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