Logo ROOT  

Reference Guide
TGeoChecker.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 01/11/01
3// CheckGeometry(), CheckOverlaps() by Mihaela Gheata
4
5/*************************************************************************
6 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7 * All rights reserved. *
8 * *
9 * For the licensing terms see $ROOTSYS/LICENSE. *
10 * For the list of contributors see $ROOTSYS/README/CREDITS. *
11 *************************************************************************/
12
13/** \class TGeoChecker
14\ingroup Geometry_classes
15
16Geometry checking package.
17
18TGeoChecker class provides several geometry checking methods. There are two
19types of tests that can be performed. One is based on random sampling or
20ray-tracing and provides a visual check on how navigation methods work for
21a given geometry. The second actually checks the validity of the geometry
22definition in terms of overlapping/extruding objects. Both types of checks
23can be done for a given branch (starting with a given volume) as well as for
24the geometry as a whole.
25
26#### TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z)
27
28This method can be called directly from the TGeoManager class and print a
29report on how a given point is classified by the modeller: which is the
30full path to the deepest node containing it, and the (under)estimation
31of the distance to the closest boundary (safety).
32
33#### TGeoChecker::RandomPoints(Int_t npoints)
34
35Can be called from TGeoVolume class. It first draws the volume and its
36content with the current visualization settings. Then randomly samples points
37in its bounding box, plotting in the geometry display only the points
38classified as belonging to visible volumes.
39
40#### TGeoChecker::RandomRays(Int_t nrays, Double_t startx, starty, startz)
41
42Can be called and acts in the same way as the previous, but instead of points,
43rays having random isotropic directions are generated from the given point.
44A raytracing algorithm propagates all rays until they exit geometry, plotting
45all segments crossing visible nodes in the same color as these.
46
47#### TGeoChecker::Test(Int_t npoints)
48
49Implementation of TGeoManager::Test(). Computes the time for the modeller
50to find out "Where am I?" for a given number of random points.
51
52#### TGeoChecker::LegoPlot(ntheta, themin, themax, nphi, phimin, phimax,...)
53
54Implementation of TGeoVolume::LegoPlot(). Draws a spherical radiation length
55lego plot for a given volume, in a given theta/phi range.
56
57#### TGeoChecker::Weight(Double_t precision)
58
59Implementation of TGeoVolume::Weight(). Estimates the total weight of a given
60volume by material sampling. Accepts as input the desired precision.
61*/
62
63#include "TVirtualPad.h"
64#include "TCanvas.h"
65#include "TStyle.h"
66#include "TFile.h"
67#include "TNtuple.h"
68#include "TH2.h"
69#include "TRandom3.h"
70#include "TPolyMarker3D.h"
71#include "TPolyLine3D.h"
72#include "TStopwatch.h"
73
74#include "TGeoVoxelFinder.h"
75#include "TGeoBBox.h"
76#include "TGeoPcon.h"
77#include "TGeoManager.h"
78#include "TGeoOverlap.h"
79#include "TGeoPainter.h"
80#include "TGeoChecker.h"
81
82#include "TBuffer3D.h"
83#include "TBuffer3DTypes.h"
84#include "TMath.h"
85
86#include <stdlib.h>
87
88// statics and globals
89
90ClassImp(TGeoChecker);
91
92////////////////////////////////////////////////////////////////////////////////
93/// Default constructor
94
95TGeoChecker::TGeoChecker()
96 :TObject(),
97 fGeoManager(NULL),
98 fVsafe(NULL),
99 fBuff1(NULL),
100 fBuff2(NULL),
101 fFullCheck(kFALSE),
102 fVal1(NULL),
103 fVal2(NULL),
104 fFlags(NULL),
105 fTimer(NULL),
106 fSelectedNode(NULL),
107 fNchecks(0),
108 fNmeshPoints(1000)
109{
110}
111
112////////////////////////////////////////////////////////////////////////////////
113/// Constructor for a given geometry
114
115TGeoChecker::TGeoChecker(TGeoManager *geom)
116 :TObject(),
117 fGeoManager(geom),
118 fVsafe(NULL),
119 fBuff1(NULL),
120 fBuff2(NULL),
121 fFullCheck(kFALSE),
122 fVal1(NULL),
123 fVal2(NULL),
124 fFlags(NULL),
125 fTimer(NULL),
126 fSelectedNode(NULL),
127 fNchecks(0),
128 fNmeshPoints(1000)
129{
130 fBuff1 = new TBuffer3D(TBuffer3DTypes::kGeneric,500,3*500,0,0,0,0);
131 fBuff2 = new TBuffer3D(TBuffer3DTypes::kGeneric,500,3*500,0,0,0,0);
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// Destructor
136
137TGeoChecker::~TGeoChecker()
138{
139 if (fBuff1) delete fBuff1;
140 if (fBuff2) delete fBuff2;
141 if (fTimer) delete fTimer;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145/// Print current operation progress.
146
147void TGeoChecker::OpProgress(const char *opname, Long64_t current, Long64_t size, TStopwatch *watch, Bool_t last, Bool_t refresh, const char *msg)
148{
149 static Long64_t icount = 0;
150 static TString oname;
151 static TString nname;
152 static Long64_t ocurrent = 0;
153 static Long64_t osize = 0;
154 static Int_t oseconds = 0;
155 static TStopwatch *owatch = 0;
156 static Bool_t oneoftwo = kFALSE;
157 static Int_t nrefresh = 0;
158 const char symbol[4] = {'=','\\','|','/'};
159 char progress[11] = " ";
160 Int_t ichar = icount%4;
161 TString message(msg);
162 message += " ";
163
164 if (!refresh) {
165 nrefresh = 0;
166 if (!size) return;
167 owatch = watch;
168 oname = opname;
169 ocurrent = TMath::Abs(current);
170 osize = TMath::Abs(size);
171 if (ocurrent > osize) ocurrent=osize;
172 } else {
173 nrefresh++;
174 if (!osize) return;
175 }
176 icount++;
177 Double_t time = 0.;
178 Int_t hours = 0;
179 Int_t minutes = 0;
180 Int_t seconds = 0;
181 if (owatch && !last) {
182 owatch->Stop();
183 time = owatch->RealTime();
184 hours = (Int_t)(time/3600.);
185 time -= 3600*hours;
186 minutes = (Int_t)(time/60.);
187 time -= 60*minutes;
188 seconds = (Int_t)time;
189 if (refresh) {
190 if (oseconds==seconds) {
191 owatch->Continue();
192 return;
193 }
194 oneoftwo = !oneoftwo;
195 }
196 oseconds = seconds;
197 }
198 if (refresh && oneoftwo) {
199 nname = oname;
200 if (fNchecks <= nrefresh) fNchecks = nrefresh+1;
201 Int_t pctdone = (Int_t)(100.*nrefresh/fNchecks);
202 oname = TString::Format(" == %3d%% ==", pctdone);
203 }
204 Double_t percent = 100.0*ocurrent/osize;
205 Int_t nchar = Int_t(percent/10);
206 if (nchar>10) nchar=10;
207 Int_t i;
208 for (i=0; i<nchar; i++) progress[i] = '=';
209 progress[nchar] = symbol[ichar];
210 for (i=nchar+1; i<10; i++) progress[i] = ' ';
211 progress[10] = '\0';
212 oname += " ";
213 oname.Remove(20);
214 if(size<10000) fprintf(stderr, "%s [%10s] %4lld ", oname.Data(), progress, ocurrent);
215 else if(size<100000) fprintf(stderr, "%s [%10s] %5lld ",oname.Data(), progress, ocurrent);
216 else fprintf(stderr, "%s [%10s] %7lld ",oname.Data(), progress, ocurrent);
217 if (time>0.) fprintf(stderr, "[%6.2f %%] TIME %.2d:%.2d:%.2d %s\r", percent, hours, minutes, seconds, message.Data());
218 else fprintf(stderr, "[%6.2f %%] %s\r", percent, message.Data());
219 if (refresh && oneoftwo) oname = nname;
220 if (owatch) owatch->Continue();
221 if (last) {
222 icount = 0;
223 owatch = 0;
224 ocurrent = 0;
225 osize = 0;
226 oseconds = 0;
227 oneoftwo = kFALSE;
228 nrefresh = 0;
229 fprintf(stderr, "\n");
230 }
231}
232
233////////////////////////////////////////////////////////////////////////////////
234/// Check pushes and pulls needed to cross the next boundary with respect to the
235/// position given by FindNextBoundary. If radius is not mentioned the full bounding
236/// box will be sampled.
237
238void TGeoChecker::CheckBoundaryErrors(Int_t ntracks, Double_t radius)
239{
240 TGeoVolume *tvol = fGeoManager->GetTopVolume();
241 Info("CheckBoundaryErrors", "Top volume is %s",tvol->GetName());
242 const TGeoShape *shape = tvol->GetShape();
243 TGeoBBox *box = (TGeoBBox *)shape;
244 Double_t dl[3];
245 Double_t ori[3];
246 Double_t xyz[3];
247 Double_t nxyz[3];
248 Double_t dir[3];
249 Double_t relp;
250 Char_t path[1024];
251 Char_t cdir[10];
252
253 // Tree part
254 TFile *f=new TFile("geobugs.root","recreate");
255 TTree *bug=new TTree("bug","Geometrical problems");
256 bug->Branch("pos",xyz,"xyz[3]/D");
257 bug->Branch("dir",dir,"dir[3]/D");
258 bug->Branch("push",&relp,"push/D");
259 bug->Branch("path",&path,"path/C");
260 bug->Branch("cdir",&cdir,"cdir/C");
261
262 dl[0] = box->GetDX();
263 dl[1] = box->GetDY();
264 dl[2] = box->GetDZ();
265 ori[0] = (box->GetOrigin())[0];
266 ori[1] = (box->GetOrigin())[1];
267 ori[2] = (box->GetOrigin())[2];
268 if (radius>0)
269 dl[0] = dl[1] = dl[2] = radius;
270
271 TH1::AddDirectory(kFALSE);
272 TH1F *hnew = new TH1F("hnew","Precision pushing",30,-20.,10.);
273 TH1F *hold = new TH1F("hold","Precision pulling", 30,-20.,10.);
274 TH2F *hplotS = new TH2F("hplotS","Problematic points",100,-dl[0],dl[0],100,-dl[1],dl[1]);
275 gStyle->SetOptStat(111111);
276
277 TGeoNode *node = 0;
278 Long_t igen=0;
279 Long_t itry=0;
280 Long_t n100 = ntracks/100;
281 Double_t rad = TMath::Sqrt(dl[0]*dl[0]+dl[1]*dl[1]);
282 printf("Random box : %f, %f, %f, %f, %f, %f\n", ori[0], ori[1], ori[2], dl[0], dl[1], dl[2]);
283 printf("Start... %i points\n", ntracks);
284 if (!fTimer) fTimer = new TStopwatch();
285 fTimer->Reset();
286 fTimer->Start();
287 while (igen<ntracks) {
288 Double_t phi1 = TMath::TwoPi()*gRandom->Rndm();
289 Double_t r = rad*gRandom->Rndm();
290 xyz[0] = ori[0] + r*TMath::Cos(phi1);
291 xyz[1] = ori[1] + r*TMath::Sin(phi1);
292 Double_t z = (1.-2.*gRandom->Rndm());
293 xyz[2] = ori[2]+dl[2]*z*TMath::Abs(z);
294 ++itry;
295 fGeoManager->SetCurrentPoint(xyz);
296 node = fGeoManager->FindNode();
297 if (!node || node==fGeoManager->GetTopNode()) continue;
298 ++igen;
299 if (n100 && !(igen%n100))
300 OpProgress("Sampling progress:",igen, ntracks, fTimer);
301 Double_t cost = 1.-2.*gRandom->Rndm();
302 Double_t sint = TMath::Sqrt((1.+cost)*(1.-cost));
303 Double_t phi = TMath::TwoPi()*gRandom->Rndm();
304 dir[0] = sint * TMath::Cos(phi);
305 dir[1] = sint * TMath::Sin(phi);
306 dir[2] = cost;
307 fGeoManager->SetCurrentDirection(dir);
308 fGeoManager->FindNextBoundary();
309 Double_t step = fGeoManager->GetStep();
310
311 relp = 1.e-21;
312 for(Int_t i=0; i<30; ++i) {
313 relp *=10.;
314 for(Int_t j=0; j<3; ++j) nxyz[j]=xyz[j]+step*(1.+relp)*dir[j];
315 if(!fGeoManager->IsSameLocation(nxyz[0],nxyz[1],nxyz[2])) {
316 hnew->Fill(i-20.);
317 if(i>15) {
318 const Double_t* norm = fGeoManager->FindNormal();
319 strncpy(path,fGeoManager->GetPath(),1024);
320 path[1023] = '\0';
321 Double_t dotp = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
322 printf("Forward error i=%d p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
323 i,xyz[0],xyz[1],xyz[2],step,dotp,path);
324 hplotS->Fill(xyz[0],xyz[1],(Double_t)i);
325 strncpy(cdir,"Forward",10);
326 bug->Fill();
327 }
328 break;
329 }
330 }
331
332 relp = -1.e-21;
333 for(Int_t i=0; i<30; ++i) {
334 relp *=10.;
335 for(Int_t j=0; j<3; ++j) nxyz[j]=xyz[j]+step*(1.+relp)*dir[j];
336 if(fGeoManager->IsSameLocation(nxyz[0],nxyz[1],nxyz[2])) {
337 hold->Fill(i-20.);
338 if(i>15) {
339 const Double_t* norm = fGeoManager->FindNormal();
340 strncpy(path,fGeoManager->GetPath(),1024);
341 path[1023] = '\0';
342 Double_t dotp = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
343 printf("Backward error i=%d p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
344 i,xyz[0],xyz[1],xyz[2],step,dotp,path);
345 strncpy(cdir,"Backward",10);
346 bug->Fill();
347 }
348 break;
349 }
350 }
351 }
352 fTimer->Stop();
353
354 if (itry) printf("CPU time/point = %5.2emus: Real time/point = %5.2emus\n",
355 1000000.*fTimer->CpuTime()/itry,1000000.*fTimer->RealTime()/itry);
356 bug->Write();
357 delete bug;
358 bug=0;
359 delete f;
360
361 CheckBoundaryReference();
362
363 if (itry) printf("Effic = %3.1f%%\n",(100.*igen)/itry);
364 TCanvas *c1 = new TCanvas("c1","Results",600,800);
365 c1->Divide(1,2);
366 c1->cd(1);
367 gPad->SetLogy();
368 hold->Draw();
369 c1->cd(2);
370 gPad->SetLogy();
371 hnew->Draw();
372 /*TCanvas *c3 = */new TCanvas("c3","Plot",600,600);
373 hplotS->Draw("cont0");
374}
375
376////////////////////////////////////////////////////////////////////////////////
377/// Check the boundary errors reference file created by CheckBoundaryErrors method.
378/// The shape for which the crossing failed is drawn with the starting point in red
379/// and the extrapolated point to boundary (+/- failing push/pull) in yellow.
380
381void TGeoChecker::CheckBoundaryReference(Int_t icheck)
382{
383 Double_t xyz[3];
384 Double_t nxyz[3];
385 Double_t dir[3];
386 Double_t lnext[3];
387 Double_t push;
388 Char_t path[1024];
389 Char_t cdir[10];
390 // Tree part
391 TFile *f=new TFile("geobugs.root","read");
392 TTree *bug=(TTree*)f->Get("bug");
393 bug->SetBranchAddress("pos",xyz);
394 bug->SetBranchAddress("dir",dir);
395 bug->SetBranchAddress("push",&push);
396 bug->SetBranchAddress("path",&path);
397 bug->SetBranchAddress("cdir",&cdir);
398
399 Int_t nentries = (Int_t)bug->GetEntries();
400 printf("nentries %d\n",nentries);
401 if (icheck<0) {
402 for (Int_t i=0;i<nentries;i++) {
403 bug->GetEntry(i);
404 printf("%-9s error push=%g p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
405 cdir,push,xyz[0],xyz[1],xyz[2],1.,1.,path);
406 }
407 } else {
408 if (icheck>=nentries) return;
409 Int_t idebug = TGeoManager::GetVerboseLevel();
410 TGeoManager::SetVerboseLevel(5);
411 bug->GetEntry(icheck);
412 printf("%-9s error push=%g p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
413 cdir,push,xyz[0],xyz[1],xyz[2],1.,1.,path);
414 fGeoManager->SetCurrentPoint(xyz);
415 fGeoManager->SetCurrentDirection(dir);
416 fGeoManager->FindNode();
417 TGeoNode *next = fGeoManager->FindNextBoundary();
418 Double_t step = fGeoManager->GetStep();
419 for (Int_t j=0; j<3; j++) nxyz[j]=xyz[j]+step*(1.+0.1*push)*dir[j];
420 Bool_t change = !fGeoManager->IsSameLocation(nxyz[0],nxyz[1],nxyz[2]);
421 printf("step=%g in: %s\n", step, fGeoManager->GetPath());
422 printf(" -> next = %s push=%g change=%d\n", next->GetName(),push, (UInt_t)change);
423 next->GetVolume()->InspectShape();
424 next->GetVolume()->DrawOnly();
425 if (next != fGeoManager->GetCurrentNode()) {
426 Int_t index1 = fGeoManager->GetCurrentVolume()->GetIndex(next);
427 if (index1>=0) fGeoManager->CdDown(index1);
428 }
429 TPolyMarker3D *pm = new TPolyMarker3D();
430 fGeoManager->MasterToLocal(xyz, lnext);
431 pm->SetNextPoint(lnext[0], lnext[1], lnext[2]);
432 pm->SetMarkerStyle(2);
433 pm->SetMarkerSize(0.2);
434 pm->SetMarkerColor(kRed);
435 pm->Draw("SAME");
436 TPolyMarker3D *pm1 = new TPolyMarker3D();
437 for (Int_t j=0; j<3; j++) nxyz[j]=xyz[j]+step*dir[j];
438 fGeoManager->MasterToLocal(nxyz, lnext);
439 pm1->SetNextPoint(lnext[0], lnext[1], lnext[2]);
440 pm1->SetMarkerStyle(2);
441 pm1->SetMarkerSize(0.2);
442 pm1->SetMarkerColor(kYellow);
443 pm1->Draw("SAME");
444 TGeoManager::SetVerboseLevel(idebug);
445 }
446 delete bug;
447 delete f;
448}
449
450////////////////////////////////////////////////////////////////////////////////
451/// Geometry checking. Optional overlap checkings (by sampling and by mesh). Optional
452/// boundary crossing check + timing per volume.
453///
454/// STAGE 1: extensive overlap checking by sampling per volume. Stdout need to be
455/// checked by user to get report, then TGeoVolume::CheckOverlaps(0.01, "s") can
456/// be called for the suspicious volumes.
457///
458/// STAGE 2: normal overlap checking using the shapes mesh - fills the list of
459/// overlaps.
460///
461/// STAGE 3: shooting NRAYS rays from VERTEX and counting the total number of
462/// crossings per volume (rays propagated from boundary to boundary until
463/// geometry exit). Timing computed and results stored in a histo.
464///
465/// STAGE 4: shooting 1 mil. random rays inside EACH volume and calling
466/// FindNextBoundary() + Safety() for each call. The timing is normalized by the
467/// number of crossings computed at stage 2 and presented as percentage.
468/// One can get a picture on which are the most "burned" volumes during
469/// transportation from geometry point of view. Another plot of the timing per
470/// volume vs. number of daughters is produced.
471///
472/// All histos are saved in the file statistics.root
473
474void TGeoChecker::CheckGeometryFull(Bool_t checkoverlaps, Bool_t checkcrossings, Int_t ntracks, const Double_t *vertex)
475{
476 Int_t nuid = fGeoManager->GetListOfUVolumes()->GetEntries();
477 if (!fTimer) fTimer = new TStopwatch();
478 Int_t i;
479 Double_t value;
480 fFlags = new Bool_t[nuid];
481 memset(fFlags, 0, nuid*sizeof(Bool_t));
482 TGeoVolume *vol;
483 TCanvas *c = new TCanvas("overlaps", "Overlaps by sampling", 800,800);
484
485// STAGE 1: Overlap checking by sampling per volume
486 if (checkoverlaps) {
487 printf("====================================================================\n");
488 printf("STAGE 1: Overlap checking by sampling within 10 microns\n");
489 printf("====================================================================\n");
490 fGeoManager->CheckOverlaps(0.001, "s");
491
492 // STAGE 2: Global overlap/extrusion checking
493 printf("====================================================================\n");
494 printf("STAGE 2: Global overlap/extrusion checking within 10 microns\n");
495 printf("====================================================================\n");
496 fGeoManager->CheckOverlaps(0.001);
497 }
498
499 if (!checkcrossings) {
500 delete [] fFlags;
501 fFlags = 0;
502 delete c;
503 return;
504 }
505
506 fVal1 = new Double_t[nuid];
507 fVal2 = new Double_t[nuid];
508 memset(fVal1, 0, nuid*sizeof(Double_t));
509 memset(fVal2, 0, nuid*sizeof(Double_t));
510 // STAGE 3: How many crossings per volume in a realistic case ?
511 // Ignore volumes with no daughters
512
513 // Generate rays from vertex in phi=[0,2*pi] theta=[0,pi]
514// Int_t ntracks = 1000000;
515 printf("====================================================================\n");
516 printf("STAGE 3: Propagating %i tracks starting from vertex\n and counting number of boundary crossings...\n", ntracks);
517 printf("====================================================================\n");
518 Int_t nbound = 0;
519 Double_t theta, phi;
520 Double_t point[3], dir[3];
521 memset(point, 0, 3*sizeof(Double_t));
522 if (vertex) memcpy(point, vertex, 3*sizeof(Double_t));
523
524 fTimer->Reset();
525 fTimer->Start();
526 for (i=0; i<ntracks; i++) {
527 phi = 2.*TMath::Pi()*gRandom->Rndm();
528 theta= TMath::ACos(1.-2.*gRandom->Rndm());
529 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
530 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
531 dir[2]=TMath::Cos(theta);
532 if ((i%100)==0) OpProgress("Transporting tracks",i, ntracks, fTimer);
533// if ((i%1000)==0) printf("... remaining tracks %i\n", ntracks-i);
534 nbound += PropagateInGeom(point,dir);
535 }
536 if (!nbound) {
537 printf("No boundary crossed\n");
538 return;
539 }
540 fTimer->Stop();
541 Double_t time1 = fTimer->CpuTime() *1.E6;
542 Double_t time2 = time1/ntracks;
543 Double_t time3 = time1/nbound;
544 OpProgress("Transporting tracks",ntracks, ntracks, fTimer, kTRUE);
545 printf("Time for crossing %i boundaries: %g [ms]\n", nbound, time1);
546 printf("Time per track for full geometry traversal: %g [ms], per crossing: %g [ms]\n", time2, time3);
547
548// STAGE 4: How much time per volume:
549
550 printf("====================================================================\n");
551 printf("STAGE 4: How much navigation time per volume per next+safety call\n");
552 printf("====================================================================\n");
553 TGeoIterator next(fGeoManager->GetTopVolume());
554 TGeoNode*current;
555 TString path;
556 vol = fGeoManager->GetTopVolume();
557 memset(fFlags, 0, nuid*sizeof(Bool_t));
558 TStopwatch timer;
559 timer.Start();
560 i = 0;
561 char volname[30];
562 strncpy(volname, vol->GetName(),15);
563 volname[15] = '\0';
564 OpProgress(volname,i++, nuid, &timer);
565 Score(vol, 1, TimingPerVolume(vol));
566 while ((current=next())) {
567 vol = current->GetVolume();
568 Int_t uid = vol->GetNumber();
569 if (fFlags[uid]) continue;
570 fFlags[uid] = kTRUE;
571 next.GetPath(path);
572 fGeoManager->cd(path.Data());
573 strncpy(volname, vol->GetName(),15);
574 volname[15] = '\0';
575 OpProgress(volname,i++, nuid, &timer);
576 Score(vol,1,TimingPerVolume(vol));
577 }
578 OpProgress("STAGE 4 completed",i, nuid, &timer, kTRUE);
579 // Draw some histos
580 Double_t time_tot_pertrack = 0.;
581 TCanvas *c1 = new TCanvas("c2","ncrossings",10,10,900,500);
582 c1->SetGrid();
583 c1->SetTopMargin(0.15);
584 TFile *f = new TFile("statistics.root", "RECREATE");
585 TH1F *h = new TH1F("h","number of boundary crossings per volume",3,0,3);
586 h->SetStats(0);
587 h->SetFillColor(38);
588 h->SetCanExtend(TH1::kAllAxes);
589
590 memset(fFlags, 0, nuid*sizeof(Bool_t));
591 for (i=0; i<nuid; i++) {
592 vol = fGeoManager->GetVolume(i);
593 if (!vol->GetNdaughters()) continue;
594 time_tot_pertrack += fVal1[i]*fVal2[i];
595 h->Fill(vol->GetName(), (Int_t)fVal1[i]);
596 }
597 time_tot_pertrack /= ntracks;
598 h->LabelsDeflate();
599 h->LabelsOption(">","X");
600 h->Draw();
601
602
603 TCanvas *c2 = new TCanvas("c3","time spent per volume in navigation",10,10,900,500);
604 c2->SetGrid();
605 c2->SetTopMargin(0.15);
606 TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0,100,100,0,15);
607 h2->SetStats(0);
608 h2->SetMarkerStyle(2);
609 TH1F *h1 = new TH1F("h1","percent of time spent per volume",3,0,3);
610 h1->SetStats(0);
611 h1->SetFillColor(38);
612 h1->SetCanExtend(TH1::kAllAxes);
613 for (i=0; i<nuid; i++) {
614 vol = fGeoManager->GetVolume(i);
615 if (!vol || !vol->GetNdaughters()) continue;
616 value = fVal1[i]*fVal2[i]/ntracks/time_tot_pertrack;
617 h1->Fill(vol->GetName(), value);
618 h2->Fill(vol->GetNdaughters(), fVal2[i]);
619 }
620 h1->LabelsDeflate();
621 h1->LabelsOption(">","X");
622 h1->Draw();
623 TCanvas *c3 = new TCanvas("c4","timing vs. ndaughters",10,10,900,500);
624 c3->SetGrid();
625 c3->SetTopMargin(0.15);
626 h2->Draw();
627 f->Write();
628 delete [] fFlags;
629 fFlags = 0;
630 delete [] fVal1;
631 fVal1 = 0;
632 delete [] fVal2;
633 fVal2 = 0;
634 delete fTimer;
635 fTimer = 0;
636 delete c;
637}
638
639////////////////////////////////////////////////////////////////////////////////
640/// Propagate from START along DIR from boundary to boundary until exiting
641/// geometry. Fill array of hits.
642
643Int_t TGeoChecker::PropagateInGeom(Double_t *start, Double_t *dir)
644{
645 fGeoManager->InitTrack(start, dir);
646 TGeoNode *current = 0;
647 Int_t nzero = 0;
648 Int_t nhits = 0;
649 while (!fGeoManager->IsOutside()) {
650 if (nzero>3) {
651 // Problems in trying to cross a boundary
652 printf("Error in trying to cross boundary of %s\n", current->GetName());
653 return nhits;
654 }
655 current = fGeoManager->FindNextBoundaryAndStep(TGeoShape::Big(), kFALSE);
656 if (!current || fGeoManager->IsOutside()) return nhits;
657 Double_t step = fGeoManager->GetStep();
658 if (step<2.*TGeoShape::Tolerance()) {
659 nzero++;
660 continue;
661 }
662 else nzero = 0;
663 // Generate the hit
664 nhits++;
665 TGeoVolume *vol = current->GetVolume();
666 Score(vol,0,1.);
667 Int_t iup = 1;
668 TGeoNode *mother = fGeoManager->GetMother(iup++);
669 while (mother && mother->GetVolume()->IsAssembly()) {
670 Score(mother->GetVolume(), 0, 1.);
671 mother = fGeoManager->GetMother(iup++);
672 }
673 }
674 return nhits;
675}
676
677////////////////////////////////////////////////////////////////////////////////
678/// Score a hit for VOL
679
680void TGeoChecker::Score(TGeoVolume *vol, Int_t ifield, Double_t value)
681{
682 Int_t uid = vol->GetNumber();
683 switch (ifield) {
684 case 0:
685 fVal1[uid] += value;
686 break;
687 case 1:
688 fVal2[uid] += value;
689 }
690}
691
692////////////////////////////////////////////////////////////////////////////////
693/// Set number of points to be generated on the shape outline when checking for overlaps.
694
695void TGeoChecker::SetNmeshPoints(Int_t npoints) {
696 fNmeshPoints = npoints;
697 if (npoints<1000) {
698 Error("SetNmeshPoints", "Cannot allow less than 1000 points for checking - set to 1000");
699 fNmeshPoints = 1000;
700 }
701}
702
703////////////////////////////////////////////////////////////////////////////////
704/// Compute timing per "FindNextBoundary" + "Safety" call. Volume must be
705/// in the current path.
706
707Double_t TGeoChecker::TimingPerVolume(TGeoVolume *vol)
708{
709 fTimer->Reset();
710 const TGeoShape *shape = vol->GetShape();
711 TGeoBBox *box = (TGeoBBox *)shape;
712 Double_t dx = box->GetDX();
713 Double_t dy = box->GetDY();
714 Double_t dz = box->GetDZ();
715 Double_t ox = (box->GetOrigin())[0];
716 Double_t oy = (box->GetOrigin())[1];
717 Double_t oz = (box->GetOrigin())[2];
718 Double_t point[3], dir[3], lpt[3], ldir[3];
719 Double_t pstep = 0.;
720 pstep = TMath::Max(pstep,dz);
721 Double_t theta, phi;
722 Int_t idaughter;
723 fTimer->Start();
724 Bool_t inside;
725 for (Int_t i=0; i<1000000; i++) {
726 lpt[0] = ox-dx+2*dx*gRandom->Rndm();
727 lpt[1] = oy-dy+2*dy*gRandom->Rndm();
728 lpt[2] = oz-dz+2*dz*gRandom->Rndm();
729 fGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
730 fGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
731 phi = 2*TMath::Pi()*gRandom->Rndm();
732 theta= TMath::ACos(1.-2.*gRandom->Rndm());
733 ldir[0]=TMath::Sin(theta)*TMath::Cos(phi);
734 ldir[1]=TMath::Sin(theta)*TMath::Sin(phi);
735 ldir[2]=TMath::Cos(theta);
736 fGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir,dir);
737 fGeoManager->SetCurrentDirection(dir);
738 fGeoManager->SetStep(pstep);
739 fGeoManager->ResetState();
740 inside = kTRUE;
741 // dist = TGeoShape::Big();
742 if (!vol->IsAssembly()) {
743 inside = vol->Contains(lpt);
744 if (!inside) {
745 // dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep);
746 // if (dist>=pstep) continue;
747 } else {
748 vol->GetShape()->DistFromInside(lpt,ldir,3,pstep);
749 }
750
751 if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside);
752 }
753 if (vol->GetNdaughters()) {
754 fGeoManager->Safety();
755 fGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE);
756 }
757 }
758 fTimer->Stop();
759 Double_t time_per_track = fTimer->CpuTime();
760 Int_t uid = vol->GetNumber();
761 Int_t ncrossings = (Int_t)fVal1[uid];
762 if (!vol->GetNdaughters())
763 printf("Time for volume %s (shape=%s): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->GetShape()->GetName(), time_per_track, vol->GetNdaughters(), ncrossings);
764 else
765 printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(), time_per_track, vol->GetNdaughters(), ncrossings);
766 return time_per_track;
767}
768
769////////////////////////////////////////////////////////////////////////////////
770/// Shoot nrays with random directions from starting point (startx, starty, startz)
771/// in the reference frame of this volume. Track each ray until exiting geometry, then
772/// shoot backwards from exiting point and compare boundary crossing points.
773
774void TGeoChecker::CheckGeometry(Int_t nrays, Double_t startx, Double_t starty, Double_t startz) const
775{
776 Int_t i, j;
777 Double_t start[3], end[3];
778 Double_t dir[3];
779 Double_t dummy[3];
780 Double_t eps = 0.;
781 Double_t *array1 = new Double_t[3*1000];
782 Double_t *array2 = new Double_t[3*1000];
783 TObjArray *pma = new TObjArray();
784 TPolyMarker3D *pm;
785 pm = new TPolyMarker3D();
786 pm->SetMarkerColor(2); // error > eps
787 pm->SetMarkerStyle(8);
788 pm->SetMarkerSize(0.4);
789 pma->AddAt(pm, 0);
790 pm = new TPolyMarker3D();
791 pm->SetMarkerColor(4); // point not found
792 pm->SetMarkerStyle(8);
793 pm->SetMarkerSize(0.4);
794 pma->AddAt(pm, 1);
795 pm = new TPolyMarker3D();
796 pm->SetMarkerColor(6); // extra point back
797 pm->SetMarkerStyle(8);
798 pm->SetMarkerSize(0.4);
799 pma->AddAt(pm, 2);
800 Int_t nelem1, nelem2;
801 Int_t dim1=1000, dim2=1000;
802 if ((startx==0) && (starty==0) && (startz==0)) eps=1E-3;
803 start[0] = startx+eps;
804 start[1] = starty+eps;
805 start[2] = startz+eps;
806 Int_t n10=nrays/10;
807 Double_t theta,phi;
808 Double_t dw, dwmin, dx, dy, dz;
809 Int_t ist1, ist2, ifound;
810 for (i=0; i<nrays; i++) {
811 if (n10) {
812 if ((i%n10) == 0) printf("%i percent\n", Int_t(100*i/nrays));
813 }
814 phi = 2*TMath::Pi()*gRandom->Rndm();
815 theta= TMath::ACos(1.-2.*gRandom->Rndm());
816 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
817 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
818 dir[2]=TMath::Cos(theta);
819 // shoot direct ray
820 nelem1=nelem2=0;
821// printf("DIRECT %i\n", i);
822 array1 = ShootRay(&start[0], dir[0], dir[1], dir[2], array1, nelem1, dim1);
823 if (!nelem1) continue;
824// for (j=0; j<nelem1; j++) printf("%i : %f %f %f\n", j, array1[3*j], array1[3*j+1], array1[3*j+2]);
825 memcpy(&end[0], &array1[3*(nelem1-1)], 3*sizeof(Double_t));
826 // shoot ray backwards
827// printf("BACK %i\n", i);
828 array2 = ShootRay(&end[0], -dir[0], -dir[1], -dir[2], array2, nelem2, dim2, &start[0]);
829 if (!nelem2) {
830 printf("#### NOTHING BACK ###########################\n");
831 for (j=0; j<nelem1; j++) {
832 pm = (TPolyMarker3D*)pma->At(0);
833 pm->SetNextPoint(array1[3*j], array1[3*j+1], array1[3*j+2]);
834 }
835 continue;
836 }
837// printf("BACKWARDS\n");
838 Int_t k=nelem2>>1;
839 for (j=0; j<k; j++) {
840 memcpy(&dummy[0], &array2[3*j], 3*sizeof(Double_t));
841 memcpy(&array2[3*j], &array2[3*(nelem2-1-j)], 3*sizeof(Double_t));
842 memcpy(&array2[3*(nelem2-1-j)], &dummy[0], 3*sizeof(Double_t));
843 }
844// for (j=0; j<nelem2; j++) printf("%i : %f ,%f ,%f \n", j, array2[3*j], array2[3*j+1], array2[3*j+2]);
845 if (nelem1!=nelem2) printf("### DIFFERENT SIZES : nelem1=%i nelem2=%i ##########\n", nelem1, nelem2);
846 ist1 = ist2 = 0;
847 // check first match
848
849 dx = array1[3*ist1]-array2[3*ist2];
850 dy = array1[3*ist1+1]-array2[3*ist2+1];
851 dz = array1[3*ist1+2]-array2[3*ist2+2];
852 dw = dx*dir[0]+dy*dir[1]+dz*dir[2];
853 fGeoManager->SetCurrentPoint(&array1[3*ist1]);
854 fGeoManager->FindNode();
855// printf("%i : %s (%g, %g, %g)\n", ist1, fGeoManager->GetPath(),
856// array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2]);
857 if (TMath::Abs(dw)<1E-4) {
858// printf(" matching %i (%g, %g, %g)\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
859 ist2++;
860 } else {
861 printf("### NOT MATCHING %i f:(%f, %f, %f) b:(%f %f %f) DCLOSE=%f\n", ist2, array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2], array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2],dw);
862 pm = (TPolyMarker3D*)pma->At(0);
863 pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
864 if (dw<0) {
865 // first boundary missed on way back
866 } else {
867 // first boundary different on way back
868 ist2++;
869 }
870 }
871
872 while ((ist1<nelem1-1) && (ist2<nelem2)) {
873 fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
874 fGeoManager->FindNode();
875// printf("%i : %s (%g, %g, %g)\n", ist1+1, fGeoManager->GetPath(),
876// array1[3*ist1+3], array1[3*ist1+4], array1[3*ist1+5]);
877
878 dx = array1[3*ist1+3]-array1[3*ist1];
879 dy = array1[3*ist1+4]-array1[3*ist1+1];
880 dz = array1[3*ist1+5]-array1[3*ist1+2];
881 // distance to next point
882 dwmin = dx+dir[0]+dy*dir[1]+dz*dir[2];
883 while (ist2<nelem2) {
884 ifound = 0;
885 dx = array2[3*ist2]-array1[3*ist1];
886 dy = array2[3*ist2+1]-array1[3*ist1+1];
887 dz = array2[3*ist2+2]-array1[3*ist1+2];
888 dw = dx+dir[0]+dy*dir[1]+dz*dir[2];
889 if (TMath::Abs(dw-dwmin)<1E-4) {
890 ist1++;
891 ist2++;
892 break;
893 }
894 if (dw<dwmin) {
895 // point found on way back. Check if close enough to ist1+1
896 ifound++;
897 dw = dwmin-dw;
898 if (dw<1E-4) {
899 // point is matching ist1+1
900// printf(" matching %i (%g, %g, %g) DCLOSE=%g\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2], dw);
901 ist2++;
902 ist1++;
903 break;
904 } else {
905 // extra boundary found on way back
906 fGeoManager->SetCurrentPoint(&array2[3*ist2]);
907 fGeoManager->FindNode();
908 pm = (TPolyMarker3D*)pma->At(2);
909 pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
910 printf("### EXTRA BOUNDARY %i : %s found at DCLOSE=%f\n", ist2, fGeoManager->GetPath(), dw);
911 ist2++;
912 continue;
913 }
914 } else {
915 if (!ifound) {
916 // point ist1+1 not found on way back
917 fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
918 fGeoManager->FindNode();
919 pm = (TPolyMarker3D*)pma->At(1);
920 pm->SetNextPoint(array2[3*ist1+3], array2[3*ist1+4], array2[3*ist1+5]);
921 printf("### BOUNDARY MISSED BACK #########################\n");
922 ist1++;
923 break;
924 } else {
925 ist1++;
926 break;
927 }
928 }
929 }
930 }
931 }
932 pm = (TPolyMarker3D*)pma->At(0);
933 pm->Draw("SAME");
934 pm = (TPolyMarker3D*)pma->At(1);
935 pm->Draw("SAME");
936 pm = (TPolyMarker3D*)pma->At(2);
937 pm->Draw("SAME");
938 if (gPad) {
939 gPad->Modified();
940 gPad->Update();
941 }
942 delete [] array1;
943 delete [] array2;
944 delete pma; // markers are drawn on the pad
945}
946
947////////////////////////////////////////////////////////////////////////////////
948/// Clean-up the mesh of pcon/pgon from useless points
949
950void TGeoChecker::CleanPoints(Double_t *points, Int_t &numPoints) const
951{
952 Int_t ipoint = 0;
953 Int_t j, k=0;
954 Double_t rsq;
955 for (Int_t i=0; i<numPoints; i++) {
956 j = 3*i;
957 rsq = points[j]*points[j]+points[j+1]*points[j+1];
958 if (rsq < 1.e-10) continue;
959 points[k] = points[j];
960 points[k+1] = points[j+1];
961 points[k+2] = points[j+2];
962 ipoint++;
963 k = 3*ipoint;
964 }
965 numPoints = ipoint;
966}
967
968////////////////////////////////////////////////////////////////////////////////
969/// Check if the 2 non-assembly volume candidates overlap/extrude. Returns overlap object.
970
971TGeoOverlap *TGeoChecker::MakeCheckOverlap(const char *name, TGeoVolume *vol1, TGeoVolume *vol2, TGeoMatrix *mat1, TGeoMatrix *mat2, Bool_t isovlp, Double_t ovlp)
972{
973 TGeoOverlap *nodeovlp = 0;
974 Int_t numPoints1 = fBuff1->NbPnts();
975 Int_t numSegs1 = fBuff1->NbSegs();
976 Int_t numPols1 = fBuff1->NbPols();
977 Int_t numPoints2 = fBuff2->NbPnts();
978 Int_t numSegs2 = fBuff2->NbSegs();
979 Int_t numPols2 = fBuff2->NbPols();
980 Int_t ip;
981 Bool_t extrude, isextrusion, isoverlapping;
982 Double_t *points1 = fBuff1->fPnts;
983 Double_t *points2 = fBuff2->fPnts;
984 Double_t local[3], local1[3];
985 Double_t point[3];
986 Double_t safety = TGeoShape::Big();
987 Double_t tolerance = TGeoShape::Tolerance();
988 if (vol1->IsAssembly() || vol2->IsAssembly()) return nodeovlp;
989 TGeoShape *shape1 = vol1->GetShape();
990 TGeoShape *shape2 = vol2->GetShape();
991 OpProgress("refresh", 0,0,NULL,kFALSE,kTRUE);
992 shape1->GetMeshNumbers(numPoints1, numSegs1, numPols1);
993 if (fBuff1->fID != (TObject*)shape1) {
994 // Fill first buffer.
995 fBuff1->SetRawSizes(TMath::Max(numPoints1,fNmeshPoints), 3*TMath::Max(numPoints1,fNmeshPoints), 0, 0, 0, 0);
996 points1 = fBuff1->fPnts;
997 if (shape1->GetPointsOnSegments(fNmeshPoints, points1)) {
998 numPoints1 = fNmeshPoints;
999 } else {
1000 shape1->SetPoints(points1);
1001 }
1002// if (shape1->InheritsFrom(TGeoPcon::Class())) {
1003// CleanPoints(points1, numPoints1);
1004// fBuff1->SetRawSizes(numPoints1, 3*numPoints1,0, 0, 0, 0);
1005// }
1006 fBuff1->fID = shape1;
1007 }
1008 shape2->GetMeshNumbers(numPoints2, numSegs2, numPols2);
1009 if (fBuff2->fID != (TObject*)shape2) {
1010 // Fill second buffer.
1011 fBuff2->SetRawSizes(TMath::Max(numPoints2,fNmeshPoints), 3*TMath::Max(numPoints2,fNmeshPoints), 0, 0, 0,0);
1012 points2 = fBuff2->fPnts;
1013 if (shape2->GetPointsOnSegments(fNmeshPoints, points2)) {
1014 numPoints2 = fNmeshPoints;
1015 } else {
1016 shape2->SetPoints(points2);
1017 }
1018// if (shape2->InheritsFrom(TGeoPcon::Class())) {
1019// CleanPoints(points2, numPoints2);
1020// fBuff1->SetRawSizes(numPoints2, 3*numPoints2,0,0,0,0);
1021// }
1022 fBuff2->fID = shape2;
1023 }
1024
1025 if (!isovlp) {
1026 // Extrusion case. Test vol2 extrude vol1.
1027 isextrusion=kFALSE;
1028 // loop all points of the daughter
1029 for (ip=0; ip<numPoints2; ip++) {
1030 memcpy(local, &points2[3*ip], 3*sizeof(Double_t));
1031 if (TMath::Abs(local[0])<tolerance && TMath::Abs(local[1])<tolerance) continue;
1032 mat2->LocalToMaster(local, point);
1033 mat1->MasterToLocal(point, local);
1034 extrude = !shape1->Contains(local);
1035 if (extrude) {
1036 safety = shape1->Safety(local, kFALSE);
1037 if (safety<ovlp) extrude=kFALSE;
1038 }
1039 if (extrude) {
1040 if (!isextrusion) {
1041 isextrusion = kTRUE;
1042 nodeovlp = new TGeoOverlap(name, vol1, vol2, mat1,mat2,kFALSE,safety);
1043 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1044 fGeoManager->AddOverlap(nodeovlp);
1045 } else {
1046 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1047 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1048 }
1049 }
1050 }
1051 // loop all points of the mother
1052 for (ip=0; ip<numPoints1; ip++) {
1053 memcpy(local, &points1[3*ip], 3*sizeof(Double_t));
1054 if (local[0]<1e-10 && local[1]<1e-10) continue;
1055 mat1->LocalToMaster(local, point);
1056 mat2->MasterToLocal(point, local1);
1057 extrude = shape2->Contains(local1);
1058 if (extrude) {
1059 // skip points on mother mesh that have no neighborhood outside mother
1060 safety = shape1->Safety(local,kTRUE);
1061 if (safety>1E-6) {
1062 extrude = kFALSE;
1063 } else {
1064 safety = shape2->Safety(local1,kTRUE);
1065 if (safety<ovlp) extrude=kFALSE;
1066 }
1067 }
1068 if (extrude) {
1069 if (!isextrusion) {
1070 isextrusion = kTRUE;
1071 nodeovlp = new TGeoOverlap(name, vol1,vol2,mat1,mat2,kFALSE,safety);
1072 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1073 fGeoManager->AddOverlap(nodeovlp);
1074 } else {
1075 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1076 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1077 }
1078 }
1079 }
1080 return nodeovlp;
1081 }
1082 // Check overlap
1083 Bool_t overlap;
1084 overlap = kFALSE;
1085 isoverlapping = kFALSE;
1086 // loop all points of first candidate
1087 for (ip=0; ip<numPoints1; ip++) {
1088 memcpy(local, &points1[3*ip], 3*sizeof(Double_t));
1089 if (local[0]<1e-10 && local[1]<1e-10) continue;
1090 mat1->LocalToMaster(local, point);
1091 mat2->MasterToLocal(point, local); // now point in local reference of second
1092 overlap = shape2->Contains(local);
1093 if (overlap) {
1094 safety = shape2->Safety(local, kTRUE);
1095 if (safety<ovlp) overlap=kFALSE;
1096 }
1097 if (overlap) {
1098 if (!isoverlapping) {
1099 isoverlapping = kTRUE;
1100 nodeovlp = new TGeoOverlap(name,vol1,vol2,mat1,mat2,kTRUE,safety);
1101 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1102 fGeoManager->AddOverlap(nodeovlp);
1103 } else {
1104 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1105 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1106 }
1107 }
1108 }
1109 // loop all points of second candidate
1110 for (ip=0; ip<numPoints2; ip++) {
1111 memcpy(local, &points2[3*ip], 3*sizeof(Double_t));
1112 if (local[0]<1e-10 && local[1]<1e-10) continue;
1113 mat2->LocalToMaster(local, point);
1114 mat1->MasterToLocal(point, local); // now point in local reference of first
1115 overlap = shape1->Contains(local);
1116 if (overlap) {
1117 safety = shape1->Safety(local, kTRUE);
1118 if (safety<ovlp) overlap=kFALSE;
1119 }
1120 if (overlap) {
1121 if (!isoverlapping) {
1122 isoverlapping = kTRUE;
1123 nodeovlp = new TGeoOverlap(name,vol1,vol2,mat1,mat2,kTRUE,safety);
1124 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1125 fGeoManager->AddOverlap(nodeovlp);
1126 } else {
1127 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1128 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1129 }
1130 }
1131 }
1132 return nodeovlp;
1133}
1134
1135////////////////////////////////////////////////////////////////////////////////
1136/// Check illegal overlaps for volume VOL within a limit OVLP by sampling npoints
1137/// inside the volume shape.
1138
1139void TGeoChecker::CheckOverlapsBySampling(TGeoVolume *vol, Double_t /* ovlp */, Int_t npoints) const
1140{
1141 Int_t nd = vol->GetNdaughters();
1142 if (nd<2) return;
1143 TGeoVoxelFinder *voxels = vol->GetVoxels();
1144 if (!voxels) return;
1145 if (voxels->NeedRebuild()) {
1146 voxels->Voxelize();
1147 vol->FindOverlaps();
1148 }
1149 TGeoBBox *box = (TGeoBBox*)vol->GetShape();
1150 TGeoShape *shape;
1151 TGeoNode *node;
1152 Double_t dx = box->GetDX();
1153 Double_t dy = box->GetDY();
1154 Double_t dz = box->GetDZ();
1155 Double_t pt[3];
1156 Double_t local[3];
1157 Int_t *check_list = 0;
1158 Int_t ncheck = 0;
1159 const Double_t *orig = box->GetOrigin();
1160 Int_t ipoint = 0;
1161 Int_t itry = 0;
1162 Int_t iovlp = 0;
1163 Int_t id=0, id0=0, id1=0;
1164 Bool_t in, incrt;
1165 Double_t safe;
1166 TString name1 = "";
1167 TString name2 = "";
1168 TGeoOverlap **flags = 0;
1169 TGeoNode *node1, *node2;
1170 Int_t novlps = 0;
1171 TGeoHMatrix mat1, mat2;
1172// Int_t tid = TGeoManager::ThreadId();
1173 TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
1174 TGeoStateInfo &td = *nav->GetCache()->GetInfo();
1175 while (ipoint < npoints) {
1176 // Shoot randomly in the bounding box.
1177 pt[0] = orig[0] - dx + 2.*dx*gRandom->Rndm();
1178 pt[1] = orig[1] - dy + 2.*dy*gRandom->Rndm();
1179 pt[2] = orig[2] - dz + 2.*dz*gRandom->Rndm();
1180 if (!vol->Contains(pt)) {
1181 itry++;
1182 if (itry>10000 && !ipoint) {
1183 Error("CheckOverlapsBySampling", "No point inside volume!!! - aborting");
1184 break;
1185 }
1186 continue;
1187 }
1188 // Check if the point is inside one or more daughters
1189 in = kFALSE;
1190 ipoint++;
1191 check_list = voxels->GetCheckList(pt, ncheck, td);
1192 if (!check_list || ncheck<2) continue;
1193 for (id=0; id<ncheck; id++) {
1194 id0 = check_list[id];
1195 node = vol->GetNode(id0);
1196 // Ignore MANY's
1197 if (node->IsOverlapping()) continue;
1198 node->GetMatrix()->MasterToLocal(pt,local);
1199 shape = node->GetVolume()->GetShape();
1200 incrt = shape->Contains(local);
1201 if (!incrt) continue;
1202 if (!in) {
1203 in = kTRUE;
1204 id1 = id0;
1205 continue;
1206 }
1207 // The point is inside 2 or more daughters, check safety
1208 safe = shape->Safety(local, kTRUE);
1209// if (safe < ovlp) continue;
1210 // We really have found an overlap -> store the point in a container
1211 iovlp++;
1212 if (!novlps) {
1213 flags = new TGeoOverlap*[nd*nd];
1214 memset(flags, 0, nd*nd*sizeof(TGeoOverlap*));
1215 }
1216 TGeoOverlap *nodeovlp = flags[nd*id1+id0];
1217 if (!nodeovlp) {
1218 novlps++;
1219 node1 = vol->GetNode(id1);
1220 name1 = node1->GetName();
1221 mat1 = node1->GetMatrix();
1222 Int_t cindex = node1->GetVolume()->GetCurrentNodeIndex();
1223 while (cindex >= 0) {
1224 node1 = node1->GetVolume()->GetNode(cindex);
1225 name1 += TString::Format("/%s", node1->GetName());
1226 mat1.Multiply(node1->GetMatrix());
1227 cindex = node1->GetVolume()->GetCurrentNodeIndex();
1228 }
1229 node2 = vol->GetNode(id0);
1230 name2 = node2->GetName();
1231 mat2 = node2->GetMatrix();
1232 cindex = node2->GetVolume()->GetCurrentNodeIndex();
1233 while (cindex >= 0) {
1234 node2 = node2->GetVolume()->GetNode(cindex);
1235 name2 += TString::Format("/%s", node2->GetName());
1236 mat2.Multiply(node2->GetMatrix());
1237 cindex = node2->GetVolume()->GetCurrentNodeIndex();
1238 }
1239 nodeovlp = new TGeoOverlap(TString::Format("Volume %s: node %s overlapping %s",
1240 vol->GetName(), name1.Data(), name2.Data()), node1->GetVolume(),node2->GetVolume(),
1241 &mat1,&mat2, kTRUE, safe);
1242 flags[nd*id1+id0] = nodeovlp;
1243 fGeoManager->AddOverlap(nodeovlp);
1244 }
1245 // Max 100 points per marker
1246 if (nodeovlp->GetPolyMarker()->GetN()<100) nodeovlp->SetNextPoint(pt[0],pt[1],pt[2]);
1247 if (nodeovlp->GetOverlap()<safe) nodeovlp->SetOverlap(safe);
1248 }
1249 }
1250 nav->GetCache()->ReleaseInfo();
1251 if (flags) delete [] flags;
1252 if (!novlps) return;
1253 Double_t capacity = vol->GetShape()->Capacity();
1254 capacity *= Double_t(iovlp)/Double_t(npoints);
1255 Double_t err = 1./TMath::Sqrt(Double_t(iovlp));
1256 Info("CheckOverlapsBySampling", "#Found %d overlaps adding-up to %g +/- %g [cm3] for daughters of %s",
1257 novlps, capacity, err*capacity, vol->GetName());
1258}
1259
1260////////////////////////////////////////////////////////////////////////////////
1261/// Compute number of overlaps combinations to check per volume
1262
1263Int_t TGeoChecker::NChecksPerVolume(TGeoVolume *vol)
1264{
1265 if (vol->GetFinder()) return 0;
1266 UInt_t nd = vol->GetNdaughters();
1267 if (!nd) return 0;
1268 Bool_t is_assembly = vol->IsAssembly();
1269 TGeoIterator next1(vol);
1270 TGeoIterator next2(vol);
1271 Int_t nchecks = 0;
1272 TGeoNode *node;
1273 UInt_t id;
1274 if (!is_assembly) {
1275 while ((node=next1())) {
1276 if (node->IsOverlapping()) {
1277 next1.Skip();
1278 continue;
1279 }
1280 if (!node->GetVolume()->IsAssembly()) {
1281 nchecks++;
1282 next1.Skip();
1283 }
1284 }
1285 }
1286 // now check if the daughters overlap with each other
1287 if (nd<2) return nchecks;
1288 TGeoVoxelFinder *vox = vol->GetVoxels();
1289 if (!vox) return nchecks;
1290
1291
1292 TGeoNode *node1, *node01, *node02;
1293 Int_t novlp;
1294 Int_t *ovlps;
1295 Int_t ko;
1296 UInt_t io;
1297 for (id=0; id<nd; id++) {
1298 node01 = vol->GetNode(id);
1299 if (node01->IsOverlapping()) continue;
1300 vox->FindOverlaps(id);
1301 ovlps = node01->GetOverlaps(novlp);
1302 if (!ovlps) continue;
1303 for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
1304 io = ovlps[ko]; // (node1, shaped, matrix1, points, fBuff1)
1305 if (io<=id) continue;
1306 node02 = vol->GetNode(io);
1307 if (node02->IsOverlapping()) continue;
1308
1309 // We have to check node against node1, but they may be assemblies
1310 if (node01->GetVolume()->IsAssembly()) {
1311 next1.Reset(node01->GetVolume());
1312 while ((node=next1())) {
1313 if (!node->GetVolume()->IsAssembly()) {
1314 if (node02->GetVolume()->IsAssembly()) {
1315 next2.Reset(node02->GetVolume());
1316 while ((node1=next2())) {
1317 if (!node1->GetVolume()->IsAssembly()) {
1318 nchecks++;
1319 next2.Skip();
1320 }
1321 }
1322 } else {
1323 nchecks++;
1324 }
1325 next1.Skip();
1326 }
1327 }
1328 } else {
1329 // node not assembly
1330 if (node02->GetVolume()->IsAssembly()) {
1331 next2.Reset(node02->GetVolume());
1332 while ((node1=next2())) {
1333 if (!node1->GetVolume()->IsAssembly()) {
1334 nchecks++;
1335 next2.Skip();
1336 }
1337 }
1338 } else {
1339 // node1 also not assembly
1340 nchecks++;
1341 }
1342 }
1343 }
1344 node01->SetOverlaps(0,0);
1345 }
1346 return nchecks;
1347}
1348
1349////////////////////////////////////////////////////////////////////////////////
1350/// Check illegal overlaps for volume VOL within a limit OVLP.
1351
1352void TGeoChecker::CheckOverlaps(const TGeoVolume *vol, Double_t ovlp, Option_t *option)
1353{
1354 if (vol->GetFinder()) return;
1355 UInt_t nd = vol->GetNdaughters();
1356 if (!nd) return;
1357 TGeoShape::SetTransform(gGeoIdentity);
1358 fNchecks = NChecksPerVolume((TGeoVolume*)vol);
1359 Bool_t sampling = kFALSE;
1360 TString opt(option);
1361 opt.ToLower();
1362 if (opt.Contains("s")) sampling = kTRUE;
1363 if (opt.Contains("f")) fFullCheck = kTRUE;
1364 else fFullCheck = kFALSE;
1365 if (sampling) {
1366 opt = opt.Strip(TString::kLeading, 's');
1367 Int_t npoints = atoi(opt.Data());
1368 if (!npoints) npoints = 1000000;
1369 CheckOverlapsBySampling((TGeoVolume*)vol, ovlp, npoints);
1370 return;
1371 }
1372 Bool_t is_assembly = vol->IsAssembly();
1373 TGeoIterator next1((TGeoVolume*)vol);
1374 TGeoIterator next2((TGeoVolume*)vol);
1375 TString path;
1376 // first, test if daughters extrude their container
1377 TGeoNode * node, *nodecheck;
1378 TGeoChecker *checker = (TGeoChecker*)this;
1379
1380// TGeoOverlap *nodeovlp = 0;
1381 UInt_t id;
1382 Int_t level;
1383// Check extrusion only for daughters of a non-assembly volume
1384 if (!is_assembly) {
1385 while ((node=next1())) {
1386 if (node->IsOverlapping()) {
1387 next1.Skip();
1388 continue;
1389 }
1390 if (!node->GetVolume()->IsAssembly()) {
1391 if (fSelectedNode) {
1392 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1393 if ((fSelectedNode != node) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1394 next1.Skip();
1395 continue;
1396 }
1397 if (node != fSelectedNode) {
1398 level = next1.GetLevel();
1399 while ((nodecheck=next1.GetNode(level--))) {
1400 if (nodecheck == fSelectedNode) break;
1401 }
1402 if (!nodecheck) {
1403 next1.Skip();
1404 continue;
1405 }
1406 }
1407 }
1408 next1.GetPath(path);
1409 checker->MakeCheckOverlap(TString::Format("%s extruded by: %s", vol->GetName(),path.Data()),
1410 (TGeoVolume*)vol,node->GetVolume(),gGeoIdentity,(TGeoMatrix*)next1.GetCurrentMatrix(),kFALSE,ovlp);
1411 next1.Skip();
1412 }
1413 }
1414 }
1415
1416 // now check if the daughters overlap with each other
1417 if (nd<2) return;
1418 TGeoVoxelFinder *vox = vol->GetVoxels();
1419 if (!vox) {
1420 Warning("CheckOverlaps", "Volume %s with %i daughters but not voxelized", vol->GetName(),nd);
1421 return;
1422 }
1423 if (vox->NeedRebuild()) {
1424 vox->Voxelize();
1425 vol->FindOverlaps();
1426 }
1427 TGeoNode *node1, *node01, *node02;
1428 TGeoHMatrix hmat1, hmat2;
1429 TString path1;
1430 Int_t novlp;
1431 Int_t *ovlps;
1432 Int_t ko;
1433 UInt_t io;
1434 for (id=0; id<nd; id++) {
1435 node01 = vol->GetNode(id);
1436 if (node01->IsOverlapping()) continue;
1437 vox->FindOverlaps(id);
1438 ovlps = node01->GetOverlaps(novlp);
1439 if (!ovlps) continue;
1440 next1.SetTopName(node01->GetName());
1441 path = node01->GetName();
1442 for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
1443 io = ovlps[ko]; // (node1, shaped, matrix1, points, fBuff1)
1444 if (io<=id) continue;
1445 node02 = vol->GetNode(io);
1446 if (node02->IsOverlapping()) continue;
1447 // Try to fasten-up things...
1448// if (!TGeoBBox::AreOverlapping((TGeoBBox*)node01->GetVolume()->GetShape(), node01->GetMatrix(),
1449// (TGeoBBox*)node02->GetVolume()->GetShape(), node02->GetMatrix())) continue;
1450 next2.SetTopName(node02->GetName());
1451 path1 = node02->GetName();
1452
1453 // We have to check node against node1, but they may be assemblies
1454 if (node01->GetVolume()->IsAssembly()) {
1455 next1.Reset(node01->GetVolume());
1456 while ((node=next1())) {
1457 if (!node->GetVolume()->IsAssembly()) {
1458 next1.GetPath(path);
1459 hmat1 = node01->GetMatrix();
1460 hmat1 *= *next1.GetCurrentMatrix();
1461 if (node02->GetVolume()->IsAssembly()) {
1462 next2.Reset(node02->GetVolume());
1463 while ((node1=next2())) {
1464 if (!node1->GetVolume()->IsAssembly()) {
1465 if (fSelectedNode) {
1466 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1467 if ((fSelectedNode != node) && (fSelectedNode != node1) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1468 next2.Skip();
1469 continue;
1470 }
1471 if ((fSelectedNode != node1) && (fSelectedNode != node)) {
1472 level = next2.GetLevel();
1473 while ((nodecheck=next2.GetNode(level--))) {
1474 if (nodecheck == fSelectedNode) break;
1475 }
1476 if (node02 == fSelectedNode) nodecheck = node02;
1477 if (!nodecheck) {
1478 level = next1.GetLevel();
1479 while ((nodecheck=next1.GetNode(level--))) {
1480 if (nodecheck == fSelectedNode) break;
1481 }
1482 }
1483 if (node01 == fSelectedNode) nodecheck = node01;
1484 if (!nodecheck) {
1485 next2.Skip();
1486 continue;
1487 }
1488 }
1489 }
1490 next2.GetPath(path1);
1491 hmat2 = node02->GetMatrix();
1492 hmat2 *= *next2.GetCurrentMatrix();
1493 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1494 node->GetVolume(),node1->GetVolume(),&hmat1,&hmat2,kTRUE,ovlp);
1495 next2.Skip();
1496 }
1497 }
1498 } else {
1499 if (fSelectedNode) {
1500 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1501 if ((fSelectedNode != node) && (fSelectedNode != node02) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1502 next1.Skip();
1503 continue;
1504 }
1505 if ((fSelectedNode != node) && (fSelectedNode != node02)) {
1506 level = next1.GetLevel();
1507 while ((nodecheck=next1.GetNode(level--))) {
1508 if (nodecheck == fSelectedNode) break;
1509 }
1510 if (node01 == fSelectedNode) nodecheck = node01;
1511 if (!nodecheck) {
1512 next1.Skip();
1513 continue;
1514 }
1515 }
1516 }
1517 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1518 node->GetVolume(),node02->GetVolume(),&hmat1,node02->GetMatrix(),kTRUE,ovlp);
1519 }
1520 next1.Skip();
1521 }
1522 }
1523 } else {
1524 // node not assembly
1525 if (node02->GetVolume()->IsAssembly()) {
1526 next2.Reset(node02->GetVolume());
1527 while ((node1=next2())) {
1528 if (!node1->GetVolume()->IsAssembly()) {
1529 if (fSelectedNode) {
1530 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1531 if ((fSelectedNode != node1) && (fSelectedNode != node01) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1532 next2.Skip();
1533 continue;
1534 }
1535 if ((fSelectedNode != node1) && (fSelectedNode != node01)) {
1536 level = next2.GetLevel();
1537 while ((nodecheck=next2.GetNode(level--))) {
1538 if (nodecheck == fSelectedNode) break;
1539 }
1540 if (node02 == fSelectedNode) nodecheck = node02;
1541 if (!nodecheck) {
1542 next2.Skip();
1543 continue;
1544 }
1545 }
1546 }
1547 next2.GetPath(path1);
1548 hmat2 = node02->GetMatrix();
1549 hmat2 *= *next2.GetCurrentMatrix();
1550 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1551 node01->GetVolume(),node1->GetVolume(),node01->GetMatrix(),&hmat2,kTRUE,ovlp);
1552 next2.Skip();
1553 }
1554 }
1555 } else {
1556 // node1 also not assembly
1557 if (fSelectedNode && (fSelectedNode != node01) && (fSelectedNode != node02)) continue;
1558 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1559 node01->GetVolume(),node02->GetVolume(),node01->GetMatrix(),node02->GetMatrix(),kTRUE,ovlp);
1560 }
1561 }
1562 }
1563 node01->SetOverlaps(0,0);
1564 }
1565}
1566
1567////////////////////////////////////////////////////////////////////////////////
1568/// Print the current list of overlaps held by the manager class.
1569
1570void TGeoChecker::PrintOverlaps() const
1571{
1572 TIter next(fGeoManager->GetListOfOverlaps());
1573 TGeoOverlap *ov;
1574 printf("=== Overlaps for %s ===\n", fGeoManager->GetName());
1575 while ((ov=(TGeoOverlap*)next())) ov->PrintInfo();
1576}
1577
1578////////////////////////////////////////////////////////////////////////////////
1579/// Draw point (x,y,z) over the picture of the daughters of the volume containing this point.
1580/// Generates a report regarding the path to the node containing this point and the distance to
1581/// the closest boundary.
1582
1583void TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z, Option_t *)
1584{
1585 Double_t point[3];
1586 Double_t local[3];
1587 point[0] = x;
1588 point[1] = y;
1589 point[2] = z;
1590 TGeoVolume *vol = fGeoManager->GetTopVolume();
1591 if (fVsafe) {
1592 TGeoNode *old = fVsafe->GetNode("SAFETY_1");
1593 if (old) fVsafe->GetNodes()->RemoveAt(vol->GetNdaughters()-1);
1594 }
1595// if (vol != fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
1596 TGeoNode *node = fGeoManager->FindNode(point[0], point[1], point[2]);
1597 fGeoManager->MasterToLocal(point, local);
1598 // get current node
1599 printf("=== Check current point : (%g, %g, %g) ===\n", point[0], point[1], point[2]);
1600 printf(" - path : %s\n", fGeoManager->GetPath());
1601 // get corresponding volume
1602 if (node) vol = node->GetVolume();
1603 // compute safety distance (distance to boundary ignored)
1604 Double_t close = fGeoManager->Safety();
1605 printf("Safety radius : %f\n", close);
1606 if (close>1E-4) {
1607 TGeoVolume *sph = fGeoManager->MakeSphere("SAFETY", vol->GetMedium(), 0, close, 0,180,0,360);
1608 sph->SetLineColor(2);
1609 sph->SetLineStyle(3);
1610 vol->AddNode(sph,1,new TGeoTranslation(local[0], local[1], local[2]));
1611 fVsafe = vol;
1612 }
1613 TPolyMarker3D *pm = new TPolyMarker3D();
1614 pm->SetMarkerColor(2);
1615 pm->SetMarkerStyle(8);
1616 pm->SetMarkerSize(0.5);
1617 pm->SetNextPoint(local[0], local[1], local[2]);
1618 if (vol->GetNdaughters()<2) fGeoManager->SetTopVisible();
1619 else fGeoManager->SetTopVisible(kFALSE);
1620 fGeoManager->SetVisLevel(1);
1621 if (!vol->IsVisible()) vol->SetVisibility(kTRUE);
1622 vol->Draw();
1623 pm->Draw("SAME");
1624 gPad->Modified();
1625 gPad->Update();
1626}
1627
1628////////////////////////////////////////////////////////////////////////////////
1629/// Test for shape navigation methods. Summary for test numbers:
1630/// - 1: DistFromInside/Outside. Sample points inside the shape. Generate
1631/// directions randomly in cos(theta). Compute DistFromInside and move the
1632/// point with bigger distance. Compute DistFromOutside back from new point.
1633/// Plot d-(d1+d2)
1634/// - 2: Safety test. Sample points inside the bounding and compute safety. Generate
1635/// directions randomly in cos(theta) and compute distance to boundary. Check if
1636/// distance to boundary is bigger than safety
1637
1638void TGeoChecker::CheckShape(TGeoShape *shape, Int_t testNo, Int_t nsamples, Option_t *option)
1639{
1640 switch (testNo) {
1641 case 1:
1642 ShapeDistances(shape, nsamples, option);
1643 break;
1644 case 2:
1645 ShapeSafety(shape, nsamples, option);
1646 break;
1647 case 3:
1648 ShapeNormal(shape, nsamples, option);
1649 break;
1650 default:
1651 Error("CheckShape", "Test number %d not existent", testNo);
1652 }
1653}
1654
1655////////////////////////////////////////////////////////////////////////////////
1656/// Test TGeoShape::DistFromInside/Outside. Sample points inside the shape. Generate
1657/// directions randomly in cos(theta). Compute d1 = DistFromInside and move the
1658/// point on the boundary. Compute DistFromOutside and propagate with d2 making sure that
1659/// the shape is not re-entered. Swap direction and call DistFromOutside that
1660/// should fall back on the same point on the boundary (at d2). Propagate back on boundary
1661/// then compute DistFromInside that should be bigger than d1.
1662/// Plot d-(d1+d2)
1663
1664void TGeoChecker::ShapeDistances(TGeoShape *shape, Int_t nsamples, Option_t *)
1665{
1666 Double_t dx = ((TGeoBBox*)shape)->GetDX();
1667 Double_t dy = ((TGeoBBox*)shape)->GetDY();
1668 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1669 Double_t dmax = 2.*TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1670 Double_t d1, d2, dmove, dnext;
1671 Int_t itot = 0;
1672 // Number of tracks shot for every point inside the shape
1673 const Int_t kNtracks = 1000;
1674 Int_t n10 = nsamples/10;
1675 Int_t i,j;
1676 Double_t point[3], pnew[3];
1677 Double_t dir[3], dnew[3];
1678 Double_t theta, phi, delta;
1679 TPolyMarker3D *pmfrominside = 0;
1680 TPolyMarker3D *pmfromoutside = 0;
1681 TH1D *hist = new TH1D("hTest1", "Residual distance from inside/outside",200,-20, 0);
1682 hist->GetXaxis()->SetTitle("delta[cm] - first bin=overflow");
1683 hist->GetYaxis()->SetTitle("count");
1684 hist->SetMarkerStyle(kFullCircle);
1685
1686 if (!fTimer) fTimer = new TStopwatch();
1687 fTimer->Reset();
1688 fTimer->Start();
1689 while (itot<nsamples) {
1690 Bool_t inside = kFALSE;
1691 while (!inside) {
1692 point[0] = gRandom->Uniform(-dx,dx);
1693 point[1] = gRandom->Uniform(-dy,dy);
1694 point[2] = gRandom->Uniform(-dz,dz);
1695 inside = shape->Contains(point);
1696 }
1697 itot++;
1698 if (n10) {
1699 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1700 }
1701 for (i=0; i<kNtracks; i++) {
1702 phi = 2*TMath::Pi()*gRandom->Rndm();
1703 theta= TMath::ACos(1.-2.*gRandom->Rndm());
1704 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
1705 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
1706 dir[2]=TMath::Cos(theta);
1707 dmove = dmax;
1708 // We have track direction, compute distance from inside
1709 d1 = shape->DistFromInside(point,dir,3);
1710 if (d1>dmove || d1<TGeoShape::Tolerance()) {
1711 // Bad distance or bbox size, to debug
1712 new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1713 shape->Draw();
1714 printf("DistFromInside: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) %f/%f(max)\n",
1715 point[0],point[1],point[2],dir[0],dir[1],dir[2], d1,dmove);
1716 pmfrominside = new TPolyMarker3D(2);
1717 pmfrominside->SetMarkerColor(kRed);
1718 pmfrominside->SetMarkerStyle(24);
1719 pmfrominside->SetMarkerSize(0.4);
1720 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1721 for (j=0; j<3; j++) pnew[j] = point[j] + d1*dir[j];
1722 pmfrominside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1723 pmfrominside->Draw();
1724 return;
1725 }
1726 // Propagate BEFORE the boundary and make sure that DistFromOutside
1727 // does not return 0 (!!!)
1728 // Check if there is a second crossing
1729 for (j=0; j<3; j++) pnew[j] = point[j] + (d1-TGeoShape::Tolerance())*dir[j];
1730 dnext = shape->DistFromOutside(pnew,dir,3);
1731 if (d1+dnext<dmax) dmove = d1+0.5*dnext;
1732 // Move point and swap direction
1733 for (j=0; j<3; j++) {
1734 pnew[j] = point[j] + dmove*dir[j];
1735 dnew[j] = -dir[j];
1736 }
1737 // Compute now distance from outside
1738 d2 = shape->DistFromOutside(pnew,dnew,3);
1739 delta = dmove-d1-d2;
1740 if (TMath::Abs(delta)>1E-6 || dnext<2.*TGeoShape::Tolerance()) {
1741 // Error->debug this
1742 new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1743 shape->Draw();
1744 printf("Error: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d2=%f dmove=%f\n",
1745 point[0],point[1],point[2],dir[0],dir[1],dir[2], d1,d2,dmove);
1746 if (dnext<2.*TGeoShape::Tolerance()) {
1747 printf(" (*)DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) dnext = %f\n",
1748 point[0]+(d1-TGeoShape::Tolerance())*dir[0],
1749 point[1]+(d1-TGeoShape::Tolerance())*dir[1],
1750 point[2]+(d1-TGeoShape::Tolerance())*dir[2], dir[0],dir[1],dir[2],dnext);
1751 } else {
1752 printf(" DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) dnext = %f\n",
1753 point[0]+d1*dir[0],point[1]+d1*dir[1], point[2]+d1*dir[2], dir[0],dir[1],dir[2],dnext);
1754 }
1755 printf(" DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) = %f\n",
1756 pnew[0],pnew[1],pnew[2],dnew[0],dnew[1],dnew[2], d2);
1757 pmfrominside = new TPolyMarker3D(2);
1758 pmfrominside->SetMarkerStyle(24);
1759 pmfrominside->SetMarkerSize(0.4);
1760 pmfrominside->SetMarkerColor(kRed);
1761 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1762 for (j=0; j<3; j++) point[j] += d1*dir[j];
1763 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1764 pmfrominside->Draw();
1765 pmfromoutside = new TPolyMarker3D(2);
1766 pmfromoutside->SetMarkerStyle(20);
1767 pmfromoutside->SetMarkerStyle(7);
1768 pmfromoutside->SetMarkerSize(0.3);
1769 pmfromoutside->SetMarkerColor(kBlue);
1770 pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1771 for (j=0; j<3; j++) pnew[j] += d2*dnew[j];
1772 if (d2<1E10) pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1773 pmfromoutside->Draw();
1774 return;
1775 }
1776 // Compute distance from inside which should be bigger than d1
1777 for (j=0; j<3; j++) pnew[j] += (d2-TGeoShape::Tolerance())*dnew[j];
1778 dnext = shape->DistFromInside(pnew,dnew,3);
1779 if (dnext<d1-TGeoShape::Tolerance() || dnext>dmax) {
1780 new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1781 shape->Draw();
1782 printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d1p=%f\n",
1783 pnew[0],pnew[1],pnew[2],dnew[0],dnew[1],dnew[2],d1,dnext);
1784 pmfrominside = new TPolyMarker3D(2);
1785 pmfrominside->SetMarkerStyle(24);
1786 pmfrominside->SetMarkerSize(0.4);
1787 pmfrominside->SetMarkerColor(kRed);
1788 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1789 for (j=0; j<3; j++) point[j] += d1*dir[j];
1790 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1791 pmfrominside->Draw();
1792 pmfromoutside = new TPolyMarker3D(2);
1793 pmfromoutside->SetMarkerStyle(20);
1794 pmfromoutside->SetMarkerStyle(7);
1795 pmfromoutside->SetMarkerSize(0.3);
1796 pmfromoutside->SetMarkerColor(kBlue);
1797 pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1798 for (j=0; j<3; j++) pnew[j] += dnext*dnew[j];
1799 if (d2<1E10) pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1800 pmfromoutside->Draw();
1801 return;
1802 }
1803 if (TMath::Abs(delta) < 1E-20) delta = 1E-30;
1804 hist->Fill(TMath::Max(TMath::Log(TMath::Abs(delta)),-20.));
1805 }
1806 }
1807 fTimer->Stop();
1808 fTimer->Print();
1809 new TCanvas("Test01", "Residuals DistFromInside/Outside", 800, 600);
1810 hist->Draw();
1811}
1812
1813////////////////////////////////////////////////////////////////////////////////
1814/// Check of validity of safe distance for a given shape.
1815/// Sample points inside the 2x bounding box and compute safety. Generate
1816/// directions randomly in cos(theta) and compute distance to boundary. Check if
1817/// distance to boundary is bigger than safety.
1818
1819void TGeoChecker::ShapeSafety(TGeoShape *shape, Int_t nsamples, Option_t *)
1820{
1821 Double_t dx = ((TGeoBBox*)shape)->GetDX();
1822 Double_t dy = ((TGeoBBox*)shape)->GetDY();
1823 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1824 // Number of tracks shot for every point inside the shape
1825 const Int_t kNtracks = 1000;
1826 Int_t n10 = nsamples/10;
1827 Int_t i;
1828 Double_t dist;
1829 Double_t point[3];
1830 Double_t dir[3];
1831 Double_t theta, phi;
1832 TPolyMarker3D *pm1 = 0;
1833 TPolyMarker3D *pm2 = 0;
1834 if (!fTimer) fTimer = new TStopwatch();
1835 fTimer->Reset();
1836 fTimer->Start();
1837 Int_t itot = 0;
1838 while (itot<nsamples) {
1839 Bool_t inside = kFALSE;
1840 point[0] = gRandom->Uniform(-2*dx,2*dx);
1841 point[1] = gRandom->Uniform(-2*dy,2*dy);
1842 point[2] = gRandom->Uniform(-2*dz,2*dz);
1843 inside = shape->Contains(point);
1844 Double_t safe = shape->Safety(point, inside);
1845 itot++;
1846 if (n10) {
1847 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1848 }
1849 for (i=0; i<kNtracks; i++) {
1850 phi = 2*TMath::Pi()*gRandom->Rndm();
1851 theta= TMath::ACos(1.-2.*gRandom->Rndm());
1852 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
1853 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
1854 dir[2]=TMath::Cos(theta);
1855 if (inside) dist = shape->DistFromInside(point,dir,3);
1856 else dist = shape->DistFromOutside(point,dir,3);
1857 if (dist<safe) {
1858 printf("Error safety (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) safe=%f dist=%f\n",
1859 point[0],point[1],point[2], dir[0], dir[1], dir[2], safe, dist);
1860 new TCanvas("shape02", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1861 shape->Draw();
1862 pm1 = new TPolyMarker3D(2);
1863 pm1->SetMarkerStyle(24);
1864 pm1->SetMarkerSize(0.4);
1865 pm1->SetMarkerColor(kRed);
1866 pm1->SetNextPoint(point[0],point[1],point[2]);
1867 pm1->SetNextPoint(point[0]+safe*dir[0],point[1]+safe*dir[1],point[2]+safe*dir[2]);
1868 pm1->Draw();
1869 pm2 = new TPolyMarker3D(1);
1870 pm2->SetMarkerStyle(7);
1871 pm2->SetMarkerSize(0.3);
1872 pm2->SetMarkerColor(kBlue);
1873 pm2->SetNextPoint(point[0]+dist*dir[0],point[1]+dist*dir[1],point[2]+dist*dir[2]);
1874 pm2->Draw();
1875 return;
1876 }
1877 }
1878 }
1879 fTimer->Stop();
1880 fTimer->Print();
1881}
1882
1883////////////////////////////////////////////////////////////////////////////////
1884/// Check of validity of the normal for a given shape.
1885/// Sample points inside the shape. Generate directions randomly in cos(theta)
1886/// and propagate to boundary. Compute normal and safety at crossing point, plot
1887/// the point and generate a random direction so that (dir) dot (norm) <0.
1888
1889void TGeoChecker::ShapeNormal(TGeoShape *shape, Int_t nsamples, Option_t *)
1890{
1891 Double_t dx = ((TGeoBBox*)shape)->GetDX();
1892 Double_t dy = ((TGeoBBox*)shape)->GetDY();
1893 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1894 Double_t dmax = 2.*TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1895 // Number of tracks shot for every point inside the shape
1896 const Int_t kNtracks = 1000;
1897 Int_t n10 = nsamples/10;
1898 Int_t itot = 0, errcnt = 0, errsame=0;
1899 Int_t i;
1900 Double_t dist, olddist, safe, dot;
1901 Double_t theta, phi, ndotd;
1902 Double_t *spoint = new Double_t[3*nsamples];
1903 Double_t *sdir = new Double_t[3*nsamples];
1904 while (itot<nsamples) {
1905 Bool_t inside = kFALSE;
1906 while (!inside) {
1907 spoint[3*itot] = gRandom->Uniform(-dx,dx);
1908 spoint[3*itot+1] = gRandom->Uniform(-dy,dy);
1909 spoint[3*itot+2] = gRandom->Uniform(-dz,dz);
1910 inside = shape->Contains(&spoint[3*itot]);
1911 }
1912 phi = 2*TMath::Pi()*gRandom->Rndm();
1913 theta= TMath::ACos(1.-2.*gRandom->Rndm());
1914 sdir[3*itot] = TMath::Sin(theta)*TMath::Cos(phi);
1915 sdir[3*itot+1] = TMath::Sin(theta)*TMath::Sin(phi);
1916 sdir[3*itot+2] = TMath::Cos(theta);
1917 itot++;
1918 }
1919 Double_t point[3],newpoint[3], oldpoint[3];
1920 Double_t dir[3], olddir[3];
1921 Double_t norm[3], newnorm[3], oldnorm[3];
1922 TCanvas *errcanvas = 0;
1923 TPolyMarker3D *pm1 = 0;
1924 TPolyMarker3D *pm2 = 0;
1925 pm2 = new TPolyMarker3D();
1926// pm2->SetMarkerStyle(24);
1927 pm2->SetMarkerSize(0.2);
1928 pm2->SetMarkerColor(kBlue);
1929 if (!fTimer) fTimer = new TStopwatch();
1930 fTimer->Reset();
1931 fTimer->Start();
1932 for (itot = 0; itot<nsamples; itot++) {
1933 if (n10) {
1934 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1935 }
1936 oldnorm[0] = oldnorm[1] = oldnorm[2] = 0.;
1937 olddist = 0.;
1938 for (Int_t j=0; j<3; j++) {
1939 oldpoint[j] = point[j] = spoint[3*itot+j];
1940 olddir[j] = dir[j] = sdir[3*itot+j];
1941 }
1942 for (i=0; i<kNtracks; i++) {
1943 if (errcnt>0) break;
1944 dist = shape->DistFromInside(point,dir,3);
1945 for (Int_t j=0; j<3; j++) {
1946 newpoint[j] = point[j] + dist*dir[j];
1947 }
1948 shape->ComputeNormal(newpoint,dir,newnorm);
1949
1950 dot = olddir[0]*oldnorm[0]+olddir[1]*oldnorm[1]+ olddir[2]*oldnorm[2];
1951 if (!shape->Contains(point) && shape->Safety(point,kFALSE)>1.E-3) {
1952 errcnt++;
1953 printf("Error point outside (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
1954 point[0],point[1],point[2], dir[0], dir[1], dir[2], dist, olddist);
1955 printf(" old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n",
1956 oldpoint[0],oldpoint[1],oldpoint[2], olddir[0], olddir[1], olddir[2]);
1957 if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1958 if (!pm1) {
1959 pm1 = new TPolyMarker3D();
1960 pm1->SetMarkerStyle(24);
1961 pm1->SetMarkerSize(0.4);
1962 pm1->SetMarkerColor(kRed);
1963 }
1964 pm1->SetNextPoint(point[0],point[1],point[2]);
1965 pm1->SetNextPoint(oldpoint[0],oldpoint[1],oldpoint[2]);
1966 break;
1967 }
1968 if ((dist<TGeoShape::Tolerance() && olddist*dot>1.E-3) || dist>dmax) {
1969 errsame++;
1970 if (errsame>1) {
1971 errcnt++;
1972 printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
1973 point[0],point[1],point[2], dir[0], dir[1], dir[2], dist, olddist);
1974 printf(" new norm: (%g, %g, %g)\n", newnorm[0], newnorm[1], newnorm[2]);
1975 printf(" old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n",
1976 oldpoint[0],oldpoint[1],oldpoint[2], olddir[0], olddir[1], olddir[2]);
1977 printf(" old norm: (%g, %g, %g)\n", oldnorm[0], oldnorm[1], oldnorm[2]);
1978 if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1979 if (!pm1) {
1980 pm1 = new TPolyMarker3D();
1981 pm1->SetMarkerStyle(24);
1982 pm1->SetMarkerSize(0.4);
1983 pm1->SetMarkerColor(kRed);
1984 }
1985 pm1->SetNextPoint(point[0],point[1],point[2]);
1986 pm1->SetNextPoint(oldpoint[0],oldpoint[1],oldpoint[2]);
1987 break;
1988 }
1989 } else errsame = 0;
1990
1991 olddist = dist;
1992 for (Int_t j=0; j<3; j++) {
1993 oldpoint[j] = point[j];
1994 point[j] += dist*dir[j];
1995 }
1996 safe = shape->Safety(point, kTRUE);
1997 if (safe>1.E-3) {
1998 errcnt++;
1999 printf("Error safety (%19.15f, %19.15f, %19.15f) safe=%g\n",
2000 point[0],point[1],point[2], safe);
2001 if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
2002 if (!pm1) {
2003 pm1 = new TPolyMarker3D();
2004 pm1->SetMarkerStyle(24);
2005 pm1->SetMarkerSize(0.4);
2006 pm1->SetMarkerColor(kRed);
2007 }
2008 pm1->SetNextPoint(point[0],point[1],point[2]);
2009 break;
2010 }
2011 // Compute normal
2012 shape->ComputeNormal(point,dir,norm);
2013 memcpy(oldnorm, norm, 3*sizeof(Double_t));
2014 memcpy(olddir, dir, 3*sizeof(Double_t));
2015 while (1) {
2016 phi = 2*TMath::Pi()*gRandom->Rndm();
2017 theta= TMath::ACos(1.-2.*gRandom->Rndm());
2018 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
2019 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
2020 dir[2]=TMath::Cos(theta);
2021 ndotd = dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2];
2022 if (ndotd<0) break; // backwards, still inside shape
2023 }
2024 if ((itot%10) == 0) pm2->SetNextPoint(point[0],point[1],point[2]);
2025 }
2026 }
2027 fTimer->Stop();
2028 fTimer->Print();
2029 if (errcanvas) {
2030 shape->Draw();
2031 pm1->Draw();
2032 }
2033
2034 new TCanvas("shape03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
2035 pm2->Draw();
2036 delete [] spoint;
2037 delete [] sdir;
2038}
2039
2040////////////////////////////////////////////////////////////////////////////////
2041/// Generate a lego plot fot the top volume, according to option.
2042
2043TH2F *TGeoChecker::LegoPlot(Int_t ntheta, Double_t themin, Double_t themax,
2044 Int_t nphi, Double_t phimin, Double_t phimax,
2045 Double_t /*rmin*/, Double_t /*rmax*/, Option_t *option)
2046{
2047 TH2F *hist = new TH2F("lego", option, nphi, phimin, phimax, ntheta, themin, themax);
2048
2049 Double_t degrad = TMath::Pi()/180.;
2050 Double_t theta, phi, step, matprop, x;
2051 Double_t start[3];
2052 Double_t dir[3];
2053 TGeoNode *startnode, *endnode;
2054 Int_t i; // loop index for phi
2055 Int_t j; // loop index for theta
2056 Int_t ntot = ntheta * nphi;
2057 Int_t n10 = ntot/10;
2058 Int_t igen = 0, iloop=0;
2059 printf("=== Lego plot sph. => nrays=%i\n", ntot);
2060 for (i=1; i<=nphi; i++) {
2061 for (j=1; j<=ntheta; j++) {
2062 igen++;
2063 if (n10) {
2064 if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/ntot));
2065 }
2066 x = 0;
2067 theta = hist->GetYaxis()->GetBinCenter(j);
2068 phi = hist->GetXaxis()->GetBinCenter(i)+1E-3;
2069 start[0] = start[1] = start[2] = 1E-3;
2070 dir[0]=TMath::Sin(theta*degrad)*TMath::Cos(phi*degrad);
2071 dir[1]=TMath::Sin(theta*degrad)*TMath::Sin(phi*degrad);
2072 dir[2]=TMath::Cos(theta*degrad);
2073 fGeoManager->InitTrack(&start[0], &dir[0]);
2074 startnode = fGeoManager->GetCurrentNode();
2075 if (fGeoManager->IsOutside()) startnode=0;
2076 if (startnode) {
2077 matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
2078 } else {
2079 matprop = 0.;
2080 }
2081 fGeoManager->FindNextBoundary();
2082// fGeoManager->IsStepEntering();
2083 // find where we end-up
2084 endnode = fGeoManager->Step();
2085 step = fGeoManager->GetStep();
2086 while (step<1E10) {
2087 // now see if we can make an other step
2088 iloop=0;
2089 while (!fGeoManager->IsEntering()) {
2090 iloop++;
2091 fGeoManager->SetStep(1E-3);
2092 step += 1E-3;
2093 endnode = fGeoManager->Step();
2094 }
2095 if (iloop>1000) printf("%i steps\n", iloop);
2096 if (matprop>0) {
2097 x += step/matprop;
2098 }
2099 if (endnode==0 && step>1E10) break;
2100 // generate an extra step to cross boundary
2101 startnode = endnode;
2102 if (startnode) {
2103 matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
2104 } else {
2105 matprop = 0.;
2106 }
2107
2108 fGeoManager->FindNextBoundary();
2109 endnode = fGeoManager->Step();
2110 step = fGeoManager->GetStep();
2111 }
2112 hist->Fill(phi, theta, x);
2113 }
2114 }
2115 return hist;
2116}
2117
2118////////////////////////////////////////////////////////////////////////////////
2119/// Draw random points in the bounding box of a volume.
2120
2121void TGeoChecker::RandomPoints(TGeoVolume *vol, Int_t npoints, Option_t *option)
2122{
2123 if (!vol) return;
2124 vol->VisibleDaughters(kTRUE);
2125 vol->Draw();
2126 TString opt = option;
2127 opt.ToLower();
2128 TObjArray *pm = new TObjArray(128);
2129 TPolyMarker3D *marker = 0;
2130 const TGeoShape *shape = vol->GetShape();
2131 TGeoBBox *box = (TGeoBBox *)shape;
2132 Double_t dx = box->GetDX();
2133 Double_t dy = box->GetDY();
2134 Double_t dz = box->GetDZ();
2135 Double_t ox = (box->GetOrigin())[0];
2136 Double_t oy = (box->GetOrigin())[1];
2137 Double_t oz = (box->GetOrigin())[2];
2138 Double_t *xyz = new Double_t[3];
2139 printf("Random box : %f, %f, %f\n", dx, dy, dz);
2140 TGeoNode *node = 0;
2141 printf("Start... %i points\n", npoints);
2142 Int_t i=0;
2143 Int_t igen=0;
2144 Int_t ic = 0;
2145 Int_t n10 = npoints/10;
2146 Double_t ratio=0;
2147 while (igen<npoints) {
2148 xyz[0] = ox-dx+2*dx*gRandom->Rndm();
2149 xyz[1] = oy-dy+2*dy*gRandom->Rndm();
2150 xyz[2] = oz-dz+2*dz*gRandom->Rndm();
2151 fGeoManager->SetCurrentPoint(xyz);
2152 igen++;
2153 if (n10) {
2154 if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/npoints));
2155 }
2156 node = fGeoManager->FindNode();
2157 if (!node) continue;
2158 if (!node->IsOnScreen()) continue;
2159 // draw only points in overlapping/non-overlapping volumes
2160 if (opt.Contains("many") && !node->IsOverlapping()) continue;
2161 if (opt.Contains("only") && node->IsOverlapping()) continue;
2162 ic = node->GetColour();
2163 if ((ic<0) || (ic>=128)) ic = 1;
2164 marker = (TPolyMarker3D*)pm->At(ic);
2165 if (!marker) {
2166 marker = new TPolyMarker3D();
2167 marker->SetMarkerColor(ic);
2168// marker->SetMarkerStyle(8);
2169// marker->SetMarkerSize(0.4);
2170 pm->AddAt(marker, ic);
2171 }
2172 marker->SetNextPoint(xyz[0], xyz[1], xyz[2]);
2173 i++;
2174 }
2175 printf("Number of visible points : %i\n", i);
2176 if (igen) ratio = (Double_t)i/(Double_t)igen;
2177 printf("efficiency : %g\n", ratio);
2178 for (Int_t m=0; m<128; m++) {
2179 marker = (TPolyMarker3D*)pm->At(m);
2180 if (marker) marker->Draw("SAME");
2181 }
2182 fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
2183 printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
2184 printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
2185 delete pm;
2186 delete [] xyz;
2187}
2188
2189////////////////////////////////////////////////////////////////////////////////
2190/// Randomly shoot nrays from point (startx,starty,startz) and plot intersections
2191/// with surfaces for current top node.
2192
2193void TGeoChecker::RandomRays(Int_t nrays, Double_t startx, Double_t starty, Double_t startz, const char *target_vol, Bool_t check_norm)
2194{
2195 TObjArray *pm = new TObjArray(128);
2196 TString starget = target_vol;
2197 TPolyLine3D *line = 0;
2198 TPolyLine3D *normline = 0;
2199 TGeoVolume *vol=fGeoManager->GetTopVolume();
2200// vol->VisibleDaughters(kTRUE);
2201
2202 Bool_t random = kFALSE;
2203 if (nrays<=0) {
2204 nrays = 100000;
2205 random = kTRUE;
2206 }
2207 Double_t start[3];
2208 Double_t dir[3];
2209 const Double_t *point = fGeoManager->GetCurrentPoint();
2210 vol->Draw();
2211 printf("Start... %i rays\n", nrays);
2212 TGeoNode *startnode, *endnode;
2213 Bool_t vis1,vis2;
2214 Int_t i=0;
2215 Int_t ipoint, inull;
2216 Int_t itot=0;
2217 Int_t n10=nrays/10;
2218 Double_t theta,phi, step, normlen;
2219 Double_t ox = ((TGeoBBox*)vol->GetShape())->GetOrigin()[0];
2220 Double_t oy = ((TGeoBBox*)vol->GetShape())->GetOrigin()[1];
2221 Double_t oz = ((TGeoBBox*)vol->GetShape())->GetOrigin()[2];
2222 Double_t dx = ((TGeoBBox*)vol->GetShape())->GetDX();
2223 Double_t dy = ((TGeoBBox*)vol->GetShape())->GetDY();
2224 Double_t dz = ((TGeoBBox*)vol->GetShape())->GetDZ();
2225 normlen = TMath::Max(dx,dy);
2226 normlen = TMath::Max(normlen,dz);
2227 normlen *= 0.05;
2228 while (itot<nrays) {
2229 itot++;
2230 inull = 0;
2231 ipoint = 0;
2232 if (n10) {
2233 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nrays));
2234 }
2235 if (random) {
2236 start[0] = ox-dx+2*dx*gRandom->Rndm();
2237 start[1] = oy-dy+2*dy*gRandom->Rndm();
2238 start[2] = oz-dz+2*dz*gRandom->Rndm();
2239 } else {
2240 start[0] = startx;
2241 start[1] = starty;
2242 start[2] = startz;
2243 }
2244 phi = 2*TMath::Pi()*gRandom->Rndm();
2245 theta= TMath::ACos(1.-2.*gRandom->Rndm());
2246 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
2247 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
2248 dir[2]=TMath::Cos(theta);
2249 startnode = fGeoManager->InitTrack(start[0],start[1],start[2], dir[0],dir[1],dir[2]);
2250 line = 0;
2251 if (fGeoManager->IsOutside()) startnode=0;
2252 vis1 = kFALSE;
2253 if (target_vol) {
2254 if (startnode && starget==startnode->GetVolume()->GetName()) vis1 = kTRUE;
2255 } else {
2256 if (startnode && startnode->IsOnScreen()) vis1 = kTRUE;
2257 }
2258 if (vis1) {
2259 line = new TPolyLine3D(2);
2260 line->SetLineColor(startnode->GetVolume()->GetLineColor());
2261 line->SetPoint(ipoint++, start[0], start[1], start[2]);
2262 i++;
2263 pm->Add(line);
2264 }
2265 while ((endnode=fGeoManager->FindNextBoundaryAndStep())) {
2266 step = fGeoManager->GetStep();
2267 if (step<TGeoShape::Tolerance()) inull++;
2268 else inull = 0;
2269 if (inull>5) break;
2270 const Double_t *normal = 0;
2271 if (check_norm) {
2272 normal = fGeoManager->FindNormalFast();
2273 if (!normal) break;
2274 }
2275 vis2 = kFALSE;
2276 if (target_vol) {
2277 if (starget==endnode->GetVolume()->GetName()) vis2 = kTRUE;
2278 } else if (endnode->IsOnScreen()) vis2 = kTRUE;
2279 if (ipoint>0) {
2280 // old visible node had an entry point -> finish segment
2281 line->SetPoint(ipoint, point[0], point[1], point[2]);
2282 if (!vis2 && check_norm) {
2283 normline = new TPolyLine3D(2);
2284 normline->SetLineColor(kBlue);
2285 normline->SetLineWidth(1);
2286 normline->SetPoint(0, point[0], point[1], point[2]);
2287 normline->SetPoint(1, point[0]+normal[0]*normlen,
2288 point[1]+normal[1]*normlen,
2289 point[2]+normal[2]*normlen);
2290 pm->Add(normline);
2291 }
2292 ipoint = 0;
2293 line = 0;
2294 }
2295 if (vis2) {
2296 // create new segment
2297 line = new TPolyLine3D(2);
2298 line->SetLineColor(endnode->GetVolume()->GetLineColor());
2299 line->SetPoint(ipoint++, point[0], point[1], point[2]);
2300 i++;
2301 if (check_norm) {
2302 normline = new TPolyLine3D(2);
2303 normline->SetLineColor(kBlue);
2304 normline->SetLineWidth(2);
2305 normline->SetPoint(0, point[0], point[1], point[2]);
2306 normline->SetPoint(1, point[0]+normal[0]*normlen,
2307 point[1]+normal[1]*normlen,
2308 point[2]+normal[2]*normlen);
2309 }
2310 pm->Add(line);
2311 if (!random) pm->Add(normline);
2312 }
2313 }
2314 }
2315 // draw all segments
2316 for (Int_t m=0; m<pm->GetEntriesFast(); m++) {
2317 line = (TPolyLine3D*)pm->At(m);
2318 if (line) line->Draw("SAME");
2319 }
2320 printf("number of segments : %i\n", i);
2321// fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
2322// printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
2323// printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
2324 delete pm;
2325}
2326
2327////////////////////////////////////////////////////////////////////////////////
2328/// shoot npoints randomly in a box of 1E-5 around current point.
2329/// return minimum distance to points outside
2330/// make sure that path to current node is updated
2331/// get the response of tgeo
2332
2333TGeoNode *TGeoChecker::SamplePoints(Int_t npoints, Double_t &dist, Double_t epsil,
2334 const char* g3path)
2335{
2336 TGeoNode *node = fGeoManager->FindNode();
2337 TGeoNode *nodegeo = 0;
2338 TGeoNode *nodeg3 = 0;
2339 TGeoNode *solg3 = 0;
2340 if (!node) {dist=-1; return 0;}
2341 Bool_t hasg3 = kFALSE;
2342 if (strlen(g3path)) hasg3 = kTRUE;
2343 TString geopath = fGeoManager->GetPath();
2344 dist = 1E10;
2345 TString common = "";
2346 // cd to common path
2347 Double_t point[3];
2348 Double_t closest[3];
2349 TGeoNode *node1 = 0;
2350 TGeoNode *node_close = 0;
2351 dist = 1E10;
2352 Double_t dist1 = 0;
2353 // initialize size of random box to epsil
2354 Double_t eps[3];
2355 eps[0] = epsil; eps[1]=epsil; eps[2]=epsil;
2356 const Double_t *pointg = fGeoManager->GetCurrentPoint();
2357 if (hasg3) {
2358 TString spath = geopath;
2359 TString name = "";
2360 Int_t index=0;
2361 while (index>=0) {
2362 index = spath.Index("/", index+1);
2363 if (index>0) {
2364 name = spath(0, index);
2365 if (strstr(g3path, name.Data())) {
2366 common = name;
2367 continue;
2368 } else break;
2369 }
2370 }
2371 // if g3 response was given, cd to common path
2372 if (strlen(common.Data())) {
2373 while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
2374 nodegeo = fGeoManager->GetCurrentNode();
2375 fGeoManager->CdUp();
2376 }
2377 fGeoManager->cd(g3path);
2378 solg3 = fGeoManager->GetCurrentNode();
2379 while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
2380 nodeg3 = fGeoManager->GetCurrentNode();
2381 fGeoManager->CdUp();
2382 }
2383 if (!nodegeo) return 0;
2384 if (!nodeg3) return 0;
2385 fGeoManager->cd(common.Data());
2386 fGeoManager->MasterToLocal(fGeoManager->GetCurrentPoint(), &point[0]);
2387 Double_t xyz[3], local[3];
2388 for (Int_t i=0; i<npoints; i++) {
2389 xyz[0] = point[0] - eps[0] + 2*eps[0]*gRandom->Rndm();
2390 xyz[1] = point[1] - eps[1] + 2*eps[1]*gRandom->Rndm();
2391 xyz[2] = point[2] - eps[2] + 2*eps[2]*gRandom->Rndm();
2392 nodeg3->MasterToLocal(&xyz[0], &local[0]);
2393 if (!nodeg3->GetVolume()->Contains(&local[0])) continue;
2394 dist1 = TMath::Sqrt((xyz[0]-point[0])*(xyz[0]-point[0])+
2395 (xyz[1]-point[1])*(xyz[1]-point[1])+(xyz[2]-point[2])*(xyz[2]-point[2]));
2396 if (dist1<dist) {
2397 // save node and closest point
2398 dist = dist1;
2399 node_close = solg3;
2400 // make the random box smaller
2401 eps[0] = TMath::Abs(point[0]-pointg[0]);
2402 eps[1] = TMath::Abs(point[1]-pointg[1]);
2403 eps[2] = TMath::Abs(point[2]-pointg[2]);
2404 }
2405 }
2406 }
2407 if (!node_close) dist = -1;
2408 return node_close;
2409 }
2410
2411 // save current point
2412 memcpy(&point[0], pointg, 3*sizeof(Double_t));
2413 for (Int_t i=0; i<npoints; i++) {
2414 // generate a random point in MARS
2415 fGeoManager->SetCurrentPoint(point[0] - eps[0] + 2*eps[0]*gRandom->Rndm(),
2416 point[1] - eps[1] + 2*eps[1]*gRandom->Rndm(),
2417 point[2] - eps[2] + 2*eps[2]*gRandom->Rndm());
2418 // check if new node is different from the old one
2419 if (node1!=node) {
2420 dist1 = TMath::Sqrt((point[0]-pointg[0])*(point[0]-pointg[0])+
2421 (point[1]-pointg[1])*(point[1]-pointg[1])+(point[2]-pointg[2])*(point[2]-pointg[2]));
2422 if (dist1<dist) {
2423 dist = dist1;
2424 node_close = node1;
2425 memcpy(&closest[0], pointg, 3*sizeof(Double_t));
2426 // make the random box smaller
2427 eps[0] = TMath::Abs(point[0]-pointg[0]);
2428 eps[1] = TMath::Abs(point[1]-pointg[1]);
2429 eps[2] = TMath::Abs(point[2]-pointg[2]);
2430 }
2431 }
2432 }
2433 // restore the original point and path
2434 fGeoManager->FindNode(point[0], point[1], point[2]); // really needed ?
2435 if (!node_close) dist=-1;
2436 return node_close;
2437}
2438
2439////////////////////////////////////////////////////////////////////////////////
2440/// Shoot one ray from start point with direction (dirx,diry,dirz). Fills input array
2441/// with points just after boundary crossings.
2442
2443Double_t *TGeoChecker::ShootRay(Double_t *start, Double_t dirx, Double_t diry, Double_t dirz, Double_t *array, Int_t &nelem, Int_t &dim, Double_t *endpoint) const
2444{
2445// Int_t array_dimension = 3*dim;
2446 nelem = 0;
2447 Int_t istep = 0;
2448 if (!dim) {
2449 printf("empty input array\n");
2450 return array;
2451 }
2452// fGeoManager->CdTop();
2453 const Double_t *point = fGeoManager->GetCurrentPoint();
2454 TGeoNode *endnode;
2455 Bool_t is_entering;
2456 Double_t step, forward;
2457 Double_t dir[3];
2458 dir[0] = dirx;
2459 dir[1] = diry;
2460 dir[2] = dirz;
2461 fGeoManager->InitTrack(start, &dir[0]);
2462 fGeoManager->GetCurrentNode();
2463// printf("Start : (%f,%f,%f)\n", point[0], point[1], point[2]);
2464 fGeoManager->FindNextBoundary();
2465 step = fGeoManager->GetStep();
2466// printf("---next : at step=%f\n", step);
2467 if (step>1E10) return array;
2468 endnode = fGeoManager->Step();
2469 is_entering = fGeoManager->IsEntering();
2470 while (step<1E10) {
2471 if (endpoint) {
2472 forward = dirx*(endpoint[0]-point[0])+diry*(endpoint[1]-point[1])+dirz*(endpoint[2]-point[2]);
2473 if (forward<1E-3) {
2474// printf("exit : Passed start point. nelem=%i\n", nelem);
2475 return array;
2476 }
2477 }
2478 if (is_entering) {
2479 if (nelem>=dim) {
2480 Double_t *temparray = new Double_t[3*(dim+20)];
2481 memcpy(temparray, array, 3*dim*sizeof(Double_t));
2482 delete [] array;
2483 array = temparray;
2484 dim += 20;
2485 }
2486 memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
2487// printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
2488 nelem++;
2489 } else {
2490 if (endnode==0 && step>1E10) {
2491// printf("exit : NULL endnode. nelem=%i\n", nelem);
2492 return array;
2493 }
2494 if (!fGeoManager->IsEntering()) {
2495// if (startnode) printf("stepping %f from (%f, %f, %f) inside %s...\n", step,point[0], point[1], point[2], startnode->GetName());
2496// else printf("stepping %f from (%f, %f, %f) OUTSIDE...\n", step,point[0], point[1], point[2]);
2497 istep = 0;
2498 }
2499 while (!fGeoManager->IsEntering()) {
2500 istep++;
2501 if (istep>1E3) {
2502// Error("ShootRay", "more than 1000 steps. Step was %f", step);
2503 nelem = 0;
2504 return array;
2505 }
2506 fGeoManager->SetStep(1E-5);
2507 endnode = fGeoManager->Step();
2508 }
2509 if (istep>0) printf("%i steps\n", istep);
2510 if (nelem>=dim) {
2511 Double_t *temparray = new Double_t[3*(dim+20)];
2512 memcpy(temparray, array, 3*dim*sizeof(Double_t));
2513 delete [] array;
2514 array = temparray;
2515 dim += 20;
2516 }
2517 memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
2518// if (endnode) printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
2519 nelem++;
2520 }
2521 fGeoManager->FindNextBoundary();
2522 step = fGeoManager->GetStep();
2523// printf("---next at step=%f\n", step);
2524 endnode = fGeoManager->Step();
2525 is_entering = fGeoManager->IsEntering();
2526 }
2527 return array;
2528// printf("exit : INFINITE step. nelem=%i\n", nelem);
2529}
2530
2531////////////////////////////////////////////////////////////////////////////////
2532/// Check time of finding "Where am I" for n points.
2533
2534void TGeoChecker::Test(Int_t npoints, Option_t *option)
2535{
2536 Bool_t recheck = !strcmp(option, "RECHECK");
2537 if (recheck) printf("RECHECK\n");
2538 const TGeoShape *shape = fGeoManager->GetTopVolume()->GetShape();
2539 Double_t dx = ((TGeoBBox*)shape)->GetDX();
2540 Double_t dy = ((TGeoBBox*)shape)->GetDY();
2541 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
2542 Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
2543 Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
2544 Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
2545 Double_t *xyz = new Double_t[3*npoints];
2546 TStopwatch *timer = new TStopwatch();
2547 printf("Random box : %f, %f, %f\n", dx, dy, dz);
2548 timer->Start(kFALSE);
2549 Int_t i;
2550 for (i=0; i<npoints; i++) {
2551 xyz[3*i] = ox-dx+2*dx*gRandom->Rndm();
2552 xyz[3*i+1] = oy-dy+2*dy*gRandom->Rndm();
2553 xyz[3*i+2] = oz-dz+2*dz*gRandom->Rndm();
2554 }
2555 timer->Stop();
2556 printf("Generation time :\n");
2557 timer->Print();
2558 timer->Reset();
2559 TGeoNode *node, *node1;
2560 printf("Start... %i points\n", npoints);
2561 timer->Start(kFALSE);
2562 for (i=0; i<npoints; i++) {
2563 fGeoManager->SetCurrentPoint(xyz+3*i);
2564 if (recheck) fGeoManager->CdTop();
2565 node = fGeoManager->FindNode();
2566 if (recheck) {
2567 node1 = fGeoManager->FindNode();
2568 if (node1 != node) {
2569 printf("Difference for x=%g y=%g z=%g\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
2570 printf(" from top : %s\n", node->GetName());
2571 printf(" redo : %s\n", fGeoManager->GetPath());
2572 }
2573 }
2574 }
2575 timer->Stop();
2576 timer->Print();
2577 delete [] xyz;
2578 delete timer;
2579}
2580
2581////////////////////////////////////////////////////////////////////////////////
2582/// Geometry overlap checker based on sampling.
2583
2584void TGeoChecker::TestOverlaps(const char* path)
2585{
2586 if (fGeoManager->GetTopVolume()!=fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
2587 printf("Checking overlaps for path :\n");
2588 if (!fGeoManager->cd(path)) return;
2589 TGeoNode *checked = fGeoManager->GetCurrentNode();
2590 checked->InspectNode();
2591 // shoot 1E4 points in the shape of the current volume
2592 Int_t npoints = 1000000;
2593 Double_t big = 1E6;
2594 Double_t xmin = big;
2595 Double_t xmax = -big;
2596 Double_t ymin = big;
2597 Double_t ymax = -big;
2598 Double_t zmin = big;
2599 Double_t zmax = -big;
2600 TObjArray *pm = new TObjArray(128);
2601 TPolyMarker3D *marker = 0;
2602 TPolyMarker3D *markthis = new TPolyMarker3D();
2603 markthis->SetMarkerColor(5);
2604 TNtuple *ntpl = new TNtuple("ntpl","random points","x:y:z");
2605 TGeoShape *shape = fGeoManager->GetCurrentNode()->GetVolume()->GetShape();
2606 Double_t *point = new Double_t[3];
2607 Double_t dx = ((TGeoBBox*)shape)->GetDX();
2608 Double_t dy = ((TGeoBBox*)shape)->GetDY();
2609 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
2610 Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
2611 Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
2612 Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
2613 Double_t *xyz = new Double_t[3*npoints];
2614 Int_t i=0;
2615 printf("Generating %i points inside %s\n", npoints, fGeoManager->GetPath());
2616 while (i<npoints) {
2617 point[0] = ox-dx+2*dx*gRandom->Rndm();
2618 point[1] = oy-dy+2*dy*gRandom->Rndm();
2619 point[2] = oz-dz+2*dz*gRandom->Rndm();
2620 if (!shape->Contains(point)) continue;
2621 // convert each point to MARS
2622// printf("local %9.3f %9.3f %9.3f\n", point[0], point[1], point[2]);
2623 fGeoManager->LocalToMaster(point, &xyz[3*i]);
2624// printf("master %9.3f %9.3f %9.3f\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
2625 xmin = TMath::Min(xmin, xyz[3*i]);
2626 xmax = TMath::Max(xmax, xyz[3*i]);
2627 ymin = TMath::Min(ymin, xyz[3*i+1]);
2628 ymax = TMath::Max(ymax, xyz[3*i+1]);
2629 zmin = TMath::Min(zmin, xyz[3*i+2]);
2630 zmax = TMath::Max(zmax, xyz[3*i+2]);
2631 i++;
2632 }
2633 delete [] point;
2634 ntpl->Fill(xmin,ymin,zmin);
2635 ntpl->Fill(xmax,ymin,zmin);
2636 ntpl->Fill(xmin,ymax,zmin);
2637 ntpl->Fill(xmax,ymax,zmin);
2638 ntpl->Fill(xmin,ymin,zmax);
2639 ntpl->Fill(xmax,ymin,zmax);
2640 ntpl->Fill(xmin,ymax,zmax);
2641 ntpl->Fill(xmax,ymax,zmax);
2642 ntpl->Draw("z:y:x");
2643
2644 // shoot the poins in the geometry
2645 TGeoNode *node;
2646 TString cpath;
2647 Int_t ic=0;
2648 TObjArray *overlaps = new TObjArray();
2649 printf("using FindNode...\n");
2650 for (Int_t j=0; j<npoints; j++) {
2651 // always start from top level (testing only)
2652 fGeoManager->CdTop();
2653 fGeoManager->SetCurrentPoint(&xyz[3*j]);
2654 node = fGeoManager->FindNode();
2655 cpath = fGeoManager->GetPath();
2656 if (cpath.Contains(path)) {
2657 markthis->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
2658 continue;
2659 }
2660 // current point is found in an overlapping node
2661 if (!node) ic=128;
2662 else ic = node->GetVolume()->GetLineColor();
2663 if (ic >= 128) ic = 0;
2664 marker = (TPolyMarker3D*)pm->At(ic);
2665 if (!marker) {
2666 marker = new TPolyMarker3D();
2667 marker->SetMarkerColor(ic);
2668 pm->AddAt(marker, ic);
2669 }
2670 // draw the overlapping point
2671 marker->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
2672 if (node) {
2673 if (overlaps->IndexOf(node) < 0) overlaps->Add(node);
2674 }
2675 }
2676 // draw all overlapping points
2677// for (Int_t m=0; m<128; m++) {
2678// marker = (TPolyMarker3D*)pm->At(m);
2679// if (marker) marker->Draw("SAME");
2680// }
2681 markthis->Draw("SAME");
2682 if (gPad) gPad->Update();
2683 // display overlaps
2684 if (overlaps->GetEntriesFast()) {
2685 printf("list of overlapping nodes :\n");
2686 for (i=0; i<overlaps->GetEntriesFast(); i++) {
2687 node = (TGeoNode*)overlaps->At(i);
2688 if (node->IsOverlapping()) printf("%s MANY\n", node->GetName());
2689 else printf("%s ONLY\n", node->GetName());
2690 }
2691 } else printf("No overlaps\n");
2692 delete ntpl;
2693 delete pm;
2694 delete [] xyz;
2695 delete overlaps;
2696}
2697
2698////////////////////////////////////////////////////////////////////////////////
2699/// Estimate weight of top level volume with a precision SIGMA(W)/W
2700/// better than PRECISION. Option can be "v" - verbose (default).
2701
2702Double_t TGeoChecker::Weight(Double_t precision, Option_t *option)
2703{
2704 TList *matlist = fGeoManager->GetListOfMaterials();
2705 Int_t nmat = matlist->GetSize();
2706 if (!nmat) return 0;
2707 Int_t *nin = new Int_t[nmat];
2708 memset(nin, 0, nmat*sizeof(Int_t));
2709 TString opt = option;
2710 opt.ToLower();
2711 Bool_t isverbose = opt.Contains("v");
2712 TGeoBBox *box = (TGeoBBox *)fGeoManager->GetTopVolume()->GetShape();
2713 Double_t dx = box->GetDX();
2714 Double_t dy = box->GetDY();
2715 Double_t dz = box->GetDZ();
2716 Double_t ox = (box->GetOrigin())[0];
2717 Double_t oy = (box->GetOrigin())[1];
2718 Double_t oz = (box->GetOrigin())[2];
2719 Double_t x,y,z;
2720 TGeoNode *node;
2721 TGeoMaterial *mat;
2722 Double_t vbox = 0.000008*dx*dy*dz; // m3
2723 Bool_t end = kFALSE;
2724 Double_t weight=0, sigma, eps, dens;
2725 Double_t eps0=1.;
2726 Int_t indmat;
2727 Int_t igen=0;
2728 Int_t iin = 0;
2729 while (!end) {
2730 x = ox-dx+2*dx*gRandom->Rndm();
2731 y = oy-dy+2*dy*gRandom->Rndm();
2732 z = oz-dz+2*dz*gRandom->Rndm();
2733 node = fGeoManager->FindNode(x,y,z);
2734 igen++;
2735 if (!node) continue;
2736 mat = node->GetVolume()->GetMedium()->GetMaterial();
2737 indmat = mat->GetIndex();
2738 if (indmat<0) continue;
2739 nin[indmat]++;
2740 iin++;
2741 if ((iin%100000)==0 || igen>1E8) {
2742 weight = 0;
2743 sigma = 0;
2744 for (indmat=0; indmat<nmat; indmat++) {
2745 mat = (TGeoMaterial*)matlist->At(indmat);
2746 dens = mat->GetDensity(); // [g/cm3]
2747 if (dens<1E-2) dens=0;
2748 dens *= 1000.; // [kg/m3]
2749 weight += dens*Double_t(nin[indmat]);
2750 sigma += dens*dens*nin[indmat];
2751 }
2752 sigma = TMath::Sqrt(sigma);
2753 eps = sigma/weight;
2754 weight *= vbox/Double_t(igen);
2755 sigma *= vbox/Double_t(igen);
2756 if (eps<precision || igen>1E8) {
2757 if (isverbose) {
2758 printf("=== Weight of %s : %g +/- %g [kg]\n",
2759 fGeoManager->GetTopVolume()->GetName(), weight, sigma);
2760 }
2761 end = kTRUE;
2762 } else {
2763 if (isverbose && eps<0.5*eps0) {
2764 printf("%8dK: %14.7g kg %g %%\n",
2765 igen/1000, weight, eps*100);
2766 eps0 = eps;
2767 }
2768 }
2769 }
2770 }
2771 delete [] nin;
2772 return weight;
2773}
2774////////////////////////////////////////////////////////////////////////////////
2775/// count voxel timing
2776
2777Double_t TGeoChecker::CheckVoxels(TGeoVolume *vol, TGeoVoxelFinder *voxels, Double_t *xyz, Int_t npoints)
2778{
2779 TStopwatch timer;
2780 Double_t time;
2781 TGeoShape *shape = vol->GetShape();
2782 TGeoNode *node;
2783 TGeoMatrix *matrix;
2784 Double_t *point;
2785 Double_t local[3];
2786 Int_t *checklist;
2787 Int_t ncheck;
2788 TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
2789 TGeoStateInfo &td = *nav->GetCache()->GetInfo();
2790 timer.Start();
2791 for (Int_t i=0; i<npoints; i++) {
2792 point = xyz + 3*i;
2793 if (!shape->Contains(point)) continue;
2794 checklist = voxels->GetCheckList(point, ncheck, td);
2795 if (!checklist) continue;
2796 if (!ncheck) continue;
2797 for (Int_t id=0; id<ncheck; id++) {
2798 node = vol->GetNode(checklist[id]);
2799 matrix = node->GetMatrix();
2800 matrix->MasterToLocal(point, &local[0]);
2801 if (node->GetVolume()->GetShape()->Contains(&local[0])) break;
2802 }
2803 }
2804 nav->GetCache()->ReleaseInfo();
2805 time = timer.CpuTime();
2806 return time;
2807}
2808
2809////////////////////////////////////////////////////////////////////////////////
2810/// Returns optimal voxelization type for volume vol.
2811/// - kFALSE - cartesian
2812/// - kTRUE - cylindrical
2813
2814Bool_t TGeoChecker::TestVoxels(TGeoVolume * /*vol*/, Int_t /*npoints*/)
2815{
2816 return kFALSE;
2817}
fFlags
unsigned int fFlags
Definition: DeclareExecutors.h:98
r
ROOT::R::TRInterface & r
Definition: Object.C:4
f
#define f(i)
Definition: RSha256.hxx:104
c
#define c(i)
Definition: RSha256.hxx:101
h
#define h(i)
Definition: RSha256.hxx:106
e
#define e(i)
Definition: RSha256.hxx:103
dummy
static RooMathCoreReg dummy
Definition: RooMathCoreReg.cxx:27
Int_t
int Int_t
Definition: RtypesCore.h:43
Char_t
char Char_t
Definition: RtypesCore.h:31
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:90
Long_t
long Long_t
Definition: RtypesCore.h:52
Double_t
double Double_t
Definition: RtypesCore.h:57
Long64_t
long long Long64_t
Definition: RtypesCore.h:71
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:89
Option_t
const char Option_t
Definition: RtypesCore.h:64
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:361
kRed
@ kRed
Definition: Rtypes.h:64
kBlue
@ kBlue
Definition: Rtypes.h:64
kYellow
@ kYellow
Definition: Rtypes.h:64
kFullCircle
@ kFullCircle
Definition: TAttMarker.h:51
id
XFontStruct * id
Definition: TGX11.cxx:108
name
char name[80]
Definition: TGX11.cxx:109
gGeoIdentity
R__EXTERN TGeoIdentity * gGeoIdentity
Definition: TGeoMatrix.h:478
xmin
float xmin
Definition: THbookFile.cxx:93
nentries
int nentries
Definition: THbookFile.cxx:89
ymin
float ymin
Definition: THbookFile.cxx:93
xmax
float xmax
Definition: THbookFile.cxx:93
ymax
float ymax
Definition: THbookFile.cxx:93
gRandom
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
Form
char * Form(const char *fmt,...)
gStyle
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
gPad
#define gPad
Definition: TVirtualPad.h:287
points
point * points
Definition: X3DBuffer.c:22
TAttFill::SetFillColor
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
TAttLine::GetLineColor
virtual Color_t GetLineColor() const
Return the line color.
Definition: TAttLine.h:33
TAttLine::SetLineWidth
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
TAxis::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:475
TBuffer3D
Generic 3D primitive description class.
Definition: TBuffer3D.h:18
TBuffer3D::NbPols
UInt_t NbPols() const
Definition: TBuffer3D.h:82
TBuffer3D::NbPnts
UInt_t NbPnts() const
Definition: TBuffer3D.h:80
TBuffer3D::NbSegs
UInt_t NbSegs() const
Definition: TBuffer3D.h:81
TBuffer3D::fID
TObject * fID
Definition: TBuffer3D.h:87
TBuffer3D::SetRawSizes
Bool_t SetRawSizes(UInt_t reqPnts, UInt_t reqPntsCapacity, UInt_t reqSegs, UInt_t reqSegsCapacity, UInt_t reqPols, UInt_t reqPolsCapacity)
Set kRaw tessellation section of buffer with supplied sizes.
Definition: TBuffer3D.cxx:359
TBuffer3D::fPnts
Double_t * fPnts
Definition: TBuffer3D.h:112
TCanvas
The Canvas class.
Definition: TCanvas.h:27
TCollection::GetSize
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
TGeoBBox
Box class.
Definition: TGeoBBox.h:18
TGeoChecker
Geometry checking package.
Definition: TGeoChecker.h:38
TGeoChecker::PropagateInGeom
Int_t PropagateInGeom(Double_t *, Double_t *)
Propagate from START along DIR from boundary to boundary until exiting geometry.
Definition: TGeoChecker.cxx:643
TGeoChecker::fNchecks
Int_t fNchecks
Selected node for overlap checking.
Definition: TGeoChecker.h:51
TGeoChecker::CheckShape
void CheckShape(TGeoShape *shape, Int_t testNo, Int_t nsamples, Option_t *option)
Test for shape navigation methods.
Definition: TGeoChecker.cxx:1638
TGeoChecker::fFlags
Bool_t * fFlags
Array of timing per volume.
Definition: TGeoChecker.h:48
TGeoChecker::fTimer
TStopwatch * fTimer
Array of flags per volume.
Definition: TGeoChecker.h:49
TGeoChecker::fVsafe
TGeoVolume * fVsafe
Definition: TGeoChecker.h:42
TGeoChecker::CheckOverlapsBySampling
void CheckOverlapsBySampling(TGeoVolume *vol, Double_t ovlp=0.1, Int_t npoints=1000000) const
Check illegal overlaps for volume VOL within a limit OVLP by sampling npoints inside the volume shape...
Definition: TGeoChecker.cxx:1139
TGeoChecker::TGeoChecker
TGeoChecker()
Default constructor.
Definition: TGeoChecker.cxx:95
TGeoChecker::~TGeoChecker
virtual ~TGeoChecker()
Destructor.
Definition: TGeoChecker.cxx:137
TGeoChecker::ShapeNormal
void ShapeNormal(TGeoShape *shape, Int_t nsamples, Option_t *option)
Check of validity of the normal for a given shape.
Definition: TGeoChecker.cxx:1889
TGeoChecker::CheckOverlaps
void CheckOverlaps(const TGeoVolume *vol, Double_t ovlp=0.1, Option_t *option="")
Check illegal overlaps for volume VOL within a limit OVLP.
Definition: TGeoChecker.cxx:1352
TGeoChecker::TestOverlaps
void TestOverlaps(const char *path)
Geometry overlap checker based on sampling.
Definition: TGeoChecker.cxx:2584
TGeoChecker::MakeCheckOverlap
TGeoOverlap * MakeCheckOverlap(const char *name, TGeoVolume *vol1, TGeoVolume *vol2, TGeoMatrix *mat1, TGeoMatrix *mat2, Bool_t isovlp, Double_t ovlp)
Check if the 2 non-assembly volume candidates overlap/extrude. Returns overlap object.
Definition: TGeoChecker.cxx:971
TGeoChecker::fVal2
Double_t * fVal2
Array of number of crossings per volume.
Definition: TGeoChecker.h:47
TGeoChecker::fGeoManager
TGeoManager * fGeoManager
Definition: TGeoChecker.h:41
TGeoChecker::ShapeDistances
void ShapeDistances(TGeoShape *shape, Int_t nsamples, Option_t *option)
Test TGeoShape::DistFromInside/Outside.
Definition: TGeoChecker.cxx:1664
TGeoChecker::NChecksPerVolume
Int_t NChecksPerVolume(TGeoVolume *vol)
Compute number of overlaps combinations to check per volume.
Definition: TGeoChecker.cxx:1263
TGeoChecker::ShapeSafety
void ShapeSafety(TGeoShape *shape, Int_t nsamples, Option_t *option)
Check of validity of safe distance for a given shape.
Definition: TGeoChecker.cxx:1819
TGeoChecker::Score
void Score(TGeoVolume *, Int_t, Double_t)
Score a hit for VOL.
Definition: TGeoChecker.cxx:680
TGeoChecker::ShootRay
Double_t * ShootRay(Double_t *start, Double_t dirx, Double_t diry, Double_t dirz, Double_t *array, Int_t &nelem, Int_t &dim, Double_t *enpoint=0) const
Shoot one ray from start point with direction (dirx,diry,dirz).
Definition: TGeoChecker.cxx:2443
TGeoChecker::CleanPoints
void CleanPoints(Double_t *points, Int_t &numPoints) const
Number of points on mesh to be checked.
Definition: TGeoChecker.cxx:950
TGeoChecker::SamplePoints
TGeoNode * SamplePoints(Int_t npoints, Double_t &dist, Double_t epsil, const char *g3path)
shoot npoints randomly in a box of 1E-5 around current point.
Definition: TGeoChecker.cxx:2333
TGeoChecker::fSelectedNode
TGeoNode * fSelectedNode
Timer.
Definition: TGeoChecker.h:50
TGeoChecker::Weight
Double_t Weight(Double_t precision=0.01, Option_t *option="v")
Estimate weight of top level volume with a precision SIGMA(W)/W better than PRECISION.
Definition: TGeoChecker.cxx:2702
TGeoChecker::CheckGeometryFull
void CheckGeometryFull(Bool_t checkoverlaps=kTRUE, Bool_t checkcrossings=kTRUE, Int_t nrays=10000, const Double_t *vertex=NULL)
Geometry checking.
Definition: TGeoChecker.cxx:474
TGeoChecker::fNmeshPoints
Int_t fNmeshPoints
Number of checks for current volume.
Definition: TGeoChecker.h:52
TGeoChecker::CheckGeometry
void CheckGeometry(Int_t nrays, Double_t startx, Double_t starty, Double_t startz) const
Shoot nrays with random directions from starting point (startx, starty, startz) in the reference fram...
Definition: TGeoChecker.cxx:774
TGeoChecker::CheckBoundaryReference
virtual void CheckBoundaryReference(Int_t icheck=-1)
Check the boundary errors reference file created by CheckBoundaryErrors method.
Definition: TGeoChecker.cxx:381
TGeoChecker::CheckPoint
void CheckPoint(Double_t x=0, Double_t y=0, Double_t z=0, Option_t *option="")
Draw point (x,y,z) over the picture of the daughters of the volume containing this point.
Definition: TGeoChecker.cxx:1583
TGeoChecker::fBuff2
TBuffer3D * fBuff2
Definition: TGeoChecker.h:44
TGeoChecker::SetNmeshPoints
void SetNmeshPoints(Int_t npoints=1000)
Set number of points to be generated on the shape outline when checking for overlaps.
Definition: TGeoChecker.cxx:695
TGeoChecker::PrintOverlaps
void PrintOverlaps() const
Print the current list of overlaps held by the manager class.
Definition: TGeoChecker.cxx:1570
TGeoChecker::Test
void Test(Int_t npoints, Option_t *option)
Check time of finding "Where am I" for n points.
Definition: TGeoChecker.cxx:2534
TGeoChecker::TimingPerVolume
Double_t TimingPerVolume(TGeoVolume *)
Compute timing per "FindNextBoundary" + "Safety" call.
Definition: TGeoChecker.cxx:707
TGeoChecker::TestVoxels
Bool_t TestVoxels(TGeoVolume *vol, Int_t npoints=1000000)
Returns optimal voxelization type for volume vol.
Definition: TGeoChecker.cxx:2814
TGeoChecker::OpProgress
void OpProgress(const char *opname, Long64_t current, Long64_t size, TStopwatch *watch=0, Bool_t last=kFALSE, Bool_t refresh=kFALSE, const char *msg="")
Print current operation progress.
Definition: TGeoChecker.cxx:147
TGeoChecker::RandomRays
void RandomRays(Int_t nrays, Double_t startx, Double_t starty, Double_t startz, const char *target_vol=0, Bool_t check_norm=kFALSE)
Randomly shoot nrays from point (startx,starty,startz) and plot intersections with surfaces for curre...
Definition: TGeoChecker.cxx:2193
TGeoChecker::fVal1
Double_t * fVal1
Definition: TGeoChecker.h:46
TGeoChecker::fFullCheck
Bool_t fFullCheck
Definition: TGeoChecker.h:45
TGeoChecker::LegoPlot
TH2F * LegoPlot(Int_t ntheta=60, Double_t themin=0., Double_t themax=180., Int_t nphi=90, Double_t phimin=0., Double_t phimax=360., Double_t rmin=0., Double_t rmax=9999999, Option_t *option="")
Generate a lego plot fot the top volume, according to option.
Definition: TGeoChecker.cxx:2043
TGeoChecker::RandomPoints
void RandomPoints(TGeoVolume *vol, Int_t npoints, Option_t *option)
Draw random points in the bounding box of a volume.
Definition: TGeoChecker.cxx:2121
TGeoChecker::CheckBoundaryErrors
virtual void CheckBoundaryErrors(Int_t ntracks=1000000, Double_t radius=-1.)
Check pushes and pulls needed to cross the next boundary with respect to the position given by FindNe...
Definition: TGeoChecker.cxx:238
TGeoChecker::CheckVoxels
Double_t CheckVoxels(TGeoVolume *vol, TGeoVoxelFinder *voxels, Double_t *xyz, Int_t npoints)
count voxel timing
Definition: TGeoChecker.cxx:2777
TGeoChecker::fBuff1
TBuffer3D * fBuff1
Definition: TGeoChecker.h:43
TGeoHMatrix
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition: TGeoMatrix.h:421
TGeoHMatrix::Multiply
void Multiply(const TGeoMatrix *right)
multiply to the right with an other transformation if right is identity matrix, just return
Definition: TGeoMatrix.cxx:2510
TGeoIterator
A geometry iterator.
Definition: TGeoNode.h:245
TGeoIterator::GetCurrentMatrix
const TGeoMatrix * GetCurrentMatrix() const
Returns global matrix for current node.
Definition: TGeoNode.cxx:1083
TGeoIterator::SetTopName
void SetTopName(const char *name)
Set the top name for path.
Definition: TGeoNode.cxx:1149
TGeoIterator::Reset
void Reset(TGeoVolume *top=0)
Resets the iterator for volume TOP.
Definition: TGeoNode.cxx:1138
TGeoIterator::GetLevel
Int_t GetLevel() const
Definition: TGeoNode.h:276
TGeoIterator::GetPath
void GetPath(TString &path) const
Returns the path for the current node.
Definition: TGeoNode.cxx:1110
TGeoIterator::GetNode
TGeoNode * GetNode(Int_t level) const
Returns current node at a given level.
Definition: TGeoNode.cxx:1099
TGeoIterator::Skip
void Skip()
Stop iterating the current branch.
Definition: TGeoNode.cxx:1158
TGeoManager
The manager class for any TGeo geometry.
Definition: TGeoManager.h:43
TGeoManager::GetMother
TGeoNode * GetMother(Int_t up=1) const
Definition: TGeoManager.h:511
TGeoManager::FindNormalFast
Double_t * FindNormalFast()
Computes fast normal to next crossed boundary, assuming that the current point is close enough to the...
Definition: TGeoManager.cxx:2579
TGeoManager::GetListOfUVolumes
TObjArray * GetListOfUVolumes() const
Definition: TGeoManager.h:495
TGeoManager::GetListOfOverlaps
TObjArray * GetListOfOverlaps()
Definition: TGeoManager.h:487
TGeoManager::CdUp
void CdUp()
Go one level up in geometry.
Definition: TGeoManager.cxx:1529
TGeoManager::cd
virtual Bool_t cd(const char *path="")
Browse the tree of nodes starting from fTopNode according to pathname.
Definition: TGeoManager.cxx:1555
TGeoManager::LocalToMaster
void LocalToMaster(const Double_t *local, Double_t *master) const
Definition: TGeoManager.h:541
TGeoManager::RestoreMasterVolume
void RestoreMasterVolume()
Restore the master volume of the geometry.
Definition: TGeoManager.cxx:2874
TGeoManager::FindNextDaughterBoundary
TGeoNode * FindNextDaughterBoundary(Double_t *point, Double_t *dir, Int_t &idaughter, Bool_t compmatrix=kFALSE)
Computes as fStep the distance to next daughter of the current volume.
Definition: TGeoManager.cxx:2546
TGeoManager::GetCurrentNavigator
TGeoNavigator * GetCurrentNavigator() const
Returns current navigator for the calling thread.
Definition: TGeoManager.cxx:739
TGeoManager::GetMasterVolume
TGeoVolume * GetMasterVolume() const
Definition: TGeoManager.h:529
TGeoManager::GetCurrentNode
TGeoNode * GetCurrentNode() const
Definition: TGeoManager.h:517
TGeoManager::SetCurrentDirection
void SetCurrentDirection(Double_t *dir)
Definition: TGeoManager.h:536
TGeoManager::SetVisLevel
void SetVisLevel(Int_t level=3)
set default level down to which visualization is performed
Definition: TGeoManager.cxx:2250
TGeoManager::FindNextBoundary
TGeoNode * FindNextBoundary(Double_t stepmax=TGeoShape::Big(), const char *path="", Bool_t frombdr=kFALSE)
Find distance to next boundary and store it in fStep.
Definition: TGeoManager.cxx:2535
TGeoManager::Step
TGeoNode * Step(Bool_t is_geom=kTRUE, Bool_t cross=kTRUE)
Make a rectilinear step of length fStep from current point (fPoint) on current direction (fDirection)...
Definition: TGeoManager.cxx:3381
TGeoManager::GetVolume
TGeoVolume * GetVolume(const char *name) const
Search for a named volume. All trailing blanks stripped.
Definition: TGeoManager.cxx:2712
TGeoManager::FindNextBoundaryAndStep
TGeoNode * FindNextBoundaryAndStep(Double_t stepmax=TGeoShape::Big(), Bool_t compsafe=kFALSE)
Compute distance to next boundary within STEPMAX.
Definition: TGeoManager.cxx:2516
TGeoManager::SetCurrentPoint
void SetCurrentPoint(Double_t *point)
Definition: TGeoManager.h:533
TGeoManager::FindNormal
Double_t * FindNormal(Bool_t forward=kTRUE)
Computes normal vector to the next surface that will be or was already crossed when propagating on a ...
Definition: TGeoManager.cxx:2590
TGeoManager::FindNode
TGeoNode * FindNode(Bool_t safe_start=kTRUE)
Returns deepest node containing current point.
Definition: TGeoManager.cxx:2562
TGeoManager::GetCurrentPoint
const Double_t * GetCurrentPoint() const
Definition: TGeoManager.h:519
TGeoManager::MakeSphere
TGeoVolume * MakeSphere(const char *name, TGeoMedium *medium, Double_t rmin, Double_t rmax, Double_t themin=0, Double_t themax=180, Double_t phimin=0, Double_t phimax=360)
Make in one step a volume pointing to a sphere shape with given medium.
Definition: TGeoManager.cxx:2939
TGeoManager::IsOutside
Bool_t IsOutside() const
Definition: TGeoManager.h:426
TGeoManager::InitTrack
TGeoNode * InitTrack(const Double_t *point, const Double_t *dir)
Initialize current point and current direction vector (normalized) in MARS.
Definition: TGeoManager.cxx:2632
TGeoManager::GetLevel
Int_t GetLevel() const
Definition: TGeoManager.h:525
TGeoManager::GetStep
Double_t GetStep() const
Definition: TGeoManager.h:406
TGeoManager::GetCurrentMatrix
TGeoHMatrix * GetCurrentMatrix() const
Definition: TGeoManager.h:514
TGeoManager::GetTopNode
TGeoNode * GetTopNode() const
Definition: TGeoManager.h:531
TGeoManager::IsSameLocation
Bool_t IsSameLocation(Double_t x, Double_t y, Double_t z, Bool_t change=kFALSE)
Checks if point (x,y,z) is still in the current node.
Definition: TGeoManager.cxx:2598
TGeoManager::GetPath
const char * GetPath() const
Get path to the current node in the form /node0/node1/...
Definition: TGeoManager.cxx:2657
TGeoManager::AddOverlap
Int_t AddOverlap(const TNamed *ovlp)
Add an illegal overlap/extrusion to the list.
Definition: TGeoManager.cxx:562
TGeoManager::SetVerboseLevel
static void SetVerboseLevel(Int_t vl)
Return current verbosity level (static function).
Definition: TGeoManager.cxx:3797
TGeoManager::SetStep
void SetStep(Double_t step)
Definition: TGeoManager.h:420
TGeoManager::GetCurrentVolume
TGeoVolume * GetCurrentVolume() const
Definition: TGeoManager.h:521
TGeoManager::GetVerboseLevel
static Int_t GetVerboseLevel()
Set verbosity level (static function).
Definition: TGeoManager.cxx:3789
TGeoManager::CdTop
void CdTop()
Make top level node the current node.
Definition: TGeoManager.cxx:1520
TGeoManager::Safety
Double_t Safety(Bool_t inside=kFALSE)
Compute safe distance from the current point.
Definition: TGeoManager.cxx:2117
TGeoManager::MasterToLocal
void MasterToLocal(const Double_t *master, Double_t *local) const
Definition: TGeoManager.h:544
TGeoManager::ResetState
void ResetState()
Reset current state flags.
Definition: TGeoManager.cxx:2554
TGeoManager::CdDown
void CdDown(Int_t index)
Make a daughter of current node current.
Definition: TGeoManager.cxx:1538
TGeoManager::GetListOfMaterials
TList * GetListOfMaterials() const
Definition: TGeoManager.h:489
TGeoManager::GetTopVolume
TGeoVolume * GetTopVolume() const
Definition: TGeoManager.h:530
TGeoManager::SetTopVisible
void SetTopVisible(Bool_t vis=kTRUE)
make top volume visible on screen
Definition: TGeoManager.cxx:2195
TGeoManager::CheckOverlaps
void CheckOverlaps(Double_t ovlp=0.1, Option_t *option="")
Check all geometry for illegal overlaps within a limit OVLP.
Definition: TGeoManager.cxx:3602
TGeoManager::IsEntering
Bool_t IsEntering() const
Definition: TGeoManager.h:422
TGeoMaterial
Base class describing materials.
Definition: TGeoMaterial.h:31
TGeoMaterial::GetIndex
Int_t GetIndex()
Retrieve material index in the list of materials.
Definition: TGeoMaterial.cxx:598
TGeoMaterial::GetRadLen
virtual Double_t GetRadLen() const
Definition: TGeoMaterial.h:110
TGeoMaterial::GetDensity
virtual Double_t GetDensity() const
Definition: TGeoMaterial.h:103
TGeoMatrix
Geometrical transformation package.
Definition: TGeoMatrix.h:41
TGeoMatrix::LocalToMasterVect
virtual void LocalToMasterVect(const Double_t *local, Double_t *master) const
convert a vector by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition: TGeoMatrix.cxx:363
TGeoMatrix::MasterToLocal
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix
Definition: TGeoMatrix.cxx:406
TGeoMatrix::LocalToMaster
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition: TGeoMatrix.cxx:339
TGeoMedium::GetMaterial
TGeoMaterial * GetMaterial() const
Definition: TGeoMedium.h:52
TGeoNavigator
Class providing navigation API for TGeo geometries.
Definition: TGeoNavigator.h:34
TGeoNavigator::GetCache
TGeoNodeCache * GetCache() const
Definition: TGeoNavigator.h:193
TGeoNodeCache::GetInfo
TGeoStateInfo * GetInfo()
Get next state info pointer.
Definition: TGeoCache.cxx:319
TGeoNodeCache::ReleaseInfo
void ReleaseInfo()
Release last used state info pointer.
Definition: TGeoCache.cxx:336
TGeoNode
A node represent a volume positioned inside another.They store links to both volumes and to the TGeoM...
Definition: TGeoNode.h:41
TGeoNode::IsOverlapping
Bool_t IsOverlapping() const
Definition: TGeoNode.h:105
TGeoNode::IsOnScreen
Bool_t IsOnScreen() const
check if this node is drawn. Assumes that this node is current
Definition: TGeoNode.cxx:273
TGeoNode::GetVolume
TGeoVolume * GetVolume() const
Definition: TGeoNode.h:97
TGeoNode::SetOverlaps
void SetOverlaps(Int_t *ovlp, Int_t novlp)
set the list of overlaps for this node (ovlp must be created with operator new)
Definition: TGeoNode.cxx:653
TGeoNode::GetMatrix
virtual TGeoMatrix * GetMatrix() const =0
TGeoNode::GetColour
Int_t GetColour() const
Definition: TGeoNode.h:88
TGeoNode::GetOverlaps
Int_t * GetOverlaps(Int_t &novlp) const
Definition: TGeoNode.h:96
TGeoNode::MasterToLocal
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
Convert the point coordinates from mother reference to local reference system.
Definition: TGeoNode.cxx:518
TGeoNode::InspectNode
void InspectNode() const
Inspect this node.
Definition: TGeoNode.cxx:282
TGeoOverlap
Base class describing geometry overlaps.
Definition: TGeoOverlap.h:41
TGeoOverlap::GetPolyMarker
TPolyMarker3D * GetPolyMarker() const
Definition: TGeoOverlap.h:72
TGeoOverlap::SetNextPoint
void SetNextPoint(Double_t x, Double_t y, Double_t z)
Set next overlapping point.
Definition: TGeoOverlap.cxx:172
TGeoOverlap::SetOverlap
void SetOverlap(Double_t ovlp)
Definition: TGeoOverlap.h:94
TGeoOverlap::PrintInfo
virtual void PrintInfo() const
Print some info.
Definition: TGeoOverlap.cxx:164
TGeoOverlap::GetOverlap
Double_t GetOverlap() const
Definition: TGeoOverlap.h:77
TGeoShape
Base abstract class for all shapes.
Definition: TGeoShape.h:26
TGeoShape::ComputeNormal
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)=0
TGeoShape::Big
static Double_t Big()
Definition: TGeoShape.h:88
TGeoShape::GetMeshNumbers
virtual void GetMeshNumbers(Int_t &, Int_t &, Int_t &) const
Definition: TGeoShape.h:125
TGeoShape::GetPointsOnSegments
virtual Bool_t GetPointsOnSegments(Int_t npoints, Double_t *array) const =0
TGeoShape::SetTransform
static void SetTransform(TGeoMatrix *matrix)
Set current transformation matrix that applies to shape.
Definition: TGeoShape.cxx:544
TGeoShape::Safety
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const =0
TGeoShape::GetName
virtual const char * GetName() const
Get the shape name.
Definition: TGeoShape.cxx:248
TGeoShape::DistFromOutside
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const =0
TGeoShape::DistFromInside
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const =0
TGeoShape::Capacity
virtual Double_t Capacity() const =0
TGeoShape::Contains
virtual Bool_t Contains(const Double_t *point) const =0
TGeoShape::Tolerance
static Double_t Tolerance()
Definition: TGeoShape.h:91
TGeoShape::SetPoints
virtual void SetPoints(Double_t *points) const =0
TGeoShape::Draw
virtual void Draw(Option_t *option="")
Draw this shape.
Definition: TGeoShape.cxx:721
TGeoTranslation
Class describing translations.
Definition: TGeoMatrix.h:122
TGeoVolume
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition: TGeoVolume.h:47
TGeoVolume::GetMedium
TGeoMedium * GetMedium() const
Definition: TGeoVolume.h:171
TGeoVolume::Contains
Bool_t Contains(const Double_t *point) const
Definition: TGeoVolume.h:107
TGeoVolume::GetMaterial
TGeoMaterial * GetMaterial() const
Definition: TGeoVolume.h:170
TGeoVolume::GetNdaughters
Int_t GetNdaughters() const
Definition: TGeoVolume.h:347
TGeoVolume::SetLineStyle
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TGeoVolume.cxx:2043
TGeoVolume::GetNodes
TObjArray * GetNodes()
Definition: TGeoVolume.h:165
TGeoVolume::VisibleDaughters
void VisibleDaughters(Bool_t vis=kTRUE)
set visibility for daughters
Definition: TGeoVolume.cxx:2291
TGeoVolume::FindOverlaps
void FindOverlaps() const
loop all nodes marked as overlaps and find overlapping brothers
Definition: TGeoVolume.cxx:2089
TGeoVolume::GetNode
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
Definition: TGeoVolume.cxx:2059
TGeoVolume::SetVisibility
virtual void SetVisibility(Bool_t vis=kTRUE)
set visibility of this volume
Definition: TGeoVolume.cxx:2197
TGeoVolume::GetIndex
Int_t GetIndex(const TGeoNode *node) const
get index number for a given daughter
Definition: TGeoVolume.cxx:1632
TGeoVolume::GetFinder
TGeoPatternFinder * GetFinder() const
Definition: TGeoVolume.h:173
TGeoVolume::GetVoxels
TGeoVoxelFinder * GetVoxels() const
Getter for optimization structure.
Definition: TGeoVolume.cxx:1686
TGeoVolume::GetNumber
Int_t GetNumber() const
Definition: TGeoVolume.h:180
TGeoVolume::GetShape
TGeoShape * GetShape() const
Definition: TGeoVolume.h:186
TGeoVolume::Draw
virtual void Draw(Option_t *option="")
draw top volume according to option
Definition: TGeoVolume.cxx:1145
TGeoVolume::GetCurrentNodeIndex
virtual Int_t GetCurrentNodeIndex() const
Definition: TGeoVolume.h:163
TGeoVolume::DrawOnly
virtual void DrawOnly(Option_t *option="")
draw only this volume
Definition: TGeoVolume.cxx:1162
TGeoVolume::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TGeoVolume.cxx:2035
TGeoVolume::IsAssembly
virtual Bool_t IsAssembly() const
Returns true if the volume is an assembly or a scaled assembly.
Definition: TGeoVolume.cxx:1704
TGeoVolume::AddNode
virtual void AddNode(TGeoVolume *vol, Int_t copy_no, TGeoMatrix *mat=0, Option_t *option="")
Add a TGeoNode to the list of nodes.
Definition: TGeoVolume.cxx:931
TGeoVolume::IsVisible
virtual Bool_t IsVisible() const
Definition: TGeoVolume.h:151
TGeoVolume::InspectShape
void InspectShape() const
Definition: TGeoVolume.h:191
TGeoVoxelFinder
Finder class handling voxels.
Definition: TGeoVoxelFinder.h:21
TGeoVoxelFinder::GetCheckList
virtual Int_t * GetCheckList(const Double_t *point, Int_t &nelem, TGeoStateInfo &td)
get the list of daughter indices for which point is inside their bbox
Definition: TGeoVoxelFinder.cxx:1206
TGeoVoxelFinder::Voxelize
virtual void Voxelize(Option_t *option="")
Voxelize attached volume according to option If the volume is an assembly, make sure the bbox is comp...
Definition: TGeoVoxelFinder.cxx:2222
TGeoVoxelFinder::FindOverlaps
virtual void FindOverlaps(Int_t inode) const
create the list of nodes for which the bboxes overlap with inode's bbox
Definition: TGeoVoxelFinder.cxx:297
TGeoVoxelFinder::NeedRebuild
Bool_t NeedRebuild() const
Definition: TGeoVoxelFinder.h:105
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
TH1::LabelsOption
virtual void LabelsOption(Option_t *option="h", Option_t *axis="X")
Set option(s) to draw axis with labels.
Definition: TH1.cxx:5222
TH1::AddDirectory
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1226
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TH1::kAllAxes
@ kAllAxes
Definition: TH1.h:73
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:317
TH1::SetCanExtend
virtual UInt_t SetCanExtend(UInt_t extendBitMask)
Make the histogram axes extendable / not extendable according to the bit mask returns the previous bi...
Definition: TH1.cxx:6291
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
TH1::LabelsDeflate
virtual void LabelsDeflate(Option_t *axis="X")
Reduce the number of bins for the axis passed in the option to the number of bins having a label.
Definition: TH1.cxx:5091
TH1::SetStats
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8446
TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
TH2::Fill
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
TList
A doubly linked list.
Definition: TList.h:44
TList::At
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:356
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TNtuple
A simple TTree restricted to a list of float variables only.
Definition: TNtuple.h:28
TNtuple::Fill
virtual Int_t Fill()
Fill a Ntuple with current values in fArgs.
Definition: TNtuple.cxx:170
TObjArray
An array of TObjects.
Definition: TObjArray.h:37
TObjArray::IndexOf
Int_t IndexOf(const TObject *obj) const
Definition: TObjArray.cxx:605
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
TObjArray::Add
void Add(TObject *obj)
Definition: TObjArray.h:74
TObjArray::GetEntries
Int_t GetEntries() const
Return the number of objects in array (i.e.
Definition: TObjArray.cxx:523
TObjArray::AddAt
virtual void AddAt(TObject *obj, Int_t idx)
Add object at position ids.
Definition: TObjArray.cxx:254
TObjArray::RemoveAt
virtual TObject * RemoveAt(Int_t idx)
Remove object at index idx.
Definition: TObjArray.cxx:694
TObjArray::At
TObject * At(Int_t idx) const
Definition: TObjArray.h:166
TObject
Mother of all ROOT objects.
Definition: TObject.h:37
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:877
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
TObject::Draw
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
TObject::Info
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:865
TPolyLine3D
A 3-dimensional polyline.
Definition: TPolyLine3D.h:32
TPolyLine3D::SetPoint
virtual void SetPoint(Int_t point, Double_t x, Double_t y, Double_t z)
Set point n to x, y, z.
Definition: TPolyLine3D.cxx:626
TPolyMarker3D
A 3D polymarker.
Definition: TPolyMarker3D.h:33
TPolyMarker3D::GetN
virtual Int_t GetN() const
Definition: TPolyMarker3D.h:58
TPolyMarker3D::SetNextPoint
virtual Int_t SetNextPoint(Double_t x, Double_t y, Double_t z)
Set point following LastPoint to x, y, z.
Definition: TPolyMarker3D.cxx:546
TPolyMarker3D::Draw
virtual void Draw(Option_t *option="")
Draws 3-D polymarker with its current attributes.
Definition: TPolyMarker3D.cxx:265
TRandom::Uniform
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:635
TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:28
TStopwatch::RealTime
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
TStopwatch::CpuTime
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
TStopwatch::Continue
void Continue()
Resume a stopped stopwatch.
Definition: TStopwatch.cxx:93
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
TStopwatch::Reset
void Reset()
Definition: TStopwatch.h:52
TStopwatch::Print
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
TString
Basic string class.
Definition: TString.h:131
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
TString::Strip
TSubString Strip(EStripType s=kTrailing, char c=' ') const
Return a substring of self stripped at beginning and/or end.
Definition: TString.cxx:1106
TString::Data
const char * Data() const
Definition: TString.h:364
TString::kLeading
@ kLeading
Definition: TString.h:262
TString::Remove
TString & Remove(Ssiz_t pos)
Definition: TString.h:668
TString::Format
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition: TString.cxx:2311
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TString::Index
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:634
TStyle::SetOptStat
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1590
TTree
A TTree represents a columnar dataset.
Definition: TTree.h:78
TTree::Fill
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4524
TTree::SetBranchAddress
virtual Int_t SetBranchAddress(const char *bname, void *add, TBranch **ptr=0)
Change branch address, dealing with clone trees properly.
Definition: TTree.cxx:8237
TTree::GetEntries
virtual Long64_t GetEntries() const
Definition: TTree.h:457
TTree::Branch
TBranch * Branch(const char *name, T *obj, Int_t bufsize=32000, Int_t splitlevel=99)
Add a new branch, and infer the data type from the type of obj being passed.
Definition: TTree.h:348
TTree::GetEntry
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition: TTree.cxx:5542
TTree::Draw
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition: TTree.h:426
TTree::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TTree.cxx:9595
pt
TPaveText * pt
Definition: entrylist_figure1.C:7
line
TLine * line
Definition: entrylistblock_figure1.C:235
box
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
c1
return c1
Definition: legend1.C:41
y
Double_t y[n]
Definition: legend1.C:17
x
Double_t x[n]
Definition: legend1.C:17
h1
TH1F * h1
Definition: legend1.C:5
c2
return c2
Definition: legend2.C:14
c3
return c3
Definition: legend3.C:15
ROOT::Math::gv_detail::dist
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
TGeant4Unit::rad
static constexpr double rad
Definition: TGeant4SystemOfUnits.h:142
TMVA::DNN::forward
void forward(const LAYERDATA &prevLayerData, LAYERDATA &currLayerData)
apply the weights (and functions) in forward direction of the DNN
Definition: NeuralNet.icc:546
TMath::ACos
Double_t ACos(Double_t)
Definition: TMath.h:658
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TMath::E
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:97
TMath::Log
Double_t Log(Double_t x)
Definition: TMath.h:750
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
TMath::Cos
Double_t Cos(Double_t)
Definition: TMath.h:631
TMath::Pi
constexpr Double_t Pi()
Definition: TMath.h:38
TMath::Sin
Double_t Sin(Double_t)
Definition: TMath.h:627
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TMath::TwoPi
constexpr Double_t TwoPi()
Definition: TMath.h:45
TGeoStateInfo
Statefull info for the current geometry level.
Definition: TGeoStateInfo.h:21
m
auto * m
Definition: textangle.C:8
dnext
#define dnext(otri1, otri2)
Definition: triangle.c:986
lnext
#define lnext(otri1, otri2)
Definition: triangle.c:942
vertex
REAL * vertex
Definition: triangle.c:512