]>
Commit | Line | Data |
---|---|---|
02e991c3 | 1 | #include "TGeoManager.h" |
2 | #include "TGeoVolume.h" | |
3 | #include "TGeoNode.h" | |
4 | #include "TGeoBBox.h" | |
5 | #include "TRandom3.h" | |
6 | #include "TMath.h" | |
7 | #include "TH1.h" | |
8 | #include "TH2.h" | |
9 | #include "TCanvas.h" | |
10 | #include "TStopwatch.h" | |
11 | #include "TFile.h" | |
12 | ||
13 | TStopwatch *timer = 0; | |
14 | ||
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 *); | |
18 | Double_t *val1 = 0; | |
19 | Double_t *val2 = 0; | |
20 | Bool_t *flags = 0; | |
21 | ||
22 | void checkgeom() | |
23 | { | |
24 | TGeoManager::Import("geometry.root"); | |
25 | Int_t nvol = gGeoManager->GetListOfVolumes()->GetEntries(); | |
26 | Int_t nuid = gGeoManager->GetListOfUVolumes()->GetEntries(); | |
27 | timer = new TStopwatch(); | |
28 | Int_t i; | |
29 | Double_t value; | |
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)); | |
36 | TGeoVolume *vol; | |
37 | ||
38 | // STAGE 1: Overlap checking by sampling per volume | |
39 | ||
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; | |
47 | flags[uid] = kTRUE; | |
48 | if (!vol->GetNdaughters()) continue; | |
49 | vol->CheckOverlaps(0.01, "s"); | |
50 | } | |
51 | ||
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); | |
57 | ||
58 | // STAGE 3: How many crossings per volume in a realistic case ? | |
59 | // Ignore volumes with no daughters | |
60 | ||
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"); | |
66 | Int_t nbound = 0; | |
67 | Double_t theta, phi; | |
68 | Double_t point[3], dir[3]; | |
69 | new TRandom3(); | |
70 | ||
71 | timer->Start(); | |
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); | |
81 | } | |
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); | |
87 | ||
88 | // STAGE 4: How much time per volume: | |
89 | ||
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()); | |
94 | TGeoNode*current; | |
95 | TString path; | |
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; | |
103 | flags[uid] = kTRUE; | |
104 | next.GetPath(path); | |
105 | gGeoManager->cd(path.Data()); | |
106 | score(vol,1,timing_per_volume(vol)); | |
107 | } | |
108 | ||
109 | // Draw some histos | |
110 | Double_t time_tot_pertrack = 0.; | |
111 | TCanvas *c1 = new TCanvas("c2","ncrossings",10,10,900,500); | |
112 | c1->SetGrid(); | |
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); | |
116 | h->SetStats(0); | |
117 | h->SetFillColor(38); | |
118 | h->SetBit(TH1::kCanRebin); | |
119 | ||
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]); | |
126 | } | |
127 | time_tot_pertrack /= ntracks; | |
128 | h->LabelsDeflate(); | |
129 | h->LabelsOption(">","X"); | |
130 | h->Draw(); | |
131 | ||
132 | ||
133 | TCanvas *c2 = new TCanvas("c3","time spent per volume in navigation",10,10,900,500); | |
134 | c2->SetGrid(); | |
135 | c2->SetTopMargin(0.15); | |
136 | TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0,100,100,0,15); | |
137 | h2->SetStats(0); | |
138 | h2->SetMarkerStyle(2); | |
139 | TH1F *h1 = new TH1F("h1","percent of time spent per volume",3,0,3); | |
140 | h1->SetStats(0); | |
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]); | |
149 | } | |
150 | h1->LabelsDeflate(); | |
151 | h1->LabelsOption(">","X"); | |
152 | h1->Draw(); | |
153 | TCanvas *c3 = new TCanvas("c4","timing vs. ndaughters",10,10,900,500); | |
154 | c3->SetGrid(); | |
155 | c3->SetTopMargin(0.15); | |
156 | h2->Draw(); | |
157 | f->Write(); | |
158 | delete [] flags; | |
159 | delete [] val1; | |
160 | delete [] val2; | |
161 | } | |
162 | ||
163 | Int_t propagate_in_geom(Double_t *start, Double_t *dir) | |
164 | { | |
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; | |
169 | Int_t nzero = 0; | |
170 | Int_t nhits = 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()) { | |
176 | nzero++; | |
177 | continue; | |
178 | } | |
179 | else nzero = 0; | |
180 | if (nzero>3) { | |
181 | // Problems in trying to cross a boundary | |
182 | printf("Error in trying to cross boundary of %s\n", current->GetName()); | |
183 | return nhits; | |
184 | } | |
185 | // Generate the hit | |
186 | nhits++; | |
187 | TGeoVolume *vol = current->GetVolume(); | |
188 | score(vol,0,1.); | |
189 | Int_t iup = 1; | |
190 | TGeoNode *mother = gGeoManager->GetMother(iup++); | |
191 | while (mother && mother->GetVolume()->IsAssembly()) { | |
192 | score(mother->GetVolume(), 0, 1.); | |
193 | mother = gGeoManager->GetMother(iup++); | |
194 | } | |
195 | } | |
196 | return nhits; | |
197 | } | |
198 | ||
199 | void score(TGeoVolume *vol, Int_t ifield, Double_t value) | |
200 | { | |
201 | // Score something for VOL | |
202 | Int_t uid = vol->GetNumber(); | |
203 | switch (ifield) { | |
204 | case 0: | |
205 | val1[uid] += value; | |
206 | break; | |
207 | case 1: | |
208 | val2[uid] += value; | |
209 | } | |
210 | } | |
211 | ||
212 | Double_t timing_per_volume(TGeoVolume *vol) | |
213 | { | |
214 | // Compute timing per "FindNextBoundary" + "Safety" call. Volume must be | |
215 | // in the current path. | |
216 | timer->Reset(); | |
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]; | |
226 | Double_t pstep = 0.; | |
227 | pstep = TMath::Max(pstep,dz); | |
228 | Double_t theta, phi; | |
229 | Int_t idaughter; | |
230 | timer->Start(); | |
231 | Double_t dist; | |
232 | Bool_t inside; | |
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(); | |
248 | inside = kTRUE; | |
249 | dist = TGeoShape::Big(); | |
250 | if (!vol->IsAssembly()) { | |
251 | inside = vol->Contains(lpt); | |
252 | if (!inside) { | |
253 | dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep); | |
254 | // if (dist>=pstep) continue; | |
255 | } else { | |
256 | vol->GetShape()->DistFromInside(lpt,ldir,3,pstep); | |
257 | } | |
258 | ||
259 | if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside); | |
260 | } | |
261 | if (vol->GetNdaughters()) { | |
262 | gGeoManager->Safety(); | |
263 | gGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE); | |
264 | } | |
265 | } | |
266 | timer->Stop(); | |
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); | |
272 | else | |
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; | |
275 | } | |
276 |