← Back to branch summary

~mbogomilov/maus/devel3

« back to all changes in this revision

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

  • Committer: Durga Rajaram
  • Date: 2014-01-14 07:07:02 UTC
  • mfrom: (659.1.80 relcand)
  • Revision ID: durga@fnal.gov-20140114070702-2l1fuj1w6rraw7xe
Tags: MAUS-v0.7.6
MAUS-v0.7.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppGlobalTrackReconstructor/MapCppGlobalTrackReconstructor.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
 /* Author: Peter Lane
 
18
 */
 
19
 
 
20
#include "src/map/MapCppGlobalTrackReconstructor/MapCppGlobalTrackReconstructor.hh"
 
21
 
 
22
// C++
 
23
#include <algorithm>
 
24
#include <map>
 
25
#include <vector>
 
26
#include <string>
 
27
 
 
28
// External
 
29
#include "TLorentzVector.h"
 
30
#include "TMinuit.h"
 
31
#include "CLHEP/Units/PhysicalConstants.h"
 
32
#include "json/json.h"
 
33
 
 
34
// Legacy/G4MICE
 
35
#include "Utils/Exception.hh"
 
36
 
 
37
// MAUS
 
38
#include "Converter/DataConverters/JsonCppSpillConverter.hh"
 
39
#include "Converter/DataConverters/CppJsonSpillConverter.hh"
 
40
#include "DataStructure/GlobalEvent.hh"
 
41
#include "DataStructure/ReconEvent.hh"
 
42
#include "DataStructure/Spill.hh"
 
43
#include "DataStructure/Global/Track.hh"
 
44
#include "DataStructure/Global/TrackPoint.hh"
 
45
#include "DataStructure/ThreeVector.hh"
 
46
#include "src/common_cpp/Optics/CovarianceMatrix.hh"
 
47
#include "src/common_cpp/Optics/LinearApproximationOpticsModel.hh"
 
48
#include "src/common_cpp/Optics/PolynomialOpticsModel.hh"
 
49
#include "Recon/Global/DataStructureHelper.hh"
 
50
#include "Recon/Global/Detector.hh"
 
51
#include "Recon/Global/MinuitTrackFitter.hh"
 
52
#include "Recon/Global/Particle.hh"
 
53
#include "Recon/Global/TrackFitter.hh"
 
54
#include "Simulation/MAUSGeant4Manager.hh"
 
55
#include "Utils/JsonWrapper.hh"
 
56
#include "Utils/CppErrorHandler.hh"
 
57
 
 
58
namespace MAUS {
 
59
 
 
60
using MAUS::recon::global::DataStructureHelper;
 
61
using MAUS::recon::global::Detector;
 
62
using MAUS::recon::global::MinuitTrackFitter;
 
63
using MAUS::recon::global::Particle;
 
64
using MAUS::recon::global::TrackFitter;
 
65
namespace GlobalDS = MAUS::DataStructure::Global;
 
66
 
 
67
MapCppGlobalTrackReconstructor::MapCppGlobalTrackReconstructor()
 
68
    : optics_model_(NULL), track_fitter_(NULL) {
 
69
}
 
70
 
 
71
MapCppGlobalTrackReconstructor::~MapCppGlobalTrackReconstructor() {
 
72
  if (optics_model_ != NULL) {
 
73
    delete optics_model_;
 
74
  }
 
75
 
 
76
  if (track_fitter_ != NULL) {
 
77
    delete track_fitter_;
 
78
  }
 
79
}
 
80
 
 
81
bool MapCppGlobalTrackReconstructor::birth(std::string configuration_string) {
 
82
  // parse the JSON document.
 
83
  try {
 
84
    configuration_ = JsonWrapper::StringToJson(configuration_string);
 
85
    Json::Value physics_processes = configuration_["physics_processes"];
 
86
    if ((physics_processes != Json::Value("mean_energy_loss"))
 
87
        && (physics_processes != Json::Value("none"))) {
 
88
      throw(Exception(Exception::nonRecoverable,
 
89
                  "The \"physics_processes\" configuration variable should be "
 
90
                  "set to \"mean_energy_loss\" or \"none\" to avoid collisions "
 
91
                  "of the test particles with walls.",
 
92
                  "MAUS::MapCppGlobalTrackReconstructor::birth()"));
 
93
    }
 
94
    DataStructureHelper::GetInstance().GetDetectorAttributes(
 
95
        configuration_, detectors_);
 
96
    SetupOpticsModel();
 
97
    SetupTrackFitter();
 
98
  } catch(Exception& exc) {
 
99
    MAUS::CppErrorHandler::getInstance()->HandleExceptionNoJson(
 
100
      exc, MapCppGlobalTrackReconstructor::kClassname);
 
101
    return false;
 
102
  } catch(std::exception& exc) {
 
103
    MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(
 
104
      exc, MapCppGlobalTrackReconstructor::kClassname);
 
105
    return false;
 
106
  }
 
107
 
 
108
  return true;  // Sucessful parsing
 
109
}
 
110
 
 
111
std::string MapCppGlobalTrackReconstructor::process(
 
112
    std::string run_data_string) {
 
113
  // parse the JSON document.
 
114
  Json::Value run_data_json
 
115
    = Json::Value(JsonWrapper::StringToJson(run_data_string));
 
116
  if (run_data_json.isNull() || run_data_json.empty()) {
 
117
    return std::string("{\"errors\":{\"bad_json_document\":"
 
118
                        "\"Failed to parse input document\"}}");
 
119
  }
 
120
 
 
121
  JsonCppSpillConverter deserialize;
 
122
  MAUS::Data * run_data = deserialize(&run_data_json);
 
123
  if (!run_data) {
 
124
    return std::string("{\"errors\":{\"failed_json_cpp_conversion\":"
 
125
                        "\"Failed to convert Json to C++ Data object\"}}");
 
126
  }
 
127
 
 
128
  const MAUS::Spill * spill = run_data->GetSpill();
 
129
  MAUS::ReconEventPArray * recon_events = spill->GetReconEvents();
 
130
  if (!recon_events) {
 
131
    return run_data_string;
 
132
  }
 
133
 
 
134
  MAUS::ReconEventPArray::const_iterator recon_event;
 
135
  for (recon_event = recon_events->begin();
 
136
      recon_event != recon_events->end();
 
137
      ++recon_event) {
 
138
    MAUS::GlobalEvent * const global_event = (*recon_event)->GetGlobalEvent();
 
139
 
 
140
    GlobalDS::TrackPArray raw_tracks;
 
141
    LoadRawTracks(global_event, raw_tracks);
 
142
    std::cerr << "Loaded " << raw_tracks.size() << " raw tracks." << std::endl;
 
143
    // Fit each raw track and store best fit track in global_event
 
144
    GlobalDS::TrackPArray::const_iterator raw_track;
 
145
    for (raw_track = raw_tracks.begin();
 
146
         raw_track != raw_tracks.end();
 
147
         ++raw_track) {
 
148
      GlobalDS::Track * best_fit_track = new GlobalDS::Track();
 
149
      best_fit_track->set_mapper_name(kClassname);
 
150
      best_fit_track->AddTrack(*raw_track);
 
151
 
 
152
      std::cout << "DEBUG MapCppGlobalTrackReconstructor::process: "
 
153
                << "raw track PID: " << (*raw_track)->get_pid() << std::endl;
 
154
      track_fitter_->Fit(*raw_track, best_fit_track, kClassname);
 
155
 
 
156
      InsertIntermediateTrackPoints(best_fit_track);
 
157
 
 
158
      global_event->add_track_recursive(best_fit_track);
 
159
    }
 
160
 
 
161
    std::vector<MAUS::DataStructure::Global::Track*>* tracks
 
162
      = global_event->get_tracks();
 
163
    std::cout << "Added " << (tracks->size()-1)<< " reconstructed tracks:"
 
164
              << std::endl;
 
165
  }
 
166
 
 
167
  // Serialize the Spill for passing on to the next map in the workflow
 
168
  CppJsonSpillConverter serialize;
 
169
  Json::FastWriter writer;
 
170
  std::string output = writer.write(*serialize(run_data));
 
171
  std::cout << "DEBUG MapCppGlobalRawTracks::process(): "
 
172
            << "Output: " << std::endl
 
173
            << output << std::endl;
 
174
 
 
175
  delete run_data;
 
176
 
 
177
  return output;
 
178
}
 
179
 
 
180
bool MapCppGlobalTrackReconstructor::death() {
 
181
  return true;  // successful
 
182
}
 
183
 
 
184
void MapCppGlobalTrackReconstructor::SetupOpticsModel() {
 
185
fprintf(stdout, "CHECKPOINT: SetupOpticsModel() 0\n");
 
186
fflush(stdout);
 
187
  Json::Value optics_model_names = JsonWrapper::GetProperty(
 
188
      configuration_,
 
189
      "global_recon_optics_models",
 
190
      JsonWrapper::arrayValue);
 
191
  Json::Value optics_model_name = JsonWrapper::GetProperty(
 
192
      configuration_,
 
193
      "global_recon_optics_model",
 
194
      JsonWrapper::stringValue);
 
195
  size_t model;
 
196
  for (model = 0; model < optics_model_names.size(); ++model) {
 
197
    if (optics_model_name == optics_model_names[model]) {
 
198
      break;  // leave the current index into optics_model_names in model
 
199
    }
 
200
  }
 
201
 
 
202
fprintf(stdout, "CHECKPOINT: SetupOpticsModel() 1\n");
 
203
fflush(stdout);
 
204
  switch (model) {
 
205
    case 0: {
 
206
      // "Differentiating"
 
207
      /* TODO(plane1@hawk.iit.edu)
 
208
      optics_model_ = new DifferentiatingOpticsModel();
 
209
      */
 
210
      throw(Exception(Exception::nonRecoverable,
 
211
                   "DifferentiatingOpticsModel is not yet implemented.",
 
212
                   "MapCppGlobalTrackReconstructor::SetupOpticsModel()"));
 
213
      break;
 
214
    }
 
215
    case 1: {
 
216
      // "Integrating"
 
217
      /* TODO(plane1@hawk.iit.edu)
 
218
      optics_model_ = new IntegratingOpticsModel();
 
219
      */
 
220
      throw(Exception(Exception::nonRecoverable,
 
221
                   "IntegratingOpticsModel is not yet implemented.",
 
222
                   "MapCppGlobalTrackReconstructor::SetupOpticsModel()()"));
 
223
      break;
 
224
    }
 
225
    case 2: {
 
226
fprintf(stdout, "CHECKPOINT SetupOpticsModel() 1.1c\n");
 
227
fflush(stdout);
 
228
     optics_model_ = new PolynomialOpticsModel(&configuration_);
 
229
      break;
 
230
    }
 
231
    case 3: {
 
232
      // "Runge Kutta"
 
233
      /* TODO(plane1@hawk.iit.edu)
 
234
      optics_model_ = new RungeKuttaOpticsModel();
 
235
      */
 
236
      throw(Exception(Exception::nonRecoverable,
 
237
                   "RungeKuttaOpticsModel is not yet implemented.",
 
238
                   "MapCppGlobalTrackReconstructor::SetupOpticsModel()()"));
 
239
      break;
 
240
    }
 
241
    case 4: {
 
242
      // "Linear Approximation"
 
243
fprintf(stdout, "CHECKPOINT SetupOpticsModel() 1.1e\n");
 
244
fflush(stdout);
 
245
     optics_model_ = new LinearApproximationOpticsModel(&configuration_);
 
246
      break;
 
247
    }
 
248
    default: {
 
249
      std::string message("Unsupported optics model \"");
 
250
      message += optics_model_name.asString();
 
251
      message += std::string("\" in recon configuration.");
 
252
      throw(Exception(Exception::nonRecoverable,
 
253
                   message.c_str(),
 
254
                   "MapCppGlobalTrackReconstructor::SetupOpticsModel()()"));
 
255
    }
 
256
  }
 
257
fprintf(stdout, "CHECKPOINT: SetupOpticsModel() 2\n");
 
258
fflush(stdout);
 
259
  optics_model_->Build();
 
260
 
 
261
fprintf(stdout, "CHECKPOINT: SetupOpticsModel() 3\n");
 
262
fflush(stdout);
 
263
}
 
264
 
 
265
void MapCppGlobalTrackReconstructor::SetupTrackFitter() {
 
266
  Json::Value fitter_names = JsonWrapper::GetProperty(
 
267
      configuration_,
 
268
      "global_recon_track_fitters",
 
269
      JsonWrapper::arrayValue);
 
270
  Json::Value fitter_name = JsonWrapper::GetProperty(
 
271
      configuration_,
 
272
      "global_recon_track_fitter",
 
273
      JsonWrapper::stringValue);
 
274
  size_t fitter;
 
275
  for (fitter = 0; fitter < fitter_names.size(); ++fitter) {
 
276
    if (fitter_name == fitter_names[fitter]) {
 
277
      break;  // leave the current index into fitter_names in fitter
 
278
    }
 
279
  }
 
280
 
 
281
  switch (fitter) {
 
282
    case 0: {
 
283
      // Minuit
 
284
      double start_plane = optics_model_->primary_plane();
 
285
      track_fitter_ = new MinuitTrackFitter(optics_model_, start_plane);
 
286
      break;
 
287
    }
 
288
    case 1: {
 
289
      // Kalman Filter
 
290
      /* TODO(plane1@hawk.iit.edu)
 
291
      track_fitter_ = new KalmanFilterTrackFitter(optics_model_);
 
292
      */
 
293
      throw(Exception(Exception::nonRecoverable,
 
294
                   "KalmanFilterTrackFitter is not yet implemented.",
 
295
                   "MapCppGlobalTrackReconstructor::SetupTrackFitter()"));
 
296
      break;
 
297
    }
 
298
    case 2: {
 
299
      // TOF
 
300
      /* TODO(plane1@hawk.iit.edu)
 
301
      track_fitter_ = new TOFTrackFitter(optics_model_);
 
302
      */
 
303
      throw(Exception(Exception::nonRecoverable,
 
304
                   "TOFTrackFitter is not yet implemented.",
 
305
                   "MapCppGlobalTrackReconstructor::SetupTrackFitter()"));
 
306
      break;
 
307
    }
 
308
    default: {
 
309
      std::string message("Unsupported track fitter \"");
 
310
      message += fitter_name.asString();
 
311
      message += std::string("\" in recon configuration.");
 
312
      throw(Exception(Exception::nonRecoverable,
 
313
                   message.c_str(),
 
314
                   "MapCppGlobalTrackReconstructor::SetupTrackFitter()"));
 
315
    }
 
316
  }
 
317
}
 
318
 
 
319
void MapCppGlobalTrackReconstructor::LoadRawTracks(
 
320
    GlobalEvent const * const global_event,
 
321
    GlobalDS::TrackPArray & tracks) const {
 
322
  GlobalDS::TrackPArray * global_tracks = global_event->get_tracks();
 
323
  GlobalDS::TrackPArray::iterator global_track;
 
324
  const std::string recon_mapper_name("MapCppGlobalRawTracks");
 
325
  for (global_track = global_tracks->begin();
 
326
       global_track != global_tracks->end();
 
327
       ++global_track) {
 
328
    if ((*global_track)->get_mapper_name() == recon_mapper_name) {
 
329
      tracks.push_back(*global_track);
 
330
    }
 
331
  }
 
332
}
 
333
 
 
334
void MapCppGlobalTrackReconstructor::InsertIntermediateTrackPoints(
 
335
    GlobalDS::Track * track) const {
 
336
  PolynomialOpticsModel * const optics_model
 
337
    = static_cast<MAUS::PolynomialOpticsModel*>(optics_model_);
 
338
  if (optics_model == NULL) {
 
339
    throw(Exception(Exception::nonRecoverable,
 
340
                "Could not reconstruct intermediate track points: the optics "
 
341
                "model being used is not yet fully supported.",
 
342
                "MAUS::MapCppGlobalTrackReconstructor::"
 
343
                "InsertIntermediateTrackPoints()"));
 
344
  }
 
345
 
 
346
  // Get the fit track points
 
347
  GlobalDS::TrackPointCPArray fit_points = track->GetTrackPoints();
 
348
 
 
349
  // Convert the first fit track point into a PhaseSpaceVector for transporting
 
350
  GlobalDS::TrackPoint const * const fit_primary_track_point = fit_points[0];
 
351
  size_t particle_event = fit_primary_track_point->get_particle_event();
 
352
  DataStructureHelper helper = DataStructureHelper::GetInstance();
 
353
  PhaseSpaceVector fit_primary
 
354
    = helper.TrackPoint2PhaseSpaceVector(*fit_primary_track_point);
 
355
 
 
356
  // Construct a list of detector z-keys (z-position rounded to nearest integer)
 
357
  std::vector<int64_t> z_keys;
 
358
  GlobalDS::TrackPointCPArray::const_iterator fit_point
 
359
    = fit_points.begin();
 
360
  for (fit_point = fit_points.begin();
 
361
       fit_point != fit_points.end();
 
362
       ++fit_point) {
 
363
    // calculate the next guess
 
364
    const double z = (*fit_point)->get_position().Z();
 
365
    int64_t z_key = (z >= 0?static_cast<int64_t>(z+.5):static_cast<int64_t>(z-.5));
 
366
    z_keys.push_back(z_key);
 
367
  }
 
368
 
 
369
  // Reconstruct the intermediate track points by transporting the fit primary
 
370
  // to all desired intermediate z-positions
 
371
  const GlobalDS::PID particle_id = track->get_pid();
 
372
  std::vector<int64_t>::const_iterator z_key = z_keys.begin();
 
373
  const std::vector<int64_t> map_positions
 
374
    = optics_model->GetAvailableMapPositions();
 
375
  std::vector<int64_t>::const_iterator map_z;
 
376
  for (map_z = map_positions.begin(); map_z != map_positions.end(); ++map_z) {
 
377
    // locate the next detector z-position
 
378
    while (z_key != z_keys.end() && (*z_key) < (*map_z)) {
 
379
      ++z_key;
 
380
    }
 
381
 
 
382
    // transport the fit primary to the desired z-position
 
383
    const PhaseSpaceVector point = optics_model_->Transport(fit_primary,
 
384
                                                            *map_z);
 
385
    std::cout << "DEBUG MapCppGlobalTrackReconstructor"
 
386
              << "::InsertIntermediateTrackPoints: track point: "
 
387
              << point << std::endl;
 
388
 
 
389
    // skip if we've already added this point to the track during fitting
 
390
    if ((*z_key) == (*map_z)) {
 
391
      continue;
 
392
    }
 
393
 
 
394
    // add the reconstructed intermediate to the track
 
395
    GlobalDS::TrackPoint track_point;
 
396
    try {
 
397
      track_point
 
398
        = helper.PhaseSpaceVector2TrackPoint(point, *map_z, particle_id);
 
399
    } catch (Exception exc) {
 
400
        std::cerr << "DEBUG MapCppGlobalTrackReconstructor"
 
401
                  << "::InsertIntermediateTrackPoints: "
 
402
                  << "something bad happened during track fitting: "
 
403
                  << exc.what() << std::endl;
 
404
    }
 
405
 
 
406
    track_point.set_mapper_name(kClassname);
 
407
    track_point.set_particle_event(particle_event);
 
408
    track->AddTrackPoint(new GlobalDS::TrackPoint(track_point));
 
409
  }
 
410
  track->SortTrackPointsByZ();
 
411
}
 
412
 
 
413
const std::string MapCppGlobalTrackReconstructor::kClassname
 
414
  = "MapCppGlobalTrackReconstructor";
 
415
 
 
416
}  // namespace MAUS
 
417