Merge pull request #314 from AndrewEdmonds11/rooutil-trig

Updates to trigger branch for RooUtil

Author: AndrewEdmonds11

README.md

@@ -3,7 +3,7 @@
 ## Table of Contents
 1. [Introduction](#Introduction)
 2. [Structure](#Structure)
-3. [How to get the version of EventNtuple](#How-to-get-the-version-of-EventNtuple)
+3. [How to Check the EventNtuple](#How-to-Check-the-EventNtuple)
 4. [How to Analyze an EventNtuple](#How-to-Analyze-an-EventNtuple)
 5. [How to Create an EventNtuple](#How-to-Create-an-EventNtuple)
 6. [Validation](#Validation)
@@ -24,13 +24,16 @@ A [list of branches is available](./doc/branches.md)
 
 The help understand what all the branches and leaves mean, we have an [```ntuplehelper```](doc/ntuplehelper.md) tool
 
-## How to get the version of EventNtuple that is in a file
-For versions of EventNtuple later than v6.3.0, we store the version number in the ROOT file. You can check the version number with the ```checkEventNtuple``` tool:
+## How to Check the EventNtuple
+You can get some basic information about the EventNtuple file useing ```checkEventNtuple``` like so:
 
 ```
 checkEventNtuple file1.root file2.root
 ```
 
+This will print the version number (for versions after v6.3.0) and the trigger branches in the file (for versions after v6.8.0)
+
+
 ## How to Analyze an EventNtuple
 To help with analyzing the EventNtuple given its complex structure, we have two sets of utilities:
 * [```RooUtil```](rooutil/README.md) for ROOT-based analyses, and

bin/checkEventNtuple

@@ -15,20 +15,35 @@ def main():
 
     foldername = "EventNtuple"
     histname = "version"
+    treename = "ntuple"
     for filename in args.filenames:
         try:
             f = ROOT.TFile(filename, "READ");
             if (f.GetListOfKeys().Contains(foldername)):
                 folder = f.Get(foldername)
-                if (folder.GetListOfKeys().Contains("version")):
-                    h = f.Get("EventNtuple/version")
+                # Check version number
+                if (folder.GetListOfKeys().Contains(histname)):
+                    h = f.Get(foldername+"/"+histname)
                     print("{} contains EventNtuple v{:02.0f}_{:02.0f}_{:02.0f}".format(filename, h.GetBinContent(1), h.GetBinContent(2), h.GetBinContent(3)))
                 else:
                     print("{} histogram does not exist in {} (it is either v06_02_00 or older)".format(histname, filename))
+
+                # Print triggers
+                if (folder.GetListOfKeys().Contains(treename)):
+                    t = f.Get(foldername+"/"+treename)
+                    print("\nIt contains the following trigger paths:")
+                    for br in t.GetListOfBranches():
+                        if "trig_" in br.GetName():
+                            print(f" - {br.GetName()}")
+                    print("For trigger definitions, please contact the Trigger group.")
+                else:
+                    print("{} tree does not exist in {}...".format(treename, filename))
             else:
                 print("{}/ folder does not exist in {}".format(foldername, filename))
         except OSError:
             print("{} threw an OSError (does it exist?)".format(filename))
 
+    print("\nFor any questions about EventNtuple, please contact the Analysis Tools group.")
+
 if __name__ == "__main__":
     main()

doc/branches.md

@@ -98,6 +98,11 @@ The branch is empty if there are no CRV hit during the event.
 | crvcoincs |  Vector branch |   information about a cluster of CRV coincidence triplets| [see CrvHitInfoReco.hh](../inc/CrvHitInfoReco.hh)
 | crvcoincsmc |  Vector branch |   information about the MC track which most likely caused the CRV coincidence triplets| [see CrvHitInfoMC.hh](../inc/CrvHitInfoMC.hh)
 | crvcoincsmcplane |  Vector branch |   information about the point where the MC trajectory crosses the xz plane of CRV-T| [see CrvPlaneInfoMC.hh](../inc/CrvPlaneInfoMC.hh)
+## Trigger Branches
+
+Each element in these branch corresponds to a different Kalman fit hypotheses to reconstruct the track:
+
+Trigger branches are prefixed with ```trig_```. To find the names of the trigger branches run ```checkEventNtuple filename.root```
 ## Deprecated Branches
 
 These branches are remnants from trkana and are deprecated.

helper/ntuplehelper.py

@@ -230,6 +230,15 @@ def list_all_branches(self, export_to_md=False):
                 tokens=explanation.split(":")
                 print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")")
 
+        if not export_to_md:
+            print("\nTrigger Branches")
+            print("================")
+            print("Trigger branches are prefixed with \"trig_\". To find the names of the trigger branches run \"checkEventNtuple filename.root\"")
+        else:
+            print("## Trigger Branches\n")
+            print("Each element in these branch corresponds to a different Kalman fit hypotheses to reconstruct the track:\n")
+            print("Trigger branches are prefixed with ```trig_```. To find the names of the trigger branches run ```checkEventNtuple filename.root```")
+
         if not export_to_md:
             print("\nDeprecated Branches")
             print("================")

inc/TrigInfo.hh

@@ -9,10 +9,10 @@ namespace mu2e
 {
   struct TrigInfo {
 
-    static const unsigned int ntrig_ = 30; 
+    static const unsigned int ntrig_ = 30;
+
+    bool _triggerArray[ntrig_];
 
-    UChar_t _triggerArray[ntrig_];
-    
     void reset() { *this = TrigInfo(); }
   };
 }

rooutil/README.md

@@ -113,8 +113,9 @@ Some branches are not contained in any of the above classes:
 * ```crvpulses``` and ```crvpulsesmc```
 * ```crvcoincsmcplane```
 * ```calohits```, ```calorecodigis```, ```calodigis```
+* ```trig_``` branches
 
-Reach out to the developers on the #analysis-tools Slack channel if you have a need to have these added somewhere.
+Reach out to the developers on the #analysis-tools Slack channel if you need to have these added somewhere.
 
 ## Accessing User-Friendly Classes
 The ```Event``` class provides access to ```Tracks``` and ```CrvCoincs```:

rooutil/examples/PlotEntranceMomentumTriggeredEvents.C

@@ -0,0 +1,64 @@
+//
+// An example of how to plot the momentum of electrons at the tracker entrance
+// This uses cut functions defined in common_cuts.hh
+//
+
+#include "EventNtuple/rooutil/inc/RooUtil.hh"
+#include "EventNtuple/rooutil/inc/common_cuts.hh"
+
+#include "TH1F.h"
+
+void PlotEntranceMomentumTriggeredEvents(std::string filename) {
+
+  // Create the histogram you want to fill
+  TH1F* hRecoMomAll = new TH1F("hRecoMomAll", "Reconstructed Momentum at Tracker Entrance", 50,95,110);
+  TH1F* hRecoMom_tprDe_highP_stopTarg = new TH1F("hRecoMom_tprDe_highP_stopTarg", "Reconstructed Momentum at Tracker Entrance", 50,95,110);
+  TH1F* hRecoMom_tprDe_highP = new TH1F("hRecoMom_tprDe_highP", "Reconstructed Momentum at Tracker Entrance", 50,95,110);
+
+  // 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
+    auto e_minus_tracks = event.GetTracks(is_e_minus);
+
+    // Loop through the e_minus tracks
+    for (auto& track : e_minus_tracks) {
+
+      // Get the track segments at the tracker entrance
+      auto trk_ent_segments = track.GetSegments([](TrackSegment& segment){ return tracker_entrance(segment) && has_reco_step(segment); });
+
+      // Loop through the tracker entrance track segments
+      for (auto& segment : trk_ent_segments) {
+
+        // Fill the histogram
+        hRecoMomAll->Fill(segment.trkseg->mom.R());
+
+        if (passes_trigger(event, "tprDe_highP_stopTarg")) {
+          hRecoMom_tprDe_highP_stopTarg->Fill(segment.trkseg->mom.R());
+        }
+        if (passes_trigger(event, "tprDe_highP_stopTarg")) {
+          hRecoMom_tprDe_highP->Fill(segment.trkseg->mom.R());
+        }
+      }
+    }
+  }
+
+  // Draw the histogram
+  hRecoMomAll->SetLineColor(kBlack);
+  hRecoMomAll->SetLineWidth(2);
+  hRecoMomAll->Draw("HIST E");
+
+  hRecoMom_tprDe_highP->SetLineColor(kBlue);
+  hRecoMom_tprDe_highP->SetLineWidth(2);
+  hRecoMom_tprDe_highP->Draw("HIST E SAME");
+
+  hRecoMom_tprDe_highP_stopTarg->SetLineColor(kRed);
+  hRecoMom_tprDe_highP_stopTarg->SetLineWidth(2);
+  hRecoMom_tprDe_highP_stopTarg->Draw("HIST E SAME");
+
+}

rooutil/examples/PrintEvents.C

@@ -249,5 +249,11 @@ void PrintEvents(std::string filename) {
         std::cout << "calodigi: " << calodigi.SiPMID_ << ", " << calodigi.t0_ << ", " << calodigi.peakpos_ << ", " << calodigi.caloRecoDigiIdx_ << std::endl;
       }
     }
+
+    // trigger branches
+    for (const auto pair : event.trigNameMap) {
+      std::cout << pair.first << ": " << event.triginfo._triggerArray[pair.second] << std::endl;
+    }
+
   }
 }

rooutil/inc/Event.hh

@@ -35,6 +35,8 @@
 
 #include "EventNtuple/inc/MVAResultInfo.hh"
 
+#include "EventNtuple/inc/TrigInfo.hh"
+
 #include "EventNtuple/rooutil/inc/Track.hh"
 #include "EventNtuple/rooutil/inc/CrvCoinc.hh"
 #include "EventNtuple/rooutil/inc/CaloCluster.hh"
@@ -43,7 +45,6 @@
 
 struct Event {
   Event(TChain* ntuple) {
-
     if (ntuple->GetBranch("evtinfo")) {
       ntuple->SetBranchAddress("evtinfo", &this->evtinfo);
     }
@@ -53,6 +54,16 @@ struct Event {
     if (ntuple->GetBranch("crvsummary")) {
       ntuple->SetBranchAddress("crvsummary", &this->crvsummary);
     }
+    const auto& branches = ntuple->GetListOfBranches();
+    int i_trig_branch = 0;
+    for (const auto& branch : *branches) {
+      std::string brname = branch->GetName();
+      if (brname.substr(0, 5) == "trig_") {
+        ntuple->SetBranchAddress(brname.c_str(), &this->triginfo._triggerArray[i_trig_branch]);
+        trigNameMap.insert({brname, i_trig_branch});
+        i_trig_branch++;
+      }
+    }
 
     if (ntuple->GetBranch("trk")) {
       ntuple->SetBranchAddress("trk", &this->trk);
@@ -363,6 +374,7 @@ struct Event {
   mu2e::HitCount* hitcount = nullptr;
   mu2e::CrvSummaryReco* crvsummary = nullptr;
   mu2e::CrvSummaryMC* crvsummarymc = nullptr;
+  mu2e::TrigInfo triginfo; // not a pointer because we give the address of array elements inside this
 
   std::vector<mu2e::TrkInfo>* trk = nullptr;
   std::vector<mu2e::TrkInfoMC>* trkmc = nullptr;
@@ -393,6 +405,9 @@ struct Event {
   std::vector<mu2e::CrvPlaneInfoMC>* crvcoincsmcplane = nullptr;
 
   std::vector<std::vector<mu2e::SimInfo>>* trkmcsim = nullptr;
+
+  // Need to keep track of trigger name to element in triginfo
+  std::map<std::string, unsigned int> trigNameMap;
 };
 
 #endif

rooutil/inc/RooUtil.hh

@@ -47,14 +47,16 @@ public:
 
   void SetVersionNumber(std::string filename) {
     TFile* file = new TFile(filename.c_str(), "READ");
-    TH1I* hVersion = (TH1I*) file->Get("EventNtuple/version");
-    if (!hVersion) {
+    TH1I* version = (TH1I*) file->Get("EventNtuple/version");
+    if (!version) {
       std::cout << "Warning: this EventNtuple file does not contain a version number. It is either v06_02_00 or older. This is just a warning..." << std::endl;
     }
     else {
-      majorVer = hVersion->GetBinContent(1);
-      minorVer = hVersion->GetBinContent(2);
-      patchVer = hVersion->GetBinContent(3);
+      hVersion = (TH1I*) version->Clone();
+      hVersion->SetDirectory(0);
+      int majorVer = hVersion->GetBinContent(1);
+      int minorVer = hVersion->GetBinContent(2);
+      int patchVer = hVersion->GetBinContent(3);
       std::cout << "EventNtuple v" << std::setw(2) << std::setfill('0') << majorVer << "_"
                 << std::setw(2) << std::setfill('0') << minorVer << "_"
                 << std::setw(2) << std::setfill('0') << patchVer << std::endl;
@@ -140,6 +142,12 @@ public:
 
     if(event->trkmcsim) { output_ntuple->Branch("trkmcsim", event->trkmcsim); }
 
+    for (const auto pair : event->trigNameMap) {
+      output_ntuple->Branch(pair.first.c_str(), &event->triginfo._triggerArray[pair.second]);
+    }
+
+    // Write out histograms from input to output
+    hVersion->Write();
   }
 
   void FillOutputEventNtuple() {
@@ -151,9 +159,7 @@ private:
   Event* event; // holds all the variables for SetBranchAddress
   bool debug;
 
-  int majorVer;
-  int minorVer;
-  int patchVer;
+  TH1I* hVersion;
 
   TTree* output_ntuple; // for output
 };