PMT Waveform Simulation

DataModel/ADCPulse.cpp

@@ -7,9 +7,9 @@
 ADCPulse::ADCPulse(int TubeId, double start_time, double peak_time,
   double baseline, double sigma_baseline, unsigned long area,
   unsigned short raw_amplitude, double calibrated_amplitude,
-  double charge) : Hit(TubeId, start_time, charge),
+  double charge, double stop_time) : Hit(TubeId, start_time, charge),
   start_time_(start_time), peak_time_(peak_time),
   baseline_(baseline), sigma_baseline_(sigma_baseline), raw_area_(area),
-  raw_amplitude_(raw_amplitude), calibrated_amplitude_(calibrated_amplitude)
+  raw_amplitude_(raw_amplitude), calibrated_amplitude_(calibrated_amplitude), stop_time_(stop_time)
 {
 }

DataModel/ADCPulse.h

@@ -23,7 +23,7 @@ class ADCPulse : public Hit {
     ADCPulse(int TubeId, double start_time, double peak_time,
       double baseline, double sigma_baseline, unsigned long raw_area,
       unsigned short raw_amplitude, double calibrated_amplitude,
-      double charge);
+      double charge, double stop_time = 0);
 
     // @brief Returns the start time (ns) of the pulse relative to the
     // start of its minibuffer
@@ -33,6 +33,10 @@ class ADCPulse : public Hit {
     // start of its minibuffer
     inline double peak_time() const { return peak_time_; }
 
+    // @brief Returns the stop time (ns) of the pulse relative to the
+    // start of its minibuffer
+    inline double stop_time() const { return stop_time_; }
+
     // @brief Returns the approximate baseline (ADC) used to calibrate the
     // pulse
     inline double baseline() const { return baseline_; }
@@ -71,16 +75,24 @@ class ADCPulse : public Hit {
       ar & raw_area_;
       ar & raw_amplitude_;
       ar & calibrated_amplitude_;
+      if (version > 0)
+	ar & stop_time_;
     }
 
   protected:
 
     double start_time_; // ns since beginning of minibuffer
     double peak_time_; // ns since beginning of minibuffer
+    double stop_time_; // ns since beginning of minibuffer
     double baseline_; // mean (ADC)
     double sigma_baseline_; // standard deviation (ADC)
     unsigned long raw_area_; // (ADC * samples)
 
     unsigned short raw_amplitude_; // ADC
     double calibrated_amplitude_; // V
+
 };
+
+// Need to increment the class version since we added time as a new variable
+// the version number ensures backward compatibility when serializing 
+BOOST_CLASS_VERSION(ADCPulse, 1)

UserTools/BackTracker/BackTracker.cpp

@@ -27,6 +27,13 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){
     Log(logmessage, v_error, verbosity);
   }
 
+  bool gotMCWaveforms = m_variables.Get("MCWaveforms", fMCWaveforms);
+  if (!gotMCWaveforms) {
+    fMCWaveforms = false;
+    logmessage = "BackTracker::Initialize: \"MCWaveforms\" not set in the config, defaulting to false";
+    Log(logmessage, v_error, verbosity);
+  }
+
 
   // Set up the pointers we're going to save. No need to 
   // delete them at Finalize, the store will handle it
@@ -42,35 +49,91 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){
 //------------------------------------------------------------------------------
 bool BackTracker::Execute()
 {
-  if (!LoadFromStores())
-    return false;
-
+  int load_status = LoadFromStores();
+  if (load_status == 0) return false;
+  if (load_status == 2) return true;
+
+
   fClusterToBestParticleID ->clear();
   fClusterToBestParticlePDG->clear();
   fClusterEfficiency       ->clear();
   fClusterPurity           ->clear();
   fClusterTotalCharge      ->clear();
 
   fParticleToTankTotalCharge.clear();
+
   SumParticleTankCharge();
 
-
-  // Loop over the clusters and do the things
-  for (std::pair<double, std::vector<MCHit>>&& apair : *fClusterMapMC) {
-    int prtId = -5;
-    int prtPdg = -5;
-    double eff = -5;
-    double pur = -5;
-    double totalCharge = 0;
-
-    MatchMCParticle(apair.second, prtId, prtPdg, eff, pur, totalCharge);
-
-    fClusterToBestParticleID ->emplace(apair.first, prtId);
-    fClusterToBestParticlePDG->emplace(apair.first, prtPdg);
-    fClusterEfficiency       ->emplace(apair.first, eff);
-    fClusterPurity           ->emplace(apair.first, pur);
-    fClusterTotalCharge      ->emplace(apair.first, totalCharge);
-  }
+  if (fMCWaveforms) { // using clusters of Hits made from simulated PMT pulses
+
+    // Produce the map from channel ID to pulse time to MCHit index
+    //std::cout << "BT: Calling MapPulsesToParentIdxs" << std::endl;
+    bool gotPulseMap = MapPulsesToParentIdxs();
+    if (!gotPulseMap) {
+      logmessage = "BackTracker: No good pulse map.";
+      Log(logmessage, v_error, verbosity);
+      return false;
+    }
+
+    // Loop over the clusters
+    for (std::pair<double, std::vector<Hit>>&& apair : *fClusterMap) {
+      // Create a vector of MCHits associated with the vector of Hits
+      std::vector<MCHit> mcHits;
+      for (auto hit : apair.second) {
+	int channel_key = hit.GetTubeId();
+	double hitTime = hit.GetTime();
+
+	// Catches if something goes wrong
+	// First make sure that we have the PMT in the outer map
+	// Then make sure we have the hit time in the inner map
+	if (fMapChannelToPulseTimeToMCHitIdx.find(channel_key) == fMapChannelToPulseTimeToMCHitIdx.end()) {
+	  std::cout << "BackTracker: No hit on this PMT: " << channel_key << std::endl;
+	  return false;
+	}
+	if (fMapChannelToPulseTimeToMCHitIdx.at(channel_key).find(hitTime) == fMapChannelToPulseTimeToMCHitIdx.at(channel_key).end()) {
+	  std::cout << "BackTracker: No hit on this PMT: " << channel_key << " at this time: " << hitTime << std::endl;
+	  return false;
+	}
+
+	std::vector<int> mcHitIdxVec = fMapChannelToPulseTimeToMCHitIdx[channel_key][hitTime];
+	for (auto mcHitIdx : mcHitIdxVec) 
+	  mcHits.push_back((fMCHitsMap->at(channel_key)).at(mcHitIdx));
+      }// end loop over cluster hits
+
+      int prtId = -5;
+      int prtPdg = -5;
+      double eff = -5;
+      double pur = -5;
+      double totalCharge = 0;
+
+      MatchMCParticle(mcHits, prtId, prtPdg, eff, pur, totalCharge);
+
+      fClusterToBestParticleID ->emplace(apair.first, prtId);
+      fClusterToBestParticlePDG->emplace(apair.first, prtPdg);
+      fClusterEfficiency       ->emplace(apair.first, eff);
+      fClusterPurity           ->emplace(apair.first, pur);
+      fClusterTotalCharge      ->emplace(apair.first, totalCharge);
+
+      m_data->Stores.at("ANNIEEvent")->Set("MapChannelToPulseTimeToMCHitIdx",  fMapChannelToPulseTimeToMCHitIdx );
+    }// end loop over the cluster map
+  } else { // using clusters of MCHits
+    // Loop over the MC clusters and do the things
+    for (std::pair<double, std::vector<MCHit>>&& apair : *fClusterMapMC) {
+      int prtId = -5;
+      int prtPdg = -5;
+      double eff = -5;
+      double pur = -5;
+      double totalCharge = 0;
+
+      MatchMCParticle(apair.second, prtId, prtPdg, eff, pur, totalCharge);
+
+      fClusterToBestParticleID ->emplace(apair.first, prtId);
+      fClusterToBestParticlePDG->emplace(apair.first, prtPdg);
+      fClusterEfficiency       ->emplace(apair.first, eff);
+      fClusterPurity           ->emplace(apair.first, pur);
+      fClusterTotalCharge      ->emplace(apair.first, totalCharge);
+    }// end loop over the MC cluster map
+  }// end if/else fMCWaveforms
 
   m_data->Stores.at("ANNIEEvent")->Set("ClusterToBestParticleID",  fClusterToBestParticleID );
   m_data->Stores.at("ANNIEEvent")->Set("ClusterToBestParticlePDG", fClusterToBestParticlePDG);
@@ -91,8 +154,8 @@ bool BackTracker::Finalise()
 //------------------------------------------------------------------------------
 void BackTracker::SumParticleTankCharge()
 {
-  for (auto mcHitsIt : *fMCHitsMap) {
-    std::vector<MCHit> mcHits = mcHitsIt.second;
+  for (auto apair : *fMCHitsMap) {
+    std::vector<MCHit> mcHits = apair.second;
     for (uint mcHitIdx = 0; mcHitIdx < mcHits.size(); ++mcHitIdx) {
 
       // technically a MCHit could have multiple parents, but they don't appear to in practice
@@ -101,8 +164,8 @@ void BackTracker::SumParticleTankCharge()
       if (parentIdxs.size() != 1) continue;
 
       int particleId = -5;
-      for (auto it : *fMCParticleIndexMap) {
-	if (it.second == parentIdxs[0]) particleId = it.first;
+      for (auto bpair : *fMCParticleIndexMap) {
+	if (bpair.second == parentIdxs[0]) particleId = bpair.first;
       }
       if (particleId == -5) continue;
 
@@ -133,8 +196,8 @@ void BackTracker::MatchMCParticle(std::vector<MCHit> const &mchits, int &prtId,
     }
 
     int particleId = -5;
-    for (auto it : *fMCParticleIndexMap) {
-      if (it.second == parentIdxs[0]) particleId = it.first;
+    for (auto apair : *fMCParticleIndexMap) {
+      if (apair.second == parentIdxs[0]) particleId = apair.first;
     }
     if (particleId == -5) continue;
 
@@ -175,39 +238,129 @@ void BackTracker::MatchMCParticle(std::vector<MCHit> const &mchits, int &prtId,
 }
 
 //------------------------------------------------------------------------------
-bool BackTracker::LoadFromStores()
+bool BackTracker::MapPulsesToParentIdxs()
 {
-  // Grab the stuff we need from the stores
-  bool goodMCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC);
-  if (!goodMCClusters) {
-    std::cerr<<"BackTracker: no ClusterMapMC in the CStore!"<<endl;
+  // Clear out the map
+  fMapChannelToPulseTimeToMCHitIdx.clear();
+
+  // Grab the pulses
+  std::map<unsigned long, std::vector< std::vector<ADCPulse>> > adcPulseMap;
+  bool goodADCPulses = m_data->Stores.at("ANNIEEvent")->Get("RecoADCHits", adcPulseMap);
+  if (!goodADCPulses) {
+    logmessage = "BackTracker: no RecoADCHits in the ANNIEEvent!";
+    Log(logmessage, v_error, verbosity);
     return false;
   }
 
+  // Loop over the ADCPulses and find MCHits that fall within the start and stop times
+  // also record the pulse time to match the MCHits to the reco Hits
+  // Pulses are indexed by PMT id then stacked in a two-deep vector
+
+  for (auto apair: adcPulseMap) {
+    int channel_key = apair.first;
+
+    std::map<double, std::vector<int>> mapHitTimeToParents;
+    bool goodPulses = false;
+    for (auto pulseVec : apair.second) {
+      for (auto pulse : pulseVec) {
+	goodPulses = true;
+	double pulseTime = pulse.peak_time();
+
+	// Record the hit index if it occurred within the pulse window
+	// If there is only one hit the no need to check the times
+	std::vector<MCHit> mcHits = fMCHitsMap->at(channel_key);
+	if (mcHits.size() == 1)
+	  mapHitTimeToParents[pulseTime].push_back(0);
+	else {
+	  for (uint mcHitIdx = 0; mcHitIdx < mcHits.size(); ++mcHitIdx) {
+	    double hitTime = mcHits[mcHitIdx].GetTime();
+	    // The hit finding has to contend with noise so allow for a 10 ns 
+	    // slew in the pulse start time (I know it seems large)
+	    if ( hitTime + 10 >= pulse.start_time() && hitTime <= pulse.stop_time())
+	      mapHitTimeToParents[pulseTime].push_back(mcHitIdx);
+	  }// end loop over MCHits
+	}
+      }// end loop over inner pulse vector
+    }// end loop over outer pulse vector
+
+    if (mapHitTimeToParents.size() == 0 && goodPulses) {
+      logmessage = "BackTracker::MapPulsesToParentIdxs: No MCHits match with this pulse! PMT channel: ";
+      logmessage += std::to_string(channel_key);
+      Log(logmessage, v_error, verbosity);
+      return false;
+    }
+
+    fMapChannelToPulseTimeToMCHitIdx.emplace(channel_key, std::move(mapHitTimeToParents));
+  }// end loop over pulse map
+
+  return true;
+}
+
+//------------------------------------------------------------------------------
+int BackTracker::LoadFromStores()
+{
+  // Grab the stuff we need from the stores
+
   bool goodAnnieEvent = m_data->Stores.count("ANNIEEvent");
   if (!goodAnnieEvent) {
-    std::cerr<<"BackTracker: no ANNIEEvent store!"<<endl;
-    return false;
+    logmessage = "BackTracker: no ANNIEEvent store!";
+    Log(logmessage, v_error, verbosity);
+    return 0;
   }
-
+
+  bool skip = false;
+  bool goodSkipStatus = m_data->Stores.at("ANNIEEvent")->Get("SkipExecute", skip);
+  if (goodSkipStatus && skip) {
+    logmessage = "BackTracker: An upstream tool told me to skip this event.";
+    Log(logmessage, v_warning, verbosity);
+    return 2;
+  }
+
   bool goodMCHits = m_data->Stores.at("ANNIEEvent")->Get("MCHits", fMCHitsMap);
   if (!goodMCHits) {
-    std::cerr<<"BackTracker: no MCHits in the ANNIEEvent!"<<endl;
-    return false;
+    logmessage = "BackTracker: no MCHits in the ANNIEEvent!";
+    Log(logmessage, v_error, verbosity);
+    return 0;
   }
+
+  if (fMCWaveforms) {
+    bool goodClusters = m_data->CStore.Get("ClusterMap", fClusterMap);
+    if (!goodClusters) {
+      logmessage = "BackTracker: no ClusterMap in the CStore!";
+      Log(logmessage, v_error, verbosity);
+      return 0;
+    }
+
+    bool goodRecoHits = m_data->Stores.at("ANNIEEvent")->Get("Hits", fRecoHitsMap);
+    if (!goodRecoHits) {
+      logmessage = "BackTracker: no Hits in the ANNIEEvent!";
+      Log(logmessage, v_error, verbosity);
+      return 0;
+    }
+
+  } else {
+    bool goodMCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC);
+    if (!goodMCClusters) {
+      logmessage = "BackTracker: no ClusterMapMC in the CStore!";
+      Log(logmessage, v_error, verbosity);
+      return 0;
+    }
+  }// end if/else fMCWaveforms
 
   bool goodMCParticles = m_data->Stores.at("ANNIEEvent")->Get("MCParticles", fMCParticles);
   if (!goodMCParticles) {
-    std::cerr<<"BackTracker: no MCParticles in the ANNIEEvent!"<<endl;
-    return false;
+    logmessage = "BackTracker: no MCParticles in the ANNIEEvent!";
+    Log(logmessage, v_error, verbosity);
+    return 0;
   }
 
   bool goodMCParticleIndexMap = m_data->Stores.at("ANNIEEvent")->Get("TrackId_to_MCParticleIndex", fMCParticleIndexMap);
   if (!goodMCParticleIndexMap) {
-    std::cerr<<"BackTracker: no TrackId_to_MCParticleIndex in the ANNIEEvent!"<<endl;
-    return false;
+    logmessage = "BackTracker: no TrackId_to_MCParticleIndex in the ANNIEEvent!";
+    Log(logmessage, v_error, verbosity);
+    return 0;
   }
 
-  return true;
+  return 1;
 }