Updated dependencies between the libraries
[u/mrichter/AliRoot.git] / macros / checkgeom.C
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