Skip to content

EventSelector modifications for new MC waveform hits #330

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marc1uk merged 5 commits into from
Nov 12, 2025
Merged

EventSelector modifications for new MC waveform hits #330

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e7807c3
Adjust EventSelector logic to handle MC waveforms
S81D Jan 6, 2025
a3c1d8c
Update README.md
S81D Jan 6, 2025
a7773fa
Update EventSelector.cpp
S81D Sep 8, 2025
3b1b36a
Update EventSelector.cpp
S81D Oct 22, 2025
eab30b8
Update EventSelector.cpp
S81D Nov 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 105 additions & 93 deletions UserTools/EventSelector/EventSelector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ bool EventSelector::Initialise(std::string configfile, DataModel &data){

fPMTMRDOffset = false;
fIsMC = true;
fMCWaveform = false;
fPMTMRDOffset = 755;
fRecoPDG = -1;

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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);
Expand All @@ -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();
Expand All @@ -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<double,std::vector<Hit>>&& apair : *m_all_clusters){
std::vector<Hit>&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<double,std::vector<MCHit>>&& apair : *m_all_clusters_MC){
Expand All @@ -690,38 +731,6 @@ bool EventSelector::EventSelectionByPMTMRDCoinc() {
}
}
}
} else {
if (m_all_clusters->size()){
double cluster_time;
for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
std::vector<Hit>&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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -1026,12 +1039,12 @@ bool EventSelector::FindPaddleChankey(double x, double y, int layer, unsigned lo

bool EventSelector::EventSelectionByRecoPDG(int recoPDG, std::vector<double> & 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<double,double> ClusterMaxPEs;
Expand All @@ -1045,52 +1058,51 @@ bool EventSelector::EventSelectionByRecoPDG(int recoPDG, std::vector<double> & c
bool found_pdg = false;

if (fabs(recoPDG)==2112){
if (fIsMC){
if (m_all_clusters_MC->size()){
for(std::pair<double,std::vector<MCHit>>&& apair : *m_all_clusters_MC){
double cluster_time = apair.first;
double charge_balance = ClusterChargeBalances.at(cluster_time);
std::vector<MCHit>&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 = "<<cluster_time<<", with CB = "<<charge_balance <<" and charge "<<charge_temp<<"!!!"<<std::endl;
} else {
std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
}
}
}
} else {
if (m_all_clusters->size()){
double cluster_time;
for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
double cluster_time = apair.first;
double charge_balance = ClusterChargeBalances.at(cluster_time);
std::vector<Hit>&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 = "<<cluster_time<<", with CB = "<<charge_balance <<"and charge "<<charge_temp<<"!!!"<<std::endl;
} else {
std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
if (!fIsMC || fMCWaveform) { // data-like
if (m_all_clusters->size()){
for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
double cluster_time = apair.first;
double charge_balance = ClusterChargeBalances.at(cluster_time);
std::vector<Hit>&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 = "<<cluster_time<<", with CB = "<<charge_balance <<"and charge "<<charge_temp<<"!!!"<<std::endl;
} else {
std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
}
}
}
} else { // MCHits
if (m_all_clusters_MC->size()){
for(std::pair<double,std::vector<MCHit>>&& apair : *m_all_clusters_MC){
double cluster_time = apair.first;
double charge_balance = ClusterChargeBalances.at(cluster_time);
std::vector<MCHit>&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 = "<<cluster_time<<", with CB = "<<charge_balance <<" and charge "<<charge_temp<<"!!!"<<std::endl;
} else {
std::cout <<"Did not pass neutron candidate cuts!!! Time = "<<cluster_time <<", CB = "<<charge_balance << " and charge "<<charge_temp<<"!!!"<<std::endl;
}
}
}
}
}
}
}

return found_pdg;
Expand Down
1 change: 1 addition & 0 deletions UserTools/EventSelector/EventSelector.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,7 @@ class EventSelector: public Tool {
bool fThroughGoing = false;
bool fEventCutStatus;
bool fIsMC;
bool fMCWaveform;
int fTriggerWord;
int fRecoPDG;
bool fTriggerExtended = false;
Expand Down
1 change: 1 addition & 0 deletions UserTools/EventSelector/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,5 +106,6 @@ ThroughGoing
TriggerWord
RecoPDG
IsMC
PMTWaveformSim
SaveStatusToStore
```