Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoChecker.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 01/11/01
3 // CheckGeometry(), CheckOverlaps() by Mihaela Gheata
4 
5 /*************************************************************************
6  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7  * All rights reserved. *
8  * *
9  * For the licensing terms see $ROOTSYS/LICENSE. *
10  * For the list of contributors see $ROOTSYS/README/CREDITS. *
11  *************************************************************************/
12 
13 /** \class TGeoChecker
14 \ingroup Geometry_classes
15 
16 Geometry checking package.
17 
18 TGeoChecker class provides several geometry checking methods. There are two
19 types of tests that can be performed. One is based on random sampling or
20 ray-tracing and provides a visual check on how navigation methods work for
21 a given geometry. The second actually checks the validity of the geometry
22 definition in terms of overlapping/extruding objects. Both types of checks
23 can be done for a given branch (starting with a given volume) as well as for
24 the geometry as a whole.
25 
26 #### TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z)
27 
28 This method can be called directly from the TGeoManager class and print a
29 report on how a given point is classified by the modeller: which is the
30 full path to the deepest node containing it, and the (under)estimation
31 of the distance to the closest boundary (safety).
32 
33 #### TGeoChecker::RandomPoints(Int_t npoints)
34 
35 Can be called from TGeoVolume class. It first draws the volume and its
36 content with the current visualization settings. Then randomly samples points
37 in its bounding box, plotting in the geometry display only the points
38 classified as belonging to visible volumes.
39 
40 #### TGeoChecker::RandomRays(Int_t nrays, Double_t startx, starty, startz)
41 
42 Can be called and acts in the same way as the previous, but instead of points,
43 rays having random isotropic directions are generated from the given point.
44 A raytracing algorithm propagates all rays until they exit geometry, plotting
45 all segments crossing visible nodes in the same color as these.
46 
47 #### TGeoChecker::Test(Int_t npoints)
48 
49 Implementation of TGeoManager::Test(). Computes the time for the modeller
50 to find out "Where am I?" for a given number of random points.
51 
52 #### TGeoChecker::LegoPlot(ntheta, themin, themax, nphi, phimin, phimax,...)
53 
54 Implementation of TGeoVolume::LegoPlot(). Draws a spherical radiation length
55 lego plot for a given volume, in a given theta/phi range.
56 
57 #### TGeoChecker::Weight(Double_t precision)
58 
59 Implementation of TGeoVolume::Weight(). Estimates the total weight of a given
60 volume by material sampling. Accepts as input the desired precision.
61 */
62 
63 #include "TVirtualPad.h"
64 #include "TCanvas.h"
65 #include "TStyle.h"
66 #include "TFile.h"
67 #include "TNtuple.h"
68 #include "TH2.h"
69 #include "TRandom3.h"
70 #include "TPolyMarker3D.h"
71 #include "TPolyLine3D.h"
72 #include "TStopwatch.h"
73 
74 #include "TGeoVoxelFinder.h"
75 #include "TGeoBBox.h"
76 #include "TGeoPcon.h"
77 #include "TGeoManager.h"
78 #include "TGeoOverlap.h"
79 #include "TGeoPainter.h"
80 #include "TGeoChecker.h"
81 
82 #include "TBuffer3D.h"
83 #include "TBuffer3DTypes.h"
84 #include "TMath.h"
85 
86 #include <stdlib.h>
87 
88 // statics and globals
89 
90 ClassImp(TGeoChecker);
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// Default constructor
94 
95 TGeoChecker::TGeoChecker()
96  :TObject(),
97  fGeoManager(NULL),
98  fVsafe(NULL),
99  fBuff1(NULL),
100  fBuff2(NULL),
101  fFullCheck(kFALSE),
102  fVal1(NULL),
103  fVal2(NULL),
104  fFlags(NULL),
105  fTimer(NULL),
106  fSelectedNode(NULL),
107  fNchecks(0),
108  fNmeshPoints(1000)
109 {
110 }
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 /// Constructor for a given geometry
114 
115 TGeoChecker::TGeoChecker(TGeoManager *geom)
116  :TObject(),
117  fGeoManager(geom),
118  fVsafe(NULL),
119  fBuff1(NULL),
120  fBuff2(NULL),
121  fFullCheck(kFALSE),
122  fVal1(NULL),
123  fVal2(NULL),
124  fFlags(NULL),
125  fTimer(NULL),
126  fSelectedNode(NULL),
127  fNchecks(0),
128  fNmeshPoints(1000)
129 {
130  fBuff1 = new TBuffer3D(TBuffer3DTypes::kGeneric,500,3*500,0,0,0,0);
131  fBuff2 = new TBuffer3D(TBuffer3DTypes::kGeneric,500,3*500,0,0,0,0);
132 }
133 
134 ////////////////////////////////////////////////////////////////////////////////
135 /// Destructor
136 
137 TGeoChecker::~TGeoChecker()
138 {
139  if (fBuff1) delete fBuff1;
140  if (fBuff2) delete fBuff2;
141  if (fTimer) delete fTimer;
142 }
143 
144 ////////////////////////////////////////////////////////////////////////////////
145 /// Print current operation progress.
146 
147 void TGeoChecker::OpProgress(const char *opname, Long64_t current, Long64_t size, TStopwatch *watch, Bool_t last, Bool_t refresh, const char *msg)
148 {
149  static Long64_t icount = 0;
150  static TString oname;
151  static TString nname;
152  static Long64_t ocurrent = 0;
153  static Long64_t osize = 0;
154  static Int_t oseconds = 0;
155  static TStopwatch *owatch = 0;
156  static Bool_t oneoftwo = kFALSE;
157  static Int_t nrefresh = 0;
158  const char symbol[4] = {'=','\\','|','/'};
159  char progress[11] = " ";
160  Int_t ichar = icount%4;
161  TString message(msg);
162  message += " ";
163 
164  if (!refresh) {
165  nrefresh = 0;
166  if (!size) return;
167  owatch = watch;
168  oname = opname;
169  ocurrent = TMath::Abs(current);
170  osize = TMath::Abs(size);
171  if (ocurrent > osize) ocurrent=osize;
172  } else {
173  nrefresh++;
174  if (!osize) return;
175  }
176  icount++;
177  Double_t time = 0.;
178  Int_t hours = 0;
179  Int_t minutes = 0;
180  Int_t seconds = 0;
181  if (owatch && !last) {
182  owatch->Stop();
183  time = owatch->RealTime();
184  hours = (Int_t)(time/3600.);
185  time -= 3600*hours;
186  minutes = (Int_t)(time/60.);
187  time -= 60*minutes;
188  seconds = (Int_t)time;
189  if (refresh) {
190  if (oseconds==seconds) {
191  owatch->Continue();
192  return;
193  }
194  oneoftwo = !oneoftwo;
195  }
196  oseconds = seconds;
197  }
198  if (refresh && oneoftwo) {
199  nname = oname;
200  if (fNchecks <= nrefresh) fNchecks = nrefresh+1;
201  Int_t pctdone = (Int_t)(100.*nrefresh/fNchecks);
202  oname = TString::Format(" == %3d%% ==", pctdone);
203  }
204  Double_t percent = 100.0*ocurrent/osize;
205  Int_t nchar = Int_t(percent/10);
206  if (nchar>10) nchar=10;
207  Int_t i;
208  for (i=0; i<nchar; i++) progress[i] = '=';
209  progress[nchar] = symbol[ichar];
210  for (i=nchar+1; i<10; i++) progress[i] = ' ';
211  progress[10] = '\0';
212  oname += " ";
213  oname.Remove(20);
214  if(size<10000) fprintf(stderr, "%s [%10s] %4lld ", oname.Data(), progress, ocurrent);
215  else if(size<100000) fprintf(stderr, "%s [%10s] %5lld ",oname.Data(), progress, ocurrent);
216  else fprintf(stderr, "%s [%10s] %7lld ",oname.Data(), progress, ocurrent);
217  if (time>0.) fprintf(stderr, "[%6.2f %%] TIME %.2d:%.2d:%.2d %s\r", percent, hours, minutes, seconds, message.Data());
218  else fprintf(stderr, "[%6.2f %%] %s\r", percent, message.Data());
219  if (refresh && oneoftwo) oname = nname;
220  if (owatch) owatch->Continue();
221  if (last) {
222  icount = 0;
223  owatch = 0;
224  ocurrent = 0;
225  osize = 0;
226  oseconds = 0;
227  oneoftwo = kFALSE;
228  nrefresh = 0;
229  fprintf(stderr, "\n");
230  }
231 }
232 
233 ////////////////////////////////////////////////////////////////////////////////
234 /// Check pushes and pulls needed to cross the next boundary with respect to the
235 /// position given by FindNextBoundary. If radius is not mentioned the full bounding
236 /// box will be sampled.
237 
238 void TGeoChecker::CheckBoundaryErrors(Int_t ntracks, Double_t radius)
239 {
240  TGeoVolume *tvol = fGeoManager->GetTopVolume();
241  Info("CheckBoundaryErrors", "Top volume is %s",tvol->GetName());
242  const TGeoShape *shape = tvol->GetShape();
243  TGeoBBox *box = (TGeoBBox *)shape;
244  Double_t dl[3];
245  Double_t ori[3];
246  Double_t xyz[3];
247  Double_t nxyz[3];
248  Double_t dir[3];
249  Double_t relp;
250  Char_t path[1024];
251  Char_t cdir[10];
252 
253  // Tree part
254  TFile *f=new TFile("geobugs.root","recreate");
255  TTree *bug=new TTree("bug","Geometrical problems");
256  bug->Branch("pos",xyz,"xyz[3]/D");
257  bug->Branch("dir",dir,"dir[3]/D");
258  bug->Branch("push",&relp,"push/D");
259  bug->Branch("path",&path,"path/C");
260  bug->Branch("cdir",&cdir,"cdir/C");
261 
262  dl[0] = box->GetDX();
263  dl[1] = box->GetDY();
264  dl[2] = box->GetDZ();
265  ori[0] = (box->GetOrigin())[0];
266  ori[1] = (box->GetOrigin())[1];
267  ori[2] = (box->GetOrigin())[2];
268  if (radius>0)
269  dl[0] = dl[1] = dl[2] = radius;
270 
271  TH1::AddDirectory(kFALSE);
272  TH1F *hnew = new TH1F("hnew","Precision pushing",30,-20.,10.);
273  TH1F *hold = new TH1F("hold","Precision pulling", 30,-20.,10.);
274  TH2F *hplotS = new TH2F("hplotS","Problematic points",100,-dl[0],dl[0],100,-dl[1],dl[1]);
275  gStyle->SetOptStat(111111);
276 
277  TGeoNode *node = 0;
278  Long_t igen=0;
279  Long_t itry=0;
280  Long_t n100 = ntracks/100;
281  Double_t rad = TMath::Sqrt(dl[0]*dl[0]+dl[1]*dl[1]);
282  printf("Random box : %f, %f, %f, %f, %f, %f\n", ori[0], ori[1], ori[2], dl[0], dl[1], dl[2]);
283  printf("Start... %i points\n", ntracks);
284  if (!fTimer) fTimer = new TStopwatch();
285  fTimer->Reset();
286  fTimer->Start();
287  while (igen<ntracks) {
288  Double_t phi1 = TMath::TwoPi()*gRandom->Rndm();
289  Double_t r = rad*gRandom->Rndm();
290  xyz[0] = ori[0] + r*TMath::Cos(phi1);
291  xyz[1] = ori[1] + r*TMath::Sin(phi1);
292  Double_t z = (1.-2.*gRandom->Rndm());
293  xyz[2] = ori[2]+dl[2]*z*TMath::Abs(z);
294  ++itry;
295  fGeoManager->SetCurrentPoint(xyz);
296  node = fGeoManager->FindNode();
297  if (!node || node==fGeoManager->GetTopNode()) continue;
298  ++igen;
299  if (n100 && !(igen%n100))
300  OpProgress("Sampling progress:",igen, ntracks, fTimer);
301  Double_t cost = 1.-2.*gRandom->Rndm();
302  Double_t sint = TMath::Sqrt((1.+cost)*(1.-cost));
303  Double_t phi = TMath::TwoPi()*gRandom->Rndm();
304  dir[0] = sint * TMath::Cos(phi);
305  dir[1] = sint * TMath::Sin(phi);
306  dir[2] = cost;
307  fGeoManager->SetCurrentDirection(dir);
308  fGeoManager->FindNextBoundary();
309  Double_t step = fGeoManager->GetStep();
310 
311  relp = 1.e-21;
312  for(Int_t i=0; i<30; ++i) {
313  relp *=10.;
314  for(Int_t j=0; j<3; ++j) nxyz[j]=xyz[j]+step*(1.+relp)*dir[j];
315  if(!fGeoManager->IsSameLocation(nxyz[0],nxyz[1],nxyz[2])) {
316  hnew->Fill(i-20.);
317  if(i>15) {
318  const Double_t* norm = fGeoManager->FindNormal();
319  strncpy(path,fGeoManager->GetPath(),1024);
320  path[1023] = '\0';
321  Double_t dotp = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
322  printf("Forward error i=%d p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
323  i,xyz[0],xyz[1],xyz[2],step,dotp,path);
324  hplotS->Fill(xyz[0],xyz[1],(Double_t)i);
325  strncpy(cdir,"Forward",10);
326  bug->Fill();
327  }
328  break;
329  }
330  }
331 
332  relp = -1.e-21;
333  for(Int_t i=0; i<30; ++i) {
334  relp *=10.;
335  for(Int_t j=0; j<3; ++j) nxyz[j]=xyz[j]+step*(1.+relp)*dir[j];
336  if(fGeoManager->IsSameLocation(nxyz[0],nxyz[1],nxyz[2])) {
337  hold->Fill(i-20.);
338  if(i>15) {
339  const Double_t* norm = fGeoManager->FindNormal();
340  strncpy(path,fGeoManager->GetPath(),1024);
341  path[1023] = '\0';
342  Double_t dotp = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
343  printf("Backward error i=%d p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
344  i,xyz[0],xyz[1],xyz[2],step,dotp,path);
345  strncpy(cdir,"Backward",10);
346  bug->Fill();
347  }
348  break;
349  }
350  }
351  }
352  fTimer->Stop();
353 
354  if (itry) printf("CPU time/point = %5.2emus: Real time/point = %5.2emus\n",
355  1000000.*fTimer->CpuTime()/itry,1000000.*fTimer->RealTime()/itry);
356  bug->Write();
357  delete bug;
358  bug=0;
359  delete f;
360 
361  CheckBoundaryReference();
362 
363  if (itry) printf("Effic = %3.1f%%\n",(100.*igen)/itry);
364  TCanvas *c1 = new TCanvas("c1","Results",600,800);
365  c1->Divide(1,2);
366  c1->cd(1);
367  gPad->SetLogy();
368  hold->Draw();
369  c1->cd(2);
370  gPad->SetLogy();
371  hnew->Draw();
372  /*TCanvas *c3 = */new TCanvas("c3","Plot",600,600);
373  hplotS->Draw("cont0");
374 }
375 
376 ////////////////////////////////////////////////////////////////////////////////
377 /// Check the boundary errors reference file created by CheckBoundaryErrors method.
378 /// The shape for which the crossing failed is drawn with the starting point in red
379 /// and the extrapolated point to boundary (+/- failing push/pull) in yellow.
380 
381 void TGeoChecker::CheckBoundaryReference(Int_t icheck)
382 {
383  Double_t xyz[3];
384  Double_t nxyz[3];
385  Double_t dir[3];
386  Double_t lnext[3];
387  Double_t push;
388  Char_t path[1024];
389  Char_t cdir[10];
390  // Tree part
391  TFile *f=new TFile("geobugs.root","read");
392  TTree *bug=(TTree*)f->Get("bug");
393  bug->SetBranchAddress("pos",xyz);
394  bug->SetBranchAddress("dir",dir);
395  bug->SetBranchAddress("push",&push);
396  bug->SetBranchAddress("path",&path);
397  bug->SetBranchAddress("cdir",&cdir);
398 
399  Int_t nentries = (Int_t)bug->GetEntries();
400  printf("nentries %d\n",nentries);
401  if (icheck<0) {
402  for (Int_t i=0;i<nentries;i++) {
403  bug->GetEntry(i);
404  printf("%-9s error push=%g p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
405  cdir,push,xyz[0],xyz[1],xyz[2],1.,1.,path);
406  }
407  } else {
408  if (icheck>=nentries) return;
409  Int_t idebug = TGeoManager::GetVerboseLevel();
410  TGeoManager::SetVerboseLevel(5);
411  bug->GetEntry(icheck);
412  printf("%-9s error push=%g p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n",
413  cdir,push,xyz[0],xyz[1],xyz[2],1.,1.,path);
414  fGeoManager->SetCurrentPoint(xyz);
415  fGeoManager->SetCurrentDirection(dir);
416  fGeoManager->FindNode();
417  TGeoNode *next = fGeoManager->FindNextBoundary();
418  Double_t step = fGeoManager->GetStep();
419  for (Int_t j=0; j<3; j++) nxyz[j]=xyz[j]+step*(1.+0.1*push)*dir[j];
420  Bool_t change = !fGeoManager->IsSameLocation(nxyz[0],nxyz[1],nxyz[2]);
421  printf("step=%g in: %s\n", step, fGeoManager->GetPath());
422  printf(" -> next = %s push=%g change=%d\n", next->GetName(),push, (UInt_t)change);
423  next->GetVolume()->InspectShape();
424  next->GetVolume()->DrawOnly();
425  if (next != fGeoManager->GetCurrentNode()) {
426  Int_t index1 = fGeoManager->GetCurrentVolume()->GetIndex(next);
427  if (index1>=0) fGeoManager->CdDown(index1);
428  }
429  TPolyMarker3D *pm = new TPolyMarker3D();
430  fGeoManager->MasterToLocal(xyz, lnext);
431  pm->SetNextPoint(lnext[0], lnext[1], lnext[2]);
432  pm->SetMarkerStyle(2);
433  pm->SetMarkerSize(0.2);
434  pm->SetMarkerColor(kRed);
435  pm->Draw("SAME");
436  TPolyMarker3D *pm1 = new TPolyMarker3D();
437  for (Int_t j=0; j<3; j++) nxyz[j]=xyz[j]+step*dir[j];
438  fGeoManager->MasterToLocal(nxyz, lnext);
439  pm1->SetNextPoint(lnext[0], lnext[1], lnext[2]);
440  pm1->SetMarkerStyle(2);
441  pm1->SetMarkerSize(0.2);
442  pm1->SetMarkerColor(kYellow);
443  pm1->Draw("SAME");
444  TGeoManager::SetVerboseLevel(idebug);
445  }
446  delete bug;
447  delete f;
448 }
449 
450 ////////////////////////////////////////////////////////////////////////////////
451 /// Geometry checking. Optional overlap checkings (by sampling and by mesh). Optional
452 /// boundary crossing check + timing per volume.
453 ///
454 /// STAGE 1: extensive overlap checking by sampling per volume. Stdout need to be
455 /// checked by user to get report, then TGeoVolume::CheckOverlaps(0.01, "s") can
456 /// be called for the suspicious volumes.
457 ///
458 /// STAGE 2: normal overlap checking using the shapes mesh - fills the list of
459 /// overlaps.
460 ///
461 /// STAGE 3: shooting NRAYS rays from VERTEX and counting the total number of
462 /// crossings per volume (rays propagated from boundary to boundary until
463 /// geometry exit). Timing computed and results stored in a histo.
464 ///
465 /// STAGE 4: shooting 1 mil. random rays inside EACH volume and calling
466 /// FindNextBoundary() + Safety() for each call. The timing is normalized by the
467 /// number of crossings computed at stage 2 and presented as percentage.
468 /// One can get a picture on which are the most "burned" volumes during
469 /// transportation from geometry point of view. Another plot of the timing per
470 /// volume vs. number of daughters is produced.
471 ///
472 /// All histos are saved in the file statistics.root
473 
474 void TGeoChecker::CheckGeometryFull(Bool_t checkoverlaps, Bool_t checkcrossings, Int_t ntracks, const Double_t *vertex)
475 {
476  Int_t nuid = fGeoManager->GetListOfUVolumes()->GetEntries();
477  if (!fTimer) fTimer = new TStopwatch();
478  Int_t i;
479  Double_t value;
480  fFlags = new Bool_t[nuid];
481  memset(fFlags, 0, nuid*sizeof(Bool_t));
482  TGeoVolume *vol;
483  TCanvas *c = new TCanvas("overlaps", "Overlaps by sampling", 800,800);
484 
485 // STAGE 1: Overlap checking by sampling per volume
486  if (checkoverlaps) {
487  printf("====================================================================\n");
488  printf("STAGE 1: Overlap checking by sampling within 10 microns\n");
489  printf("====================================================================\n");
490  fGeoManager->CheckOverlaps(0.001, "s");
491 
492  // STAGE 2: Global overlap/extrusion checking
493  printf("====================================================================\n");
494  printf("STAGE 2: Global overlap/extrusion checking within 10 microns\n");
495  printf("====================================================================\n");
496  fGeoManager->CheckOverlaps(0.001);
497  }
498 
499  if (!checkcrossings) {
500  delete [] fFlags;
501  fFlags = 0;
502  delete c;
503  return;
504  }
505 
506  fVal1 = new Double_t[nuid];
507  fVal2 = new Double_t[nuid];
508  memset(fVal1, 0, nuid*sizeof(Double_t));
509  memset(fVal2, 0, nuid*sizeof(Double_t));
510  // STAGE 3: How many crossings per volume in a realistic case ?
511  // Ignore volumes with no daughters
512 
513  // Generate rays from vertex in phi=[0,2*pi] theta=[0,pi]
514 // Int_t ntracks = 1000000;
515  printf("====================================================================\n");
516  printf("STAGE 3: Propagating %i tracks starting from vertex\n and counting number of boundary crossings...\n", ntracks);
517  printf("====================================================================\n");
518  Int_t nbound = 0;
519  Double_t theta, phi;
520  Double_t point[3], dir[3];
521  memset(point, 0, 3*sizeof(Double_t));
522  if (vertex) memcpy(point, vertex, 3*sizeof(Double_t));
523 
524  fTimer->Reset();
525  fTimer->Start();
526  for (i=0; i<ntracks; i++) {
527  phi = 2.*TMath::Pi()*gRandom->Rndm();
528  theta= TMath::ACos(1.-2.*gRandom->Rndm());
529  dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
530  dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
531  dir[2]=TMath::Cos(theta);
532  if ((i%100)==0) OpProgress("Transporting tracks",i, ntracks, fTimer);
533 // if ((i%1000)==0) printf("... remaining tracks %i\n", ntracks-i);
534  nbound += PropagateInGeom(point,dir);
535  }
536  if (!nbound) {
537  printf("No boundary crossed\n");
538  return;
539  }
540  fTimer->Stop();
541  Double_t time1 = fTimer->CpuTime() *1.E6;
542  Double_t time2 = time1/ntracks;
543  Double_t time3 = time1/nbound;
544  OpProgress("Transporting tracks",ntracks, ntracks, fTimer, kTRUE);
545  printf("Time for crossing %i boundaries: %g [ms]\n", nbound, time1);
546  printf("Time per track for full geometry traversal: %g [ms], per crossing: %g [ms]\n", time2, time3);
547 
548 // STAGE 4: How much time per volume:
549 
550  printf("====================================================================\n");
551  printf("STAGE 4: How much navigation time per volume per next+safety call\n");
552  printf("====================================================================\n");
553  TGeoIterator next(fGeoManager->GetTopVolume());
554  TGeoNode*current;
555  TString path;
556  vol = fGeoManager->GetTopVolume();
557  memset(fFlags, 0, nuid*sizeof(Bool_t));
558  TStopwatch timer;
559  timer.Start();
560  i = 0;
561  char volname[30];
562  strncpy(volname, vol->GetName(),15);
563  volname[15] = '\0';
564  OpProgress(volname,i++, nuid, &timer);
565  Score(vol, 1, TimingPerVolume(vol));
566  while ((current=next())) {
567  vol = current->GetVolume();
568  Int_t uid = vol->GetNumber();
569  if (fFlags[uid]) continue;
570  fFlags[uid] = kTRUE;
571  next.GetPath(path);
572  fGeoManager->cd(path.Data());
573  strncpy(volname, vol->GetName(),15);
574  volname[15] = '\0';
575  OpProgress(volname,i++, nuid, &timer);
576  Score(vol,1,TimingPerVolume(vol));
577  }
578  OpProgress("STAGE 4 completed",i, nuid, &timer, kTRUE);
579  // Draw some histos
580  Double_t time_tot_pertrack = 0.;
581  TCanvas *c1 = new TCanvas("c2","ncrossings",10,10,900,500);
582  c1->SetGrid();
583  c1->SetTopMargin(0.15);
584  TFile *f = new TFile("statistics.root", "RECREATE");
585  TH1F *h = new TH1F("h","number of boundary crossings per volume",3,0,3);
586  h->SetStats(0);
587  h->SetFillColor(38);
588  h->SetCanExtend(TH1::kAllAxes);
589 
590  memset(fFlags, 0, nuid*sizeof(Bool_t));
591  for (i=0; i<nuid; i++) {
592  vol = fGeoManager->GetVolume(i);
593  if (!vol->GetNdaughters()) continue;
594  time_tot_pertrack += fVal1[i]*fVal2[i];
595  h->Fill(vol->GetName(), (Int_t)fVal1[i]);
596  }
597  time_tot_pertrack /= ntracks;
598  h->LabelsDeflate();
599  h->LabelsOption(">","X");
600  h->Draw();
601 
602 
603  TCanvas *c2 = new TCanvas("c3","time spent per volume in navigation",10,10,900,500);
604  c2->SetGrid();
605  c2->SetTopMargin(0.15);
606  TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0,100,100,0,15);
607  h2->SetStats(0);
608  h2->SetMarkerStyle(2);
609  TH1F *h1 = new TH1F("h1","percent of time spent per volume",3,0,3);
610  h1->SetStats(0);
611  h1->SetFillColor(38);
612  h1->SetCanExtend(TH1::kAllAxes);
613  for (i=0; i<nuid; i++) {
614  vol = fGeoManager->GetVolume(i);
615  if (!vol || !vol->GetNdaughters()) continue;
616  value = fVal1[i]*fVal2[i]/ntracks/time_tot_pertrack;
617  h1->Fill(vol->GetName(), value);
618  h2->Fill(vol->GetNdaughters(), fVal2[i]);
619  }
620  h1->LabelsDeflate();
621  h1->LabelsOption(">","X");
622  h1->Draw();
623  TCanvas *c3 = new TCanvas("c4","timing vs. ndaughters",10,10,900,500);
624  c3->SetGrid();
625  c3->SetTopMargin(0.15);
626  h2->Draw();
627  f->Write();
628  delete [] fFlags;
629  fFlags = 0;
630  delete [] fVal1;
631  fVal1 = 0;
632  delete [] fVal2;
633  fVal2 = 0;
634  delete fTimer;
635  fTimer = 0;
636  delete c;
637 }
638 
639 ////////////////////////////////////////////////////////////////////////////////
640 /// Propagate from START along DIR from boundary to boundary until exiting
641 /// geometry. Fill array of hits.
642 
643 Int_t TGeoChecker::PropagateInGeom(Double_t *start, Double_t *dir)
644 {
645  fGeoManager->InitTrack(start, dir);
646  TGeoNode *current = 0;
647  Int_t nzero = 0;
648  Int_t nhits = 0;
649  while (!fGeoManager->IsOutside()) {
650  if (nzero>3) {
651  // Problems in trying to cross a boundary
652  printf("Error in trying to cross boundary of %s\n", current->GetName());
653  return nhits;
654  }
655  current = fGeoManager->FindNextBoundaryAndStep(TGeoShape::Big(), kFALSE);
656  if (!current || fGeoManager->IsOutside()) return nhits;
657  Double_t step = fGeoManager->GetStep();
658  if (step<2.*TGeoShape::Tolerance()) {
659  nzero++;
660  continue;
661  }
662  else nzero = 0;
663  // Generate the hit
664  nhits++;
665  TGeoVolume *vol = current->GetVolume();
666  Score(vol,0,1.);
667  Int_t iup = 1;
668  TGeoNode *mother = fGeoManager->GetMother(iup++);
669  while (mother && mother->GetVolume()->IsAssembly()) {
670  Score(mother->GetVolume(), 0, 1.);
671  mother = fGeoManager->GetMother(iup++);
672  }
673  }
674  return nhits;
675 }
676 
677 ////////////////////////////////////////////////////////////////////////////////
678 /// Score a hit for VOL
679 
680 void TGeoChecker::Score(TGeoVolume *vol, Int_t ifield, Double_t value)
681 {
682  Int_t uid = vol->GetNumber();
683  switch (ifield) {
684  case 0:
685  fVal1[uid] += value;
686  break;
687  case 1:
688  fVal2[uid] += value;
689  }
690 }
691 
692 ////////////////////////////////////////////////////////////////////////////////
693 /// Set number of points to be generated on the shape outline when checking for overlaps.
694 
695 void TGeoChecker::SetNmeshPoints(Int_t npoints) {
696  fNmeshPoints = npoints;
697  if (npoints<1000) {
698  Error("SetNmeshPoints", "Cannot allow less than 1000 points for checking - set to 1000");
699  fNmeshPoints = 1000;
700  }
701 }
702 
703 ////////////////////////////////////////////////////////////////////////////////
704 /// Compute timing per "FindNextBoundary" + "Safety" call. Volume must be
705 /// in the current path.
706 
707 Double_t TGeoChecker::TimingPerVolume(TGeoVolume *vol)
708 {
709  fTimer->Reset();
710  const TGeoShape *shape = vol->GetShape();
711  TGeoBBox *box = (TGeoBBox *)shape;
712  Double_t dx = box->GetDX();
713  Double_t dy = box->GetDY();
714  Double_t dz = box->GetDZ();
715  Double_t ox = (box->GetOrigin())[0];
716  Double_t oy = (box->GetOrigin())[1];
717  Double_t oz = (box->GetOrigin())[2];
718  Double_t point[3], dir[3], lpt[3], ldir[3];
719  Double_t pstep = 0.;
720  pstep = TMath::Max(pstep,dz);
721  Double_t theta, phi;
722  Int_t idaughter;
723  fTimer->Start();
724  Bool_t inside;
725  for (Int_t i=0; i<1000000; i++) {
726  lpt[0] = ox-dx+2*dx*gRandom->Rndm();
727  lpt[1] = oy-dy+2*dy*gRandom->Rndm();
728  lpt[2] = oz-dz+2*dz*gRandom->Rndm();
729  fGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
730  fGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
731  phi = 2*TMath::Pi()*gRandom->Rndm();
732  theta= TMath::ACos(1.-2.*gRandom->Rndm());
733  ldir[0]=TMath::Sin(theta)*TMath::Cos(phi);
734  ldir[1]=TMath::Sin(theta)*TMath::Sin(phi);
735  ldir[2]=TMath::Cos(theta);
736  fGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir,dir);
737  fGeoManager->SetCurrentDirection(dir);
738  fGeoManager->SetStep(pstep);
739  fGeoManager->ResetState();
740  inside = kTRUE;
741  // dist = TGeoShape::Big();
742  if (!vol->IsAssembly()) {
743  inside = vol->Contains(lpt);
744  if (!inside) {
745  // dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep);
746  // if (dist>=pstep) continue;
747  } else {
748  vol->GetShape()->DistFromInside(lpt,ldir,3,pstep);
749  }
750 
751  if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside);
752  }
753  if (vol->GetNdaughters()) {
754  fGeoManager->Safety();
755  fGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE);
756  }
757  }
758  fTimer->Stop();
759  Double_t time_per_track = fTimer->CpuTime();
760  Int_t uid = vol->GetNumber();
761  Int_t ncrossings = (Int_t)fVal1[uid];
762  if (!vol->GetNdaughters())
763  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);
764  else
765  printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(), time_per_track, vol->GetNdaughters(), ncrossings);
766  return time_per_track;
767 }
768 
769 ////////////////////////////////////////////////////////////////////////////////
770 /// Shoot nrays with random directions from starting point (startx, starty, startz)
771 /// in the reference frame of this volume. Track each ray until exiting geometry, then
772 /// shoot backwards from exiting point and compare boundary crossing points.
773 
774 void TGeoChecker::CheckGeometry(Int_t nrays, Double_t startx, Double_t starty, Double_t startz) const
775 {
776  Int_t i, j;
777  Double_t start[3], end[3];
778  Double_t dir[3];
779  Double_t dummy[3];
780  Double_t eps = 0.;
781  Double_t *array1 = new Double_t[3*1000];
782  Double_t *array2 = new Double_t[3*1000];
783  TObjArray *pma = new TObjArray();
784  TPolyMarker3D *pm;
785  pm = new TPolyMarker3D();
786  pm->SetMarkerColor(2); // error > eps
787  pm->SetMarkerStyle(8);
788  pm->SetMarkerSize(0.4);
789  pma->AddAt(pm, 0);
790  pm = new TPolyMarker3D();
791  pm->SetMarkerColor(4); // point not found
792  pm->SetMarkerStyle(8);
793  pm->SetMarkerSize(0.4);
794  pma->AddAt(pm, 1);
795  pm = new TPolyMarker3D();
796  pm->SetMarkerColor(6); // extra point back
797  pm->SetMarkerStyle(8);
798  pm->SetMarkerSize(0.4);
799  pma->AddAt(pm, 2);
800  Int_t nelem1, nelem2;
801  Int_t dim1=1000, dim2=1000;
802  if ((startx==0) && (starty==0) && (startz==0)) eps=1E-3;
803  start[0] = startx+eps;
804  start[1] = starty+eps;
805  start[2] = startz+eps;
806  Int_t n10=nrays/10;
807  Double_t theta,phi;
808  Double_t dw, dwmin, dx, dy, dz;
809  Int_t ist1, ist2, ifound;
810  for (i=0; i<nrays; i++) {
811  if (n10) {
812  if ((i%n10) == 0) printf("%i percent\n", Int_t(100*i/nrays));
813  }
814  phi = 2*TMath::Pi()*gRandom->Rndm();
815  theta= TMath::ACos(1.-2.*gRandom->Rndm());
816  dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
817  dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
818  dir[2]=TMath::Cos(theta);
819  // shoot direct ray
820  nelem1=nelem2=0;
821 // printf("DIRECT %i\n", i);
822  array1 = ShootRay(&start[0], dir[0], dir[1], dir[2], array1, nelem1, dim1);
823  if (!nelem1) continue;
824 // for (j=0; j<nelem1; j++) printf("%i : %f %f %f\n", j, array1[3*j], array1[3*j+1], array1[3*j+2]);
825  memcpy(&end[0], &array1[3*(nelem1-1)], 3*sizeof(Double_t));
826  // shoot ray backwards
827 // printf("BACK %i\n", i);
828  array2 = ShootRay(&end[0], -dir[0], -dir[1], -dir[2], array2, nelem2, dim2, &start[0]);
829  if (!nelem2) {
830  printf("#### NOTHING BACK ###########################\n");
831  for (j=0; j<nelem1; j++) {
832  pm = (TPolyMarker3D*)pma->At(0);
833  pm->SetNextPoint(array1[3*j], array1[3*j+1], array1[3*j+2]);
834  }
835  continue;
836  }
837 // printf("BACKWARDS\n");
838  Int_t k=nelem2>>1;
839  for (j=0; j<k; j++) {
840  memcpy(&dummy[0], &array2[3*j], 3*sizeof(Double_t));
841  memcpy(&array2[3*j], &array2[3*(nelem2-1-j)], 3*sizeof(Double_t));
842  memcpy(&array2[3*(nelem2-1-j)], &dummy[0], 3*sizeof(Double_t));
843  }
844 // for (j=0; j<nelem2; j++) printf("%i : %f ,%f ,%f \n", j, array2[3*j], array2[3*j+1], array2[3*j+2]);
845  if (nelem1!=nelem2) printf("### DIFFERENT SIZES : nelem1=%i nelem2=%i ##########\n", nelem1, nelem2);
846  ist1 = ist2 = 0;
847  // check first match
848 
849  dx = array1[3*ist1]-array2[3*ist2];
850  dy = array1[3*ist1+1]-array2[3*ist2+1];
851  dz = array1[3*ist1+2]-array2[3*ist2+2];
852  dw = dx*dir[0]+dy*dir[1]+dz*dir[2];
853  fGeoManager->SetCurrentPoint(&array1[3*ist1]);
854  fGeoManager->FindNode();
855 // printf("%i : %s (%g, %g, %g)\n", ist1, fGeoManager->GetPath(),
856 // array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2]);
857  if (TMath::Abs(dw)<1E-4) {
858 // printf(" matching %i (%g, %g, %g)\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
859  ist2++;
860  } else {
861  printf("### NOT MATCHING %i f:(%f, %f, %f) b:(%f %f %f) DCLOSE=%f\n", ist2, array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2], array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2],dw);
862  pm = (TPolyMarker3D*)pma->At(0);
863  pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
864  if (dw<0) {
865  // first boundary missed on way back
866  } else {
867  // first boundary different on way back
868  ist2++;
869  }
870  }
871 
872  while ((ist1<nelem1-1) && (ist2<nelem2)) {
873  fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
874  fGeoManager->FindNode();
875 // printf("%i : %s (%g, %g, %g)\n", ist1+1, fGeoManager->GetPath(),
876 // array1[3*ist1+3], array1[3*ist1+4], array1[3*ist1+5]);
877 
878  dx = array1[3*ist1+3]-array1[3*ist1];
879  dy = array1[3*ist1+4]-array1[3*ist1+1];
880  dz = array1[3*ist1+5]-array1[3*ist1+2];
881  // distance to next point
882  dwmin = dx+dir[0]+dy*dir[1]+dz*dir[2];
883  while (ist2<nelem2) {
884  ifound = 0;
885  dx = array2[3*ist2]-array1[3*ist1];
886  dy = array2[3*ist2+1]-array1[3*ist1+1];
887  dz = array2[3*ist2+2]-array1[3*ist1+2];
888  dw = dx+dir[0]+dy*dir[1]+dz*dir[2];
889  if (TMath::Abs(dw-dwmin)<1E-4) {
890  ist1++;
891  ist2++;
892  break;
893  }
894  if (dw<dwmin) {
895  // point found on way back. Check if close enough to ist1+1
896  ifound++;
897  dw = dwmin-dw;
898  if (dw<1E-4) {
899  // point is matching ist1+1
900 // printf(" matching %i (%g, %g, %g) DCLOSE=%g\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2], dw);
901  ist2++;
902  ist1++;
903  break;
904  } else {
905  // extra boundary found on way back
906  fGeoManager->SetCurrentPoint(&array2[3*ist2]);
907  fGeoManager->FindNode();
908  pm = (TPolyMarker3D*)pma->At(2);
909  pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
910  printf("### EXTRA BOUNDARY %i : %s found at DCLOSE=%f\n", ist2, fGeoManager->GetPath(), dw);
911  ist2++;
912  continue;
913  }
914  } else {
915  if (!ifound) {
916  // point ist1+1 not found on way back
917  fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
918  fGeoManager->FindNode();
919  pm = (TPolyMarker3D*)pma->At(1);
920  pm->SetNextPoint(array2[3*ist1+3], array2[3*ist1+4], array2[3*ist1+5]);
921  printf("### BOUNDARY MISSED BACK #########################\n");
922  ist1++;
923  break;
924  } else {
925  ist1++;
926  break;
927  }
928  }
929  }
930  }
931  }
932  pm = (TPolyMarker3D*)pma->At(0);
933  pm->Draw("SAME");
934  pm = (TPolyMarker3D*)pma->At(1);
935  pm->Draw("SAME");
936  pm = (TPolyMarker3D*)pma->At(2);
937  pm->Draw("SAME");
938  if (gPad) {
939  gPad->Modified();
940  gPad->Update();
941  }
942  delete [] array1;
943  delete [] array2;
944 }
945 
946 ////////////////////////////////////////////////////////////////////////////////
947 /// Clean-up the mesh of pcon/pgon from useless points
948 
949 void TGeoChecker::CleanPoints(Double_t *points, Int_t &numPoints) const
950 {
951  Int_t ipoint = 0;
952  Int_t j, k=0;
953  Double_t rsq;
954  for (Int_t i=0; i<numPoints; i++) {
955  j = 3*i;
956  rsq = points[j]*points[j]+points[j+1]*points[j+1];
957  if (rsq < 1.e-10) continue;
958  points[k] = points[j];
959  points[k+1] = points[j+1];
960  points[k+2] = points[j+2];
961  ipoint++;
962  k = 3*ipoint;
963  }
964  numPoints = ipoint;
965 }
966 
967 ////////////////////////////////////////////////////////////////////////////////
968 /// Check if the 2 non-assembly volume candidates overlap/extrude. Returns overlap object.
969 
970 TGeoOverlap *TGeoChecker::MakeCheckOverlap(const char *name, TGeoVolume *vol1, TGeoVolume *vol2, TGeoMatrix *mat1, TGeoMatrix *mat2, Bool_t isovlp, Double_t ovlp)
971 {
972  TGeoOverlap *nodeovlp = 0;
973  Int_t numPoints1 = fBuff1->NbPnts();
974  Int_t numSegs1 = fBuff1->NbSegs();
975  Int_t numPols1 = fBuff1->NbPols();
976  Int_t numPoints2 = fBuff2->NbPnts();
977  Int_t numSegs2 = fBuff2->NbSegs();
978  Int_t numPols2 = fBuff2->NbPols();
979  Int_t ip;
980  Bool_t extrude, isextrusion, isoverlapping;
981  Double_t *points1 = fBuff1->fPnts;
982  Double_t *points2 = fBuff2->fPnts;
983  Double_t local[3], local1[3];
984  Double_t point[3];
985  Double_t safety = TGeoShape::Big();
986  Double_t tolerance = TGeoShape::Tolerance();
987  if (vol1->IsAssembly() || vol2->IsAssembly()) return nodeovlp;
988  TGeoShape *shape1 = vol1->GetShape();
989  TGeoShape *shape2 = vol2->GetShape();
990  OpProgress("refresh", 0,0,NULL,kFALSE,kTRUE);
991  shape1->GetMeshNumbers(numPoints1, numSegs1, numPols1);
992  if (fBuff1->fID != (TObject*)shape1) {
993  // Fill first buffer.
994  fBuff1->SetRawSizes(TMath::Max(numPoints1,fNmeshPoints), 3*TMath::Max(numPoints1,fNmeshPoints), 0, 0, 0, 0);
995  points1 = fBuff1->fPnts;
996  if (shape1->GetPointsOnSegments(fNmeshPoints, points1)) {
997  numPoints1 = fNmeshPoints;
998  } else {
999  shape1->SetPoints(points1);
1000  }
1001 // if (shape1->InheritsFrom(TGeoPcon::Class())) {
1002 // CleanPoints(points1, numPoints1);
1003 // fBuff1->SetRawSizes(numPoints1, 3*numPoints1,0, 0, 0, 0);
1004 // }
1005  fBuff1->fID = shape1;
1006  }
1007  shape2->GetMeshNumbers(numPoints2, numSegs2, numPols2);
1008  if (fBuff2->fID != (TObject*)shape2) {
1009  // Fill second buffer.
1010  fBuff2->SetRawSizes(TMath::Max(numPoints2,fNmeshPoints), 3*TMath::Max(numPoints2,fNmeshPoints), 0, 0, 0,0);
1011  points2 = fBuff2->fPnts;
1012  if (shape2->GetPointsOnSegments(fNmeshPoints, points2)) {
1013  numPoints2 = fNmeshPoints;
1014  } else {
1015  shape2->SetPoints(points2);
1016  }
1017 // if (shape2->InheritsFrom(TGeoPcon::Class())) {
1018 // CleanPoints(points2, numPoints2);
1019 // fBuff1->SetRawSizes(numPoints2, 3*numPoints2,0,0,0,0);
1020 // }
1021  fBuff2->fID = shape2;
1022  }
1023 
1024  if (!isovlp) {
1025  // Extrusion case. Test vol2 extrude vol1.
1026  isextrusion=kFALSE;
1027  // loop all points of the daughter
1028  for (ip=0; ip<numPoints2; ip++) {
1029  memcpy(local, &points2[3*ip], 3*sizeof(Double_t));
1030  if (TMath::Abs(local[0])<tolerance && TMath::Abs(local[1])<tolerance) continue;
1031  mat2->LocalToMaster(local, point);
1032  mat1->MasterToLocal(point, local);
1033  extrude = !shape1->Contains(local);
1034  if (extrude) {
1035  safety = shape1->Safety(local, kFALSE);
1036  if (safety<ovlp) extrude=kFALSE;
1037  }
1038  if (extrude) {
1039  if (!isextrusion) {
1040  isextrusion = kTRUE;
1041  nodeovlp = new TGeoOverlap(name, vol1, vol2, mat1,mat2,kFALSE,safety);
1042  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1043  fGeoManager->AddOverlap(nodeovlp);
1044  } else {
1045  if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1046  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1047  }
1048  }
1049  }
1050  // loop all points of the mother
1051  for (ip=0; ip<numPoints1; ip++) {
1052  memcpy(local, &points1[3*ip], 3*sizeof(Double_t));
1053  if (local[0]<1e-10 && local[1]<1e-10) continue;
1054  mat1->LocalToMaster(local, point);
1055  mat2->MasterToLocal(point, local1);
1056  extrude = shape2->Contains(local1);
1057  if (extrude) {
1058  // skip points on mother mesh that have no neighborhood outside mother
1059  safety = shape1->Safety(local,kTRUE);
1060  if (safety>1E-6) {
1061  extrude = kFALSE;
1062  } else {
1063  safety = shape2->Safety(local1,kTRUE);
1064  if (safety<ovlp) extrude=kFALSE;
1065  }
1066  }
1067  if (extrude) {
1068  if (!isextrusion) {
1069  isextrusion = kTRUE;
1070  nodeovlp = new TGeoOverlap(name, vol1,vol2,mat1,mat2,kFALSE,safety);
1071  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1072  fGeoManager->AddOverlap(nodeovlp);
1073  } else {
1074  if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1075  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1076  }
1077  }
1078  }
1079  return nodeovlp;
1080  }
1081  // Check overlap
1082  Bool_t overlap;
1083  overlap = kFALSE;
1084  isoverlapping = kFALSE;
1085  // loop all points of first candidate
1086  for (ip=0; ip<numPoints1; ip++) {
1087  memcpy(local, &points1[3*ip], 3*sizeof(Double_t));
1088  if (local[0]<1e-10 && local[1]<1e-10) continue;
1089  mat1->LocalToMaster(local, point);
1090  mat2->MasterToLocal(point, local); // now point in local reference of second
1091  overlap = shape2->Contains(local);
1092  if (overlap) {
1093  safety = shape2->Safety(local, kTRUE);
1094  if (safety<ovlp) overlap=kFALSE;
1095  }
1096  if (overlap) {
1097  if (!isoverlapping) {
1098  isoverlapping = kTRUE;
1099  nodeovlp = new TGeoOverlap(name,vol1,vol2,mat1,mat2,kTRUE,safety);
1100  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1101  fGeoManager->AddOverlap(nodeovlp);
1102  } else {
1103  if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1104  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1105  }
1106  }
1107  }
1108  // loop all points of second candidate
1109  for (ip=0; ip<numPoints2; ip++) {
1110  memcpy(local, &points2[3*ip], 3*sizeof(Double_t));
1111  if (local[0]<1e-10 && local[1]<1e-10) continue;
1112  mat2->LocalToMaster(local, point);
1113  mat1->MasterToLocal(point, local); // now point in local reference of first
1114  overlap = shape1->Contains(local);
1115  if (overlap) {
1116  safety = shape1->Safety(local, kTRUE);
1117  if (safety<ovlp) overlap=kFALSE;
1118  }
1119  if (overlap) {
1120  if (!isoverlapping) {
1121  isoverlapping = kTRUE;
1122  nodeovlp = new TGeoOverlap(name,vol1,vol2,mat1,mat2,kTRUE,safety);
1123  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1124  fGeoManager->AddOverlap(nodeovlp);
1125  } else {
1126  if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
1127  nodeovlp->SetNextPoint(point[0],point[1],point[2]);
1128  }
1129  }
1130  }
1131  return nodeovlp;
1132 }
1133 
1134 ////////////////////////////////////////////////////////////////////////////////
1135 /// Check illegal overlaps for volume VOL within a limit OVLP by sampling npoints
1136 /// inside the volume shape.
1137 
1138 void TGeoChecker::CheckOverlapsBySampling(TGeoVolume *vol, Double_t /* ovlp */, Int_t npoints) const
1139 {
1140  Int_t nd = vol->GetNdaughters();
1141  if (nd<2) return;
1142  TGeoVoxelFinder *voxels = vol->GetVoxels();
1143  if (!voxels) return;
1144  if (voxels->NeedRebuild()) {
1145  voxels->Voxelize();
1146  vol->FindOverlaps();
1147  }
1148  TGeoBBox *box = (TGeoBBox*)vol->GetShape();
1149  TGeoShape *shape;
1150  TGeoNode *node;
1151  Double_t dx = box->GetDX();
1152  Double_t dy = box->GetDY();
1153  Double_t dz = box->GetDZ();
1154  Double_t pt[3];
1155  Double_t local[3];
1156  Int_t *check_list = 0;
1157  Int_t ncheck = 0;
1158  const Double_t *orig = box->GetOrigin();
1159  Int_t ipoint = 0;
1160  Int_t itry = 0;
1161  Int_t iovlp = 0;
1162  Int_t id=0, id0=0, id1=0;
1163  Bool_t in, incrt;
1164  Double_t safe;
1165  TString name1 = "";
1166  TString name2 = "";
1167  TGeoOverlap **flags = 0;
1168  TGeoNode *node1, *node2;
1169  Int_t novlps = 0;
1170  TGeoHMatrix mat1, mat2;
1171 // Int_t tid = TGeoManager::ThreadId();
1172  TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
1173  TGeoStateInfo &td = *nav->GetCache()->GetInfo();
1174  while (ipoint < npoints) {
1175  // Shoot randomly in the bounding box.
1176  pt[0] = orig[0] - dx + 2.*dx*gRandom->Rndm();
1177  pt[1] = orig[1] - dy + 2.*dy*gRandom->Rndm();
1178  pt[2] = orig[2] - dz + 2.*dz*gRandom->Rndm();
1179  if (!vol->Contains(pt)) {
1180  itry++;
1181  if (itry>10000 && !ipoint) {
1182  Error("CheckOverlapsBySampling", "No point inside volume!!! - aborting");
1183  break;
1184  }
1185  continue;
1186  }
1187  // Check if the point is inside one or more daughters
1188  in = kFALSE;
1189  ipoint++;
1190  check_list = voxels->GetCheckList(pt, ncheck, td);
1191  if (!check_list || ncheck<2) continue;
1192  for (id=0; id<ncheck; id++) {
1193  id0 = check_list[id];
1194  node = vol->GetNode(id0);
1195  // Ignore MANY's
1196  if (node->IsOverlapping()) continue;
1197  node->GetMatrix()->MasterToLocal(pt,local);
1198  shape = node->GetVolume()->GetShape();
1199  incrt = shape->Contains(local);
1200  if (!incrt) continue;
1201  if (!in) {
1202  in = kTRUE;
1203  id1 = id0;
1204  continue;
1205  }
1206  // The point is inside 2 or more daughters, check safety
1207  safe = shape->Safety(local, kTRUE);
1208 // if (safe < ovlp) continue;
1209  // We really have found an overlap -> store the point in a container
1210  iovlp++;
1211  if (!novlps) {
1212  flags = new TGeoOverlap*[nd*nd];
1213  memset(flags, 0, nd*nd*sizeof(TGeoOverlap*));
1214  }
1215  TGeoOverlap *nodeovlp = flags[nd*id1+id0];
1216  if (!nodeovlp) {
1217  novlps++;
1218  node1 = vol->GetNode(id1);
1219  name1 = node1->GetName();
1220  mat1 = node1->GetMatrix();
1221  Int_t cindex = node1->GetVolume()->GetCurrentNodeIndex();
1222  while (cindex >= 0) {
1223  node1 = node1->GetVolume()->GetNode(cindex);
1224  name1 += TString::Format("/%s", node1->GetName());
1225  mat1.Multiply(node1->GetMatrix());
1226  cindex = node1->GetVolume()->GetCurrentNodeIndex();
1227  }
1228  node2 = vol->GetNode(id0);
1229  name2 = node2->GetName();
1230  mat2 = node2->GetMatrix();
1231  cindex = node2->GetVolume()->GetCurrentNodeIndex();
1232  while (cindex >= 0) {
1233  node2 = node2->GetVolume()->GetNode(cindex);
1234  name2 += TString::Format("/%s", node2->GetName());
1235  mat2.Multiply(node2->GetMatrix());
1236  cindex = node2->GetVolume()->GetCurrentNodeIndex();
1237  }
1238  nodeovlp = new TGeoOverlap(TString::Format("Volume %s: node %s overlapping %s",
1239  vol->GetName(), name1.Data(), name2.Data()), node1->GetVolume(),node2->GetVolume(),
1240  &mat1,&mat2, kTRUE, safe);
1241  flags[nd*id1+id0] = nodeovlp;
1242  fGeoManager->AddOverlap(nodeovlp);
1243  }
1244  // Max 100 points per marker
1245  if (nodeovlp->GetPolyMarker()->GetN()<100) nodeovlp->SetNextPoint(pt[0],pt[1],pt[2]);
1246  if (nodeovlp->GetOverlap()<safe) nodeovlp->SetOverlap(safe);
1247  }
1248  }
1249  nav->GetCache()->ReleaseInfo();
1250  if (flags) delete [] flags;
1251  if (!novlps) return;
1252  Double_t capacity = vol->GetShape()->Capacity();
1253  capacity *= Double_t(iovlp)/Double_t(npoints);
1254  Double_t err = 1./TMath::Sqrt(Double_t(iovlp));
1255  Info("CheckOverlapsBySampling", "#Found %d overlaps adding-up to %g +/- %g [cm3] for daughters of %s",
1256  novlps, capacity, err*capacity, vol->GetName());
1257 }
1258 
1259 ////////////////////////////////////////////////////////////////////////////////
1260 /// Compute number of overlaps combinations to check per volume
1261 
1262 Int_t TGeoChecker::NChecksPerVolume(TGeoVolume *vol)
1263 {
1264  if (vol->GetFinder()) return 0;
1265  UInt_t nd = vol->GetNdaughters();
1266  if (!nd) return 0;
1267  Bool_t is_assembly = vol->IsAssembly();
1268  TGeoIterator next1(vol);
1269  TGeoIterator next2(vol);
1270  Int_t nchecks = 0;
1271  TGeoNode *node;
1272  UInt_t id;
1273  if (!is_assembly) {
1274  while ((node=next1())) {
1275  if (node->IsOverlapping()) {
1276  next1.Skip();
1277  continue;
1278  }
1279  if (!node->GetVolume()->IsAssembly()) {
1280  nchecks++;
1281  next1.Skip();
1282  }
1283  }
1284  }
1285  // now check if the daughters overlap with each other
1286  if (nd<2) return nchecks;
1287  TGeoVoxelFinder *vox = vol->GetVoxels();
1288  if (!vox) return nchecks;
1289 
1290 
1291  TGeoNode *node1, *node01, *node02;
1292  Int_t novlp;
1293  Int_t *ovlps;
1294  Int_t ko;
1295  UInt_t io;
1296  for (id=0; id<nd; id++) {
1297  node01 = vol->GetNode(id);
1298  if (node01->IsOverlapping()) continue;
1299  vox->FindOverlaps(id);
1300  ovlps = node01->GetOverlaps(novlp);
1301  if (!ovlps) continue;
1302  for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
1303  io = ovlps[ko]; // (node1, shaped, matrix1, points, fBuff1)
1304  if (io<=id) continue;
1305  node02 = vol->GetNode(io);
1306  if (node02->IsOverlapping()) continue;
1307 
1308  // We have to check node against node1, but they may be assemblies
1309  if (node01->GetVolume()->IsAssembly()) {
1310  next1.Reset(node01->GetVolume());
1311  while ((node=next1())) {
1312  if (!node->GetVolume()->IsAssembly()) {
1313  if (node02->GetVolume()->IsAssembly()) {
1314  next2.Reset(node02->GetVolume());
1315  while ((node1=next2())) {
1316  if (!node1->GetVolume()->IsAssembly()) {
1317  nchecks++;
1318  next2.Skip();
1319  }
1320  }
1321  } else {
1322  nchecks++;
1323  }
1324  next1.Skip();
1325  }
1326  }
1327  } else {
1328  // node not assembly
1329  if (node02->GetVolume()->IsAssembly()) {
1330  next2.Reset(node02->GetVolume());
1331  while ((node1=next2())) {
1332  if (!node1->GetVolume()->IsAssembly()) {
1333  nchecks++;
1334  next2.Skip();
1335  }
1336  }
1337  } else {
1338  // node1 also not assembly
1339  nchecks++;
1340  }
1341  }
1342  }
1343  node01->SetOverlaps(0,0);
1344  }
1345  return nchecks;
1346 }
1347 
1348 ////////////////////////////////////////////////////////////////////////////////
1349 /// Check illegal overlaps for volume VOL within a limit OVLP.
1350 
1351 void TGeoChecker::CheckOverlaps(const TGeoVolume *vol, Double_t ovlp, Option_t *option)
1352 {
1353  if (vol->GetFinder()) return;
1354  UInt_t nd = vol->GetNdaughters();
1355  if (!nd) return;
1356  TGeoShape::SetTransform(gGeoIdentity);
1357  fNchecks = NChecksPerVolume((TGeoVolume*)vol);
1358  Bool_t sampling = kFALSE;
1359  TString opt(option);
1360  opt.ToLower();
1361  if (opt.Contains("s")) sampling = kTRUE;
1362  if (opt.Contains("f")) fFullCheck = kTRUE;
1363  else fFullCheck = kFALSE;
1364  if (sampling) {
1365  opt = opt.Strip(TString::kLeading, 's');
1366  Int_t npoints = atoi(opt.Data());
1367  if (!npoints) npoints = 1000000;
1368  CheckOverlapsBySampling((TGeoVolume*)vol, ovlp, npoints);
1369  return;
1370  }
1371  Bool_t is_assembly = vol->IsAssembly();
1372  TGeoIterator next1((TGeoVolume*)vol);
1373  TGeoIterator next2((TGeoVolume*)vol);
1374  TString path;
1375  // first, test if daughters extrude their container
1376  TGeoNode * node, *nodecheck;
1377  TGeoChecker *checker = (TGeoChecker*)this;
1378 
1379 // TGeoOverlap *nodeovlp = 0;
1380  UInt_t id;
1381  Int_t level;
1382 // Check extrusion only for daughters of a non-assembly volume
1383  if (!is_assembly) {
1384  while ((node=next1())) {
1385  if (node->IsOverlapping()) {
1386  next1.Skip();
1387  continue;
1388  }
1389  if (!node->GetVolume()->IsAssembly()) {
1390  if (fSelectedNode) {
1391  // We have to check only overlaps of the selected node (or real daughters if an assembly)
1392  if ((fSelectedNode != node) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1393  next1.Skip();
1394  continue;
1395  }
1396  if (node != fSelectedNode) {
1397  level = next1.GetLevel();
1398  while ((nodecheck=next1.GetNode(level--))) {
1399  if (nodecheck == fSelectedNode) break;
1400  }
1401  if (!nodecheck) {
1402  next1.Skip();
1403  continue;
1404  }
1405  }
1406  }
1407  next1.GetPath(path);
1408  checker->MakeCheckOverlap(TString::Format("%s extruded by: %s", vol->GetName(),path.Data()),
1409  (TGeoVolume*)vol,node->GetVolume(),gGeoIdentity,(TGeoMatrix*)next1.GetCurrentMatrix(),kFALSE,ovlp);
1410  next1.Skip();
1411  }
1412  }
1413  }
1414 
1415  // now check if the daughters overlap with each other
1416  if (nd<2) return;
1417  TGeoVoxelFinder *vox = vol->GetVoxels();
1418  if (!vox) {
1419  Warning("CheckOverlaps", "Volume %s with %i daughters but not voxelized", vol->GetName(),nd);
1420  return;
1421  }
1422  if (vox->NeedRebuild()) {
1423  vox->Voxelize();
1424  vol->FindOverlaps();
1425  }
1426  TGeoNode *node1, *node01, *node02;
1427  TGeoHMatrix hmat1, hmat2;
1428  TString path1;
1429  Int_t novlp;
1430  Int_t *ovlps;
1431  Int_t ko;
1432  UInt_t io;
1433  for (id=0; id<nd; id++) {
1434  node01 = vol->GetNode(id);
1435  if (node01->IsOverlapping()) continue;
1436  vox->FindOverlaps(id);
1437  ovlps = node01->GetOverlaps(novlp);
1438  if (!ovlps) continue;
1439  next1.SetTopName(node01->GetName());
1440  path = node01->GetName();
1441  for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
1442  io = ovlps[ko]; // (node1, shaped, matrix1, points, fBuff1)
1443  if (io<=id) continue;
1444  node02 = vol->GetNode(io);
1445  if (node02->IsOverlapping()) continue;
1446  // Try to fasten-up things...
1447 // if (!TGeoBBox::AreOverlapping((TGeoBBox*)node01->GetVolume()->GetShape(), node01->GetMatrix(),
1448 // (TGeoBBox*)node02->GetVolume()->GetShape(), node02->GetMatrix())) continue;
1449  next2.SetTopName(node02->GetName());
1450  path1 = node02->GetName();
1451 
1452  // We have to check node against node1, but they may be assemblies
1453  if (node01->GetVolume()->IsAssembly()) {
1454  next1.Reset(node01->GetVolume());
1455  while ((node=next1())) {
1456  if (!node->GetVolume()->IsAssembly()) {
1457  next1.GetPath(path);
1458  hmat1 = node01->GetMatrix();
1459  hmat1 *= *next1.GetCurrentMatrix();
1460  if (node02->GetVolume()->IsAssembly()) {
1461  next2.Reset(node02->GetVolume());
1462  while ((node1=next2())) {
1463  if (!node1->GetVolume()->IsAssembly()) {
1464  if (fSelectedNode) {
1465  // We have to check only overlaps of the selected node (or real daughters if an assembly)
1466  if ((fSelectedNode != node) && (fSelectedNode != node1) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1467  next2.Skip();
1468  continue;
1469  }
1470  if ((fSelectedNode != node1) && (fSelectedNode != node)) {
1471  level = next2.GetLevel();
1472  while ((nodecheck=next2.GetNode(level--))) {
1473  if (nodecheck == fSelectedNode) break;
1474  }
1475  if (node02 == fSelectedNode) nodecheck = node02;
1476  if (!nodecheck) {
1477  level = next1.GetLevel();
1478  while ((nodecheck=next1.GetNode(level--))) {
1479  if (nodecheck == fSelectedNode) break;
1480  }
1481  }
1482  if (node01 == fSelectedNode) nodecheck = node01;
1483  if (!nodecheck) {
1484  next2.Skip();
1485  continue;
1486  }
1487  }
1488  }
1489  next2.GetPath(path1);
1490  hmat2 = node02->GetMatrix();
1491  hmat2 *= *next2.GetCurrentMatrix();
1492  checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1493  node->GetVolume(),node1->GetVolume(),&hmat1,&hmat2,kTRUE,ovlp);
1494  next2.Skip();
1495  }
1496  }
1497  } else {
1498  if (fSelectedNode) {
1499  // We have to check only overlaps of the selected node (or real daughters if an assembly)
1500  if ((fSelectedNode != node) && (fSelectedNode != node02) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1501  next1.Skip();
1502  continue;
1503  }
1504  if ((fSelectedNode != node) && (fSelectedNode != node02)) {
1505  level = next1.GetLevel();
1506  while ((nodecheck=next1.GetNode(level--))) {
1507  if (nodecheck == fSelectedNode) break;
1508  }
1509  if (node01 == fSelectedNode) nodecheck = node01;
1510  if (!nodecheck) {
1511  next1.Skip();
1512  continue;
1513  }
1514  }
1515  }
1516  checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1517  node->GetVolume(),node02->GetVolume(),&hmat1,node02->GetMatrix(),kTRUE,ovlp);
1518  }
1519  next1.Skip();
1520  }
1521  }
1522  } else {
1523  // node not assembly
1524  if (node02->GetVolume()->IsAssembly()) {
1525  next2.Reset(node02->GetVolume());
1526  while ((node1=next2())) {
1527  if (!node1->GetVolume()->IsAssembly()) {
1528  if (fSelectedNode) {
1529  // We have to check only overlaps of the selected node (or real daughters if an assembly)
1530  if ((fSelectedNode != node1) && (fSelectedNode != node01) && (!fSelectedNode->GetVolume()->IsAssembly())) {
1531  next2.Skip();
1532  continue;
1533  }
1534  if ((fSelectedNode != node1) && (fSelectedNode != node01)) {
1535  level = next2.GetLevel();
1536  while ((nodecheck=next2.GetNode(level--))) {
1537  if (nodecheck == fSelectedNode) break;
1538  }
1539  if (node02 == fSelectedNode) nodecheck = node02;
1540  if (!nodecheck) {
1541  next2.Skip();
1542  continue;
1543  }
1544  }
1545  }
1546  next2.GetPath(path1);
1547  hmat2 = node02->GetMatrix();
1548  hmat2 *= *next2.GetCurrentMatrix();
1549  checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1550  node01->GetVolume(),node1->GetVolume(),node01->GetMatrix(),&hmat2,kTRUE,ovlp);
1551  next2.Skip();
1552  }
1553  }
1554  } else {
1555  // node1 also not assembly
1556  if (fSelectedNode && (fSelectedNode != node01) && (fSelectedNode != node02)) continue;
1557  checker->MakeCheckOverlap(TString::Format("%s/%s overlapping %s/%s", vol->GetName(),path.Data(),vol->GetName(),path1.Data()),
1558  node01->GetVolume(),node02->GetVolume(),node01->GetMatrix(),node02->GetMatrix(),kTRUE,ovlp);
1559  }
1560  }
1561  }
1562  node01->SetOverlaps(0,0);
1563  }
1564 }
1565 
1566 ////////////////////////////////////////////////////////////////////////////////
1567 /// Print the current list of overlaps held by the manager class.
1568 
1569 void TGeoChecker::PrintOverlaps() const
1570 {
1571  TIter next(fGeoManager->GetListOfOverlaps());
1572  TGeoOverlap *ov;
1573  printf("=== Overlaps for %s ===\n", fGeoManager->GetName());
1574  while ((ov=(TGeoOverlap*)next())) ov->PrintInfo();
1575 }
1576 
1577 ////////////////////////////////////////////////////////////////////////////////
1578 /// Draw point (x,y,z) over the picture of the daughters of the volume containing this point.
1579 /// Generates a report regarding the path to the node containing this point and the distance to
1580 /// the closest boundary.
1581 
1582 void TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z, Option_t *)
1583 {
1584  Double_t point[3];
1585  Double_t local[3];
1586  point[0] = x;
1587  point[1] = y;
1588  point[2] = z;
1589  TGeoVolume *vol = fGeoManager->GetTopVolume();
1590  if (fVsafe) {
1591  TGeoNode *old = fVsafe->GetNode("SAFETY_1");
1592  if (old) fVsafe->GetNodes()->RemoveAt(vol->GetNdaughters()-1);
1593  }
1594 // if (vol != fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
1595  TGeoNode *node = fGeoManager->FindNode(point[0], point[1], point[2]);
1596  fGeoManager->MasterToLocal(point, local);
1597  // get current node
1598  printf("=== Check current point : (%g, %g, %g) ===\n", point[0], point[1], point[2]);
1599  printf(" - path : %s\n", fGeoManager->GetPath());
1600  // get corresponding volume
1601  if (node) vol = node->GetVolume();
1602  // compute safety distance (distance to boundary ignored)
1603  Double_t close = fGeoManager->Safety();
1604  printf("Safety radius : %f\n", close);
1605  if (close>1E-4) {
1606  TGeoVolume *sph = fGeoManager->MakeSphere("SAFETY", vol->GetMedium(), 0, close, 0,180,0,360);
1607  sph->SetLineColor(2);
1608  sph->SetLineStyle(3);
1609  vol->AddNode(sph,1,new TGeoTranslation(local[0], local[1], local[2]));
1610  fVsafe = vol;
1611  }
1612  TPolyMarker3D *pm = new TPolyMarker3D();
1613  pm->SetMarkerColor(2);
1614  pm->SetMarkerStyle(8);
1615  pm->SetMarkerSize(0.5);
1616  pm->SetNextPoint(local[0], local[1], local[2]);
1617  if (vol->GetNdaughters()<2) fGeoManager->SetTopVisible();
1618  else fGeoManager->SetTopVisible(kFALSE);
1619  fGeoManager->SetVisLevel(1);
1620  if (!vol->IsVisible()) vol->SetVisibility(kTRUE);
1621  vol->Draw();
1622  pm->Draw("SAME");
1623  gPad->Modified();
1624  gPad->Update();
1625 }
1626 
1627 ////////////////////////////////////////////////////////////////////////////////
1628 /// Test for shape navigation methods. Summary for test numbers:
1629 /// - 1: DistFromInside/Outside. Sample points inside the shape. Generate
1630 /// directions randomly in cos(theta). Compute DistFromInside and move the
1631 /// point with bigger distance. Compute DistFromOutside back from new point.
1632 /// Plot d-(d1+d2)
1633 /// - 2: Safety test. Sample points inside the bounding and compute safety. Generate
1634 /// directions randomly in cos(theta) and compute distance to boundary. Check if
1635 /// distance to boundary is bigger than safety
1636 
1637 void TGeoChecker::CheckShape(TGeoShape *shape, Int_t testNo, Int_t nsamples, Option_t *option)
1638 {
1639  switch (testNo) {
1640  case 1:
1641  ShapeDistances(shape, nsamples, option);
1642  break;
1643  case 2:
1644  ShapeSafety(shape, nsamples, option);
1645  break;
1646  case 3:
1647  ShapeNormal(shape, nsamples, option);
1648  break;
1649  default:
1650  Error("CheckShape", "Test number %d not existent", testNo);
1651  }
1652 }
1653 
1654 ////////////////////////////////////////////////////////////////////////////////
1655 /// Test TGeoShape::DistFromInside/Outside. Sample points inside the shape. Generate
1656 /// directions randomly in cos(theta). Compute d1 = DistFromInside and move the
1657 /// point on the boundary. Compute DistFromOutside and propagate with d2 making sure that
1658 /// the shape is not re-entered. Swap direction and call DistFromOutside that
1659 /// should fall back on the same point on the boundary (at d2). Propagate back on boundary
1660 /// then compute DistFromInside that should be bigger than d1.
1661 /// Plot d-(d1+d2)
1662 
1663 void TGeoChecker::ShapeDistances(TGeoShape *shape, Int_t nsamples, Option_t *)
1664 {
1665  Double_t dx = ((TGeoBBox*)shape)->GetDX();
1666  Double_t dy = ((TGeoBBox*)shape)->GetDY();
1667  Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1668  Double_t dmax = 2.*TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1669  Double_t d1, d2, dmove, dnext;
1670  Int_t itot = 0;
1671  // Number of tracks shot for every point inside the shape
1672  const Int_t kNtracks = 1000;
1673  Int_t n10 = nsamples/10;
1674  Int_t i,j;
1675  Double_t point[3], pnew[3];
1676  Double_t dir[3], dnew[3];
1677  Double_t theta, phi, delta;
1678  TPolyMarker3D *pmfrominside = 0;
1679  TPolyMarker3D *pmfromoutside = 0;
1680  TH1D *hist = new TH1D("hTest1", "Residual distance from inside/outside",200,-20, 0);
1681  hist->GetXaxis()->SetTitle("delta[cm] - first bin=overflow");
1682  hist->GetYaxis()->SetTitle("count");
1683  hist->SetMarkerStyle(kFullCircle);
1684 
1685  if (!fTimer) fTimer = new TStopwatch();
1686  fTimer->Reset();
1687  fTimer->Start();
1688  while (itot<nsamples) {
1689  Bool_t inside = kFALSE;
1690  while (!inside) {
1691  point[0] = gRandom->Uniform(-dx,dx);
1692  point[1] = gRandom->Uniform(-dy,dy);
1693  point[2] = gRandom->Uniform(-dz,dz);
1694  inside = shape->Contains(point);
1695  }
1696  itot++;
1697  if (n10) {
1698  if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1699  }
1700  for (i=0; i<kNtracks; i++) {
1701  phi = 2*TMath::Pi()*gRandom->Rndm();
1702  theta= TMath::ACos(1.-2.*gRandom->Rndm());
1703  dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
1704  dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
1705  dir[2]=TMath::Cos(theta);
1706  dmove = dmax;
1707  // We have track direction, compute distance from inside
1708  d1 = shape->DistFromInside(point,dir,3);
1709  if (d1>dmove || d1<TGeoShape::Tolerance()) {
1710  // Bad distance or bbox size, to debug
1711  new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1712  shape->Draw();
1713  printf("DistFromInside: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) %f/%f(max)\n",
1714  point[0],point[1],point[2],dir[0],dir[1],dir[2], d1,dmove);
1715  pmfrominside = new TPolyMarker3D(2);
1716  pmfrominside->SetMarkerColor(kRed);
1717  pmfrominside->SetMarkerStyle(24);
1718  pmfrominside->SetMarkerSize(0.4);
1719  pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1720  for (j=0; j<3; j++) pnew[j] = point[j] + d1*dir[j];
1721  pmfrominside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1722  pmfrominside->Draw();
1723  return;
1724  }
1725  // Propagate BEFORE the boundary and make sure that DistFromOutside
1726  // does not return 0 (!!!)
1727  // Check if there is a second crossing
1728  for (j=0; j<3; j++) pnew[j] = point[j] + (d1-TGeoShape::Tolerance())*dir[j];
1729  dnext = shape->DistFromOutside(pnew,dir,3);
1730  if (d1+dnext<dmax) dmove = d1+0.5*dnext;
1731  // Move point and swap direction
1732  for (j=0; j<3; j++) {
1733  pnew[j] = point[j] + dmove*dir[j];
1734  dnew[j] = -dir[j];
1735  }
1736  // Compute now distance from outside
1737  d2 = shape->DistFromOutside(pnew,dnew,3);
1738  delta = dmove-d1-d2;
1739  if (TMath::Abs(delta)>1E-6 || dnext<2.*TGeoShape::Tolerance()) {
1740  // Error->debug this
1741  new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1742  shape->Draw();
1743  printf("Error: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d2=%f dmove=%f\n",
1744  point[0],point[1],point[2],dir[0],dir[1],dir[2], d1,d2,dmove);
1745  if (dnext<2.*TGeoShape::Tolerance()) {
1746  printf(" (*)DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) dnext = %f\n",
1747  point[0]+(d1-TGeoShape::Tolerance())*dir[0],
1748  point[1]+(d1-TGeoShape::Tolerance())*dir[1],
1749  point[2]+(d1-TGeoShape::Tolerance())*dir[2], dir[0],dir[1],dir[2],dnext);
1750  } else {
1751  printf(" DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) dnext = %f\n",
1752  point[0]+d1*dir[0],point[1]+d1*dir[1], point[2]+d1*dir[2], dir[0],dir[1],dir[2],dnext);
1753  }
1754  printf(" DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) = %f\n",
1755  pnew[0],pnew[1],pnew[2],dnew[0],dnew[1],dnew[2], d2);
1756  pmfrominside = new TPolyMarker3D(2);
1757  pmfrominside->SetMarkerStyle(24);
1758  pmfrominside->SetMarkerSize(0.4);
1759  pmfrominside->SetMarkerColor(kRed);
1760  pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1761  for (j=0; j<3; j++) point[j] += d1*dir[j];
1762  pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1763  pmfrominside->Draw();
1764  pmfromoutside = new TPolyMarker3D(2);
1765  pmfromoutside->SetMarkerStyle(20);
1766  pmfromoutside->SetMarkerStyle(7);
1767  pmfromoutside->SetMarkerSize(0.3);
1768  pmfromoutside->SetMarkerColor(kBlue);
1769  pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1770  for (j=0; j<3; j++) pnew[j] += d2*dnew[j];
1771  if (d2<1E10) pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1772  pmfromoutside->Draw();
1773  return;
1774  }
1775  // Compute distance from inside which should be bigger than d1
1776  for (j=0; j<3; j++) pnew[j] += (d2-TGeoShape::Tolerance())*dnew[j];
1777  dnext = shape->DistFromInside(pnew,dnew,3);
1778  if (dnext<d1-TGeoShape::Tolerance() || dnext>dmax) {
1779  new TCanvas("shape01", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1780  shape->Draw();
1781  printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d1p=%f\n",
1782  pnew[0],pnew[1],pnew[2],dnew[0],dnew[1],dnew[2],d1,dnext);
1783  pmfrominside = new TPolyMarker3D(2);
1784  pmfrominside->SetMarkerStyle(24);
1785  pmfrominside->SetMarkerSize(0.4);
1786  pmfrominside->SetMarkerColor(kRed);
1787  pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1788  for (j=0; j<3; j++) point[j] += d1*dir[j];
1789  pmfrominside->SetNextPoint(point[0],point[1],point[2]);
1790  pmfrominside->Draw();
1791  pmfromoutside = new TPolyMarker3D(2);
1792  pmfromoutside->SetMarkerStyle(20);
1793  pmfromoutside->SetMarkerStyle(7);
1794  pmfromoutside->SetMarkerSize(0.3);
1795  pmfromoutside->SetMarkerColor(kBlue);
1796  pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1797  for (j=0; j<3; j++) pnew[j] += dnext*dnew[j];
1798  if (d2<1E10) pmfromoutside->SetNextPoint(pnew[0],pnew[1],pnew[2]);
1799  pmfromoutside->Draw();
1800  return;
1801  }
1802  if (TMath::Abs(delta) < 1E-20) delta = 1E-30;
1803  hist->Fill(TMath::Max(TMath::Log(TMath::Abs(delta)),-20.));
1804  }
1805  }
1806  fTimer->Stop();
1807  fTimer->Print();
1808  new TCanvas("Test01", "Residuals DistFromInside/Outside", 800, 600);
1809  hist->Draw();
1810 }
1811 
1812 ////////////////////////////////////////////////////////////////////////////////
1813 /// Check of validity of safe distance for a given shape.
1814 /// Sample points inside the 2x bounding box and compute safety. Generate
1815 /// directions randomly in cos(theta) and compute distance to boundary. Check if
1816 /// distance to boundary is bigger than safety.
1817 
1818 void TGeoChecker::ShapeSafety(TGeoShape *shape, Int_t nsamples, Option_t *)
1819 {
1820  Double_t dx = ((TGeoBBox*)shape)->GetDX();
1821  Double_t dy = ((TGeoBBox*)shape)->GetDY();
1822  Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1823  // Number of tracks shot for every point inside the shape
1824  const Int_t kNtracks = 1000;
1825  Int_t n10 = nsamples/10;
1826  Int_t i;
1827  Double_t dist;
1828  Double_t point[3];
1829  Double_t dir[3];
1830  Double_t theta, phi;
1831  TPolyMarker3D *pm1 = 0;
1832  TPolyMarker3D *pm2 = 0;
1833  if (!fTimer) fTimer = new TStopwatch();
1834  fTimer->Reset();
1835  fTimer->Start();
1836  Int_t itot = 0;
1837  while (itot<nsamples) {
1838  Bool_t inside = kFALSE;
1839  point[0] = gRandom->Uniform(-2*dx,2*dx);
1840  point[1] = gRandom->Uniform(-2*dy,2*dy);
1841  point[2] = gRandom->Uniform(-2*dz,2*dz);
1842  inside = shape->Contains(point);
1843  Double_t safe = shape->Safety(point, inside);
1844  itot++;
1845  if (n10) {
1846  if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1847  }
1848  for (i=0; i<kNtracks; i++) {
1849  phi = 2*TMath::Pi()*gRandom->Rndm();
1850  theta= TMath::ACos(1.-2.*gRandom->Rndm());
1851  dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
1852  dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
1853  dir[2]=TMath::Cos(theta);
1854  if (inside) dist = shape->DistFromInside(point,dir,3);
1855  else dist = shape->DistFromOutside(point,dir,3);
1856  if (dist<safe) {
1857  printf("Error safety (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) safe=%f dist=%f\n",
1858  point[0],point[1],point[2], dir[0], dir[1], dir[2], safe, dist);
1859  new TCanvas("shape02", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1860  shape->Draw();
1861  pm1 = new TPolyMarker3D(2);
1862  pm1->SetMarkerStyle(24);
1863  pm1->SetMarkerSize(0.4);
1864  pm1->SetMarkerColor(kRed);
1865  pm1->SetNextPoint(point[0],point[1],point[2]);
1866  pm1->SetNextPoint(point[0]+safe*dir[0],point[1]+safe*dir[1],point[2]+safe*dir[2]);
1867  pm1->Draw();
1868  pm2 = new TPolyMarker3D(1);
1869  pm2->SetMarkerStyle(7);
1870  pm2->SetMarkerSize(0.3);
1871  pm2->SetMarkerColor(kBlue);
1872  pm2->SetNextPoint(point[0]+dist*dir[0],point[1]+dist*dir[1],point[2]+dist*dir[2]);
1873  pm2->Draw();
1874  return;
1875  }
1876  }
1877  }
1878  fTimer->Stop();
1879  fTimer->Print();
1880 }
1881 
1882 ////////////////////////////////////////////////////////////////////////////////
1883 /// Check of validity of the normal for a given shape.
1884 /// Sample points inside the shape. Generate directions randomly in cos(theta)
1885 /// and propagate to boundary. Compute normal and safety at crossing point, plot
1886 /// the point and generate a random direction so that (dir) dot (norm) <0.
1887 
1888 void TGeoChecker::ShapeNormal(TGeoShape *shape, Int_t nsamples, Option_t *)
1889 {
1890  Double_t dx = ((TGeoBBox*)shape)->GetDX();
1891  Double_t dy = ((TGeoBBox*)shape)->GetDY();
1892  Double_t dz = ((TGeoBBox*)shape)->GetDZ();
1893  Double_t dmax = 2.*TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1894  // Number of tracks shot for every point inside the shape
1895  const Int_t kNtracks = 1000;
1896  Int_t n10 = nsamples/10;
1897  Int_t itot = 0, errcnt = 0, errsame=0;
1898  Int_t i;
1899  Double_t dist, olddist, safe, dot;
1900  Double_t theta, phi, ndotd;
1901  Double_t *spoint = new Double_t[3*nsamples];
1902  Double_t *sdir = new Double_t[3*nsamples];
1903  while (itot<nsamples) {
1904  Bool_t inside = kFALSE;
1905  while (!inside) {
1906  spoint[3*itot] = gRandom->Uniform(-dx,dx);
1907  spoint[3*itot+1] = gRandom->Uniform(-dy,dy);
1908  spoint[3*itot+2] = gRandom->Uniform(-dz,dz);
1909  inside = shape->Contains(&spoint[3*itot]);
1910  }
1911  phi = 2*TMath::Pi()*gRandom->Rndm();
1912  theta= TMath::ACos(1.-2.*gRandom->Rndm());
1913  sdir[3*itot] = TMath::Sin(theta)*TMath::Cos(phi);
1914  sdir[3*itot+1] = TMath::Sin(theta)*TMath::Sin(phi);
1915  sdir[3*itot+2] = TMath::Cos(theta);
1916  itot++;
1917  }
1918  Double_t point[3],newpoint[3], oldpoint[3];
1919  Double_t dir[3], olddir[3];
1920  Double_t norm[3], newnorm[3], oldnorm[3];
1921  TCanvas *errcanvas = 0;
1922  TPolyMarker3D *pm1 = 0;
1923  TPolyMarker3D *pm2 = 0;
1924  pm2 = new TPolyMarker3D();
1925 // pm2->SetMarkerStyle(24);
1926  pm2->SetMarkerSize(0.2);
1927  pm2->SetMarkerColor(kBlue);
1928  if (!fTimer) fTimer = new TStopwatch();
1929  fTimer->Reset();
1930  fTimer->Start();
1931  for (itot = 0; itot<nsamples; itot++) {
1932  if (n10) {
1933  if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nsamples));
1934  }
1935  oldnorm[0] = oldnorm[1] = oldnorm[2] = 0.;
1936  olddist = 0.;
1937  for (Int_t j=0; j<3; j++) {
1938  oldpoint[j] = point[j] = spoint[3*itot+j];
1939  olddir[j] = dir[j] = sdir[3*itot+j];
1940  }
1941  for (i=0; i<kNtracks; i++) {
1942  if (errcnt>0) break;
1943  dist = shape->DistFromInside(point,dir,3);
1944  for (Int_t j=0; j<3; j++) {
1945  newpoint[j] = point[j] + dist*dir[j];
1946  }
1947  shape->ComputeNormal(newpoint,dir,newnorm);
1948 
1949  dot = olddir[0]*oldnorm[0]+olddir[1]*oldnorm[1]+ olddir[2]*oldnorm[2];
1950  if (!shape->Contains(point) && shape->Safety(point,kFALSE)>1.E-3) {
1951  errcnt++;
1952  printf("Error point outside (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
1953  point[0],point[1],point[2], dir[0], dir[1], dir[2], dist, olddist);
1954  printf(" old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n",
1955  oldpoint[0],oldpoint[1],oldpoint[2], olddir[0], olddir[1], olddir[2]);
1956  if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1957  if (!pm1) {
1958  pm1 = new TPolyMarker3D();
1959  pm1->SetMarkerStyle(24);
1960  pm1->SetMarkerSize(0.4);
1961  pm1->SetMarkerColor(kRed);
1962  }
1963  pm1->SetNextPoint(point[0],point[1],point[2]);
1964  pm1->SetNextPoint(oldpoint[0],oldpoint[1],oldpoint[2]);
1965  break;
1966  }
1967  if ((dist<TGeoShape::Tolerance() && olddist*dot>1.E-3) || dist>dmax) {
1968  errsame++;
1969  if (errsame>1) {
1970  errcnt++;
1971  printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
1972  point[0],point[1],point[2], dir[0], dir[1], dir[2], dist, olddist);
1973  printf(" new norm: (%g, %g, %g)\n", newnorm[0], newnorm[1], newnorm[2]);
1974  printf(" old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n",
1975  oldpoint[0],oldpoint[1],oldpoint[2], olddir[0], olddir[1], olddir[2]);
1976  printf(" old norm: (%g, %g, %g)\n", oldnorm[0], oldnorm[1], oldnorm[2]);
1977  if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
1978  if (!pm1) {
1979  pm1 = new TPolyMarker3D();
1980  pm1->SetMarkerStyle(24);
1981  pm1->SetMarkerSize(0.4);
1982  pm1->SetMarkerColor(kRed);
1983  }
1984  pm1->SetNextPoint(point[0],point[1],point[2]);
1985  pm1->SetNextPoint(oldpoint[0],oldpoint[1],oldpoint[2]);
1986  break;
1987  }
1988  } else errsame = 0;
1989 
1990  olddist = dist;
1991  for (Int_t j=0; j<3; j++) {
1992  oldpoint[j] = point[j];
1993  point[j] += dist*dir[j];
1994  }
1995  safe = shape->Safety(point, kTRUE);
1996  if (safe>1.E-3) {
1997  errcnt++;
1998  printf("Error safety (%19.15f, %19.15f, %19.15f) safe=%g\n",
1999  point[0],point[1],point[2], safe);
2000  if (!errcanvas) errcanvas = new TCanvas("shape_err03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
2001  if (!pm1) {
2002  pm1 = new TPolyMarker3D();
2003  pm1->SetMarkerStyle(24);
2004  pm1->SetMarkerSize(0.4);
2005  pm1->SetMarkerColor(kRed);
2006  }
2007  pm1->SetNextPoint(point[0],point[1],point[2]);
2008  break;
2009  }
2010  // Compute normal
2011  shape->ComputeNormal(point,dir,norm);
2012  memcpy(oldnorm, norm, 3*sizeof(Double_t));
2013  memcpy(olddir, dir, 3*sizeof(Double_t));
2014  while (1) {
2015  phi = 2*TMath::Pi()*gRandom->Rndm();
2016  theta= TMath::ACos(1.-2.*gRandom->Rndm());
2017  dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
2018  dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
2019  dir[2]=TMath::Cos(theta);
2020  ndotd = dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2];
2021  if (ndotd<0) break; // backwards, still inside shape
2022  }
2023  if ((itot%10) == 0) pm2->SetNextPoint(point[0],point[1],point[2]);
2024  }
2025  }
2026  fTimer->Stop();
2027  fTimer->Print();
2028  if (errcanvas) {
2029  shape->Draw();
2030  pm1->Draw();
2031  }
2032 
2033  new TCanvas("shape03", Form("Shape %s (%s)",shape->GetName(),shape->ClassName()), 1000, 800);
2034  pm2->Draw();
2035  delete [] spoint;
2036  delete [] sdir;
2037 }
2038 
2039 ////////////////////////////////////////////////////////////////////////////////
2040 /// Generate a lego plot fot the top volume, according to option.
2041 
2042 TH2F *TGeoChecker::LegoPlot(Int_t ntheta, Double_t themin, Double_t themax,
2043  Int_t nphi, Double_t phimin, Double_t phimax,
2044  Double_t /*rmin*/, Double_t /*rmax*/, Option_t *option)
2045 {
2046  TH2F *hist = new TH2F("lego", option, nphi, phimin, phimax, ntheta, themin, themax);
2047 
2048  Double_t degrad = TMath::Pi()/180.;
2049  Double_t theta, phi, step, matprop, x;
2050  Double_t start[3];
2051  Double_t dir[3];
2052  TGeoNode *startnode, *endnode;
2053  Int_t i; // loop index for phi
2054  Int_t j; // loop index for theta
2055  Int_t ntot = ntheta * nphi;
2056  Int_t n10 = ntot/10;
2057  Int_t igen = 0, iloop=0;
2058  printf("=== Lego plot sph. => nrays=%i\n", ntot);
2059  for (i=1; i<=nphi; i++) {
2060  for (j=1; j<=ntheta; j++) {
2061  igen++;
2062  if (n10) {
2063  if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/ntot));
2064  }
2065  x = 0;
2066  theta = hist->GetYaxis()->GetBinCenter(j);
2067  phi = hist->GetXaxis()->GetBinCenter(i)+1E-3;
2068  start[0] = start[1] = start[2] = 1E-3;
2069  dir[0]=TMath::Sin(theta*degrad)*TMath::Cos(phi*degrad);
2070  dir[1]=TMath::Sin(theta*degrad)*TMath::Sin(phi*degrad);
2071  dir[2]=TMath::Cos(theta*degrad);
2072  fGeoManager->InitTrack(&start[0], &dir[0]);
2073  startnode = fGeoManager->GetCurrentNode();
2074  if (fGeoManager->IsOutside()) startnode=0;
2075  if (startnode) {
2076  matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
2077  } else {
2078  matprop = 0.;
2079  }
2080  fGeoManager->FindNextBoundary();
2081 // fGeoManager->IsStepEntering();
2082  // find where we end-up
2083  endnode = fGeoManager->Step();
2084  step = fGeoManager->GetStep();
2085  while (step<1E10) {
2086  // now see if we can make an other step
2087  iloop=0;
2088  while (!fGeoManager->IsEntering()) {
2089  iloop++;
2090  fGeoManager->SetStep(1E-3);
2091  step += 1E-3;
2092  endnode = fGeoManager->Step();
2093  }
2094  if (iloop>1000) printf("%i steps\n", iloop);
2095  if (matprop>0) {
2096  x += step/matprop;
2097  }
2098  if (endnode==0 && step>1E10) break;
2099  // generate an extra step to cross boundary
2100  startnode = endnode;
2101  if (startnode) {
2102  matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
2103  } else {
2104  matprop = 0.;
2105  }
2106 
2107  fGeoManager->FindNextBoundary();
2108  endnode = fGeoManager->Step();
2109  step = fGeoManager->GetStep();
2110  }
2111  hist->Fill(phi, theta, x);
2112  }
2113  }
2114  return hist;
2115 }
2116 
2117 ////////////////////////////////////////////////////////////////////////////////
2118 /// Draw random points in the bounding box of a volume.
2119 
2120 void TGeoChecker::RandomPoints(TGeoVolume *vol, Int_t npoints, Option_t *option)
2121 {
2122  if (!vol) return;
2123  vol->VisibleDaughters(kTRUE);
2124  vol->Draw();
2125  TString opt = option;
2126  opt.ToLower();
2127  TObjArray *pm = new TObjArray(128);
2128  TPolyMarker3D *marker = 0;
2129  const TGeoShape *shape = vol->GetShape();
2130  TGeoBBox *box = (TGeoBBox *)shape;
2131  Double_t dx = box->GetDX();
2132  Double_t dy = box->GetDY();
2133  Double_t dz = box->GetDZ();
2134  Double_t ox = (box->GetOrigin())[0];
2135  Double_t oy = (box->GetOrigin())[1];
2136  Double_t oz = (box->GetOrigin())[2];
2137  Double_t *xyz = new Double_t[3];
2138  printf("Random box : %f, %f, %f\n", dx, dy, dz);
2139  TGeoNode *node = 0;
2140  printf("Start... %i points\n", npoints);
2141  Int_t i=0;
2142  Int_t igen=0;
2143  Int_t ic = 0;
2144  Int_t n10 = npoints/10;
2145  Double_t ratio=0;
2146  while (igen<npoints) {
2147  xyz[0] = ox-dx+2*dx*gRandom->Rndm();
2148  xyz[1] = oy-dy+2*dy*gRandom->Rndm();
2149  xyz[2] = oz-dz+2*dz*gRandom->Rndm();
2150  fGeoManager->SetCurrentPoint(xyz);
2151  igen++;
2152  if (n10) {
2153  if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/npoints));
2154  }
2155  node = fGeoManager->FindNode();
2156  if (!node) continue;
2157  if (!node->IsOnScreen()) continue;
2158  // draw only points in overlapping/non-overlapping volumes
2159  if (opt.Contains("many") && !node->IsOverlapping()) continue;
2160  if (opt.Contains("only") && node->IsOverlapping()) continue;
2161  ic = node->GetColour();
2162  if ((ic<0) || (ic>=128)) ic = 1;
2163  marker = (TPolyMarker3D*)pm->At(ic);
2164  if (!marker) {
2165  marker = new TPolyMarker3D();
2166  marker->SetMarkerColor(ic);
2167 // marker->SetMarkerStyle(8);
2168 // marker->SetMarkerSize(0.4);
2169  pm->AddAt(marker, ic);
2170  }
2171  marker->SetNextPoint(xyz[0], xyz[1], xyz[2]);
2172  i++;
2173  }
2174  printf("Number of visible points : %i\n", i);
2175  if (igen) ratio = (Double_t)i/(Double_t)igen;
2176  printf("efficiency : %g\n", ratio);
2177  for (Int_t m=0; m<128; m++) {
2178  marker = (TPolyMarker3D*)pm->At(m);
2179  if (marker) marker->Draw("SAME");
2180  }
2181  fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
2182  printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
2183  printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
2184  delete pm;
2185  delete [] xyz;
2186 }
2187 
2188 ////////////////////////////////////////////////////////////////////////////////
2189 /// Randomly shoot nrays from point (startx,starty,startz) and plot intersections
2190 /// with surfaces for current top node.
2191 
2192 void TGeoChecker::RandomRays(Int_t nrays, Double_t startx, Double_t starty, Double_t startz, const char *target_vol, Bool_t check_norm)
2193 {
2194  TObjArray *pm = new TObjArray(128);
2195  TString starget = target_vol;
2196  TPolyLine3D *line = 0;
2197  TPolyLine3D *normline = 0;
2198  TGeoVolume *vol=fGeoManager->GetTopVolume();
2199 // vol->VisibleDaughters(kTRUE);
2200 
2201  Bool_t random = kFALSE;
2202  if (nrays<=0) {
2203  nrays = 100000;
2204  random = kTRUE;
2205  }
2206  Double_t start[3];
2207  Double_t dir[3];
2208  const Double_t *point = fGeoManager->GetCurrentPoint();
2209  vol->Draw();
2210  printf("Start... %i rays\n", nrays);
2211  TGeoNode *startnode, *endnode;
2212  Bool_t vis1,vis2;
2213  Int_t i=0;
2214  Int_t ipoint, inull;
2215  Int_t itot=0;
2216  Int_t n10=nrays/10;
2217  Double_t theta,phi, step, normlen;
2218  Double_t ox = ((TGeoBBox*)vol->GetShape())->GetOrigin()[0];
2219  Double_t oy = ((TGeoBBox*)vol->GetShape())->GetOrigin()[1];
2220  Double_t oz = ((TGeoBBox*)vol->GetShape())->GetOrigin()[2];
2221  Double_t dx = ((TGeoBBox*)vol->GetShape())->GetDX();
2222  Double_t dy = ((TGeoBBox*)vol->GetShape())->GetDY();
2223  Double_t dz = ((TGeoBBox*)vol->GetShape())->GetDZ();
2224  normlen = TMath::Max(dx,dy);
2225  normlen = TMath::Max(normlen,dz);
2226  normlen *= 0.05;
2227  while (itot<nrays) {
2228  itot++;
2229  inull = 0;
2230  ipoint = 0;
2231  if (n10) {
2232  if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nrays));
2233  }
2234  if (random) {
2235  start[0] = ox-dx+2*dx*gRandom->Rndm();
2236  start[1] = oy-dy+2*dy*gRandom->Rndm();
2237  start[2] = oz-dz+2*dz*gRandom->Rndm();
2238  } else {
2239  start[0] = startx;
2240  start[1] = starty;
2241  start[2] = startz;
2242  }
2243  phi = 2*TMath::Pi()*gRandom->Rndm();
2244  theta= TMath::ACos(1.-2.*gRandom->Rndm());
2245  dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
2246  dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
2247  dir[2]=TMath::Cos(theta);
2248  startnode = fGeoManager->InitTrack(start[0],start[1],start[2], dir[0],dir[1],dir[2]);
2249  line = 0;
2250  if (fGeoManager->IsOutside()) startnode=0;
2251  vis1 = kFALSE;
2252  if (target_vol) {
2253  if (startnode && starget==startnode->GetVolume()->GetName()) vis1 = kTRUE;
2254  } else {
2255  if (startnode && startnode->IsOnScreen()) vis1 = kTRUE;
2256  }
2257  if (vis1) {
2258  line = new TPolyLine3D(2);
2259  line->SetLineColor(startnode->GetVolume()->GetLineColor());
2260  line->SetPoint(ipoint++, start[0], start[1], start[2]);
2261  i++;
2262  pm->Add(line);
2263  }
2264  while ((endnode=fGeoManager->FindNextBoundaryAndStep())) {
2265  step = fGeoManager->GetStep();
2266  if (step<TGeoShape::Tolerance()) inull++;
2267  else inull = 0;
2268  if (inull>5) break;
2269  const Double_t *normal = 0;
2270  if (check_norm) {
2271  normal = fGeoManager->FindNormalFast();
2272  if (!normal) break;
2273  }
2274  vis2 = kFALSE;
2275  if (target_vol) {
2276  if (starget==endnode->GetVolume()->GetName()) vis2 = kTRUE;
2277  } else if (endnode->IsOnScreen()) vis2 = kTRUE;
2278  if (ipoint>0) {
2279  // old visible node had an entry point -> finish segment
2280  line->SetPoint(ipoint, point[0], point[1], point[2]);
2281  if (!vis2 && check_norm) {
2282  normline = new TPolyLine3D(2);
2283  normline->SetLineColor(kBlue);
2284  normline->SetLineWidth(1);
2285  normline->SetPoint(0, point[0], point[1], point[2]);
2286  normline->SetPoint(1, point[0]+normal[0]*normlen,
2287  point[1]+normal[1]*normlen,
2288  point[2]+normal[2]*normlen);
2289  pm->Add(normline);
2290  }
2291  ipoint = 0;
2292  line = 0;
2293  }
2294  if (vis2) {
2295  // create new segment
2296  line = new TPolyLine3D(2);
2297  line->SetLineColor(endnode->GetVolume()->GetLineColor());
2298  line->SetPoint(ipoint++, point[0], point[1], point[2]);
2299  i++;
2300  if (check_norm) {
2301  normline = new TPolyLine3D(2);
2302  normline->SetLineColor(kBlue);
2303  normline->SetLineWidth(2);
2304  normline->SetPoint(0, point[0], point[1], point[2]);
2305  normline->SetPoint(1, point[0]+normal[0]*normlen,
2306  point[1]+normal[1]*normlen,
2307  point[2]+normal[2]*normlen);
2308  }
2309  pm->Add(line);
2310  if (!random) pm->Add(normline);
2311  }
2312  }
2313  }
2314  // draw all segments
2315  for (Int_t m=0; m<pm->GetEntriesFast(); m++) {
2316  line = (TPolyLine3D*)pm->At(m);
2317  if (line) line->Draw("SAME");
2318  }
2319  printf("number of segments : %i\n", i);
2320 // fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
2321 // printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
2322 // printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
2323  delete pm;
2324 }
2325 
2326 ////////////////////////////////////////////////////////////////////////////////
2327 /// shoot npoints randomly in a box of 1E-5 around current point.
2328 /// return minimum distance to points outside
2329 /// make sure that path to current node is updated
2330 /// get the response of tgeo
2331 
2332 TGeoNode *TGeoChecker::SamplePoints(Int_t npoints, Double_t &dist, Double_t epsil,
2333  const char* g3path)
2334 {
2335  TGeoNode *node = fGeoManager->FindNode();
2336  TGeoNode *nodegeo = 0;
2337  TGeoNode *nodeg3 = 0;
2338  TGeoNode *solg3 = 0;
2339  if (!node) {dist=-1; return 0;}
2340  Bool_t hasg3 = kFALSE;
2341  if (strlen(g3path)) hasg3 = kTRUE;
2342  TString geopath = fGeoManager->GetPath();
2343  dist = 1E10;
2344  TString common = "";
2345  // cd to common path
2346  Double_t point[3];
2347  Double_t closest[3];
2348  TGeoNode *node1 = 0;
2349  TGeoNode *node_close = 0;
2350  dist = 1E10;
2351  Double_t dist1 = 0;
2352  // initialize size of random box to epsil
2353  Double_t eps[3];
2354  eps[0] = epsil; eps[1]=epsil; eps[2]=epsil;
2355  const Double_t *pointg = fGeoManager->GetCurrentPoint();
2356  if (hasg3) {
2357  TString spath = geopath;
2358  TString name = "";
2359  Int_t index=0;
2360  while (index>=0) {
2361  index = spath.Index("/", index+1);
2362  if (index>0) {
2363  name = spath(0, index);
2364  if (strstr(g3path, name.Data())) {
2365  common = name;
2366  continue;
2367  } else break;
2368  }
2369  }
2370  // if g3 response was given, cd to common path
2371  if (strlen(common.Data())) {
2372  while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
2373  nodegeo = fGeoManager->GetCurrentNode();
2374  fGeoManager->CdUp();
2375  }
2376  fGeoManager->cd(g3path);
2377  solg3 = fGeoManager->GetCurrentNode();
2378  while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
2379  nodeg3 = fGeoManager->GetCurrentNode();
2380  fGeoManager->CdUp();
2381  }
2382  if (!nodegeo) return 0;
2383  if (!nodeg3) return 0;
2384  fGeoManager->cd(common.Data());
2385  fGeoManager->MasterToLocal(fGeoManager->GetCurrentPoint(), &point[0]);
2386  Double_t xyz[3], local[3];
2387  for (Int_t i=0; i<npoints; i++) {
2388  xyz[0] = point[0] - eps[0] + 2*eps[0]*gRandom->Rndm();
2389  xyz[1] = point[1] - eps[1] + 2*eps[1]*gRandom->Rndm();
2390  xyz[2] = point[2] - eps[2] + 2*eps[2]*gRandom->Rndm();
2391  nodeg3->MasterToLocal(&xyz[0], &local[0]);
2392  if (!nodeg3->GetVolume()->Contains(&local[0])) continue;
2393  dist1 = TMath::Sqrt((xyz[0]-point[0])*(xyz[0]-point[0])+
2394  (xyz[1]-point[1])*(xyz[1]-point[1])+(xyz[2]-point[2])*(xyz[2]-point[2]));
2395  if (dist1<dist) {
2396  // save node and closest point
2397  dist = dist1;
2398  node_close = solg3;
2399  // make the random box smaller
2400  eps[0] = TMath::Abs(point[0]-pointg[0]);
2401  eps[1] = TMath::Abs(point[1]-pointg[1]);
2402  eps[2] = TMath::Abs(point[2]-pointg[2]);
2403  }
2404  }
2405  }
2406  if (!node_close) dist = -1;
2407  return node_close;
2408  }
2409 
2410  // save current point
2411  memcpy(&point[0], pointg, 3*sizeof(Double_t));
2412  for (Int_t i=0; i<npoints; i++) {
2413  // generate a random point in MARS
2414  fGeoManager->SetCurrentPoint(point[0] - eps[0] + 2*eps[0]*gRandom->Rndm(),
2415  point[1] - eps[1] + 2*eps[1]*gRandom->Rndm(),
2416  point[2] - eps[2] + 2*eps[2]*gRandom->Rndm());
2417  // check if new node is different from the old one
2418  if (node1!=node) {
2419  dist1 = TMath::Sqrt((point[0]-pointg[0])*(point[0]-pointg[0])+
2420  (point[1]-pointg[1])*(point[1]-pointg[1])+(point[2]-pointg[2])*(point[2]-pointg[2]));
2421  if (dist1<dist) {
2422  dist = dist1;
2423  node_close = node1;
2424  memcpy(&closest[0], pointg, 3*sizeof(Double_t));
2425  // make the random box smaller
2426  eps[0] = TMath::Abs(point[0]-pointg[0]);
2427  eps[1] = TMath::Abs(point[1]-pointg[1]);
2428  eps[2] = TMath::Abs(point[2]-pointg[2]);
2429  }
2430  }
2431  }
2432  // restore the original point and path
2433  fGeoManager->FindNode(point[0], point[1], point[2]); // really needed ?
2434  if (!node_close) dist=-1;
2435  return node_close;
2436 }
2437 
2438 ////////////////////////////////////////////////////////////////////////////////
2439 /// Shoot one ray from start point with direction (dirx,diry,dirz). Fills input array
2440 /// with points just after boundary crossings.
2441 
2442 Double_t *TGeoChecker::ShootRay(Double_t *start, Double_t dirx, Double_t diry, Double_t dirz, Double_t *array, Int_t &nelem, Int_t &dim, Double_t *endpoint) const
2443 {
2444 // Int_t array_dimension = 3*dim;
2445  nelem = 0;
2446  Int_t istep = 0;
2447  if (!dim) {
2448  printf("empty input array\n");
2449  return array;
2450  }
2451 // fGeoManager->CdTop();
2452  const Double_t *point = fGeoManager->GetCurrentPoint();
2453  TGeoNode *endnode;
2454  Bool_t is_entering;
2455  Double_t step, forward;
2456  Double_t dir[3];
2457  dir[0] = dirx;
2458  dir[1] = diry;
2459  dir[2] = dirz;
2460  fGeoManager->InitTrack(start, &dir[0]);
2461  fGeoManager->GetCurrentNode();
2462 // printf("Start : (%f,%f,%f)\n", point[0], point[1], point[2]);
2463  fGeoManager->FindNextBoundary();
2464  step = fGeoManager->GetStep();
2465 // printf("---next : at step=%f\n", step);
2466  if (step>1E10) return array;
2467  endnode = fGeoManager->Step();
2468  is_entering = fGeoManager->IsEntering();
2469  while (step<1E10) {
2470  if (endpoint) {
2471  forward = dirx*(endpoint[0]-point[0])+diry*(endpoint[1]-point[1])+dirz*(endpoint[2]-point[2]);
2472  if (forward<1E-3) {
2473 // printf("exit : Passed start point. nelem=%i\n", nelem);
2474  return array;
2475  }
2476  }
2477  if (is_entering) {
2478  if (nelem>=dim) {
2479  Double_t *temparray = new Double_t[3*(dim+20)];
2480  memcpy(temparray, array, 3*dim*sizeof(Double_t));
2481  delete [] array;
2482  array = temparray;
2483  dim += 20;
2484  }
2485  memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
2486 // printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
2487  nelem++;
2488  } else {
2489  if (endnode==0 && step>1E10) {
2490 // printf("exit : NULL endnode. nelem=%i\n", nelem);
2491  return array;
2492  }
2493  if (!fGeoManager->IsEntering()) {
2494 // if (startnode) printf("stepping %f from (%f, %f, %f) inside %s...\n", step,point[0], point[1], point[2], startnode->GetName());
2495 // else printf("stepping %f from (%f, %f, %f) OUTSIDE...\n", step,point[0], point[1], point[2]);
2496  istep = 0;
2497  }
2498  while (!fGeoManager->IsEntering()) {
2499  istep++;
2500  if (istep>1E3) {
2501 // Error("ShootRay", "more than 1000 steps. Step was %f", step);
2502  nelem = 0;
2503  return array;
2504  }
2505  fGeoManager->SetStep(1E-5);
2506  endnode = fGeoManager->Step();
2507  }
2508  if (istep>0) printf("%i steps\n", istep);
2509  if (nelem>=dim) {
2510  Double_t *temparray = new Double_t[3*(dim+20)];
2511  memcpy(temparray, array, 3*dim*sizeof(Double_t));
2512  delete [] array;
2513  array = temparray;
2514  dim += 20;
2515  }
2516  memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
2517 // if (endnode) printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
2518  nelem++;
2519  }
2520  fGeoManager->FindNextBoundary();
2521  step = fGeoManager->GetStep();
2522 // printf("---next at step=%f\n", step);
2523  endnode = fGeoManager->Step();
2524  is_entering = fGeoManager->IsEntering();
2525  }
2526  return array;
2527 // printf("exit : INFINITE step. nelem=%i\n", nelem);
2528 }
2529 
2530 ////////////////////////////////////////////////////////////////////////////////
2531 /// Check time of finding "Where am I" for n points.
2532 
2533 void TGeoChecker::Test(Int_t npoints, Option_t *option)
2534 {
2535  Bool_t recheck = !strcmp(option, "RECHECK");
2536  if (recheck) printf("RECHECK\n");
2537  const TGeoShape *shape = fGeoManager->GetTopVolume()->GetShape();
2538  Double_t dx = ((TGeoBBox*)shape)->GetDX();
2539  Double_t dy = ((TGeoBBox*)shape)->GetDY();
2540  Double_t dz = ((TGeoBBox*)shape)->GetDZ();
2541  Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
2542  Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
2543  Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
2544  Double_t *xyz = new Double_t[3*npoints];
2545  TStopwatch *timer = new TStopwatch();
2546  printf("Random box : %f, %f, %f\n", dx, dy, dz);
2547  timer->Start(kFALSE);
2548  Int_t i;
2549  for (i=0; i<npoints; i++) {
2550  xyz[3*i] = ox-dx+2*dx*gRandom->Rndm();
2551  xyz[3*i+1] = oy-dy+2*dy*gRandom->Rndm();
2552  xyz[3*i+2] = oz-dz+2*dz*gRandom->Rndm();
2553  }
2554  timer->Stop();
2555  printf("Generation time :\n");
2556  timer->Print();
2557  timer->Reset();
2558  TGeoNode *node, *node1;
2559  printf("Start... %i points\n", npoints);
2560  timer->Start(kFALSE);
2561  for (i=0; i<npoints; i++) {
2562  fGeoManager->SetCurrentPoint(xyz+3*i);
2563  if (recheck) fGeoManager->CdTop();
2564  node = fGeoManager->FindNode();
2565  if (recheck) {
2566  node1 = fGeoManager->FindNode();
2567  if (node1 != node) {
2568  printf("Difference for x=%g y=%g z=%g\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
2569  printf(" from top : %s\n", node->GetName());
2570  printf(" redo : %s\n", fGeoManager->GetPath());
2571  }
2572  }
2573  }
2574  timer->Stop();
2575  timer->Print();
2576  delete [] xyz;
2577  delete timer;
2578 }
2579 
2580 ////////////////////////////////////////////////////////////////////////////////
2581 /// Geometry overlap checker based on sampling.
2582 
2583 void TGeoChecker::TestOverlaps(const char* path)
2584 {
2585  if (fGeoManager->GetTopVolume()!=fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
2586  printf("Checking overlaps for path :\n");
2587  if (!fGeoManager->cd(path)) return;
2588  TGeoNode *checked = fGeoManager->GetCurrentNode();
2589  checked->InspectNode();
2590  // shoot 1E4 points in the shape of the current volume
2591  Int_t npoints = 1000000;
2592  Double_t big = 1E6;
2593  Double_t xmin = big;
2594  Double_t xmax = -big;
2595  Double_t ymin = big;
2596  Double_t ymax = -big;
2597  Double_t zmin = big;
2598  Double_t zmax = -big;
2599  TObjArray *pm = new TObjArray(128);
2600  TPolyMarker3D *marker = 0;
2601  TPolyMarker3D *markthis = new TPolyMarker3D();
2602  markthis->SetMarkerColor(5);
2603  TNtuple *ntpl = new TNtuple("ntpl","random points","x:y:z");
2604  TGeoShape *shape = fGeoManager->GetCurrentNode()->GetVolume()->GetShape();
2605  Double_t *point = new Double_t[3];
2606  Double_t dx = ((TGeoBBox*)shape)->GetDX();
2607  Double_t dy = ((TGeoBBox*)shape)->GetDY();
2608  Double_t dz = ((TGeoBBox*)shape)->GetDZ();
2609  Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
2610  Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
2611  Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
2612  Double_t *xyz = new Double_t[3*npoints];
2613  Int_t i=0;
2614  printf("Generating %i points inside %s\n", npoints, fGeoManager->GetPath());
2615  while (i<npoints) {
2616  point[0] = ox-dx+2*dx*gRandom->Rndm();
2617  point[1] = oy-dy+2*dy*gRandom->Rndm();
2618  point[2] = oz-dz+2*dz*gRandom->Rndm();
2619  if (!shape->Contains(point)) continue;
2620  // convert each point to MARS
2621 // printf("local %9.3f %9.3f %9.3f\n", point[0], point[1], point[2]);
2622  fGeoManager->LocalToMaster(point, &xyz[3*i]);
2623 // printf("master %9.3f %9.3f %9.3f\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
2624  xmin = TMath::Min(xmin, xyz[3*i]);
2625  xmax = TMath::Max(xmax, xyz[3*i]);
2626  ymin = TMath::Min(ymin, xyz[3*i+1]);
2627  ymax = TMath::Max(ymax, xyz[3*i+1]);
2628  zmin = TMath::Min(zmin, xyz[3*i+2]);
2629  zmax = TMath::Max(zmax, xyz[3*i+2]);
2630  i++;
2631  }
2632  delete [] point;
2633  ntpl->Fill(xmin,ymin,zmin);
2634  ntpl->Fill(xmax,ymin,zmin);
2635  ntpl->Fill(xmin,ymax,zmin);
2636  ntpl->Fill(xmax,ymax,zmin);
2637  ntpl->Fill(xmin,ymin,zmax);
2638  ntpl->Fill(xmax,ymin,zmax);
2639  ntpl->Fill(xmin,ymax,zmax);
2640  ntpl->Fill(xmax,ymax,zmax);
2641  ntpl->Draw("z:y:x");
2642 
2643  // shoot the poins in the geometry
2644  TGeoNode *node;
2645  TString cpath;
2646  Int_t ic=0;
2647  TObjArray *overlaps = new TObjArray();
2648  printf("using FindNode...\n");
2649  for (Int_t j=0; j<npoints; j++) {
2650  // always start from top level (testing only)
2651  fGeoManager->CdTop();
2652  fGeoManager->SetCurrentPoint(&xyz[3*j]);
2653  node = fGeoManager->FindNode();
2654  cpath = fGeoManager->GetPath();
2655  if (cpath.Contains(path)) {
2656  markthis->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
2657  continue;
2658  }
2659  // current point is found in an overlapping node
2660  if (!node) ic=128;
2661  else ic = node->GetVolume()->GetLineColor();
2662  if (ic >= 128) ic = 0;
2663  marker = (TPolyMarker3D*)pm->At(ic);
2664  if (!marker) {
2665  marker = new TPolyMarker3D();
2666  marker->SetMarkerColor(ic);
2667  pm->AddAt(marker, ic);
2668  }
2669  // draw the overlapping point
2670  marker->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
2671  if (node) {
2672  if (overlaps->IndexOf(node) < 0) overlaps->Add(node);
2673  }
2674  }
2675  // draw all overlapping points
2676 // for (Int_t m=0; m<128; m++) {
2677 // marker = (TPolyMarker3D*)pm->At(m);
2678 // if (marker) marker->Draw("SAME");
2679 // }
2680  markthis->Draw("SAME");
2681  if (gPad) gPad->Update();
2682  // display overlaps
2683  if (overlaps->GetEntriesFast()) {
2684  printf("list of overlapping nodes :\n");
2685  for (i=0; i<overlaps->GetEntriesFast(); i++) {
2686  node = (TGeoNode*)overlaps->At(i);
2687  if (node->IsOverlapping()) printf("%s MANY\n", node->GetName());
2688  else printf("%s ONLY\n", node->GetName());
2689  }
2690  } else printf("No overlaps\n");
2691  delete ntpl;
2692  delete pm;
2693  delete [] xyz;
2694  delete overlaps;
2695 }
2696 
2697 ////////////////////////////////////////////////////////////////////////////////
2698 /// Estimate weight of top level volume with a precision SIGMA(W)/W
2699 /// better than PRECISION. Option can be "v" - verbose (default).
2700 
2701 Double_t TGeoChecker::Weight(Double_t precision, Option_t *option)
2702 {
2703  TList *matlist = fGeoManager->GetListOfMaterials();
2704  Int_t nmat = matlist->GetSize();
2705  if (!nmat) return 0;
2706  Int_t *nin = new Int_t[nmat];
2707  memset(nin, 0, nmat*sizeof(Int_t));
2708  TString opt = option;
2709  opt.ToLower();
2710  Bool_t isverbose = opt.Contains("v");
2711  TGeoBBox *box = (TGeoBBox *)fGeoManager->GetTopVolume()->GetShape();
2712  Double_t dx = box->GetDX();
2713  Double_t dy = box->GetDY();
2714  Double_t dz = box->GetDZ();
2715  Double_t ox = (box->GetOrigin())[0];
2716  Double_t oy = (box->GetOrigin())[1];
2717  Double_t oz = (box->GetOrigin())[2];
2718  Double_t x,y,z;
2719  TGeoNode *node;
2720  TGeoMaterial *mat;
2721  Double_t vbox = 0.000008*dx*dy*dz; // m3
2722  Bool_t end = kFALSE;
2723  Double_t weight=0, sigma, eps, dens;
2724  Double_t eps0=1.;
2725  Int_t indmat;
2726  Int_t igen=0;
2727  Int_t iin = 0;
2728  while (!end) {
2729  x = ox-dx+2*dx*gRandom->Rndm();
2730  y = oy-dy+2*dy*gRandom->Rndm();
2731  z = oz-dz+2*dz*gRandom->Rndm();
2732  node = fGeoManager->FindNode(x,y,z);
2733  igen++;
2734  if (!node) continue;
2735  mat = node->GetVolume()->GetMedium()->GetMaterial();
2736  indmat = mat->GetIndex();
2737  if (indmat<0) continue;
2738  nin[indmat]++;
2739  iin++;
2740  if ((iin%100000)==0 || igen>1E8) {
2741  weight = 0;
2742  sigma = 0;
2743  for (indmat=0; indmat<nmat; indmat++) {
2744  mat = (TGeoMaterial*)matlist->At(indmat);
2745  dens = mat->GetDensity(); // [g/cm3]
2746  if (dens<1E-2) dens=0;
2747  dens *= 1000.; // [kg/m3]
2748  weight += dens*Double_t(nin[indmat]);
2749  sigma += dens*dens*nin[indmat];
2750  }
2751  sigma = TMath::Sqrt(sigma);
2752  eps = sigma/weight;
2753  weight *= vbox/Double_t(igen);
2754  sigma *= vbox/Double_t(igen);
2755  if (eps<precision || igen>1E8) {
2756  if (isverbose) {
2757  printf("=== Weight of %s : %g +/- %g [kg]\n",
2758  fGeoManager->GetTopVolume()->GetName(), weight, sigma);
2759  }
2760  end = kTRUE;
2761  } else {
2762  if (isverbose && eps<0.5*eps0) {
2763  printf("%8dK: %14.7g kg %g %%\n",
2764  igen/1000, weight, eps*100);
2765  eps0 = eps;
2766  }
2767  }
2768  }
2769  }
2770  delete [] nin;
2771  return weight;
2772 }
2773 ////////////////////////////////////////////////////////////////////////////////
2774 /// count voxel timing
2775 
2776 Double_t TGeoChecker::CheckVoxels(TGeoVolume *vol, TGeoVoxelFinder *voxels, Double_t *xyz, Int_t npoints)
2777 {
2778  TStopwatch timer;
2779  Double_t time;
2780  TGeoShape *shape = vol->GetShape();
2781  TGeoNode *node;
2782  TGeoMatrix *matrix;
2783  Double_t *point;
2784  Double_t local[3];
2785  Int_t *checklist;
2786  Int_t ncheck;
2787  TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
2788  TGeoStateInfo &td = *nav->GetCache()->GetInfo();
2789  timer.Start();
2790  for (Int_t i=0; i<npoints; i++) {
2791  point = xyz + 3*i;
2792  if (!shape->Contains(point)) continue;
2793  checklist = voxels->GetCheckList(point, ncheck, td);
2794  if (!checklist) continue;
2795  if (!ncheck) continue;
2796  for (Int_t id=0; id<ncheck; id++) {
2797  node = vol->GetNode(checklist[id]);
2798  matrix = node->GetMatrix();
2799  matrix->MasterToLocal(point, &local[0]);
2800  if (node->GetVolume()->GetShape()->Contains(&local[0])) break;
2801  }
2802  }
2803  nav->GetCache()->ReleaseInfo();
2804  time = timer.CpuTime();
2805  return time;
2806 }
2807 
2808 ////////////////////////////////////////////////////////////////////////////////
2809 /// Returns optimal voxelization type for volume vol.
2810 /// - kFALSE - cartesian
2811 /// - kTRUE - cylindrical
2812 
2813 Bool_t TGeoChecker::TestVoxels(TGeoVolume * /*vol*/, Int_t /*npoints*/)
2814 {
2815  return kFALSE;
2816 }