Add mcsteps_virtualdetector branch to RooUtil

Author: AndrewEdmonds11

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/PrintEvents.C

@@ -263,5 +263,10 @@ 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();
     }

validation/create_val_file_rooutil.C

@@ -505,6 +505,25 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) {
   TH1F* h_calodigis_caloRecoDigiIdx_ = new TH1F("h_calodigis_caloRecoDigiIdx_", "", 100,0,100);
 
   TH1F* h_triginfo = new TH1F("h_triginfo", "", mu2e::TrigInfo::ntrig_,0,mu2e::TrigInfo::ntrig_);
+
+  TH1F* h_mcsteps_virtualdetector_vid = new TH1F("h_mcsteps_virtualdetector_vid", "", 150,0,150);
+  TH1F* h_mcsteps_virtualdetector_sid = new TH1F("h_mcsteps_virtualdetector_sid", "", 100,0,100);
+  TH1F* h_mcsteps_virtualdetector_iinter = new TH1F("h_mcsteps_virtualdetector_iinter", "", 100,0,100);
+  TH1F* h_mcsteps_virtualdetector_time = new TH1F("h_mcsteps_virtualdetector_time", "", 100,0,2000);
+  TH1F* h_mcsteps_virtualdetector_de = new TH1F("h_mcsteps_virtualdetector_de", "", 100,0,200);
+  TH1F* h_mcsteps_virtualdetector_dp = new TH1F("h_mcsteps_virtualdetector_dp", "", 100,0,200);
+  TH1F* h_mcsteps_virtualdetector_early = new TH1F("h_mcsteps_virtualdetector_early", "", 2,0,2);
+  TH1F* h_mcsteps_virtualdetector_late = new TH1F("h_mcsteps_virtualdetector_late", "", 2,0,2);
+  TH1F* h_mcsteps_virtualdetector_mom_x = new TH1F("h_mcsteps_virtualdetector_mom_x", "", 100,-200,200);
+  TH1F* h_mcsteps_virtualdetector_mom_y = new TH1F("h_mcsteps_virtualdetector_mom_y", "", 100,-200,200);
+  TH1F* h_mcsteps_virtualdetector_mom_z = new TH1F("h_mcsteps_virtualdetector_mom_z", "", 320,-200,200);
+  TH1F* h_mcsteps_virtualdetector_pos_x = new TH1F("h_mcsteps_virtualdetector_pos_x", "", 100,-1000,1000);
+  TH1F* h_mcsteps_virtualdetector_pos_y = new TH1F("h_mcsteps_virtualdetector_pos_y", "", 100,-1000,1000);
+  TH1F* h_mcsteps_virtualdetector_pos_z = new TH1F("h_mcsteps_virtualdetector_pos_z", "", 320,-1600,1600);
+  TH1F* h_mcsteps_virtualdetector_pdg = new TH1F("h_mcsteps_virtualdetector_pdg", "", 100,0,100);
+  TH1F* h_mcsteps_virtualdetector_startCode = new TH1F("h_mcsteps_virtualdetector_startCode", "", 100,0,100);
+  TH1F* h_mcsteps_virtualdetector_stopCode = new TH1F("h_mcsteps_virtualdetector_stopCode", "", 100,0,100);
+
   // Grab the first event to name the bins
   const auto& evt = util.GetEvent(0);
   for (const auto pair : evt.trigger.NameToIndexMap()) {
@@ -1153,6 +1172,30 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) {
       h_triginfo->AddBinContent(i_trig+1, event.triginfo._triggerArray[i_trig]);
     }
 
+    // Creating mcsteps_virtualdetector histogram
+    std::cout << "Creating mcsteps_virtualdetector histogram" << std::endl;
+    if (event.mcsteps_virtualdetector != nullptr) {
+      for (const auto& vdstep : *(event.mcsteps_virtualdetector)) {
+        h_mcsteps_virtualdetector_vid->Fill(mcsteps_virtualdetector.vid);
+        h_mcsteps_virtualdetector_sid->Fill(mcsteps_virtualdetector.sid);
+        h_mcsteps_virtualdetector_iinter->Fill(mcsteps_virtualdetector.iinter);
+        h_mcsteps_virtualdetector_time->Fill(mcsteps_virtualdetector.time);
+        h_mcsteps_virtualdetector_de->Fill(mcsteps_virtualdetector.de);
+        h_mcsteps_virtualdetector_dp->Fill(mcsteps_virtualdetector.dp);
+        h_mcsteps_virtualdetector_early->Fill(mcsteps_virtualdetector.early);
+        h_mcsteps_virtualdetector_late->Fill(mcsteps_virtualdetector.late);
+        h_mcsteps_virtualdetector_mom_x->Fill(mcsteps_virtualdetector.mom.x());
+        h_mcsteps_virtualdetector_mom_y->Fill(mcsteps_virtualdetector.mom.y());
+        h_mcsteps_virtualdetector_mom_z->Fill(mcsteps_virtualdetector.mom.z());
+        h_mcsteps_virtualdetector_pos_x->Fill(mcsteps_virtualdetector.pos.x());
+        h_mcsteps_virtualdetector_pos_y->Fill(mcsteps_virtualdetector.pos.y());
+        h_mcsteps_virtualdetector_pos_z->Fill(mcsteps_virtualdetector.pos.z());
+        h_mcsteps_virtualdetector_pdg->Fill(mcsteps_virtualdetector.pdg);
+        h_mcsteps_virtualdetector_startCode->Fill(mcsteps_virtualdetector.startCode);
+        h_mcsteps_virtualdetector_stopCode->Fill(mcsteps_virtualdetector.stopCode);
+
+      }
+    }
   }
 
   file->Write();