24#include "BDSOutputROOTEventTrajectory.hh"
27#include "G4VPhysicalVolume.hh"
29#include "BDSHitEnergyDeposition.hh"
30#include "BDSAuxiliaryNavigator.hh"
31#include "BDSPhysicalVolumeInfoRegistry.hh"
32#include "BDSPhysicalVolumeInfo.hh"
33#include "BDSTrajectory.hh"
34#include "BDSTrajectoryOptions.hh"
41BDSOutputROOTEventTrajectory::BDSOutputROOTEventTrajectory():
42 auxNavigator(nullptr),
46BDSOutputROOTEventTrajectory::~BDSOutputROOTEventTrajectory()
56 if (!traj->GetParent())
59 if (traj->GetParent()->GetTrackID() == 1)
60 {
return traj->GetParentStepIndex();}
62 {
return findPrimaryStepIndex(traj->GetParent());}
67 bool storeStepPointLast,
69 const std::map<G4Material*, short int>& materialToID)
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;
84 for (
auto trajFlag : trajectories->trajectories)
97 for (
auto trajFlag : trajectories->trajectories)
101 if (trajFlag.second && parent)
106 if (parent->GetTrajIndex() != -1)
108 auto trajStartPos = traj->
GetPoint(0)->GetPosition();
112 if(parent->
GetPoint(i)->GetPosition() == trajStartPos)
127 for (
auto trajFlag : trajectories->trajectories)
132 if (!trajFlag.second)
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());
144 IndividualTrajectory itj;
145 if (storeStepPointsN > 0)
148 G4int nPoints = std::min(nSteps, storeStepPointsN);
149 for (
int i = 0; i < nPoints; ++i)
152 if (storeStepPointLast && (nPoints < nSteps))
162 filters.push_back(trajectories->filtersMatched.at(traj));
164 XYZ.push_back(itj.XYZ);
165 modelIndicies.push_back(itj.modelIndex);
168 {PXPYPZ.push_back(itj.PXPYPZ);}
174 preProcessTypes.push_back(itj.preProcessType);
175 preProcessSubTypes.push_back(itj.preProcessSubType);
176 postProcessTypes.push_back(itj.postProcessType);
177 postProcessSubTypes.push_back(itj.postProcessSubType);
180 preWeights.push_back(itj.preWeight);
181 postWeights.push_back(itj.postWeight);
182 energyDeposit.push_back(itj.energyDeposit);
185 {T.push_back(itj.T);}
191 {materialID.push_back(itj.materialID);}
193 if (!itj.xyz.empty())
195 xyz.push_back(itj.xyz);
196 pxpypz.push_back(itj.pxpypz);
199 if (!itj.charge.empty())
201 charge.push_back(itj.charge);
203 mass.push_back(itj.mass);
207 if (!itj.isIon.empty())
209 isIon.push_back(itj.isIon);
210 ionA.push_back(itj.ionA);
211 ionZ.push_back(itj.ionZ);
216 primaryStepIndex.push_back(findPrimaryStepIndex(traj));
219 trackID_trackIndex.insert(std::pair<int,int>(traj->GetTrackID(),n));
227 for (
auto iT = trajVec.begin(); iT != trajVec.end(); ++iT)
232 trackID_trackIndex.insert(std::pair<int, int>(traj->GetTrackID(),trackIndex));
234 std::cout << trajVec.size() <<
" " << parentID.size() <<
" " << parentIndex.size() <<
" "
235 << traj->GetTrackID() <<
" " << traj->GetParentID() <<
" " << trackIndex << std::endl;
238 auto processPair = findParentProcess(trackIndex);
239 trackIndex_trackProcess.insert(std::pair<
int,std::pair<int,int>>(trackIndex,processPair));
242 if(processPair.first != -1)
244 int mi = modelIndicies[processPair.first][processPair.second];
245 trackIndex_modelIndex.insert(std::pair<int,int>(trackIndex, mi));
247 {modelIndex_trackIndex.at(mi).push_back(trackIndex);}
248 catch(
const std::exception&)
250 modelIndex_trackIndex.insert(std::pair<
int,std::vector<int>>(mi,std::vector<int>()));
251 modelIndex_trackIndex.at(mi).push_back(trackIndex);
263 const std::map<G4Material*, short int>& materialToID)
const
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));
280 {itj.modelIndex.push_back(-1);}
292 itj.PXPYPZ.emplace_back(TVector3(mom.getX(),
295 itj.S.push_back(point->
GetPreS() / CLHEP::m);
299 itj.materialID.push_back(materialToID.at(point->
GetMaterial()));
301 if (point->extraLocal)
305 itj.xyz.emplace_back(TVector3(localPos.getX(),
308 itj.pxpypz.emplace_back(TVector3(localMom.getX(),
313 if (point->extraLink)
315 itj.charge.push_back((
int) (point->
GetCharge() / (G4double)CLHEP::eplus));
317 itj.mass.push_back(point->
GetMass() / CLHEP::GeV);
318 itj.rigidity.push_back(point->
GetRigidity() / (CLHEP::tesla*CLHEP::m));
323 itj.isIon.push_back(point->
GetIsIon());
324 itj.ionA.push_back(point->
GetIonA());
325 itj.ionZ.push_back(point->
GetIonZ());
332 G4cout << phc->GetSize() << G4endl;
345 parentStepIndex.clear();
346 primaryStepIndex.clear();
348 preProcessTypes.clear();
349 preProcessSubTypes.clear();
350 postProcessTypes.clear();
351 postProcessSubTypes.clear();
354 energyDeposit.clear();
371 modelIndicies.clear();
372 trackID_trackIndex.clear();
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;
402 PXPYPZ = other->PXPYPZ;
403 modelIndicies = other->modelIndicies;
404 trackID_trackIndex = other->trackID_trackIndex;
417 materialID = other->materialID;
421std::pair<int,int> BDSOutputROOTEventTrajectory::findParentProcess(
int trackIndex)
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;
432 {
return std::pair<int,int>(-1,-1);}
433 int sin = parentStepIndex.at(tid);
434 std::cout << sin << std::endl;
441 pid = parentID.at(tid);
442 pin = parentIndex.at(tid);
443 sin = parentStepIndex.at(tid);
445 std::cout << tid <<
" " << pid <<
" " << pin <<
" " << sin <<
" " << std::endl;
448 return std::pair<int,int>(pin,sin);
452std::vector<BDSOutputROOTEventTrajectoryPoint> BDSOutputROOTEventTrajectory::trackInteractions(
int trackIDIn)
const
455 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
457 std::cout <<
"No such track ID" << std::endl;
458 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
461 int ti = trackID_trackIndex.at(trackIDIn);
463 std::vector<BDSOutputROOTEventTrajectoryPoint> tpv;
465 int nstep = (int)XYZ[ti].size();
467 {
return std::vector<BDSOutputROOTEventTrajectoryPoint>();}
469 if (postProcessTypes[ti].empty())
471 std::cout <<
"Processes not stored for this file - not possible." << std::endl;
472 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
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();
481 bool useMat = !materialID.empty();
483 for (
int i = 0; i < nstep; ++i)
485 int ppt = postProcessTypes[ti][i];
486 int ppst = postProcessSubTypes[ti][i];
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)
500 postProcessTypes[ti][i],
501 postProcessSubTypes[ti][i],
503 energyDeposit[ti][i],
505 usePXPYPZ ? PXPYPZ[ti][i] : TVector3(),
506 modelIndicies[ti][i],
508 useLocal ?
xyz[ti][i] : TVector3(),
509 useLocal ?
pxpypz[ti][i] : TVector3(),
510 useLinks ?
charge[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,
519 useMat ? materialID[ti][i] : -1,
530 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
532 std::cout <<
"No such track ID" << std::endl;
535 int ti = trackID_trackIndex.at(trackIDIn);
536 int pid = (int)parentID[ti];
537 int chosenTrackID = trackIDIn;
540 if (parentID[trackID_trackIndex.at(pid)] > 0)
541 {chosenTrackID = pid;}
542 ti = trackID_trackIndex.at(pid);
543 pid = (int)parentID[ti];
545 return parentProcessPoint(chosenTrackID);
551 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
553 std::cout <<
"No such track ID" << std::endl;
557 int ti = trackID_trackIndex.at(trackIDIn);
558 int pti = (int)parentID[ti];
562 std::cout <<
"Track is a parent" << std::endl;
566 int si = (int)parentStepIndex.at(ti);
569 {pi = trackID_trackIndex.at(pti);}
573 if (si > (
int)XYZ[pi].size())
575 std::cout <<
"Not all step points are stored. Parent step index is outside points stored." << std::endl;
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();
585 bool useMat = !materialID.empty();
591 postProcessTypes[pi][si],
592 postProcessSubTypes[pi][si],
594 energyDeposit[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,
605 useLinks ?
mass[pi][si] : 0,
606 useIon ?
isIon[pi][si] :
false,
607 useIon ?
ionA[pi][si] : 0,
608 useIon ?
ionZ[pi][si] : 0,
610 useMat ? materialID[pi][si] : -1,
615std::vector<BDSOutputROOTEventTrajectoryPoint> BDSOutputROOTEventTrajectory::processHistory(
int trackIDIn)
const
618 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
620 std::cout <<
"No such track ID" << std::endl;
621 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
624 if (postProcessTypes.empty())
626 std::cout <<
"Processes not stored for this file - not possible." << std::endl;
627 return std::vector<BDSOutputROOTEventTrajectoryPoint>();
630 int ti = trackID_trackIndex.at(trackIDIn);
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();
638 bool useMat = !materialID.empty();
640 std::vector<BDSOutputROOTEventTrajectoryPoint> tpv;
643 unsigned int pi = parentIndex.at(ti);
644 unsigned int psi = parentStepIndex.at(ti);
645 if (psi > (
unsigned int)XYZ[pi].size())
647 std::cout <<
"Not all points stored - defaulting to creation point." << std::endl;
655 postProcessTypes[pi][psi],
656 postProcessSubTypes[pi][psi],
657 postWeights[pi][psi],
658 energyDeposit[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,
669 useLinks ?
mass[ti][psi] : 0,
670 useIon ?
isIon[ti][psi] :
false,
671 useIon ?
ionA[ti][psi] : 0,
672 useIon ?
ionZ[ti][psi] : 0,
674 useMat ? materialID[ti][psi] : -1,
679 std::reverse(tpv.begin(),tpv.end());
683void BDSOutputROOTEventTrajectory::printTrajectoryInfoByTrackID(
int trackIDIn)
const
686 int storageIndex = 0;
687 auto search = trackID_trackIndex.find(trackIDIn);
688 if (search == trackID_trackIndex.end())
690 std::cout <<
"No such track ID" << std::endl;
694 {storageIndex = search->second;}
695 printTrajectoryInfo(storageIndex);
698void BDSOutputROOTEventTrajectory::printTrajectoryInfo(
int storageIndex)
const
700 int i = storageIndex;
703 if (i > (
int)partID.size() || i < 0)
704 {std::cout <<
"Invalid index" << std::endl;
return;}
707 {std::cout <<
"Index chosen is greater than maximum index of: " << n-1 << std::endl;
return;}
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;
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;
726 for (
size_t j=0; j<XYZ[i].size(); ++j)
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;
739bool BDSOutputROOTEventTrajectory::parentIsPrimary(
int trackIDIn)
const
742 if (trackID_trackIndex.find(trackIDIn) == trackID_trackIndex.end())
744 std::cout <<
"No such track ID" << std::endl;
748 unsigned int storageIndex = (
unsigned int)trackID_trackIndex.at(trackIDIn);
749 unsigned int parentStorageIndex = parentIndex[storageIndex];
750 return parentID[parentStorageIndex] == 0;
755 for (
int i=0; i< (int)t.preProcessTypes.size();++i)
757 for (
int j=0; j< (int)t.preProcessTypes[i].size(); ++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]
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.