diff --git a/UserTools/EventSelector/EventSelector.cpp b/UserTools/EventSelector/EventSelector.cpp index 97efba7a6..e06fb99c5 100644 --- a/UserTools/EventSelector/EventSelector.cpp +++ b/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>&& apair : *m_all_clusters){ + std::vector&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>&& apair : *m_all_clusters_MC){ @@ -690,38 +731,6 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() { } } } - } else { - if (m_all_clusters->size()){ - double cluster_time; - for(std::pair>&& apair : *m_all_clusters){ - std::vector&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 & 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 ClusterMaxPEs; @@ -1045,52 +1058,51 @@ bool EventSelector::EventSelectionByRecoPDG(int recoPDG, std::vector & c bool found_pdg = false; if (fabs(recoPDG)==2112){ - if (fIsMC){ - if (m_all_clusters_MC->size()){ - for(std::pair>&& apair : *m_all_clusters_MC){ - double cluster_time = apair.first; - double charge_balance = ClusterChargeBalances.at(cluster_time); - std::vector&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 = "<size()){ - double cluster_time; - for(std::pair>&& apair : *m_all_clusters){ - double cluster_time = apair.first; - double charge_balance = ClusterChargeBalances.at(cluster_time); - std::vector&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 = "<size()){ + for(std::pair>&& apair : *m_all_clusters){ + double cluster_time = apair.first; + double charge_balance = ClusterChargeBalances.at(cluster_time); + std::vector&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 = "<size()){ + for(std::pair>&& apair : *m_all_clusters_MC){ + double cluster_time = apair.first; + double charge_balance = ClusterChargeBalances.at(cluster_time); + std::vector&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 = "<