PMTWaveformSim tool (re-make of PR #321)

UserTools/AssignBunchTimingMC/AssignBunchTimingMC.cpp

@@ -17,6 +17,7 @@ bool AssignBunchTimingMC::Initialise(std::string configfile, DataModel &data){
     bool got_count     = m_variables.Get("bunchcount", fbunchcount);
     bool got_sample    = m_variables.Get("sampletype", fsample);
     bool got_trigger   = m_variables.Get("prompttriggertime", ftriggertime);
+    bool got_waveform  = m_variables.Get("PMTWaveformSim", fPMTWaveformSim);
 
     if (!got_verbosity) {
         verbosity = 0;
@@ -62,6 +63,13 @@ bool AssignBunchTimingMC::Initialise(std::string configfile, DataModel &data){
         Log(logmessage, v_warning, verbosity);
     }
 
+    if (!got_waveform) {
+        fPMTWaveformSim = false;     // assume they are using the standard, parametric MC Hits
+        logmessage = ("Warning (AssignBunchTimingMC): \"PMTWaveformSim\" not "
+            "set in the config file. Using default MC Hits instead of waveform hits");
+        Log(logmessage, v_warning, verbosity);
+    }
+
 
     if (verbosity >= v_message) {
         std::cout<<"------------------------------------"<<"\n";
@@ -70,6 +78,7 @@ bool AssignBunchTimingMC::Initialise(std::string configfile, DataModel &data){
         std::cout<<"bunch width       = "<<fbunchwidth<<" ns"<<"\n";
         std::cout<<"bunch interval    = "<<fbunchinterval<<" ns"<<"\n";
         std::cout<<"number of bunches = "<<fbunchcount<<"\n";
+        std::cout<<"PMTWaveformSim    = "<<fPMTWaveformSim<<"\n";
         std::cout<<"sample type       = "<<(fsample == 0 ? "(0) Tank" : "(1) World")<<"\n";
         std::cout<<"trigger time      = "<<(ftriggertime == 0 ? "(0) prompt trigger starts when first particle arrives (default WCSim)" 
                                                    : "(1) prompt trigger starts at beam dump (modified WCSim)")<<"\n";
@@ -90,6 +99,16 @@ bool AssignBunchTimingMC::Execute()
         std::cout << "AssignBunchTimingMC: Executing tool..." << std::endl;
     }
 
+    // An upstream tool may opt to skip this execution stage
+    // For example the PMTWaveformSim tool will skip events with no MCHits or if
+    // no waveforms are produced.
+    bool skip = false;
+    bool got_skip_status = m_data->Stores["ANNIEEvent"]->Get("SkipExecute", skip);
+    if (got_skip_status && skip) {
+        Log("AssignBunchTimingMC: An upstream tool told me to skip this event.",v_warning,verbosity);
+        return true;
+    } 
+
     if (!LoadStores())      // Load info from store
         return false;
     if (verbosity >= v_debug) {
@@ -103,15 +122,35 @@ bool AssignBunchTimingMC::Execute()
         std::cout << "AssignBunchTimingMC: BNB timing successful" << std::endl;
     }
 
-    for (std::pair<double, std::vector<MCHit>>&& apair : *fClusterMapMC) {
-        double totalHitTime = 0;
-        int hitCount = 0;
-        int totalHits = apair.second.size();
+    if (fPMTWaveformSim) {   // PMTWaveformSim data-like clusters (read from ClusterMap)
+        if (verbosity >= v_debug) {
+            std::cout << "AssignBunchTimingMC: Reading from ClusterMap (PMTWaveformSim)" << std::endl;
+        }
+        for (std::pair<double, std::vector<Hit>>&& apair : *fClusterMap) {
+            double totalHitTime = 0;
+            int hitCount = 0;
+            int totalHits = apair.second.size();
+
+            CalculateClusterAndBunchTimesPMTWaveformSim(apair.second, totalHitTime, hitCount, totalHits);
+
+            // store the cluster time in a map (e.g., keyed by cluster identifier)
+            fbunchTimes->emplace(apair.first, bunchTime);
+        }
+
+    } else {   // otherwise default to reading the ClusterMapMC
+        if (verbosity >= v_debug) {
+            std::cout << "AssignBunchTimingMC: Reading from ClusterMapMC (default MCHits)" << std::endl;
+        }
+        for (std::pair<double, std::vector<MCHit>>&& apair : *fClusterMapMC) {
+            double totalHitTime = 0;
+            int hitCount = 0;
+            int totalHits = apair.second.size();
 
-        CalculateClusterAndBunchTimes(apair.second, totalHitTime, hitCount, totalHits);
+            CalculateClusterAndBunchTimes(apair.second, totalHitTime, hitCount, totalHits);
 
-        // store the cluster time in a map (e.g., keyed by cluster identifier)
-        fbunchTimes->emplace(apair.first, bunchTime);
+            // store the cluster time in a map (e.g., keyed by cluster identifier)
+            fbunchTimes->emplace(apair.first, bunchTime);
+        }
     }
 
     if (verbosity >= v_debug) {
@@ -142,23 +181,25 @@ bool AssignBunchTimingMC::LoadStores()
 {
     // grab necessary information from Stores
 
-    bool get_MCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC);
-    if (!get_MCClusters) {
-        Log("AssignBunchTimingMC: no ClusterMapMC in the CStore! Are you sure you ran the ClusterFinder tool?", v_error, verbosity);
-        return false;
+    if (fPMTWaveformSim) {
+        bool get_Clusters = m_data->CStore.Get("ClusterMap", fClusterMap);
+        if (!get_Clusters) {
+            Log("AssignBunchTimingMC: no ClusterMap in the CStore! Are you sure you ran the ClusterFinder (and the PMTWaveformSim) tool?", v_error, verbosity);
+            return false;
+        }
+    } else {
+        bool get_MCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC);
+        if (!get_MCClusters) {
+            Log("AssignBunchTimingMC: no ClusterMapMC in the CStore! Are you sure you ran the ClusterFinder tool?", v_error, verbosity);
+            return false;
+        }
     }
 
     bool get_AnnieEvent = m_data->Stores.count("ANNIEEvent");
     if (!get_AnnieEvent) {
         Log("AssignBunchTimingMC: no ANNIEEvent store!", v_error, verbosity);
         return false;
     }
-
-    bool get_MCHits = m_data->Stores.at("ANNIEEvent")->Get("MCHits", fMCHitsMap);
-    if (!get_MCHits) {
-        Log("AssignBunchTimingMC: no MCHits in the ANNIEEvent!", v_error, verbosity);
-        return false;
-    }
 
     bool get_neutrino_vtxt = m_data->Stores["GenieInfo"]->Get("NuIntxVtx_T",TrueNuIntxVtx_T);   // grab neutrino interaction time
     if (!get_neutrino_vtxt) {
@@ -177,7 +218,7 @@ void AssignBunchTimingMC::BNBtiming()
 {
     // Determined from GENIE samples (as of Dec 2024)
     const double tank_time = 33.0;    // Tank neutrino arrival time: 33ns
-    const double world_time = 33.0;   // WORLD neutrino arrival time: 33ns  (As of Dec 2024, World samples have not yet been re-produced fully)
+    const double world_time = 33.0;   // WORLD neutrino arrival time: 33ns
 
     if (ftriggertime == 0) {
         new_nu_time = (fsample == 0) ? (TrueNuIntxVtx_T - tank_time) : (TrueNuIntxVtx_T - world_time);
@@ -240,6 +281,37 @@ void AssignBunchTimingMC::CalculateClusterAndBunchTimes(std::vector<MCHit> const
 }
 
 
+void AssignBunchTimingMC::CalculateClusterAndBunchTimesPMTWaveformSim(std::vector<Hit> const &hits, double &totalHitTime, int &hitCount, int &totalHits)
+{
+
+    // loop over the hits to get their times
+    for (auto hit : hits) { 
+        double hitTime = hit.GetTime();
+        totalHitTime += hitTime;
+        hitCount++;
+        if (verbosity >= v_debug) {
+            std::string logmessage = "AssignBunchTimingMC: (" + std::to_string(hitCount) + "/" + std::to_string(totalHits) + ") PMTWaveformSim Hit time = " + std::to_string(hitTime);
+            Log(logmessage, v_debug, verbosity);
+        }
+    }
+
+    // find nominal cluster time (average hit time)
+    double clusterTime = (hitCount > 0) ? totalHitTime / hitCount : -9999;
+    if (verbosity >= v_debug) {
+        std::string logmessage = "AssignBunchTimingMC: ClusterTime = " + std::to_string(clusterTime);
+        Log(logmessage, v_debug, verbosity);
+    }
+
+    // calculate BunchTime
+    bunchTime = fbunchinterval * bunchNumber + clusterTime + jitter + new_nu_time ;
+    if (verbosity >= v_debug) {
+        std::string logmessage = "AssignBunchTimingMC: bunchTime = " + std::to_string(bunchTime);
+        Log(logmessage, v_debug, verbosity);
+    }
+
+}
+
+
 //------------------------------------------------------------------------------
 
 // done

UserTools/AssignBunchTimingMC/AssignBunchTimingMC.h

@@ -30,11 +30,12 @@ class AssignBunchTimingMC: public Tool {
         bool LoadStores(); ///< Loads all relevant information from the store, away from the Execute function
         void BNBtiming(); ///< Calculates the appropriate BNB timing to apply to the clusters
         void CalculateClusterAndBunchTimes(std::vector<MCHit> const &mchits, double &totalHitTime, int &hitCount, int &totalHits); ///< Loops through the MCHits, finds the cluster times (avg hit time), and calculates the new bunch timing
+        void CalculateClusterAndBunchTimesPMTWaveformSim(std::vector<Hit> const &hits, double &totalHitTime, int &hitCount, int &totalHits); ///< Same as above but for the data-like PMTWaveformSim clusters and hits
 
     private:
 
-        std::map<unsigned long, std::vector<MCHit>> *fMCHitsMap    = nullptr;  ///< All of the MCHits keyed by channel number
         std::map<double, std::vector<MCHit>>        *fClusterMapMC = nullptr;  ///< All MC clusters
+        std::map<double, std::vector<Hit>>          *fClusterMap = nullptr;    ///< All MC clusters (data-like hits from the PMTWaveformSim tool)
         double TrueNuIntxVtx_T;                                                ///< true neutrino interaction time in ns, from GenieInfo store
 
         std::map<double, double> *fbunchTimes                      = nullptr;  ///< Bunch-realistic timing from the cluster times; 
@@ -47,6 +48,7 @@ class AssignBunchTimingMC: public Tool {
         int fbunchcount;          ///< number of BNB bunches per spill
         int fsample;              ///< GENIE Tank or WORLD samples
         int ftriggertime;         ///< whether the samples used the default WCSim prompt trigger = 0 (when particles enter the volume), or the adjusted prompt trigger based on the start of the beam dump
+        bool fPMTWaveformSim;     ///< whether to use the PMTWaveform data-like hits or the defaul MCHits
 
         double new_nu_time;       ///< offset needed to make the cluster times beam realistic
         int bunchNumber;          ///< randomly assigned bunch number

UserTools/AssignBunchTimingMC/README.md

@@ -15,12 +15,14 @@
 
 verbosity 0
 
+PMTWaveformSim 1       # Whether we're using the data-like hits output from PMTWaveformSim
+
 # BNB properties taken from: MicroBooNE https://doi.org/10.1103/PhysRevD.108.052010
 bunchwidth 1.308       # BNB instrinic bunch spread [ns]
 bunchinterval 18.936   # BNB bunch spacings [ns]
 bunchcount 81          # number of BNB bunches per spill
 
-sampletype 0           # Tank (0) or World (1) genie samples you are running over
+sampletype 1           # Tank (0) or World (1) genie samples you are running over
 prompttriggertime 1    # WCSim prompt trigger settings: (0 = default, t0 = 0 when a particle enters the volume)
                        #                                (1 = modified, t0 = 0 when the neutrino beam dump begins)
 ```

UserTools/ClusterClassifiers/ClusterClassifiers.cpp

@@ -36,6 +36,14 @@ bool ClusterClassifiers::Initialise(std::string configfile, DataModel &data){
 bool ClusterClassifiers::Execute(){
 
   //We're gonna make ourselves a couple cluster classifier maps boyeeee
+  bool skip = false;
+  bool goodSkipStatus = m_data->Stores.at("ANNIEEvent")->Get("SkipExecute", skip);
+  if (goodSkipStatus && skip) {
+    logmessage = "ClusterClassifier: An upstream tool told me to skip this event.";
+    Log(logmessage, v_warning, verbosity);
+    return true;
+  }
+
   if(verbosity>4) std::cout << "ClusterClassifiers tool: Accessing cluster map in CStore" << std::endl;
   bool get_clusters = false;
   m_data->CStore.Get("ClusterMap",m_all_clusters);

UserTools/ClusterFinder/ClusterFinder.cpp

@@ -164,6 +164,16 @@ bool ClusterFinder::Execute(){
   //----------------------------------------------------------------------------
   //---------------get the members of the ANNIEEvent----------------------------
   //----------------------------------------------------------------------------
+
+  // An upstream tool may opt to skip this execution stage
+  // For example the PMTWaveformSim tool will skip events with no MCHits or if
+  // no waveforms are produced.
+  bool skip = false;
+  bool got_skip_status = m_data->Stores["ANNIEEvent"]->Get("SkipExecute", skip);
+  if (got_skip_status && skip) {
+    Log("ClusterFinder: An upstream tool told me to skip this event.",v_warning,verbose);
+    return true;
+  } 
 
   m_data->Stores["ANNIEEvent"]->Get("EventNumber", evnum);
   //m_data->Stores["ANNIEEvent"]->Get("BeamStatus", BeamStatus);

UserTools/Factory/Factory.cpp

@@ -174,6 +174,7 @@ if (tool=="BackTracker") ret=new BackTracker;
 if (tool=="PrintDQ") ret=new PrintDQ;
 if (tool=="AssignBunchTimingMC") ret=new AssignBunchTimingMC;
 if (tool=="FitRWMWaveform") ret=new FitRWMWaveform;
+if (tool=="PMTWaveformSim") ret=new PMTWaveformSim;
 if (tool=="LAPPDWaveformDisplay") ret=new LAPPDWaveformDisplay;
 return ret;
 }

UserTools/LoadGeometry/LoadGeometry.cpp

@@ -77,6 +77,7 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
   TankPMTCrateSpaceToChannelNumMap = new std::map<std::vector<int>,int>;
   ChannelNumToTankPMTSPEChargeMap = new std::map<int,double>;
   ChannelNumToTankPMTTimingOffsetMap = new std::map<unsigned long,double>;
+  ChannelNumToTankPMTTimingSigmaMap = new std::map<unsigned long,double>;
   ChannelNumToTankPMTCrateSpaceMap = new std::map<int,std::vector<int>>;
   AuxCrateSpaceToChannelNumMap = new std::map<std::vector<int>,int>;
   AuxChannelNumToTypeMap = new std::map<int,std::string>;
@@ -112,6 +113,7 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
   m_data->CStore.Set("ChannelNumToTankPMTCrateSpaceMap",ChannelNumToTankPMTCrateSpaceMap);
   m_data->CStore.Set("ChannelNumToTankPMTSPEChargeMap",ChannelNumToTankPMTSPEChargeMap);
   m_data->CStore.Set("ChannelNumToTankPMTTimingOffsetMap",ChannelNumToTankPMTTimingOffsetMap);
+  m_data->CStore.Set("ChannelNumToTankPMTTimingSigmaMap",ChannelNumToTankPMTTimingSigmaMap);
   m_data->CStore.Set("AuxCrateSpaceToChannelNumMap",AuxCrateSpaceToChannelNumMap);
   m_data->CStore.Set("AuxChannelNumToCrateSpaceMap",AuxChannelNumToCrateSpaceMap);
   m_data->CStore.Set("AuxChannelNumToTypeMap",AuxChannelNumToTypeMap);
@@ -919,21 +921,25 @@ void LoadGeometry::LoadTankPMTGains(){
   return;
 }
 
+// load in both the timing offsets and uncertainties from the 2023 laser campaign
 void LoadGeometry::LoadTankPMTTimingOffsets(){
   ifstream myfile(fTankPMTTimingOffsetFile.c_str());
   std::string line;
   if (myfile.is_open()){
-    //Loop over lines, collect all detector data (should only be one line here)
+    // Timing offset file has columns: [0] chankey [1] PMT_location [2] offset value (ns) [3] sigma (ns) [4] notes
     while(getline(myfile,line)){
       if(verbosity>3) std::cout << line << std::endl; //has our stuff;
       if(line.find("#")!=std::string::npos) continue;
       std::vector<std::string> DataEntries;
       boost::split(DataEntries,line, boost::is_any_of(","), boost::token_compress_on);
       int channelkey = -9999;
-      double TimingOffset = -9999.;
+      double TimingOffset = -9999;
+      double TimingSigma = -9999;
       channelkey = std::stoul(DataEntries.at(0));
       TimingOffset= std::stod(DataEntries.at(2));
+      TimingSigma = std::stod(DataEntries.at(3));
       ChannelNumToTankPMTTimingOffsetMap->emplace(channelkey,TimingOffset);
+      ChannelNumToTankPMTTimingSigmaMap->emplace(channelkey,TimingSigma);
     }
   }
   return;

UserTools/LoadGeometry/LoadGeometry.h

@@ -69,6 +69,7 @@ class LoadGeometry: public Tool {
   std::map<int,std::vector<int>>* AuxChannelNumToCrateSpaceMap;
   std::map<int,double>* ChannelNumToTankPMTSPEChargeMap;
   std::map<unsigned long,double>* ChannelNumToTankPMTTimingOffsetMap;
+  std::map<unsigned long,double>* ChannelNumToTankPMTTimingSigmaMap;
   std::map<int,std::string>* AuxChannelNumToTypeMap;
   std::map<std::vector<unsigned int>,int>* LAPPDCrateSpaceToChannelNumMap;