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