11 #include <simulation/modules/MaterialScan.h>
12 #include <framework/core/ModuleParam.templateDetails.h>
13 #include <framework/gearbox/Unit.h>
14 #include <framework/utilities/Utils.h>
23 #include <G4EventManager.hh>
24 #include <G4UserEventAction.hh>
25 #include <G4UserStackingAction.hh>
26 #include <G4UserTrackingAction.hh>
27 #include <G4UserSteppingAction.hh>
28 #include <G4RayShooter.hh>
31 #include <boost/format.hpp>
32 #include <boost/algorithm/string.hpp>
46 bool MaterialScanBase::checkStep(
const G4Step* step)
48 double stlen = step->GetStepLength();
49 G4StepPoint* preStepPoint = step->GetPreStepPoint();
50 G4Region* region = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
51 if (stlen < c_zeroTolerance) {
56 if (m_zeroSteps > c_maxZeroStepsNudge) {
57 if (m_zeroSteps > c_maxZeroStepsKill) {
58 B2ERROR(
"Track is stuck at " << preStepPoint->GetPosition() <<
" in volume '"
59 << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
60 <<
" (" << region->GetName() <<
"): "
61 << m_zeroSteps <<
" consecutive steps with length less then "
62 << c_zeroTolerance <<
" mm, killing it");
63 step->GetTrack()->SetTrackStatus(fStopAndKill);
65 B2WARNING(
"Track is stuck at " << preStepPoint->GetPosition() <<
" in volume '"
66 << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
67 <<
" (" << region->GetName() <<
"): "
68 << m_zeroSteps <<
" consecutive steps with length less then "
69 << c_zeroTolerance <<
" mm, nudging it along");
70 G4ThreeVector pos = step->GetTrack()->GetPosition();
71 G4ThreeVector dir = step->GetTrack()->GetMomentumDirection();
72 step->GetTrack()->SetPosition(pos + c_zeroTolerance * dir);
80 MaterialScan2D::MaterialScan2D(TFile* rootFile,
const std::string& name,
const std::string& axisLabel,
const ScanParams& params):
121 getRay(origin, direction);
129 std::unique_ptr<TH2D>& hist =
m_regions[name];
133 hist.reset(
new TH2D(name.c_str(), (name +
";" +
m_axisLabel).c_str(),
149 G4StepPoint* preStepPoint = step->GetPreStepPoint();
150 G4Region* region = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
151 G4Material* material = preStepPoint->GetMaterial();
152 double stlen = step->GetStepLength();
153 double x0 = stlen / (material->GetRadlen());
154 double lambda = stlen / (material->GetNuclearInterLength());
155 B2DEBUG(20,
"Step in at " << preStepPoint->GetPosition() <<
" in volume '"
156 << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
157 <<
" (" << region->GetName() <<
")"
158 <<
" with length=" << stlen <<
" mm");
166 G4Track* track = step->GetTrack();
167 track->SetTrackStatus(fStopAndKill);
169 x0 *= remaining / stlen;
170 lambda *= remaining / stlen;
179 string x0_total =
"All_Regions_x0";
180 string lambda_total =
"All_Regions_lambda";
181 string x0_name = region->GetName() +
"_x0";
182 string lambda_name = region->GetName() +
"_lambda";
185 x0_total =
"All_Materials_x0";
186 lambda_total =
"All_Materials_lambda";
187 x0_name = material->GetName() +
"_x0";
188 lambda_name = material->GetName() +
"_lambda";
205 direction.set(sin(theta) * cos(phi),
206 sin(theta) * sin(phi),
230 B2INFO(
"Set maximum depth to " <<
m_maxDepth);
245 double theta = acos(gRandom->Uniform(cos(
m_opening / 2.), 1));
246 double phi = gRandom->Uniform(0, 2 * M_PI);
247 auto orthogonal = direction.orthogonal();
248 orthogonal.rotate(phi, direction);
249 direction.rotate(theta, orthogonal);
251 B2DEBUG(10,
"Create Ray " <<
m_curRay <<
" from " << origin <<
" direction " <<
252 direction <<
" [theta=" << (direction.theta() / M_PI * 180) <<
", phi=" << (direction.phi() / M_PI * 180) <<
"]");
258 std::unique_ptr<TH1D>& hist =
m_regions[name];
276 const double dx = steplength / samples;
277 const double dv = value / samples /
m_count / h->GetBinWidth(1);
278 for (
int i = 0; i < samples; ++i) {
287 G4StepPoint* preStepPoint = step->GetPreStepPoint();
288 G4Region* region = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
289 G4Material* material = preStepPoint->GetMaterial();
291 double stlen = step->GetStepLength();
292 double x0 = stlen / (material->GetRadlen());
293 double lambda = stlen / (material->GetNuclearInterLength());
294 B2DEBUG(20,
"Step in at " << preStepPoint->GetPosition() <<
" in volume '"
295 << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
296 <<
" (" << region->GetName() <<
")"
297 <<
" with length=" << stlen <<
" mm");
304 G4Track* track = step->GetTrack();
305 track->SetTrackStatus(fStopAndKill);
307 x0 *= remaining / stlen;
308 lambda *= remaining / stlen;
315 string x0_total =
"All_Regions_x0";
316 string lambda_total =
"All_Regions_lambda";
317 string x0_name = region->GetName() +
"_x0";
318 string lambda_name = region->GetName() +
"_lambda";
321 x0_total =
"All_Materials_x0";
322 lambda_total =
"All_Materials_lambda";
323 x0_name = material->GetName() +
"_x0";
324 lambda_name = material->GetName() +
"_lambda";
338 MaterialScanModule::MaterialScanModule(): m_rootFile(0), m_sphericalOrigin(3, 0), m_doCosTheta(false)
341 setDescription(
"This Module is intended to scan the material budget of the "
342 "geometry. Currently, there are two different kinds of scans "
343 "available: Spherical and Planar scan. \n"
344 "Spherical scan will shoot rays from the origin of the "
345 "detector and scan along the polar and azimuth angle.\n"
346 "Planar scan will shoot rays perpendicular to a given "
350 m_spherical.ignoredMaterials.push_back(
"Vacuum");
351 m_spherical.ignoredMaterials.push_back(
"Air");
352 m_spherical.ignoredMaterials.push_back(
"G4_AIR");
353 m_planar.ignoredMaterials = m_spherical.ignoredMaterials;
354 m_rayIgnoredMaterials = m_spherical.ignoredMaterials;
356 addParam(
"Filename", m_filename,
357 "The filename where the material scan will be stored",
358 string(
"MaterialScan.root"));
359 addParam(
"spherical", m_doSpherical,
360 "Do a spherical scan, that is shooting rays from the origin with "
361 "varying angles",
true);
362 addParam(
"spherical.origin", m_sphericalOrigin,
363 "Origin for the spherical scan", m_sphericalOrigin);
364 addParam(
"spherical.nTheta", m_spherical.nU,
365 "Number of rays in theta", 200);
366 addParam(
"spherical.minTheta", m_spherical.minU,
367 "Theta start angle", 17.);
368 addParam(
"spherical.maxTheta", m_spherical.maxU,
369 "Theta stop angle", 150.);
370 addParam(
"spherical.nPhi", m_spherical.nV,
371 "Number of rays in phi", 200);
372 addParam(
"spherical.minPhi", m_spherical.minV,
373 "Phi start angle", 0.);
374 addParam(
"spherical.maxPhi", m_spherical.maxV,
375 "Phi stop angle", 360.);
376 addParam(
"spherical.maxDepth", m_spherical.maxDepth,
377 "Maximum scan depth in cm. The ray will be killed after having "
378 "reached the maximum Depth. <=0 means no Limit.", -1.0);
379 addParam(
"spherical.ignored", m_spherical.ignoredMaterials,
380 "Names of Materials which should be ignored when doing the scan", m_spherical.ignoredMaterials);
381 addParam(
"spherical.splitByMaterials", m_spherical.splitByMaterials,
382 "If True, split output by material names instead of by regions",
false);
383 addParam(
"spherical.cosTheta", m_doCosTheta,
384 "If True, perform the spherical scan uniform in cos(theta) instead of theta",
false);
386 addParam(
"planar", m_doPlanar,
387 "Do a plane scan, that is shooting parallel rays from a defined "
389 addParam(
"planar.plane", m_planeName,
390 "Plane to use for scanning, available are all two letter "
391 "combinations of X,Y and Z like XY and XZ and custom",
string(
"ZX"));
392 addParam(
"planar.nU", m_planar.nU,
393 "Number of rays in U", 200);
394 addParam(
"planar.minU", m_planar.minU,
396 addParam(
"planar.maxU", m_planar.maxU,
398 addParam(
"planar.nV", m_planar.nV,
399 "Number of rays in V", 200);
400 addParam(
"planar.minV", m_planar.minV,
402 addParam(
"planar.maxV", m_planar.maxV,
404 addParam(
"planar.maxDepth", m_planar.maxDepth,
405 "Maximum scan depth in cm. The ray will be killed after having "
406 "reached the maximum Depth. <=0 means no Limit.", -1.0);
407 addParam(
"planar.custom", m_customPlane,
408 "Parameters of the plane when choosing custom. This is supposed to "
409 "be a list of 9 values: The first three are the coordinates of the "
410 "plane origin. The second three are the direction of U in the "
411 "r-phi coordinates. The last three are the direction of V in the "
412 "coordinates parallel to the detector axis (beamline). ", m_customPlane);
413 addParam(
"planar.ignored", m_planar.ignoredMaterials,
414 "Names of Materials which should be ignored when doing the scan", m_planar.ignoredMaterials);
415 addParam(
"planar.splitByMaterials", m_planar.splitByMaterials,
416 "If True, split output by material names instead of by regions",
false);
418 addParam(
"ray", m_doRay,
"Do a ray scan: Shoot from the given origin in a "
419 "given direction and record the material encountered along the way",
false);
420 addParam(
"ray.origin", m_rayOrigin,
"Origin for the ray", m_rayOrigin);
421 addParam(
"ray.theta", m_rayTheta,
"Theta angle for the ray in degrees", m_rayTheta);
422 addParam(
"ray.phi", m_rayPhi,
"Phi angle of the ray in degrees", m_rayPhi);
423 addParam(
"ray.direction", m_rayDirection,
"Set the ray direction as vector "
424 "[x,y,z]. You cannot set this and theta/phi at the same time",
426 addParam(
"ray.opening", m_rayOpening,
"Opening angle for the ray in degrees. "
427 "If bigger than 0 then multiple rays (``ray.count``) will be shot "
428 "randomly inside the opening angle uniformly distributed in area "
429 "(cosine of the angle)", m_rayOpening);
430 addParam(
"ray.count", m_rayCount,
"Count of rays if the opening angle is >0",
432 addParam(
"ray.sampleDepth", m_raySampleDepth,
"Distance in which to sample the "
433 "material", m_raySampleDepth);
434 addParam(
"ray.maxDepth", m_rayMaxDepth,
"Max distance for the ray before it's "
435 "stopped, 0 for no limit. The actual depth can be smaller if the "
436 "simulation volume ends earlier", m_rayMaxDepth);
437 addParam(
"ray.ignored", m_rayIgnoredMaterials,
"Names of Materials which "
438 "should be ignored when doing the scan", m_rayIgnoredMaterials);
445 B2ERROR(
"Cannot open rootfile for writing" <<
LogVar(
"rootfile",
m_filename));
454 B2ERROR(
"planar.custom: At least 9 values needed to define custom plane, only " <<
m_customPlane.size() <<
" provided");
457 B2ERROR(
"planar.plane: Only custom or any two letter combinations of X,Y and Z are recognized");
459 B2ERROR(
"planar.plane: " <<
m_planeName <<
" not valid, cannot use the same axis twice");
464 B2ERROR(
"spherical.origin: Three values are needed to define a point, only " <<
m_sphericalOrigin.size() <<
" given.");
472 if (getParam<double>(
"ray.theta").isSetInSteering() || getParam<double>(
"ray.phi").isSetInSteering()) {
473 B2ERROR(
"Cannot set ray.theta/ray.phi and ray.direction at the same time, please only set one");
482 G4EventManager* eventManager = G4EventManager::GetEventManager();
485 G4UserEventAction* vanillaEventAction = eventManager->GetUserEventAction();
486 G4UserStackingAction* vanillaStackingAction = eventManager->GetUserStackingAction();
487 G4UserTrackingAction* vanillaTrackingAction = eventManager->GetUserTrackingAction();
488 G4UserSteppingAction* vanillaSteppingAction = eventManager->GetUserSteppingAction();
491 eventManager->SetUserAction((G4UserEventAction*)0);
492 eventManager->SetUserAction((G4UserStackingAction*)0);
493 eventManager->SetUserAction((G4UserTrackingAction*)0);
498 vector<MaterialScanBase*> scans;
504 G4ThreeVector origin(0, 0, 0);
535 eventManager->SetUserAction(scan);
537 G4RayShooter rayShooter;
538 G4ThreeVector origin;
539 G4ThreeVector direction;
540 int maxRays = scan->getNRays();
545 while (scan->createNext(origin, direction)) {
546 G4Event*
event =
new G4Event(++curRay);
547 rayShooter.Shoot(
event, origin, direction);
548 eventManager->ProcessOneEvent(
event);
552 int donePercent = 100 * curRay / maxRays;
553 if (donePercent > lastPercent) {
555 double eta = perRay * (maxRays - curRay);
556 B2INFO(boost::format(
"%s Scan: %3d%%, %.3f ms per ray, ETA: %.2f seconds")
557 % scan->getName() % donePercent
559 lastPercent = donePercent;
573 eventManager->SetUserAction(vanillaEventAction);
574 eventManager->SetUserAction(vanillaStackingAction);
575 eventManager->SetUserAction(vanillaTrackingAction);
576 eventManager->SetUserAction(vanillaSteppingAction);
582 case 'x':
return G4ThreeVector(1, 0, 0);
583 case 'y':
return G4ThreeVector(0, 1, 0);
584 case 'z':
return G4ThreeVector(0, 0, 1);
586 B2FATAL(
"Unknown axis: " << name);
588 return G4ThreeVector(0, 0, 0);