PMT timing offsets, new hit timing method, and changes to ClusterFinder to fix charge omission seen for MC Hits

UserTools/ClusterClassifiers/ClusterClassifiers.cpp

@@ -222,8 +222,8 @@ double ClusterClassifiers::CalculateChargeBalance(std::vector<Hit> cluster_hits)
     total_Q+= tube_charge;
     total_QSquared += (tube_charge * tube_charge);
   }
-  //FIXME: Need a method to have the 123 be equal to the number of operating detectors
-  double charge_balance  = sqrt((total_QSquared)/(total_Q*total_Q) - (1./123.));
+  //FIXME: Need a method to have the 1/N be equal to the number of operating detectors
+  double charge_balance  = sqrt((total_QSquared)/(total_Q*total_Q) - (1./121.));
   if(verbosity>4) std::cout << "ClusterClassifiers Tool: Calculated charge balance of " << charge_balance << std::endl;
   return charge_balance;
 }
@@ -249,8 +249,8 @@ double ClusterClassifiers::CalculateChargeBalanceMC(std::vector<MCHit> cluster_h
      total_Q+= tube_charge;
      total_QSquared += (tube_charge * tube_charge);
    }
-   //FIXME: Need a method to have the 123 be equal to the number of operating detectors
-   double charge_balance  = sqrt((total_QSquared)/(total_Q*total_Q) - (1./123.));
+   //FIXME: Need a method to have the 1/N be equal to the number of operating detectors
+   double charge_balance  = sqrt((total_QSquared)/(total_Q*total_Q) - (1./121.));
    if(verbosity>4) std::cout << "ClusterClassifiers Tool: Calculated charge balance of " << charge_balance << std::endl;
    return charge_balance;
  }

UserTools/ClusterFinder/ClusterFinder.cpp

@@ -211,39 +211,84 @@ bool ClusterFinder::Execute(){
         std::vector<double> datalike_hits;
         std::vector<double> datalike_hits_charge;
         for (MCHit &ahit : ThisPMTHits){
-          //std::cout <<"Key: "<<detectorkey<<", charge "<<ahit.GetCharge()<<", time "<<ahit.GetTime()<<std::endl;
-          //if (ahit.GetTime() > 2000.) std::cout <<"Found hit later than 2us! Hit time : "<<ahit.GetTime()<<", chankey: "<<chankey<<std::endl;
           if (ahit.GetTime() < end_of_window_time_cut*AcqTimeWindow) {
             //Make MC more like data --> combine multiple photons if they are within a 10ns range
-            //hit times can only be recorded with 2ns precision --> possible times are 0ns, 2ns, 4ns, ...
-            hits_2ns_res.push_back(2*(int(ahit.GetTime())/2.)+(int(ahit.GetTime())%2));
+            hits_2ns_res.push_back(ahit.GetTime());
             hits_2ns_res_charge.push_back(ahit.GetCharge());
           }
         }
 	//Combine multiple MC hits to one pulse
         std::sort(hits_2ns_res.begin(),hits_2ns_res.end());
-	for (int i_hit=0; i_hit < (int) hits_2ns_res.size(); i_hit++){
-          double hit1 = hits_2ns_res.at(i_hit);
-          if (datalike_hits.size()==0) {
-            datalike_hits.push_back(hit1);
-            datalike_hits_charge.push_back(hits_2ns_res_charge.at(i_hit));
-          }
-          else {
+	std::vector<double> temp_times;
+	double temp_charges = 0.0;
+	double mid_time;
+	if (verbose > 0){
+		std::cout << " " << std::endl;
+		std::cout << hits_2ns_res.size() << " total photon hits(s)" << std::endl;
+	}
+
+	if (hits_2ns_res.size() > 0){
+          for (int i_hit=0; i_hit < (int) hits_2ns_res.size(); i_hit++){
+            double hit1 = hits_2ns_res.at(i_hit);
+            if (temp_times.size()==0) {
+              temp_times.push_back(hit1);
+              temp_charges += hits_2ns_res_charge.at(i_hit);
+            }
+
+	else {
             bool new_pulse = false;
-            for (int j_hit=0; j_hit < (int) datalike_hits.size(); j_hit++){
-              if (fabs(datalike_hits.at(j_hit)-hit1)<10.) {
+        	if (fabs(temp_times[0]-hit1)<10.) {
                 new_pulse=false;
-                datalike_hits_charge.at(j_hit)+=hits_2ns_res_charge.at(i_hit);
-                break;
+                temp_charges+=hits_2ns_res_charge.at(i_hit);
+		temp_times.push_back(hit1);
               } else new_pulse=true;
-            }
+		  
             if (new_pulse) {
-              datalike_hits.push_back(hit1);      //Only count as a new pulse if it was 10ns away from every other pulse
-              datalike_hits_charge.push_back(hits_2ns_res_charge.at(i_hit));
+		// following the DigitBuilder tool --> take median photon hit time as the hit time of the "pulse"
+		if (temp_times.size() % 2 == 0){
+                  mid_time = (temp_times.at(temp_times.size()/2 - 1) + temp_times.at(temp_times.size()/2))/2;
+                } else{
+                  mid_time = temp_times.at(temp_times.size()/2);
+                }
+
+              datalike_hits.push_back(mid_time);
+              datalike_hits_charge.push_back(temp_charges);
+              temp_times.clear();
+              temp_charges = 0;
+              temp_times.push_back(hit1);
+              temp_charges += hits_2ns_res_charge.at(i_hit);
+
             }
           }
         }
 
+	if (temp_times.size() % 2 == 0){
+            mid_time = (temp_times.at(temp_times.size()/2 - 1) + temp_times.at(temp_times.size()/2))/2;
+        } else{
+            mid_time = temp_times.at(temp_times.size()/2);
+        }
+
+        datalike_hits.push_back(mid_time);
+        datalike_hits_charge.push_back(temp_charges);
+
+        if (verbose > 0){
+          std::cout << " " << std::endl;
+          std::cout << datalike_hits.size() << " MC pulse(s) identified from the raw photon hit(s)" << std::endl;
+          std::cout << "Pulse time(s):" << std::endl;
+          for (int ih=0; ih < (int) datalike_hits.size(); ih++){
+            double junk = datalike_hits.at(ih);
+            std::cout << junk << " ";
+          }
+          std::cout << " " << std::endl;
+          std::cout << "Pulse charge(s):" << std::endl;
+          for (int ih=0; ih < (int) datalike_hits_charge.size(); ih++){
+            double junk2 = datalike_hits_charge.at(ih);
+            std::cout << junk2 << " ";
+          }
+          std::cout << " " << std::endl;
+        }
+        }
+
         for (int i_hit = 0; i_hit < (int) datalike_hits.size(); i_hit++){
           //v_hittimes.push_back(ahit.GetTime()); // fill a vector with all hit times (unsorted)
           v_hittimes.push_back(datalike_hits.at(i_hit));
@@ -295,7 +340,7 @@ bool ClusterFinder::Execute(){
 
   // Now sort the hit time array, fill the highest time in a new array until the old array is empty
   do {
-    double max_time = 0;
+    double max_time = -9999;
     int i_max_time = 0;
     for (std::vector<double>::iterator it = v_hittimes.begin(); it != v_hittimes.end(); ++it) {
       if (*it > max_time) {
@@ -321,12 +366,21 @@ bool ClusterFinder::Execute(){
     }
     thiswindow_Nhits = 0;   
     v_mini_hits.clear();
-    for (double j_time = *it; j_time < *it + ClusterFindingWindow; j_time+=2){  // loops through times in the window and check if there's a hit at this time
+    for (double j_time = *it; j_time < *it + ClusterFindingWindow; j_time+=1){  // loops through times in the window and check if there's a hit at this time
       for(std::vector<double>::iterator it2 = v_hittimes_sorted.begin(); it2 != v_hittimes_sorted.end(); ++it2) {
-        if (*it2 == j_time) {
-          thiswindow_Nhits++;
-          v_mini_hits.push_back(*it2);
+        if(HitStoreName=="MCHits"){
+          if (static_cast<int>(j_time) == static_cast<int>(*it2)) {     // accept all hit times (some may be smeared to negative values)
+            thiswindow_Nhits++;
+            v_mini_hits.push_back(*it2);
+          }
         }
+        if(HitStoreName=="Hits"){
+          if (*it2 > 0 && static_cast<int>(j_time) == static_cast<int>(*it2)) {  // reject hit times in the data that are 0
+            thiswindow_Nhits++;
+            v_mini_hits.push_back(*it2);
+          }
+        }
+
       }
     }
     if (!v_mini_hits.empty()) {

UserTools/EventSelector/EventSelector.cpp

@@ -14,7 +14,7 @@ bool EventSelector::Initialise(std::string configfile, DataModel &data){
 
   fPMTMRDOffset = false;
   fIsMC = true;
-  fPMTMRDOffset = 745;
+  fPMTMRDOffset = 755;
   fRecoPDG = -1;
 
   //Get the tool configuration variables
@@ -695,8 +695,15 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
         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;
+          // 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);
@@ -769,7 +776,7 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
     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){
+    if (time_diff > pmtmrd_coinc_min && time_diff < pmtmrd_coinc_max){
       coincidence = true;
       vector_mrd_coincidence.push_back(i_mrd);
     }

UserTools/LoadGeometry/LoadGeometry.cpp

@@ -20,6 +20,7 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
   m_variables.Get("FACCMRDGeoFile", fFACCMRDGeoFile);
   m_variables.Get("TankPMTGeoFile", fTankPMTGeoFile);
   m_variables.Get("TankPMTGainFile", fTankPMTGainFile);
+  m_variables.Get("TankPMTTimingOffsetFile", fTankPMTTimingOffsetFile);
   m_variables.Get("AuxiliaryChannelFile", fAuxChannelFile);
   m_variables.Get("LAPPDGeoFile", fLAPPDGeoFile);
   m_variables.Get("DetectorGeoFile", fDetectorGeoFile);
@@ -58,6 +59,12 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
     return false;
   }
 
+  if(!this->FileExists(fTankPMTTimingOffsetFile)){
+    Log("LoadGeometry Tool: File for Tank PMT Timing offsets does not exist!",v_error,verbosity);
+    if (verbosity > 0) std::cout << "Filepath was... " << fTankPMTTimingOffsetFile << std::endl;
+    return false;
+  }
+
   if(!this->FileExists(fAuxChannelFile)){
     Log("LoadGeometry Tool: File for Auxiliary Channels does not exist!",v_error,verbosity);
     if (verbosity > 0) std::cout << "Filepath was... " << fAuxChannelFile << std::endl;
@@ -69,6 +76,7 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
   MRDChannelNumToCrateSpaceMap = new std::map<int,std::vector<int>>;
   TankPMTCrateSpaceToChannelNumMap = new std::map<std::vector<int>,int>;
   ChannelNumToTankPMTSPEChargeMap = new std::map<int,double>;
+  ChannelNumToTankPMTTimingOffsetMap = 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>;
@@ -87,6 +95,9 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
   //Load TankPMT charge to PE conversion
   this->LoadTankPMTGains();
 
+  //Load TankPMT timing offsets
+  this->LoadTankPMTTimingOffsets();
+
   //Load auxiliary and spare channels
   this->LoadAuxiliaryChannels();
 
@@ -100,6 +111,7 @@ bool LoadGeometry::Initialise(std::string configfile, DataModel &data){
   m_data->CStore.Set("TankPMTCrateSpaceToChannelNumMap",TankPMTCrateSpaceToChannelNumMap);
   m_data->CStore.Set("ChannelNumToTankPMTCrateSpaceMap",ChannelNumToTankPMTCrateSpaceMap);
   m_data->CStore.Set("ChannelNumToTankPMTSPEChargeMap",ChannelNumToTankPMTSPEChargeMap);
+  m_data->CStore.Set("ChannelNumToTankPMTTimingOffsetMap",ChannelNumToTankPMTTimingOffsetMap);
   m_data->CStore.Set("AuxCrateSpaceToChannelNumMap",AuxCrateSpaceToChannelNumMap);
   m_data->CStore.Set("AuxChannelNumToCrateSpaceMap",AuxChannelNumToCrateSpaceMap);
   m_data->CStore.Set("AuxChannelNumToTypeMap",AuxChannelNumToTypeMap);
@@ -906,3 +918,23 @@ void LoadGeometry::LoadTankPMTGains(){
   }
   return;
 }
+
+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)
+    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.;
+      channelkey = std::stoul(DataEntries.at(0));
+      TimingOffset= std::stod(DataEntries.at(2));
+      ChannelNumToTankPMTTimingOffsetMap->emplace(channelkey,TimingOffset);
+    }
+  }
+  return;
+}

UserTools/LoadGeometry/LoadGeometry.h

@@ -36,6 +36,7 @@ class LoadGeometry: public Tool {
   bool ParseAuxChannelDataEntry(std::vector<std::string> SpecLine,
                                std::vector<std::string> AuxChannelLegendEntries);
   void LoadTankPMTGains();
+  void LoadTankPMTTimingOffsets();
 
   Geometry* AnnieGeometry;
 
@@ -49,6 +50,7 @@ class LoadGeometry: public Tool {
   std::string fFACCMRDGeoFile;
   std::string fTankPMTGeoFile;
   std::string fTankPMTGainFile;
+  std::string fTankPMTTimingOffsetFile; 
   std::string fAuxChannelFile;
   std::string fLAPPDGeoFile;
   std::string fDetectorGeoFile;
@@ -66,6 +68,7 @@ class LoadGeometry: public Tool {
   std::map<int,std::vector<int>>* ChannelNumToTankPMTCrateSpaceMap;
   std::map<int,std::vector<int>>* AuxChannelNumToCrateSpaceMap;
   std::map<int,double>* ChannelNumToTankPMTSPEChargeMap;
+  std::map<unsigned long,double>* ChannelNumToTankPMTTimingOffsetMap;
   std::map<int,std::string>* AuxChannelNumToTypeMap;
   std::map<std::vector<unsigned int>,int>* LAPPDCrateSpaceToChannelNumMap;
 

UserTools/LoadWCSim/LoadWCSim.cpp

@@ -418,7 +418,7 @@ bool LoadWCSim::Execute(){
 			int genieentry = firsttrigt->GetHeader()->GetGenieEntryNum();
 			if(verbosity>1) cout<<"Genie file is "<<geniefilename<<", genie event num was "<<genieentry<<endl;
 			m_data->CStore.Set("GenieFile",geniefilename);
-			m_data->CStore.Set("GenieEntry",std::to_string(genieentry));
+			m_data->CStore.Set("GenieEntry",genieentry);
 
 			for(int trigi=0; trigi<WCSimEntry->wcsimrootevent->GetNumberOfEvents(); trigi++){