BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSOutputROOTEventTrajectory.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 <algorithm>
20#include <bitset>
21#include <iostream>
22#include <iomanip>
23
24#include "BDSOutputROOTEventTrajectory.hh"
25
26#ifndef __ROOTBUILD__
27#include "G4VPhysicalVolume.hh"
28
29#include "BDSHitEnergyDeposition.hh"
30#include "BDSAuxiliaryNavigator.hh"
31#include "BDSPhysicalVolumeInfoRegistry.hh"
32#include "BDSPhysicalVolumeInfo.hh"
33#include "BDSTrajectory.hh"
34#include "BDSTrajectoryOptions.hh"
35
36#include <cmath>
37#include <map>
38#endif
39
41BDSOutputROOTEventTrajectory::BDSOutputROOTEventTrajectory():
42 auxNavigator(nullptr),
43 n(0)
44{;}
45
46BDSOutputROOTEventTrajectory::~BDSOutputROOTEventTrajectory()
47{
48#ifndef __ROOTBUILD__
49 delete auxNavigator;
50#endif
51}
52
53#ifndef __ROOTBUILD__
54int findPrimaryStepIndex(BDSTrajectory* traj)
55{
56 if (!traj->GetParent())
57 {return -1;}
58
59 if (traj->GetParent()->GetTrackID() == 1)
60 {return traj->GetParentStepIndex();}
61 else
62 {return findPrimaryStepIndex(traj->GetParent());}
63}
64
65void BDSOutputROOTEventTrajectory::Fill(const BDSTrajectoriesToStore* trajectories,
66 int storeStepPointsN,
67 bool storeStepPointLast,
68 const BDS::TrajectoryOptions& storageOptions,
69 const std::map<G4Material*, short int>& materialToID)
70{
71 if (!auxNavigator)
72 {// navigator for converting coordinates to curvilinear coordinate system
74 }
75
76 G4bool stEK = storageOptions.storeLinks || storageOptions.storeKineticEnergy;
77 G4bool stMo = storageOptions.storeMomentumVector;
78 G4bool stPr = storageOptions.storeProcesses;
79 G4bool stTi = storageOptions.storeTime;
80 G4bool stMa = storageOptions.storeMaterial;
81
82 // assign trajectory indices
83 int idx = 0;
84 for (auto trajFlag : trajectories->trajectories)
85 {
86 BDSTrajectory* traj = trajFlag.first;
87 if (trajFlag.second) // ie we want to save this trajectory
88 {
89 traj->SetTrajIndex(idx);
90 idx++;
91 }
92 else // we don't want to save it
93 {traj->SetTrajIndex(-1);}
94 }
95
96 // assign parent (and step) indices
97 for (auto trajFlag : trajectories->trajectories)
98 {
99 BDSTrajectory* traj = trajFlag.first;
100 BDSTrajectory* parent = traj->GetParent();
101 if (trajFlag.second && parent)
102 { // to store and not primary
103 traj->SetParentIndex(parent->GetTrajIndex());
104
105 // search for parent step index
106 if (parent->GetTrajIndex() != -1)
107 {
108 auto trajStartPos = traj->GetPoint(0)->GetPosition();
109 traj->SetParentStepIndex(-1);
110 for (int i = 0; i < parent->GetPointEntries(); ++i)
111 {
112 if(parent->GetPoint(i)->GetPosition() == trajStartPos)
113 {
114 traj->SetParentStepIndex(i);
115 break;
116 }
117 }
118 }
119 else
120 {parent->SetParentStepIndex(-1);}
121 }
122 else
123 {traj->SetParentIndex(-1);}
124 }
125
126 n = 0;
127 for (auto trajFlag : trajectories->trajectories)
128 {
129 BDSTrajectory* traj = trajFlag.first;
130
131 // check if the trajectory is to be stored
132 if (!trajFlag.second) // ie false, then continue and don't store
133 {continue;}
134
135 partID.push_back((int) traj->GetPDGEncoding());
136 trackID.push_back((unsigned int) std::abs(traj->GetTrackID()));
137 parentID.push_back((unsigned int) std::abs(traj->GetParentID()));
138 parentIndex.push_back((unsigned int) std::abs(traj->GetParentIndex()));
139 parentStepIndex.push_back((unsigned int) std::abs(traj->GetParentStepIndex()));
140 depth.push_back((int) traj->GetDepth());
141
142 // now we convert the geant4 type based BDSTrajectory information into
143 // basic C++ and ROOT types for the output
144 IndividualTrajectory itj;
145 if (storeStepPointsN > 0)
146 {// store specific number of step points along the trajectory
147 G4int nSteps = traj->GetPointEntries();
148 G4int nPoints = std::min(nSteps, storeStepPointsN);
149 for (int i = 0; i < nPoints; ++i)
150 {FillIndividualTrajectory(itj, traj, i, materialToID);}
151 // optionally include the last point if required and not already stored
152 if (storeStepPointLast && (nPoints < nSteps))
153 {FillIndividualTrajectory(itj, traj, nSteps-1, materialToID);}
154 }
155 else
156 {// store all points as usual
157 for (int i = 0; i < traj->GetPointEntries(); ++i)
158 {FillIndividualTrajectory(itj, traj, i, materialToID);}
159 }
160
161 // record the filters that were matched for this trajectory
162 filters.push_back(trajectories->filtersMatched.at(traj));
163
164 XYZ.push_back(itj.XYZ);
165 modelIndicies.push_back(itj.modelIndex);
166
167 if (stMo)
168 {PXPYPZ.push_back(itj.PXPYPZ);}
169
170 S.push_back(itj.S);
171
172 if (stPr)
173 {
174 preProcessTypes.push_back(itj.preProcessType);
175 preProcessSubTypes.push_back(itj.preProcessSubType);
176 postProcessTypes.push_back(itj.postProcessType);
177 postProcessSubTypes.push_back(itj.postProcessSubType);
178 }
179
180 preWeights.push_back(itj.preWeight);
181 postWeights.push_back(itj.postWeight);
182 energyDeposit.push_back(itj.energyDeposit);
183
184 if (stTi)
185 {T.push_back(itj.T);}
186
187 if (stEK)
188 {kineticEnergy.push_back(itj.kineticEnergy);}
189
190 if (stMa)
191 {materialID.push_back(itj.materialID);}
192
193 if (!itj.xyz.empty())
194 {
195 xyz.push_back(itj.xyz);
196 pxpypz.push_back(itj.pxpypz);
197 }
198
199 if (!itj.charge.empty())
200 {
201 charge.push_back(itj.charge);
202 turnsTaken.push_back(itj.turn);
203 mass.push_back(itj.mass);
204 rigidity.push_back(itj.rigidity);
205 }
206
207 if (!itj.isIon.empty())
208 {
209 isIon.push_back(itj.isIon);
210 ionA.push_back(itj.ionA);
211 ionZ.push_back(itj.ionZ);
212 nElectrons.push_back(itj.nElectrons);
213 }
214
215 // recursively search for primary interaction step
216 primaryStepIndex.push_back(findPrimaryStepIndex(traj));
217
218 // geant4 trackID to trackIndex in this table
219 trackID_trackIndex.insert(std::pair<int,int>(traj->GetTrackID(),n));
220
221 n++;
222 }
223
224#if 0
225 // Fill maps for later analysis
226 int trackIndex = 0;
227 for (auto iT = trajVec.begin(); iT != trajVec.end(); ++iT)
228 {
229 BDSTrajectory* traj = *iT;
230
231 // map of trackID to trackIndex
232 trackID_trackIndex.insert(std::pair<int, int>(traj->GetTrackID(),trackIndex));
233
234 std::cout << trajVec.size() << " " << parentID.size() << " " << parentIndex.size() << " "
235 << traj->GetTrackID() << " " << traj->GetParentID() << " " << trackIndex << std::endl;
236
237 // map of trackIndex to trackProcess
238 auto processPair = findParentProcess(trackIndex);
239 trackIndex_trackProcess.insert(std::pair<int,std::pair<int,int>>(trackIndex,processPair));
240
241 // map of modelIndex to trackProcess
242 if(processPair.first != -1)
243 {
244 int mi = modelIndicies[processPair.first][processPair.second];
245 trackIndex_modelIndex.insert(std::pair<int,int>(trackIndex, mi));
246 try
247 {modelIndex_trackIndex.at(mi).push_back(trackIndex);}
248 catch(const std::exception&)
249 {
250 modelIndex_trackIndex.insert(std::pair<int,std::vector<int>>(mi,std::vector<int>()));
251 modelIndex_trackIndex.at(mi).push_back(trackIndex);
252 }
253 }
254
255 ++trackIndex;
256 }
257#endif
258}
259
261 BDSTrajectory* traj,
262 int i,
263 const std::map<G4Material*, short int>& materialToID) const
264{
265 BDSTrajectoryPoint* point = dynamic_cast<BDSTrajectoryPoint*>(traj->GetPoint(i));
266 if (!point)
267 {return;}
268
269 // Position
270 G4ThreeVector pos = point->GetPosition();
271 itj.XYZ.emplace_back(TVector3(pos.getX() / CLHEP::m,
272 pos.getY() / CLHEP::m,
273 pos.getZ() / CLHEP::m));
274
275 G4VPhysicalVolume* vol = auxNavigator->LocateGlobalPointAndSetup(pos,nullptr,true,true,true);
277 if (theInfo)
278 {itj.modelIndex.push_back(theInfo->GetBeamlineIndex());}
279 else
280 {itj.modelIndex.push_back(-1);}
281
282 // Process types
283 itj.preProcessType.push_back(point->GetPreProcessType());
284 itj.preProcessSubType.push_back(point->GetPreProcessSubType());
285 itj.postProcessType.push_back(point->GetPostProcessType());
286 itj.postProcessSubType.push_back(point->GetPostProcessSubType());
287
288 itj.preWeight.push_back(point->GetPreWeight());
289 itj.postWeight.push_back(point->GetPostWeight());
290 itj.energyDeposit.push_back(point->GetEnergyDeposit() / CLHEP::GeV);
291 G4ThreeVector mom = point->GetPreMomentum() / CLHEP::GeV;
292 itj.PXPYPZ.emplace_back(TVector3(mom.getX(),
293 mom.getY(),
294 mom.getZ()));
295 itj.S.push_back(point->GetPreS() / CLHEP::m);
296 itj.T.push_back(point->GetPreGlobalTime() / CLHEP::ns);
297 itj.kineticEnergy.push_back(point->GetKineticEnergy() / CLHEP::GeV);
298
299 itj.materialID.push_back(materialToID.at(point->GetMaterial()));
300
301 if (point->extraLocal)
302 {
303 G4ThreeVector localPos = point->GetPositionLocal() / CLHEP::m;
304 G4ThreeVector localMom = point->GetMomentumLocal() / CLHEP::GeV;
305 itj.xyz.emplace_back(TVector3(localPos.getX(),
306 localPos.getY(),
307 localPos.getZ()));
308 itj.pxpypz.emplace_back(TVector3(localMom.getX(),
309 localMom.getY(),
310 localMom.getZ()));
311 }
312
313 if (point->extraLink)
314 {
315 itj.charge.push_back((int) (point->GetCharge() / (G4double)CLHEP::eplus));
316 itj.turn.push_back(point->GetTurnsTaken());
317 itj.mass.push_back(point->GetMass() / CLHEP::GeV);
318 itj.rigidity.push_back(point->GetRigidity() / (CLHEP::tesla*CLHEP::m));
319 }
320
321 if (point->extraIon)
322 {
323 itj.isIon.push_back(point->GetIsIon());
324 itj.ionA.push_back(point->GetIonA());
325 itj.ionZ.push_back(point->GetIonZ());
326 itj.nElectrons.push_back(point->GetNElectrons());
327 }
328}
329
330void BDSOutputROOTEventTrajectory::Fill(const BDSHitsCollectionEnergyDeposition* phc)
331{
332 G4cout << phc->GetSize() << G4endl;
333}
334
335#endif
336
338{
339 n = 0;
340 filters.clear();
341 partID.clear();
342 trackID.clear();
343 parentID.clear();
344 parentIndex.clear();
345 parentStepIndex.clear();
346 primaryStepIndex.clear();
347 depth.clear();
348 preProcessTypes.clear();
349 preProcessSubTypes.clear();
350 postProcessTypes.clear();
351 postProcessSubTypes.clear();
352 preWeights.clear();
353 postWeights.clear();
354 energyDeposit.clear();
355 XYZ.clear();
356 S.clear();
357 PXPYPZ.clear();
358 T.clear();
359 xyz.clear();
360 pxpypz.clear();
361 charge.clear();
362 kineticEnergy.clear();
363 turnsTaken.clear();
364 mass.clear();
365 rigidity.clear();
366 isIon.clear();
367 ionA.clear();
368 ionZ.clear();
369 nElectrons.clear();
370 materialID.clear();
371 modelIndicies.clear();
372 trackID_trackIndex.clear();
373
374 // trackIndex_trackProcess.clear();
375 // trackIndex_modelIndex.clear();
376 // modelIndex_trackIndex.clear();
377}
378
379void BDSOutputROOTEventTrajectory::Fill(const BDSOutputROOTEventTrajectory* other)
380{
381 if (!other)
382 {return;}
383
384 n = 0;
385 filters = other->filters;
386 partID = other->partID;
387 trackID = other->trackID;
388 parentID = other->parentID;
389 parentIndex = other->parentIndex;
390 parentStepIndex = other->parentStepIndex;
391 primaryStepIndex = other->primaryStepIndex;
392 depth = other->depth;
393 preProcessTypes = other->preProcessTypes;
394 preProcessSubTypes = other->preProcessSubTypes;
395 postProcessTypes = other->postProcessTypes;
396 postProcessSubTypes = other->postProcessSubTypes;
397 preWeights = other->preWeights;
398 postWeights = other->postWeights;
399 energyDeposit = other->energyDeposit;
400 XYZ = other->XYZ;
401 S = other->S;
402 PXPYPZ = other->PXPYPZ;
403 modelIndicies = other->modelIndicies;
404 trackID_trackIndex = other->trackID_trackIndex;
405 T = other->T;
406
407 xyz = other->xyz;
408 pxpypz = other->pxpypz;
409 charge = other->charge;
411 turnsTaken = other->turnsTaken;
412 rigidity = other->rigidity;
413 isIon = other->isIon;
414 ionA = other->ionA;
415 ionZ = other->ionZ;
416 nElectrons = other->nElectrons;
417 materialID = other->materialID;
418}
419
420#if 0
421std::pair<int,int> BDSOutputROOTEventTrajectory::findParentProcess(int trackIndex)
422{
423 std::cout << "BDSOutputROOTEventTrajectory::findParentProcess> "
424 << trackIndex << " " << parentID.size() << " " << parentIndex.size() << std::endl;
425 int tid = trackIndex;
426 int pid = parentID.at(tid);
427 std::cout << pid << std::endl;
428 int pin = parentIndex.at(tid);
429 std::cout << pin << std::endl;
430
431 if (pin == -1)
432 {return std::pair<int,int>(-1,-1);}
433 int sin = parentStepIndex.at(tid);
434 std::cout << sin << std::endl;
435
436 while (pid > 0)
437 {
438 if(pin == 0)
439 {break;}
440 tid = pin;
441 pid = parentID.at(tid);
442 pin = parentIndex.at(tid);
443 sin = parentStepIndex.at(tid);
444
445 std::cout << tid << " " << pid << " " << pin << " " << sin << " " << std::endl;
446 }
447
448 return std::pair<int,int>(pin,sin);
449}
450#endif
451
452std::vector<BDSOutputROOTEventTrajectoryPoint> BDSOutputROOTEventTrajectory::trackInteractions(int trackIDIn) const
453{
454 // prevent a bad access
455 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
456 {
457 std::cout << "No such track ID" << std::endl;
458 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
459 }
460
461 int ti = trackID_trackIndex.at(trackIDIn); // get track index
462
463 std::vector<BDSOutputROOTEventTrajectoryPoint> tpv; // trajectory point vector - result
464
465 int nstep = (int)XYZ[ti].size();
466 if (nstep == 0) // no points or it wasn't stored
467 {return std::vector<BDSOutputROOTEventTrajectoryPoint>();}
468
469 if (postProcessTypes[ti].empty()) // and implicitly nstep>0
470 {// we require processes for this function
471 std::cout << "Processes not stored for this file - not possible." << std::endl;
472 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
473 }
474
475 bool usePXPYPZ = !PXPYPZ[ti].empty();
476 bool useT = !T.empty();
477 bool useIon = !isIon.empty();
478 bool useLocal = !xyz.empty();
479 bool useLinks = !charge.empty();
480 bool useEK = !kineticEnergy.empty();
481 bool useMat = !materialID.empty();
482
483 for (int i = 0; i < nstep; ++i)
484 {
485 int ppt = postProcessTypes[ti][i];
486 int ppst = postProcessSubTypes[ti][i];
487 // this is a hard coded version of BDSTrajectoryPoint::IsScatteringPoint which is
488 // only available if we link to Geant4 for the enums (we don't here).
489 // -1 = undefined, 1 = G4ProcessType::fTransportation, 10 = G4ProcessTypes::fParallel
490 // 0 = G4ProcessType::fNotDefined for crystal channeling (exclude the thousands of points)
491 // 401 = G4TransportationProcessSubType::STEP_LIMITER which is categorised under G4ProcessType::fGeneral
492 bool notGeneral = ppt != 7 && ppst != 401;
493 bool changeInEnergy = energyDeposit[ti][i] > 1e-9;
494 if ( (ppt != -1 && ppt != 1 && ppt != 10 && ppt != 0 && notGeneral) || changeInEnergy)
495 {
497 trackID[ti],
498 parentID[ti],
499 parentIndex[ti],
500 postProcessTypes[ti][i],
501 postProcessSubTypes[ti][i],
502 postWeights[ti][i],
503 energyDeposit[ti][i],
504 XYZ[ti][i],
505 usePXPYPZ ? PXPYPZ[ti][i] : TVector3(),
506 modelIndicies[ti][i],
507 useT ? T[ti][i] : 0,
508 useLocal ? xyz[ti][i] : TVector3(),
509 useLocal ? pxpypz[ti][i] : TVector3(),
510 useLinks ? charge[ti][i] : 0,
511 useEK ? kineticEnergy[ti][i] : 0,
512 useLinks ? turnsTaken[ti][i] : 0,
513 useLinks ? rigidity[ti][i] : 0,
514 useLinks ? mass[ti][i] : 0,
515 useIon ? isIon[ti][i] : false,
516 useIon ? ionA[ti][i] : 0,
517 useIon ? ionZ[ti][i] : 0,
518 useIon ? nElectrons[ti][i] : 0,
519 useMat ? materialID[ti][i] : -1,
520 i);
521 tpv.push_back(p);
522 }
523 }
524 return tpv;
525}
526
527BDSOutputROOTEventTrajectoryPoint BDSOutputROOTEventTrajectory::primaryProcessPoint(int trackIDIn) const
528{
529 // prevent a bad access
530 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
531 {
532 std::cout << "No such track ID" << std::endl;
534 }
535 int ti = trackID_trackIndex.at(trackIDIn); // get track index
536 int pid = (int)parentID[ti]; // parent trackID
537 int chosenTrackID = trackIDIn;
538 while (pid != 0)
539 {
540 if (parentID[trackID_trackIndex.at(pid)] > 0)
541 {chosenTrackID = pid;}
542 ti = trackID_trackIndex.at(pid);
543 pid = (int)parentID[ti];
544 }
545 return parentProcessPoint(chosenTrackID);
546}
547
548BDSOutputROOTEventTrajectoryPoint BDSOutputROOTEventTrajectory::parentProcessPoint(int trackIDIn) const
549{
550 // prevent a bad access
551 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
552 {
553 std::cout << "No such track ID" << std::endl;
555 }
556
557 int ti = trackID_trackIndex.at(trackIDIn); // get track index
558 int pti = (int)parentID[ti]; // parent trackID
559
560 if (pti == 0)
561 {
562 std::cout << "Track is a parent" << std::endl;
564 }
565
566 int si = (int)parentStepIndex.at(ti); // get primary index
567 int pi = 0;
568 if (pti > 0)
569 {pi = trackID_trackIndex.at(pti);} // parent track storage index
570 else
571 {pi = ti;}
572
573 if (si > (int)XYZ[pi].size())
574 {// evidently not all step points are stored
575 std::cout << "Not all step points are stored. Parent step index is outside points stored." << std::endl;
577 }
578
579 bool usePXPYPZ = !PXPYPZ[ti].empty();
580 bool useT = !T.empty();
581 bool useIon = !isIon.empty();
582 bool useLocal = !xyz.empty();
583 bool useLinks = !charge.empty();
584 bool useEK = !kineticEnergy.empty();
585 bool useMat = !materialID.empty();
586
588 trackID[pi],
589 parentID[pi],
590 parentIndex[pi],
591 postProcessTypes[pi][si],
592 postProcessSubTypes[pi][si],
593 postWeights[pi][si],
594 energyDeposit[pi][si],
595 XYZ[pi][si],
596 usePXPYPZ ? PXPYPZ[pi][si] : TVector3(),
597 modelIndicies[pi][si],
598 useT ? T[pi][si] : 0,
599 useLocal ? xyz[pi][si] : TVector3(),
600 useLocal ? pxpypz[pi][si] : TVector3(),
601 useLinks ? charge[pi][si] : 0,
602 useEK ? kineticEnergy[pi][si] : 0,
603 useLinks ? turnsTaken[pi][si] : 0,
604 useLinks ? rigidity[pi][si] : 0,
605 useLinks ? mass[pi][si] : 0,
606 useIon ? isIon[pi][si] : false,
607 useIon ? ionA[pi][si] : 0,
608 useIon ? ionZ[pi][si] : 0,
609 useIon ? nElectrons[pi][si] : 0,
610 useMat ? materialID[pi][si] : -1,
611 si);
612 return p;
613}
614
615std::vector<BDSOutputROOTEventTrajectoryPoint> BDSOutputROOTEventTrajectory::processHistory(int trackIDIn) const
616{
617 // prevent a bad access
618 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
619 {
620 std::cout << "No such track ID" << std::endl;
621 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
622 }
623
624 if (postProcessTypes.empty())
625 {
626 std::cout << "Processes not stored for this file - not possible." << std::endl;
627 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
628 }
629
630 int ti = trackID_trackIndex.at(trackIDIn);
631
632 bool usePXPYPZ = !PXPYPZ.empty();
633 bool useT = !T.empty();
634 bool useIon = !isIon.empty();
635 bool useLocal = !xyz.empty();
636 bool useLinks = !charge.empty();
637 bool useEK = !kineticEnergy.empty();
638 bool useMat = !materialID.empty();
639
640 std::vector<BDSOutputROOTEventTrajectoryPoint> tpv; // trajectory point vector
641 while (ti != 0)
642 {
643 unsigned int pi = parentIndex.at(ti);
644 unsigned int psi = parentStepIndex.at(ti);
645 if (psi > (unsigned int)XYZ[pi].size())
646 {
647 std::cout << "Not all points stored - defaulting to creation point." << std::endl;
648 psi = 0;
649 }
650
652 trackID[pi],
653 parentID[pi],
654 parentIndex[pi],
655 postProcessTypes[pi][psi],
656 postProcessSubTypes[pi][psi],
657 postWeights[pi][psi],
658 energyDeposit[pi][psi],
659 XYZ[pi][psi],
660 usePXPYPZ ? PXPYPZ[ti][psi] : TVector3(),
661 modelIndicies[ti][psi],
662 useT ? T[ti][psi] : 0,
663 useLocal ? xyz[ti][psi] : TVector3(),
664 useLocal ? pxpypz[ti][psi] : TVector3(),
665 useLinks ? charge[ti][psi] : 0,
666 useEK ? kineticEnergy[ti][psi] : 0,
667 useLinks ? turnsTaken[ti][psi] : 0,
668 useLinks ? rigidity[ti][psi] : 0,
669 useLinks ? mass[ti][psi] : 0,
670 useIon ? isIon[ti][psi] : false,
671 useIon ? ionA[ti][psi] : 0,
672 useIon ? ionZ[ti][psi] : 0,
673 useIon ? nElectrons[ti][psi] : 0,
674 useMat ? materialID[ti][psi] : -1,
675 (int)psi);
676 tpv.push_back(p);
677 ti = (int)pi;
678 }
679 std::reverse(tpv.begin(),tpv.end());
680 return tpv;
681}
682
683void BDSOutputROOTEventTrajectory::printTrajectoryInfoByTrackID(int trackIDIn) const
684{
685 // prevent a bad access
686 int storageIndex = 0;
687 auto search = trackID_trackIndex.find(trackIDIn);
688 if (search == trackID_trackIndex.end())
689 {
690 std::cout << "No such track ID" << std::endl;
691 return;
692 }
693 else
694 {storageIndex = search->second;}
695 printTrajectoryInfo(storageIndex);
696}
697
698void BDSOutputROOTEventTrajectory::printTrajectoryInfo(int storageIndex) const
699{
700 int i = storageIndex; // shortcut
701 int wdt = 11; // width of columns for print out
702
703 if (i > (int)partID.size() || i < 0)
704 {std::cout << "Invalid index" << std::endl; return;}
705
706 if (i+1 > n) // safety
707 {std::cout << "Index chosen is greater than maximum index of: " << n-1 << std::endl; return;}
708
709 // print out regarding the trajectory generally
710 std::cout << "Storage index " << std::setw(wdt) << i
711 << ", PDG ID " << std::setw(wdt) << partID[i]
712 << ", Track ID " << std::setw(wdt) << trackID[i]
713 << ", Parent ID " << std::setw(wdt) << parentID[i] << std::endl;
714 std::cout << "Created by Track ID " << parentID[i] << ", with storage index " << parentIndex[i] << ", and at step index " << parentStepIndex[i] << std::endl;
715
716 // print out regarding each step of the trajectory
717 std::cout << std::setw(wdt) << "step ind" << " "
718 << std::setw(wdt) << "prePrcT" << " " << std::setw(wdt) << "prePrcST" << " "
719 << std::setw(wdt) << "pstPrcT" << " " << std::setw(wdt) << "pstPrcST" << " "
720 << std::setw(wdt) << "X" << " " << std::setw(wdt) << "Y" << " "
721 << std::setw(wdt) << "Z" << " " << std::setw(wdt) << "E" << " "
722 << std::setw(wdt) << "p" << " " << std::setw(wdt) << "p_x" << " "
723 << std::setw(wdt) << "p_y" << " " << std::setw(wdt) << "p_z" << " "
724 << std::setw(wdt) << "t" << std::endl;
725
726 for (size_t j=0; j<XYZ[i].size(); ++j)
727 {
728 std::cout << std::setw(wdt) << j << " "
729 << std::setw(wdt) << preProcessTypes[i][j] << " " << std::setw(wdt) << preProcessSubTypes[i][j] << " "
730 << std::setw(wdt) << postProcessTypes[i][j] << " " << std::setw(wdt) << postProcessSubTypes[i][j] << " "
731 << std::setw(wdt) << XYZ[i][j].X() << " " << std::setw(wdt) << XYZ[i][j].Y() << " "
732 << std::setw(wdt) << XYZ[i][j].Z() << " " << std::setw(wdt) << energyDeposit[i][j] << " "
733 << std::setw(wdt) << PXPYPZ[i][j].Mag() << " " << std::setw(wdt) << PXPYPZ[i][j].X() << " "
734 << std::setw(wdt) << PXPYPZ[i][j].Y() << " " << std::setw(wdt) << PXPYPZ[i][j].Z() << " "
735 << std::setw(wdt) << T[i][j] << std::endl;
736 }
737}
738
739bool BDSOutputROOTEventTrajectory::parentIsPrimary(int trackIDIn) const
740{
741 // prevent a bad access
742 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
743 {
744 std::cout << "No such track ID" << std::endl;
745 return false;
746 }
747
748 unsigned int storageIndex = (unsigned int)trackID_trackIndex.at(trackIDIn);
749 unsigned int parentStorageIndex = parentIndex[storageIndex];
750 return parentID[parentStorageIndex] == 0;
751}
752
753std::ostream& operator<< (std::ostream& out, BDSOutputROOTEventTrajectory const &t)
754{
755 for (int i=0; i< (int)t.preProcessTypes.size();++i)
756 {
757 for (int j=0; j< (int)t.preProcessTypes[i].size(); ++j)
758 {
759 //if(t.preProcessTypes[i][j] != 1 && t.preProcessTypes[i][j] != 7)
760 //{
761 out << i << " " << j
762 << " " << t.preProcessTypes[i][j] << " " << t.preProcessSubTypes[i][j]
763 << " " << t.postProcessTypes[i][j] << " " << t.postProcessSubTypes[i][j]
764 << " " << t.preWeights[i][j] << " " << t.postWeights[i][j]
765 << " " << t.energyDeposit[i][j] << " " << t.T[i][j]
766 << std::endl;
767 //}
768 }
769 }
770 return out;
771}
Extra G4Navigator to get coordinate transforms.
G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=nullptr, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true, G4bool useCurvilinear=true) const
A wrapper for the underlying static navigator instance located within this class.
Structure to record a trajectory point.
Structure to record a trajectory.
std::vector< std::vector< int > > nElectrons
Ion trajectory information.
std::vector< std::vector< int > > charge
Link trajectory information.
std::vector< std::vector< TVector3 > > pxpypz
Local coordinates.
std::vector< std::vector< bool > > isIon
Ion trajectory information.
void Flush()
add comment to avoid warning (no need to make persistent, see issue #191)
void FillIndividualTrajectory(IndividualTrajectory &itj, BDSTrajectory *traj, int i, const std::map< G4Material *, short int > &materialToID) const
BDSAuxiliaryNavigator * auxNavigator
Required to find beamline index careful including in streamer.
std::vector< std::vector< double > > mass
Link trajectory information.
std::vector< std::vector< int > > ionZ
Ion trajectory information.
std::vector< std::vector< TVector3 > > xyz
Local coordinates.
std::vector< std::vector< int > > turnsTaken
Link trajectory information.
std::vector< std::vector< double > > rigidity
Link trajectory information.
std::vector< std::vector< int > > ionA
Ion trajectory information.
std::vector< std::vector< double > > kineticEnergy
Link trajectory information.
BDSPhysicalVolumeInfo * GetInfo(G4VPhysicalVolume *logicalVolume, G4bool isTunnel=false)
static BDSPhysicalVolumeInfoRegistry * Instance()
Singleton accessor.
A class holding any information pertaining to a particular physical volume in a BDSIM geant4 model.
G4int GetBeamlineIndex() const
Accessor.
Double map of trajectories to bitset of which filters matched whether to store them.
A Point in a trajectory with extra information.
G4double GetKineticEnergy() const
Accessor for the extra information links.
G4int GetCharge() const
Accessor for the extra information links.
G4double GetPreGlobalTime() const
Accessor.
G4int GetPostProcessType() const
Accessor.
G4ThreeVector GetMomentumLocal() const
Accessor.
G4int GetPreProcessSubType() const
Accessor.
G4bool GetIsIon() const
Accessor for the extra information ions.
G4double GetPreWeight() const
Accessor.
G4double GetRigidity() const
Accessor for the extra information links.
G4int GetTurnsTaken() const
Accessor for the extra information links.
G4double GetEnergyDeposit() const
Accessor.
G4double GetMass() const
Accessor for the extra information links.
G4double GetPreS() const
Accessor.
G4int GetPreProcessType() const
Accessor.
G4int GetIonZ() const
Accessor for the extra information ions.
G4ThreeVector GetPositionLocal() const
Accessor for the extra information local.
G4double GetPostWeight() const
Accessor.
G4int GetPostProcessSubType() const
Accessor.
G4int GetNElectrons() const
Accessor for the extra information ions.
G4int GetIonA() const
Accessor for the extra information ions.
G4Material * GetMaterial() const
Accessor.
G4ThreeVector GetPreMomentum() const
Accessor.
Trajectory information from track including last scatter etc.
virtual int GetPointEntries() const
Get number of trajectory points in this trajectory.
void SetTrajIndex(G4int trajIndexIn)
void SetParentStepIndex(G4int parentStepIndexIn)
The index of the step along the parent trajectory from which this one was created.
virtual G4VTrajectoryPoint * GetPoint(G4int i) const
Access a point - use this class's container.
void SetParentIndex(G4int parentIndexIn)
G4int GetDepth() const
Depth in the tree. Will be filled later once all trajectories are created and sorted.
Temporary structure for an individual trajectory used to convert types.