/* 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) {
}
MapCppGlobalPID::~MapCppGlobalPID() {
}
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();
_histFile = new TFile(PDF_file.c_str(), "READ");
Json::Value pid_bounds = _configJSON["pid_bounds"];
_minA = pid_bounds["minA"].asDouble();
_maxA = pid_bounds["maxA"].asDouble();
_XminB = pid_bounds["XminB"].asDouble();
_XmaxB = pid_bounds["XmaxB"].asDouble();
_YminB = pid_bounds["YminB"].asDouble();
_YmaxB = pid_bounds["YmaxB"].asDouble();
_XminC = pid_bounds["XminC"].asDouble();
_XmaxC = pid_bounds["XmaxC"].asDouble();
_YminC = pid_bounds["YminC"].asDouble();
_YmaxC = pid_bounds["YmaxC"].asDouble();
_minD = pid_bounds["minD"].asDouble();
_maxD = pid_bounds["maxD"].asDouble();
_minE = pid_bounds["minE"].asDouble();
_maxE = pid_bounds["maxE"].asDouble();
_XminF = pid_bounds["XminF"].asDouble();
_XmaxF = pid_bounds["XmaxF"].asDouble();
_YminF = pid_bounds["YminF"].asDouble();
_YmaxF = pid_bounds["YmaxF"].asDouble();
_minG = pid_bounds["minG"].asDouble();
_maxG = pid_bounds["maxG"].asDouble();
_XminH = pid_bounds["XminH"].asDouble();
_XmaxH = pid_bounds["XmaxH"].asDouble();
_YminH = pid_bounds["YminH"].asDouble();
_YmaxH = pid_bounds["YmaxH"].asDouble();
_XminI = pid_bounds["XminI"].asDouble();
_XmaxI = pid_bounds["XmaxI"].asDouble();
_YminI = pid_bounds["YminI"].asDouble();
_YmaxI = pid_bounds["YmaxI"].asDouble();
_XminJ = pid_bounds["XminJ"].asDouble();
_XmaxJ = pid_bounds["XmaxJ"].asDouble();
_YminJ = pid_bounds["YminJ"].asDouble();
_YmaxJ = pid_bounds["YmaxJ"].asDouble();
_minComA = pid_bounds["minComA"].asDouble();
_maxComA = pid_bounds["maxComA"].asDouble();
_XminComB = pid_bounds["XminComB"].asDouble();
_XmaxComB = pid_bounds["XmaxComB"].asDouble();
_YminComB = pid_bounds["YminComB"].asDouble();
_YmaxComB = pid_bounds["YmaxComB"].asDouble();
_minComC = pid_bounds["minComC"].asDouble();
_maxComC = pid_bounds["maxComC"].asDouble();
_minComD = pid_bounds["minComD"].asDouble();
_maxComD = pid_bounds["maxComD"].asDouble();
_XminComE = pid_bounds["XminComE"].asDouble();
_XmaxComE = pid_bounds["XmaxComE"].asDouble();
_YminComE = pid_bounds["YminComE"].asDouble();
_YmaxComE = pid_bounds["YmaxComE"].asDouble();
_minComF = pid_bounds["minComF"].asDouble();
_maxComF = pid_bounds["maxComF"].asDouble();
_XminComG = pid_bounds["XminComG"].asDouble();
_XmaxComG = pid_bounds["XmaxComG"].asDouble();
_YminComG = pid_bounds["YminComG"].asDouble();
_YmaxComG = pid_bounds["YmaxComG"].asDouble();
_XminComH = pid_bounds["XminComH"].asDouble();
_XmaxComH = pid_bounds["XmaxComH"].asDouble();
_YminComH = pid_bounds["YminComH"].asDouble();
_YmaxComH = pid_bounds["YmaxComH"].asDouble();
_XminComI = pid_bounds["XminComI"].asDouble();
_XmaxComI = pid_bounds["XmaxComI"].asDouble();
_YminComI = pid_bounds["YminComI"].asDouble();
_YmaxComI = pid_bounds["YmaxComI"].asDouble();
_pid_config = _configJSON["pid_config"].asString();
_pid_mode = _configJSON["pid_mode"].asString();
_custom_pid_set = _configJSON["custom_pid_set"].asString();
_pid_confidence_level = _configJSON["pid_confidence_level"].asInt();
_pid_beamline_polarity = _configJSON["pid_beamline_polarity"].asString();
_pid_beam_setting = _configJSON["pid_beam_setting"].asString();
_pid_track_selection = _configJSON["pid_track_selection"].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, "
<< "MapCppGlobalPID::_birth" << std::endl;
}
if (_pid_config == "step_4") {
if (_pid_mode == "online") {
for (unsigned int i =0; i < _hypotheses.size(); ++i) {
// vector of pid vars
_pid_vars.push_back(new MAUS::recon::global::PIDVarB(_histFile,
_hypotheses[i],
_XminB, _XmaxB,
_YminB, _YmaxB));
_pid_vars.push_back(new MAUS::recon::global::PIDVarC(_histFile,
_hypotheses[i],
_XminC, _XmaxC,
_YminC, _YmaxC));
_pid_vars.push_back(new MAUS::recon::global::PIDVarF(_histFile,
_hypotheses[i],
_XminF, _XmaxF,
_YminF, _YmaxF));
_pid_vars.push_back(new MAUS::recon::global::PIDVarH(_histFile,
_hypotheses[i],
_XminH, _XmaxH,
_YminH, _YmaxH));
_pid_vars.push_back(new MAUS::recon::global::PIDVarI(_histFile,
_hypotheses[i],
_XminI, _XmaxI,
_YminI, _YmaxI));
_pid_vars.push_back(new MAUS::recon::global::PIDVarJ(_histFile,
_hypotheses[i],
_XminJ, _XmaxJ,
_YminJ, _YmaxJ));
}
} else if (_pid_mode == "offline") {
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],
_minA, _maxA));
_pid_vars.push_back(new MAUS::recon::global::PIDVarB(_histFile,
_hypotheses[i],
_XminB, _XmaxB,
_YminB, _YmaxB));
_pid_vars.push_back(new MAUS::recon::global::PIDVarC(_histFile,
_hypotheses[i],
_XminC, _XmaxC,
_YminC, _YmaxC));
_pid_vars.push_back(new MAUS::recon::global::PIDVarD(_histFile,
_hypotheses[i],
_minD, _maxD));
_pid_vars.push_back(new MAUS::recon::global::PIDVarE(_histFile,
_hypotheses[i],
_minE, _maxE));
_pid_vars.push_back(new MAUS::recon::global::PIDVarF(_histFile,
_hypotheses[i],
_XminF, _XmaxF,
_YminF, _YmaxF));
_pid_vars.push_back(new MAUS::recon::global::PIDVarG(_histFile,
_hypotheses[i],
_minG, _maxG));
_pid_vars.push_back(new MAUS::recon::global::PIDVarH(_histFile,
_hypotheses[i],
_XminH, _XmaxH,
_YminH, _YmaxH));
_pid_vars.push_back(new MAUS::recon::global::PIDVarI(_histFile,
_hypotheses[i],
_XminI, _XmaxI,
_YminI, _YmaxI));
_pid_vars.push_back(new MAUS::recon::global::PIDVarJ(_histFile,
_hypotheses[i],
_XminJ, _XmaxJ,
_YminJ, _YmaxJ));
}
} else if (_pid_mode == "custom") {
std::istringstream ss(_custom_pid_set);
std::string token;
std::vector input_pid_vars;
while(std::getline(ss, token, ' ')) {
input_pid_vars.push_back(token);
}
if (input_pid_vars.size() == 0) {
Squeak::mout(Squeak::warning) << "No PID variables in custom_pid_set,"
<< " MapCppGlobalPID::birth" << std::endl;
} else {
for (unsigned int i =0; i < _hypotheses.size(); ++i) {
// vector of pid vars
if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarA") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarA(_histFile,
_hypotheses[i],
_minA, _maxA));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarB") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarB(_histFile,
_hypotheses[i],
_XminB, _XmaxB,
_YminB, _YmaxB));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarC") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarC(_histFile,
_hypotheses[i],
_XminC, _XmaxC,
_YminC, _YmaxC));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarD") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarD(_histFile,
_hypotheses[i],
_minD, _maxD));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarE") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarE(_histFile,
_hypotheses[i],
_minE, _maxE));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarF") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarF(_histFile,
_hypotheses[i],
_XminF, _XmaxF,
_YminF, _YmaxF));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarG") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarG(_histFile,
_hypotheses[i],
_minG, _maxG));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarH") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarH(_histFile,
_hypotheses[i],
_XminH, _XmaxH,
_YminH, _YmaxH));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarI") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarI(_histFile,
_hypotheses[i],
_XminI, _XmaxI,
_YminI, _YmaxI));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"PIDVarJ") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::PIDVarJ(_histFile,
_hypotheses[i],
_XminJ, _XmaxJ,
_YminJ, _YmaxJ));
} else {
Squeak::mout(Squeak::warning) << "No valid PID variables given in "
<< "custom_pid_set, "
<< "MapCppGlobalPID::birth" << std::endl;
}
}
}
} else {
Squeak::mout(Squeak::warning) << "Invalid pid_mode, "
<< " MapCppGlobalPID::birth" << std::endl;
}
} else if (_pid_config == "commissioning") {
if (_pid_mode == "online") {
for (unsigned int i =0; i < _hypotheses.size(); ++i) {
// vector of pid vars
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarB(_histFile,
_hypotheses[i],
_XminComB, _XmaxComB,
_YminComB, _YmaxComB));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarE(_histFile,
_hypotheses[i],
_XminComE, _XmaxComE,
_YminComE, _YmaxComE));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarG(_histFile,
_hypotheses[i],
_XminComG, _XmaxComG,
_YminComG, _YmaxComG));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarH(_histFile,
_hypotheses[i],
_XminComH, _XmaxComH,
_YminComH, _YmaxComH));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarI(_histFile,
_hypotheses[i],
_XminComI, _XmaxComI,
_YminComI, _YmaxComI));
}
} else if (_pid_mode == "offline") {
for (unsigned int i =0; i < _hypotheses.size(); ++i) {
// vector of pid vars
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarA(_histFile,
_hypotheses[i],
_minComA, _maxComA));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarB(_histFile,
_hypotheses[i],
_XminComB, _XmaxComB,
_YminComB, _YmaxComB));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarC(_histFile,
_hypotheses[i],
_minComC, _maxComC));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarD(_histFile,
_hypotheses[i],
_minComD, _maxComD));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarF(_histFile,
_hypotheses[i],
_minComF, _maxComF));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarE(_histFile,
_hypotheses[i],
_XminComE, _XmaxComE,
_YminComE, _YmaxComE));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarG(_histFile,
_hypotheses[i],
_XminComG, _XmaxComG,
_YminComG, _YmaxComG));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarH(_histFile,
_hypotheses[i],
_XminComH, _XmaxComH,
_YminComH, _YmaxComH));
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarI(_histFile,
_hypotheses[i],
_XminComI, _XmaxComI,
_YminComI, _YmaxComI));
}
} else if (_pid_mode == "custom") {
std::istringstream ss(_custom_pid_set);
std::string token;
std::vector input_pid_vars;
while(std::getline(ss, token, ' ')) {
input_pid_vars.push_back(token);
}
if (input_pid_vars.size() == 0) {
Squeak::mout(Squeak::warning) << "No PID variables in custom_pid_set,"
<< " MapCppGlobalPID::birth" << std::endl;
} else {
for (unsigned int i =0; i < _hypotheses.size(); ++i) {
// vector of pid vars
if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarA") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarA(_histFile,
_hypotheses[i],
_minComA, _maxComA));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarB") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarB(_histFile,
_hypotheses[i],
_XminComB, _XmaxComB,
_YminComB, _YmaxComB));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarC") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarC(_histFile,
_hypotheses[i],
_minComC, _maxComC));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarD") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarD(_histFile,
_hypotheses[i],
_minComD, _maxComD));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarE") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarE(_histFile,
_hypotheses[i],
_XminComE, _XmaxComE,
_YminComE, _YmaxComE));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarF") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarF(_histFile,
_hypotheses[i],
_minComF, _maxComF));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarG") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarG(_histFile,
_hypotheses[i],
_XminComG, _XmaxComG,
_YminComG, _YmaxComG));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarH") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarH(_histFile,
_hypotheses[i],
_XminComH, _XmaxComH,
_YminComH, _YmaxComH));
} else if (std::find(input_pid_vars.begin(), input_pid_vars.end(),
"ComPIDVarI") != input_pid_vars.end()) {
_pid_vars.push_back(new MAUS::recon::global::ComPIDVarI(_histFile,
_hypotheses[i],
_XminComI, _XmaxComI,
_YminComI, _YmaxComI));
} else {
Squeak::mout(Squeak::warning) << "No valid PID variables given in "
<< "custom_pid_set, "
<< "MapCppGlobalPID::birth" << std::endl;
}
}
}
} else {
Squeak::mout(Squeak::warning) << "Invalid pid_mode, "
<< " MapCppGlobalPID::birth" << std::endl;
}
} else {
Squeak::mout(Squeak::warning) << "Invalid pid_config, "
<< " MapCppGlobalPID::birth" << std::endl;
}
_configCheck = true;
}
void MapCppGlobalPID::_death() {
_pid_vars.clear();
}
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 ) {
if (_spill->GetReconEvents()->at(event_i)->GetGlobalEvent()) {
MAUS::GlobalEvent* global_event =
_spill->GetReconEvents()->at(event_i)->GetGlobalEvent();
std::vector *GlobalTrackArray =
global_event->get_tracks();
std::istringstream track_selector(_pid_track_selection);
std::string token;
std::vector track_types;
while (std::getline(track_selector, token, ' ')) {
track_types.push_back(token);
}
if (track_types.size() == 0) {
Squeak::mout(Squeak::warning) << "No tracks have been selected to be ID'd,"
<< " MapCppGlobalPID::_process" << std::endl;
}
for (unsigned int q = 0; q < track_types.size(); ++q) {
if (track_types[q] == "all") {
US_PID(GlobalTrackArray, global_event);
DS_PID(GlobalTrackArray, global_event);
Through_US_PID(GlobalTrackArray, global_event);
Through_DS_PID(GlobalTrackArray, global_event);
Through_PID(GlobalTrackArray, global_event);
} else if (track_types[q] == "US") {
US_PID(GlobalTrackArray, global_event);
} else if (track_types[q] == "DS") {
DS_PID(GlobalTrackArray, global_event);
} else if (track_types[q] == "Through_US") {
Through_US_PID(GlobalTrackArray, global_event);
} else if (track_types[q] == "Through_DS") {
Through_DS_PID(GlobalTrackArray, global_event);
} else if (track_types[q] == "Through") {
Through_PID(GlobalTrackArray, global_event);
}
}
}
}
}
}
double MapCppGlobalPID::ConfidenceLevel(double LL_x, double sum_L) const {
return (exp(LL_x)/sum_L)*100;
}
void MapCppGlobalPID::US_PID(std::vector*
GlobalTrackArray,
MAUS::GlobalEvent* global_event) const {
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") {
perform_pid("MapCppGlobalPID-Final_US",
"US_PID", track, _pid_vars,
_pid_beamline_polarity,
_pid_beam_setting, global_event);
}
}
}
void MapCppGlobalPID::DS_PID(std::vector*
GlobalTrackArray,
MAUS::GlobalEvent* global_event) const {
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-DS") {
perform_pid("MapCppGlobalPID-Final_DS",
"DS_PID", track, _pid_vars,
_pid_beamline_polarity,
_pid_beam_setting, global_event);
}
}
}
void MapCppGlobalPID::Through_US_PID(std::vector*
GlobalTrackArray,
MAUS::GlobalEvent* global_event) const {
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") {
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));
perform_pid("MapCppGlobalPID-Final_Through_US",
"Through_US_PID", US_track, _pid_vars,
_pid_beamline_polarity,
_pid_beam_setting, global_event);
}
}
}
}
}
void MapCppGlobalPID::Through_DS_PID(std::vector*
GlobalTrackArray,
MAUS::GlobalEvent* global_event) const {
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") {
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-DS") {
MAUS::DataStructure::Global::Track* DS_track =
const_cast
(const_tracks.at(track_j));
perform_pid("MapCppGlobalPID-Final_Through_DS",
"Through_DS_PID", DS_track, _pid_vars,
_pid_beamline_polarity,
_pid_beam_setting, global_event);
}
}
}
}
}
void MapCppGlobalPID::Through_PID(std::vector*
GlobalTrackArray,
MAUS::GlobalEvent* global_event) const {
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") {
perform_pid("MapCppGlobalPID-Final_Through",
"Through_PID", track, _pid_vars,
_pid_beamline_polarity,
_pid_beam_setting, global_event);
}
}
}
void MapCppGlobalPID::perform_pid(std::string output_track_name,
std::string function_name,
MAUS::DataStructure::Global::Track* track,
std::vector _pid_vars,
std::string _pid_beamline_polarity,
std::string _pid_beam_setting,
MAUS::GlobalEvent* global_event) const {
MAUS::DataStructure::Global::Track* pidtrack = track->Clone();
global_event->add_track_recursive(pidtrack);
std::string candidate_name = "MapCppGlobalPID-Candidate_" + function_name;
pidtrack->set_mapper_name(candidate_name);
pidtrack->set_pid(MAUS::DataStructure::Global::kNoPID);
// doubles to hold cumulative log likelihoods for each hypothesis
double logL_mu = 0;
double logL_e = 0;
double logL_pi = 0;
/// Bools to check if log likelihood for PID variable returned an
/// "error value" logL, 1 for an unsuitable track, 2 for a track
/// that fails the PID due to a user defined cut on the variable
/// value. Only used for the purposes of efficiency calculations,
/// does not affect PID performance
bool logL_1 = false;
bool logL_2 = false;
// std::cerr << "_pid_vars.size() : " << _pid_vars.size() << std::endl;
for (size_t pid_var_count = 0; pid_var_count < _pid_vars.size();
++pid_var_count) {
// TODO Pidcott: Need to do this separately for +ve and -ve
// particles for now. Figure out a better way to organise this.
if (_pid_beamline_polarity == "positive") {
logL_1 = false;
logL_2 = false;
std::string hyp = _pid_vars[pid_var_count]->Get_hyp();
if (hyp == (_pid_beam_setting + "_mu_plus")) {
if (_pid_vars[pid_var_count]->logL(track) == 1) {
logL_1 = true;
continue;
} else if (_pid_vars[pid_var_count]->logL(track) == 2) {
logL_2 = true;
continue;
} else {
logL_mu += _pid_vars[pid_var_count]->logL(track);
logL_1 = false;
logL_2 = false;
}
} else if (hyp == (_pid_beam_setting + "_e_plus")) {
if (_pid_vars[pid_var_count]->logL(track) == 1) {
logL_1 = true;
continue;
} else if (_pid_vars[pid_var_count]->logL(track) == 2) {
logL_2 = true;
continue;
} else {
logL_e += _pid_vars[pid_var_count]->logL(track);
logL_1 = false;
logL_2 = false;
}
} else if (hyp == (_pid_beam_setting + "_pi_plus")) {
if (_pid_vars[pid_var_count]->logL(track) == 1) {
logL_1 = true;
continue;
} else if (_pid_vars[pid_var_count]->logL(track) == 2) {
logL_2 = true;
continue;
} else {
logL_pi += _pid_vars[pid_var_count]->logL(track);
logL_1 = false;
logL_2 = false;
}
} else {
Squeak::mout(Squeak::debug) << "Unrecognised particle hypothesis,"
<< " MapCppGlobalPID::perform_pid" << std::endl;
}
} else if (_pid_beamline_polarity == "negative") {
logL_1 = false;
logL_2 = false;
std::string hyp = _pid_vars[pid_var_count]->Get_hyp();
if (hyp == (_pid_beam_setting + "_mu_minus")) {
if (_pid_vars[pid_var_count]->logL(track) == 1) {
logL_1 = true;
continue;
} else if (_pid_vars[pid_var_count]->logL(track) == 2) {
logL_2 = true;
continue;
} else {
logL_mu += _pid_vars[pid_var_count]->logL(track);
logL_1 = false;
logL_2 = false;
}
} else if (hyp == (_pid_beam_setting + "_e_minus")) {
if (_pid_vars[pid_var_count]->logL(track) == 1) {
logL_1 = true;
continue;
} else if (_pid_vars[pid_var_count]->logL(track) == 2) {
logL_2 = true;
continue;
} else {
logL_e += _pid_vars[pid_var_count]->logL(track);
logL_1 = false;
logL_2 = false;
}
} else if (hyp == (_pid_beam_setting + "_pi_minus")) {
if (_pid_vars[pid_var_count]->logL(track) == 1) {
logL_1 = true;
continue;
} else if (_pid_vars[pid_var_count]->logL(track) == 2) {
logL_2 = true;
continue;
} else {
logL_pi += _pid_vars[pid_var_count]->logL(track);
logL_1 = false;
logL_2 = false;
}
} else {
Squeak::mout(Squeak::debug) << "Unrecognised particle hypothesis,"
<< " MapCppGlobalPID::perform_pid" << std::endl;
}
}
}
if ((logL_mu == 0 || logL_pi == 0 ||
logL_e == 0) && logL_1 == true) {
logL_mu = 1;
logL_pi = 1;
logL_e = 1;
} else if ((logL_mu == 0 || logL_pi == 0 ||
logL_e == 0) && logL_2 == true) {
logL_mu = 2;
logL_pi = 2;
logL_e = 2;
}
// make pid_ll_values
if (_pid_beamline_polarity == "positive") {
MAUS::DataStructure::Global::PIDLogLPair mu_LL(-13, logL_mu);
MAUS::DataStructure::Global::PIDLogLPair e_LL(-11, logL_e);
MAUS::DataStructure::Global::PIDLogLPair pi_LL(211, logL_pi);
pidtrack->AddPIDLogLValues(mu_LL);
pidtrack->AddPIDLogLValues(pi_LL);
pidtrack->AddPIDLogLValues(e_LL);
} else if (_pid_beamline_polarity == "negative") {
MAUS::DataStructure::Global::PIDLogLPair mu_LL(13, logL_mu);
MAUS::DataStructure::Global::PIDLogLPair e_LL(11, logL_e);
MAUS::DataStructure::Global::PIDLogLPair pi_LL(-211, logL_pi);
pidtrack->AddPIDLogLValues(mu_LL);
pidtrack->AddPIDLogLValues(pi_LL);
pidtrack->AddPIDLogLValues(e_LL);
}
/*std::vector
pid_ll_values = pidtrack->get_pid_logL_values();
for (std::vector::const_iterator
i = pid_ll_values.begin(); i != pid_ll_values.end(); ++i) {
std::cerr << i->first << "\t" << i->second << std::endl;
}*/
if (logL_mu < 0 || logL_pi < 0 || logL_e < 0) {
// calculate CLs
double sum_exp_LLs = exp(logL_mu) + exp(logL_e) + exp(logL_pi);
// std::cerr << "sum exp LLs: " << sum_exp_LLs << std::endl;
double CL_mu = ConfidenceLevel(logL_mu, sum_exp_LLs);
double CL_e = ConfidenceLevel(logL_e, sum_exp_LLs);
double CL_pi = ConfidenceLevel(logL_pi, sum_exp_LLs);
// std::cerr << "CL_mu " << CL_mu << std::endl;
// std::cerr << "CL_e " << CL_e << std::endl;
// std::cerr << "CL_pi " << CL_pi << std::endl;
// compare CLs and select winning hypothesis. set g.o.f. of track to CL
if (CL_mu - CL_e > _pid_confidence_level &&
CL_mu - CL_pi > _pid_confidence_level) {
if (_pid_beamline_polarity == "positive") {
pidtrack->set_pid(MAUS::DataStructure::Global::kMuPlus);
pidtrack->set_goodness_of_fit(CL_mu);
if (track->get_pid() == MAUS::DataStructure::Global::kMuPlus) {
MAUS::DataStructure::Global::Track* final_pidtrack = pidtrack->Clone();
global_event->add_track_recursive(final_pidtrack);
final_pidtrack->set_mapper_name(output_track_name);
}
} else if (_pid_beamline_polarity == "negative") {
pidtrack->set_pid(MAUS::DataStructure::Global::kMuMinus);
pidtrack->set_goodness_of_fit(CL_mu);
if (track->get_pid() == MAUS::DataStructure::Global::kMuMinus) {
MAUS::DataStructure::Global::Track* final_pidtrack = pidtrack->Clone();
global_event->add_track_recursive(final_pidtrack);
final_pidtrack->set_mapper_name(output_track_name);
}
}
} else if (CL_pi - CL_e > _pid_confidence_level &&
CL_pi - CL_mu > _pid_confidence_level) {
if (_pid_beamline_polarity == "positive") {
pidtrack->set_pid(MAUS::DataStructure::Global::kPiPlus);
pidtrack->set_goodness_of_fit(CL_pi);
if (track->get_pid() == MAUS::DataStructure::Global::kPiPlus) {
MAUS::DataStructure::Global::Track* final_pidtrack = pidtrack->Clone();
global_event->add_track_recursive(final_pidtrack);
final_pidtrack->set_mapper_name(output_track_name);
}
} else if (_pid_beamline_polarity == "negative") {
pidtrack->set_pid(MAUS::DataStructure::Global::kPiMinus);
pidtrack->set_goodness_of_fit(CL_pi);
if (track->get_pid() == MAUS::DataStructure::Global::kPiMinus) {
MAUS::DataStructure::Global::Track* final_pidtrack = pidtrack->Clone();
global_event->add_track_recursive(final_pidtrack);
final_pidtrack->set_mapper_name(output_track_name);
}
}
} else if (CL_e - CL_mu > _pid_confidence_level &&
CL_e - CL_pi > _pid_confidence_level) {
if (_pid_beamline_polarity == "positive") {
pidtrack->set_pid(MAUS::DataStructure::Global::kEPlus);
pidtrack->set_goodness_of_fit(CL_e);
if (track->get_pid() == MAUS::DataStructure::Global::kEPlus) {
MAUS::DataStructure::Global::Track* final_pidtrack = pidtrack->Clone();
global_event->add_track_recursive(final_pidtrack);
final_pidtrack->set_mapper_name(output_track_name);
}
} else if (_pid_beamline_polarity == "negative") {
pidtrack->set_pid(MAUS::DataStructure::Global::kEMinus);
pidtrack->set_goodness_of_fit(CL_e);
if (track->get_pid() == MAUS::DataStructure::Global::kEMinus) {
MAUS::DataStructure::Global::Track* final_pidtrack = pidtrack->Clone();
global_event->add_track_recursive(final_pidtrack);
final_pidtrack->set_mapper_name(output_track_name);
}
}
} else {
Squeak::mout(Squeak::debug) << "PID for track could not be" <<
" determined." << std::endl;
}
}
}
} // ~MAUS