[u/mrichter/AliRoot.git] / macros / checkgeom.C

#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TGeoNode.h"
#include "TGeoBBox.h"
#include "TRandom3.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "TFile.h"

TStopwatch *timer = 0;

Int_t propagate_in_geom(Double_t *, Double_t *);
void score(TGeoVolume *, Int_t, Double_t);
Double_t timing_per_volume(TGeoVolume *);
Double_t *val1 = 0;
Double_t *val2 = 0;
Bool_t   *flags = 0;

void checkgeom()
{
   TGeoManager::Import("geometry.root");
   Int_t nvol = gGeoManager->GetListOfVolumes()->GetEntries();
   Int_t nuid = gGeoManager->GetListOfUVolumes()->GetEntries();
   timer = new TStopwatch();
   Int_t i;
   Double_t value;
   flags = new Bool_t[nuid];
   val1 = new Double_t[nuid];
   val2 = new Double_t[nuid];
   memset(flags, 0, nuid*sizeof(Bool_t));
   memset(val1, 0, nuid*sizeof(Double_t));
   memset(val2, 0, nuid*sizeof(Double_t));
   TGeoVolume *vol;

// STAGE 1: Overlap checking by sampling per volume

   printf("====================================================================\n");
   printf("STAGE 1: Overlap checking by sampling per volume\n");
   printf("====================================================================\n");
   for (i=0; i<nvol; i++) {
      vol = (TGeoVolume*)gGeoManager->GetListOfVolumes()->At(i);
      Int_t uid = vol->GetNumber();
      if (flags[uid]) continue;
      flags[uid] = kTRUE;
      if (!vol->GetNdaughters()) continue;
      vol->CheckOverlaps(0.01, "s"); 
   }

// STAGE 2: Global overlap checking
   printf("====================================================================\n");
   printf("STAGE 2: Global overlap checking\n");
   printf("====================================================================\n");
   gGeoManager->CheckOverlaps(0.01);
   
// STAGE 3: How many crossings per volume in a realistic case ?
   // Ignore volumes with no daughters

   // Generate rays from vertex in phi=[0,2*pi] theta=[0,pi]
   Int_t ntracks = 1000000;
   printf("====================================================================\n");
   printf("STAGE 3: Propagating %i tracks starting from vertex\n and conting number of boundary crossings...\n", ntracks);
   printf("====================================================================\n");
   Int_t nbound = 0;
   Double_t theta, phi;
   Double_t point[3], dir[3];
   new TRandom3();

   timer->Start();
   for (i=0; i<ntracks; i++) {
      phi = 2.*TMath::Pi()*gRandom->Rndm();
      theta= TMath::ACos(1.-2.*gRandom->Rndm());
      point[0] = point[1] = point[2] = 0;
      dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
      dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
      dir[2]=TMath::Cos(theta);
      if ((i%1000)==0) printf("... remaining tracks %i\n", ntracks-i);
      nbound += propagate_in_geom(point,dir);
   }
   Double_t time1 = timer->CpuTime() *1.E6;
   Double_t time2 = time1/ntracks;
   Double_t time3 = time1/nbound;
   printf("Time for crossing %i boundaries: %g [ms]\n", nbound, time1);
   printf("Time per track for full geometry traversal: %g [ms], per crossing: %g [ms]\n", time2, time3);

// STAGE 4: How much time per volume:
   
   printf("====================================================================\n");
   printf("STAGE 4: How much navigation time per volume per next+safety call\n");
   printf("====================================================================\n");
   TGeoIterator next(gGeoManager->GetTopVolume());
   TGeoNode*current;
   TString path;
   vol = gGeoManager->GetTopVolume();
   memset(flags, 0, nuid*sizeof(Bool_t));
   score(vol, 1, timing_per_volume(vol)); 
   while ((current=next())) {
      vol = current->GetVolume();
      Int_t uid = vol->GetNumber();
      if (flags[uid]) continue;
      flags[uid] = kTRUE;
      next.GetPath(path);
      gGeoManager->cd(path.Data());
      score(vol,1,timing_per_volume(vol));
   }   

   // Draw some histos
   Double_t time_tot_pertrack = 0.;
   TCanvas *c1 = new TCanvas("c2","ncrossings",10,10,900,500);
   c1->SetGrid();
   c1->SetTopMargin(0.15);
   TFile *f = new TFile("histos.root", "RECREATE");
   TH1F *h = new TH1F("h","number of boundary crossings per volume",3,0,3);
   h->SetStats(0);
   h->SetFillColor(38);
   h->SetBit(TH1::kCanRebin);
   
   memset(flags, 0, nuid*sizeof(Bool_t));
   for (i=0; i<nuid; i++) {
      vol = gGeoManager->GetVolume(i);
      if (!vol->GetNdaughters()) continue;
      time_tot_pertrack += val1[i]*val2[i];
      h->Fill(vol->GetName(), (Int_t)val1[i]);
   }
   time_tot_pertrack /= ntracks;
   h->LabelsDeflate();
   h->LabelsOption(">","X");
   h->Draw();   


   TCanvas *c2 = new TCanvas("c3","time spent per volume in navigation",10,10,900,500);
   c2->SetGrid();
   c2->SetTopMargin(0.15);
   TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0,100,100,0,15);
   h2->SetStats(0);
   h2->SetMarkerStyle(2);
   TH1F *h1 = new TH1F("h1","percent of time spent per volume",3,0,3);
   h1->SetStats(0);
   h1->SetFillColor(38);
   h1->SetBit(TH1::kCanRebin);
   for (i=0; i<nuid; i++) {
      vol = gGeoManager->GetVolume(i);
      if (!vol->GetNdaughters()) continue;
      value = val1[i]*val2[i]/ntracks/time_tot_pertrack;
      h1->Fill(vol->GetName(), value);
      h2->Fill(vol->GetNdaughters(), val2[i]);
   }     
   h1->LabelsDeflate();
   h1->LabelsOption(">","X");
   h1->Draw();
   TCanvas *c3 = new TCanvas("c4","timing vs. ndaughters",10,10,900,500);
   c3->SetGrid();
   c3->SetTopMargin(0.15);
   h2->Draw();   
   f->Write();
   delete [] flags;
   delete [] val1;
   delete [] val2;
}

Int_t propagate_in_geom(Double_t *start, Double_t *dir)
{
// Propagate from START along DIR from boundary to boundary until exiting 
// geometry. Fill array of hits.
   gGeoManager->InitTrack(start, dir);
   TGeoNode *current = 0;
   Int_t nzero = 0;
   Int_t nhits = 0;
   while (!gGeoManager->IsOutside()) {
      current = gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big(), kFALSE);
      if (!current || gGeoManager->IsOutside()) return nhits;
      Double_t step = gGeoManager->GetStep();
      if (step<2.*TGeoShape::Tolerance()) {
         nzero++;
         continue;
      } 
      else nzero = 0;
      if (nzero>3) {
      // Problems in trying to cross a boundary
         printf("Error in trying to cross boundary of %s\n", current->GetName());
         return nhits;
      }
      // Generate the hit
      nhits++;
      TGeoVolume *vol = current->GetVolume();
      score(vol,0,1.);
      Int_t iup = 1;
      TGeoNode *mother = gGeoManager->GetMother(iup++);
      while (mother && mother->GetVolume()->IsAssembly()) {
         score(mother->GetVolume(), 0, 1.);
         mother = gGeoManager->GetMother(iup++);
      }   
   }
   return nhits;
}      

void score(TGeoVolume *vol, Int_t ifield, Double_t value)
{
// Score something for VOL
   Int_t uid = vol->GetNumber();
   switch (ifield) {
      case 0:
         val1[uid] += value;
         break;
      case 1:
         val2[uid] += value;
   }
}   

Double_t timing_per_volume(TGeoVolume *vol)
{
// Compute timing per "FindNextBoundary" + "Safety" call. Volume must be
// in the current path.
   timer->Reset();
   const TGeoShape *shape = vol->GetShape();
   TGeoBBox *box = (TGeoBBox *)shape;
   Double_t dx = box->GetDX();
   Double_t dy = box->GetDY();
   Double_t dz = box->GetDZ();
   Double_t ox = (box->GetOrigin())[0];
   Double_t oy = (box->GetOrigin())[1];
   Double_t oz = (box->GetOrigin())[2];
   Double_t point[3], dir[3], lpt[3], ldir[3];
   Double_t pstep = 0.;
   pstep = TMath::Max(pstep,dz);
   Double_t theta, phi;
   Int_t idaughter;
   timer->Start();
   Double_t dist;
   Bool_t inside;
   for (Int_t i=0; i<1000000; i++) {
      lpt[0] = ox-dx+2*dx*gRandom->Rndm();
      lpt[1] = oy-dy+2*dy*gRandom->Rndm();
      lpt[2] = oz-dz+2*dz*gRandom->Rndm();
      gGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
      gGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
      phi = 2*TMath::Pi()*gRandom->Rndm();
      theta= TMath::ACos(1.-2.*gRandom->Rndm());
      ldir[0]=TMath::Sin(theta)*TMath::Cos(phi);
      ldir[1]=TMath::Sin(theta)*TMath::Sin(phi);
      ldir[2]=TMath::Cos(theta);
      gGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir,dir);
      gGeoManager->SetCurrentDirection(dir);
      gGeoManager->SetStep(pstep);
      gGeoManager->ResetState();
      inside = kTRUE;
      dist = TGeoShape::Big();
      if (!vol->IsAssembly()) {
         inside = vol->Contains(lpt);
         if (!inside) {
            dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep); 
//            if (dist>=pstep) continue;
         } else {   
            vol->GetShape()->DistFromInside(lpt,ldir,3,pstep);
         }   
            
         if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside);
      }   
      if (vol->GetNdaughters()) {
         gGeoManager->Safety();
         gGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE);
      }   
   }
   timer->Stop();
   Double_t time_per_track = timer->CpuTime();
   Int_t uid = vol->GetNumber();
   Int_t ncrossings = (Int_t)val1[uid];
   if (!vol->GetNdaughters())
      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);
   else
      printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(), time_per_track, vol->GetNdaughters(), ncrossings);
   return time_per_track;
}
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+2dxgRandom->Rndm();
235	lpt[1] = oy-dy+2dygRandom->Rndm();
236	lpt[2] = oz-dz+2dzgRandom->Rndm();
237	gGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
238	gGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
239	phi = 2TMath::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