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