Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
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_painter
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 <cstdlib>
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 for (Int_t i=0; i<1000000; i++) {
725 lpt[0] = ox-dx+2*dx*gRandom->Rndm();
726 lpt[1] = oy-dy+2*dy*gRandom->Rndm();
727 lpt[2] = oz-dz+2*dz*gRandom->Rndm();
728 fGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
729 fGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
730 phi = 2*TMath::Pi()*gRandom->Rndm();
731 theta= TMath::ACos(1.-2.*gRandom->Rndm());
732 ldir[0]=TMath::Sin(theta)*TMath::Cos(phi);
733 ldir[1]=TMath::Sin(theta)*TMath::Sin(phi);
734 ldir[2]=TMath::Cos(theta);
735 fGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir,dir);
736 fGeoManager->SetCurrentDirection(dir);
737 fGeoManager->SetStep(pstep);
738 fGeoManager->ResetState();
739 // dist = TGeoShape::Big();
740 if (!vol->IsAssembly()) {
741 Bool_t inside = vol->Contains(lpt);
742 if (!inside) {
743 // dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep);
744 // if (dist>=pstep) continue;
745 } else {
746 vol->GetShape()->DistFromInside(lpt,ldir,3,pstep);
747 }
748
749 if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside);
750 }
751 if (vol->GetNdaughters()) {
752 fGeoManager->Safety();
753 fGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE);
754 }
755 }
756 fTimer->Stop();
757 Double_t time_per_track = fTimer->CpuTime();
758 Int_t uid = vol->GetNumber();
759 Int_t ncrossings = (Int_t)fVal1[uid];
760 if (!vol->GetNdaughters())
761 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);
762 else
763 printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(), time_per_track, vol->GetNdaughters(), ncrossings);
764 return time_per_track;
765}
766
767////////////////////////////////////////////////////////////////////////////////
768/// Shoot nrays with random directions from starting point (startx, starty, startz)
769/// in the reference frame of this volume. Track each ray until exiting geometry, then
770/// shoot backwards from exiting point and compare boundary crossing points.
771
772void TGeoChecker::CheckGeometry(Int_t nrays, Double_t startx, Double_t starty, Double_t startz) const
773{
774 Int_t i, j;
775 Double_t start[3], end[3];
776 Double_t dir[3];
777 Double_t dummy[3];
778 Double_t eps = 0.;
779 Double_t *array1 = new Double_t[3*1000];
780 Double_t *array2 = new Double_t[3*1000];
781 TObjArray *pma = new TObjArray();
782 TPolyMarker3D *pm;
783 pm = new TPolyMarker3D();
784 pm->SetMarkerColor(2); // error > eps
785 pm->SetMarkerStyle(8);
786 pm->SetMarkerSize(0.4);
787 pma->AddAt(pm, 0);
788 pm = new TPolyMarker3D();
789 pm->SetMarkerColor(4); // point not found
790 pm->SetMarkerStyle(8);
791 pm->SetMarkerSize(0.4);
792 pma->AddAt(pm, 1);
793 pm = new TPolyMarker3D();
794 pm->SetMarkerColor(6); // extra point back
795 pm->SetMarkerStyle(8);
796 pm->SetMarkerSize(0.4);
797 pma->AddAt(pm, 2);
798 Int_t nelem1, nelem2;
799 Int_t dim1=1000, dim2=1000;
800 if ((startx==0) && (starty==0) && (startz==0)) eps=1E-3;
801 start[0] = startx+eps;
802 start[1] = starty+eps;
803 start[2] = startz+eps;
804 Int_t n10=nrays/10;
805 Double_t theta,phi;
806 Double_t dw, dwmin, dx, dy, dz;
807 Int_t ist1, ist2, ifound;
808 for (i=0; i<nrays; i++) {
809 if (n10) {
810 if ((i%n10) == 0) printf("%i percent\n", Int_t(100*i/nrays));
811 }
812 phi = 2*TMath::Pi()*gRandom->Rndm();
813 theta= TMath::ACos(1.-2.*gRandom->Rndm());
814 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
815 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
816 dir[2]=TMath::Cos(theta);
817 // shoot direct ray
818 nelem1=nelem2=0;
819// printf("DIRECT %i\n", i);
820 array1 = ShootRay(&start[0], dir[0], dir[1], dir[2], array1, nelem1, dim1);
821 if (!nelem1) continue;
822// for (j=0; j<nelem1; j++) printf("%i : %f %f %f\n", j, array1[3*j], array1[3*j+1], array1[3*j+2]);
823 memcpy(&end[0], &array1[3*(nelem1-1)], 3*sizeof(Double_t));
824 // shoot ray backwards
825// printf("BACK %i\n", i);
826 array2 = ShootRay(&end[0], -dir[0], -dir[1], -dir[2], array2, nelem2, dim2, &start[0]);
827 if (!nelem2) {
828 printf("#### NOTHING BACK ###########################\n");
829 for (j=0; j<nelem1; j++) {
830 pm = (TPolyMarker3D*)pma->At(0);
831 pm->SetNextPoint(array1[3*j], array1[3*j+1], array1[3*j+2]);
832 }
833 continue;
834 }
835// printf("BACKWARDS\n");
836 Int_t k=nelem2>>1;
837 for (j=0; j<k; j++) {
838 memcpy(&dummy[0], &array2[3*j], 3*sizeof(Double_t));
839 memcpy(&array2[3*j], &array2[3*(nelem2-1-j)], 3*sizeof(Double_t));
840 memcpy(&array2[3*(nelem2-1-j)], &dummy[0], 3*sizeof(Double_t));
841 }
842// for (j=0; j<nelem2; j++) printf("%i : %f ,%f ,%f \n", j, array2[3*j], array2[3*j+1], array2[3*j+2]);
843 if (nelem1!=nelem2) printf("### DIFFERENT SIZES : nelem1=%i nelem2=%i ##########\n", nelem1, nelem2);
844 ist1 = ist2 = 0;
845 // check first match
846
847 dx = array1[3*ist1]-array2[3*ist2];
848 dy = array1[3*ist1+1]-array2[3*ist2+1];
849 dz = array1[3*ist1+2]-array2[3*ist2+2];
850 dw = dx*dir[0]+dy*dir[1]+dz*dir[2];
851 fGeoManager->SetCurrentPoint(&array1[3*ist1]);
852 fGeoManager->FindNode();
853// printf("%i : %s (%g, %g, %g)\n", ist1, fGeoManager->GetPath(),
854// array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2]);
855 if (TMath::Abs(dw)<1E-4) {
856// printf(" matching %i (%g, %g, %g)\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
857 ist2++;
858 } else {
859 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);
860 pm = (TPolyMarker3D*)pma->At(0);
861 pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
862 if (dw<0) {
863 // first boundary missed on way back
864 } else {
865 // first boundary different on way back
866 ist2++;
867 }
868 }
869
870 while ((ist1<nelem1-1) && (ist2<nelem2)) {
871 fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
872 fGeoManager->FindNode();
873// printf("%i : %s (%g, %g, %g)\n", ist1+1, fGeoManager->GetPath(),
874// array1[3*ist1+3], array1[3*ist1+4], array1[3*ist1+5]);
875
876 dx = array1[3*ist1+3]-array1[3*ist1];
877 dy = array1[3*ist1+4]-array1[3*ist1+1];
878 dz = array1[3*ist1+5]-array1[3*ist1+2];
879 // distance to next point
880 dwmin = dx+dir[0]+dy*dir[1]+dz*dir[2];
881 while (ist2<nelem2) {
882 ifound = 0;
883 dx = array2[3*ist2]-array1[3*ist1];
884 dy = array2[3*ist2+1]-array1[3*ist1+1];
885 dz = array2[3*ist2+2]-array1[3*ist1+2];
886 dw = dx+dir[0]+dy*dir[1]+dz*dir[2];
887 if (TMath::Abs(dw-dwmin)<1E-4) {
888 ist1++;
889 ist2++;
890 break;
891 }
892 if (dw<dwmin) {
893 // point found on way back. Check if close enough to ist1+1
894 ifound++;
895 dw = dwmin-dw;
896 if (dw<1E-4) {
897 // point is matching ist1+1
898// printf(" matching %i (%g, %g, %g) DCLOSE=%g\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2], dw);
899 ist2++;
900 ist1++;
901 break;
902 } else {
903 // extra boundary found on way back
904 fGeoManager->SetCurrentPoint(&array2[3*ist2]);
905 fGeoManager->FindNode();
906 pm = (TPolyMarker3D*)pma->At(2);
907 pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
908 printf("### EXTRA BOUNDARY %i : %s found at DCLOSE=%f\n", ist2, fGeoManager->GetPath(), dw);
909 ist2++;
910 continue;
911 }
912 } else {
913 if (!ifound) {
914 // point ist1+1 not found on way back
915 fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
916 fGeoManager->FindNode();
917 pm = (TPolyMarker3D*)pma->At(1);
918 pm->SetNextPoint(array2[3*ist1+3], array2[3*ist1+4], array2[3*ist1+5]);
919 printf("### BOUNDARY MISSED BACK #########################\n");
920 ist1++;
921 break;
922 } else {
923 ist1++;
924 break;
925 }
926 }
927 }
928 }
929 }
930 pm = (TPolyMarker3D*)pma->At(0);
931 pm->Draw("SAME");
932 pm = (TPolyMarker3D*)pma->At(1);
933 pm->Draw("SAME");
934 pm = (TPolyMarker3D*)pma->At(2);
935 pm->Draw("SAME");
936 if (gPad) {
937 gPad->Modified();
938 gPad->Update();
939 }
940 delete [] array1;
941 delete [] array2;
942 delete pma; // markers are drawn on the pad
943}
944
945////////////////////////////////////////////////////////////////////////////////
946/// Clean-up the mesh of pcon/pgon from useless points
947
948void TGeoChecker::CleanPoints(Double_t *points, Int_t &numPoints) const
949{
950 Int_t ipoint = 0;
951 Int_t j, k=0;
952 Double_t rsq;
953 for (Int_t i=0; i<numPoints; i++) {
954 j = 3*i;
955 rsq = points[j]*points[j]+points[j+1]*points[j+1];
956 if (rsq < 1.e-10) continue;
957 points[k] = points[j];
958 points[k+1] = points[j+1];
959 points[k+2] = points[j+2];
960 ipoint++;
961 k = 3*ipoint;
962 }
963 numPoints = ipoint;
964}
965
966////////////////////////////////////////////////////////////////////////////////
967/// Check if the 2 non-assembly volume candidates overlap/extrude. Returns overlap object.
968
969TGeoOverlap *TGeoChecker::MakeCheckOverlap(const char *name, TGeoVolume *vol1, TGeoVolume *vol2, TGeoMatrix *mat1, TGeoMatrix *mat2, Bool_t isovlp, Double_t ovlp)
970{
971 TGeoOverlap *nodeovlp = 0;
972 Int_t numPoints1 = fBuff1->NbPnts();
973 Int_t numSegs1 = fBuff1->NbSegs();
974 Int_t numPols1 = fBuff1->NbPols();
975 Int_t numPoints2 = fBuff2->NbPnts();
976 Int_t numSegs2 = fBuff2->NbSegs();
977 Int_t numPols2 = fBuff2->NbPols();
978 Int_t ip;
979 Bool_t extrude, isextrusion, isoverlapping;
980 Double_t *points1 = fBuff1->fPnts;
981 Double_t *points2 = fBuff2->fPnts;
982 Double_t local[3], local1[3];
983 Double_t point[3];
984 Double_t safety = TGeoShape::Big();
985 Double_t tolerance = TGeoShape::Tolerance();
986 if (vol1->IsAssembly() || vol2->IsAssembly()) return nodeovlp;
987 TGeoShape *shape1 = vol1->GetShape();
988 TGeoShape *shape2 = vol2->GetShape();
989 OpProgress("refresh", 0,0,NULL,kFALSE,kTRUE);
990 shape1->GetMeshNumbers(numPoints1, numSegs1, numPols1);
991 if (fBuff1->fID != (TObject*)shape1) {
992 // Fill first buffer.
993 fBuff1->SetRawSizes(TMath::Max(numPoints1,fNmeshPoints), 3*TMath::Max(numPoints1,fNmeshPoints), 0, 0, 0, 0);
994 points1 = fBuff1->fPnts;
995 if (shape1->GetPointsOnSegments(fNmeshPoints, points1)) {
996 numPoints1 = fNmeshPoints;
997 } else {
998 shape1->SetPoints(points1);
999 }
1000// if (shape1->InheritsFrom(TGeoPcon::Class())) {
1001// CleanPoints(points1, numPoints1);
1002// fBuff1->SetRawSizes(numPoints1, 3*numPoints1,0, 0, 0, 0);
1003// }
1004 fBuff1->fID = shape1;
1005 }
1006 shape2->GetMeshNumbers(numPoints2, numSegs2, numPols2);
1007 if (fBuff2->fID != (TObject*)shape2) {
1008 // Fill second buffer.
1009 fBuff2->SetRawSizes(TMath::Max(numPoints2,fNmeshPoints), 3*TMath::Max(numPoints2,fNmeshPoints), 0, 0, 0,0);
1010 points2 = fBuff2->fPnts;
1011 if (shape2->GetPointsOnSegments(fNmeshPoints, points2)) {
1012 numPoints2 = fNmeshPoints;
1013 } else {
1014 shape2->SetPoints(points2);
1015 }
1016// if (shape2->InheritsFrom(TGeoPcon::Class())) {
1017// CleanPoints(points2, numPoints2);
1018// fBuff1->SetRawSizes(numPoints2, 3*numPoints2,0,0,0,0);
1019// }
1020 fBuff2->fID = shape2;
1021 }
1022
1023 if (!isovlp) {
1024 // Extrusion case. Test vol2 extrude vol1.
1025 isextrusion=kFALSE;
1026 // loop all points of the daughter
1027 for (ip=0; ip<numPoints2; ip++) {
1028 memcpy(local, &points2[3*ip], 3*sizeof(Double_t));
1029 if (TMath::Abs(local[0])<tolerance && TMath::Abs(local[1])<tolerance) continue;
1030 mat2->LocalToMaster(local, point);
1031 mat1->MasterToLocal(point, local);
1032 extrude = !shape1->Contains(local);
1033 if (extrude) {
1034 safety = shape1->Safety(local, kFALSE);
1035 if (safety<ovlp) extrude=kFALSE;
1036 }
1037 if (extrude) {
1038 if (!isextrusion) {
1039 isextrusion = kTRUE;
1040 nodeovlp = new TGeoOverlap(name, vol1, vol2, mat1,mat2,kFALSE,safety);
1041 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1042 fGeoManager->AddOverlap(nodeovlp);
1043 } else {
1044 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1045 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1046 }
1047 }
1048 }
1049 // loop all points of the mother
1050 for (ip=0; ip<numPoints1; ip++) {
1051 memcpy(local, &points1[3*ip], 3*sizeof(Double_t));
1052 if (local[0]<1e-10 && local[1]<1e-10) continue;
1053 mat1->LocalToMaster(local, point);
1054 mat2->MasterToLocal(point, local1);
1055 extrude = shape2->Contains(local1);
1056 if (extrude) {
1057 // skip points on mother mesh that have no neighborhood outside mother
1058 safety = shape1->Safety(local,kTRUE);
1059 if (safety>1E-6) {
1060 extrude = kFALSE;
1061 } else {
1062 safety = shape2->Safety(local1,kTRUE);
1063 if (safety<ovlp) extrude=kFALSE;
1064 }
1065 }
1066 if (extrude) {
1067 if (!isextrusion) {
1068 isextrusion = kTRUE;
1069 nodeovlp = new TGeoOverlap(name, vol1,vol2,mat1,mat2,kFALSE,safety);
1070 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1071 fGeoManager->AddOverlap(nodeovlp);
1072 } else {
1073 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1074 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1075 }
1076 }
1077 }
1078 return nodeovlp;
1079 }
1080 // Check overlap
1081 isoverlapping = kFALSE;
1082 // loop all points of first candidate
1083 for (ip=0; ip<numPoints1; ip++) {
1084 memcpy(local, &points1[3*ip], 3*sizeof(Double_t));
1085 if (local[0]<1e-10 && local[1]<1e-10) continue;
1086 mat1->LocalToMaster(local, point);
1087 mat2->MasterToLocal(point, local); // now point in local reference of second
1088 Bool_t overlap = shape2->Contains(local);
1089 if (overlap) {
1090 safety = shape2->Safety(local, kTRUE);
1091 if (safety<ovlp) overlap=kFALSE;
1092 }
1093 if (overlap) {
1094 if (!isoverlapping) {
1095 isoverlapping = kTRUE;
1096 nodeovlp = new TGeoOverlap(name,vol1,vol2,mat1,mat2,kTRUE,safety);
1097 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1098 fGeoManager->AddOverlap(nodeovlp);
1099 } else {
1100 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1101 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1102 }
1103 }
1104 }
1105 // loop all points of second candidate
1106 for (ip=0; ip<numPoints2; ip++) {
1107 memcpy(local, &points2[3*ip], 3*sizeof(Double_t));
1108 if (local[0]<1e-10 && local[1]<1e-10) continue;
1109 mat2->LocalToMaster(local, point);
1110 mat1->MasterToLocal(point, local); // now point in local reference of first
1111 Bool_t overlap = shape1->Contains(local);
1112 if (overlap) {
1113 safety = shape1->Safety(local, kTRUE);
1114 if (safety<ovlp) overlap=kFALSE;
1115 }
1116 if (overlap) {
1117 if (!isoverlapping) {
1118 isoverlapping = kTRUE;
1119 nodeovlp = new TGeoOverlap(name,vol1,vol2,mat1,mat2,kTRUE,safety);
1120 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1121 fGeoManager->AddOverlap(nodeovlp);
1122 } else {
1123 if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1124 nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1125 }
1126 }
1127 }
1128 return nodeovlp;
1129}
1130
1131////////////////////////////////////////////////////////////////////////////////
1132/// Check illegal overlaps for volume VOL within a limit OVLP by sampling npoints
1133/// inside the volume shape.
1134
1135void TGeoChecker::CheckOverlapsBySampling(TGeoVolume *vol, Double_t /* ovlp */, Int_t npoints) const
1136{
1137 Int_t nd = vol->GetNdaughters();
1138 if (nd<2) return;
1139 TGeoVoxelFinder *voxels = vol->GetVoxels();
1140 if (!voxels) return;
1141 if (voxels->NeedRebuild()) {
1142 voxels->Voxelize();
1143 vol->FindOverlaps();
1144 }
1145 TGeoBBox *box = (TGeoBBox*)vol->GetShape();
1146 TGeoShape *shape;
1147 TGeoNode *node;
1148 Double_t dx = box->GetDX();
1149 Double_t dy = box->GetDY();
1150 Double_t dz = box->GetDZ();
1151 Double_t pt[3];
1152 Double_t local[3];
1153 Int_t *check_list = nullptr;
1154 Int_t ncheck = 0;
1155 const Double_t *orig = box->GetOrigin();
1156 Int_t ipoint = 0;
1157 Int_t itry = 0;
1158 Int_t iovlp = 0;
1159 Int_t id=0, id0=0, id1=0;
1160 Bool_t in, incrt;
1161 Double_t safe;
1162 TString name1 = "";
1163 TString name2 = "";
1164 TGeoOverlap **flags = nullptr;
1165 TGeoNode *node1, *node2;
1166 Int_t novlps = 0;
1167 TGeoHMatrix mat1, mat2;
1168// Int_t tid = TGeoManager::ThreadId();
1169 TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
1170 TGeoStateInfo &td = *nav->GetCache()->GetInfo();
1171 while (ipoint < npoints) {
1172 // Shoot randomly in the bounding box.
1173 pt[0] = orig[0] - dx + 2.*dx*gRandom->Rndm();
1174 pt[1] = orig[1] - dy + 2.*dy*gRandom->Rndm();
1175 pt[2] = orig[2] - dz + 2.*dz*gRandom->Rndm();
1176 if (!vol->Contains(pt)) {
1177 itry++;
1178 if (itry>10000 && !ipoint) {
1179 Error("CheckOverlapsBySampling", "No point inside volume!!! - aborting");
1180 break;
1181 }
1182 continue;
1183 }
1184 // Check if the point is inside one or more daughters
1185 in = kFALSE;
1186 ipoint++;
1187 check_list = voxels->GetCheckList(pt, ncheck, td);
1188 if (!check_list || ncheck<2) continue;
1189 for (id=0; id<ncheck; id++) {
1190 id0 = check_list[id];
1191 node = vol->GetNode(id0);
1192 // Ignore MANY's
1193 if (node->IsOverlapping()) continue;
1194 node->GetMatrix()->MasterToLocal(pt,local);
1195 shape = node->GetVolume()->GetShape();
1196 incrt = shape->Contains(local);
1197 if (!incrt) continue;
1198 if (!in) {
1199 in = kTRUE;
1200 id1 = id0;
1201 continue;
1202 }
1203 // The point is inside 2 or more daughters, check safety
1204 safe = shape->Safety(local, kTRUE);
1205// if (safe < ovlp) continue;
1206 // We really have found an overlap -> store the point in a container
1207 iovlp++;
1208 if (!novlps) {
1209 if (flags) delete [] flags; // should never happen
1210 flags = new TGeoOverlap*[nd*nd];
1211 memset(flags, 0, nd*nd*sizeof(TGeoOverlap*));
1212 }
1213 TGeoOverlap *nodeovlp = flags[nd*id1+id0];
1214 if (!nodeovlp) {
1215 novlps++;
1216 node1 = vol->GetNode(id1);
1217 name1 = node1->GetName();
1218 mat1 = node1->GetMatrix();
1219 Int_t cindex = node1->GetVolume()->GetCurrentNodeIndex();
1220 while (cindex >= 0) {
1221 node1 = node1->GetVolume()->GetNode(cindex);
1222 name1 += TString::Format("/%s", node1->GetName());
1223 mat1.Multiply(node1->GetMatrix());
1224 cindex = node1->GetVolume()->GetCurrentNodeIndex();
1225 }
1226 node2 = vol->GetNode(id0);
1227 name2 = node2->GetName();
1228 mat2 = node2->GetMatrix();
1229 cindex = node2->GetVolume()->GetCurrentNodeIndex();
1230 while (cindex >= 0) {
1231 node2 = node2->GetVolume()->GetNode(cindex);
1232 name2 += TString::Format("/%s", node2->GetName());
1233 mat2.Multiply(node2->GetMatrix());
1234 cindex = node2->GetVolume()->GetCurrentNodeIndex();
1235 }
1236 nodeovlp = new TGeoOverlap(TString::Format("Volume %s: node %s overlapping %s",
1237 vol->GetName(), name1.Data(), name2.Data()), node1->GetVolume(),node2->GetVolume(),
1238 &mat1,&mat2, kTRUE, safe);
1239 flags[nd*id1+id0] = nodeovlp;
1240 fGeoManager->AddOverlap(nodeovlp);
1241 }
1242 // Max 100 points per marker
1243 if (nodeovlp->GetPolyMarker()->GetN()<100) nodeovlp->SetNextPoint(pt[0],pt[1],pt[2]);
1244 if (nodeovlp->GetOverlap()<safe) nodeovlp->SetOverlap(safe);
1245 }
1246 }
1247 nav->GetCache()->ReleaseInfo();
1248 if (flags) delete [] flags;
1249 if (!novlps) return;
1250 Double_t capacity = vol->GetShape()->Capacity();
1251 capacity *= Double_t(iovlp)/Double_t(npoints);
1252 Double_t err = 1./TMath::Sqrt(Double_t(iovlp));
1253 Info("CheckOverlapsBySampling", "#Found %d overlaps adding-up to %g +/- %g [cm3] for daughters of %s",
1254 novlps, capacity, err*capacity, vol->GetName());
1255}
1256
1257////////////////////////////////////////////////////////////////////////////////
1258/// Compute number of overlaps combinations to check per volume
1259
1260Int_t TGeoChecker::NChecksPerVolume(TGeoVolume *vol)
1261{
1262 if (vol->GetFinder()) return 0;
1263 UInt_t nd = vol->GetNdaughters();
1264 if (!nd) return 0;
1265 Bool_t is_assembly = vol->IsAssembly();
1266 TGeoIterator next1(vol);
1267 TGeoIterator next2(vol);
1268 Int_t nchecks = 0;
1269 TGeoNode *node;
1270 UInt_t id;
1271 if (!is_assembly) {
1272 while ((node=next1())) {
1273 if (node->IsOverlapping()) {
1274 next1.Skip();
1275 continue;
1276 }
1277 if (!node->GetVolume()->IsAssembly()) {
1278 nchecks++;
1279 next1.Skip();
1280 }
1281 }
1282 }
1283 // now check if the daughters overlap with each other
1284 if (nd<2) return nchecks;
1285 TGeoVoxelFinder *vox = vol->GetVoxels();
1286 if (!vox) return nchecks;
1287
1288
1289 TGeoNode *node1, *node01, *node02;
1290 Int_t novlp;
1291 Int_t *ovlps;
1292 Int_t ko;
1293 UInt_t io;
1294 for (id=0; id<nd; id++) {
1295 node01 = vol->GetNode(id);
1296 if (node01->IsOverlapping()) continue;
1297 vox->FindOverlaps(id);
1298 ovlps = node01->GetOverlaps(novlp);
1299 if (!ovlps) continue;
1300 for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
1301 io = ovlps[ko]; // (node1, shaped, matrix1, points, fBuff1)
1302 if (io<=id) continue;
1303 node02 = vol->GetNode(io);
1304 if (node02->IsOverlapping()) continue;
1305
1306 // We have to check node against node1, but they may be assemblies
1307 if (node01->GetVolume()->IsAssembly()) {
1308 next1.Reset(node01->GetVolume());
1309 while ((node=next1())) {
1310 if (!node->GetVolume()->IsAssembly()) {
1311 if (node02->GetVolume()->IsAssembly()) {
1312 next2.Reset(node02->GetVolume());
1313 while ((node1=next2())) {
1314 if (!node1->GetVolume()->IsAssembly()) {
1315 nchecks++;
1316 next2.Skip();
1317 }
1318 }
1319 } else {
1320 nchecks++;
1321 }
1322 next1.Skip();
1323 }
1324 }
1325 } else {
1326 // node not assembly
1327 if (node02->GetVolume()->IsAssembly()) {
1328 next2.Reset(node02->GetVolume());
1329 while ((node1=next2())) {
1330 if (!node1->GetVolume()->IsAssembly()) {
1331 nchecks++;
1332 next2.Skip();
1333 }
1334 }
1335 } else {
1336 // node1 also not assembly
1337 nchecks++;
1338 }
1339 }
1340 }
1341 node01->SetOverlaps(0,0);
1342 }
1343 return nchecks;
1344}
1345
1346////////////////////////////////////////////////////////////////////////////////
1347/// Check illegal overlaps for volume VOL within a limit OVLP.
1348
1349void TGeoChecker::CheckOverlaps(const TGeoVolume *vol, Double_t ovlp, Option_t *option)
1350{
1351 if (vol->GetFinder()) return;
1352 UInt_t nd = vol->GetNdaughters();
1353 if (!nd) return;
1354 TGeoShape::SetTransform(gGeoIdentity);
1355 fNchecks = NChecksPerVolume((TGeoVolume*)vol);
1356 Bool_t sampling = kFALSE;
1357 TString opt(option);
1358 opt.ToLower();
1359 if (opt.Contains("s")) sampling = kTRUE;
1360 if (opt.Contains("f")) fFullCheck = kTRUE;
1361 else fFullCheck = kFALSE;
1362 if (sampling) {
1363 opt = opt.Strip(TString::kLeading, 's');
1364 Int_t npoints = atoi(opt.Data());
1365 if (!npoints) npoints = 1000000;
1366 CheckOverlapsBySampling((TGeoVolume*)vol, ovlp, npoints);
1367 return;
1368 }
1369 Bool_t is_assembly = vol->IsAssembly();
1370 TGeoIterator next1((TGeoVolume*)vol);
1371 TGeoIterator next2((TGeoVolume*)vol);
1372 TString path;
1373 // first, test if daughters extrude their container
1374 TGeoNode * node, *nodecheck;
1375 TGeoChecker *checker = (TGeoChecker*)this;
1376
1377// TGeoOverlap *nodeovlp = 0;
1378 UInt_t id;
1379 Int_t level;
1380// Check extrusion only for daughters of a non-assembly volume
1381 if (!is_assembly) {
1382 while ((node=next1())) {
1383 if (node->IsOverlapping()) {
1384 next1.Skip();
1385 continue;
1386 }
1387 if (!node->GetVolume()->IsAssembly()) {
1388 if (fSelectedNode) {
1389 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1390 if ((fSelectedNode != node) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1391 next1.Skip();
1392 continue;
1393 }
1394 if (node != fSelectedNode) {
1395 level = next1.GetLevel();
1396 while ((nodecheck=next1.GetNode(level--))) {
1397 if (nodecheck == fSelectedNode) break;
1398 }
1399 if (!nodecheck) {
1400 next1.Skip();
1401 continue;
1402 }
1403 }
1404 }
1405 next1.GetPath(path);
1406 checker->MakeCheckOverlap(TString::Format("%s extruded by: %s", vol->GetName(),path.Data()),
1407 (TGeoVolume*)vol,node->GetVolume(),gGeoIdentity,(TGeoMatrix*)next1.GetCurrentMatrix(),kFALSE,ovlp);
1408 next1.Skip();
1409 }
1410 }
1411 }
1412
1413 // now check if the daughters overlap with each other
1414 if (nd<2) return;
1415 TGeoVoxelFinder *vox = vol->GetVoxels();
1416 if (!vox) {
1417 Warning("CheckOverlaps", "Volume %s with %i daughters but not voxelized", vol->GetName(),nd);
1418 return;
1419 }
1420 if (vox->NeedRebuild()) {
1421 vox->Voxelize();
1422 vol->FindOverlaps();
1423 }
1424 TGeoNode *node1, *node01, *node02;
1425 TGeoHMatrix hmat1, hmat2;
1426 TString path1;
1427 Int_t novlp;
1428 Int_t *ovlps;
1429 Int_t ko;
1430 UInt_t io;
1431 for (id=0; id<nd; id++) {
1432 node01 = vol->GetNode(id);
1433 if (node01->IsOverlapping()) continue;
1434 vox->FindOverlaps(id);
1435 ovlps = node01->GetOverlaps(novlp);
1436 if (!ovlps) continue;
1437 next1.SetTopName(node01->GetName());
1438 path = node01->GetName();
1439 for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
1440 io = ovlps[ko]; // (node1, shaped, matrix1, points, fBuff1)
1441 if (io<=id) continue;
1442 node02 = vol->GetNode(io);
1443 if (node02->IsOverlapping()) continue;
1444 // Try to fasten-up things...
1445// if (!TGeoBBox::AreOverlapping((TGeoBBox*)node01->GetVolume()->GetShape(), node01->GetMatrix(),
1446// (TGeoBBox*)node02->GetVolume()->GetShape(), node02->GetMatrix())) continue;
1447 next2.SetTopName(node02->GetName());
1448 path1 = node02->GetName();
1449
1450 // We have to check node against node1, but they may be assemblies
1451 if (node01->GetVolume()->IsAssembly()) {
1452 next1.Reset(node01->GetVolume());
1453 while ((node=next1())) {
1454 if (!node->GetVolume()->IsAssembly()) {
1455 next1.GetPath(path);
1456 hmat1 = node01->GetMatrix();
1457 hmat1 *= *next1.GetCurrentMatrix();
1458 if (node02->GetVolume()->IsAssembly()) {
1459 next2.Reset(node02->GetVolume());
1460 while ((node1=next2())) {
1461 if (!node1->GetVolume()->IsAssembly()) {
1462 if (fSelectedNode) {
1463 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1464 if ((fSelectedNode != node) && (fSelectedNode != node1) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1465 next2.Skip();
1466 continue;
1467 }
1468 if ((fSelectedNode != node1) && (fSelectedNode != node)) {
1469 level = next2.GetLevel();
1470 while ((nodecheck=next2.GetNode(level--))) {
1471 if (nodecheck == fSelectedNode) break;
1472 }
1473 if (node02 == fSelectedNode) nodecheck = node02;
1474 if (!nodecheck) {
1475 level = next1.GetLevel();
1476 while ((nodecheck=next1.GetNode(level--))) {
1477 if (nodecheck == fSelectedNode) break;
1478 }
1479 }
1480 if (node01 == fSelectedNode) nodecheck = node01;
1481 if (!nodecheck) {
1482 next2.Skip();
1483 continue;
1484 }
1485 }
1486 }
1487 next2.GetPath(path1);
1488 hmat2 = node02->GetMatrix();
1489 hmat2 *= *next2.GetCurrentMatrix();
1490 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1491 node->GetVolume(),node1->GetVolume(),&hmat1,&hmat2,kTRUE,ovlp);
1492 next2.Skip();
1493 }
1494 }
1495 } else {
1496 if (fSelectedNode) {
1497 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1498 if ((fSelectedNode != node) && (fSelectedNode != node02) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1499 next1.Skip();
1500 continue;
1501 }
1502 if ((fSelectedNode != node) && (fSelectedNode != node02)) {
1503 level = next1.GetLevel();
1504 while ((nodecheck=next1.GetNode(level--))) {
1505 if (nodecheck == fSelectedNode) break;
1506 }
1507 if (node01 == fSelectedNode) nodecheck = node01;
1508 if (!nodecheck) {
1509 next1.Skip();
1510 continue;
1511 }
1512 }
1513 }
1514 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1515 node->GetVolume(),node02->GetVolume(),&hmat1,node02->GetMatrix(),kTRUE,ovlp);
1516 }
1517 next1.Skip();
1518 }
1519 }
1520 } else {
1521 // node not assembly
1522 if (node02->GetVolume()->IsAssembly()) {
1523 next2.Reset(node02->GetVolume());
1524 while ((node1=next2())) {
1525 if (!node1->GetVolume()->IsAssembly()) {
1526 if (fSelectedNode) {
1527 // We have to check only overlaps of the selected node (or real daughters if an assembly)
1528 if ((fSelectedNode != node1) && (fSelectedNode != node01) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1529 next2.Skip();
1530 continue;
1531 }
1532 if ((fSelectedNode != node1) && (fSelectedNode != node01)) {
1533 level = next2.GetLevel();
1534 while ((nodecheck=next2.GetNode(level--))) {
1535 if (nodecheck == fSelectedNode) break;
1536 }
1537 if (node02 == fSelectedNode) nodecheck = node02;
1538 if (!nodecheck) {
1539 next2.Skip();
1540 continue;
1541 }
1542 }
1543 }
1544 next2.GetPath(path1);
1545 hmat2 = node02->GetMatrix();
1546 hmat2 *= *next2.GetCurrentMatrix();
1547 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1548 node01->GetVolume(),node1->GetVolume(),node01->GetMatrix(),&hmat2,kTRUE,ovlp);
1549 next2.Skip();
1550 }
1551 }
1552 } else {
1553 // node1 also not assembly
1554 if (fSelectedNode && (fSelectedNode != node01) && (fSelectedNode != node02)) continue;
1555 checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1556 node01->GetVolume(),node02->GetVolume(),node01->GetMatrix(),node02->GetMatrix(),kTRUE,ovlp);
1557 }
1558 }
1559 }
1560 node01->SetOverlaps(0,0);
1561 }
1562}
1563
1564////////////////////////////////////////////////////////////////////////////////
1565/// Print the current list of overlaps held by the manager class.
1566
1567void TGeoChecker::PrintOverlaps() const
1568{
1569 TIter next(fGeoManager->GetListOfOverlaps());
1570 TGeoOverlap *ov;
1571 printf("=== Overlaps for %s ===\n", fGeoManager->GetName());
1572 while ((ov=(TGeoOverlap*)next())) ov->PrintInfo();
1573}
1574
1575////////////////////////////////////////////////////////////////////////////////
1576/// Draw point (x,y,z) over the picture of the daughters of the volume containing this point.
1577/// Generates a report regarding the path to the node containing this point and the distance to
1578/// the closest boundary.
1579
1580void TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z, Option_t *)
1581{
1582 Double_t point[3];
1583 Double_t local[3];
1584 point[0] = x;
1585 point[1] = y;
1586 point[2] = z;
1587 TGeoVolume *vol = fGeoManager->GetTopVolume();
1588 if (fVsafe) {
1589 TGeoNode *old = fVsafe->GetNode("SAFETY_1");
1590 if (old) fVsafe->GetNodes()->RemoveAt(vol->GetNdaughters()-1);
1591 }
1592// if (vol != fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
1593 TGeoNode *node = fGeoManager->FindNode(point[0], point[1], point[2]);
1594 fGeoManager->MasterToLocal(point, local);
1595 // get current node
1596 printf("=== Check current point : (%g, %g, %g) ===\n", point[0], point[1], point[2]);
1597 printf(" - path : %s\n", fGeoManager->GetPath());
1598 // get corresponding volume
1599 if (node) vol = node->GetVolume();
1600 // compute safety distance (distance to boundary ignored)
1601 Double_t close = fGeoManager->Safety();
1602 printf("Safety radius : %f\n", close);
1603 if (close>1E-4) {
1604 TGeoVolume *sph = fGeoManager->MakeSphere("SAFETY", vol->GetMedium(), 0, close, 0,180,0,360);
1605 sph->SetLineColor(2);
1606 sph->SetLineStyle(3);
1607 vol->AddNode(sph,1,new TGeoTranslation(local[0], local[1], local[2]));
1608 fVsafe = vol;
1609 }
1610 TPolyMarker3D *pm = new TPolyMarker3D();
1611 pm->SetMarkerColor(2);
1612 pm->SetMarkerStyle(8);
1613 pm->SetMarkerSize(0.5);
1614 pm->SetNextPoint(local[0], local[1], local[2]);
1615 if (vol->GetNdaughters()<2) fGeoManager->SetTopVisible();
1616 else fGeoManager->SetTopVisible(kFALSE);
1617 fGeoManager->SetVisLevel(1);
1618 if (!vol->IsVisible()) vol->SetVisibility(kTRUE);
1619 vol->Draw();
1620 pm->Draw("SAME");
1621 gPad->Modified();
1622 gPad->Update();
1623}
1624
1625////////////////////////////////////////////////////////////////////////////////
1626/// Test for shape navigation methods. Summary for test numbers:
1627/// - 1: DistFromInside/Outside. Sample points inside the shape. Generate
1628/// directions randomly in cos(theta). Compute DistFromInside and move the
1629/// point with bigger distance. Compute DistFromOutside back from new point.
1630/// Plot d-(d1+d2)
1631/// - 2: Safety test. Sample points inside the bounding and compute safety. Generate
1632/// directions randomly in cos(theta) and compute distance to boundary. Check if
1633/// distance to boundary is bigger than safety
1634
1635void TGeoChecker::CheckShape(TGeoShape *shape, Int_t testNo, Int_t nsamples, Option_t *option)
1636{
1637 switch (testNo) {
1638 case 1:
1639 ShapeDistances(shape, nsamples, option);
1640 break;
1641 case 2:
1642 ShapeSafety(shape, nsamples, option);
1643 break;
1644 case 3:
1645 ShapeNormal(shape, nsamples, option);
1646 break;
1647 default:
1648 Error("CheckShape", "Test number %d not existent", testNo);
1649 }
1650}
1651
1652////////////////////////////////////////////////////////////////////////////////
1653/// Test TGeoShape::DistFromInside/Outside. Sample points inside the shape. Generate
1654/// directions randomly in cos(theta). Compute d1 = DistFromInside and move the
1655/// point on the boundary. Compute DistFromOutside and propagate with d2 making sure that
1656/// the shape is not re-entered. Swap direction and call DistFromOutside that
1657/// should fall back on the same point on the boundary (at d2). Propagate back on boundary
1658/// then compute DistFromInside that should be bigger than d1.
1659/// Plot d-(d1+d2)
1660
1661void TGeoChecker::ShapeDistances(TGeoShape *shape, Int_t nsamples, Option_t *)
1662{
1663 Double_t dx = ((TGeoBBox*)shape)->GetDX();
1664 Double_t dy = ((TGeoBBox*)shape)->GetDY();
1665 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1666 Double_t dmax = 2.*TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1667 Double_t d1, d2, dmove, dnext;
1668 Int_t itot = 0;
1669 // Number of tracks shot for every point inside the shape
1670 const Int_t kNtracks = 1000;
1671 Int_t n10 = nsamples/10;
1672 Int_t i,j;
1673 Double_t point[3], pnew[3];
1674 Double_t dir[3], dnew[3];
1675 Double_t theta, phi, delta;
1676 TPolyMarker3D *pmfrominside = 0;
1677 TPolyMarker3D *pmfromoutside = 0;
1678 TH1D *hist = new TH1D("hTest1", "Residual distance from inside/outside",200,-20, 0);
1679 hist->GetXaxis()->SetTitle("delta[cm] - first bin=overflow");
1680 hist->GetYaxis()->SetTitle("count");
1681 hist->SetMarkerStyle(kFullCircle);
1682
1683 if (!fTimer) fTimer = new TStopwatch();
1684 fTimer->Reset();
1685 fTimer->Start();
1686 while (itot<nsamples) {
1687 Bool_t inside = kFALSE;
1688 while (!inside) {
1689 point[0] = gRandom->Uniform(-dx,dx);
1690 point[1] = gRandom->Uniform(-dy,dy);
1691 point[2] = gRandom->Uniform(-dz,dz);
1692 inside = shape->Contains(point);
1693 }
1694 itot++;
1695 if (n10) {
1696 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1697 }
1698 for (i=0; i<kNtracks; i++) {
1699 phi = 2*TMath::Pi()*gRandom->Rndm();
1700 theta= TMath::ACos(1.-2.*gRandom->Rndm());
1701 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
1702 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
1703 dir[2]=TMath::Cos(theta);
1704 dmove = dmax;
1705 // We have track direction, compute distance from inside
1706 d1 = shape->DistFromInside(point,dir,3);
1707 if (d1>dmove || d1<TGeoShape::Tolerance()) {
1708 // Bad distance or bbox size, to debug
1709 new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1710 shape->Draw();
1711 printf("DistFromInside: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) %f/%f(max)\n",
1712 point[0],point[1],point[2],dir[0],dir[1],dir[2], d1,dmove);
1713 pmfrominside = new TPolyMarker3D(2);
1714 pmfrominside->SetMarkerColor(kRed);
1715 pmfrominside->SetMarkerStyle(24);
1716 pmfrominside->SetMarkerSize(0.4);
1717 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1718 for (j=0; j<3; j++) pnew[j] = point[j] + d1*dir[j];
1719 pmfrominside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1720 pmfrominside->Draw();
1721 return;
1722 }
1723 // Propagate BEFORE the boundary and make sure that DistFromOutside
1724 // does not return 0 (!!!)
1725 // Check if there is a second crossing
1726 for (j=0; j<3; j++) pnew[j] = point[j] + (d1-TGeoShape::Tolerance())*dir[j];
1727 dnext = shape->DistFromOutside(pnew,dir,3);
1728 if (d1+dnext<dmax) dmove = d1+0.5*dnext;
1729 // Move point and swap direction
1730 for (j=0; j<3; j++) {
1731 pnew[j] = point[j] + dmove*dir[j];
1732 dnew[j] = -dir[j];
1733 }
1734 // Compute now distance from outside
1735 d2 = shape->DistFromOutside(pnew,dnew,3);
1736 delta = dmove-d1-d2;
1737 if (TMath::Abs(delta)>1E-6 || dnext<2.*TGeoShape::Tolerance()) {
1738 // Error->debug this
1739 new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1740 shape->Draw();
1741 printf("Error: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d2=%f dmove=%f\n",
1742 point[0],point[1],point[2],dir[0],dir[1],dir[2], d1,d2,dmove);
1743 if (dnext<2.*TGeoShape::Tolerance()) {
1744 printf(" (*)DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) dnext = %f\n",
1745 point[0]+(d1-TGeoShape::Tolerance())*dir[0],
1746 point[1]+(d1-TGeoShape::Tolerance())*dir[1],
1747 point[2]+(d1-TGeoShape::Tolerance())*dir[2], dir[0],dir[1],dir[2],dnext);
1748 } else {
1749 printf(" DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) dnext = %f\n",
1750 point[0]+d1*dir[0],point[1]+d1*dir[1], point[2]+d1*dir[2], dir[0],dir[1],dir[2],dnext);
1751 }
1752 printf(" DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) = %f\n",
1753 pnew[0],pnew[1],pnew[2],dnew[0],dnew[1],dnew[2], d2);
1754 pmfrominside = new TPolyMarker3D(2);
1755 pmfrominside->SetMarkerStyle(24);
1756 pmfrominside->SetMarkerSize(0.4);
1757 pmfrominside->SetMarkerColor(kRed);
1758 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1759 for (j=0; j<3; j++) point[j] += d1*dir[j];
1760 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1761 pmfrominside->Draw();
1762 pmfromoutside = new TPolyMarker3D(2);
1763 pmfromoutside->SetMarkerStyle(20);
1764 pmfromoutside->SetMarkerStyle(7);
1765 pmfromoutside->SetMarkerSize(0.3);
1766 pmfromoutside->SetMarkerColor(kBlue);
1767 pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1768 for (j=0; j<3; j++) pnew[j] += d2*dnew[j];
1769 if (d2<1E10) pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1770 pmfromoutside->Draw();
1771 return;
1772 }
1773 // Compute distance from inside which should be bigger than d1
1774 for (j=0; j<3; j++) pnew[j] += (d2-TGeoShape::Tolerance())*dnew[j];
1775 dnext = shape->DistFromInside(pnew,dnew,3);
1776 if (dnext<d1-TGeoShape::Tolerance() || dnext>dmax) {
1777 new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1778 shape->Draw();
1779 printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d1p=%f\n",
1780 pnew[0],pnew[1],pnew[2],dnew[0],dnew[1],dnew[2],d1,dnext);
1781 pmfrominside = new TPolyMarker3D(2);
1782 pmfrominside->SetMarkerStyle(24);
1783 pmfrominside->SetMarkerSize(0.4);
1784 pmfrominside->SetMarkerColor(kRed);
1785 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1786 for (j=0; j<3; j++) point[j] += d1*dir[j];
1787 pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1788 pmfrominside->Draw();
1789 pmfromoutside = new TPolyMarker3D(2);
1790 pmfromoutside->SetMarkerStyle(20);
1791 pmfromoutside->SetMarkerStyle(7);
1792 pmfromoutside->SetMarkerSize(0.3);
1793 pmfromoutside->SetMarkerColor(kBlue);
1794 pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1795 for (j=0; j<3; j++) pnew[j] += dnext*dnew[j];
1796 if (d2<1E10) pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1797 pmfromoutside->Draw();
1798 return;
1799 }
1800 if (TMath::Abs(delta) < 1E-20) delta = 1E-30;
1801 hist->Fill(TMath::Max(TMath::Log(TMath::Abs(delta)),-20.));
1802 }
1803 }
1804 fTimer->Stop();
1805 fTimer->Print();
1806 new TCanvas("Test01", "Residuals DistFromInside/Outside", 800, 600);
1807 hist->Draw();
1808}
1809
1810////////////////////////////////////////////////////////////////////////////////
1811/// Check of validity of safe distance for a given shape.
1812/// Sample points inside the 2x bounding box and compute safety. Generate
1813/// directions randomly in cos(theta) and compute distance to boundary. Check if
1814/// distance to boundary is bigger than safety.
1815
1816void TGeoChecker::ShapeSafety(TGeoShape *shape, Int_t nsamples, Option_t *)
1817{
1818 Double_t dx = ((TGeoBBox*)shape)->GetDX();
1819 Double_t dy = ((TGeoBBox*)shape)->GetDY();
1820 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1821 // Number of tracks shot for every point inside the shape
1822 const Int_t kNtracks = 1000;
1823 Int_t n10 = nsamples/10;
1824 Int_t i;
1825 Double_t dist;
1826 Double_t point[3];
1827 Double_t dir[3];
1828 Double_t theta, phi;
1829 TPolyMarker3D *pm1 = 0;
1830 TPolyMarker3D *pm2 = 0;
1831 if (!fTimer) fTimer = new TStopwatch();
1832 fTimer->Reset();
1833 fTimer->Start();
1834 Int_t itot = 0;
1835 while (itot<nsamples) {
1836 Bool_t inside = kFALSE;
1837 point[0] = gRandom->Uniform(-2*dx,2*dx);
1838 point[1] = gRandom->Uniform(-2*dy,2*dy);
1839 point[2] = gRandom->Uniform(-2*dz,2*dz);
1840 inside = shape->Contains(point);
1841 Double_t safe = shape->Safety(point, inside);
1842 itot++;
1843 if (n10) {
1844 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1845 }
1846 for (i=0; i<kNtracks; i++) {
1847 phi = 2*TMath::Pi()*gRandom->Rndm();
1848 theta= TMath::ACos(1.-2.*gRandom->Rndm());
1849 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
1850 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
1851 dir[2]=TMath::Cos(theta);
1852 if (inside) dist = shape->DistFromInside(point,dir,3);
1853 else dist = shape->DistFromOutside(point,dir,3);
1854 if (dist<safe) {
1855 printf("Error safety (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) safe=%f dist=%f\n",
1856 point[0],point[1],point[2], dir[0], dir[1], dir[2], safe, dist);
1857 new TCanvas("shape02", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1858 shape->Draw();
1859 pm1 = new TPolyMarker3D(2);
1860 pm1->SetMarkerStyle(24);
1861 pm1->SetMarkerSize(0.4);
1862 pm1->SetMarkerColor(kRed);
1863 pm1->SetNextPoint(point[0],point[1],point[2]);
1864 pm1->SetNextPoint(point[0]+safe*dir[0],point[1]+safe*dir[1],point[2]+safe*dir[2]);
1865 pm1->Draw();
1866 pm2 = new TPolyMarker3D(1);
1867 pm2->SetMarkerStyle(7);
1868 pm2->SetMarkerSize(0.3);
1869 pm2->SetMarkerColor(kBlue);
1870 pm2->SetNextPoint(point[0]+dist*dir[0],point[1]+dist*dir[1],point[2]+dist*dir[2]);
1871 pm2->Draw();
1872 return;
1873 }
1874 }
1875 }
1876 fTimer->Stop();
1877 fTimer->Print();
1878}
1879
1880////////////////////////////////////////////////////////////////////////////////
1881/// Check of validity of the normal for a given shape.
1882/// Sample points inside the shape. Generate directions randomly in cos(theta)
1883/// and propagate to boundary. Compute normal and safety at crossing point, plot
1884/// the point and generate a random direction so that (dir) dot (norm) <0.
1885
1886void TGeoChecker::ShapeNormal(TGeoShape *shape, Int_t nsamples, Option_t *)
1887{
1888 Double_t dx = ((TGeoBBox*)shape)->GetDX();
1889 Double_t dy = ((TGeoBBox*)shape)->GetDY();
1890 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1891 Double_t dmax = 2.*TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1892 // Number of tracks shot for every point inside the shape
1893 const Int_t kNtracks = 1000;
1894 Int_t n10 = nsamples/10;
1895 Int_t itot = 0, errcnt = 0, errsame=0;
1896 Int_t i;
1897 Double_t dist, olddist, safe, dot;
1898 Double_t theta, phi, ndotd;
1899 Double_t *spoint = new Double_t[3*nsamples];
1900 Double_t *sdir = new Double_t[3*nsamples];
1901 while (itot<nsamples) {
1902 Bool_t inside = kFALSE;
1903 while (!inside) {
1904 spoint[3*itot] = gRandom->Uniform(-dx,dx);
1905 spoint[3*itot+1] = gRandom->Uniform(-dy,dy);
1906 spoint[3*itot+2] = gRandom->Uniform(-dz,dz);
1907 inside = shape->Contains(&spoint[3*itot]);
1908 }
1909 phi = 2*TMath::Pi()*gRandom->Rndm();
1910 theta= TMath::ACos(1.-2.*gRandom->Rndm());
1911 sdir[3*itot] = TMath::Sin(theta)*TMath::Cos(phi);
1912 sdir[3*itot+1] = TMath::Sin(theta)*TMath::Sin(phi);
1913 sdir[3*itot+2] = TMath::Cos(theta);
1914 itot++;
1915 }
1916 Double_t point[3],newpoint[3], oldpoint[3];
1917 Double_t dir[3], olddir[3];
1918 Double_t norm[3], newnorm[3], oldnorm[3];
1919 TCanvas *errcanvas = 0;
1920 TPolyMarker3D *pm1 = 0;
1921 TPolyMarker3D *pm2 = 0;
1922 pm2 = new TPolyMarker3D();
1923// pm2->SetMarkerStyle(24);
1924 pm2->SetMarkerSize(0.2);
1925 pm2->SetMarkerColor(kBlue);
1926 if (!fTimer) fTimer = new TStopwatch();
1927 fTimer->Reset();
1928 fTimer->Start();
1929 for (itot = 0; itot<nsamples; itot++) {
1930 if (n10) {
1931 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1932 }
1933 oldnorm[0] = oldnorm[1] = oldnorm[2] = 0.;
1934 olddist = 0.;
1935 for (Int_t j=0; j<3; j++) {
1936 oldpoint[j] = point[j] = spoint[3*itot+j];
1937 olddir[j] = dir[j] = sdir[3*itot+j];
1938 }
1939 for (i=0; i<kNtracks; i++) {
1940 if (errcnt>0) break;
1941 dist = shape->DistFromInside(point,dir,3);
1942 for (Int_t j=0; j<3; j++) {
1943 newpoint[j] = point[j] + dist*dir[j];
1944 }
1945 shape->ComputeNormal(newpoint,dir,newnorm);
1946
1947 dot = olddir[0]*oldnorm[0]+olddir[1]*oldnorm[1]+ olddir[2]*oldnorm[2];
1948 if (!shape->Contains(point) && shape->Safety(point,kFALSE)>1.E-3) {
1949 errcnt++;
1950 printf("Error point outside (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
1951 point[0],point[1],point[2], dir[0], dir[1], dir[2], dist, olddist);
1952 printf(" old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n",
1953 oldpoint[0],oldpoint[1],oldpoint[2], olddir[0], olddir[1], olddir[2]);
1954 if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1955 if (!pm1) {
1956 pm1 = new TPolyMarker3D();
1957 pm1->SetMarkerStyle(24);
1958 pm1->SetMarkerSize(0.4);
1959 pm1->SetMarkerColor(kRed);
1960 }
1961 pm1->SetNextPoint(point[0],point[1],point[2]);
1962 pm1->SetNextPoint(oldpoint[0],oldpoint[1],oldpoint[2]);
1963 break;
1964 }
1965 if ((dist<TGeoShape::Tolerance() && olddist*dot>1.E-3) || dist>dmax) {
1966 errsame++;
1967 if (errsame>1) {
1968 errcnt++;
1969 printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
1970 point[0],point[1],point[2], dir[0], dir[1], dir[2], dist, olddist);
1971 printf(" new norm: (%g, %g, %g)\n", newnorm[0], newnorm[1], newnorm[2]);
1972 printf(" old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n",
1973 oldpoint[0],oldpoint[1],oldpoint[2], olddir[0], olddir[1], olddir[2]);
1974 printf(" old norm: (%g, %g, %g)\n", oldnorm[0], oldnorm[1], oldnorm[2]);
1975 if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1976 if (!pm1) {
1977 pm1 = new TPolyMarker3D();
1978 pm1->SetMarkerStyle(24);
1979 pm1->SetMarkerSize(0.4);
1980 pm1->SetMarkerColor(kRed);
1981 }
1982 pm1->SetNextPoint(point[0],point[1],point[2]);
1983 pm1->SetNextPoint(oldpoint[0],oldpoint[1],oldpoint[2]);
1984 break;
1985 }
1986 } else errsame = 0;
1987
1988 olddist = dist;
1989 for (Int_t j=0; j<3; j++) {
1990 oldpoint[j] = point[j];
1991 point[j] += dist*dir[j];
1992 }
1993 safe = shape->Safety(point, kTRUE);
1994 if (safe>1.E-3) {
1995 errcnt++;
1996 printf("Error safety (%19.15f, %19.15f, %19.15f) safe=%g\n",
1997 point[0],point[1],point[2], safe);
1998 if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1999 if (!pm1) {
2000 pm1 = new TPolyMarker3D();
2001 pm1->SetMarkerStyle(24);
2002 pm1->SetMarkerSize(0.4);
2003 pm1->SetMarkerColor(kRed);
2004 }
2005 pm1->SetNextPoint(point[0],point[1],point[2]);
2006 break;
2007 }
2008 // Compute normal
2009 shape->ComputeNormal(point,dir,norm);
2010 memcpy(oldnorm, norm, 3*sizeof(Double_t));
2011 memcpy(olddir, dir, 3*sizeof(Double_t));
2012 while (1) {
2013 phi = 2*TMath::Pi()*gRandom->Rndm();
2014 theta= TMath::ACos(1.-2.*gRandom->Rndm());
2015 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
2016 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
2017 dir[2]=TMath::Cos(theta);
2018 ndotd = dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2];
2019 if (ndotd<0) break; // backwards, still inside shape
2020 }
2021 if ((itot%10) == 0) pm2->SetNextPoint(point[0],point[1],point[2]);
2022 }
2023 }
2024 fTimer->Stop();
2025 fTimer->Print();
2026 if (errcanvas) {
2027 shape->Draw();
2028 pm1->Draw();
2029 }
2030
2031 new TCanvas("shape03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
2032 pm2->Draw();
2033 delete [] spoint;
2034 delete [] sdir;
2035}
2036
2037////////////////////////////////////////////////////////////////////////////////
2038/// Generate a lego plot fot the top volume, according to option.
2039
2040TH2F *TGeoChecker::LegoPlot(Int_t ntheta, Double_t themin, Double_t themax,
2041 Int_t nphi, Double_t phimin, Double_t phimax,
2042 Double_t /*rmin*/, Double_t /*rmax*/, Option_t *option)
2043{
2044 TH2F *hist = new TH2F("lego", option, nphi, phimin, phimax, ntheta, themin, themax);
2045
2046 Double_t degrad = TMath::Pi()/180.;
2047 Double_t theta, phi, step, matprop, x;
2048 Double_t start[3];
2049 Double_t dir[3];
2050 TGeoNode *startnode, *endnode;
2051 Int_t i; // loop index for phi
2052 Int_t j; // loop index for theta
2053 Int_t ntot = ntheta * nphi;
2054 Int_t n10 = ntot/10;
2055 Int_t igen = 0, iloop=0;
2056 printf("=== Lego plot sph. => nrays=%i\n", ntot);
2057 for (i=1; i<=nphi; i++) {
2058 for (j=1; j<=ntheta; j++) {
2059 igen++;
2060 if (n10) {
2061 if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/ntot));
2062 }
2063 x = 0;
2064 theta = hist->GetYaxis()->GetBinCenter(j);
2065 phi = hist->GetXaxis()->GetBinCenter(i)+1E-3;
2066 start[0] = start[1] = start[2] = 1E-3;
2067 dir[0]=TMath::Sin(theta*degrad)*TMath::Cos(phi*degrad);
2068 dir[1]=TMath::Sin(theta*degrad)*TMath::Sin(phi*degrad);
2069 dir[2]=TMath::Cos(theta*degrad);
2070 fGeoManager->InitTrack(&start[0], &dir[0]);
2071 startnode = fGeoManager->GetCurrentNode();
2072 if (fGeoManager->IsOutside()) startnode=0;
2073 if (startnode) {
2074 matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
2075 } else {
2076 matprop = 0.;
2077 }
2078 fGeoManager->FindNextBoundary();
2079// fGeoManager->IsStepEntering();
2080 // find where we end-up
2081 endnode = fGeoManager->Step();
2082 step = fGeoManager->GetStep();
2083 while (step<1E10) {
2084 // now see if we can make an other step
2085 iloop=0;
2086 while (!fGeoManager->IsEntering()) {
2087 iloop++;
2088 fGeoManager->SetStep(1E-3);
2089 step += 1E-3;
2090 endnode = fGeoManager->Step();
2091 }
2092 if (iloop>1000) printf("%i steps\n", iloop);
2093 if (matprop>0) {
2094 x += step/matprop;
2095 }
2096 if (endnode==0 && step>1E10) break;
2097 // generate an extra step to cross boundary
2098 startnode = endnode;
2099 if (startnode) {
2100 matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
2101 } else {
2102 matprop = 0.;
2103 }
2104
2105 fGeoManager->FindNextBoundary();
2106 endnode = fGeoManager->Step();
2107 step = fGeoManager->GetStep();
2108 }
2109 hist->Fill(phi, theta, x);
2110 }
2111 }
2112 return hist;
2113}
2114
2115////////////////////////////////////////////////////////////////////////////////
2116/// Draw random points in the bounding box of a volume.
2117
2118void TGeoChecker::RandomPoints(TGeoVolume *vol, Int_t npoints, Option_t *option)
2119{
2120 if (!vol) return;
2121 vol->VisibleDaughters(kTRUE);
2122 vol->Draw();
2123 TString opt = option;
2124 opt.ToLower();
2125 TObjArray *pm = new TObjArray(128);
2126 TPolyMarker3D *marker = 0;
2127 const TGeoShape *shape = vol->GetShape();
2128 TGeoBBox *box = (TGeoBBox *)shape;
2129 Double_t dx = box->GetDX();
2130 Double_t dy = box->GetDY();
2131 Double_t dz = box->GetDZ();
2132 Double_t ox = (box->GetOrigin())[0];
2133 Double_t oy = (box->GetOrigin())[1];
2134 Double_t oz = (box->GetOrigin())[2];
2135 Double_t *xyz = new Double_t[3];
2136 printf("Random box : %f, %f, %f\n", dx, dy, dz);
2137 TGeoNode *node = 0;
2138 printf("Start... %i points\n", npoints);
2139 Int_t i=0;
2140 Int_t igen=0;
2141 Int_t ic = 0;
2142 Int_t n10 = npoints/10;
2143 Double_t ratio=0;
2144 while (igen<npoints) {
2145 xyz[0] = ox-dx+2*dx*gRandom->Rndm();
2146 xyz[1] = oy-dy+2*dy*gRandom->Rndm();
2147 xyz[2] = oz-dz+2*dz*gRandom->Rndm();
2148 fGeoManager->SetCurrentPoint(xyz);
2149 igen++;
2150 if (n10) {
2151 if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/npoints));
2152 }
2153 node = fGeoManager->FindNode();
2154 if (!node) continue;
2155 if (!node->IsOnScreen()) continue;
2156 // draw only points in overlapping/non-overlapping volumes
2157 if (opt.Contains("many") && !node->IsOverlapping()) continue;
2158 if (opt.Contains("only") && node->IsOverlapping()) continue;
2159 ic = node->GetColour();
2160 if ((ic<0) || (ic>=128)) ic = 1;
2161 marker = (TPolyMarker3D*)pm->At(ic);
2162 if (!marker) {
2163 marker = new TPolyMarker3D();
2164 marker->SetMarkerColor(ic);
2165// marker->SetMarkerStyle(8);
2166// marker->SetMarkerSize(0.4);
2167 pm->AddAt(marker, ic);
2168 }
2169 marker->SetNextPoint(xyz[0], xyz[1], xyz[2]);
2170 i++;
2171 }
2172 printf("Number of visible points : %i\n", i);
2173 if (igen) ratio = (Double_t)i/(Double_t)igen;
2174 printf("efficiency : %g\n", ratio);
2175 for (Int_t m=0; m<128; m++) {
2176 marker = (TPolyMarker3D*)pm->At(m);
2177 if (marker) marker->Draw("SAME");
2178 }
2179 fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
2180 printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
2181 printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
2182 delete pm;
2183 delete [] xyz;
2184}
2185
2186////////////////////////////////////////////////////////////////////////////////
2187/// Randomly shoot nrays from point (startx,starty,startz) and plot intersections
2188/// with surfaces for current top node.
2189
2190void TGeoChecker::RandomRays(Int_t nrays, Double_t startx, Double_t starty, Double_t startz, const char *target_vol, Bool_t check_norm)
2191{
2192 TObjArray *pm = new TObjArray(128);
2193 TString starget = target_vol;
2194 TPolyLine3D *line = 0;
2195 TPolyLine3D *normline = 0;
2196 TGeoVolume *vol=fGeoManager->GetTopVolume();
2197// vol->VisibleDaughters(kTRUE);
2198
2199 Bool_t random = kFALSE;
2200 if (nrays<=0) {
2201 nrays = 100000;
2202 random = kTRUE;
2203 }
2204 Double_t start[3];
2205 Double_t dir[3];
2206 const Double_t *point = fGeoManager->GetCurrentPoint();
2207 vol->Draw();
2208 printf("Start... %i rays\n", nrays);
2209 TGeoNode *startnode, *endnode;
2210 Bool_t vis1,vis2;
2211 Int_t i=0;
2212 Int_t ipoint, inull;
2213 Int_t itot=0;
2214 Int_t n10=nrays/10;
2215 Double_t theta,phi, step, normlen;
2216 Double_t ox = ((TGeoBBox*)vol->GetShape())->GetOrigin()[0];
2217 Double_t oy = ((TGeoBBox*)vol->GetShape())->GetOrigin()[1];
2218 Double_t oz = ((TGeoBBox*)vol->GetShape())->GetOrigin()[2];
2219 Double_t dx = ((TGeoBBox*)vol->GetShape())->GetDX();
2220 Double_t dy = ((TGeoBBox*)vol->GetShape())->GetDY();
2221 Double_t dz = ((TGeoBBox*)vol->GetShape())->GetDZ();
2222 normlen = TMath::Max(dx,dy);
2223 normlen = TMath::Max(normlen,dz);
2224 normlen *= 0.05;
2225 while (itot<nrays) {
2226 itot++;
2227 inull = 0;
2228 ipoint = 0;
2229 if (n10) {
2230 if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nrays));
2231 }
2232 if (random) {
2233 start[0] = ox-dx+2*dx*gRandom->Rndm();
2234 start[1] = oy-dy+2*dy*gRandom->Rndm();
2235 start[2] = oz-dz+2*dz*gRandom->Rndm();
2236 } else {
2237 start[0] = startx;
2238 start[1] = starty;
2239 start[2] = startz;
2240 }
2241 phi = 2*TMath::Pi()*gRandom->Rndm();
2242 theta= TMath::ACos(1.-2.*gRandom->Rndm());
2243 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
2244 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
2245 dir[2]=TMath::Cos(theta);
2246 startnode = fGeoManager->InitTrack(start[0],start[1],start[2], dir[0],dir[1],dir[2]);
2247 line = 0;
2248 if (fGeoManager->IsOutside()) startnode=0;
2249 vis1 = kFALSE;
2250 if (target_vol) {
2251 if (startnode && starget==startnode->GetVolume()->GetName()) vis1 = kTRUE;
2252 } else {
2253 if (startnode && startnode->IsOnScreen()) vis1 = kTRUE;
2254 }
2255 if (vis1) {
2256 line = new TPolyLine3D(2);
2257 line->SetLineColor(startnode->GetVolume()->GetLineColor());
2258 line->SetPoint(ipoint++, start[0], start[1], start[2]);
2259 i++;
2260 pm->Add(line);
2261 }
2262 while ((endnode=fGeoManager->FindNextBoundaryAndStep())) {
2263 step = fGeoManager->GetStep();
2264 if (step<TGeoShape::Tolerance()) inull++;
2265 else inull = 0;
2266 if (inull>5) break;
2267 const Double_t *normal = 0;
2268 if (check_norm) {
2269 normal = fGeoManager->FindNormalFast();
2270 if (!normal) break;
2271 }
2272 vis2 = kFALSE;
2273 if (target_vol) {
2274 if (starget==endnode->GetVolume()->GetName()) vis2 = kTRUE;
2275 } else if (endnode->IsOnScreen()) vis2 = kTRUE;
2276 if (ipoint>0) {
2277 // old visible node had an entry point -> finish segment
2278 line->SetPoint(ipoint, point[0], point[1], point[2]);
2279 if (!vis2 && check_norm) {
2280 normline = new TPolyLine3D(2);
2281 normline->SetLineColor(kBlue);
2282 normline->SetLineWidth(1);
2283 normline->SetPoint(0, point[0], point[1], point[2]);
2284 normline->SetPoint(1, point[0]+normal[0]*normlen,
2285 point[1]+normal[1]*normlen,
2286 point[2]+normal[2]*normlen);
2287 pm->Add(normline);
2288 }
2289 ipoint = 0;
2290 line = 0;
2291 }
2292 if (vis2) {
2293 // create new segment
2294 line = new TPolyLine3D(2);
2295 line->SetLineColor(endnode->GetVolume()->GetLineColor());
2296 line->SetPoint(ipoint++, point[0], point[1], point[2]);
2297 i++;
2298 if (check_norm) {
2299 normline = new TPolyLine3D(2);
2300 normline->SetLineColor(kBlue);
2301 normline->SetLineWidth(2);
2302 normline->SetPoint(0, point[0], point[1], point[2]);
2303 normline->SetPoint(1, point[0]+normal[0]*normlen,
2304 point[1]+normal[1]*normlen,
2305 point[2]+normal[2]*normlen);
2306 }
2307 pm->Add(line);
2308 if (!random) pm->Add(normline);
2309 }
2310 }
2311 }
2312 // draw all segments
2313 for (Int_t m=0; m<pm->GetEntriesFast(); m++) {
2314 line = (TPolyLine3D*)pm->At(m);
2315 if (line) line->Draw("SAME");
2316 }
2317 printf("number of segments : %i\n", i);
2318// fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
2319// printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
2320// printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
2321 delete pm;
2322}
2323
2324////////////////////////////////////////////////////////////////////////////////
2325/// shoot npoints randomly in a box of 1E-5 around current point.
2326/// return minimum distance to points outside
2327/// make sure that path to current node is updated
2328/// get the response of tgeo
2329
2330TGeoNode *TGeoChecker::SamplePoints(Int_t npoints, Double_t &dist, Double_t epsil,
2331 const char* g3path)
2332{
2333 TGeoNode *node = fGeoManager->FindNode();
2334 TGeoNode *nodegeo = 0;
2335 TGeoNode *nodeg3 = 0;
2336 TGeoNode *solg3 = 0;
2337 if (!node) {dist=-1; return 0;}
2338 Bool_t hasg3 = kFALSE;
2339 if (strlen(g3path)) hasg3 = kTRUE;
2340 TString geopath = fGeoManager->GetPath();
2341 dist = 1E10;
2342 TString common = "";
2343 // cd to common path
2344 Double_t point[3];
2345 Double_t closest[3];
2346 TGeoNode *node1 = 0;
2347 TGeoNode *node_close = 0;
2348 dist = 1E10;
2349 Double_t dist1 = 0;
2350 // initialize size of random box to epsil
2351 Double_t eps[3];
2352 eps[0] = epsil; eps[1]=epsil; eps[2]=epsil;
2353 const Double_t *pointg = fGeoManager->GetCurrentPoint();
2354 if (hasg3) {
2355 TString spath = geopath;
2356 TString name = "";
2357 Int_t index=0;
2358 while (index>=0) {
2359 index = spath.Index("/", index+1);
2360 if (index>0) {
2361 name = spath(0, index);
2362 if (strstr(g3path, name.Data())) {
2363 common = name;
2364 continue;
2365 } else break;
2366 }
2367 }
2368 // if g3 response was given, cd to common path
2369 if (strlen(common.Data())) {
2370 while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
2371 nodegeo = fGeoManager->GetCurrentNode();
2372 fGeoManager->CdUp();
2373 }
2374 fGeoManager->cd(g3path);
2375 solg3 = fGeoManager->GetCurrentNode();
2376 while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
2377 nodeg3 = fGeoManager->GetCurrentNode();
2378 fGeoManager->CdUp();
2379 }
2380 if (!nodegeo) return 0;
2381 if (!nodeg3) return 0;
2382 fGeoManager->cd(common.Data());
2383 fGeoManager->MasterToLocal(fGeoManager->GetCurrentPoint(), &point[0]);
2384 Double_t xyz[3], local[3];
2385 for (Int_t i=0; i<npoints; i++) {
2386 xyz[0] = point[0] - eps[0] + 2*eps[0]*gRandom->Rndm();
2387 xyz[1] = point[1] - eps[1] + 2*eps[1]*gRandom->Rndm();
2388 xyz[2] = point[2] - eps[2] + 2*eps[2]*gRandom->Rndm();
2389 nodeg3->MasterToLocal(&xyz[0], &local[0]);
2390 if (!nodeg3->GetVolume()->Contains(&local[0])) continue;
2391 dist1 = TMath::Sqrt((xyz[0]-point[0])*(xyz[0]-point[0])+
2392 (xyz[1]-point[1])*(xyz[1]-point[1])+(xyz[2]-point[2])*(xyz[2]-point[2]));
2393 if (dist1<dist) {
2394 // save node and closest point
2395 dist = dist1;
2396 node_close = solg3;
2397 // make the random box smaller
2398 eps[0] = TMath::Abs(point[0]-pointg[0]);
2399 eps[1] = TMath::Abs(point[1]-pointg[1]);
2400 eps[2] = TMath::Abs(point[2]-pointg[2]);
2401 }
2402 }
2403 }
2404 if (!node_close) dist = -1;
2405 return node_close;
2406 }
2407
2408 // save current point
2409 memcpy(&point[0], pointg, 3*sizeof(Double_t));
2410 for (Int_t i=0; i<npoints; i++) {
2411 // generate a random point in MARS
2412 fGeoManager->SetCurrentPoint(point[0] - eps[0] + 2*eps[0]*gRandom->Rndm(),
2413 point[1] - eps[1] + 2*eps[1]*gRandom->Rndm(),
2414 point[2] - eps[2] + 2*eps[2]*gRandom->Rndm());
2415 // check if new node is different from the old one
2416 if (node1!=node) {
2417 dist1 = TMath::Sqrt((point[0]-pointg[0])*(point[0]-pointg[0])+
2418 (point[1]-pointg[1])*(point[1]-pointg[1])+(point[2]-pointg[2])*(point[2]-pointg[2]));
2419 if (dist1<dist) {
2420 dist = dist1;
2421 node_close = node1;
2422 memcpy(&closest[0], pointg, 3*sizeof(Double_t));
2423 // make the random box smaller
2424 eps[0] = TMath::Abs(point[0]-pointg[0]);
2425 eps[1] = TMath::Abs(point[1]-pointg[1]);
2426 eps[2] = TMath::Abs(point[2]-pointg[2]);
2427 }
2428 }
2429 }
2430 // restore the original point and path
2431 fGeoManager->FindNode(point[0], point[1], point[2]); // really needed ?
2432 if (!node_close) dist=-1;
2433 return node_close;
2434}
2435
2436////////////////////////////////////////////////////////////////////////////////
2437/// Shoot one ray from start point with direction (dirx,diry,dirz). Fills input array
2438/// with points just after boundary crossings.
2439
2440Double_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
2441{
2442// Int_t array_dimension = 3*dim;
2443 nelem = 0;
2444 Int_t istep = 0;
2445 if (!dim) {
2446 printf("empty input array\n");
2447 return array;
2448 }
2449// fGeoManager->CdTop();
2450 const Double_t *point = fGeoManager->GetCurrentPoint();
2451 TGeoNode *endnode;
2452 Bool_t is_entering;
2453 Double_t step, forward;
2454 Double_t dir[3];
2455 dir[0] = dirx;
2456 dir[1] = diry;
2457 dir[2] = dirz;
2458 fGeoManager->InitTrack(start, &dir[0]);
2459 fGeoManager->GetCurrentNode();
2460// printf("Start : (%f,%f,%f)\n", point[0], point[1], point[2]);
2461 fGeoManager->FindNextBoundary();
2462 step = fGeoManager->GetStep();
2463// printf("---next : at step=%f\n", step);
2464 if (step>1E10) return array;
2465 endnode = fGeoManager->Step();
2466 is_entering = fGeoManager->IsEntering();
2467 while (step<1E10) {
2468 if (endpoint) {
2469 forward = dirx*(endpoint[0]-point[0])+diry*(endpoint[1]-point[1])+dirz*(endpoint[2]-point[2]);
2470 if (forward<1E-3) {
2471// printf("exit : Passed start point. nelem=%i\n", nelem);
2472 return array;
2473 }
2474 }
2475 if (is_entering) {
2476 if (nelem>=dim) {
2477 Double_t *temparray = new Double_t[3*(dim+20)];
2478 memcpy(temparray, array, 3*dim*sizeof(Double_t));
2479 delete [] array;
2480 array = temparray;
2481 dim += 20;
2482 }
2483 memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
2484// printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
2485 nelem++;
2486 } else {
2487 if (endnode==0 && step>1E10) {
2488// printf("exit : NULL endnode. nelem=%i\n", nelem);
2489 return array;
2490 }
2491 if (!fGeoManager->IsEntering()) {
2492// if (startnode) printf("stepping %f from (%f, %f, %f) inside %s...\n", step,point[0], point[1], point[2], startnode->GetName());
2493// else printf("stepping %f from (%f, %f, %f) OUTSIDE...\n", step,point[0], point[1], point[2]);
2494 istep = 0;
2495 }
2496 while (!fGeoManager->IsEntering()) {
2497 istep++;
2498 if (istep>1E3) {
2499// Error("ShootRay", "more than 1000 steps. Step was %f", step);
2500 nelem = 0;
2501 return array;
2502 }
2503 fGeoManager->SetStep(1E-5);
2504 endnode = fGeoManager->Step();
2505 }
2506 if (istep>0) printf("%i steps\n", istep);
2507 if (nelem>=dim) {
2508 Double_t *temparray = new Double_t[3*(dim+20)];
2509 memcpy(temparray, array, 3*dim*sizeof(Double_t));
2510 delete [] array;
2511 array = temparray;
2512 dim += 20;
2513 }
2514 memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
2515// if (endnode) printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
2516 nelem++;
2517 }
2518 fGeoManager->FindNextBoundary();
2519 step = fGeoManager->GetStep();
2520// printf("---next at step=%f\n", step);
2521 endnode = fGeoManager->Step();
2522 is_entering = fGeoManager->IsEntering();
2523 }
2524 return array;
2525// printf("exit : INFINITE step. nelem=%i\n", nelem);
2526}
2527
2528////////////////////////////////////////////////////////////////////////////////
2529/// Check time of finding "Where am I" for n points.
2530
2531void TGeoChecker::Test(Int_t npoints, Option_t *option)
2532{
2533 Bool_t recheck = !strcmp(option, "RECHECK");
2534 if (recheck) printf("RECHECK\n");
2535 const TGeoShape *shape = fGeoManager->GetTopVolume()->GetShape();
2536 Double_t dx = ((TGeoBBox*)shape)->GetDX();
2537 Double_t dy = ((TGeoBBox*)shape)->GetDY();
2538 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
2539 Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
2540 Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
2541 Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
2542 Double_t *xyz = new Double_t[3*npoints];
2543 TStopwatch *timer = new TStopwatch();
2544 printf("Random box : %f, %f, %f\n", dx, dy, dz);
2545 timer->Start(kFALSE);
2546 Int_t i;
2547 for (i=0; i<npoints; i++) {
2548 xyz[3*i] = ox-dx+2*dx*gRandom->Rndm();
2549 xyz[3*i+1] = oy-dy+2*dy*gRandom->Rndm();
2550 xyz[3*i+2] = oz-dz+2*dz*gRandom->Rndm();
2551 }
2552 timer->Stop();
2553 printf("Generation time :\n");
2554 timer->Print();
2555 timer->Reset();
2556 TGeoNode *node, *node1;
2557 printf("Start... %i points\n", npoints);
2558 timer->Start(kFALSE);
2559 for (i=0; i<npoints; i++) {
2560 fGeoManager->SetCurrentPoint(xyz+3*i);
2561 if (recheck) fGeoManager->CdTop();
2562 node = fGeoManager->FindNode();
2563 if (recheck) {
2564 node1 = fGeoManager->FindNode();
2565 if (node1 != node) {
2566 printf("Difference for x=%g y=%g z=%g\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
2567 printf(" from top : %s\n", node->GetName());
2568 printf(" redo : %s\n", fGeoManager->GetPath());
2569 }
2570 }
2571 }
2572 timer->Stop();
2573 timer->Print();
2574 delete [] xyz;
2575 delete timer;
2576}
2577
2578////////////////////////////////////////////////////////////////////////////////
2579/// Geometry overlap checker based on sampling.
2580
2581void TGeoChecker::TestOverlaps(const char* path)
2582{
2583 if (fGeoManager->GetTopVolume()!=fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
2584 printf("Checking overlaps for path :\n");
2585 if (!fGeoManager->cd(path)) return;
2586 TGeoNode *checked = fGeoManager->GetCurrentNode();
2587 checked->InspectNode();
2588 // shoot 1E4 points in the shape of the current volume
2589 Int_t npoints = 1000000;
2590 Double_t big = 1E6;
2591 Double_t xmin = big;
2592 Double_t xmax = -big;
2593 Double_t ymin = big;
2594 Double_t ymax = -big;
2595 Double_t zmin = big;
2596 Double_t zmax = -big;
2597 TObjArray *pm = new TObjArray(128);
2598 TPolyMarker3D *marker = 0;
2599 TPolyMarker3D *markthis = new TPolyMarker3D();
2600 markthis->SetMarkerColor(5);
2601 TNtuple *ntpl = new TNtuple("ntpl","random points","x:y:z");
2602 TGeoShape *shape = fGeoManager->GetCurrentNode()->GetVolume()->GetShape();
2603 Double_t *point = new Double_t[3];
2604 Double_t dx = ((TGeoBBox*)shape)->GetDX();
2605 Double_t dy = ((TGeoBBox*)shape)->GetDY();
2606 Double_t dz = ((TGeoBBox*)shape)->GetDZ();
2607 Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
2608 Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
2609 Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
2610 Double_t *xyz = new Double_t[3*npoints];
2611 Int_t i=0;
2612 printf("Generating %i points inside %s\n", npoints, fGeoManager->GetPath());
2613 while (i<npoints) {
2614 point[0] = ox-dx+2*dx*gRandom->Rndm();
2615 point[1] = oy-dy+2*dy*gRandom->Rndm();
2616 point[2] = oz-dz+2*dz*gRandom->Rndm();
2617 if (!shape->Contains(point)) continue;
2618 // convert each point to MARS
2619// printf("local %9.3f %9.3f %9.3f\n", point[0], point[1], point[2]);
2620 fGeoManager->LocalToMaster(point, &xyz[3*i]);
2621// printf("master %9.3f %9.3f %9.3f\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
2622 xmin = TMath::Min(xmin, xyz[3*i]);
2623 xmax = TMath::Max(xmax, xyz[3*i]);
2624 ymin = TMath::Min(ymin, xyz[3*i+1]);
2625 ymax = TMath::Max(ymax, xyz[3*i+1]);
2626 zmin = TMath::Min(zmin, xyz[3*i+2]);
2627 zmax = TMath::Max(zmax, xyz[3*i+2]);
2628 i++;
2629 }
2630 delete [] point;
2631 ntpl->Fill(xmin,ymin,zmin);
2632 ntpl->Fill(xmax,ymin,zmin);
2633 ntpl->Fill(xmin,ymax,zmin);
2634 ntpl->Fill(xmax,ymax,zmin);
2635 ntpl->Fill(xmin,ymin,zmax);
2636 ntpl->Fill(xmax,ymin,zmax);
2637 ntpl->Fill(xmin,ymax,zmax);
2638 ntpl->Fill(xmax,ymax,zmax);
2639 ntpl->Draw("z:y:x");
2640
2641 // shoot the poins in the geometry
2642 TGeoNode *node;
2643 TString cpath;
2644 Int_t ic=0;
2645 TObjArray *overlaps = new TObjArray();
2646 printf("using FindNode...\n");
2647 for (Int_t j=0; j<npoints; j++) {
2648 // always start from top level (testing only)
2649 fGeoManager->CdTop();
2650 fGeoManager->SetCurrentPoint(&xyz[3*j]);
2651 node = fGeoManager->FindNode();
2652 cpath = fGeoManager->GetPath();
2653 if (cpath.Contains(path)) {
2654 markthis->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
2655 continue;
2656 }
2657 // current point is found in an overlapping node
2658 if (!node) ic=128;
2659 else ic = node->GetVolume()->GetLineColor();
2660 if (ic >= 128) ic = 0;
2661 marker = (TPolyMarker3D*)pm->At(ic);
2662 if (!marker) {
2663 marker = new TPolyMarker3D();
2664 marker->SetMarkerColor(ic);
2665 pm->AddAt(marker, ic);
2666 }
2667 // draw the overlapping point
2668 marker->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
2669 if (node) {
2670 if (overlaps->IndexOf(node) < 0) overlaps->Add(node);
2671 }
2672 }
2673 // draw all overlapping points
2674// for (Int_t m=0; m<128; m++) {
2675// marker = (TPolyMarker3D*)pm->At(m);
2676// if (marker) marker->Draw("SAME");
2677// }
2678 markthis->Draw("SAME");
2679 if (gPad) gPad->Update();
2680 // display overlaps
2681 if (overlaps->GetEntriesFast()) {
2682 printf("list of overlapping nodes :\n");
2683 for (i=0; i<overlaps->GetEntriesFast(); i++) {
2684 node = (TGeoNode*)overlaps->At(i);
2685 if (node->IsOverlapping()) printf("%s MANY\n", node->GetName());
2686 else printf("%s ONLY\n", node->GetName());
2687 }
2688 } else printf("No overlaps\n");
2689 delete ntpl;
2690 delete pm;
2691 delete [] xyz;
2692 delete overlaps;
2693}
2694
2695////////////////////////////////////////////////////////////////////////////////
2696/// Estimate weight of top level volume with a precision SIGMA(W)/W
2697/// better than PRECISION. Option can be "v" - verbose (default).
2698
2699Double_t TGeoChecker::Weight(Double_t precision, Option_t *option)
2700{
2701 TList *matlist = fGeoManager->GetListOfMaterials();
2702 Int_t nmat = matlist->GetSize();
2703 if (!nmat) return 0;
2704 Int_t *nin = new Int_t[nmat];
2705 memset(nin, 0, nmat*sizeof(Int_t));
2706 TString opt = option;
2707 opt.ToLower();
2708 Bool_t isverbose = opt.Contains("v");
2709 TGeoBBox *box = (TGeoBBox *)fGeoManager->GetTopVolume()->GetShape();
2710 Double_t dx = box->GetDX();
2711 Double_t dy = box->GetDY();
2712 Double_t dz = box->GetDZ();
2713 Double_t ox = (box->GetOrigin())[0];
2714 Double_t oy = (box->GetOrigin())[1];
2715 Double_t oz = (box->GetOrigin())[2];
2716 Double_t x,y,z;
2717 TGeoNode *node;
2718 TGeoMaterial *mat;
2719 Double_t vbox = 0.000008*dx*dy*dz; // m3
2720 Bool_t end = kFALSE;
2721 Double_t weight=0, sigma, eps, dens;
2722 Double_t eps0=1.;
2723 Int_t indmat;
2724 Int_t igen=0;
2725 Int_t iin = 0;
2726 while (!end) {
2727 x = ox-dx+2*dx*gRandom->Rndm();
2728 y = oy-dy+2*dy*gRandom->Rndm();
2729 z = oz-dz+2*dz*gRandom->Rndm();
2730 node = fGeoManager->FindNode(x,y,z);
2731 igen++;
2732 if (!node) continue;
2733 mat = node->GetVolume()->GetMedium()->GetMaterial();
2734 indmat = mat->GetIndex();
2735 if (indmat<0) continue;
2736 nin[indmat]++;
2737 iin++;
2738 if ((iin%100000)==0 || igen>1E8) {
2739 weight = 0;
2740 sigma = 0;
2741 for (indmat=0; indmat<nmat; indmat++) {
2742 mat = (TGeoMaterial*)matlist->At(indmat);
2743 dens = mat->GetDensity(); // [g/cm3]
2744 if (dens<1E-2) dens=0;
2745 dens *= 1000.; // [kg/m3]
2746 weight += dens*Double_t(nin[indmat]);
2747 sigma += dens*dens*nin[indmat];
2748 }
2749 sigma = TMath::Sqrt(sigma);
2750 eps = sigma/weight;
2751 weight *= vbox/Double_t(igen);
2752 sigma *= vbox/Double_t(igen);
2753 if (eps<precision || igen>1E8) {
2754 if (isverbose) {
2755 printf("=== Weight of %s : %g +/- %g [kg]\n",
2756 fGeoManager->GetTopVolume()->GetName(), weight, sigma);
2757 }
2758 end = kTRUE;
2759 } else {
2760 if (isverbose && eps<0.5*eps0) {
2761 printf("%8dK: %14.7g kg %g %%\n",
2762 igen/1000, weight, eps*100);
2763 eps0 = eps;
2764 }
2765 }
2766 }
2767 }
2768 delete [] nin;
2769 return weight;
2770}
2771////////////////////////////////////////////////////////////////////////////////
2772/// count voxel timing
2773
2774Double_t TGeoChecker::CheckVoxels(TGeoVolume *vol, TGeoVoxelFinder *voxels, Double_t *xyz, Int_t npoints)
2775{
2776 TStopwatch timer;
2777 Double_t time;
2778 TGeoShape *shape = vol->GetShape();
2779 TGeoNode *node;
2780 TGeoMatrix *matrix;
2781 Double_t *point;
2782 Double_t local[3];
2783 Int_t *checklist;
2784 Int_t ncheck;
2785 TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
2786 TGeoStateInfo &td = *nav->GetCache()->GetInfo();
2787 timer.Start();
2788 for (Int_t i=0; i<npoints; i++) {
2789 point = xyz + 3*i;
2790 if (!shape->Contains(point)) continue;
2791 checklist = voxels->GetCheckList(point, ncheck, td);
2792 if (!checklist) continue;
2793 if (!ncheck) continue;
2794 for (Int_t id=0; id<ncheck; id++) {
2795 node = vol->GetNode(checklist[id]);
2796 matrix = node->GetMatrix();
2797 matrix->MasterToLocal(point, &local[0]);
2798 if (node->GetVolume()->GetShape()->Contains(&local[0])) break;
2799 }
2800 }
2801 nav->GetCache()->ReleaseInfo();
2802 time = timer.CpuTime();
2803 return time;
2804}
2805
2806////////////////////////////////////////////////////////////////////////////////
2807/// Returns optimal voxelization type for volume vol.
2808/// - kFALSE - cartesian
2809/// - kTRUE - cylindrical
2810
2811Bool_t TGeoChecker::TestVoxels(TGeoVolume * /*vol*/, Int_t /*npoints*/)
2812{
2813 return kFALSE;
2814}
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
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Int_t
int Int_t
Definition RtypesCore.h:45
Char_t
char Char_t
Definition RtypesCore.h:37
kFALSE
const Bool_t kFALSE
Definition RtypesCore.h:101
Long_t
long Long_t
Definition RtypesCore.h:54
Double_t
double Double_t
Definition RtypesCore.h:59
Long64_t
long long Long64_t
Definition RtypesCore.h:80
kTRUE
const Bool_t kTRUE
Definition RtypesCore.h:100
Option_t
const char Option_t
Definition RtypesCore.h:66
ClassImp
#define ClassImp(name)
Definition Rtypes.h:364
kRed
@ kRed
Definition Rtypes.h:66
kBlue
@ kBlue
Definition Rtypes.h:66
kYellow
@ kYellow
Definition Rtypes.h:66
kFullCircle
@ kFullCircle
Definition TAttMarker.h:51
id
XFontStruct * id
Definition TGX11.cxx:109
name
char name[80]
Definition TGX11.cxx:110
gGeoIdentity
R__EXTERN TGeoIdentity * gGeoIdentity
Definition TGeoMatrix.h:478
xmin
float xmin
Definition THbookFile.cxx:95
nentries
int nentries
Definition THbookFile.cxx:91
ymin
float ymin
Definition THbookFile.cxx:95
xmax
float xmax
Definition THbookFile.cxx:95
ymax
float ymax
Definition THbookFile.cxx:95
gRandom
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
Form
char * Form(const char *fmt,...)
gStyle
R__EXTERN TStyle * gStyle
Definition TStyle.h:413
gPad
#define gPad
Definition TVirtualPad.h:288
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:478
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:23
TCollection::GetSize
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition TCollection.h:184
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition TFile.h:54
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:1635
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:1135
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:1886
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:1349
TGeoChecker::TestOverlaps
void TestOverlaps(const char *path)
Geometry overlap checker based on sampling.
Definition TGeoChecker.cxx:2581
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:969
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:1661
TGeoChecker::NChecksPerVolume
Int_t NChecksPerVolume(TGeoVolume *vol)
Compute number of overlaps combinations to check per volume.
Definition TGeoChecker.cxx:1260
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:1816
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:2440
TGeoChecker::CleanPoints
void CleanPoints(Double_t *points, Int_t &numPoints) const
Number of points on mesh to be checked.
Definition TGeoChecker.cxx:948
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:2330
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:2699
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:772
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:1580
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:1567
TGeoChecker::Test
void Test(Int_t npoints, Option_t *option)
Check time of finding "Where am I" for n points.
Definition TGeoChecker.cxx:2531
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:2811
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:2190
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:2040
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:2118
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:2774
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:45
TGeoManager::GetMother
TGeoNode * GetMother(Int_t up=1) const
Definition TGeoManager.h:513
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:2589
TGeoManager::GetListOfUVolumes
TObjArray * GetListOfUVolumes() const
Definition TGeoManager.h:497
TGeoManager::GetListOfOverlaps
TObjArray * GetListOfOverlaps()
Definition TGeoManager.h:489
TGeoManager::CdUp
void CdUp()
Go one level up in geometry.
Definition TGeoManager.cxx:1539
TGeoManager::cd
virtual Bool_t cd(const char *path="")
Browse the tree of nodes starting from fTopNode according to pathname.
Definition TGeoManager.cxx:1565
TGeoManager::LocalToMaster
void LocalToMaster(const Double_t *local, Double_t *master) const
Definition TGeoManager.h:543
TGeoManager::RestoreMasterVolume
void RestoreMasterVolume()
Restore the master volume of the geometry.
Definition TGeoManager.cxx:2884
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:2556
TGeoManager::GetCurrentNavigator
TGeoNavigator * GetCurrentNavigator() const
Returns current navigator for the calling thread.
Definition TGeoManager.cxx:738
TGeoManager::GetMasterVolume
TGeoVolume * GetMasterVolume() const
Definition TGeoManager.h:531
TGeoManager::GetCurrentNode
TGeoNode * GetCurrentNode() const
Definition TGeoManager.h:519
TGeoManager::SetCurrentDirection
void SetCurrentDirection(Double_t *dir)
Definition TGeoManager.h:538
TGeoManager::SetVisLevel
void SetVisLevel(Int_t level=3)
set default level down to which visualization is performed
Definition TGeoManager.cxx:2260
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:2545
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:3391
TGeoManager::GetVolume
TGeoVolume * GetVolume(const char *name) const
Search for a named volume. All trailing blanks stripped.
Definition TGeoManager.cxx:2722
TGeoManager::FindNextBoundaryAndStep
TGeoNode * FindNextBoundaryAndStep(Double_t stepmax=TGeoShape::Big(), Bool_t compsafe=kFALSE)
Compute distance to next boundary within STEPMAX.
Definition TGeoManager.cxx:2526
TGeoManager::SetCurrentPoint
void SetCurrentPoint(Double_t *point)
Definition TGeoManager.h:535
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:2600
TGeoManager::FindNode
TGeoNode * FindNode(Bool_t safe_start=kTRUE)
Returns deepest node containing current point.
Definition TGeoManager.cxx:2572
TGeoManager::GetCurrentPoint
const Double_t * GetCurrentPoint() const
Definition TGeoManager.h:521
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:2949
TGeoManager::IsOutside
Bool_t IsOutside() const
Definition TGeoManager.h:428
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:2642
TGeoManager::GetLevel
Int_t GetLevel() const
Definition TGeoManager.h:527
TGeoManager::GetStep
Double_t GetStep() const
Definition TGeoManager.h:408
TGeoManager::GetCurrentMatrix
TGeoHMatrix * GetCurrentMatrix() const
Definition TGeoManager.h:516
TGeoManager::GetTopNode
TGeoNode * GetTopNode() const
Definition TGeoManager.h:533
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:2608
TGeoManager::GetPath
const char * GetPath() const
Get path to the current node in the form /node0/node1/...
Definition TGeoManager.cxx:2667
TGeoManager::AddOverlap
Int_t AddOverlap(const TNamed *ovlp)
Add an illegal overlap/extrusion to the list.
Definition TGeoManager.cxx:561
TGeoManager::SetVerboseLevel
static void SetVerboseLevel(Int_t vl)
Return current verbosity level (static function).
Definition TGeoManager.cxx:3807
TGeoManager::SetStep
void SetStep(Double_t step)
Definition TGeoManager.h:422
TGeoManager::GetCurrentVolume
TGeoVolume * GetCurrentVolume() const
Definition TGeoManager.h:523
TGeoManager::GetVerboseLevel
static Int_t GetVerboseLevel()
Set verbosity level (static function).
Definition TGeoManager.cxx:3799
TGeoManager::CdTop
void CdTop()
Make top level node the current node.
Definition TGeoManager.cxx:1530
TGeoManager::Safety
Double_t Safety(Bool_t inside=kFALSE)
Compute safe distance from the current point.
Definition TGeoManager.cxx:2127
TGeoManager::MasterToLocal
void MasterToLocal(const Double_t *master, Double_t *local) const
Definition TGeoManager.h:546
TGeoManager::ResetState
void ResetState()
Reset current state flags.
Definition TGeoManager.cxx:2564
TGeoManager::CdDown
void CdDown(Int_t index)
Make a daughter of current node current.
Definition TGeoManager.cxx:1548
TGeoManager::GetListOfMaterials
TList * GetListOfMaterials() const
Definition TGeoManager.h:491
TGeoManager::GetTopVolume
TGeoVolume * GetTopVolume() const
Definition TGeoManager.h:532
TGeoManager::SetTopVisible
void SetTopVisible(Bool_t vis=kTRUE)
make top volume visible on screen
Definition TGeoManager.cxx:2205
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:3612
TGeoManager::IsEntering
Bool_t IsEntering() const
Definition TGeoManager.h:424
TGeoMaterial
Base class describing materials.
Definition TGeoMaterial.h:36
TGeoMaterial::GetIndex
Int_t GetIndex()
Retrieve material index in the list of materials.
Definition TGeoMaterial.cxx:623
TGeoMaterial::GetRadLen
virtual Double_t GetRadLen() const
Definition TGeoMaterial.h:115
TGeoMaterial::GetDensity
virtual Double_t GetDensity() const
Definition TGeoMaterial.h:108
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:49
TGeoVolume::GetMedium
TGeoMedium * GetMedium() const
Definition TGeoVolume.h:174
TGeoVolume::Contains
Bool_t Contains(const Double_t *point) const
Definition TGeoVolume.h:109
TGeoVolume::GetMaterial
TGeoMaterial * GetMaterial() const
Definition TGeoVolume.h:173
TGeoVolume::GetNdaughters
Int_t GetNdaughters() const
Definition TGeoVolume.h:351
TGeoVolume::SetLineStyle
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TGeoVolume.cxx:2045
TGeoVolume::GetNodes
TObjArray * GetNodes()
Definition TGeoVolume.h:168
TGeoVolume::AddNode
virtual TGeoNode * AddNode(TGeoVolume *vol, Int_t copy_no, TGeoMatrix *mat=0, Option_t *option="")
Add a TGeoNode to the list of nodes.
Definition TGeoVolume.cxx:932
TGeoVolume::VisibleDaughters
void VisibleDaughters(Bool_t vis=kTRUE)
set visibility for daughters
Definition TGeoVolume.cxx:2293
TGeoVolume::FindOverlaps
void FindOverlaps() const
loop all nodes marked as overlaps and find overlapping brothers
Definition TGeoVolume.cxx:2091
TGeoVolume::GetNode
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
Definition TGeoVolume.cxx:2061
TGeoVolume::SetVisibility
virtual void SetVisibility(Bool_t vis=kTRUE)
set visibility of this volume
Definition TGeoVolume.cxx:2199
TGeoVolume::GetIndex
Int_t GetIndex(const TGeoNode *node) const
get index number for a given daughter
Definition TGeoVolume.cxx:1634
TGeoVolume::GetFinder
TGeoPatternFinder * GetFinder() const
Definition TGeoVolume.h:176
TGeoVolume::GetVoxels
TGeoVoxelFinder * GetVoxels() const
Getter for optimization structure.
Definition TGeoVolume.cxx:1688
TGeoVolume::GetNumber
Int_t GetNumber() const
Definition TGeoVolume.h:183
TGeoVolume::GetShape
TGeoShape * GetShape() const
Definition TGeoVolume.h:189
TGeoVolume::Draw
virtual void Draw(Option_t *option="")
draw top volume according to option
Definition TGeoVolume.cxx:1147
TGeoVolume::GetCurrentNodeIndex
virtual Int_t GetCurrentNodeIndex() const
Definition TGeoVolume.h:166
TGeoVolume::DrawOnly
virtual void DrawOnly(Option_t *option="")
draw only this volume
Definition TGeoVolume.cxx:1164
TGeoVolume::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TGeoVolume.cxx:2037
TGeoVolume::IsAssembly
virtual Bool_t IsAssembly() const
Returns true if the volume is an assembly or a scaled assembly.
Definition TGeoVolume.cxx:1706
TGeoVolume::IsVisible
virtual Bool_t IsVisible() const
Definition TGeoVolume.h:154
TGeoVolume::InspectShape
void InspectShape() const
Definition TGeoVolume.h:194
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:1205
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:2221
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:296
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:618
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:575
TH1::kAllAxes
@ kAllAxes
Definition TH1.h:75
TH1::LabelsOption
virtual void LabelsOption(Option_t *option="h", Option_t *axis="X")
Sort bins with labels or set option(s) to draw axis with labels.
Definition TH1.cxx:5315
TH1::AddDirectory
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition TH1.cxx:1283
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3351
TH1::GetYaxis
TAxis * GetYaxis()
Definition TH1.h:321
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:6598
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3074
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:5178
TH1::SetStats
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition TH1.cxx:8820
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:358
TList
A doubly linked list.
Definition TList.h:38
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:357
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:169
TObjArray
An array of TObjects.
Definition TObjArray.h:31
TObjArray::IndexOf
Int_t IndexOf(const TObject *obj) const
Definition TObjArray.cxx:605
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
Definition TObjArray.h:58
TObjArray::Add
void Add(TObject *obj)
Definition TObjArray.h:68
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:164
TObject
Mother of all ROOT objects.
Definition TObject.h:41
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:200
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:949
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:963
TObject::Draw
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:267
TObject::Info
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:937
TPolyLine3D
A 3-dimensional polyline.
Definition TPolyLine3D.h:33
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:621
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:542
TPolyMarker3D::Draw
virtual void Draw(Option_t *option="")
Draws 3-D polymarker with its current attributes.
Definition TPolyMarker3D.cxx:261
TRandom::Uniform
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:672
TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom.cxx:552
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:136
TString::ToLower
void ToLower()
Change string to lower-case.
Definition TString.cxx:1150
TString::Strip
TSubString Strip(EStripType s=kTrailing, char c=' ') const
Return a substring of self stripped at beginning and/or end.
Definition TString.cxx:1131
TString::Data
const char * Data() const
Definition TString.h:369
TString::kLeading
@ kLeading
Definition TString.h:267
TString::Remove
TString & Remove(Ssiz_t pos)
Definition TString.h:673
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:2336
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
TString::Index
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition TString.h:639
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:1589
TTree
A TTree represents a columnar dataset.
Definition TTree.h:79
TTree::Fill
virtual Int_t Fill()
Fill all branches.
Definition TTree.cxx:4594
TTree::GetEntry
virtual Int_t GetEntry(Long64_t entry, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition TTree.cxx:5622
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:8356
TTree::GetEntries
virtual Long64_t GetEntries() const
Definition TTree.h:460
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:350
TTree::Draw
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition TTree.h:428
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:9693
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
TMath::ACos
Double_t ACos(Double_t)
Definition TMath.h:619
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:208
TMath::Log
Double_t Log(Double_t x)
Definition TMath.h:710
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition TMath.h:641
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition TMathBase.h:176
TMath::Cos
Double_t Cos(Double_t)
Definition TMath.h:593
TMath::Pi
constexpr Double_t Pi()
Definition TMath.h:37
TMath::Sin
Double_t Sin(Double_t)
Definition TMath.h:589
TMath::Abs
Short_t Abs(Short_t d)
Definition TMathBase.h:120
TMath::TwoPi
constexpr Double_t TwoPi()
Definition TMath.h:44
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:987
lnext
#define lnext(otri1, otri2)
Definition triangle.c:943
vertex
REAL * vertex
Definition triangle.c:513