EventSelector modifications for new MC waveform hits (#330)

* Adjust EventSelector logic to handle MC waveforms
Another downstream tool modification to adjust how the tool uses MC hits that are coming from the waveform generator.
* Update README.md
Added config variable that tells EventSelector to use MC waveform hits rather than the default MC hits.
* Update EventSelector.cpp
Added skip event for PMTWaveformSim
* Update EventSelector.cpp
Cleaned up logic
* Update EventSelector.cpp
- code condensed by merging data and MC logic together
- removed the event skipping check

Author: S81D

UserTools/EventSelector/EventSelector.cpp

@@ -14,6 +14,7 @@ bool EventSelector::Initialise(std::string configfile, DataModel &data){
 
   fPMTMRDOffset = false;
   fIsMC = true;
+  fMCWaveform = false;
   fPMTMRDOffset = 755;
   fRecoPDG = -1;
 
@@ -50,6 +51,7 @@ bool EventSelector::Initialise(std::string configfile, DataModel &data){
   }
   m_variables.Get("SaveStatusToStore", fSaveStatusToStore);
   m_variables.Get("IsMC",fIsMC);
+  m_variables.Get("PMTWaveformSim",fMCWaveform);
   m_variables.Get("RecoPDG",fRecoPDG);
   m_variables.Get("TriggerExtendedWindow",fTriggerExtended);
   m_variables.Get("BeamOK",fBeamOK);
@@ -94,7 +96,7 @@ bool EventSelector::Execute(){
   	Log("EventSelector Tool: No ANNIEEvent store!",v_error,verbosity); 
   	return false;
   };
-
+	
   // ANNIE Event number
   m_data->Stores.at("ANNIEEvent")->Get("EventNumber",fEventNumber);
 
@@ -631,12 +633,12 @@ bool EventSelector::EventSelectionByMCProjectedMRDHit() {
 
 bool EventSelector::EventSelectionByPMTMRDCoinc() {
 
-  if (fIsMC){
-    bool has_clustered_pmt = m_data->CStore.Get("ClusterMapMC",m_all_clusters_MC);
-    if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMapMC from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
+  if (!fIsMC || fMCWaveform) {
+	  bool has_clustered_pmt = m_data->CStore.Get("ClusterMap",m_all_clusters);
+      if (not has_clustered_pmt) { Log("EventSelector Tool: MCWaveform - Error retrieving ClusterMap from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
   } else {
-    bool has_clustered_pmt = m_data->CStore.Get("ClusterMap",m_all_clusters);
-    if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMap from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
+	  bool has_clustered_pmt = m_data->CStore.Get("ClusterMapMC",m_all_clusters_MC);
+        if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMapMC from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
   }
 
   bool has_clustered_mrd = m_data->CStore.Get("MrdTimeClusters",MrdTimeClusters);
@@ -649,8 +651,12 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
   }
 
   int pmt_cluster_size;
-  if (fIsMC) pmt_cluster_size = (int) m_all_clusters_MC->size();
-  else pmt_cluster_size = (int) m_all_clusters->size();
+  if (!fIsMC || fMCWaveform) {
+	  pmt_cluster_size = (int) m_all_clusters->size();
+  } else {
+	  pmt_cluster_size = (int) m_all_clusters_MC->size();
+  }
+	  
   m_data->Stores["RecoEvent"]->Set("NumPMTClusters",pmt_cluster_size);
   vec_pmtclusters_charge->clear();
   vec_pmtclusters_time->clear();
@@ -667,7 +673,42 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
 
   pmt_time = -1;
 
-  if (fIsMC){
+  // MC Waveform or Data
+  if (!fIsMC || fMCWaveform) {
+    if (m_all_clusters->size()){
+    double cluster_time;
+    for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
+      std::vector<Hit>&Hits = apair.second;
+      double time_temp = 0;
+      double charge_temp = 0;
+      for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){
+        time_temp+=Hits.at(i_hit).GetTime();
+        int tube = Hits.at(i_hit).GetTubeId();
+        // check if PMT is present in the map before accessing it
+	  auto it = ChannelNumToTankPMTSPEChargeMap->find(tube);
+	  if (it != ChannelNumToTankPMTSPEChargeMap->end()) {
+	    double charge_pe = Hits.at(i_hit).GetCharge() / it->second;
+	    charge_temp += charge_pe;
+	  } else {
+	    std::cerr << "PMT channel with hit not found in ChannelNumToTankPMTSPEChargeMap. Skipping this hit." << std::endl;
+	    continue;
+	  }
+      }
+      if (Hits.size()>0) time_temp/=Hits.size();
+      vec_pmtclusters_charge->push_back(charge_temp);
+      vec_pmtclusters_time->push_back(time_temp);
+      if (time_temp > 2000.) continue;	//not a prompt event
+      if (charge_temp > max_charge){
+        max_charge = charge_temp;
+        prompt_cluster = true;
+        pmt_time = time_temp;
+        n_hits = int(Hits.size());
+      }
+    }
+  }
+
+  // MC (parametric)
+  } else {
     if (m_all_clusters_MC->size()){
       double cluster_time;
       for(std::pair<double,std::vector<MCHit>>&& apair : *m_all_clusters_MC){
@@ -690,38 +731,6 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
         }
       }
     }
-  } else {
-    if (m_all_clusters->size()){
-      double cluster_time;
-      for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
-        std::vector<Hit>&Hits = apair.second;
-        double time_temp = 0;
-        double charge_temp = 0;
-        for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){
-          time_temp+=Hits.at(i_hit).GetTime();
-          int tube = Hits.at(i_hit).GetTubeId();
-          // check if PMT is present in the map before accessing it
-	  auto it = ChannelNumToTankPMTSPEChargeMap->find(tube);
-	  if (it != ChannelNumToTankPMTSPEChargeMap->end()) {
-	  	double charge_pe = Hits.at(i_hit).GetCharge() / it->second;
-	  	charge_temp += charge_pe;
-	  } else {
-	  	std::cerr << "PMT channel with hit not found in ChannelNumToTankPMTSPEChargeMap. Skipping this hit." << std::endl;
-	  	continue;
-	  }
-        }
-        if (Hits.size()>0) time_temp/=Hits.size();
-        vec_pmtclusters_charge->push_back(charge_temp);
-        vec_pmtclusters_time->push_back(time_temp);
-        if (time_temp > 2000.) continue;	//not a prompt event
-        if (charge_temp > max_charge){
-          max_charge = charge_temp;
-          prompt_cluster = true;
-          pmt_time = time_temp;
-          n_hits = int(Hits.size());
-        }
-      }
-    }
   }
 
   if (verbosity > 1) {
@@ -762,10 +771,14 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
   }
   m_data->Stores["RecoEvent"]->Set("MRDClustersTime",vec_mrdclusters_time, true);
 
-  if (fIsMC){
-    if (MrdTimeClusters.size() == 0 || m_all_clusters_MC->size() == 0) return false;
+  if (fIsMC) {
+    if (MrdTimeClusters.size() == 0 || (fMCWaveform ? m_all_clusters->size() == 0 : m_all_clusters_MC->size() == 0)) {
+        return false;
+    }
   } else {
-    if (MrdTimeClusters.size() == 0 || m_all_clusters->size() == 0) return false;
+      if (MrdTimeClusters.size() == 0 || m_all_clusters->size() == 0) {
+          return false;
+      }
   }
 
   pmtmrd_coinc_min = fPMTMRDOffset - 50;
@@ -1026,12 +1039,12 @@ bool EventSelector::FindPaddleChankey(double x, double y, int layer, unsigned lo
 
 bool EventSelector::EventSelectionByRecoPDG(int recoPDG, std::vector<double> & cluster_reco_pdg){
 
-  if (fIsMC){
-    bool has_clustered_pmt = m_data->CStore.Get("ClusterMapMC",m_all_clusters_MC);
-    if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMapMC from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
+  if (!fIsMC || fMCWaveform) {
+	bool has_clustered_pmt = m_data->CStore.Get("ClusterMap",m_all_clusters);
+	if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMap from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
   } else {
-    bool has_clustered_pmt = m_data->CStore.Get("ClusterMap",m_all_clusters);
-    if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMap from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
+	  bool has_clustered_pmt = m_data->CStore.Get("ClusterMapMC",m_all_clusters_MC);
+      if (not has_clustered_pmt) { Log("EventSelector Tool: Error retrieving ClusterMapMC from CStore, did you run ClusterFinder beforehand?",v_error,verbosity); return false; }
   }
 
   std::map<double,double> ClusterMaxPEs;
@@ -1045,52 +1058,51 @@ bool EventSelector::EventSelectionByRecoPDG(int recoPDG, std::vector<double> & c
   bool found_pdg = false;
 
   if (fabs(recoPDG)==2112){
-    if (fIsMC){
-      if (m_all_clusters_MC->size()){
-        for(std::pair<double,std::vector<MCHit>>&& apair : *m_all_clusters_MC){
-          double cluster_time = apair.first;
-          double charge_balance = ClusterChargeBalances.at(cluster_time);
-          std::vector<MCHit>&MCHits = apair.second;
-          double time_temp = 0;
-          double charge_temp = 0;
-          for (unsigned int i_hit = 0; i_hit < MCHits.size(); i_hit++){
-            time_temp+=MCHits.at(i_hit).GetTime();
-            charge_temp+=MCHits.at(i_hit).GetCharge();
-          }
-          if (cluster_time > 10000 && charge_balance < 0.4 && charge_temp < 120) {
-            cluster_reco_pdg.push_back(cluster_time);
-            found_pdg = true;
-            std::cout <<"Found neutron candidate for cluster at time = "<<cluster_time<<", with CB = "<<charge_balance <<" and charge "<<charge_temp<<"!!!"<<std::endl;
-          } else {
-            std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
-          }
-        }
-      }
-    } else {
-      if (m_all_clusters->size()){
-        double cluster_time;
-        for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
-          double cluster_time = apair.first;
-          double charge_balance = ClusterChargeBalances.at(cluster_time);
-          std::vector<Hit>&Hits = apair.second;
-          double time_temp = 0;
-          double charge_temp = 0;
-          for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){
-            time_temp+=Hits.at(i_hit).GetTime();
-            int tube = Hits.at(i_hit).GetTubeId();
-            double charge_pe = Hits.at(i_hit).GetCharge()/ChannelNumToTankPMTSPEChargeMap->at(tube);
-            charge_temp+=charge_pe;
-          }
-          if (cluster_time > 10000 && charge_balance < 0.4 && charge_temp < 120) {
-            cluster_reco_pdg.push_back(cluster_time);
-            found_pdg = true;
-            std::cout <<"Found neutron candidate for cluster at time = "<<cluster_time<<", with CB = "<<charge_balance <<"and charge "<<charge_temp<<"!!!"<<std::endl;
-          } else {
-            std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
+	if (!fIsMC || fMCWaveform) {    // data-like
+		if (m_all_clusters->size()){
+	        for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
+	          double cluster_time = apair.first;
+	          double charge_balance = ClusterChargeBalances.at(cluster_time);
+	          std::vector<Hit>&Hits = apair.second;
+	          double time_temp = 0;
+	          double charge_temp = 0;
+	          for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){
+	            time_temp+=Hits.at(i_hit).GetTime();
+	            int tube = Hits.at(i_hit).GetTubeId();
+	            double charge_pe = Hits.at(i_hit).GetCharge()/ChannelNumToTankPMTSPEChargeMap->at(tube);
+	            charge_temp+=charge_pe;
+	          }
+	          if (cluster_time > 10000 && charge_balance < 0.4 && charge_temp < 120) {
+	            cluster_reco_pdg.push_back(cluster_time);
+	            found_pdg = true;
+	            std::cout <<"Found neutron candidate for cluster at time = "<<cluster_time<<", with CB = "<<charge_balance <<"and charge "<<charge_temp<<"!!!"<<std::endl;
+	          } else {
+	            std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
+	          }
+        	}
+     	 }
+	  } else {     // MCHits
+		if (m_all_clusters_MC->size()){
+          for(std::pair<double,std::vector<MCHit>>&& apair : *m_all_clusters_MC){
+            double cluster_time = apair.first;
+            double charge_balance = ClusterChargeBalances.at(cluster_time);
+            std::vector<MCHit>&MCHits = apair.second;
+            double time_temp = 0;
+            double charge_temp = 0;
+            for (unsigned int i_hit = 0; i_hit < MCHits.size(); i_hit++){
+              time_temp+=MCHits.at(i_hit).GetTime();
+              charge_temp+=MCHits.at(i_hit).GetCharge();
+            }
+            if (cluster_time > 10000 && charge_balance < 0.4 && charge_temp < 120) {
+              cluster_reco_pdg.push_back(cluster_time);
+              found_pdg = true;
+              std::cout <<"Found neutron candidate for cluster at time = "<<cluster_time<<", with CB = "<<charge_balance <<" and charge "<<charge_temp<<"!!!"<<std::endl;
+            } else {
+              std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
+            }
           }
         }
-      }
-    }
+	  }
   }
 
   return found_pdg;

UserTools/EventSelector/EventSelector.h

@@ -249,6 +249,7 @@ class EventSelector: public Tool {
   bool fThroughGoing = false;
   bool fEventCutStatus;
   bool fIsMC; 
+  bool fMCWaveform;
   int fTriggerWord;
   int fRecoPDG;
   bool fTriggerExtended = false;

UserTools/EventSelector/README.md

@@ -106,5 +106,6 @@ ThroughGoing
 TriggerWord
 RecoPDG
 IsMC
+PMTWaveformSim
 SaveStatusToStore
 ```