Neutron Multiplicity toolchain (#245)

* added framework for simple reco

* added Factory & Unity

* added muon and neutron identification cuts in ClusterClassifiers

* add FindNeutrons tool

* added FindNeutrons

* new tool

* updated Makefile & Setup.sh

* fixed compilation errors in FindNeutrons & ClusterClassifiers

* more compilation fixes

* moved neutron finding algorithms to FindNeutrons tool

* SimpleReconstruction & FindNeutrons progress

* wip SimpleReconstruction

* first version SimpleReco + FindNeutrons

* first working version of simple reco

* added extended window and beamok flag to EventSelector

* fix for variables SimpleReco/NeutronMult tools

* add tree framework to NeutronMultiplicity tool

* fill more variables in output tree

* added BoostStore read-in and store mode for NeutronMultiplicity

* Update README of FindNeutrons

Added README of FindNeutrons tool with relevant information

* Update README of SimpleReconstruction

Added README with basic properties of SimpleReconstruction tool

* Update README of NeutronMultiplicity tool

Added README with description of properties of NeutronMultiplicity tool.

* moved some variables to RecoEvent instead of ANNIEEvent

* Update visible_energy_vertex.C

Removed unnecessary parts of the code (which were commented out)

* Update neutrino_selection.C

Added verbose configuration variable

* Update reconstruction_migration.C

Added verbose configuration variable

* Update and rename plot_neutrons_mc_new.C to plot_neutrons_mc.C

rename + add configuration variable verbose

* Update and rename plot_neutrons_dirt.C to plot_neutrons_data_dirt.C

Rename + add verbose configuration variable

* Update and rename plot_neutrons_data_new.C to plot_neutrons_data_beam.C

Rename + add verbose configuration variable

* Create README.md

Added README for the root-scripts

* Update README.md

Extended description

* Update README.md

Added text to README file

* Update README.md

Added additional sentence about ReadFromBoostStore version of the toolchain

* wip NeutronMultiplicity

* re-add LoadGenieEvent & PrintGenieEvent

* Add FindParticlePdgs to MCRecoEventLoader

* more MRD variables in SimpleReco tool

* add MC version of NeutronMultiplicity toolchain

* Add Particles object in LoadWCSim

* verbosity changes

* add primary/secondary information for neutron captures

* more verbosity cleanup

* add efficiency map in FindNeutrons

* add effiency map in NeutronMultiplicity

* Add transverse momentum in SimpleReconstruction

* Add neutron efficiency scripts (MC)

* bugfix in neutronmultiplicity efficiency correction

* add scripts for running on many files

* Reduced verbosity of NeutronMultiplicity tool

* Added NeutronMultiplicity script for filtered ANNIEEvent BoostStore input files

* config file madness

* removed duplicate filelist file

* Update README.md

update readme of NeutronMultiplicity tool

* Update README.md

update README of FindNeutrons tool

* Add NHits method for neutron candidates

* add event offset to LoadGenieEvent tool

* toolchain modifications for different neutron cuts

UserTools/ClusterClassifiers/ClusterClassifiers.cpp

@@ -62,6 +62,9 @@ bool ClusterClassifiers::Execute(){
   std::map<double,double> ClusterMaxPEs;
   std::map<double,Position> ClusterChargePoints;
   std::map<double,double> ClusterChargeBalances;
+  std::map<double,double> ClusterTotalPEs;
+  std::vector<double> ClusterTimes;
+  std::map<double,int> ClusterNHits;
 
   if (isData){
     for (std::pair<double,std::vector<Hit>>&& cluster_pair : *m_all_clusters) {
@@ -74,6 +77,10 @@ bool ClusterClassifiers::Execute(){
       ClusterChargeBalances.emplace(cluster_time,ChargeBalance);
       double max_PE = this->CalculateMaxPE(cluster_hits);
       ClusterMaxPEs.emplace(cluster_time,max_PE);
+      double total_PE = this->CalculateTotalPE(cluster_hits);
+      ClusterTotalPEs.emplace(cluster_time,total_PE);
+      ClusterTimes.push_back(cluster_time);
+      ClusterNHits.emplace(cluster_time,int(cluster_hits.size()));
     }
   } else {
      for (std::pair<double,std::vector<MCHit>>&& cluster_pair : *m_all_clusters_MC) {
@@ -87,6 +94,10 @@ bool ClusterClassifiers::Execute(){
        ClusterChargeBalances.emplace(cluster_time,ChargeBalance);
        double max_PE = this->CalculateMaxPEMC(cluster_hits,cluster_detkey);
        ClusterMaxPEs.emplace(cluster_time,max_PE);
+       double total_PE = this->CalculateTotalPEMC(cluster_hits,cluster_detkey);
+       ClusterTotalPEs.emplace(cluster_time,total_PE);
+       ClusterTimes.push_back(cluster_time);
+       ClusterNHits.emplace(cluster_time,int(cluster_hits.size()));
      }
   }
 
@@ -95,6 +106,19 @@ bool ClusterClassifiers::Execute(){
   m_data->Stores.at("ANNIEEvent")->Set("ClusterChargePoints", ClusterChargePoints);
   m_data->Stores.at("ANNIEEvent")->Set("ClusterChargeBalances", ClusterChargeBalances);
   m_data->Stores.at("ANNIEEvent")->Set("ClusterMaxPEs", ClusterMaxPEs);
+  m_data->Stores.at("ANNIEEvent")->Set("ClusterTotalPEs", ClusterTotalPEs);
+  m_data->Stores.at("ANNIEEvent")->Set("ClusterNHits", ClusterNHits);
+
+  //identify prompt muon candidate
+  bool found_prompt_muon = this->IdentifyPromptMuonCluster(ClusterTotalPEs);
+
+  //store indices of muon and neutron clusters to ANNIEEvent Store
+  m_data->Stores.at("RecoEvent")->Set("ClusterIndexPromptMuon", prompt_muon_index);
+  if (found_prompt_muon){
+    m_data->Stores.at("RecoEvent")->Set("PromptMuonTotalPE", ClusterTotalPEs.at(ClusterTimes.at(prompt_muon_index)));
+    m_data->Stores.at("RecoEvent")->Set("PromptMuonTime", ClusterTimes.at(prompt_muon_index));
+  }
+
   return true;
 }
 
@@ -272,3 +296,94 @@ double ClusterClassifiers::CalculateMaxPEMC(std::vector<MCHit> cluster_hits, std
    if(verbosity>4) std::cout << "ClusterClassifiers Tool: Calculated max PE hit of " << max_PE << std::endl;
    return max_PE;
  }
+
+bool ClusterClassifiers::IdentifyPromptMuonCluster(std::map<double,double> cluster_totalq){
+
+  bool return_prompt = false;
+  int tmp_cluster_id = 0;
+  int cluster_muon = -1;
+  double max_pe = 0;
+  for (std::map<double,double>::iterator it = cluster_totalq.begin(); it!= cluster_totalq.end(); it++){
+    //Check if the cluster charge is higher than the current maximum
+    //Only consider clusters in the prompt window (time < 2000 ns)
+    if (it->second > max_pe && it->first < 2000){
+      max_pe = it->second;
+      cluster_muon = tmp_cluster_id;	//For now, associate the cluster with the highest charge with the muon
+      return_prompt = true;
+    }
+    tmp_cluster_id ++;
+  }
+
+  prompt_muon_index = cluster_muon; 
+
+  return return_prompt;
+}
+
+bool ClusterClassifiers::IdentifyDelayedNeutronClusters(std::map<double,double> cluster_cb, std::map<double,double> cluster_totalq){
+
+  bool return_delayed = false;
+  int tmp_cluster_id = 0;
+  std::vector<int> cluster_neutron;
+  for (std::map<double,double>::iterator it = cluster_totalq.begin(); it!= cluster_totalq.end(); it++){
+    //Check if the cluster is in the delayed window and has a time > 10 us (exclude afterpulses)
+    //The window should should be extended in the future after relevant exlusion cuts for afterpulsing have been implemented
+    //check if the charge balance cut for neutrons is passed -> consider a neutron candidate
+    //Improve the neutron selection cuts in the future, probably cutting more signal than necessary at the moment
+    if (it->first > 10000){
+      double current_cb = cluster_cb.at(it->first);
+      double current_q = cluster_totalq.at(it->first);
+      if (current_cb < 0.4 && current_q < 150 && (current_cb <= (1. - current_q/150.)*0.5)){
+        cluster_neutron.push_back(tmp_cluster_id);
+        return_delayed = true;
+      }
+    }
+    tmp_cluster_id ++;
+  }
+
+  delayed_neutron_index = cluster_neutron;
+
+  return return_delayed; 
+
+}
+
+double ClusterClassifiers::CalculateTotalPE(std::vector<Hit> cluster_hits)
+{ 
+  double total_PE = 0;
+  for (int i = 0; i < (int) cluster_hits.size(); i++){
+    Hit ahit = cluster_hits.at(i); 
+    double hit_charge = ahit.GetCharge();
+    int channel_key = ahit.GetTubeId();
+    double hit_PE = 0;
+    std::map<int, double>::iterator it = ChannelKeyToSPEMap.find(channel_key);
+    if(it != ChannelKeyToSPEMap.end()){ //Charge to SPE conversion is available
+      hit_PE  = hit_charge / ChannelKeyToSPEMap.at(channel_key);
+      total_PE += hit_PE;
+    }
+  }
+  if(verbosity>4) std::cout << "ClusterClassifiers Tool: Calculated total PE of " << total_PE << std::endl;
+  return total_PE;
+}
+
+double ClusterClassifiers::CalculateTotalPEMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys)
+ {
+   double total_PE = 0;
+   for (int i = 0; i < (int) cluster_hits.size(); i++){
+     MCHit ahit = cluster_hits.at(i); 
+     double hit_charge = ahit.GetCharge();
+     //int channel_key = ahit.GetTubeId();
+     unsigned long detkey = cluster_detkeys.at(i);
+     int channel_key = (int) detkey;
+     int tubeid = ahit.GetTubeId();
+     unsigned long utubeid = (unsigned long) tubeid;
+     int wcsimid = channelkey_to_pmtid.at(utubeid);
+     unsigned long detkey_data = pmtid_to_channelkey[wcsimid];
+     int channel_key_data = (int) detkey_data;
+     std::map<int, double>::iterator it = ChannelKeyToSPEMap.find(channel_key_data);
+     if(it != ChannelKeyToSPEMap.end()){ //Charge to SPE conversion is available
+       total_PE += hit_charge;
+     }
+   }
+   if(verbosity>4) std::cout << "ClusterClassifiers Tool: Calculated total PE of " << total_PE << std::endl;
+   return total_PE;
+ }
+

UserTools/ClusterClassifiers/ClusterClassifiers.h

@@ -27,12 +27,16 @@ class ClusterClassifiers: public Tool {
   bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools.
   bool Execute(); ///< Execute function used to perform Tool purpose.
   bool Finalise(); ///< Finalise function used to clean up resources.
-  Position CalculateChargePoint(std::vector<Hit> cluster_hits);
-  double CalculateChargeBalance(std::vector<Hit> cluster_hits);
-  double CalculateMaxPE(std::vector<Hit> cluster_hits);
-  Position CalculateChargePointMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys);
-  double CalculateChargeBalanceMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys);
-  double CalculateMaxPEMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys);
+  Position CalculateChargePoint(std::vector<Hit> cluster_hits); ///< function to calculate the charge points for clusters (data)
+  double CalculateChargeBalance(std::vector<Hit> cluster_hits); ///< function to calculate the charge balance values for clusters (data)
+  double CalculateMaxPE(std::vector<Hit> cluster_hits); ///< function to identify the MaxPE value recorded by a single PMT in a cluster (data)
+  Position CalculateChargePointMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys); ///< function to calculate the charge points for clusters (MC)
+  double CalculateChargeBalanceMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys); ///< function to calculate the charge balance values for clusters (MC)
+  double CalculateMaxPEMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys); ///< function to calculate the MaxPE value recorded by a single PMT in a cluster (MC)
+  double CalculateTotalPE(std::vector<Hit> cluster_hits); ///< function to calculate the total PE recorded in a cluster (data)
+  double CalculateTotalPEMC(std::vector<MCHit> cluster_hits, std::vector<unsigned long> cluster_detkeys); ///< function to calculate the total PE recorded in a cluster (MC)
+  bool IdentifyPromptMuonCluster(std::map<double,double> cluster_totalq); ///< function to identify the prompt muon cluster
+  bool IdentifyDelayedNeutronClusters(std::map<double,double> cluster_cb, std::map<double,double> cluster_totalq); ///< function to identify delayed neutron clusters
 
  private:
 
@@ -48,6 +52,9 @@ class ClusterClassifiers: public Tool {
   std::map<int,unsigned long> pmtid_to_channelkey;
   std::map<unsigned long,int> channelkey_to_pmtid;
 
+  int prompt_muon_index;	//cluster index for prompt muon candidate
+  std::vector<int> delayed_neutron_index;	//cluster indices for delayed neutron candidates
+
   /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged.
   int verbosity;
   int v_error=0;

UserTools/ClusterFinder/ClusterFinder.cpp

@@ -542,7 +542,7 @@ bool ClusterFinder::Execute(){
 
   } else {
 
-    std::cout <<"ClusterFinder: RecoADCHits Store does not exist and is not read out"<<std::endl;
+    if (verbose > 0) std::cout <<"ClusterFinder: RecoADCHits Store does not exist and is not read out"<<std::endl;
   }
 
   return true;

UserTools/EventSelector/EventSelector.cpp

@@ -46,6 +46,8 @@ bool EventSelector::Initialise(std::string configfile, DataModel &data){
   m_variables.Get("SaveStatusToStore", fSaveStatusToStore);
   m_variables.Get("IsMC",fIsMC);
   m_variables.Get("RecoPDG",fRecoPDG);
+  m_variables.Get("TriggerExtendedWindow",fTriggerExtended);
+  m_variables.Get("BeamOK",fBeamOK);
   m_variables.Get("CutConfiguration",fCutConfigurationName);
 
   if (!fIsMC){fMCFVCut = false; fMCPMTVolCut = false; fMCMRDCut = false; fMCPiKCut = false; fMCIsMuonCut = false; fMCIsElectronCut = false; fMCIsSingleRingCut = false; fMCIsMultiRingCut = false; fMCProjectedMRDHit = false; fMCEnergyCut = false; fPromptTrigOnly = false;}
@@ -211,6 +213,12 @@ bool EventSelector::Execute(){
   m_data->Stores.at("RecoEvent")->Set("RecoPDGVector",cluster_reco_pdg);
   m_data->Stores.at("RecoEvent")->Set("PDG",fRecoPDG);
 
+  bool passExtendedCut = this->EventSelectionByTriggerExtended();
+  m_data->Stores.at("RecoEvent")->Set("TriggerExtended",passExtendedCut);
+
+  bool passBeamOKCut = this->EventSelectionByBeamOK();
+  m_data->Stores.at("RecoEvent")->Set("BeamOK",passBeamOKCut);
+
   // Fill the EventSelection mask for the cuts that are supposed to be applied
   if (fMCPiKCut){
     fEventApplied |= EventSelector::kFlagMCPiK; 
@@ -329,6 +337,16 @@ bool EventSelector::Execute(){
     if (!passTriggerCut) fEventFlagged |= EventSelector::kFlagTrigger;
   }
 
+  if (fTriggerExtended){
+    fEventApplied |= EventSelector::kFlagExtended;
+    if (!passExtendedCut) fEventFlagged |= EventSelector::kFlagExtended;
+  }
+
+  if (fBeamOK){
+    fEventApplied |= EventSelector::kFlagBeamOK;
+    if (!passBeamOKCut) fEventFlagged |= EventSelector::kFlagBeamOK;
+  }
+
   if (fRecoPDG != -1){
     fEventApplied |= EventSelector::kFlagRecoPDG;
     if (!passRecoPDGCut) fEventFlagged |= EventSelector::kFlagRecoPDG;
@@ -338,7 +356,10 @@ bool EventSelector::Execute(){
   if(fEventCutStatus){  
     Log("EventSelector Tool: Event is clean according to current event selection.",v_message,verbosity);
   }
-  if(fSaveStatusToStore) m_data->Stores.at("RecoEvent")->Set("EventCutStatus", fEventCutStatus);
+  if(fSaveStatusToStore) {
+    m_data->Stores.at("RecoEvent")->Set("EventCutStatus", fEventCutStatus);
+    m_data->Stores.at("ANNIEEvent")->Set("EventCutStatus", fEventCutStatus);
+  }
   m_data->Stores.at("RecoEvent")->Set("EventFlagApplied", fEventApplied);
   m_data->Stores.at("RecoEvent")->Set("EventFlagged", fEventFlagged);
 
@@ -349,6 +370,9 @@ bool EventSelector::Execute(){
 
   if (verbosity > 1) std::cout <<"EventCutStatus: "<<fEventCutStatus<<std::endl;
 
+  //std::cout << "EventSelector tool: Bit representation: fEventApplied: " << std::bitset<32>(fEventApplied) << ", fEventFlagged: " << std::bitset<32>(fEventFlagged) << std::endl;
+  //std::cout <<"EventCutStatus: "<<fEventCutStatus<<std::endl;
+
 
   return true;
 }
@@ -546,8 +570,7 @@ bool EventSelector::EventSelectionByMCTruthMRD() {
   double mrdEndZ = fGeometry->GetMrdEnd()*100-168.1;
   double mrdHeightY = fGeometry->GetMrdHeight()*100;
   double mrdWidthX = fGeometry->GetMrdWidth()*100;
-  std::cout <<"mrdStartZ: "<<mrdStartZ<<", mrdEndZ: "<<mrdEndZ<<", mrdHeightY: "<<mrdHeightY<<", mrdWidthX: "<<mrdWidthX<<std::endl;                                                                                     
-  Log("EventSelector tool: Read in MuonStop (X,Y,Z) = ("+std::to_string(muonStopX)+","+std::to_string(muonStopY)+","+std::to_string(muonStopZ)+")");
+  Log("EventSelector tool: Read in MuonStop (X,Y,Z) = ("+std::to_string(muonStopX)+","+std::to_string(muonStopY)+","+std::to_string(muonStopZ)+")",2,verbosity);
   if(muonStopZ<mrdStartZ || muonStopZ>mrdEndZ
   	|| muonStopX<-1.0*mrdWidthX || muonStopX>mrdWidthX
   	|| muonStopY<-1.0*mrdHeightY || muonStopY>mrdHeightY) {
@@ -626,10 +649,10 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
   m_data->Stores["RecoEvent"]->Set("NumPMTClusters",pmt_cluster_size);
   vec_pmtclusters_charge->clear();
   vec_pmtclusters_time->clear();
-  m_data->Stores["RecoEvent"]->Set("PMTClustersCharge",vec_pmtclusters_charge,false);
-  m_data->Stores["RecoEvent"]->Set("PMTClustersTime",vec_pmtclusters_time,false);
+  m_data->Stores["RecoEvent"]->Set("PMTClustersCharge",vec_pmtclusters_charge,true);
+  m_data->Stores["RecoEvent"]->Set("PMTClustersTime",vec_pmtclusters_time,true);
   vec_mrdclusters_time->clear();
-  m_data->Stores["RecoEvent"]->Set("MRDClustersTime",vec_mrdclusters_time);
+  m_data->Stores["RecoEvent"]->Set("MRDClustersTime",vec_mrdclusters_time,true);
   if (verbosity > 1) std::cout <<"pmt_cluster_size: "<<pmt_cluster_size<<", mrd cluster size: "<<MrdTimeClusters.size()<<std::endl;
 
   bool prompt_cluster = false;
@@ -694,8 +717,8 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
     std::cout <<"Number of PMT hits in muon cluster: "<<n_hits<<std::endl;
   }
 
-  m_data->Stores["RecoEvent"]->Set("PMTClustersCharge",vec_pmtclusters_charge,false);
-  m_data->Stores["RecoEvent"]->Set("PMTClustersTime",vec_pmtclusters_time,false);
+  m_data->Stores["RecoEvent"]->Set("PMTClustersCharge",vec_pmtclusters_charge,true);
+  m_data->Stores["RecoEvent"]->Set("PMTClustersTime",vec_pmtclusters_time,true);
 
   std::vector<double> mrd_meantimes;
   if (verbosity > 1) std::cout <<"MrdTimeClusters.size(): "<<MrdTimeClusters.size()<<std::endl;
@@ -725,7 +748,7 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
   for (int i=0; i<(int)mrd_meantimes.size(); i++){
     vec_mrdclusters_time->push_back(mrd_meantimes.at(i));
   }
-  m_data->Stores["RecoEvent"]->Set("MRDClustersTime",vec_mrdclusters_time);
+  m_data->Stores["RecoEvent"]->Set("MRDClustersTime",vec_mrdclusters_time, true);
 
   if (fIsMC){
     if (MrdTimeClusters.size() == 0 || m_all_clusters_MC->size() == 0) return false;
@@ -736,16 +759,23 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
   pmtmrd_coinc_min = fPMTMRDOffset - 50;
   pmtmrd_coinc_max = fPMTMRDOffset + 50;
 
+  std::vector<int> vector_mrd_coincidence;
+
   bool coincidence = false;
   for (int i_mrd = 0; i_mrd < int(mrd_meantimes.size()); i_mrd++){
     double time_diff = mrd_meantimes.at(i_mrd) - pmt_time;
     if (verbosity > 0) std::cout <<"MRD time: "<<mrd_meantimes.at(i_mrd)<<", PMT time: "<<pmt_time<<", difference: "<<time_diff<<std::endl;
     Log("EventSelector tool: MRD/Tank coincidene candidate "+std::to_string(i_mrd)+ " has time difference: "+std::to_string(time_diff),v_message,verbosity);
     if (verbosity > 1) std::cout <<"max_charge: "<<max_charge<<", n_hits: "<<n_hits<<std::endl;
+    Log("EventSelector tool: MRD/Tank coincidene candidate "+std::to_string(i_mrd)+ " has time difference: "+std::to_string(time_diff),1,verbosity);
+
     if (time_diff > pmtmrd_coinc_min && time_diff < pmtmrd_coinc_max && max_charge > 200 && n_hits >= 20){
       coincidence = true;
+      vector_mrd_coincidence.push_back(i_mrd);
     }
   }
+
+  m_data->Stores["RecoEvent"]->Set("MRDCoincidenceCluster",vector_mrd_coincidence);
 
   return coincidence;
 
@@ -1054,4 +1084,20 @@ bool EventSelector::EventSelectionByRecoPDG(int recoPDG, std::vector<double> & c
   return found_pdg;
 
 }
+bool EventSelector::EventSelectionByTriggerExtended(){
+
+  int fExtended = 0;
+  m_data->Stores["ANNIEEvent"]->Get("TriggerExtended",fExtended);
+  if (fExtended == 1) return true;
+  else return false;
+
+}
+
+bool EventSelector::EventSelectionByBeamOK(){
 
+  BeamStatus beamstat;
+  m_data->Stores["ANNIEEvent"]->Get("BeamStatus",beamstat);
+  if (beamstat.ok()) return true;
+  else return false;
+
+}

UserTools/EventSelector/EventSelector.h

@@ -19,6 +19,8 @@
 #include "ANNIEGeometry.h"
 #include "TMath.h"
 
+#include "BeamStatus.h"
+
 class EventSelector: public Tool {
 
 
@@ -51,7 +53,9 @@ class EventSelector: public Tool {
    kFlagVeto        = 0x20000, //131072
    kFlagTrigger      = 0x40000,
    kFlagThroughGoing     = 0x80000,
-   kFlagRecoPDG        = 0x1000000,
+   kFlagRecoPDG        = 0x100000,
+   kFlagExtended = 0x200000,
+   kFlagBeamOK        = 0x400000,
   } EventFlags_t;
 
  private:
@@ -176,6 +180,16 @@ class EventSelector: public Tool {
   //
   bool EventSelectionByRecoPDG(int recoPDG, std::vector<double> & vector_reco_pdg);
 
+  /// \brief Event selection for extended acquisition windows
+  //
+  /// This event selection criterion flags events with an extended trigger window (70us)
+  bool EventSelectionByTriggerExtended();
+
+  /// \brief Event selection for beam status
+  ///
+  /// This event selection criterion flags events for which beam status was ok
+  bool EventSelectionByBeamOK();
+
   /// \brief MC entry number
   uint64_t fMCEventNum;
 
@@ -237,6 +251,8 @@ class EventSelector: public Tool {
   bool fIsMC; 
   int fTriggerWord;
   int fRecoPDG;
+  bool fTriggerExtended = false;
+  bool fBeamOK = false;
   std::string fCutConfigurationName;  
 
   bool get_mrd = false;