BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSOutputROOTEventCollimator.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 "BDSOutputROOTEventCollimator.hh"
20#include "BDSOutputROOTParticleData.hh"
21
22#ifndef __ROOTBUILD__
23#include "BDSHitCollimator.hh"
24#include "BDSHitEnergyDeposition.hh"
25#include "BDSUtilities.hh"
26
27#include "globals.hh"
28#include "G4TwoVector.hh"
29
30#include "CLHEP/Units/SystemOfUnits.h"
31
32#include <utility>
33#endif
34
35#include <set>
36#include <vector>
37
39
41
42BDSOutputROOTParticleData* BDSOutputROOTEventCollimator::particleTable = nullptr;
43
45 primaryInteracted(false),
46 primaryStopped(false),
47 n(0),
48 totalEnergyDeposited(0)
49{;}
50
51BDSOutputROOTEventCollimator::~BDSOutputROOTEventCollimator()
52{;}
53
55{
56 primaryInteracted = false;
57 primaryStopped = false;
58 n = 0;
59 energy.clear();
60 energyDeposited.clear();
61 xIn.clear();
62 yIn.clear();
63 zIn.clear();
64 xpIn.clear();
65 ypIn.clear();
66 zpIn.clear();
67 T.clear();
68 weight.clear();
69 partID.clear();
70 parentID.clear();
71 turn.clear();
73 impactParameterX.clear();
74 impactParameterY.clear();
75 isIon.clear();
76 ionZ.clear();
77 ionA.clear();
78 turnSet.clear();
79 charge.clear();
80 kineticEnergy.clear();
81 mass.clear();
82 rigidity.clear();
84}
85
86#ifndef __ROOTBUILD__
87void BDSOutputROOTEventCollimator::Fill(const BDSHitCollimator* hit,
89 const std::pair<G4double, G4double>& differences,
90 G4bool storeHits)
91{
92 if (storeHits)
93 {
94 n++;
95 energy.push_back((float) (hit->totalEnergy / CLHEP::GeV));
96 const G4ThreeVector& pos = hit->preStepPosition;
97 const G4ThreeVector& mom = hit->preStepMomentum;
98 xIn.push_back((float) (pos.x() / CLHEP::m));
99 yIn.push_back((float) (pos.y() / CLHEP::m));
100 zIn.push_back((float) (pos.z() / CLHEP::m));
101 xpIn.push_back((float) (mom.x() / CLHEP::rad));
102 ypIn.push_back((float) (mom.y() / CLHEP::rad));
103 zpIn.push_back((float) (mom.z() / CLHEP::rad));
104
105 // calculate impact parameters - note done in output units (as is info object)
106 G4double impactX = std::abs(xIn.back() - info.offsetX);
107 G4double impactY = std::abs(yIn.back() - info.offsetY);
108 G4double impactZ = zIn.back();
109
110 // interpolate aperture to that point
111 G4double zFromStart = -0.5 * info.length - impactZ;
112 if (zFromStart < 0)
113 {zFromStart = 0;} // sometimes rounding problems
114 G4double fraction = zFromStart / info.length;
115 G4double xAperAtZ = info.xSizeIn + differences.first * fraction;
116 G4double yAperAtZ = info.ySizeIn + differences.second * fraction;
117
118 // impact parameter is absolute
119 impactX = impactX - xAperAtZ;
120 impactY = impactY - yAperAtZ;
121
122 if (BDS::IsFinite(info.tilt))
123 {
124 G4TwoVector impactPos(impactX, impactY);
125 impactPos.rotate(info.tilt);
126 impactX = impactPos.x();
127 impactY = impactPos.y();
128 }
129 impactParameterX.push_back((float) impactX);
130 impactParameterY.push_back((float) impactY);
131 }
132
134 if (eHit)
135 {
136 G4double eDep = eHit->GetEnergy() / CLHEP::GeV;
137 G4double eW = eHit->GetEnergyWeighted() / CLHEP::GeV;
138 G4double w = eHit->GetWeight();
139
141
142 primaryInteracted = primaryInteracted || eHit->GetParentID() == 0;
143
144 if (!storeHits)
145 {return;} // skip the rest of this function to avoid storing extra bits
146
147 energyDeposited.push_back((float)eDep);
148 T.push_back((float)(eHit->GetGlobalTime() / CLHEP::ns));
149 weight.push_back((float)w);
150 partID.push_back(eHit->GetPartID());
151 parentID.push_back(eHit->GetParentID());
152 G4int tn = eHit->GetTurnsTaken();
153 turn.push_back(tn);
154 // Mark if this is the first primary impact of this turn as there may
155 // be multiple 'hits' inside the collimator for a primary passing through
156 // but we're really interested in the first one. Judge by whether we have
157 // an hits for this turn already. We assume (safely) that the first entry
158 // for a given turn will be the primary.
159 if (parentID.back() == 0)
160 {firstPrimaryHitThisTurn.push_back(turnSet.find(tn) == turnSet.end());}
161 else
162 {firstPrimaryHitThisTurn.push_back(false);} // can't be true
163 turnSet.insert(tn);
164 }
165}
166
167void BDSOutputROOTEventCollimator::FillExtras(G4bool fillIonInfo,
168 G4bool fillLinks)
169{
170 if (!particleTable)
171 {return;}
172
173 if (!(fillIonInfo || fillLinks))
174 {return;}
175
176 for (int i = 0; i < n; ++i)
177 {// loop over all existing entries in the branch vectors
178 auto& pid = partID[i];
179 if (particleTable->IsIon(pid))
180 {
181 auto& ionInfo = particleTable->GetIonInfo(pid);
182 if (fillIonInfo)
183 {
184 isIon.push_back(true);
185 ionA.push_back(ionInfo.a);
186 ionZ.push_back(ionInfo.z);
187 }
188 if (fillLinks)
189 {
190 charge.push_back(ionInfo.charge);
191 mass.push_back((float)ionInfo.mass);
192 rigidity.push_back((float)ionInfo.rigidity(energy[i], ionInfo.charge));
193 kineticEnergy.push_back((float)particleTable->KineticEnergy(pid, energy[i]));
194 }
195 }
196 else
197 {// particle
198 auto& pInfo = particleTable->GetParticleInfo(pid);
199 if (fillIonInfo)
200 {
201 isIon.push_back(false);
202 ionA.push_back(0);
203 ionZ.push_back(0);
204 }
205 if (fillLinks)
206 {
207 charge.push_back(pInfo.charge);
208 mass.push_back((float)pInfo.mass);
209 rigidity.push_back((float)pInfo.rigidity(energy[i], pInfo.charge));
210 kineticEnergy.push_back((float)particleTable->KineticEnergy(pid, energy[i]));
211 }
212 }
213 }
214}
215
216#endif
217
218void BDSOutputROOTEventCollimator::Fill(const BDSOutputROOTEventCollimator* other)
219{
220 if (!other)
221 {return;}
222 primaryInteracted = other->primaryInteracted;
223 primaryStopped = other->primaryStopped;
224 n = other->n;
225 energy = other->energy;
226 energyDeposited = other->energyDeposited;
227
228 xIn = other->xIn;
229 yIn = other->yIn;
230 zIn = other->zIn;
231 xpIn = other->xpIn;
232 ypIn = other->ypIn;
233 zpIn = other->zpIn;
234 T = other->T;
235
236 weight = other->weight;
237 partID = other->partID;
238 parentID = other->parentID;
239 turn = other->turn;
240
242 impactParameterX = other->impactParameterX;
243 impactParameterY = other->impactParameterY;
244
245 isIon = other->isIon;
246 ionA = other->ionA;
247 ionZ = other->ionZ;
248 turnSet = other->turnSet;
249 charge = other->charge;
250
252 mass = other->mass;
253 rigidity = other->rigidity;
254}
Snapshot of information for particle passing through a collimator.
G4double totalEnergy
Total energy of particle.
G4ThreeVector preStepMomentum
Local pre step point momentum.
BDSHitEnergyDeposition * energyDepositionHit
G4ThreeVector preStepPosition
Local pre step point (z from centre of object).
Information recorded for a single piece of energy deposition.
G4int GetTurnsTaken() const
Accessor for extra piece of information.
G4int GetPartID() const
Accessor for extra piece of information.
G4int GetParentID() const
Accessor for extra piece of information.
G4double GetGlobalTime() const
Accessor for extra piece of information.
Data stored for each collimator in the model.
Data stored for each collimator per event.
std::vector< float > energy
Total energy of particle for each hit.
std::set< int > turnSet
Different length set of turn number.
std::vector< int > charge
These are not filled by default.
std::vector< float > mass
These are not filled by default.
virtual void Flush()
Flush this instance.
std::vector< float > kineticEnergy
These are not filled by default.
BDSOutputROOTEventCollimator()
Default constructor for ROOT.
double totalEnergyDeposited
Sum of energy deposits including weights.
std::vector< float > rigidity
These are not filled by default.
Geant4 particle data for particles used in simulation.
bool IsIon(const int &pdgID) const
This doesn't count a proton (even with electrons) as an ion.
double KineticEnergy(const int &pdgID, const double &totalEnergy) const
Calculate kinetic energy of particle given PDG ID and total energy in GeV.
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())