New feature: Keep record of where a track originates

This allows to see which volumes have a (source) influence on which other volumes.

Author: sawenzel

include/MCStepLogger/SimpleStepAnalysis.h

@@ -58,6 +58,11 @@ class SimpleStepAnalysis : public MCAnalysis
   TH1I* histNStepsPerPDG;
   TH1I* histNStepsPerVolSorted;
 
+  // track production origins per module and volume
+  TH1I* histOriginPerMod;
+  TH1I* histOriginPerVol;
+  TH1I* histOriginPerVolSorted;
+
   // accumulated number of secondaries produced per volume
   TH1I* histNSecondariesPerVol;
   // accumulated number of secondaries produces per module

include/MCStepLogger/StepInfo.h

@@ -60,7 +60,9 @@ struct StepLookups {
   std::vector<int> tracktopdg;
   std::vector<int> tracktoparent; // when parent is -1 we mean primary
   std::vector<int> stepcounterpertrack;
-  std::vector<float> tracktoenergy; // mapping of trackID to the track energy
+
+  std::vector<int> trackorigin;      // the volumeID where the track originates
+  std::vector<float> tracktoenergy;  // mapping of trackID to the track energy
   std::vector<bool> crossedboundary; // if track every crossed a geometry boundary
   std::vector<bool> producedsecondary; // if track ever produced a secondary
 
@@ -117,6 +119,14 @@ struct StepLookups {
     }
   }
 
+  void setTrackOrigin(int trackindex, int volID)
+  {
+    if (trackindex >= trackorigin.size()) {
+      trackorigin.resize(trackindex + 1, -1);
+    }
+    trackorigin[trackindex] = volID;
+  }
+
   bool initSensitiveVolLookup(std::string const& filename);
 
   void insertParent(int trackindex, int parent)

src/SimpleStepAnalysis.cxx

@@ -34,6 +34,11 @@ void SimpleStepAnalysis::initialize()
   histNStepsPerMod = getHistogram<TH1I>("nStepsPerMod", 1, 2., 1.);
   // accumulated number of steps per volume
   histNStepsPerVol = getHistogram<TH1I>("nStepsPerVol", 1, 2., 1.);
+
+  histOriginPerMod = getHistogram<TH1I>("OriginsPerMod", 1, 2., 1.);
+  histOriginPerVol = getHistogram<TH1I>("OriginsPerVol", 1, 2., 1.);
+  histOriginPerVolSorted = getHistogram<TH1I>("OriginsPerVolSorted", 1, 2., 1.);
+
   // accumulated number of steps per pdg
   histNStepsPerPDG = getHistogram<TH1I>("nStepsPerPDG", 1, 2., 1.);
   histNStepsPerVolSorted = getHistogram<TH1I>("nStepsPerVolSorted", 1, 2., 1.);
@@ -115,9 +120,7 @@ void SimpleStepAnalysis::analyze(const std::vector<StepInfo>* const steps, const
 
   // loop over all steps in an event
   for (const auto& step : *steps) {
-    int currentTrackID = step.trackID;
-    bool newtrack = (currentTrackID != oldTrackID);
-    if (newtrack) oldTrackID = currentTrackID;
+
 
     // prepare for PDG ids and volume names
     mAnalysisManager->getLookupPDG(step.trackID, pdgId);
@@ -133,6 +136,12 @@ void SimpleStepAnalysis::analyze(const std::vector<StepInfo>* const steps, const
       continue;
     }
 
+    int currentTrackID = step.trackID;
+    bool newtrack = (currentTrackID != oldTrackID);
+    if (newtrack) {
+      oldTrackID = currentTrackID;
+    }
+
     if (keepsteps) {
       stepptr = &step;
       steptree->Fill();
@@ -145,6 +154,16 @@ void SimpleStepAnalysis::analyze(const std::vector<StepInfo>* const steps, const
       histTrackEnergySpectrum->Fill(log10f(step.E));
       histTrackPDGSpectrum->Fill(pdgasstring.c_str(),1);
       histTrackProdProcess->Fill(step.getProdProcessAsString(),1);
+
+      auto originid = mAnalysisManager->getLookups()->trackorigin[step.trackID];
+      std::string originVolName;
+      // to store the module name
+      std::string originModName;
+      mAnalysisManager->getLookupVolName(originid, originVolName);
+      mAnalysisManager->getLookupModName(originid, originModName);
+
+      histOriginPerMod->Fill(originModName.c_str(), 1);
+      histOriginPerVol->Fill(originVolName.c_str(), 1);
     }
 
     // record number of steps per module
@@ -169,6 +188,11 @@ void SimpleStepAnalysis::finalize()
   histNStepsPerVolSorted->SetName("nStepsPerVolSorted");
   histNStepsPerVolSorted->LabelsOption(">", "X");
 
+  *histOriginPerVolSorted = *histOriginPerVol;
+  histOriginPerVolSorted->SetName("OriginPerVolSorted");
+  histOriginPerVolSorted->LabelsOption(">", "X");
+  histOriginPerVolSorted->SetBins(30, 0, 30);
+
   std::cerr << "MOD have " << histNStepsPerMod->GetEntries() << " entries \n";
 
   *histTrackPDGSpectrumSorted = *histTrackPDGSpectrum;

src/StepInfo.cxx

@@ -118,7 +118,11 @@ StepInfo::StepInfo(TVirtualMC* mc)
   if (mc->IsTrackExiting()) {
     lookupstructures.setCrossedBoundary(trackID, true);
   }
-  
+
+  if (mc->IsNewTrack()) {
+    lookupstructures.setTrackOrigin(trackID, volId);
+  }
+
   prodprocess = stack->GetCurrentTrack()->GetUniqueID();
 
   TArrayI procs;