BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSLinkDetectorConstruction.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSAcceleratorModel.hh"
20#include "BDSBeamline.hh"
21#include "BDSBeamlineElement.hh"
22#include "BDSCollimatorJaw.hh"
23#include "BDSComponentFactory.hh"
24#include "BDSCrystalInfo.hh"
25#include "BDSDebug.hh"
26#include "BDSException.hh"
27#include "BDSExtent.hh"
28#include "BDSExtentGlobal.hh"
29#include "BDSGlobalConstants.hh"
30#include "BDSLinkComponent.hh"
31#include "BDSLinkDetectorConstruction.hh"
32#include "BDSLinkOpaqueBox.hh"
33#include "BDSLinkPrimaryGeneratorAction.hh"
34#include "BDSLinkRegistry.hh"
35#include "BDSMaterials.hh"
36#include "BDSParallelWorldSampler.hh"
37#include "BDSParser.hh"
38#include "BDSSampler.hh"
39#include "BDSSamplerInfo.hh"
40#include "BDSSamplerPlane.hh"
41#include "BDSSamplerRegistry.hh"
42#include "BDSSamplerType.hh"
43#include "BDSSDManager.hh"
44#include "BDSTiltOffset.hh"
45
46#include "parser/element.h"
47#include "parser/elementtype.h"
48
49#include "G4Box.hh"
50#include "G4PVPlacement.hh"
51#include "G4String.hh"
52#include "G4ThreeVector.hh"
53#include "G4Types.hh"
54#include "G4Version.hh"
55#include "G4VisAttributes.hh"
56#if G4VERSION_NUMBER > 1039
57#include "G4ChannelingOptrMultiParticleChangeCrossSection.hh"
58#endif
59
60#include <set>
61#include <vector>
62
64
66 worldSolid(nullptr),
67 worldPV(nullptr),
68 linkBeamline(nullptr),
69 linkRegistry(nullptr),
70 primaryGeneratorAction(nullptr),
71 designParticle(nullptr),
72#if G4VERSION_NUMBER > 1039
73 crystalBiasing(nullptr),
74#endif
75 samplerWorldID(-1)
76{
77 linkRegistry = new BDSLinkRegistry();
78 BDSSDManager::Instance()->SetLinkRegistry(linkRegistry);
79}
80
81BDSLinkDetectorConstruction::~BDSLinkDetectorConstruction()
82{
83 delete linkBeamline;
84 delete linkRegistry;
85#if G4VERSION_NUMBER > 1039
86 delete crystalBiasing;
87#endif
88}
89
90G4VPhysicalVolume* BDSLinkDetectorConstruction::Construct()
91{
93
94 auto componentFactory = new BDSComponentFactory(designParticle, nullptr, false);
95 auto beamline = BDSParser::Instance()->GetBeamline();
96
97 std::vector<BDSLinkOpaqueBox*> opaqueBoxes = {};
98 linkBeamline = new BDSBeamline();
99
100 auto acceleratorModel = BDSAcceleratorModel::Instance();
101
102 for (const auto& element : beamline)
103 {
104 GMAD::ElementType eType = element.type;
105
106 if (eType == GMAD::ElementType::_LINE || eType == GMAD::ElementType::_REV_LINE)
107 {continue;}
108
109 std::set<GMAD::ElementType> acceptedTypes = {GMAD::ElementType::_ECOL,
110 GMAD::ElementType::_RCOL,
111 GMAD::ElementType::_JCOL,
112 GMAD::ElementType::_CRYSTALCOL,
113 GMAD::ElementType::_ELEMENT};
114 auto search = acceptedTypes.find(eType);
115 if (search == acceptedTypes.end())
116 {throw BDSException(G4String("Unsupported element type for link = " + GMAD::typestr(eType)));}
117
118 // Only need first argument, the rest pertain to beamlines.
119 BDSAcceleratorComponent* component = componentFactory->CreateComponent(&element,
120 nullptr,
121 nullptr,
122 0);
123
124 BDSTiltOffset* to = new BDSTiltOffset(element.offsetX * CLHEP::m,
125 element.offsetY * CLHEP::m,
126 element.tilt * CLHEP::rad);
127 auto extentTiltOffset = component->GetExtent().TiltOffset(to);
128 G4double encompassingRadius = extentTiltOffset.TransverseBoundingRadius();
129 BDSLinkOpaqueBox* opaqueBox = new BDSLinkOpaqueBox(component, to, encompassingRadius);
130
131 delete to; // opaqueBox doesn't own it
132 opaqueBoxes.push_back(opaqueBox);
133
134 BDSLinkComponent* comp = new BDSLinkComponent(opaqueBox->GetName(),
135 opaqueBox,
136 opaqueBox->GetExtent().DZ());
137 acceleratorModel->RegisterLinkComponent(comp); // memory management
138 nameToElementIndex[element.name] = (G4int)linkBeamline->size();
139 linkBeamline->AddComponent(comp);
140 }
141
142 // update world extents and world solid
144
145 G4LogicalVolume* worldLV = new G4LogicalVolume(worldSolid,
146 BDSMaterials::Instance()->GetMaterial("G4_Galactic"),
147 "world_lv");
148
149 G4VisAttributes* debugWorldVis = new G4VisAttributes(*(BDSGlobalConstants::Instance()->ContainerVisAttr()));
150 debugWorldVis->SetForceWireframe(true);//just wireframe so we can see inside it
151 worldLV->SetVisAttributes(debugWorldVis);
152 worldLV->SetUserLimits(globalConstants->DefaultUserLimits());
153
154 worldPV = new G4PVPlacement(nullptr,
155 G4ThreeVector(),
156 worldLV,
157 "world_pv",
158 nullptr,
159 false,
160 0,
161 true);
162
163 // place any defined link elements in input
164 for (auto element : *linkBeamline)
165 {
166 BDSLinkComponent* lc = dynamic_cast<BDSLinkComponent*>(element->GetAcceleratorComponent());
167 BDSSamplerInfo* samplerInfo = element->GetSamplerInfo();
168 G4String samplerName = samplerInfo ? samplerInfo->name : "unknown";
169 G4String name = lc ? lc->LinkName() : samplerName;
170 G4int linkID = PlaceOneComponent(element, name);
171 nameToElementIndex[name] = linkID;
172 }
173
174 delete componentFactory;
175
176 return worldPV;
177}
178
179G4int BDSLinkDetectorConstruction::AddLinkCollimatorJaw(const std::string& collimatorName,
180 const std::string& materialName,
181 G4double length,
182 G4double halfApertureLeft,
183 G4double halfApertureRight,
184 G4double rotation,
185 G4double xOffset,
186 G4double yOffset,
187 G4double jawTiltLeft,
188 G4double jawTiltright,
189 G4bool buildLeftJaw,
190 G4bool buildRightJaw,
191 G4bool isACrystalIn,
192 G4double crystalAngle,
193 G4bool /*sampleIn*/)
194{
195 auto componentFactory = new BDSComponentFactory(designParticle, nullptr, false);
196
197 std::map<std::string, std::string> collimatorToCrystal =
198 {
199 {"cry.mio.b1", "stf75"}, // b1 h
200 {"cry.mio.b2", "tcp76"}, // b2 h
201 {"tcpv.a6l7.b1", "qmp34"}, // b1 v
202 {"tcpv.a6r7.b2", "qmp53"}, // b2 v
203 {"tcpch.a4l7.b1", "stf75"},// b1 h new name
204 {"tcpcv.a6l7.b1", "tcp76"} // b2 v new name
205 };
206 G4String collimatorLower = BDS::LowerCase(collimatorName);
207 auto searchC = collimatorToCrystal.find(collimatorLower); // will use later if needed
208 G4bool isACrystal = searchC != collimatorToCrystal.end();
209 if (!isACrystal && isACrystalIn)
210 {throw BDSException("BDSLinkDetectorConstruction", "no matching crystal name found but it is flagged as a crystal in input");}
211 G4cout << "XYZ isACrystal " << isACrystal << G4endl;
212 if (isACrystal)
213 {G4cout << "crystal name " << searchC->first << " " << searchC->second << G4endl;}
214
215 std::map<std::string, std::string> sixtrackToBDSIM =
216 {
217 {"CU", "Cu"},
218 {"W", "W"},
219 {"C", "G4_GRAPHITE_POROUS"},
220 {"Si", "Si"},
221 {"SI", "Si"}
222 };
223 std::string g4material;
224 auto search = sixtrackToBDSIM.find(materialName);
225 if (search != sixtrackToBDSIM.end())
226 {g4material = search->second;}
227 else
228 {g4material = materialName;}
229
230 // build component
232 el.type = GMAD::ElementType::_JCOL;
233 el.name = collimatorName;
234 el.material = g4material;
235 el.l = length / CLHEP::m;
236 el.xsizeLeft = halfApertureLeft / CLHEP::m;
237 el.xsizeRight = halfApertureRight / CLHEP::m;
238 el.ysize = 0.2; // half size
239 el.tilt = rotation / CLHEP::rad;
240 el.offsetX = xOffset / CLHEP::m;
241 el.offsetY = yOffset / CLHEP::m;
242 el.horizontalWidth = 2.0; // m
243 el.jawTiltLeft = jawTiltLeft; // rad
244 el.jawTiltRight = jawTiltright; // rad
245
246 // if we don't want to build a jaw, then we set it to outside the width.
247 if (!buildLeftJaw)
248 {el.xsizeLeft = el.horizontalWidth * 1.2;}
249 if (!buildRightJaw)
250 {el.xsizeRight = el.horizontalWidth * 1.2;}
251
252 if (isACrystal)
253 {
254 // find the bending angle of this particular crystal
255 // so we can add half of that on for the BDSIM convention of the 0 about the centre for crystals
256 G4String crystalNameC = searchC->second;
257 G4cout << "crystal name " << crystalNameC << G4endl;
258 BDSCrystalInfo* ci = componentFactory->PrepareCrystalInfo(crystalNameC);
259 crystalAngle *= CLHEP::rad;
260
261 // crucial - crystal only responds to xsize - not xsizeLeft
262 el.xsize = el.xsizeLeft;
263
264 el.type = GMAD::ElementType::_CRYSTALCOL;
265 el.apertureType = "circularvacuum";
266 el.aper1 = 0.5; // m
267 // need a small margin in length as crystal may have angled face and be rotated
268 el.l += 10e-6; // TODO - confirm margin with sixtrack interface backtracking on input side
269 if (collimatorName.find('2') != std::string::npos) // b2
270 {
271 el.crystalLeft = crystalNameC;
272 el.crystalAngleYAxisLeft = crystalAngle + 0.5 * ci->bendingAngleYAxis;
273 }
274 else
275 {
276 el.crystalRight = crystalNameC;
277 el.crystalAngleYAxisRight = crystalAngle + 0.5 * ci->bendingAngleYAxis;
278 }
279
280 G4cout << "XYZKEY Crystal angle " << crystalAngle << G4endl;
281 G4cout << "xsizeLeft " << el.xsizeLeft << G4endl;
282 G4cout << "xsizeRight " << el.xsizeRight << G4endl;
283 G4cout << "l crystal angle " << el.crystalAngleYAxisLeft << G4endl;
284 G4cout << "r crystal angle " << el.crystalAngleYAxisRight << G4endl;
285 G4cout << "rotation (tilt) " << el.tilt << G4endl;
286 delete ci; // no longer needed
287 }
288 else
289 {el.region = "r1";} // stricter range cuts for default collimators
290
291 BDSAcceleratorComponent* component = nullptr;
292 try
293 {component = componentFactory->CreateComponent(&el, nullptr, nullptr, 0);}
294 catch (const BDSException& e)
295 {
296 G4cout << e.what() << G4endl;
297 G4cout << "Replacing component " << el.name << " with drift" << G4endl;
298 // well it didn't work (maybe ridiculous unphysical gap - so replace it with a drift
299 el.type = GMAD::ElementType::_DRIFT;
300 el.apertureType = "circularvacuum";
301 component = componentFactory->CreateComponent(&el, nullptr, nullptr, 0);
302 }
303
304 // wrap in box
305 BDSTiltOffset* to = new BDSTiltOffset(el.offsetX * CLHEP::m,
306 el.offsetY * CLHEP::m,
307 el.tilt * CLHEP::rad);
308 auto extentTiltOffset = component->GetExtent().TiltOffset(to);
309 G4double encompassingRadius = extentTiltOffset.TransverseBoundingRadius();
310 BDSLinkOpaqueBox* opaqueBox = new BDSLinkOpaqueBox(component, to, encompassingRadius);
311
312 // add to beam line
313 BDSLinkComponent* comp = new BDSLinkComponent(opaqueBox->GetName(),
314 opaqueBox,
315 opaqueBox->GetExtent().DZ());
316 BDSAcceleratorModel::Instance()->RegisterLinkComponent(comp);
317 BDSSamplerInfo* samplerInfo = new BDSSamplerInfo(comp->GetName() + "_out", BDSSamplerType::plane);
318 linkBeamline->AddComponent(comp, nullptr, samplerInfo);
319
320 // update world extents and world solid
322
323 // place that one element
324 G4int linkID = PlaceOneComponent(linkBeamline->back(), collimatorName);
325 nameToElementIndex[collimatorName] = linkID;
326 linkIDToBeamlineIndex[linkID] = (G4int)linkBeamline->size() - 1;
327
328 // update crystal biasing
329 BuildPhysicsBias();
330
331 return linkID;
332}
333
335{
336 BDSExtentGlobal we = linkBeamline->GetExtentGlobal();
337 we = we.ExpandToEncompass(BDSExtentGlobal(BDSExtent(10*CLHEP::m, 10*CLHEP::m, 10*CLHEP::m))); // minimum size
338 G4ThreeVector worldExtentAbs = we.GetMaximumExtentAbsolute();
339 worldExtentAbs *= 1.2;
340
341 if (!worldSolid)
342 {
343 worldSolid = new G4Box("world_solid",
344 worldExtentAbs.x(),
345 worldExtentAbs.y(),
346 worldExtentAbs.z());
347 }
348 else
349 {
350 worldSolid->SetXHalfLength(worldExtentAbs.x());
351 worldSolid->SetYHalfLength(worldExtentAbs.y());
352 worldSolid->SetZHalfLength(worldExtentAbs.z());
353 }
354 worldExtent = BDSExtent(worldExtentAbs);
355 if (primaryGeneratorAction)
356 {primaryGeneratorAction->SetWorldExtent(worldExtent);}
357}
358
360 const G4String& originalName)
361{
362 G4String placementName = element->GetPlacementName() + "_pv";
363 G4int copyNumber = element->GetCopyNo();
364 std::set<G4VPhysicalVolume*> pvs = element->PlaceElement(placementName, worldPV, false, copyNumber, true);
365
366 auto lc = dynamic_cast<BDSLinkComponent*>(element->GetAcceleratorComponent());
367 if (!lc)
368 {return -1;}
369 BDSLinkOpaqueBox* el = lc->Component();
370 G4Transform3D* placementTransform = element->GetPlacementTransform();
371 G4Transform3D elCentreToStart = el->TransformToStart();
372 G4Transform3D globalToStart = elCentreToStart * (*placementTransform);
373 G4int linkID = linkRegistry->Register(el, globalToStart);
374
375 G4ThreeVector zOffset = G4ThreeVector(0,0,BDSGlobalConstants::Instance()->LengthSafety()+BDSSamplerPlane::ChordLength());
376 G4Transform3D samplerPosition = globalToStart * G4Transform3D(G4RotationMatrix(), globalToStart.getRotation()*zOffset);
377
378 if (element->GetSamplerType() == BDSSamplerType::plane && samplerWorldID >= 0)
379 {
380 auto samplerWorldRaw = GetParallelWorld(samplerWorldID);
381 auto samplerWorld = dynamic_cast<BDSParallelWorldSampler*>(samplerWorldRaw);
382 if (!samplerWorld)
383 {return -1;}
384
385 BDSSamplerPlane* sampler = samplerWorld->GeneralPlane();
386 G4String samplerName = originalName + "_in";
387 G4double sStart = element->GetSPositionStart();
388
389 G4int samplerID = BDSSamplerRegistry::Instance()->RegisterSampler(samplerName, sampler, samplerPosition, sStart, element);
390
391 G4LogicalVolume* samplerWorldLV = samplerWorld->WorldLV();
392 new G4PVPlacement(samplerPosition,
393 sampler->GetContainerLogicalVolume(),
394 samplerName + "_pv",
395 samplerWorldLV,
396 false,
397 samplerID,
398 false);
399 }
400 return linkID;
401}
402
403void BDSLinkDetectorConstruction::BuildPhysicsBias()
404{
405#if G4VERSION_NUMBER > 1039
406 if (!crystalBiasing) // cache it because we may have to dynamically add later
407 {crystalBiasing = new G4ChannelingOptrMultiParticleChangeCrossSection();}
408
409 // crystal biasing necessary for implementation of variable density
410 std::set<G4LogicalVolume*>* crystals = BDSAcceleratorModel::Instance()->VolumeSet("crystals");
411 if (!crystals->empty())
412 {
413 G4cout << __METHOD_NAME__ << "Using crystal biasing: true" << G4endl; // to match print out style further down
414 for (auto crystal : *crystals)
415 {
416 // if it hasn't been added already - g4 will whine for double registration
417 if (!crystalBiasing->GetBiasingOperator(crystal))
418 {crystalBiasing->AttachTo(crystal);}
419 }
420 }
421#endif
422}
Abstract class that represents a component of an accelerator.
virtual G4String GetName() const
The name of the component without modification.
std::set< G4LogicalVolume * > * VolumeSet(const G4String &name)
Returns pointer to a set of logical volumes. If no set by that name exits, create it.
A class that holds a fully constructed BDSAcceleratorComponent as well as any information relevant to...
BDSSamplerType GetSamplerType() const
Accessor.
std::set< G4VPhysicalVolume * > PlaceElement(const G4String &pvName, G4VPhysicalVolume *containerPV, G4bool useCLPlacementTransform, G4int copyNumber, G4bool checkOverlaps) const
G4String GetPlacementName() const
Accessor.
G4Transform3D * GetPlacementTransform() const
Accessor.
BDSAcceleratorComponent * GetAcceleratorComponent() const
Accessor.
G4int GetCopyNo() const
Accessor.
G4double GetSPositionStart() const
Accessor.
A vector of BDSBeamlineElement instances - a beamline.
Definition: BDSBeamline.hh:61
BeamlineVector::size_type size() const
Get the number of elements.
Definition: BDSBeamline.hh:148
void AddComponent(BDSAcceleratorComponent *component, BDSTiltOffset *tiltOffset=nullptr, BDSSamplerInfo *samplerInfo=nullptr)
Definition: BDSBeamline.cc:124
Factory to produce all types of BDSAcceleratorComponents.
Holder for all information required to create a crystal.
G4double bendingAngleYAxis
Bending angle about Y axis.
General exception with possible name of object and message.
Definition: BDSException.hh:35
const char * what() const noexcept override
Override message in std::exception.
Definition: BDSException.hh:55
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.
G4ThreeVector GetMaximumExtentAbsolute() const
Get the maximum extent absolute in each dimension.
Holder for +- extents in 3 dimensions.
Definition: BDSExtent.hh:39
G4double DZ() const
The difference in a dimension.
Definition: BDSExtent.hh:85
BDSExtent TiltOffset(const BDSTiltOffset *tiltOffset) const
Provide a new copy of this extent with both a tilt and an offset applied.
Definition: BDSExtent.cc:105
G4double TransverseBoundingRadius() const
Return a radius that would encompass the maximum x,y extent.
Definition: BDSExtent.cc:183
G4LogicalVolume * GetContainerLogicalVolume() const
Accessor - see member for more info.
BDSExtent GetExtent() const
Accessor - see member for more info.
virtual G4String GetName() const
Accessor - see member for more info.
A class that holds global options and constants.
static BDSGlobalConstants * Instance()
Access method.
A BDSAcceleratorComponent wrapper for BDSLinkOpaqueBox.
BDSLinkOpaqueBox * Component() const
Accessor.
G4String LinkName() const
Accessor.
G4int AddLinkCollimatorJaw(const std::string &collimatorName, const std::string &materialName, G4double length, G4double halfApertureLeft, G4double halfApertureRight, G4double rotation, G4double xOffset, G4double yOffset, G4double jawTiltLeft=0.0, G4double jawTiltRight=0.0, G4bool buildLeftJaw=true, G4bool buildRightJaw=true, G4bool isACrystal=false, G4double crystalAngle=0, G4bool sampleIn=false)
Interface to append a collimator of jaw style to the linking.
G4int samplerWorldID
Cache of the index to which parallel world the sampler one is.
G4int PlaceOneComponent(const BDSBeamlineElement *element, const G4String &originalName)
Place a beam line element in the world.
const BDSParticleDefinition * designParticle
Particle definition all components are built w.r.t. Includes rigidity etc.
BDSLinkDetectorConstruction()
Default constructor.
std::map< std::string, G4int > nameToElementIndex
Build up a copy here too.
std::map< G4int, G4int > linkIDToBeamlineIndex
Special linkID to linkBeamline index.
Wrapper box for an accelerator component.
void SetWorldExtent(const BDSExtent worldExtentIn)
Set the world extent that particle coordinates will be checked against.
Registry / map of components for tracker linkage.
static BDSMaterials * Instance()
Singleton pattern access.
Definition: BDSMaterials.cc:39
A parallel world for sampler planes.
static BDSParser * Instance()
Access method.
Definition: BDSParser.cc:28
Wrapper for particle definition.
void SetLinkRegistry(BDSLinkRegistry *registry)
If samplerLink member exists, set the registry to look up links for that SD.
All info required to build a sampler but not place it.
Rectangular sampler with fixed thickness but variable x,y.
static G4double ChordLength()
Access the sampler plane length in other classes.
static BDSSamplerRegistry * Instance()
Accessor for registry.
G4int RegisterSampler(const G4String &name, BDSSampler *sampler, const G4Transform3D &transform=G4Transform3D(), G4double S=-1000, const BDSBeamlineElement *element=nullptr, BDSSamplerType type=BDSSamplerType::plane, G4double radius=0)
A holder for any placement offsets and rotations for a BDSAcceleratorComponent.
const FastList< Element > & GetBeamline() const
Definition: parser.cc:903
G4String LowerCase(const G4String &str)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
ElementType
types of elements
Definition: elementtype.h:28
std::string typestr(ElementType type)
conversion from enum to string
Definition: elementtype.cc:29
Element class.
Definition: element.h:43
double aper1
beampipe information, new aperture model
Definition: element.h:107
double xsizeRight
individual collimator jaw half widths
Definition: element.h:127
double jawTiltRight
jaw collimator jaw tilts (angle in x-z plane)
Definition: element.h:128
std::string region
region with range cuts
Definition: element.h:209
double offsetX
offset X
Definition: element.h:129
double ysize
collimator aperture or laser spotsize for laser
Definition: element.h:125
double tilt
tilt
Definition: element.h:124
double l
length in metres
Definition: element.h:49
ElementType type
element enum
Definition: element.h:44
double offsetY
offset Y
Definition: element.h:130
std::string apertureType
beampipe information, new aperture model
Definition: element.h:111