/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus * * MAUS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * MAUS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MAUS. If not, see . * */ #include #include #include "src/common_cpp/Utils/JsonWrapper.hh" #include "src/common_cpp/Utils/CppErrorHandler.hh" #include "Utils/Squeak.hh" #include "src/common_cpp/Utils/Exception.hh" #include "Interface/dataCards.hh" #include "src/common_cpp/DataStructure/Spill.hh" #include "src/common_cpp/DataStructure/ReconEvent.hh" #include "src/common_cpp/JsonCppProcessors/SpillProcessor.hh" #include "src/common_cpp/Converter/DataConverters/JsonCppSpillConverter.hh" #include "src/reduce/ReduceCppGlobalPID/ReduceCppGlobalPID.hh" namespace MAUS { ReduceCppGlobalPID::~ReduceCppGlobalPID() { } void ReduceCppGlobalPID::_birth(const std::string& argJsonConfigDocument) { _configCheck = false; _classname = "ReduceCppGlobalPID"; // JsonCpp setup - check file parses correctly, if not return false Json::Value _configJSON; try { _configJSON = JsonWrapper::StringToJson(argJsonConfigDocument); _pid_beam_setting = _configJSON["pid_beam_setting"].asString(); _unique_identifier = _configJSON["unique_identifier"].asString(); if (_pid_beam_setting.empty() || !_configJSON.isMember("pid_beam_setting")) { Squeak::mout(Squeak::error) << "Config did not contain a valid " << "pid_beam_setting, which is required for PDF production, " << "ReduceCppGlobalPID::birth" << std::endl; } if (_unique_identifier.empty() || !_configJSON.isMember("unique_identifier")) { Squeak::mout(Squeak::error) << "Config did not contain a valid " << "unique_identifier, which is required for PDF production, " << "ReduceCppGlobalPID::birth" << std::endl; } _pid_config = _configJSON["pid_config"].asString(); _pid_beamline_polarity = _configJSON["pid_beamline_polarity"].asString(); // vector of hypotheses if (_pid_beamline_polarity == "positive") { _hypotheses.push_back((_pid_beam_setting + "_mu_plus")); _hypotheses.push_back((_pid_beam_setting + "_e_plus")); _hypotheses.push_back((_pid_beam_setting + "_pi_plus")); } else if (_pid_beamline_polarity == "negative") { _hypotheses.push_back((_pid_beam_setting + "_mu_minus")); _hypotheses.push_back((_pid_beam_setting + "_e_minus")); _hypotheses.push_back((_pid_beam_setting + "_pi_minus")); } else { Squeak::mout(Squeak::warning) << "Invalid pid_beamline_polarity " << "set in ConfigurationDefaults, " << "ReduceCppGlobalPID::_birth" << std::endl; } if (_pid_config == "step_4") { // vector of mu pid vars _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarA(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarB(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarC(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarD(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarE(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarF(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarG(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarH(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarI(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::PIDVarJ(_hypotheses[0], _unique_identifier)); // vector of e pid vars _e_pid_vars.push_back(new MAUS::recon::global::PIDVarA(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarB(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarC(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarD(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarE(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarF(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarG(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarH(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarI(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::PIDVarJ(_hypotheses[1], _unique_identifier)); // vector of pi pid vars _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarA(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarB(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarC(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarD(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarE(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarF(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarG(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarH(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarI(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::PIDVarJ(_hypotheses[2], _unique_identifier)); } else if (_pid_config == "commissioning") { // vector of mu pid vars _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarA(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarB(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarC(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarD(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarE(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarF(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarG(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarH(_hypotheses[0], _unique_identifier)); _mu_pid_vars.push_back(new MAUS::recon::global::ComPIDVarI(_hypotheses[0], _unique_identifier)); // vector of e pid vars _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarA(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarB(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarC(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarD(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarE(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarF(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarG(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarH(_hypotheses[1], _unique_identifier)); _e_pid_vars.push_back(new MAUS::recon::global::ComPIDVarI(_hypotheses[1], _unique_identifier)); // vector of pi pid vars _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarA(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarB(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarC(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarD(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarE(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarF(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarG(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarH(_hypotheses[2], _unique_identifier)); _pi_pid_vars.push_back(new MAUS::recon::global::ComPIDVarI(_hypotheses[2], _unique_identifier)); } else { Squeak::mout(Squeak::warning) << "Invalid pid_config, " << " ReduceCppGlobalPID::birth" << std::endl; } _configCheck = true; } catch (Exceptions::Exception& exc) { MAUS::CppErrorHandler::getInstance()->HandleExceptionNoJson(exc, _classname); } catch (std::exception& exc) { MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(exc, _classname); } } void ReduceCppGlobalPID::_process(MAUS::Data* data_cpp) { if (data_cpp == NULL) throw Exceptions::Exception(Exceptions::recoverable, "Data was NULL", "ReduceCppMCProp::_process"); if (data_cpp->GetSpill() == NULL) throw Exceptions::Exception(Exceptions::recoverable, "Spill was NULL", "ReduceCppMCProp::_process"); if (data_cpp->GetSpill()->GetDaqEventType() != "physics_event") { } if (!_configCheck) { throw Exceptions::Exception(Exceptions::recoverable, "Birth was not called successfully", "ReduceCppGlobalPID::process"); } /*if (data_cpp->GetSpill()->GetReconEvents() == NULL) throw Exceptions::Exception(Exceptions::recoverable, "ReconEvents were NULL", "ReduceCppGlobalPID::_process");*/ _spill = data_cpp->GetSpill(); if (_spill) { if ( _spill->GetReconEvents() ) { for ( unsigned int event_i = 0; event_i < _spill->GetReconEvents()->size(); ++event_i) { if (_spill->GetReconEvents()->at(event_i)->GetGlobalEvent()) { MAUS::GlobalEvent* global_event = _spill->GetReconEvents()->at(event_i)->GetGlobalEvent(); std::vector *GlobalTrackArray = global_event->get_tracks(); bool through_track_check = false; for (unsigned int track_i = 0; track_i < GlobalTrackArray->size(); ++track_i) { MAUS::DataStructure::Global::Track* track = GlobalTrackArray->at(track_i); if (track->get_mapper_name() != "MapCppGlobalTrackMatching_Through") continue; through_track_check = true; if (_pid_config == "step_4") { // get constituent tracks std::vector const_tracks = track->GetConstituentTracks(); for (unsigned int track_j = 0; track_j < const_tracks.size(); ++track_j) { if (const_tracks.at(track_j)->get_mapper_name() == "MapCppGlobalTrackMatching_US") { MAUS::DataStructure::Global::Track* US_track = const_cast (const_tracks.at(track_j)); if ((_pid_beamline_polarity == "positive" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -13) || (_pid_beamline_polarity == "negative" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 13)) { for (size_t pid_var_count = 0; pid_var_count < _mu_pid_vars.size(); ++pid_var_count) { _mu_pid_vars[pid_var_count]->Fill_Hist(US_track); } // loop over US mu_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -11) || (_pid_beamline_polarity == "negative" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 11)) { for (size_t pid_var_count = 0; pid_var_count < _e_pid_vars.size(); ++pid_var_count) { _e_pid_vars[pid_var_count]->Fill_Hist(US_track); } // loop over US e_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 211) || (_pid_beamline_polarity == "negative" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -211)) { for (size_t pid_var_count = 0; pid_var_count < _pi_pid_vars.size(); ++pid_var_count) { _pi_pid_vars[pid_var_count]->Fill_Hist(US_track); } // loop over US pi_pid_vars } // US mc pid check } else if (const_tracks.at(track_j)->get_mapper_name() == "MapCppGlobalTrackMatching_DS") { MAUS::DataStructure::Global::Track* DS_track = const_cast (const_tracks.at(track_j)); if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -13) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 13)) { for (size_t pid_var_count = 0; pid_var_count < _mu_pid_vars.size(); ++pid_var_count) { _mu_pid_vars[pid_var_count]->Fill_Hist(DS_track); } // loop over DS mu_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -11) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 11)) { for (size_t pid_var_count = 0; pid_var_count < _e_pid_vars.size(); ++pid_var_count) { _e_pid_vars[pid_var_count]->Fill_Hist(DS_track); } // loop over DS e_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 211) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -211)) { for (size_t pid_var_count = 0; pid_var_count < _pi_pid_vars.size(); ++pid_var_count) { _pi_pid_vars[pid_var_count]->Fill_Hist(DS_track); } // loop over DS pi_pid_vars } // DS mc pid check } // get DS TM track } // loop over constituent tracks } else if (_pid_config == "commissioning") { if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -13) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 13)) { for (size_t pid_var_count = 0; pid_var_count < _mu_pid_vars.size(); ++pid_var_count) { _mu_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over mu_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -11) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 11)) { for (size_t pid_var_count = 0; pid_var_count < _e_pid_vars.size(); ++pid_var_count) { _e_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over e_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 211) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -211)) { for (size_t pid_var_count = 0; pid_var_count < _pi_pid_vars.size(); ++pid_var_count) { _pi_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over pi_pid_vars } // check straight through track mc pid at DS tracker ref plane } // change behaviour based on straight or helical tracks } // loop over global tracks if (through_track_check == false) { for (unsigned int track_i = 0; track_i < GlobalTrackArray->size(); ++track_i) { MAUS::DataStructure::Global::Track* track = GlobalTrackArray->at(track_i); if (track->get_mapper_name() == "MapCppGlobalTrackMatching_US") { if ((_pid_beamline_polarity == "positive" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -13) || (_pid_beamline_polarity == "negative" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 13)) { for (size_t pid_var_count = 0; pid_var_count < _mu_pid_vars.size(); ++pid_var_count) { _mu_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over US mu_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -11) || (_pid_beamline_polarity == "negative" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 11)) { for (size_t pid_var_count = 0; pid_var_count < _e_pid_vars.size(); ++pid_var_count) { _e_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over US e_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 211) || (_pid_beamline_polarity == "negative" && _mc_pid_US_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -211)) { for (size_t pid_var_count = 0; pid_var_count < _pi_pid_vars.size(); ++pid_var_count) { _pi_pid_vars[pid_var_count]->Fill_Hist(track); } } } // US tracks if (track->get_mapper_name() == "MapCppGlobalTrackMatching_DS") { if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -13) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 13)) { for (size_t pid_var_count = 0; pid_var_count < _mu_pid_vars.size(); ++pid_var_count) { _mu_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over US mu_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -11) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 11)) { for (size_t pid_var_count = 0; pid_var_count < _e_pid_vars.size(); ++pid_var_count) { _e_pid_vars[pid_var_count]->Fill_Hist(track); } // loop over US e_pid_vars } else if ((_pid_beamline_polarity == "positive" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == 211) || (_pid_beamline_polarity == "negative" && _mc_pid_DS_tracker_ref(_spill->GetMCEvents()->at(event_i)) == -211)) { for (size_t pid_var_count = 0; pid_var_count < _pi_pid_vars.size(); ++pid_var_count) { _pi_pid_vars[pid_var_count]->Fill_Hist(track); } } } // DS tracks } } // temporary measure for when through tracks haven't been formed } // get global event } // get recon event } // get recon event array } else { Squeak::mout(Squeak::error) << "Failed to import spill from data\n"; } // check for spill } void ReduceCppGlobalPID::_death() { if (_configCheck) { // _pid_vars.clear(); for (size_t pid_var_count = 0; pid_var_count < _mu_pid_vars.size(); ++pid_var_count) { delete _mu_pid_vars[pid_var_count]; } for (size_t pid_var_count = 0; pid_var_count < _e_pid_vars.size(); ++pid_var_count) { delete _e_pid_vars[pid_var_count]; } for (size_t pid_var_count = 0; pid_var_count < _pi_pid_vars.size(); ++pid_var_count) { delete _pi_pid_vars[pid_var_count]; } } } int ReduceCppGlobalPID::_mc_pid_US_tracker_ref(MAUS::MCEvent* mc_event) { int us_pid = 0; MAUS::SciFiHitArray* scifihits = mc_event->GetSciFiHits(); std::vector::iterator scifihit; for ( size_t s = 0; s < scifihits->size(); ++s ) { MAUS::SciFiHit scifihit = scifihits->at(s); if (scifihit.GetChannelId()->GetTrackerNumber() == 0 && scifihit.GetChannelId()->GetStationNumber() == 1 && scifihit.GetChannelId()->GetPlaneNumber() == 0) { us_pid = scifihit.GetParticleId(); break; } } return us_pid; } int ReduceCppGlobalPID::_mc_pid_DS_tracker_ref(MAUS::MCEvent* mc_event) { int ds_pid = 0; MAUS::SciFiHitArray* scifihits = mc_event->GetSciFiHits(); std::vector::iterator scifihit; for ( size_t s = 0; s < scifihits->size(); ++s ) { MAUS::SciFiHit scifihit = scifihits->at(s); if (scifihit.GetChannelId()->GetTrackerNumber() == 1 && scifihit.GetChannelId()->GetStationNumber() == 1 && scifihit.GetChannelId()->GetPlaneNumber() == 0) { ds_pid = scifihit.GetParticleId(); break; } } return ds_pid; } PyMODINIT_FUNC init_ReduceCppGlobalPID(void) { PyWrapReduceBase::PyWrapReduceBaseModInit( "ReduceCppGlobalPID", "", "", "", ""); } } // ~namespace MAUS