1 #include "TGeoManager.h"
2 #include "TGeoVolume.h"
10 #include "TStopwatch.h"
13 TStopwatch *timer = 0;
15 Int_t propagate_in_geom(Double_t *, Double_t *);
16 void score(TGeoVolume *, Int_t, Double_t);
17 Double_t timing_per_volume(TGeoVolume *);
24 TGeoManager::Import("geometry.root");
25 Int_t nvol = gGeoManager->GetListOfVolumes()->GetEntries();
26 Int_t nuid = gGeoManager->GetListOfUVolumes()->GetEntries();
27 timer = new TStopwatch();
30 flags = new Bool_t[nuid];
31 val1 = new Double_t[nuid];
32 val2 = new Double_t[nuid];
33 memset(flags, 0, nuid*sizeof(Bool_t));
34 memset(val1, 0, nuid*sizeof(Double_t));
35 memset(val2, 0, nuid*sizeof(Double_t));
38 // STAGE 1: Overlap checking by sampling per volume
40 printf("====================================================================\n");
41 printf("STAGE 1: Overlap checking by sampling per volume\n");
42 printf("====================================================================\n");
43 for (i=0; i<nvol; i++) {
44 vol = (TGeoVolume*)gGeoManager->GetListOfVolumes()->At(i);
45 Int_t uid = vol->GetNumber();
46 if (flags[uid]) continue;
48 if (!vol->GetNdaughters()) continue;
49 vol->CheckOverlaps(0.01, "s");
52 // STAGE 2: Global overlap checking
53 printf("====================================================================\n");
54 printf("STAGE 2: Global overlap checking\n");
55 printf("====================================================================\n");
56 gGeoManager->CheckOverlaps(0.01);
58 // STAGE 3: How many crossings per volume in a realistic case ?
59 // Ignore volumes with no daughters
61 // Generate rays from vertex in phi=[0,2*pi] theta=[0,pi]
62 Int_t ntracks = 1000000;
63 printf("====================================================================\n");
64 printf("STAGE 3: Propagating %i tracks starting from vertex\n and conting number of boundary crossings...\n", ntracks);
65 printf("====================================================================\n");
68 Double_t point[3], dir[3];
72 for (i=0; i<ntracks; i++) {
73 phi = 2.*TMath::Pi()*gRandom->Rndm();
74 theta= TMath::ACos(1.-2.*gRandom->Rndm());
75 point[0] = point[1] = point[2] = 0;
76 dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
77 dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
78 dir[2]=TMath::Cos(theta);
79 if ((i%1000)==0) printf("... remaining tracks %i\n", ntracks-i);
80 nbound += propagate_in_geom(point,dir);
82 Double_t time1 = timer->CpuTime() *1.E6;
83 Double_t time2 = time1/ntracks;
84 Double_t time3 = time1/nbound;
85 printf("Time for crossing %i boundaries: %g [ms]\n", nbound, time1);
86 printf("Time per track for full geometry traversal: %g [ms], per crossing: %g [ms]\n", time2, time3);
88 // STAGE 4: How much time per volume:
90 printf("====================================================================\n");
91 printf("STAGE 4: How much navigation time per volume per next+safety call\n");
92 printf("====================================================================\n");
93 TGeoIterator next(gGeoManager->GetTopVolume());
96 vol = gGeoManager->GetTopVolume();
97 memset(flags, 0, nuid*sizeof(Bool_t));
98 score(vol, 1, timing_per_volume(vol));
99 while ((current=next())) {
100 vol = current->GetVolume();
101 Int_t uid = vol->GetNumber();
102 if (flags[uid]) continue;
105 gGeoManager->cd(path.Data());
106 score(vol,1,timing_per_volume(vol));
110 Double_t time_tot_pertrack = 0.;
111 TCanvas *c1 = new TCanvas("c2","ncrossings",10,10,900,500);
113 c1->SetTopMargin(0.15);
114 TFile *f = new TFile("histos.root", "RECREATE");
115 TH1F *h = new TH1F("h","number of boundary crossings per volume",3,0,3);
118 h->SetBit(TH1::kCanRebin);
120 memset(flags, 0, nuid*sizeof(Bool_t));
121 for (i=0; i<nuid; i++) {
122 vol = gGeoManager->GetVolume(i);
123 if (!vol->GetNdaughters()) continue;
124 time_tot_pertrack += val1[i]*val2[i];
125 h->Fill(vol->GetName(), (Int_t)val1[i]);
127 time_tot_pertrack /= ntracks;
129 h->LabelsOption(">","X");
133 TCanvas *c2 = new TCanvas("c3","time spent per volume in navigation",10,10,900,500);
135 c2->SetTopMargin(0.15);
136 TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0,100,100,0,15);
138 h2->SetMarkerStyle(2);
139 TH1F *h1 = new TH1F("h1","percent of time spent per volume",3,0,3);
141 h1->SetFillColor(38);
142 h1->SetBit(TH1::kCanRebin);
143 for (i=0; i<nuid; i++) {
144 vol = gGeoManager->GetVolume(i);
145 if (!vol->GetNdaughters()) continue;
146 value = val1[i]*val2[i]/ntracks/time_tot_pertrack;
147 h1->Fill(vol->GetName(), value);
148 h2->Fill(vol->GetNdaughters(), val2[i]);
151 h1->LabelsOption(">","X");
153 TCanvas *c3 = new TCanvas("c4","timing vs. ndaughters",10,10,900,500);
155 c3->SetTopMargin(0.15);
163 Int_t propagate_in_geom(Double_t *start, Double_t *dir)
165 // Propagate from START along DIR from boundary to boundary until exiting
166 // geometry. Fill array of hits.
167 gGeoManager->InitTrack(start, dir);
168 TGeoNode *current = 0;
171 while (!gGeoManager->IsOutside()) {
172 current = gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big(), kFALSE);
173 if (!current || gGeoManager->IsOutside()) return nhits;
174 Double_t step = gGeoManager->GetStep();
175 if (step<2.*TGeoShape::Tolerance()) {
181 // Problems in trying to cross a boundary
182 printf("Error in trying to cross boundary of %s\n", current->GetName());
187 TGeoVolume *vol = current->GetVolume();
190 TGeoNode *mother = gGeoManager->GetMother(iup++);
191 while (mother && mother->GetVolume()->IsAssembly()) {
192 score(mother->GetVolume(), 0, 1.);
193 mother = gGeoManager->GetMother(iup++);
199 void score(TGeoVolume *vol, Int_t ifield, Double_t value)
201 // Score something for VOL
202 Int_t uid = vol->GetNumber();
212 Double_t timing_per_volume(TGeoVolume *vol)
214 // Compute timing per "FindNextBoundary" + "Safety" call. Volume must be
215 // in the current path.
217 const TGeoShape *shape = vol->GetShape();
218 TGeoBBox *box = (TGeoBBox *)shape;
219 Double_t dx = box->GetDX();
220 Double_t dy = box->GetDY();
221 Double_t dz = box->GetDZ();
222 Double_t ox = (box->GetOrigin())[0];
223 Double_t oy = (box->GetOrigin())[1];
224 Double_t oz = (box->GetOrigin())[2];
225 Double_t point[3], dir[3], lpt[3], ldir[3];
227 pstep = TMath::Max(pstep,dz);
233 for (Int_t i=0; i<1000000; i++) {
234 lpt[0] = ox-dx+2*dx*gRandom->Rndm();
235 lpt[1] = oy-dy+2*dy*gRandom->Rndm();
236 lpt[2] = oz-dz+2*dz*gRandom->Rndm();
237 gGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
238 gGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
239 phi = 2*TMath::Pi()*gRandom->Rndm();
240 theta= TMath::ACos(1.-2.*gRandom->Rndm());
241 ldir[0]=TMath::Sin(theta)*TMath::Cos(phi);
242 ldir[1]=TMath::Sin(theta)*TMath::Sin(phi);
243 ldir[2]=TMath::Cos(theta);
244 gGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir,dir);
245 gGeoManager->SetCurrentDirection(dir);
246 gGeoManager->SetStep(pstep);
247 gGeoManager->ResetState();
249 dist = TGeoShape::Big();
250 if (!vol->IsAssembly()) {
251 inside = vol->Contains(lpt);
253 dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep);
254 // if (dist>=pstep) continue;
256 vol->GetShape()->DistFromInside(lpt,ldir,3,pstep);
259 if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside);
261 if (vol->GetNdaughters()) {
262 gGeoManager->Safety();
263 gGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE);
267 Double_t time_per_track = timer->CpuTime();
268 Int_t uid = vol->GetNumber();
269 Int_t ncrossings = (Int_t)val1[uid];
270 if (!vol->GetNdaughters())
271 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);
273 printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(), time_per_track, vol->GetNdaughters(), ncrossings);
274 return time_per_track;