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