BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSOutputROOTEventModel.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 "BDSOutputROOTEventModel.hh"
20
21#ifndef __ROOTBUILD__
22#include "BDSAcceleratorModel.hh"
23#include "BDSBeamline.hh"
24#include "BDSBeamlineElement.hh"
25#include "BDSBeamPipeInfo.hh"
26#include "BDSMagnet.hh"
27#include "BDSMagnetStrength.hh"
28#include "BDSPhysicalVolumeInfoRegistry.hh"
29#include "BDSSamplerRegistry.hh"
30
31#include "G4RotationMatrix.hh"
32#include "G4String.hh"
33#include "G4Types.hh"
34
35#include <algorithm>
36#include <iterator>
37#include <map>
38#include <set>
39#include <sstream>
40#include <string>
41#include <vector>
42#include <utility>
43#endif
44
46
48 n(0),
49 storeCollimatorInfo(false),
50 storeCavityInfo(false),
51 nCavities(0),
52 nCollimators(0)
53{
54 Flush();
55}
56
58{;}
59
60int BDSOutputROOTEventModel::findNearestElement(const TVector3& point) const
61{
62 // we make a vector of pairs of the distance between the mid of each element
63 // and the specified point along with the index. We then sort this vector by
64 // the distance (first part of pair only). Could use a map as it's sorted but
65 // can't guarantee access with a double or float as a key due to binary representation.
66 std::vector<std::pair<float, int> > distanceAndIndex;
67 distanceAndIndex.reserve(midRefPos.size());
68 for (int i = 0; i < (int)midRefPos.size(); i++)
69 {distanceAndIndex.emplace_back(std::make_pair((midRefPos[i]-point).Mag(), i));}
70
71 struct customLess
72 {bool operator()(std::pair<float, int>& a, std::pair<float, int>& b) const {return a.first < b.first;};};
73 std::sort(distanceAndIndex.begin(), distanceAndIndex.end(), customLess());
74 return distanceAndIndex[0].second;
75}
76
78{
79 n = 0;
80 componentName.clear();
81 placementName.clear();
82 componentType.clear();
83 length.clear();
84 angle.clear();
85 staPos.clear();
86 midPos.clear();
87 endPos.clear();
88 staRot.clear();
89 midRot.clear();
90 endRot.clear();
91 staRefPos.clear();
92 midRefPos.clear();
93 endRefPos.clear();
94 staRefRot.clear();
95 midRefRot.clear();
96 endRefRot.clear();
97 tilt.clear();
98 offsetX.clear();
99 offsetY.clear();
100 staS.clear();
101 midS.clear();
102 endS.clear();
103 beamPipeType.clear();
104 beamPipeAper1.clear();
105 beamPipeAper2.clear();
106 beamPipeAper3.clear();
107 beamPipeAper4.clear();
108 material.clear();
109 k1.clear();
110 k2.clear();
111 k3.clear();
112 k4.clear();
113 k5.clear();
114 k6.clear();
115 k7.clear();
116 k8.clear();
117 k9.clear();
118 k10.clear();
119 k11.clear();
120 k12.clear();
121 k1s.clear();
122 k2s.clear();
123 k3s.clear();
124 k4s.clear();
125 k5s.clear();
126 k6s.clear();
127 k7s.clear();
128 k8s.clear();
129 k9s.clear();
130 k10s.clear();
131 k11s.clear();
132 k12s.clear();
133 ks.clear();
134 hkick.clear();
135 vkick.clear();
136 bField.clear();
137 eField.clear();
138 e1.clear();
139 e2.clear();
140 fint.clear();
141 fintx.clear();
142 fintk2.clear();
143 fintxk2.clear();
144 pvName.clear();
145 pvNameWPointer.clear();
146
147 storeCavityInfo = false;
148 cavityIndices.clear();
149 cavityIndicesByName.clear();
150 nCavities = 0;
151 cavityInfo.clear();
153 storeCollimatorInfo = false;
154 collimatorIndices.clear();
156 nCollimators = 0;
157 collimatorInfo.clear();
159 scoringMeshTranslation.clear();
160 scoringMeshRotation.clear();
161 scoringMeshName.clear();
162
163 materialIDToName.clear();
164 materialNameToID.clear();
165
166 samplerNamesUnique.clear();
167 samplerSPosition.clear();
168 samplerCNamesUnique.clear();
169 samplerSNamesUnique.clear();
170 samplerCRadius.clear();
171 samplerSRadius.clear();
172}
173
174#ifndef __ROOTBUILD__
176 G4bool storeCavityInfoIn):
177 n(0),
178 storeCollimatorInfo(storeCollimatorInfoIn),
179 storeCavityInfo(storeCavityInfoIn),
180 nCavities(0),
181 nCollimators(0)
182{;}
183
184TRotation BDSOutputROOTEventModel::ConvertToROOT(const G4RotationMatrix* rm) const
185{
186 TRotation rr = TRotation();
187 rr.SetToIdentity();
188 if (!rm)
189 {return rr;}
190
191 double rotAngle;
192 CLHEP::Hep3Vector axis;
193
194 rm->getAngleAxis(rotAngle, axis);
195 rr.Rotate(rotAngle, TVector3(axis.x(), axis.y(), axis.z()));
196 return rr;
197}
198
199TRotation BDSOutputROOTEventModel::ConvertToROOT(const G4RotationMatrix& rm) const
200{
201 return ConvertToROOT(&rm);
202}
203
204TVector3 BDSOutputROOTEventModel::ConvertToROOT(const G4ThreeVector& v) const
205{
206 return TVector3(v.x() / CLHEP::m, v.y() / CLHEP::m, v.z() / CLHEP::m);
207}
208
209void BDSOutputROOTEventModel::Fill(const std::vector<G4int>& collimatorIndicesIn,
210 const std::map<G4String, G4int>& collimatorIndicesByNameIn,
211 const std::vector<BDSOutputROOTEventCollimatorInfo>& collimatorInfoIn,
212 const std::vector<G4String>& collimatorBranchNamesIn,
213 const std::vector<G4int>& cavityIndicesIn,
214 const std::map<G4String, G4int>& cavityIndicesByNameIn,
215 const std::vector<BDSOutputROOTEventCavityInfo>& cavityInfoIn,
216 const std::vector<G4String>& cavityBranchNamesIn,
217 const std::map<G4String, G4Transform3D>* scorerMeshPlacements,
218 const std::map<short int, G4String>* materialIDToNameUnique,
219 G4bool storeTrajectory)
220{
222 for (const auto& nameSPos : sr->GetUniquePlaneNamesAndSPosition())
223 {
224 samplerNamesUnique.push_back(std::string(nameSPos.first) + ".");
225 samplerSPosition.push_back((double) nameSPos.second / CLHEP::m);
226 }
227 for (const auto& name : sr->GetUniqueNamesCylinder())
228 {samplerCNamesUnique.push_back(std::string(name) + ".");}
229 for (const auto& name : sr->GetUniqueNamesSphere())
230 {samplerSNamesUnique.push_back(std::string(name) + ".");}
231 samplerCRadius = sr->GetUniqueNameToRadiusCylinder();
232 samplerSRadius = sr->GetUniqueNameToRadiusSphere();
233
234 for (const auto& name : collimatorBranchNamesIn)
235 {collimatorBranchNamesUnique.push_back(std::string(name) + ".");}
236 for (const auto& name : cavityBranchNamesIn)
237 {cavityBranchNamesUnique.push_back(std::string(name) + ".");}
238
239 if (scorerMeshPlacements)
240 {
241 for (const auto& kv : *scorerMeshPlacements)
242 {
243 const G4String& name = kv.first;
244 scoringMeshTranslation[name] = ConvertToROOT(kv.second.getTranslation());
245 scoringMeshRotation[name] = ConvertToROOT(kv.second.getRotation());
246 scoringMeshName.emplace_back(name);
247 }
248 }
249
250 const BDSBeamline* beamline = BDSAcceleratorModel::Instance()->BeamlineMain();
251 if (!beamline)
252 {return;} // in case of generatePrimariesOnly there is no model - return
253
255 {
256 for (const auto value : collimatorIndicesIn)
257 {collimatorIndices.push_back((int)value);}
258
259 nCollimators = (int)collimatorIndices.size();
260
261 for (const auto& kv : collimatorIndicesByNameIn)
262 {collimatorIndicesByName[(std::string)kv.first] = (int)kv.second;}
263
264 collimatorInfo = collimatorInfoIn;
265 }
266
267 if (storeCavityInfo)
268 {
269 for (const auto value : cavityIndicesIn)
270 {cavityIndices.push_back((int)value);}
271
272 nCavities = (int)cavityIndices.size();
273
274 for (const auto& kv : cavityIndicesByNameIn)
275 {cavityIndicesByName[(std::string)kv.first] = (int)kv.second;}
276
277 cavityInfo = cavityInfoIn;
278 }
279
280 if (materialIDToNameUnique && storeTrajectory)
281 {
282 for (const auto& kv : *materialIDToNameUnique)
283 {
284 materialIDToName[kv.first] = (std::string)kv.second;
285 materialNameToID[(std::string)kv.second] = kv.first;
286 }
287 }
288
289 n = (int)beamline->size();
290
291 for (auto i = beamline->begin(); i != beamline->end(); ++i)
292 {
293 componentName.push_back((*i)->GetName());
294 placementName.push_back((*i)->GetPlacementName());
295 componentType.push_back((*i)->GetType());
296 length.push_back((float)((*i)->GetAcceleratorComponent()->GetArcLength() / CLHEP::m));
297 angle.push_back((float)((*i)->GetAcceleratorComponent()->GetAngle() / CLHEP::radian));
298 staPos.push_back(ConvertToROOT((*i)->GetPositionStart()));
299 midPos.push_back(ConvertToROOT((*i)->GetPositionMiddle()));
300 endPos.push_back(ConvertToROOT((*i)->GetPositionEnd()));
301 staRot.push_back(ConvertToROOT((*i)->GetRotationStart()));
302 midRot.push_back(ConvertToROOT((*i)->GetRotationMiddle()));
303 endRot.push_back(ConvertToROOT((*i)->GetRotationEnd()));
304 staRefPos.emplace_back(ConvertToROOT((*i)->GetReferencePositionStart()));
305 midRefPos.emplace_back(ConvertToROOT((*i)->GetReferencePositionMiddle()));
306 endRefPos.emplace_back(ConvertToROOT((*i)->GetReferencePositionEnd()));
307 staRefRot.push_back(ConvertToROOT((*i)->GetReferenceRotationStart()));
308 midRefRot.push_back(ConvertToROOT((*i)->GetReferenceRotationMiddle()));
309 endRefRot.push_back(ConvertToROOT((*i)->GetReferenceRotationEnd()));
310
311 // tilt and offset
312 BDSTiltOffset* to = (*i)->GetTiltOffset();
313 if (to)
314 {
315 tilt.push_back((float)(to->GetTilt() / CLHEP::rad));
316 offsetX.push_back((float)(to->GetXOffset() / CLHEP::m));
317 offsetY.push_back((float)(to->GetYOffset() / CLHEP::m));
318 }
319 else
320 {
321 tilt.push_back(0);
322 offsetX.push_back(0);
323 offsetY.push_back(0);
324 }
325
326 // S positions
327 staS.push_back((float)((*i)->GetSPositionStart() / CLHEP::m));
328 midS.push_back((float)((*i)->GetSPositionMiddle() / CLHEP::m));
329 endS.push_back((float)((*i)->GetSPositionEnd() / CLHEP::m));
330
331 // beam pipe
332 BDSBeamPipeInfo* beampipeinfo = (*i)->GetBeamPipeInfo();
333 beamPipeType.push_back(beampipeinfo ? beampipeinfo->beamPipeType.ToString() : "");
334 beamPipeAper1.push_back(beampipeinfo ? beampipeinfo->aper1 / CLHEP::m : 0);
335 beamPipeAper2.push_back(beampipeinfo ? beampipeinfo->aper2 / CLHEP::m : 0);
336 beamPipeAper3.push_back(beampipeinfo ? beampipeinfo->aper3 / CLHEP::m : 0);
337 beamPipeAper4.push_back(beampipeinfo ? beampipeinfo->aper4 / CLHEP::m : 0);
338
339 // associated material if any
340 const auto accComp = (*i)->GetAcceleratorComponent();
341 material.push_back(accComp->Material());
342
343 // helper shortcuts to all the member vectors
344 std::vector<std::vector<float>*> localNorm = {&k1,&k2,&k3,&k4,&k5,&k6,&k7,&k8,&k9,&k10,&k11,&k12};
345 std::vector<std::vector<float>*> localSkew = {&k1s,&k2s,&k3s,&k4s,&k5s,&k6s,&k7s,&k8s,&k9s,&k10s,&k11s,&k12s};
346
347 // helper lambda to avoid duplication
348 auto fillzero = [&]
349 {
350 for (int j = 0; j < (int)localNorm.size(); j++)
351 {localNorm[j]->push_back(0);}
352 for (int j = 0; j < (int)localSkew.size(); j++)
353 {localSkew[j]->push_back(0);}
354 ks.push_back(0);
355 hkick.push_back(0);
356 vkick.push_back(0);
357 bField.push_back(0);
358 eField.push_back(0);
359 e1.push_back(0);
360 e2.push_back(0);
361 hgap.push_back(0);
362 fint.push_back(0);
363 fintx.push_back(0);
364 fintk2.push_back(0);
365 fintxk2.push_back(0);
366 };
367
368 // do this bit first as we test for magnet strengths later and then do a 'continue' in the for loop
370 std::vector<std::string> localPVNames;
371 std::vector<std::string> localPVNamesWPointer;
372 if (setOfPVs)
373 {
374 auto name = [](G4VPhysicalVolume* pv){return pv->GetName();};
375 std::transform(setOfPVs->begin(), setOfPVs->end(), std::back_inserter(localPVNames), name);
376 auto fullName = [](G4VPhysicalVolume* pv)
377 {
378 std::stringstream ss;
379 ss << static_cast<const void*>(pv);
380 std::string addressName = ss.str();
381 return pv->GetName() + addressName;
382 };
383 std::transform(setOfPVs->begin(), setOfPVs->end(), std::back_inserter(localPVNamesWPointer), fullName);
384 }
385 // always push it back even if an empty vector
386 pvName.push_back(localPVNames);
387 pvNameWPointer.push_back(localPVNamesWPointer);
388
389 // fill magnet strength data
390 // NOTE - has a 'continue'
391 if (BDSMagnet* mag = dynamic_cast<BDSMagnet*>(accComp))
392 {
393 const BDSMagnetStrength* ms = mag->MagnetStrength();
394 if (!ms)
395 {
396 fillzero();
397 continue;
398 }
399 // assume localNorm and normComponents are same size
400 std::vector<G4double> normComponents = ms->NormalComponents();
401 for (int j = 0; j < (int)localNorm.size(); j++)
402 {localNorm[j]->push_back((float)normComponents[j]);}
403 std::vector<G4double> skewComponents = ms->SkewComponents();
404 for (int j = 0; j < (int)localSkew.size(); j++)
405 {localSkew[j]->push_back((float)skewComponents[j]);}
406
407 ks.push_back((float)((*ms)["ks"]/BDSMagnetStrength::Unit("ks")));
408 hkick.push_back( (float)((*ms)["hkick"]/BDSMagnetStrength::Unit("hkick")));
409 vkick.push_back( (float)((*ms)["vkick"]/BDSMagnetStrength::Unit("vkick")));
410 bField.push_back((float)((*ms)["field"]/BDSMagnetStrength::Unit("field")));
411 eField.push_back((float)((*ms)["efield"]/BDSMagnetStrength::Unit("efield")));
412 e1.push_back((float)((*ms)["e1"]/BDSMagnetStrength::Unit("e1")));
413 e2.push_back((float)((*ms)["e2"]/BDSMagnetStrength::Unit("e2")));
414 hgap.push_back((float)((*ms)["hgap"]/BDSMagnetStrength::Unit("hgap")));
415 fint.push_back((float)(*ms)["fint"]);
416 fintx.push_back((float)(*ms)["fintx"]);
417 fintk2.push_back((float)(*ms)["fintk2"]);
418 fintxk2.push_back((float)(*ms)["fintxk2"]);
419 }
420 else
421 {// not a magnet
422 fillzero();
423 }
424 }
425}
426#endif
const BDSBeamline * BeamlineMain() const
Accessor.
Holder class for all information required to describe a beam pipe model.
G4double aper3
Public member for direct access.
G4double aper1
Public member for direct access.
G4double aper4
Public member for direct access.
G4double aper2
Public member for direct access.
BDSBeamPipeType beamPipeType
Public member for direct access.
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
iterator begin()
Iterator mechanics.
Definition: BDSBeamline.hh:167
iterator end()
Iterator mechanics.
Definition: BDSBeamline.hh:168
Efficient storage of magnet strengths.
std::vector< G4double > NormalComponents() const
Accessor for all normal components - k1 - k12.
std::vector< G4double > SkewComponents() const
Accessor for all skew components - k1 - k12.
static G4double Unit(const G4String &key)
Access a unit factor for a given key.
Abstract base class that implements features common to all magnets.
Definition: BDSMagnet.hh:45
Information stored per model representing accelerator.
std::vector< BDSOutputROOTEventCollimatorInfo > collimatorInfo
Collimator information explicitly.
bool storeCollimatorInfo
Whether optional collimator information was stored.
std::vector< float > vkick
Vertical fractional momentum kick.
std::vector< float > eField
E field in V/m.
std::vector< std::string > collimatorBranchNamesUnique
Vector of all collimator branch names in event tree used to load data.
int nCollimators
Number of collimators in beam line.
int nCavities
Number of cavities in beam line.
virtual ~BDSOutputROOTEventModel()
Destructor.
int findNearestElement(const TVector3 &point) const
Find element index closest to point.
std::vector< float > bField
B field in T.
std::vector< float > hkick
Horizontal fractional momentum kick.
std::vector< int > collimatorIndices
std::vector< std::string > material
Material associated with element if any.
std::vector< BDSOutputROOTEventCavityInfo > cavityInfo
cavity information explicitly.
void Flush()
Initialise all members.
virtual void Fill(const std::vector< G4int > &collimatorIndicesIn={}, const std::map< G4String, G4int > &collimatorIndicesByNameIn={}, const std::vector< BDSOutputROOTEventCollimatorInfo > &collimatorInfoIn={}, const std::vector< G4String > &collimatorBranchNamesIn={}, const std::vector< G4int > &cavityIndicesIn={}, const std::map< G4String, G4int > &cavityIndicesByNameIn={}, const std::vector< BDSOutputROOTEventCavityInfo > &cavityInfoIn={}, const std::vector< G4String > &cavityBranchNamesIn={}, const std::map< G4String, G4Transform3D > *scorerMeshPlacements=nullptr, const std::map< short int, G4String > *materialIDToNameUnique=nullptr, G4bool storeTrajectory=false)
Fill root output.
std::vector< std::string > cavityBranchNamesUnique
Vector of all cavity branch names in event tree used to load data.
BDSOutputROOTEventModel()
Default constructor.
std::map< std::string, int > collimatorIndicesByName
Similar cache but by name of collimator as built by BDSIM.
std::vector< float > ks
Solenoid strength.
bool storeCavityInfo
Whether optional cavity information was stored.
TRotation ConvertToROOT(const G4RotationMatrix *rm) const
Utility function.
std::map< std::string, int > cavityIndicesByName
Similar cache but by name of cavity as built by BDSIM.
const std::set< G4VPhysicalVolume * > * PVsForBeamlineElement(BDSBeamlineElement *element) const
Access a set of volumes registered for the placement of a beamline element.
static BDSPhysicalVolumeInfoRegistry * Instance()
Singleton accessor.
static BDSSamplerRegistry * Instance()
Accessor for registry.
A holder for any placement offsets and rotations for a BDSAcceleratorComponent.
G4double GetYOffset() const
Accessor.
G4double GetXOffset() const
Accessor.
G4double GetTilt() const
Accessor.