Merge pull request #355 from AndrewEdmonds11/rooutil-calohitmc

Add CaloHitMC to RooUtil analyzer class

Author: AndrewEdmonds11

rooutil/README.md

@@ -106,7 +106,7 @@ The ```EventNtupleTimeCluster``` class contains all information related to a sin
 The ```CaloCluster``` class contains all information related to a single calorimeter cluster
 
 * single objects: ```calocluster```, ```caloclustermc```
-* vector objects: ```calohits```, ```calomcsim```
+* vector objects: ```calohits```, ```calomcsim```, ```calohitsmc```
 
 Examples: [PlotCaloClusterEnergy.C](./examples/PlotCaloClusterEnergy.C), [PlotCaloClusterEnergy_RecoVsTrue.C](./examples/PlotCaloClusterEnergy_RecoVsTrue.C), [PlotCaloClusterAndHits.C](./examples/PlotCaloClusterAndHits.C), [PlotCaloCluster_SimParticles.C](./examples/PlotCaloCluster_SimParticles.C)
 
@@ -125,10 +125,12 @@ Some branches are not contained in any of the above classes:
 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```:
+The ```Event``` class provides access to ```Tracks```, ```CaloClusters``` and ```CrvCoincs```:
 
 * ```CountTracks()```: counts the number of tracks in the event
 * ```GetTracks()```: gets the tracks (passes you a copy)
+* ```CountCaloClusters()```: counts the number of calo clusters in the event
+* ```GetCaloClusters()```: gets the calo clusters (passes you a copy)
 * ```CountCrvCoincs()```: counts the number of CRV coincidences in the event
 * ```GetCrvCoincs()```: gets the CRV coincidences
 
@@ -137,6 +139,10 @@ The ```Track``` class provides access to the ```TrackSegments```, the ```TrackHi
 * ```GetHits()```, ```CountHits()```
 * ```GetMCParticles()```, ```CountMCParticles()```
 
+The ```CaloCluster``` class provides access to the ```CaloHits``` and the ```MCParticles```:
+* ```GetHits()```, ```CountHits()```
+* ```GetMCParticles()```, ```CountMCParticles()```
+
 All of these can be passed a cut function (see below) to count / select just a subset of the objects
 
 ## Cut Functions

rooutil/examples/PlotCaloHits_RecoVsTrue.C

@@ -0,0 +1,48 @@
+//
+// 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 "TH2F.h"
+
+using namespace rooutil;
+
+void PlotCaloHits_RecoVsTrue(std::string filename) {
+
+  // Set up RooUtil
+  RooUtil util(filename);
+
+  TH2F* hCaloHits_RecoVsTrue = new TH2F("hCaloHits_RecoVsTrue", "Reco vs True (Calo Hits)", 110,0,110, 110,0,110);
+  hCaloHits_RecoVsTrue->SetXTitle("Reco EDep [MeV]");
+  hCaloHits_RecoVsTrue->SetYTitle("True EDep [MeV]");
+
+  // Loop through the events
+  int n_events = util.GetNEvents();
+  for (int i_event = 0; i_event < n_events; ++i_event) {
+    // Get the next event
+    auto& event = util.GetEvent(i_event);
+
+        // Get the e_minus tracks from the event
+    auto calo_clusters = event.GetCaloClusters();
+
+    // Loop through the calo clusters
+    int i_cluster = 0;
+    for (auto& calo_cluster : calo_clusters) {
+
+      // Get the calohits for this cluster
+      auto calohits = calo_cluster.GetHits([](CaloHit& hit){ return has_mc_hit(hit) && has_reco_hit(hit); });
+
+      // Loop through the calo hits
+      for (auto& calohit : calohits) {
+        hCaloHits_RecoVsTrue->Fill(calohit.reco->eDep_, calohit.mc->eDep);
+      }
+    }
+  }
+
+  TCanvas* c1 = new TCanvas();
+  hCaloHits_RecoVsTrue->Draw("COLZ");
+}

rooutil/examples/PrintEvents.C

@@ -287,6 +287,13 @@ void PrintEvents(std::string filename) {
       }
     }
 
+    // calohitsmc branch
+    if (event.calohitsmc != nullptr) {
+      for (const auto& calohitmc : *(event.calohitsmc)) {
+        std::cout << "calohitmc: " << calohitmc.nsim << ", " << calohitmc.eDep << ", " << calohitmc.eDepG4 << ", " << calohitmc.clusterIdx_ << std::endl;
+      }
+    }
+
 
     // trigger branches
     for (const auto& pair : event.trigger.NameToIndexMap()) {

rooutil/inc/CaloHit.hh

@@ -2,16 +2,16 @@
 #define CaloHit_hh_
 
 #include <functional>
-#include "EventNtuple/inc/TrkStrawHitInfo.hh"
-#include "EventNtuple/inc/TrkStrawHitInfoMC.hh"
-#include "EventNtuple/inc/TrkStrawHitCalibInfo.hh"
+#include "EventNtuple/inc/CaloHitInfo.hh"
+#include "EventNtuple/inc/CaloHitInfoMC.hh"
 
 namespace rooutil {
   struct CaloHit {
     CaloHit() { }
 
     // Pointers to the data
     mu2e::CaloHitInfo* reco = nullptr;
+    mu2e::CaloHitInfoMC* mc = nullptr;
   };
 
   typedef std::function<bool(CaloHit&)> CaloHitCut;

rooutil/inc/Event.hh

@@ -185,6 +185,10 @@ namespace rooutil {
               CaloHit calohit;
               calohit.reco = &(calohits->at(calohit_Idx)); // passing the addresses of the underlying structs
 
+              if (calohitsmc != nullptr) { // 1-to-1 matching of calohits and calohitsmc
+                calohit.mc = &(calohitsmc->at(calohit_Idx));
+              }
+
               calo_cluster.hits.emplace_back(calohit);
             }
           }

rooutil/inc/common_cuts.hh

@@ -317,4 +317,15 @@ bool has_reco_cluster(CaloCluster& cluster) { // CaloCluster has a reco object
   else { return false; }
 }
 
+//+ Calo Hit Cuts - Other
+bool has_mc_hit(CaloHit& hit) { // calo hit has MC truth info
+  if (hit.mc != nullptr) { return true; }
+  else { return false; }
+}
+bool has_reco_hit(CaloHit& hit) { // calo hit has reco info
+  if (hit.reco != nullptr) { return true; }
+  else { return false; }
+}
+
+
 #endif