19#include "BDSAcceleratorComponent.hh"
20#include "BDSAcceleratorComponentRegistry.hh"
21#include "BDSAcceleratorModel.hh"
22#include "BDSApertureInfo.hh"
23#include "BDSAuxiliaryNavigator.hh"
24#include "BDSBeamline.hh"
25#include "BDSBeamlineBLMBuilder.hh"
26#include "BDSBeamlineEndPieceBuilder.hh"
27#include "BDSBeamlineElement.hh"
28#include "BDSBeamlinePlacementBuilder.hh"
29#include "BDSBeamlineSet.hh"
30#include "BDSBeamPipeInfo.hh"
32#include "BDSBLMRegistry.hh"
33#include "BDSBOptrMultiParticleChangeCrossSection.hh"
34#include "BDSComponentFactory.hh"
35#include "BDSComponentFactoryUser.hh"
36#include "BDSCurvilinearBuilder.hh"
38#include "BDSDetectorConstruction.hh"
39#include "BDSException.hh"
40#include "BDSExtent.hh"
41#include "BDSFieldBuilder.hh"
42#include "BDSFieldObjects.hh"
43#include "BDSFieldQuery.hh"
44#include "BDSFieldQueryInfo.hh"
45#include "BDSFieldLoaderQueryPoints.hh"
47#include "BDSGeometryComponent.hh"
48#include "BDSGeometryExternal.hh"
49#include "BDSGeometryFactory.hh"
50#include "BDSGlobalConstants.hh"
51#include "BDSHistBinMapper.hh"
52#include "BDSIntegratorSet.hh"
54#include "BDSMaterials.hh"
55#include "BDSParser.hh"
56#include "BDSPhysicalVolumeInfo.hh"
57#include "BDSPhysicalVolumeInfoRegistry.hh"
58#include "BDSRegion.hh"
59#include "BDSSamplerInfo.hh"
60#include "BDSSamplerType.hh"
61#include "BDSScorerFactory.hh"
62#include "BDSScorerInfo.hh"
63#include "BDSScorerMeshInfo.hh"
64#include "BDSScoringMeshBox.hh"
65#include "BDSScoringMeshCylinder.hh"
66#include "BDSSDEnergyDeposition.hh"
67#include "BDSSDManager.hh"
68#include "BDSSDType.hh"
69#include "BDSSurvey.hh"
70#include "BDSTeleporter.hh"
71#include "BDSTrajectoryPoint.hh"
72#include "BDSTunnelBuilder.hh"
73#include "BDSUtilities.hh"
74#include "BDSWarning.hh"
76#include "parser/blmplacement.h"
77#include "parser/element.h"
78#include "parser/fastlist.h"
79#include "parser/options.h"
80#include "parser/physicsbiasing.h"
81#include "parser/placement.h"
82#include "parser/samplerplacement.h"
83#include "parser/scorermesh.h"
86#include "G4AffineTransform.hh"
88#include "G4LogicalVolume.hh"
89#include "G4Material.hh"
90#include "G4ProductionCuts.hh"
91#include "G4PVPlacement.hh"
92#include "G4VPrimitiveScorer.hh"
94#include "G4ScoringManager.hh"
96#include "G4Transform3D.hh"
97#include "G4Version.hh"
98#include "G4VisAttributes.hh"
99#include "G4VPhysicalVolume.hh"
100#if G4VERSION_NUMBER > 1039
101#include "G4ChannelingOptrMultiParticleChangeCrossSection.hh"
104#ifdef BDSCHECKUSERLIMITS
105#include "G4UserLimits.hh"
108#include "CLHEP/Units/SystemOfUnits.h"
109#include "CLHEP/Vector/EulerAngles.h"
123 placementBL(nullptr),
124 designParticle(nullptr),
125 userComponentFactory(userComponentFactoryIn),
127 buildPlacementFieldsWorld(false),
128 worldLogicalVolume(nullptr)
131 verbose = globals->Verbose();
132 checkOverlaps = globals->CheckOverlaps();
133 circular = globals->Circular();
138 acceleratorModel = BDSAcceleratorModel::Instance();
139 canSampleAngledFaces =
true;
141 if ( (integratorSetType == BDSIntegratorSetType::bdsimtwo)
142 || (integratorSetType == BDSIntegratorSetType::geant4)
143#
if G4VERSION_NUMBER > 1039
144 || (integratorSetType == BDSIntegratorSetType::geant4dp)
148 canSampleAngledFaces = globals->SampleElementsWithPoleface();
151 UpdateSamplerDiameterAndCountSamplers();
152 PrepareExtraSamplerSDs();
153 CountPlacementFields();
160 G4double maxBendingRatio = 1e-9;
161 for (
const auto& blElement : beamline)
164 auto st = BDS::DetermineSamplerType(blElement.samplerType);
165 if (st != BDSSamplerType::none)
168 G4double length = blElement.l;
169 G4double angle = blElement.angle;
172 G4double ratio = angle / length;
173 maxBendingRatio = std::max(maxBendingRatio, ratio);
177 G4double curvilinearRadius = 0.5*globals->CurvilinearDiameter();
178 G4double tolerance = 0.9;
179 if (maxBendingRatio > 0.4*tolerance)
181 G4double curvilinearRadiusBends = (tolerance / maxBendingRatio)*CLHEP::m;
182 if (curvilinearRadiusBends < curvilinearRadius)
184 G4cout << __METHOD_NAME__ <<
"Reducing curvilinear diameter from " << 2*curvilinearRadius / CLHEP::m
185 <<
"m to " << 2*curvilinearRadiusBends / CLHEP::m <<
"m" << G4endl;
189 G4double sd = globals->SamplerDiameter();
190 if (curvilinearRadius*2 < sd)
192 G4cout << __METHOD_NAME__ <<
"Reducing sampler diameter from " << sd / CLHEP::m <<
"m to the same" << G4endl;
211 for (
const auto& placement : placements)
213 if (!placement.fieldAll.empty() || !placement.bdsimElement.empty())
216 buildPlacementFieldsWorld = nFields > 0;
222 {G4cout << __METHOD_NAME__ <<
"starting accelerator geometry construction\n" << G4endl;}
240 delete componentFactory;
258 {G4cout << __METHOD_NAME__ <<
"detector Construction done" << G4endl;}
263 G4cout << G4endl << __METHOD_NAME__ <<
"printing material table" << G4endl;
264 G4cout << *(G4Material::GetMaterialTable()) << G4endl << G4endl;
265 if (
verbose || debug) {G4cout <<
"Finished listing materials, returning physiWorld" << G4endl;}
267#ifdef BDSCHECKUSERLIMITS
268 PrintUserLimitsSummary(worldPV);
273BDSDetectorConstruction::~BDSDetectorConstruction()
275#if G4VERSION_NUMBER > 1009
280 for (
auto q : fieldQueries)
290 G4cout <<
"New region defined: " << G4endl << *reg << G4endl;
293 delete defaultRegion;
298 std::map<G4String, BDSApertureInfo*> apertures;
307 apertures[a.name] = ap;
316 {G4cout <<
"parsing the beamline element list..."<< G4endl;}
327 G4cout <<
"Registry size "
329 G4cout <<
"Parser beam line size "
335 if (!
circular && mainBeamline.massWorld)
338 {BDS::Warning(
"Total sum of all element angles is approximately 2*pi but the circular option was not specified, this simulation may run indefinitely");}
347 for (
const auto& placement : placements)
349 if (placement.sequence.empty())
364 G4String beamlineName = placement.name +
"_" + placement.sequence;
379 if (sType == BDSSamplerType::none)
388 const G4String& name,
389 const G4Transform3D& initialTransform,
391 G4bool beamlineIsCircular,
392 G4bool isPlacementBeamline)
394 if (beamLine.
empty())
397 if (userComponentFactory)
402 if (beamlineIsCircular)
407 G4cerr <<
"The first element in the beam line is unsuitable for a circular "
408 <<
"model as the first element will " << G4endl <<
"overlap with the "
409 <<
"teleporter and terminator - the necessary mechanics for a circular "
410 <<
"model in Geant4" << G4endl;
411 throw BDSException(__METHOD_NAME__,
"check construction for circular machine");
415 if (beamLine.
size() <= 1)
416 {
throw BDSException(__METHOD_NAME__,
"BDSIM requires the sequence defined with the use command to have at least one element.");}
418 for (
auto elementIt = beamLine.
begin(); elementIt != beamLine.
end(); ++elementIt)
422 auto prevIt = elementIt;
423 while (prevIt != beamLine.
begin())
428 prevElement = &(*prevIt);
434 auto nextIt = elementIt;
436 G4double nextElementInputFace = 0;
437 while (nextIt != beamLine.
end())
441 nextElement = &(*nextIt);
443 nextElementInputFace = nextElement->
e1;
455 G4bool forceNoSamplerOnThisElement =
false;
457 {forceNoSamplerOnThisElement =
true;}
459 {forceNoSamplerOnThisElement =
true;}
461 {forceNoSamplerOnThisElement =
true;}
471 if (beamlineIsCircular && !massWorld->
empty())
474 G4cout << __METHOD_NAME__ <<
"Circular machine - creating terminator & teleporter" << G4endl;
476 G4double teleporterLength = 0;
479 auto hasBeamPipeInfo = [](
BDSBeamlineElement* ble) {
return ble->GetBeamPipeInfo() !=
nullptr;};
480 auto firstElementWithBPInfo = std::find_if(massWorld->
begin(), massWorld->
end(), hasBeamPipeInfo);
481 auto lastElementWithBPInfo = std::find_if(massWorld->
rbegin(), massWorld->
rend(), hasBeamPipeInfo);
483 G4double firstbeamPipeMaxExtent = (*firstElementWithBPInfo)->GetBeamPipeInfo()->Extent().MaximumAbsTransverse();
484 G4double lastbeamPipeMaxExtent = (*lastElementWithBPInfo)->GetBeamPipeInfo()->Extent().MaximumAbsTransverse();
487 G4double teleporterHorizontalWidth = 2 * std::max(firstbeamPipeMaxExtent, lastbeamPipeMaxExtent);
497 teleporterHorizontalWidth,
498 teleporterTransform);
509 if (isPlacementBeamline)
512 survey->
Write(massWorld);
515 delete theComponentFactory;
518 G4cout << __METHOD_NAME__ <<
"\"" << name <<
"\" " << G4endl << *massWorld;
529 beamlineSet.massWorld = massWorld;
530 beamlineSet.curvilinearWorld = clBeamline;
531 beamlineSet.curvilinearBridgeWorld = clBridgeBeamline;
532 beamlineSet.endPieces = endPieces;
550 std::vector<BDSExtentGlobal> extents;
567 const auto& extras = BDSAcceleratorModel::Instance()->
ExtraBeamlines();
569 for (
const auto& bl : extras)
570 {bl.second.GetExtentGlobals(extents);}
579 G4ThreeVector worldR;
581 for (
const auto& ext : extents)
583 for (G4int i = 0; i < 3; i++)
584 {worldR[i] = std::max(worldR[i], ext.GetMaximumExtentAbsolute()[i]);}
587 G4String worldName =
"World";
588 G4VSolid* worldSolid =
nullptr;
589 G4LogicalVolume* worldLV =
nullptr;
592 if (!worldGeometryFile.empty())
595 {BDS::Warning(__METHOD_NAME__,
"conflicting options - world material option specified but material will be taken from world GDML file");}
605 &namedWorldVacuumVolumes,
607 BDSSDType::energydepworldcontents);
613 for (
auto volume : worldVacuumLogicalVolumes)
614 {allWorldVolumes.erase(volume);}
619 G4bool worldContainsAllBeamlines = worldExtentGlobal.
Encompasses(extents);
621 G4cout <<
"External world geometry: \"" << worldGeometryFile <<
"\"" << G4endl;
622 G4cout <<
"Loaded world extent: \n" <<
worldExtent << G4endl;
625 if (!worldContainsAllBeamlines)
627 G4String message =
"Beamlines cannot be constructed, beamline extents are larger than \n";
628 message +=
"the extents of the external world";
644 worldLV->SetSensitiveDetector(BDSSDManager::Instance()->WorldComplete());
651 G4cout << __METHOD_NAME__ <<
"world extent absolute: " << worldR << G4endl;
654 margin = std::max(margin, 2*CLHEP::m);
655 worldR += G4ThreeVector(margin,margin,margin);
657 G4cout << __METHOD_NAME__ <<
"with " << margin <<
"m margin, it becomes in all dimensions: " << worldR << G4endl;
659 G4cout << __METHOD_NAME__ <<
"World dimensions: " << worldR / CLHEP::m <<
" m" << G4endl;
664 worldSolid =
new G4Box(worldName +
"_solid", worldR.x(), worldR.y(), worldR.z());
667 worldLV =
new G4LogicalVolume(worldSolid,
673 {worldLV->SetSensitiveDetector(BDSSDManager::Instance()->WorldComplete());}
679 debugWorldVis->SetForceWireframe(
true);
680 worldLV->SetVisAttributes(debugWorldVis);
686 G4VPhysicalVolume* worldPV =
new G4PVPlacement(
nullptr,
709 worldLogicalVolume = worldLV;
740 const auto& extras = BDSAcceleratorModel::Instance()->
ExtraBeamlines();
741 for (
auto const& bl : extras)
750 G4VPhysicalVolume* containerPV,
751 G4bool checkOverlaps,
754 G4bool useCLPlacementTransform,
755 G4bool useIncrementalCopyNumbers)
761 for (
auto element : *beamline)
764 if (
dynamic_cast<BDSGap*
>(element->GetAcceleratorComponent()))
769 auto accComp = element->GetAcceleratorComponent();
770 const G4String regionName = accComp->GetRegion();
771 if(!regionName.empty())
773 G4Region* region = BDSAcceleratorModel::Instance()->
Region(regionName);
774 auto contLV = accComp->GetContainerLogicalVolume();
775 contLV->SetRegion(region);
776 region->AddRootLogicalVolume(contLV);
781 element->GetAcceleratorComponent()->AttachSensitiveDetectors();
784 G4int copyNumber = useIncrementalCopyNumbers ? i : element->GetCopyNo();
785 G4String placementName = element->GetPlacementName() +
"_pv";
786 std::set<G4VPhysicalVolume*> pvs = element->PlaceElement(placementName, containerPV, useCLPlacementTransform,
793 element->GetSPositionMiddle(),
807 const G4String& objectTypeForErrorMsg)
809 G4Transform3D result;
817 G4RotationMatrix rm = G4RotationMatrix();
820 G4ThreeVector axis = G4ThreeVector(placement.
axisX,
823 rm = G4RotationMatrix(axis, placement.
angle*CLHEP::rad);
831 CLHEP::HepEulerAngles ang = CLHEP::HepEulerAngles(placement.
phi*CLHEP::rad,
832 placement.
theta*CLHEP::rad,
833 placement.
psi*CLHEP::rad);
834 rm = G4RotationMatrix(ang);
844 {
throw BDSException(__METHOD_NAME__,
"no valid beam line yet for " + objectTypeForErrorMsg +
" w.r.t. a beam line.");}
849 G4cerr << __METHOD_NAME__ <<
"No element named \"" << placement.
referenceElement <<
"\" (instance #"
851 << placement.
name <<
"\"" << G4endl;
852 G4cout <<
"Note, this may be because the element is a bend and split into " << G4endl;
853 G4cout <<
"multiple sections with unique names. Run the visualiser to get " << G4endl;
854 G4cout <<
"the name of the segment, or place w.r.t. the element before / after." << G4endl;
855 throw BDSException(
"Error in "+objectTypeForErrorMsg+
".");
860 G4String message = objectTypeForErrorMsg +
" \"" + placement.
name +
"\" is placed using";
861 message +=
" a referenceElement but the z offset is\n non zero. Note, s should be used to";
862 message +=
" offset the placement in this case and z will\n have no effect.";
863 BDS::Warning(message);
866 sCoordinate += placement.
s * CLHEP::m;
870 G4ThreeVector offset = G4ThreeVector();
875 placement.
x * CLHEP::m + offset.x(),
876 placement.
y * CLHEP::m + offset.y());
877 G4Transform3D localRotation(rm, G4ThreeVector());
878 result = beamlinePart * localRotation;
883 {
throw BDSException(__METHOD_NAME__,
"no valid beam line yet placement w.r.t. a beam line.");}
884 G4ThreeVector offset = G4ThreeVector();
889 placement.
x * CLHEP::m + offset.x(),
890 placement.
y * CLHEP::m + offset.y());
892 G4Transform3D localRotation(rm, G4ThreeVector());
893 result = beamlinePart * localRotation;
895 {*S = placement.
s*CLHEP::m;}
899 G4ThreeVector translation = G4ThreeVector(placement.
x*CLHEP::m,
900 placement.
y*CLHEP::m,
901 placement.
z*CLHEP::m);
904 result = G4Transform3D(rm, translation);
952 G4ThreeVector result = G4ThreeVector();
953 G4String side = G4String(placement.
side);
958 G4double pathLength = placement.
s*CLHEP::m;
959 std::pair<G4double, G4double> extentZ = placementExtent.
ExtentZ();
960 G4double sLow = pathLength + extentZ.first;
961 G4double sHigh = pathLength + extentZ.second;
965 if (end != beamLine->
end())
971 for (
auto iter = start; iter != end; ++iter)
982 result.setY(sectionMaxExtent.
YPos() + placementExtent.
YPos() + ls);
983 G4double xOffset = sectionMaxExtent.
XPos() - 0.5*sectionMaxExtent.
DX();
984 result.setX(xOffset);
986 else if (side ==
"bottom")
988 result.setY(sectionMaxExtent.
YNeg() + placementExtent.
YNeg() - ls);
989 G4double xOffset = sectionMaxExtent.
XPos() - 0.5*sectionMaxExtent.
DX();
990 result.setX(xOffset);
992 else if (side ==
"left")
993 {result.setX(sectionMaxExtent.
XPos() + placementExtent.
XPos() + ls);}
994 else if (side ==
"right")
995 {result.setX(sectionMaxExtent.
XNeg() + placementExtent.
XNeg() - ls);}
997 {
throw BDSException(std::string(
"Unknown side in placement: " + side));}
1004 if (sp.apertureModel.empty())
1009 sp.aper1 * CLHEP::m,
1010 sp.aper2 * CLHEP::m,
1011 sp.aper3 * CLHEP::m,
1012 sp.aper4 * CLHEP::m,
1014 apertureExtent = aperture.
Extent();
1017 {apertureExtent =
BDSExtent(sp.aper1*CLHEP::m, sp.aper1*CLHEP::m, sp.aper2*CLHEP::m);}
1019 {apertureExtent =
BDSExtent(sp.aper1*CLHEP::m, sp.aper1*CLHEP::m, sp.aper1*CLHEP::m);}
1026 apertureExtent = aperture->
Extent();
1031 apertureExtent.
YNeg(), apertureExtent.
YPos(),
1032 1*CLHEP::um, 1*CLHEP::um);
1040 for (
const auto& samplerPlacement : samplerPlacements)
1053 return recipe.Extent();
1060 for (
const auto& mesh : scorerMeshes)
1070#if G4VERSION_NUMBER > 1009
1073 const std::list<std::string>& defaultBias,
1074 const G4String& elementName)
1077 if (biasList.empty() && defaultBias.empty())
1080 std::list<std::string> biasesAll = biasList.empty() ? defaultBias : biasList;
1083 std::set<std::string> biasNamesSorted = {biasesAll.begin(), biasesAll.end()};
1084 G4String biasSetKey;
1085 G4String biasSetPrintOut;
1086 for (
const auto& n : biasNamesSorted)
1088 biasSetKey += n +
"_";
1089 biasSetPrintOut += n +
" ";
1091 biasSetKey =
BDS::StrStrip(biasSetKey,
'_', BDS::StringStripType::trailing);
1093 auto exists = biasSetObjects.find(biasSetKey);
1094 if (exists != biasSetObjects.end())
1095 {
return exists->second;}
1098 G4cout <<
"Bias> Creating unique set of bias objects ( " << biasSetPrintOut <<
")" << G4endl;
1102 for (std::string
const & bs : biasesAll)
1104 auto result = biasObjectList.
find(bs);
1105 if (result == biasObjectList.end())
1106 {
throw BDSException(
"Error: bias named \"" + bs +
"\" not found for element named \"" + elementName +
"\"");}
1110 {G4cout << __METHOD_NAME__ <<
"bias loop : " << bs <<
" " << pb.
particle <<
" " << pb.
process << G4endl;}
1115 if (pb.
flag.size() != pb.processList.size())
1116 {
throw BDSException(__METHOD_NAME__,
"number of flag entries in \"" + pb.
name +
"\" doesn't match number of processes");}
1117 if (pb.
factor.size() != pb.processList.size())
1118 {
throw BDSException(__METHOD_NAME__,
"number of factor entries in \"" + pb.
name +
"\" doesn't match number of processes");}
1119 for (
unsigned int p = 0; p < pb.processList.size(); ++p)
1124 biasSetObjects[biasSetKey] = eg;
1131#if G4VERSION_NUMBER > 1009
1134 {G4cout << __METHOD_NAME__ <<
"registry=" << registry << G4endl;}
1136#if G4VERSION_NUMBER > 1039
1139 std::set<G4LogicalVolume*>* crystals = BDSAcceleratorModel::Instance()->
VolumeSet(
"crystals");
1140 if (!crystals->empty())
1142 G4cout << __METHOD_NAME__ <<
"Using crystal biasing: true" << G4endl;
1143 auto crystalBiasing =
new G4ChannelingOptrMultiParticleChangeCrossSection();
1144 for (
auto crystal : *crystals)
1145 {crystalBiasing->AttachTo(crystal);}
1150 for (
auto blm : blmSet)
1152 G4String biasNamesS = blm->Bias();
1153 if (biasNamesS.empty())
1156 std::list<std::string> biasNames = {biasNamesV.begin(), biasNamesV.end()};
1157 std::list<std::string> emptyDefaultBias;
1159 for (
auto lv : blm->GetAllLogicalVolumes())
1160 {biasForBLM->AttachTo(lv);}
1161 biasForBLM->AttachTo(blm->GetContainerLogicalVolume());
1167 auto defaultBiasVacuumList = std::list<std::string>(defaultBiasVacuumVector.begin(), defaultBiasVacuumVector.end());
1170 auto defaultBiasMaterialList = std::list<std::string>(defaultBiasMaterialVector.begin(), defaultBiasMaterialVector.end());
1171 G4String biasForWorldVolume = g->BiasForWorldVolume();
1173 auto biasForWorldVolumeList = std::list<std::string>(biasForWorldVolumeVector.begin(), biasForWorldVolumeVector.end());
1174 G4String biasForWorldContents = g->BiasForWorldContents();
1176 auto biasForWorldContentsList = std::list<std::string>(biasForWorldContentsVector.begin(), biasForWorldContentsVector.end());
1177 G4String biasForWorldVacuum = g->BiasForWorldVacuum();
1179 auto biasForWorldVacuumList = std::list<std::string>(biasForWorldVacuumVector.begin(), biasForWorldVacuumVector.end());
1181 G4bool useDefaultBiasVacuum = !defaultBiasVacuum.empty();
1182 G4bool useDefaultBiasMaterial = !defaultBiasMaterial.empty();
1184 G4bool useBiasForWorldVolume = !biasForWorldVolume.
empty();
1185 G4bool useBiasForWorldContents = !biasForWorldContents.empty();
1186 G4bool useBiasForWorldVacuum = !biasForWorldVacuum.empty();
1187 G4bool biasesDefined = !biasObjectList.empty();
1189 G4bool overallUseBiasing = useDefaultBiasVacuum || useDefaultBiasMaterial || biasesDefined || useBiasForWorldVolume || useBiasForWorldContents;
1190 G4cout << __METHOD_NAME__ <<
"Using generic biasing: " << BDS::BoolToString(overallUseBiasing) << G4endl;
1191 if (!overallUseBiasing)
1196 for (
auto const & item : allAcceleratorComponents)
1199 {G4cout << __METHOD_NAME__ <<
"checking component named: " << item.first << G4endl;}
1204 G4String accName = accCom->
GetName();
1209 if (!vacuumBiasList.empty() && vacuumLVs.empty())
1210 {BDS::Warning(
"biasVacuum set for component \"" + accName +
"\" but there's no 'vacuum' volume for it and it can't be biased.\nRemove biasVacuum or name it with the namedVacuumVolumes parameter.");}
1211 if ((!vacuumBiasList.empty() || useDefaultBiasVacuum) && !vacuumLVs.empty())
1214 for (
auto lv : vacuumLVs)
1217 {G4cout << __METHOD_NAME__ <<
"vacuum volume name: " << lv <<
" " << lv->GetName() << G4endl;}
1218 egVacuum->AttachTo(lv);
1224 if (!materialBiasList.empty() || useDefaultBiasMaterial)
1229 {G4cout << __METHOD_NAME__ <<
"# of logical volumes for biasing under 'material': " << allLVs.size() << G4endl;}
1230 for (
auto lv : allLVs)
1233 {G4cout << __METHOD_NAME__ <<
"Biasing 'material' logical volume: " << lv <<
" " << lv->GetName() << G4endl;}
1234 egMaterial->AttachTo(lv);
1239 if (useBiasForWorldContents)
1241 std::list<std::string> emptyList;
1244 {egWC->AttachTo(lv);}
1246 if (useBiasForWorldVolume)
1248 std::list<std::string> emptyList;
1250 egWV->AttachTo(worldLogicalVolume);
1252 if (useBiasForWorldVacuum)
1254 std::list<std::string> emptyList;
1256 for (
auto lv : worldVacuumLogicalVolumes)
1257 {egWVac->AttachTo(lv);}
1273 while ((*element).type == GMAD::ElementType::_LINE)
1276 if ((*element).type == GMAD::ElementType::_RBEND)
1292 if (scoringMeshes.empty())
1295 G4ScoringManager* scManager = G4ScoringManager::GetScoringManager();
1296 scManager->SetVerboseLevel(1);
1299 std::map<G4String, BDSScorerInfo> scorerRecipes;
1300 for (
const auto& scorer : scorers)
1303 scorerRecipes.insert(std::make_pair(si.
name, si));
1308 for (
const auto& mesh : scoringMeshes)
1314 G4String meshName = meshRecipe.name;
1324 G4String geometryType =
BDS::LowerCase(G4String(mesh.geometryType));
1326 if (geometryType ==
"box")
1329 mapper = scorerBox->Mapper();
1331 else if (geometryType ==
"cylindrical")
1334 mapper = scorerCylindrical->Mapper();
1338 G4String msg =
"mesh geometry type \"" + geometryType +
"\" is not correct. The possible options are \"box\" and \"cylindrical\"";
1343 std::vector<G4String> meshPrimitiveScorerNames;
1344 std::vector<G4double> meshPrimitiveScorerUnits;
1345 std::vector<G4String> scorerNames;
1348 for (
const auto& word : words)
1350 auto search = scorerRecipes.find(word);
1351 if (search == scorerRecipes.end())
1352 {
throw BDSException(__METHOD_NAME__,
"scorerQuantity \"" + word +
"\" for mesh \"" + meshName +
"\" not found.");}
1354 G4double psUnit = 1.0;
1355 G4VPrimitiveScorer* ps = scorerFactory.
CreateScorer(&(search->second), mapper, &psUnit, worldLV);
1360 G4String uniqueName = meshName +
"/" + ps->GetName();
1361 meshPrimitiveScorerNames.push_back(uniqueName);
1362 meshPrimitiveScorerUnits.push_back(psUnit);
1365 if (geometryType ==
"box")
1366 {scorerBox->SetPrimitiveScorer(ps);}
1367 else if (geometryType ==
"cylindrical")
1368 {scorerCylindrical->SetPrimitiveScorer(ps);}
1375 if (geometryType ==
"box")
1376 {scManager->RegisterScoringMesh(scorerBox);}
1377 else if (geometryType ==
"cylindrical")
1378 {scManager->RegisterScoringMesh(scorerCylindrical);}
1389 std::vector<BDSFieldQueryInfo*> result;
1391 for (
const auto& def : parserQueries)
1393 if (!def.queryMagneticField && !def.queryElectricField)
1394 {
throw BDSException(__METHOD_NAME__,
"neither \"queryMagneticField\" nor \"queryElectricField\" are true (=1) - one must be turned on.");}
1396 if (!def.pointsFile.empty())
1398 std::vector<G4String> columnNames;
1401 G4String(def.outfileMagnetic),
1402 G4String(def.outfileElectric),
1403 G4bool(def.queryMagneticField),
1404 G4bool(def.queryElectricField),
1407 G4bool(def.overwriteExistingFiles),
1408 G4String(def.fieldObject),
1409 def.checkParameters,
1411 def.drawZeroValuePoints,
1418 auto rot = globalTransform3D.getRotation();
1419 rot = rot.inverse();
1420 G4AffineTransform globalTransform(rot, globalTransform3D.getTranslation());
1423 G4String(def.outfileMagnetic),
1424 G4String(def.outfileElectric),
1425 G4bool(def.queryMagneticField),
1426 G4bool(def.queryElectricField),
1427 {def.nx, def.xmin*CLHEP::m, def.xmax*CLHEP::m},
1428 {def.ny, def.ymin*CLHEP::m, def.ymax*CLHEP::m},
1429 {def.nz, def.zmin*CLHEP::m, def.zmax*CLHEP::m},
1430 {def.nt, def.tmin*CLHEP::ns, def.tmax*CLHEP::ns},
1432 G4bool(def.overwriteExistingFiles),
1433 G4String(def.fieldObject),
1434 G4bool(def.printTransform),
1435 def.checkParameters,
1437 def.drawZeroValuePoints,
1445#ifdef BDSCHECKUSERLIMITS
1446void BDSDetectorConstruction::PrintUserLimitsSummary(
const G4VPhysicalVolume* world)
const
1448 G4cout <<
"USERLIMITS START" << G4endl;
1450 PrintUserLimitsPV(world, globalMinEK);
1451 G4cout <<
"USERLIMITS END" << G4endl;
1454void BDSDetectorConstruction::PrintUserLimitsPV(
const G4VPhysicalVolume* aPV, G4double globalMinEK)
const
1456 const auto lv = aPV->GetLogicalVolume();
1457 const auto ul = lv->GetUserLimits();
1461 G4double ekUL = ul->GetUserMinEkine(dummyTrack);
1463 if (ekUL < globalMinEK)
1464 {G4cout << lv->GetName() <<
" Ek Min: " << ekUL <<
" MeV < global: " << globalMinEK <<
" MeV" << G4endl;}
1467 {G4cout << lv->GetName() <<
" no G4UserLimits" << G4endl;}
1468 for (G4int i = 0; i < (G4int)lv->GetNoDaughters(); i++)
1469 {PrintUserLimitsPV(lv->GetDaughter(i), globalMinEK);}
A registry of constructed BDSAcceleratorComponent instances that can be searched.
static BDSAcceleratorComponentRegistry * Instance()
Singleton accessor.
std::map< G4String, BDSAcceleratorComponent * > AllComponentsIncludingUnique() const
void PrintNumberOfEachType() const
Print out the number of each type of component registered.
size_t size() const
Size of registry.
Abstract class that represents a component of an accelerator.
virtual void Initialise()
virtual G4String GetName() const
The name of the component without modification.
virtual std::set< G4LogicalVolume * > GetAcceleratorVacuumLogicalVolumes() const
Access the 'vacuum' volume(s) in this component if any. Default is empty set.
virtual std::set< G4LogicalVolume * > GetAcceleratorMaterialLogicalVolumes() const
Return a set of logical volumes excluding the ones in the 'vacuum' set.
std::list< std::string > GetBiasVacuumList() const
Access the bias list copied from parser.
G4String GetType() const
Get a string describing the type of the component.
std::list< std::string > GetBiasMaterialList() const
Access the bias list copied from parser.
void RegisterFields(std::vector< BDSFieldObjects * > &fieldsIn)
Register all field objects.
void RegisterBeamlineSetExtra(const G4String &name, const BDSBeamlineSet &setIn)
Register a set of beam lines to be managed and cleared up at the end of the simulation.
void RegisterWorldLV(G4LogicalVolume *worldIn)
Register constituent of world.
std::set< G4LogicalVolume * > * VolumeSet(const G4String &name)
Returns pointer to a set of logical volumes. If no set by that name exits, create it.
BDSApertureInfo * Aperture(G4String name) const
const std::map< G4String, BDSBeamlineSet > & ExtraBeamlines() const
Accessor.
BDSBeamline * TunnelBeamline() const
Access the beam line containing all the tunnel segments.
void RegisterApertures(const std::map< G4String, BDSApertureInfo * > &aperturesIn)
Register a map of apertures with associated names.
void RegisterScorerPlacement(const G4String &meshName, const G4Transform3D &placement)
Register a copy of the scorer placement so it can be used in the output.
const BDSBeamlineSet & BeamlineSetMain() const
Accessor.
const BDSBeamline * BeamlineMain() const
Accessor.
void RegisterTunnelBeamline(BDSBeamline *beamlineIn)
Register the beam line containing all the tunnel segments.
BDSBeamline * PlacementBeamline() const
Access the beam line of arbitrary placements.
void RegisterWorldPV(G4VPhysicalVolume *worldIn)
Register constituent of world.
void RegisterPlacementBeamline(BDSBeamline *placementBLIn)
Register the beam line of arbitrary placements.
void RegisterBeamlineSetMain(const BDSBeamlineSet &setIn)
Register the main beam line set.
void RegisterScorerHistogramDefinition(const BDSScorerHistogramDef &def)
G4LogicalVolume * WorldLV() const
Access the logical volume of the world.
void RegisterBLMs(BDSBeamline *blmsIn)
Register a beam line of blm placements.
G4Region * Region(const G4String &name) const
Access region information. Will exit if not found.
void RegisterWorldSolid(G4VSolid *worldIn)
Register constituent of world.
void RegisterRegion(BDSRegion *region)
Register a region.
Holder class for all information required to describe an aperture.
static void AttachWorldVolumeToNavigator(G4VPhysicalVolume *worldPVIn)
Setup the navigator w.r.t. to a world volume - typically real world.
static BDSBLMRegistry * Instance()
Accessor for registry.
Multi-particle cross-section changer.
A class that holds a fully constructed BDSAcceleratorComponent as well as any information relevant to...
G4double GetSPositionEnd() const
Accessor.
G4double GetSPositionMiddle() const
Accessor.
Simple struct to return a beamline plus associated beam lines.
void GetExtentGlobals(std::vector< BDSExtentGlobal > &extents) const
Append global extents of all beamlines to supplied vector.
A vector of BDSBeamlineElement instances - a beamline.
const BDSBeamlineElement * GetElement(G4String acceleratorComponentName, G4int i=0) const
Get the ith placement of an element in the beam line. Returns null pointer if not found.
BDSBeamlineElement * back() const
Return a reference to the last element.
G4bool ElementAnglesSumToCircle() const
Check if the sum of the angle of all elements is two pi.
reverse_iterator rend()
Iterator mechanics.
const_iterator FindFromS(G4double S) const
Returns an iterator to the beamline element at s.
G4bool empty() const
Iterator mechanics.
G4Transform3D GetGlobalEuclideanTransform(G4double s, G4double x=0, G4double y=0, G4int *indexOfFoundElement=nullptr) const
BDSExtentGlobal GetExtentGlobal() const
Get the global extents for this beamline.
reverse_iterator rbegin()
Iterator mechanics.
iterator begin()
Iterator mechanics.
void AddComponent(BDSAcceleratorComponent *component, BDSTiltOffset *tiltOffset=nullptr, BDSSamplerInfo *samplerInfo=nullptr)
iterator end()
Iterator mechanics.
G4double GetTotalArcLength() const
Get the total ARC length for the beamline - ie the maximum s position.
Factory for user specified accelerator components.
void SetDesignParticle(const BDSParticleDefinition *designParticleIn)
Update values for nominal rigidity and relativisitic beta of the beam particle.
Factory to produce all types of BDSAcceleratorComponents.
BDSAcceleratorComponent * CreateTeleporter(const G4double teleporterLength, const G4double teleporterWidth, const G4Transform3D &transformIn)
BDSAcceleratorComponent * CreateTerminator(G4double witdth)
BDSAcceleratorComponent * CreateComponent(GMAD::Element const *elementIn, GMAD::Element const *prevElementIn, GMAD::Element const *nextElementIn, G4double currentArcLength=0)
static BDSTiltOffset * CreateTiltOffset(GMAD::Element const *el)
Create the tilt and offset information object by inspecting the parser element.
Factory for simple parallel geometry for curvilinear coordinates.
BDSBeamline * BuildCurvilinearBeamLine1To1(BDSBeamline const *const beamline, const G4bool circular)
Build a beam line of curvilinear geometry based on another beam line.
BDSBeamline * BuildCurvilinearBridgeBeamLine(BDSBeamline const *const beamline)
Build bridging volumes to join the curvilinear ones.
G4bool circular
Whether or not we're building a circular machine.
BDSBeamlineSet BuildBeamline(const GMAD::FastList< GMAD::Element > &beamLine, const G4String &name, const G4Transform3D &initialTransform=G4Transform3D(), G4double initialS=0.0, G4bool beamlineIsCircular=false, G4bool isPlacementBeamline=false)
void UpdateSamplerDiameterAndCountSamplers()
BDSBOptrMultiParticleChangeCrossSection * BuildCrossSectionBias(const std::list< std::string > &biasList, const std::list< std::string > &defaultBias, const G4String &elementName)
Function that creates physics biasing cross section.
void InitialiseRegions()
Create and set parameters for various G4Regions.
std::vector< BDSBOptrMultiParticleChangeCrossSection * > biasObjects
List of bias objects - for memory management.
BDSExtent CalculateExtentOfScorerMesh(const GMAD::ScorerMesh &sm) const
Calculate local extent of scorer mesh in 3D.
static void PlaceBeamlineInWorld(BDSBeamline *beamline, G4VPhysicalVolume *containerPV, G4bool checkOverlaps=false, G4bool setRegions=false, G4bool registerInfo=false, G4bool useCLPlacementTransform=false, G4bool useIncrementalCopyNumbers=false)
BDSAcceleratorModel * acceleratorModel
Accelerator model pointer.
static G4ThreeVector SideToLocalOffset(const GMAD::Placement &placement, const BDSBeamline *beamLine, const BDSExtent &placementExtent)
Attach component with extent2 to component with extent1 with placement.
G4bool verbose
Variable copied from global constants.
BDSBeamline * placementBL
static std::vector< BDSFieldQueryInfo * > PrepareFieldQueries(const BDSBeamline *mainBeamline)
Prepare field queries from parser information.
G4VPhysicalVolume * BuildWorld()
void InitialiseApertures()
Create all aperture definitions from parser and store in BDSAcceleratorModel.
void CountPlacementFields()
Count number of fields required for placements.
G4bool canSampleAngledFaces
Whether the integrator set permits sampling elements with angled faces.
BDSExtent worldExtent
Record of the world extent.
std::set< G4LogicalVolume * > worldContentsLogicalVolumes
Cache of possibly loaded logical volumes from a world geometry file - used for biasing.
static BDSSamplerInfo * BuildSamplerInfo(const GMAD::Element *element)
BDSExtent CalculateExtentOfSamplerPlacement(const GMAD::SamplerPlacement &sp) const
Calculate local extent of custom user sampler.
void ConstructScoringMeshes()
Construct scoring meshes.
BDSExtentGlobal CalculateExtentOfScorerMeshes(const BDSBeamline *bl) const
G4bool checkOverlaps
Variable copied from global constants.
G4int nSamplers
Count of number of samplers to be built.
void ComponentPlacement(G4VPhysicalVolume *worldPV)
Place beam line, tunnel beam line, end pieces and placements in world.
void BuildBeamlines()
Build the main beam line and then any other required beam lines.
const BDSParticleDefinition * designParticle
Particle definition all components are built w.r.t. Includes rigidity etc.
void BuildTunnel()
Build the tunnel around the already constructed flat beam line.
virtual void ConstructSDandField()
Construct sensitive detectors and fields.
void PrepareExtraSamplerSDs()
virtual G4VPhysicalVolume * Construct()
static G4Transform3D CreatePlacementTransform(const GMAD::Placement &placement, const BDSBeamline *beamLine, G4double *S=nullptr, BDSExtent *placementExtent=nullptr, const G4String &objectTypeForErrorMsg="placement")
G4bool UnsuitableFirstElement(std::list< GMAD::Element >::const_iterator element)
BDSExtentGlobal CalculateExtentOfSamplerPlacements(const BDSBeamline *beamLine) const
General exception with possible name of object and message.
Holder for +- extents in 3 dimensions with a rotation and translation.
BDSExtentGlobal ExpandToEncompass(const BDSExtentGlobal &other) const
Return a copy of this extent but expanded to encompass another global extent.
G4bool Encompasses(const BDSExtentGlobal &otherExtent)
Return whether the extent encompasses another extent.
Holder for +- extents in 3 dimensions.
G4double XPos() const
Accessor.
G4double XNeg() const
Accessor.
G4double YNeg() const
Accessor.
std::pair< G4double, G4double > ExtentZ() const
Accessor.
G4double DX() const
The difference in a dimension.
G4double YPos() const
Accessor.
static BDSFieldBuilder * Instance()
Singleton pattern accessor.
Holder class for all information required for a field query.
static void AttachWorldVolumeToNavigator(G4VPhysicalVolume *worldPVIn)
Setup the navigator w.r.t. to a world volume.
Gap in accelerator beamline.
G4LogicalVolume * GetContainerLogicalVolume() const
Accessor - see member for more info.
BDSExtent GetExtent() const
Accessor - see member for more info.
virtual void AttachSensitiveDetectors()
Attach a sensitive detector class to all registered sensitive volumes in this component.
virtual std::set< G4LogicalVolume * > GetAllLogicalVolumes() const
Access all logical volumes belonging to this component.
G4VSolid * GetContainerSolid() const
Accessor - see member for more info.
A loaded piece of externally provided geometry.
const std::set< G4LogicalVolume * > & VacuumVolumes() const
Access the vacuum volumes.
static BDSGeometryFactory * Instance()
Singleton accessor.
BDSGeometryExternal * BuildGeometry(G4String componentName, const G4String &formatAndFilePath, std::map< G4String, G4Colour * > *colourMapping=nullptr, G4bool autoColour=true, G4double suggestedLength=0, G4double suggestedHorizontalWidth=0, std::vector< G4String > *namedVacuumVolumes=nullptr, G4bool makeSensitive=true, BDSSDType sensitivityType=BDSSDType::energydep, G4bool stripOuterVolumeAndMakeAssembly=false, G4UserLimits *userLimitsToAttachToAllLVs=nullptr, G4bool dontReloadGeometry=false)
A class that holds global options and constants.
void SetCurvilinearDiameterShrunkForBends()
Setter.
static BDSGlobalConstants * Instance()
Access method.
void SetCurvilinearDiameter(G4double curvilinearDiameterIn)
Setter.
void SetSamplerDiameter(G4double samplerDiameterIn)
Setter.
Mapping from axis indices to 1D index.
A class that hold multiple accelerator components.
static BDSMaterials * Instance()
Singleton pattern access.
G4Material * GetMaterial(G4String material) const
Get material by name.
const GMAD::FastList< GMAD::Element > & GetSequence(const std::string &name)
Return sequence.
std::vector< GMAD::Scorer > GetScorers() const
Return scorer list.
static BDSParser * Instance()
Access method.
const GMAD::Options & GetOptions() const
Return options.
std::vector< GMAD::Query > GetQuery() const
Query list.
std::vector< GMAD::SamplerPlacement > GetSamplerPlacements() const
Return sampler placement list.
const GMAD::FastList< GMAD::PhysicsBiasing > & GetBiasing() const
Return biasing list.
std::vector< GMAD::ScorerMesh > GetScorerMesh() const
Return scorermesh list.
std::vector< GMAD::Placement > GetPlacements() const
Return the vector of placement objects.
void RegisterExcludedPV(G4VPhysicalVolume *physicalVolume)
void RegisterInfo(G4VPhysicalVolume *physicalVolume, BDSPhysicalVolumeInfo *info, G4bool isReadOutVolume=false, G4bool isTunnel=false)
static BDSPhysicalVolumeInfoRegistry * Instance()
Singleton accessor.
A class holding any information pertaining to a particular physical volume in a BDSIM geant4 model.
Range cuts for a region. Help with defaults.
void ConstructSamplerSDsForParticleSets(const std::map< int, std::set< int > > &samplerFilterIDtoPDGSetIn)
void RegisterPrimitiveScorerNames(const std::vector< G4String > &namesIn, const std::vector< G4double > *units=nullptr)
Loop over a vector and apply above single name function.
All info required to build a sampler but not place it.
Create primitive scorers on demand.
G4VPrimitiveScorer * CreateScorer(const BDSScorerInfo *info, const BDSHistBinMapper *mapper, G4double *unit=nullptr, G4LogicalVolume *worldLV=nullptr)
Main function to create a scorer.
Definition for a scorer histogram.
Recipe class for scorer. Checks values.
G4String name
Scorer name.
Recipe class for a scoring mesh.
Wrapper for G4ScoringBox to allow full access to placement.
Wrapper for G4ScoringCylinder to allow full access to placement.
A class of functions to output Geant4/Mokka/BDSIM parameters for the beamline.
void Write(BDSBeamlineElement *beamlineElement)
write line
A holder for any placement offsets and rotations for a BDSAcceleratorComponent.
static G4double dEThresholdForScattering
A class that constructs tunnel segments around a beamline.
BDSBeamline * BuildTunnelSections(const BDSBeamline *flatBeamLine) const
bool empty() const
Whether the list is empty.
FastListConstIterator find(std::string name, unsigned int count=1) const
int size() const
size of list
const FastList< Element > & GetBeamline() const
Physics biasing class for parser.
std::vector< double > factor
factors corresponding to process
std::vector< PhysicsBiasingType > flag
flag which particles are biased
std::string particle
particle name
std::string process
geant4 process: single string, but can have multiple processes separated with a space
Placement class for parser.
bool axisAngle
Flag to use the axis angle construction of rotation.
std::string name
Name of this placement.
std::string side
which side to attach to: top, bottom, left, right.
double theta
Euler angle for rotation.
double sideOffset
Gap between side and component.
double psi
Euler angle for rotation.
double s
Curvilinear s position to place w.r.t..
int referenceElementNumber
Index of repetition of element if there are multiple uses.
double phi
Euler angle for rotation.
std::string referenceElement
Name of reference element w.r.t. to place to.
Query structure class for parser.
Sampler placement class for parser.
std::string name
Name of this samplerplacement.
std::string samplerType
Plane, Cylinder, Sphere.
ScorerMesh class for parser.
BDSBeamline * BuildEndPieceBeamline(const BDSBeamline *beamline, const G4bool circularMachine)
G4Transform3D CalculateTeleporterDelta(const BDSBeamline *beamline, G4double &teleporterLength)
std::vector< BDSFourVector< G4double > > LoadFieldQueryPoints(const G4String &fileName, std::vector< G4String > *columnNamesIn)
BDSExtent MaximumCombinedExtent(const BDSExtent &first, const BDSExtent &second)
Returns the extent which is the greatest extent in all six directions.
BDSBeamline * BuildPlacementGeometry(const std::vector< GMAD::Placement > &placements, const BDSBeamline *parentBeamLine, BDSComponentFactory *componentFactory)
G4String LowerCase(const G4String &str)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
std::vector< G4String > SplitOnWhiteSpace(const G4String &input)
Split a string on whitespace and return a vector of these 'words'.
BDSBeamline * BuildBLMs(const std::vector< GMAD::BLMPlacement > &blmPlacements, const BDSBeamline *parentBeamLine)
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())
G4String StrStrip(const G4String &str, char ch, StringStripType stripType=StringStripType::both)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
std::string samplerType
element has a sampler of this type (default "none")
double e1
input pole face rotation for bends
std::string samplerName
name of sampler (default empty)