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