← Back to branch summary

~chris-rogers/maus/emr_mc_digitization

« back to all changes in this revision

Viewing changes to src/map/MapCppGlobalPID/MapCppGlobalPID.cc

  • Committer: Chris Rogers
  • Date: 2014-04-16 11:48:45 UTC
  • mfrom: (707 merge)
  • mto: This revision was merged to the branch mainline in revision 711.
  • Revision ID: chris.rogers@stfc.ac.uk-20140416114845-h3u3q7pdcxkxvovs
Update to trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppGlobalPID/MapCppGlobalPID.cc
 
1
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
 
2
 *
 
3
 * MAUS is free software: you can redistribute it and/or modify
 
4
 * it under the terms of the GNU General Public License as published by
 
5
 * the Free Software Foundation, either version 3 of the License, or
 
6
 * (at your option) any later version.
 
7
 *
 
8
 * MAUS is distributed in the hope that it will be useful,
 
9
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
 * GNU General Public License for more details.
 
12
 *
 
13
 * You should have received a copy of the GNU General Public License
 
14
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
 
15
 *
 
16
 */
 
17
 
 
18
#include "src/map/MapCppGlobalPID/MapCppGlobalPID.hh"
 
19
 
 
20
#include "Interface/Squeak.hh"
 
21
#include "src/common_cpp/DataStructure/Data.hh"
 
22
#include "src/common_cpp/Converter/DataConverters/JsonCppSpillConverter.hh"
 
23
#include "src/common_cpp/Converter/DataConverters/CppJsonSpillConverter.hh"
 
24
 
 
25
 
 
26
namespace MAUS {
 
27
  MapCppGlobalPID::MapCppGlobalPID() {
 
28
    _classname = "MapCppGlobalPID";
 
29
  }
 
30
 
 
31
  bool MapCppGlobalPID::birth(std::string argJsonConfigDocument) {
 
32
    // Check if the JSON document can be parsed, else return error only.
 
33
    bool parsingSuccessful = _reader.parse(argJsonConfigDocument, _configJSON);
 
34
    if (!parsingSuccessful) {
 
35
      _configCheck = false;
 
36
      return false;
 
37
    }
 
38
 
 
39
    char* pMAUS_ROOT_DIR = getenv("MAUS_ROOT_DIR");
 
40
    if (!pMAUS_ROOT_DIR) {
 
41
      Squeak::mout(Squeak::error)
 
42
        << "Could not find the $MAUS_ROOT_DIR env variable." << std::endl;
 
43
      Squeak::mout(Squeak::error)
 
44
        << "Did you try running: source env.sh ?" << std::endl;
 
45
      return false;
 
46
    }
 
47
 
 
48
    _configCheck = true;
 
49
 
 
50
    _hypotheses.clear();
 
51
    _pid_vars.clear();
 
52
 
 
53
    PDF_file = _configJSON["PID_PDFs_file"].asString();
 
54
 
 
55
    _histFile = new TFile(PDF_file.c_str(), "READ");
 
56
 
 
57
    // vector of hypotheses
 
58
    // TODO(Pidcott) find a more elegant way of accessing hypotheses
 
59
    _hypotheses.push_back("200MeV_mu_plus");
 
60
    _hypotheses.push_back("200MeV_e_plus");
 
61
    _hypotheses.push_back("200MeV_pi_plus");
 
62
 
 
63
    for (unsigned int i =0; i < _hypotheses.size(); ++i) {
 
64
      // vector of pid vars
 
65
      _pid_vars.push_back(new MAUS::recon::global::PIDVarA(_histFile,
 
66
                                                           _hypotheses[i]));
 
67
      // _pid_vars.push_back(new MAUS::recon::global::PIDVarB(histFile,
 
68
      //                                                      _hypotheses[i]));
 
69
      // etc.
 
70
      }
 
71
 
 
72
    return true;
 
73
  }
 
74
 
 
75
  bool MapCppGlobalPID::death() {
 
76
    return true;
 
77
  }
 
78
 
 
79
  std::string MapCppGlobalPID::process(std::string document) const {
 
80
    Json::FastWriter writer;
 
81
    Json::Value root;
 
82
 
 
83
    if (document.empty()) {
 
84
      Json::Value errors;
 
85
      std::stringstream ss;
 
86
      ss << _classname << " says: Empty document passed to process";
 
87
      errors["bad_json_document"] = ss.str();
 
88
      root["errors"] = errors;
 
89
      return writer.write(root);
 
90
    }
 
91
 
 
92
    if (!_configCheck) {
 
93
      Json::Value errors;
 
94
      std::stringstream ss;
 
95
      ss << _classname << " says: process was not passed a valid configuration";
 
96
      errors["bad_json_document"] = ss.str();
 
97
      root["errors"] = errors;
 
98
      return writer.write(root);
 
99
    }
 
100
 
 
101
    // Prepare converters, spill and json objects
 
102
    JsonCppSpillConverter json2cppconverter;
 
103
    CppJsonSpillConverter cpp2jsonconverter;
 
104
    Json::Value *data_json = NULL;
 
105
    MAUS::Data *data_cpp = NULL;
 
106
 
 
107
    // Read string and convert to a Json object
 
108
    try {
 
109
      Json::Value imported_json = JsonWrapper::StringToJson(document);
 
110
      data_json = new Json::Value(imported_json);
 
111
    } catch (...) {
 
112
      Json::Value errors;
 
113
      std::stringstream ss;
 
114
      ss << _classname << " says: Bad json document";
 
115
      errors["bad_json_document"] = ss.str();
 
116
      root["errors"] = errors;
 
117
      return writer.write(root);
 
118
    }
 
119
 
 
120
    if (!data_json || data_json->isNull()) {
 
121
      if (data_json) delete data_json;
 
122
      return std::string("{\"errors\":{\"bad_json_document\":")+
 
123
        std::string("\"Failed to parse input document\"}}");
 
124
    }
 
125
 
 
126
    if (data_json->empty()) {
 
127
      delete data_json;
 
128
      return std::string("{\"errors\":{\"bad_json_document\":")+
 
129
        std::string("\"Failed to parse input document\"}}");
 
130
    }
 
131
 
 
132
    std::string maus_event = JsonWrapper::GetProperty(
 
133
        *data_json, "maus_event_type",
 
134
        JsonWrapper::stringValue).asString();
 
135
 
 
136
    if ( maus_event.compare("Spill") != 0 ) {
 
137
      Squeak::mout(Squeak::error) << "Line of json document did not contain "
 
138
          << "a Spill" << std::endl;
 
139
      delete data_json;
 
140
      return document;
 
141
    }
 
142
 
 
143
    std::string daq_event = JsonWrapper::GetProperty(
 
144
        *data_json, "daq_event_type",
 
145
        JsonWrapper::stringValue).asString();
 
146
 
 
147
    if ( daq_event.compare("physics_event") != 0 ) {
 
148
      Squeak::mout(Squeak::error) << "daq_event_type did not return a "
 
149
                                  << "physics event" << std::endl;
 
150
      delete data_json;
 
151
      return document;
 
152
    }
 
153
 
 
154
    // Convert Json into MAUS::Spill object.  In future, this will all
 
155
    // be done for me, and process will take/return whichever object we
 
156
    // prefer.
 
157
    try {
 
158
      data_cpp = json2cppconverter(data_json);
 
159
      delete data_json;
 
160
    } catch (...) {
 
161
      Squeak::mout(Squeak::error) << "Missing required branch daq_event_type"
 
162
          << " converting json->cpp, MapCppGlobalPID" << std::endl;
 
163
    }
 
164
 
 
165
    if (!data_cpp) {
 
166
      return std::string("{\"errors\":{\"failed_json_cpp_conversion\":")+
 
167
        std::string("\"Failed to convert Json to Cpp Spill object\"}}");
 
168
    }
 
169
 
 
170
    const MAUS::Spill* _spill = data_cpp->GetSpill();
 
171
 
 
172
    if ( _spill->GetReconEvents() ) {
 
173
      for ( unsigned int event_i = 0;
 
174
            event_i < _spill->GetReconEvents()->size(); ++event_i ) {
 
175
        MAUS::GlobalEvent* global_event =
 
176
          _spill->GetReconEvents()->at(event_i)->GetGlobalEvent();
 
177
        std::vector<MAUS::DataStructure::Global::Track*> *GlobalTrackArray =
 
178
          global_event->get_tracks();
 
179
        for (unsigned int track_i = 0; track_i < GlobalTrackArray->size();
 
180
             ++track_i) {
 
181
          MAUS::DataStructure::Global::Track* track =
 
182
            GlobalTrackArray->at(track_i);
 
183
          // doubles to hold cumulative log likelihoods for each hypothesis
 
184
          double logL_200MeV_mu_plus = 0;
 
185
          double logL_200MeV_e_plus = 0;
 
186
          double logL_200MeV_pi_plus = 0;
 
187
          for (size_t pid_var_count = 0; pid_var_count < _pid_vars.size();
 
188
               ++pid_var_count) {
 
189
            std::string hyp = _pid_vars[pid_var_count]->Get_hyp();
 
190
            if (hyp == "200MeV_mu_plus") {
 
191
              if (_pid_vars[pid_var_count]->logL(track) == 1) {
 
192
                continue;
 
193
              } else {
 
194
                logL_200MeV_mu_plus += _pid_vars[pid_var_count]->logL(track);
 
195
              }
 
196
            } else if (hyp == "200MeV_e_plus") {
 
197
              if (_pid_vars[pid_var_count]->logL(track) == 1) {
 
198
                continue;
 
199
              } else {
 
200
                logL_200MeV_e_plus += _pid_vars[pid_var_count]->logL(track);
 
201
              }
 
202
            } else if (hyp == "200MeV_pi_plus") {
 
203
              if (_pid_vars[pid_var_count]->logL(track) == 1) {
 
204
                continue;
 
205
              } else {
 
206
                logL_200MeV_pi_plus += _pid_vars[pid_var_count]->logL(track);
 
207
              }
 
208
            } else {
 
209
              Squeak::mout(Squeak::error) << "Unrecognised particle hypothesis,"
 
210
                  << " MapCppGlobalPID::process" << std::endl;
 
211
            }
 
212
          }
 
213
          if ((logL_200MeV_mu_plus - logL_200MeV_e_plus > 0.5) &&
 
214
              (logL_200MeV_mu_plus - logL_200MeV_pi_plus > 0.5)) {
 
215
            track->set_pid(MAUS::DataStructure::Global::kMuPlus);
 
216
          } else if ((logL_200MeV_e_plus - logL_200MeV_mu_plus > 0.5) &&
 
217
                     (logL_200MeV_e_plus - logL_200MeV_pi_plus > 0.5)) {
 
218
            track->set_pid(MAUS::DataStructure::Global::kEPlus);
 
219
          } else if ((logL_200MeV_pi_plus - logL_200MeV_mu_plus > 0.5) &&
 
220
                     (logL_200MeV_pi_plus - logL_200MeV_e_plus > 0.5)) {
 
221
            track->set_pid(MAUS::DataStructure::Global::kPiPlus);
 
222
          } else {
 
223
            Squeak::mout(Squeak::error) << "PID for track could not be" <<
 
224
              " determined." << std::endl;
 
225
            continue;
 
226
          }
 
227
          Squeak::mout(Squeak::error) << "PID of track : " <<
 
228
            track->get_pid() << std::endl;
 
229
        }
 
230
      }
 
231
      }
 
232
 
 
233
    data_json = cpp2jsonconverter(data_cpp);
 
234
 
 
235
    if (!data_json) {
 
236
      delete data_cpp;
 
237
      return std::string("{\"errors\":{\"failed_cpp_json_conversion\":")+
 
238
        std::string("\"Failed to convert Cpp to Json Spill object\"}}");
 
239
    }
 
240
 
 
241
    std::string output_document = JsonWrapper::JsonToString(*data_json);
 
242
    delete data_json;
 
243
    delete data_cpp;
 
244
    return output_document;
 
245
  }
 
246
} // ~MAUS