STEER/STEER: Install AliRawDataHeaderSim.h for HMPID AMORE module
[u/mrichter/AliRoot.git] / macros / checkgeom.C
CommitLineData
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
13TStopwatch *timer = 0;
14
15Int_t propagate_in_geom(Double_t *, Double_t *);
16void score(TGeoVolume *, Int_t, Double_t);
17Double_t timing_per_volume(TGeoVolume *);
18Double_t *val1 = 0;
19Double_t *val2 = 0;
20Bool_t *flags = 0;
21
22void 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
163Int_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
199void 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
212Double_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