/* This file is part of MAUS: http://micewww.pp.rl.ac.uk/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 "src/map/MapCppGlobalPID/MapCppGlobalPID.hh" #include "Interface/Squeak.hh" #include "src/common_cpp/DataStructure/Data.hh" #include "src/common_cpp/API/PyWrapMapBase.hh" #include "src/common_cpp/Converter/DataConverters/JsonCppSpillConverter.hh" #include "src/common_cpp/Converter/DataConverters/CppJsonSpillConverter.hh" namespace MAUS { PyMODINIT_FUNC init_MapCppGlobalPID(void) { PyWrapMapBase::PyWrapMapBaseModInit ("MapCppGlobalPID", "", "", "", ""); } MapCppGlobalPID::MapCppGlobalPID() : MapBase("MapCppGlobalPID"), _configCheck(false) { } void MapCppGlobalPID::_birth(const std::string& argJsonConfigDocument) { // Check if the JSON document can be parsed, else return error only. _configCheck = false; bool parsingSuccessful = _reader.parse(argJsonConfigDocument, _configJSON); if (!parsingSuccessful) { throw MAUS::Exception(Exception::recoverable, "Failed to parse Json configuration file", "MapCppGlobalPID::_birth"); } char* pMAUS_ROOT_DIR = getenv("MAUS_ROOT_DIR"); if (!pMAUS_ROOT_DIR) { throw MAUS::Exception(Exception::recoverable, std::string("Could not find the $MAUS_ROOT_DIR env variable. ")+\ std::string("Did you try running: source env.sh?"), "MapCppGlobalPID::_birth"); } _hypotheses.clear(); _pid_vars.clear(); PDF_file = _configJSON["PID_PDFs_file"].asString(); // PDF_file = "/home/celeste/MICE/MAUS/1389a/src/map/MapCppGlobalPID/PIDhists.root"; std::cerr << PDF_file << std::endl; _histFile = new TFile(PDF_file.c_str(), "READ"); // vector of hypotheses // TODO(Pidcott) find a more elegant way of accessing hypotheses _hypotheses.push_back("200MeV_mu_plus"); _hypotheses.push_back("200MeV_e_plus"); _hypotheses.push_back("200MeV_pi_plus"); for (unsigned int i =0; i < _hypotheses.size(); ++i) { // vector of pid vars _pid_vars.push_back(new MAUS::recon::global::PIDVarA(_histFile, _hypotheses[i])); _pid_vars.push_back(new MAUS::recon::global::PIDVarB(_histFile, _hypotheses[i])); _pid_vars.push_back(new MAUS::recon::global::PIDVarC(_histFile, _hypotheses[i])); // etc. } _configCheck = true; } void MapCppGlobalPID::_death() { } void MapCppGlobalPID::_process(MAUS::Data* data) const { MAUS::Data* data_cpp = data; if (!data_cpp) { throw Exception(Exception::recoverable, "Data was NULL", "MapCppGlobalPID::process"); } if (!_configCheck) { throw Exception(Exception::recoverable, "Birth was not called successfully", "MapCppGlobalPID::process"); } const MAUS::Spill* _spill = data_cpp->GetSpill(); if ( _spill->GetReconEvents() ) { for ( unsigned int event_i = 0; event_i < _spill->GetReconEvents()->size(); ++event_i ) { MAUS::GlobalEvent* global_event = _spill->GetReconEvents()->at(event_i)->GetGlobalEvent(); std::vector *GlobalTrackArray = global_event->get_tracks(); 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") continue; // doubles to hold cumulative log likelihoods for each hypothesis double logL_200MeV_mu_plus = 0; double logL_200MeV_e_plus = 0; double logL_200MeV_pi_plus = 0; for (size_t pid_var_count = 0; pid_var_count < _pid_vars.size(); ++pid_var_count) { std::string hyp = _pid_vars[pid_var_count]->Get_hyp(); if (hyp == "200MeV_mu_plus") { if (_pid_vars[pid_var_count]->logL(track) == 1) { continue; } else { logL_200MeV_mu_plus += _pid_vars[pid_var_count]->logL(track); } } else if (hyp == "200MeV_e_plus") { if (_pid_vars[pid_var_count]->logL(track) == 1) { continue; } else { logL_200MeV_e_plus += _pid_vars[pid_var_count]->logL(track); } } else if (hyp == "200MeV_pi_plus") { if (_pid_vars[pid_var_count]->logL(track) == 1) { continue; } else { logL_200MeV_pi_plus += _pid_vars[pid_var_count]->logL(track); } } else { Squeak::mout(Squeak::debug) << "Unrecognised particle hypothesis," << " MapCppGlobalPID::process" << std::endl; } } if ((logL_200MeV_mu_plus - logL_200MeV_e_plus > 0.5) && (logL_200MeV_mu_plus - logL_200MeV_pi_plus > 0.5)) { track->set_pid(MAUS::DataStructure::Global::kMuPlus); } else if ((logL_200MeV_e_plus - logL_200MeV_mu_plus > 0.5) && (logL_200MeV_e_plus - logL_200MeV_pi_plus > 0.5)) { track->set_pid(MAUS::DataStructure::Global::kEPlus); } else if ((logL_200MeV_pi_plus - logL_200MeV_mu_plus > 0.5) && (logL_200MeV_pi_plus - logL_200MeV_e_plus > 0.5)) { track->set_pid(MAUS::DataStructure::Global::kPiPlus); } else { Squeak::mout(Squeak::debug) << "PID for track could not be" << " determined." << std::endl; continue; } } } } } } // ~MAUS