Add branch for StepPointMCs

doc/developers.md

@@ -62,8 +62,9 @@ mu2e -c EventNtuple/fcl/from_mcs-mockdata.fcl -s test-art-file.art -n 10
 
 8. Add to ```ntuplehelper```
    * instructions [here](./ntuplehelper.md#Adding-a-branch)
-9. Add to validation, RooUtil, and PyUtil
-   * [instructions to be completed... for the time being contact Andy and Sophie for this]
+9. Add to validation, RooUtil, and pyutils
+   * instructions for validation and RooUtil are [here](../rooutil/README.md#For-Developers). Contact Andy for more info if needed
+   * contact Sophie for pyutils
 10. Run validation
     * both [test_fcls.sh](../validation/README.md#Validating-EventNtuple-Runs), and
     * and [valCompare](../validation/README.md#Validating-EventNtuple-Contents) steps

fcl/README.md

@@ -22,6 +22,7 @@ where ```tier``` is the data tier of the input dataset, ```type``` is the type o
 | from_mcs-mockdata_separateTrkBranches.fcl | mock datasets | example on how to separate the tracks into separate branches again|
 | from_mcs-mockdata_selectorExample.fcl | mock datasets | example on how to use a selector to select certain types of tracks before putting them into the EventNtuple |
 | from_mcs-mixed_trkQualCompare.fcl | reconstructed mixed (i.e. primary+background hits) datasets | shows how to output result of more than one TrkQual |
+| from_mcs-primary_addVDSteps.fcl | reconstructed primary (i.e. no background hits) datasets | shows how to add the branch for virtual detector steps |
 | from_mcs-DeMCalib.fcl | reconstructed primary or mixed datasets | only writes one track per event |
 | from_mcs-OffSpill.fcl | off spill datasets | only contains ```CentralHelix``` tracks (i.e. field-on cosmics) |
 | from_dig-mockdata.fcl | mock datasets (digis) | runs reconstruction and creates EventNtuple in one job |

fcl/from_mcs-primary_addVDSteps.fcl

@@ -0,0 +1,3 @@
+#include "EventNtuple/fcl/from_mcs-primary.fcl"
+
+physics.analyzers.EventNtuple.StepPointMCTags : [ "compressRecoMCs:virtualdetector" ]

fcl/prolog.fcl

@@ -319,6 +319,7 @@ EventNtupleMaker : {
   ExtraMCStepCollectionTags : []
   SurfaceStepCollectionTag : "compressRecoMCs"
   FitType : LoopHelix
+  StepPointMCTags : [ ]
 }
 # instance for processing trigger (digitization) output from simulation
 EventNtupleMakerTTMC: {

helper/ntuplehelper.py

@@ -3,14 +3,15 @@
 class nthelper:
 
     single_object_branches = ['evtinfo', 'evtinfomc', 'hitcount', 'tcnt', 'crvsummary', 'crvsummarymc']
-    vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'trkqual', 'trkpid']
+    vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'trkqual', 'trkpid', 'mcsteps']
     vector_vector_object_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trkmcsim', 'trkhits', 'trkhitsmc', 'trkmats', 'trkhitcalibs', 'trkmcsci', 'trkmcssi', 'trksegsmc' ]
 
     evt_branches = ['evtinfo','evtinfomc','hitcount','tcnt']
     trk_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'trkqual', 'trkpid']
     trksegs_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trksegsmc']
     straw_branches = ['trkhits', 'trkmats', 'trkhitsmc', 'trkhitcalibs']
-    mc_branches = ['trkmcsim']
+    trk_mc_branches = [ 'trkmcsim' ]
+    general_mc_branches = [ 'mcsteps' ]
     calo_branches = ['caloclusters', 'calohits', 'calorecodigis', 'calodigis']
     crv_branches = ['crvsummary','crvsummarymc','crvcoincs','crvcoincsmc','crvcoincsmcplane']
     deprecated_branches = ['trkmcsci','trkmcssi']
@@ -58,7 +59,8 @@ class nthelper:
                            "trkqual" : "MVAResultInfo",
                            "trkpid" : "MVAResultInfo",
                            "helices" : "HelixInfo",
-                           "trksegsmc" : "SurfaceStepInfo"
+                           "trksegsmc" : "SurfaceStepInfo",
+                           "mcsteps" : "MCStepInfo"
                           }
 
     #
@@ -170,16 +172,34 @@ def list_all_branches(self, export_to_md=False):
                 print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")")
 
         if not export_to_md:
-            print("\nMonte Carlo Branches")
+            print("\nTrk Monte Carlo Branches")
             print("================")
         else:
-            print("## Monte Carlo Branches\n")
+            print("## Trk Monte Carlo Branches\n")
             print("These branches contain 4 elements per event corresponding to different Kalman fit hypotheses (see Track branches).\n")
             print("Within each Kalman fit element, there is a vector containing Monte Carlo truth information about the particle making the track and its parent particles.\n")
             print("The vector is sorted in reverse chronological order, such that the last element is the initial particle simulated in GEANT4, and each element before correspond to one of its daughter particles.\n")
             print("| branch | structure | explanation | leaf information |")
             print("|--------|-----------|-------------|------------------|")
-        for branch in self.mc_branches:
+        for branch in self.trk_mc_branches:
+            explanation = self.get_branch_explanation(branch)
+            struct = self.branch_struct_dict[branch]
+            struct_file = struct + ".hh";
+            if not export_to_md:
+                print(explanation)
+            else:
+                tokens=explanation.split(":")
+                print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")")
+
+        if not export_to_md:
+            print("\nGeneral Monte Carlo Branches")
+            print("================")
+        else:
+            print("## General Monte Carlo Branches\n")
+            print("These branches contain MC information in the event and are not required to be related to tracks")
+            print("| branch | structure | explanation | leaf information |")
+            print("|--------|-----------|-------------|------------------|")
+        for branch in self.general_mc_branches:
             explanation = self.get_branch_explanation(branch)
             struct = self.branch_struct_dict[branch]
             struct_file = struct + ".hh";

inc/InfoMCStructHelper.hh

@@ -75,7 +75,7 @@ namespace mu2e {
       void fillExtraMCStepInfos(KalSeedMC const& kseedmc, StepPointMCCollection const& mcsteps,
                                 std::vector<MCStepInfos>& mcsics, std::vector<MCStepSummaryInfo>& mcssis);
       void fillSurfaceStepInfos(KalSeedMC const& kseedmc, SurfaceStepCollection const& surfsteps,std::vector<SurfaceStepInfo>& ssic);
-
+      void fillStepPointMCInfo(StepPointMCCollection const& mcsteps, MCStepInfos& mcstepinfos);
   };
 }
 

inc/MCStepInfo.hh

@@ -14,12 +14,16 @@ namespace mu2e {
     float time = 0;  // time of this step WRT MC primary proton (ns)
     float de = 0; // energy deposit through this step (MeV)
     float dp = 0; // momentum magnitude change throw this step (MeV/c)
-    bool early = false;
-    bool late = false; // flag if this is the earliest or latest step
+    bool early = false; // flag if this is the earliest step
+    bool late = false; // flag if this is the latest step
     XYZVectorF mom; // particle momentum at the start of this step
     XYZVectorF pos;  // particle position at the start of this step
     void reset() {*this = MCStepInfo(); }
     bool valid() { return vid>=0; }
+    int pdg = -1; // true PDG code
+    int startCode = -1; // creation process code
+    int stopCode = -1; // stop process code
+
   };
 
   using MCStepInfos = std::vector<MCStepInfo>;

rooutil/README.md

@@ -114,6 +114,7 @@ Some branches are not contained in any of the above classes:
 * ```crvcoincsmcplane```
 * ```calohits```, ```calorecodigis```, ```calodigis```
 * ```trig_``` branches
+* ```mcsteps_virtualdetector```
 
 Reach out to the developers on the #analysis-tools Slack channel if you need to have these added somewhere.
 
@@ -224,23 +225,22 @@ Beware: there are a lot of debug messages.
 If a new branch is added to the EventNtuple, then the following needs to be done to so that RooUtil can access the branch:
 
 1. In [Event.hh](inc/Event.hh) add the pointers at the bottom of the file
-2. In [Event.hh](inc/Event.hh) constructor, set the branch address
-  - make sure to test that the branch exists
-3. In [Event.hh](inc/Event.hh) add the #include to the underlying object
+2. In [Event.hh](inc/Event.hh) constructor, set the branch address using the ```CheckForBranch()``` function
+3. In [Event.hh](inc/Event.hh) add the #include to the underlying info struct
 4. In [RooUtil.hh](inc/RooUtil.hh) add it to the ```CreateOutputEventNtuple()``` function
 5. Add to validation places:
   - [PrintEvents.C](examples/PrintEvents.C): at least the first and last leaf in the struct
   - [create_val_file_rooutil.C](../../validation/create_val_file_rooutil.C)
      - copy contents of struct into place where histograms are defined
      - then copy and replace "type " with "TH1F* h_branchname_" (e.g. "float " with "TH1F* h_trkcalohit_")
-     - then copy in " = new TH1F("h_branchname_", "", 100,0,100); //"
+     - then copy in " = new TH1F("h_branchname_", "", 100,0,100); //" after the histogram name
      - delete line after "//"
      - add in leaf names to histname (and make separate x, y, z histograms for XYZVectorF leaves)
      - copy histogram lines into main loop
-     - search and replace "TH1F* " with ""
-     - search and replace " = new TH1F("h_" with "->Fill("
-     - search and replace "", "", 100,0,100);" with ");"
-     - search and replace "(branchname_" with "(branchname."
+     - search and replace ```TH1F* ``` (note space) with nothing
+     - search and replace ``` = new TH1F("h_``` with ```->Fill(```
+     - search and replace ```", "", 100,0,100);``` with ```);```
+     - search and replace ```(branchname_``` with ```(branchname.``` in the ```Fill()``` commands
      - update histogram ranges
 6. (Optional) If appropriate, add branches to other classes (e.g. Track.hh) and to ```Event::Update()```
    - be sure to test it with an example script

rooutil/examples/PlotVDSteps.C

@@ -0,0 +1,52 @@
+//
+// An example of how to plot reco vs true
+// This uses cut functions defined in common_cuts.hh
+//
+
+#include "EventNtuple/rooutil/inc/RooUtil.hh"
+#include "EventNtuple/rooutil/inc/common_cuts.hh"
+
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+
+using namespace rooutil;
+void PlotVDSteps(std::string filename) {
+
+  // Create the histogram you want to fill
+  TH2F* hVD13_XY = new TH2F("hVD13_XY", "XY position at VD13 (tracker entrance)", 100,-500,500, 100,-500,500);
+  hVD13_XY->SetXTitle("X (det. coords) [mm]");
+  hVD13_XY->SetYTitle("Y (det. coords) [mm]");
+
+  TH1F* hVD13_Mom = new TH1F("hVD13_Mom", "Momentum at VD13 (tracker entrance)", 120,0,120);
+  hVD13_Mom->SetXTitle("Momentum [MeV/c]");
+
+  // Set up RooUtil
+  RooUtil util(filename);
+
+  // Loop through the events
+  for (int i_event = 0; i_event < util.GetNEvents(); ++i_event) {
+    // Get the next event
+    auto& event = util.GetEvent(i_event);
+
+    // Get the e_minus tracks from the event
+    if (event.mcsteps_virtualdetector != nullptr) {
+      for (const auto& vdstep : *(event.mcsteps_virtualdetector)) {
+        if (vdstep.vid == mu2e::VirtualDetectorId::TT_FrontHollow) {
+          // Fill the histogram
+          hVD13_XY->Fill(vdstep.pos.x(), vdstep.pos.y());
+
+          hVD13_Mom->Fill(vdstep.mom.R());
+        }
+      }
+    }
+  }
+
+  // Draw the histogram
+  TCanvas* c = new TCanvas();
+  c->Divide(2);
+  c->cd(1);
+  hVD13_XY->Draw("COLZ");
+  c->cd(2);
+  hVD13_Mom->Draw("HIST E");
+}

rooutil/examples/PrintEvents.C

@@ -263,5 +263,11 @@ void PrintEvents(std::string filename) {
       std::cout << pair.first << ": " << event.trigger.Fired(pair.second) << std::endl;
     }
 
+    // mcsteps_virtualdetector branch
+    if (event.mcsteps_virtualdetector != nullptr) {
+      for (const auto& vdstep : *(event.mcsteps_virtualdetector)) {
+        std::cout << "vdstep: " << vdstep.vid << ", " << vdstep.time << ", " << vdstep.mom << ", " << vdstep.pos << ", " << vdstep.pdg << ", " << vdstep.startCode << ", " << vdstep.stopCode << std::endl;
+      }
+    }
   }
 }

rooutil/inc/Event.hh

@@ -37,6 +37,7 @@
 #include "EventNtuple/inc/MVAResultInfo.hh"
 
 #include "EventNtuple/inc/TrigInfo.hh"
+#include "EventNtuple/inc/MCStepInfo.hh"
 
 #include "EventNtuple/rooutil/inc/Track.hh"
 #include "EventNtuple/rooutil/inc/CrvCoinc.hh"
@@ -87,6 +88,7 @@ namespace rooutil {
       CheckForBranch(ntuple, "calorecodigis", &this->calorecodigis);
       CheckForBranch(ntuple, "calodigis", &this->calodigis);
 
+      CheckForBranch(ntuple, "mcsteps_virtualdetector", &this->mcsteps_virtualdetector);
     }
 
     // Check if a branch exists in the TChain, and optionally set its address
@@ -342,6 +344,7 @@ namespace rooutil {
     std::vector<mu2e::CrvPlaneInfoMC>* crvcoincsmcplane = nullptr;
 
     std::vector<std::vector<mu2e::SimInfo>>* trkmcsim = nullptr;
+    std::vector<mu2e::MCStepInfo>* mcsteps_virtualdetector = nullptr; // TODO: EventNtuple could have other mcsteps_* branches but for the time being just hardcode for the virtualdetector ones
   };
 } // namespace rooutil
 #endif

rooutil/inc/RooUtil.hh

@@ -148,6 +148,8 @@ namespace rooutil {
         output_ntuple->Branch(("trig_" + pair.first).c_str(), &event->triginfo._triggerArray[pair.second]);
       }
 
+      if (event->mcsteps_virtualdetector) { output_ntuple->Branch("mcsteps_virtualdetector", event->mcsteps_virtualdetector); }
+
       // Write out histograms from input to output
       hVersion->Write();
     }

rooutil/inc/common_cuts.hh

@@ -10,6 +10,7 @@
 #include "Offline/MCDataProducts/inc/ProcessCode.hh"
 #include "Offline/MCDataProducts/inc/GenId.hh"
 #include "Offline/MCDataProducts/inc/MCRelationship.hh"
+#include "Offline/DataProducts/inc/VirtualDetectorId.hh"
 
 #include "EventNtuple/rooutil/inc/Track.hh"
 #include "EventNtuple/rooutil/inc/TrackSegment.hh"

src/EventNtupleMaker_module.cc

@@ -182,9 +182,10 @@ namespace mu2e {
         fhicl::Atom<art::InputTag> primaryParticleTag{Name("PrimaryParticleTag"), Comment("Tag for PrimaryParticle"), art::InputTag()};
         fhicl::Atom<art::InputTag> kalSeedMCTag{Name("KalSeedMCAssns"), Comment("Tag for KalSeedMCAssn"), art::InputTag()};
         // extra MC
-        fhicl::OptionalSequence<art::InputTag> extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out")};
+        fhicl::OptionalSequence<art::InputTag> extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out. This will only write out steps associated with the KalSeed")};
         // passive elements and Virtual Detector MC information
         fhicl::OptionalAtom<art::InputTag> SurfaceStepsTag{Name("SurfaceStepCollectionTag"), Comment("SurfaceStep Collection")};
+        fhicl::OptionalSequence<art::InputTag> stepPointMCTags{Name("StepPointMCTags"), Comment("Input tags for any other StepPointMCCollections you want written out. This will write out all steps in the collection")};
         // Calo MC
         fhicl::Atom<bool> fillCaloMC{ Name("FillCaloMC"),Comment("Fill CaloMC information"), true};
         fhicl::Atom<art::InputTag> caloClusterMCTag{Name("CaloClusterMCTag"), Comment("Tag for CaloClusterMCCollection") ,art::InputTag()};
@@ -249,7 +250,7 @@ namespace mu2e {
       std::map<size_t,unsigned> _tmap; // map between path and trigger ID.  ID should come from trigger itself FIXME!
       // MC truth (fcl parameters)
       bool _fillmc;
-      // MC truth inputs
+      // MC steps (associated with KalSeed)
       std::vector<art::InputTag> _extraMCStepTags;
       std::vector<art::Handle<StepPointMCCollection>> _extraMCStepCollections;
       std::map<BranchIndex, std::map<StepCollIndex, std::vector<MCStepInfos>>> _extraMCStepInfos;
@@ -258,6 +259,11 @@ namespace mu2e {
       art::InputTag _surfaceStepsTag;
       std::map<BranchIndex, std::vector<std::vector<SurfaceStepInfo>>> _surfaceStepInfos;
       art::Handle<SurfaceStepCollection> _surfaceStepsHandle;
+      // MC steps (all)
+      std::vector<art::InputTag> _stepPointMCTags;
+      std::vector<art::Handle<StepPointMCCollection>> _stepPointMCCollections;
+      std::map<StepCollIndex, MCStepInfos> _stepPointMCInfos;
+
       //
       art::Handle<PrimaryParticle> _pph;
       art::Handle<KalSeedMCAssns> _ksmcah;
@@ -441,6 +447,11 @@ namespace mu2e {
           _surfaceStepInfos[i_branch] = std::vector<std::vector<SurfaceStepInfo>>();
         }
       }
+      if(_conf.stepPointMCTags(_stepPointMCTags)){
+        for (StepCollIndex i_stepPointMCTag = 0; i_stepPointMCTag < _stepPointMCTags.size(); ++i_stepPointMCTag) {
+          _stepPointMCInfos[i_stepPointMCTag] = MCStepInfos();
+        }
+      }
     }
   }
 
@@ -566,6 +577,18 @@ namespace mu2e {
     }
     // helix info
     if(_conf.helices()) _ntuple->Branch("helices.",&_hinfos,_buffsize,_splitlevel);
+
+    // all MC information
+    if (_fillmc) {
+      if(_conf.stepPointMCTags(_stepPointMCTags)){
+        for(size_t icoll=0;icoll < _stepPointMCTags.size(); ++icoll) {
+          auto const& tag = _stepPointMCTags[icoll];
+          auto inst = tag.instance();
+          std::string branchname  = "mcsteps_" + inst + ".";
+          _ntuple->Branch(branchname.c_str(),&_stepPointMCInfos[icoll],_buffsize,_splitlevel);
+        }
+      }
+    }
   }
 
   void EventNtupleMaker::beginSubRun(const art::SubRun & subrun ) {
@@ -675,6 +698,16 @@ namespace mu2e {
     }
     event.getByLabel(_surfaceStepsTag,_surfaceStepsHandle);
 
+    // find extra MCStep collections
+    _stepPointMCCollections.clear();
+    for(size_t icoll = 0; icoll < _stepPointMCTags.size(); icoll++){
+      auto const& tag = _stepPointMCTags[icoll];
+      art::Handle<StepPointMCCollection> mcstepch;
+      event.getByLabel(tag,mcstepch);
+      _stepPointMCCollections.push_back(mcstepch);
+    }
+
+
     // loop through all track types
     unsigned ntrks(0);
     for (BranchIndex i_branch = 0; i_branch < _allBranches.size(); ++i_branch) {
@@ -709,6 +742,10 @@ namespace mu2e {
         _extraMCStepSummaryInfos.at(i_branch).at(i_extraMCStepTag).clear();
       }
       _surfaceStepInfos.at(i_branch).clear();
+      for (StepCollIndex i_stepPointMCTag = 0; i_stepPointMCTag < _stepPointMCTags.size(); ++i_stepPointMCTag) {
+        _stepPointMCInfos.at(i_stepPointMCTag).clear();
+      }
+
 
       if(_fillcalomc) { _allMCTCHIs.at(i_branch).clear(); }
 
@@ -875,6 +912,17 @@ namespace mu2e {
       }
 
     }
+
+    // fill MC information unrelated to any reconstructed object
+    if (_fillmc) {
+      // fill MCStepCollection branch
+      for(size_t icoll = 0; icoll < _stepPointMCTags.size(); icoll++){
+        auto const& mcsteps = *_stepPointMCCollections[icoll];
+        auto& mcstepinfos = _stepPointMCInfos[icoll];
+        _infoMCStructHelper.fillStepPointMCInfo(mcsteps, mcstepinfos);
+      }
+    }
+
     // fill this row in the TTree
     bool fill = true; // default to fliling event
     if(_hastrks && ntrks == 0) { // if we require tracks (_hastrks) but we have none, don't write