← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

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

merging in changes in merge branch

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerRecon/MapCppTrackerRecon.cc
24
24
 
25
25
namespace MAUS {
26
26
 
27
 
// void print_state(Kalman::State& state);
28
 
 
29
 
void print_helical(SciFiHelicalPRTrack* helical);
30
 
void print_straight(SciFiStraightPRTrack* straight);
31
 
 
32
 
 
33
 
SciFiClusterPArray find_clusters(SciFiClusterPArray cluster_array, int tracker);
34
 
SciFiSpacePointPArray find_spacepoints(SciFiSpacePointPArray spacepoint_array, int tracker);
35
 
 
36
 
 
37
27
PyMODINIT_FUNC init_MapCppTrackerRecon(void) {
38
28
  PyWrapMapBase<MapCppTrackerRecon>::PyWrapMapBaseModInit
39
29
                                        ("MapCppTrackerRecon", "", "", "", "");
40
30
}
41
31
 
42
 
MapCppTrackerRecon::MapCppTrackerRecon() : _up_straight_pr_on(true),
 
32
 
 
33
MapCppTrackerRecon::MapCppTrackerRecon() : MapBase<Data>("MapCppTrackerRecon"),
 
34
                                           _clusters_on(true),
 
35
                                           _spacepoints_on(true),
 
36
                                           _up_straight_pr_on(true),
43
37
                                           _down_straight_pr_on(true),
44
38
                                           _up_helical_pr_on(true),
45
39
                                           _down_helical_pr_on(true),
46
 
                                           MapBase<Data>("MapCppTrackerRecon"),
47
 
#ifdef KALMAN_TEST
48
 
  _spacepoint_helical_track_fitter(NULL),
49
 
  _spacepoint_straight_track_fitter(NULL) {
50
 
#else
51
 
  _helical_track_fitter(NULL),
52
 
  _straight_track_fitter(NULL) {
53
 
#endif
 
40
                                           _kalman_on(true),
 
41
                                           _patrec_on(true),
 
42
                                           _patrec_debug_on(false),
 
43
                                           _helical_track_fitter(NULL),
 
44
                                           _straight_track_fitter(NULL) {
54
45
}
55
46
 
 
47
 
56
48
MapCppTrackerRecon::~MapCppTrackerRecon() {
57
49
}
58
50
 
 
51
 
59
52
void MapCppTrackerRecon::_birth(const std::string& argJsonConfigDocument) {
 
53
  // Pull out the global settings
60
54
  if (!Globals::HasInstance()) {
61
55
    GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
62
56
  }
66
60
  int user_up_straight_pr_on   = (*json)["SciFiPRStraightTkUSOn"].asInt();
67
61
  int user_down_straight_pr_on = (*json)["SciFiPRStraightTkDSOn"].asInt();
68
62
  _kalman_on          = (*json)["SciFiKalmanOn"].asBool();
 
63
  _clusters_on        = (*json)["SciFiClusterReconOn"].asBool();
 
64
  _spacepoints_on     = (*json)["SciFiSpacepointReconOn"].asBool();
69
65
  _patrec_on          = (*json)["SciFiPatRecOn"].asBool();
 
66
  _patrec_debug_on    = (*json)["SciFiPatRecDebugOn"].asBool();
70
67
  _size_exception     = (*json)["SciFiClustExcept"].asInt();
71
68
  _min_npe            = (*json)["SciFiNPECut"].asDouble();
72
69
  _use_mcs            = (*json)["SciFiKalman_use_MCS"].asBool();
73
70
  _use_eloss          = (*json)["SciFiKalman_use_Eloss"].asBool();
74
71
  _use_patrec_seed    = (*json)["SciFiSeedPatRec"].asBool();
75
72
  _correct_pz         = (*json)["SciFiKalmanCorrectPz"].asBool();
76
 
// Values used to set the track rating:
 
73
  // Values used to set the track rating:
77
74
  _excellent_num_trackpoints     = (*json)["SciFiExcellentNumTrackpoints"].asInt();
78
75
  _good_num_trackpoints          = (*json)["SciFiGoodNumTrackpoints"].asInt();
79
76
  _poor_num_trackpoints          = (*json)["SciFiPoorNumTrackpoints"].asInt();
90
87
  logfile.open(path);
91
88
  logfile.close();
92
89
 
 
90
  // Build the geometery helper instance
93
91
  MiceModule* module = Globals::GetReconstructionMiceModules();
94
92
  std::vector<const MiceModule*> modules =
95
93
    module->findModulesByPropertyString("SensitiveDetector", "SciFi");
96
94
  _geometry_helper = SciFiGeometryHelper(modules);
97
95
  _geometry_helper.Build();
98
96
 
 
97
  // Set up cluster reconstruction
99
98
  _cluster_recon = SciFiClusterRec(_size_exception, _min_npe, _geometry_helper.GeometryMap());
100
99
 
 
100
  // Set up spacepoint reconstruction
101
101
  _spacepoint_recon = SciFiSpacePointRec();
102
102
 
 
103
  // Setup Pattern Recognition
103
104
  double up_field = _geometry_helper.GetFieldValue(0);
104
105
  double down_field = _geometry_helper.GetFieldValue(1);
105
 
 
106
106
  if (user_up_helical_pr_on == 2)
107
107
    _up_helical_pr_on = true;
108
108
  else if (user_up_helical_pr_on == 1)
147
147
  _pattern_recognition.set_down_straight_pr_on(_down_straight_pr_on);
148
148
  _pattern_recognition.set_bz_t1(up_field);
149
149
  _pattern_recognition.set_bz_t2(down_field);
 
150
  if (_patrec_debug_on) _pattern_recognition.setup_debug();
150
151
 
151
152
  if (_use_patrec_seed) {
152
153
    _seed_value = -1.0;
154
155
    _seed_value = (*json)["SciFiSeedCovariance"].asDouble();
155
156
  }
156
157
 
157
 
#ifdef KALMAN_TEST
158
 
  HelicalPropagator* spacepoint_helical_prop = new HelicalPropagator(&_geometry_helper);
159
 
  spacepoint_helical_prop->SetCorrectPz(_correct_pz);
160
 
  spacepoint_helical_prop->SetIncludeMCS(_use_mcs);
161
 
  spacepoint_helical_prop->SetSubtractELoss(_use_eloss);
162
 
  _spacepoint_helical_track_fitter = new Kalman::TrackFit(spacepoint_helical_prop,
163
 
      new SciFiSpacepointMeasurements<5>(0.1));
164
 
 
165
 
  StraightPropagator* spacepoint_straight_prop = new StraightPropagator(&_geometry_helper);
166
 
  spacepoint_straight_prop->SetIncludeMCS(_use_mcs);
167
 
  _spacepoint_straight_track_fitter = new Kalman::TrackFit(spacepoint_straight_prop,
168
 
      new SciFiSpacepointMeasurements<4>(0.1));
169
 
 
170
 
  _spacepoint_recon_plane = 0;
171
 
#else
 
158
  // Set up final track fit (Kalman filter)
172
159
  HelicalPropagator* helical_prop = new HelicalPropagator(&_geometry_helper);
173
160
  helical_prop->SetCorrectPz(_correct_pz);
174
161
  helical_prop->SetIncludeMCS(_use_mcs);
180
167
  straight_prop->SetIncludeMCS(_use_mcs);
181
168
  _straight_measurement = new SciFiStraightMeasurements(&_geometry_helper);
182
169
  _straight_track_fitter = new Kalman::TrackFit(straight_prop, _straight_measurement);
183
 
#endif
184
170
}
185
171
 
 
172
 
186
173
void MapCppTrackerRecon::_death() {
187
 
#ifdef KALMAN_TEST
188
 
  if (_spacepoint_helical_track_fitter) {
189
 
    delete _spacepoint_helical_track_fitter;
190
 
    _spacepoint_helical_track_fitter = NULL;
191
 
  }
192
 
  if (_spacepoint_straight_track_fitter) {
193
 
    delete _spacepoint_straight_track_fitter;
194
 
    _spacepoint_straight_track_fitter = NULL;
195
 
  }
196
 
#else
197
174
  if (_helical_track_fitter) {
198
175
    delete _helical_track_fitter;
199
176
    _helical_track_fitter = NULL;
202
179
    delete _straight_track_fitter;
203
180
    _straight_track_fitter = NULL;
204
181
  }
205
 
#endif
206
182
}
207
183
 
 
184
 
208
185
void MapCppTrackerRecon::_process(Data* data) const {
209
186
  Spill& spill = *(data->GetSpill());
210
187
 
219
196
        continue;
220
197
      }
221
198
 
222
 
      // Build Clusters.
223
 
      if (event->digits().size()) {
224
 
        _cluster_recon.process(*event);
 
199
      if ( _clusters_on ) {
 
200
      // Clear any exising higher level data
 
201
        event->clear_clusters();
 
202
        // Build Clusters.
 
203
        if (event->digits().size()) {
 
204
          _cluster_recon.process(*event);
 
205
        }
225
206
      }
226
 
      // Build SpacePoints.
227
 
      if (event->clusters().size()) {
228
 
        _spacepoint_recon.process(*event);
229
 
        set_spacepoint_global_output(event->spacepoints());
 
207
      if ( _spacepoints_on ) {
 
208
        // Build SpacePoints.
 
209
        event->clear_spacepoints();
 
210
        if (event->clusters().size()) {
 
211
          _spacepoint_recon.process(*event);
 
212
          set_spacepoint_global_output(event->spacepoints());
 
213
        }
230
214
      }
231
215
      // Pattern Recognition.
232
216
      if (_patrec_on && event->spacepoints().size()) {
 
217
        event->clear_seeds();
 
218
        event->clear_stracks();
 
219
        event->clear_htracks();
233
220
        _pattern_recognition.process(*event);
 
221
        set_straight_prtrack_global_output(event->straightprtracks());
234
222
        extrapolate_helical_reference(*event);
235
223
        extrapolate_straight_reference(*event);
236
224
      }
237
225
      // Kalman Track Fit.
238
226
      if (_kalman_on) {
 
227
        event->clear_scifitracks();
239
228
        track_fit(*event);
240
229
      }
241
230
//      print_event_info(*event);
245
234
  }
246
235
}
247
236
 
 
237
 
248
238
void MapCppTrackerRecon::track_fit(SciFiEvent &evt) const {
249
 
#ifdef KALMAN_TEST
250
 
  std::cerr << "KALMAN_TEST ACTIVE\n";
251
 
  if (_patrec_on) {
252
 
    std::cerr << "PATREC ACTIVE\n";
253
 
    size_t number_helical_tracks = evt.helicalprtracks().size();
254
 
    size_t number_straight_tracks = evt.straightprtracks().size();
255
 
 
256
 
    for (size_t track_i = 0; track_i < number_helical_tracks; track_i++) {
257
 
      SciFiHelicalPRTrack* helical = evt.helicalprtracks().at(track_i);
258
 
      SciFiSpacePointPArray spacepoints = helical->get_spacepoints_pointers();
259
 
 
260
 
      Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
261
 
                                                                          _spacepoint_recon_plane);
262
 
 
263
 
      Kalman::State seed = ComputeSeed(helical, &_geometry_helper, false, _seed_value);
264
 
 
265
 
      std::cerr << Kalman::print_state(seed, "Helical Seed");
266
 
      seed.GetCovariance().Print();
267
 
 
268
 
      std::cerr << "Helical Track\n";
269
 
      print_helical(helical);
270
 
 
271
 
      _spacepoint_helical_track_fitter->SetSeed(seed);
272
 
      _spacepoint_helical_track_fitter->SetData(data_track);
273
 
 
274
 
      _spacepoint_helical_track_fitter->Filter(false);
275
 
      _spacepoint_helical_track_fitter->Smooth(false);
276
 
 
277
 
      Kalman::Track predicted = _spacepoint_helical_track_fitter->GetPredicted();
278
 
      Kalman::Track filtered = _spacepoint_helical_track_fitter->GetFiltered();
279
 
      Kalman::Track smoothed = _spacepoint_helical_track_fitter->GetSmoothed();
280
 
      std::cerr << "Data Track\n";
281
 
      std::cerr << Kalman::print_track(data_track);
282
 
      std::cerr << "Measured Predicted\n";
283
 
      SciFiSpacepointMeasurements<5>* temp_measurement = new SciFiSpacepointMeasurements<5>(1.0);
284
 
      for (unsigned int k = 0; k < predicted.GetLength(); ++k) {
285
 
        Kalman::State temp_state = temp_measurement->Measure(predicted[k]);
286
 
        std::cerr << Kalman::print_state(temp_state);
287
 
      }
288
 
      std::cerr << "Predicted Track\n";
289
 
      std::cerr << Kalman::print_track(predicted);
290
 
      std::cerr << "Filtered Track\n";
291
 
      std::cerr << Kalman::print_track(filtered);
292
 
      std::cerr << "Smoothed Track\n";
293
 
      std::cerr << Kalman::print_track(smoothed);
294
 
      delete temp_measurement;
295
 
 
296
 
      SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter,
297
 
                                                                       &_geometry_helper, helical);
298
 
  //    SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter->GetFiltered(),
299
 
  //                                                             &_geometry_helper);
300
 
      evt.add_scifitrack(track);
301
 
    }
302
 
    for (size_t track_i = 0; track_i < number_straight_tracks; track_i++) {
303
 
      SciFiStraightPRTrack* straight = evt.straightprtracks().at(track_i);
304
 
      SciFiSpacePointPArray spacepoints = straight->get_spacepoints_pointers();
305
 
 
306
 
      Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
307
 
                                                                          _spacepoint_recon_plane);
308
 
      Kalman::State seed = ComputeSeed(straight, &_geometry_helper, _seed_value);
309
 
 
310
 
      std::cerr << Kalman::print_state(seed, "Helical Seed");
311
 
      std::cerr << "Straight Track\n";
312
 
      print_straight(straight);
313
 
 
314
 
      _spacepoint_straight_track_fitter->SetSeed(seed);
315
 
      _spacepoint_straight_track_fitter->SetData(data_track);
316
 
 
317
 
      _spacepoint_straight_track_fitter->Filter(false);
318
 
      _spacepoint_straight_track_fitter->Smooth(false);
319
 
 
320
 
      Kalman::Track predicted = _spacepoint_straight_track_fitter->GetPredicted();
321
 
      Kalman::Track filtered = _spacepoint_straight_track_fitter->GetFiltered();
322
 
      Kalman::Track smoothed = _spacepoint_straight_track_fitter->GetSmoothed();
323
 
 
324
 
      std::cerr << "Data Track\n";
325
 
      std::cerr << Kalman::print_track(data_track);
326
 
      std::cerr << "Predicted Track\n";
327
 
      std::cerr << Kalman::print_track(predicted);
328
 
      std::cerr << "Filtered Track\n";
329
 
      std::cerr << Kalman::print_track(filtered);
330
 
      std::cerr << "Smoothed Track\n";
331
 
      std::cerr << Kalman::print_track(smoothed);
332
 
 
333
 
      SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_straight_track_fitter,
334
 
                                                                      &_geometry_helper, straight);
335
 
      evt.add_scifitrack(track);
336
 
    }
337
 
  } else { // PATREC Not Active
338
 
    std::cerr << "PATREC NOT ACTIVE\n";
339
 
 
340
 
    SciFiSpacePointPArray event_spacepoints = evt.spacepoints();
341
 
    for (int tracker = 0; tracker < 2; ++tracker) {
342
 
 
343
 
      SciFiSpacePointPArray spacepoints = find_spacepoints(event_spacepoints, tracker);
344
 
      std::cerr << " USING " << spacepoints.size() << " SPACEPOINTS from TRACKER : "
345
 
                                                                                << tracker << "\n";
346
 
      if (spacepoints.size() != 5) break;
347
 
 
348
 
      int id = ((4*3 + 1) + _spacepoint_recon_plane) * (tracker == 0 ? -1 : 1);
349
 
      Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
350
 
                                                                          _spacepoint_recon_plane);
351
 
      Kalman::State seed(5);
352
 
      seed.SetId(id);
353
 
 
354
 
      bool seed_found = false;
355
 
      for (SciFiSpacePointPArray::iterator it_s = spacepoints.begin();
356
 
                                                               it_s != spacepoints.end(); ++it_s) {
357
 
        if ((*it_s)->get_station() == 5) {
358
 
 
359
 
          ThreeVector pos = (*it_s)->get_true_position();
360
 
          ThreeVector mom = (*it_s)->get_true_momentum();
361
 
          std::cerr << pos.x() << ", " << pos.y() << ", " << pos.z() << " ;  " << mom.x()
362
 
                                                     << ", " << mom.y() << ", " << mom.z() << "\n";
363
 
 
364
 
          TMatrixD state(5, 1);
365
 
          state(0, 0) = pos.x();
366
 
          state(1, 0) = mom.x();
367
 
          state(2, 0) = pos.y();
368
 
          state(3, 0) = mom.y();
369
 
          state(4, 0) = 1.0 / mom.z();
370
 
 
371
 
          seed.SetVector(state);
372
 
          seed.SetPosition(pos.z());
373
 
          seed_found = true;
374
 
          break;
375
 
        }
376
 
      }
377
 
      if (!seed_found) continue;
378
 
 
379
 
      TMatrixD seed_cov(5, 5);
380
 
      seed_cov(0, 0) = _seed_value;
381
 
      seed_cov(1, 1) = _seed_value;
382
 
      seed_cov(2, 2) = _seed_value;
383
 
      seed_cov(3, 3) = _seed_value;
384
 
      seed_cov(4, 4) = _seed_value;
385
 
      seed.SetCovariance(seed_cov);
386
 
 
387
 
      std::cerr << Kalman::print_state(seed, "Helical Seed");
388
 
 
389
 
      _spacepoint_helical_track_fitter->SetSeed(seed);
390
 
      _spacepoint_helical_track_fitter->SetData(data_track);
391
 
 
392
 
      _spacepoint_helical_track_fitter->Filter(false);
393
 
      _spacepoint_helical_track_fitter->Smooth(false);
394
 
 
395
 
      Kalman::Track predicted = _spacepoint_helical_track_fitter->GetPredicted();
396
 
      Kalman::Track filtered = _spacepoint_helical_track_fitter->GetFiltered();
397
 
      Kalman::Track smoothed = _spacepoint_helical_track_fitter->GetSmoothed();
398
 
      std::cerr << "Data Track\n";
399
 
      std::cerr << Kalman::print_track(data_track);
400
 
      std::cerr << "Measured Predicted\n";
401
 
      SciFiSpacepointMeasurements<5>* temp_measurement = new SciFiSpacepointMeasurements<5>(1.0);
402
 
      for (unsigned int k = 0; k < predicted.GetLength(); ++k) {
403
 
        Kalman::State temp_state = temp_measurement->Measure(predicted[k]);
404
 
        std::cerr << Kalman::print_state(temp_state);
405
 
      }
406
 
      std::cerr << "Predicted Track\n";
407
 
      std::cerr << Kalman::print_track(predicted);
408
 
      std::cerr << "Filtered Track\n";
409
 
      std::cerr << Kalman::print_track(filtered);
410
 
      std::cerr << "Smoothed Track\n";
411
 
      std::cerr << Kalman::print_track(smoothed);
412
 
      delete temp_measurement;
413
 
 
414
 
      SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter,
415
 
                                                                       &_geometry_helper, helical);
416
 
  //    SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter->GetFiltered(),
417
 
  //                                                             &_geometry_helper);
418
 
      evt.add_scifitrack(track);
419
 
    }
420
 
  }
421
 
#else
422
 
  if (_patrec_on) {
423
 
    size_t number_helical_tracks = evt.helicalprtracks().size();
424
 
    size_t number_straight_tracks = evt.straightprtracks().size();
425
 
    for (size_t track_i = 0; track_i < number_helical_tracks; track_i++) {
426
 
      SciFiHelicalPRTrack* helical = evt.helicalprtracks().at(track_i);
427
 
 
428
 
      Kalman::Track data_track = BuildTrack(helical, &_geometry_helper);
429
 
      Kalman::State seed = ComputeSeed(helical, &_geometry_helper, _use_eloss, _seed_value);
430
 
 
431
 
      _helical_track_fitter->SetData(data_track);
432
 
      _helical_track_fitter->SetSeed(seed);
433
 
 
434
 
      _helical_track_fitter->Filter(false);
435
 
      _helical_track_fitter->Smooth(false);
436
 
 
437
 
      SciFiTrack* track = ConvertToSciFiTrack(_helical_track_fitter, &_geometry_helper, helical);
438
 
      calculate_track_rating(track);
439
 
 
440
 
      evt.add_scifitrack(track);
441
 
    }
442
 
    for (size_t track_i = 0; track_i < number_straight_tracks; track_i++) {
443
 
      SciFiStraightPRTrack* straight = evt.straightprtracks().at(track_i);
444
 
      Kalman::Track data_track = BuildTrack(straight, &_geometry_helper);
445
 
 
446
 
      Kalman::State seed = ComputeSeed(straight, &_geometry_helper, _seed_value);
447
 
 
448
 
      _straight_track_fitter->SetData(data_track);
449
 
      _straight_track_fitter->SetSeed(seed);
450
 
 
451
 
      _straight_track_fitter->Filter(false);
452
 
      _straight_track_fitter->Smooth(false);
453
 
 
454
 
      SciFiTrack* track = ConvertToSciFiTrack(_straight_track_fitter, &_geometry_helper, straight);
455
 
      calculate_track_rating(track);
456
 
 
457
 
      evt.add_scifitrack(track);
458
 
    }
459
 
  } else { // PatRec Not Switched on. => Assume EVERYTHING is a single track.
460
 
    throw Exception(Exception::nonRecoverable,
461
 
        "No pattern recognition, required to seed Kalman",
462
 
        "MapCppTrackerRecon::track_fit()" );
463
 
 
464
 
    /*
465
 
    SciFiClusterPArray event_clusters = evt.clusters();
466
 
    for (int tracker = 0; tracker < 2; ++tracker) {
467
 
 
468
 
      SciFiClusterPArray clusters = find_clusters(event_clusters, tracker);
469
 
      if (clusters.size() != 15) break;
470
 
 
471
 
      Kalman::Track data_track = BuildTrack(clusters, &_geometry_helper);
472
 
      Kalman::State seed(5);
473
 
 
474
 
      bool seed_found = false;
475
 
      for (SciFiClusterPArray::iterator it_c = clusters.begin(); it_c != clusters.end(); ++it_c) {
476
 
        if (((*it_c)->get_station() == 5) && ((*it_c)->get_plane() == 2)) {
477
 
 
478
 
          ThreeVector pos = (*it_c)->get_true_position();
479
 
          ThreeVector mom = (*it_c)->get_true_momentum();
480
 
 
481
 
          TMatrixD state(5, 1);
482
 
          state(0, 0) = pos.x();
483
 
          state(1, 0) = mom.x();
484
 
          state(2, 0) = pos.y();
485
 
          state(3, 0) = mom.y();
486
 
          state(4, 0) = 1.0 / mom.z();
487
 
 
488
 
          seed.SetVector(state);
489
 
          seed.SetPosition(pos.z());
490
 
          seed_found = true;
491
 
          break;
492
 
        }
493
 
      }
494
 
      if (!seed_found) continue;
495
 
 
496
 
      TMatrixD seed_cov(5, 5);
497
 
      seed_cov(0, 0) = _seed_value;
498
 
      seed_cov(1, 1) = _seed_value;
499
 
      seed_cov(2, 2) = _seed_value;
500
 
      seed_cov(3, 3) = _seed_value;
501
 
      seed_cov(4, 4) = _seed_value;
502
 
      seed.SetCovariance(seed_cov);
503
 
 
504
 
      _helical_track_fitter->SetData(data_track);
505
 
      _helical_track_fitter->SetSeed(seed);
506
 
 
507
 
      _helical_track_fitter->Filter(false);
508
 
      _helical_track_fitter->Smooth(false);
509
 
 
510
 
      SciFiTrack* track = ConvertToSciFiTrack(_helical_track_fitter, &_geometry_helper);
511
 
 
512
 
      evt.add_scifitrack(track);
513
 
    }
514
 
    */
515
 
  }
516
 
#endif
 
239
  size_t number_helical_tracks = evt.helicalprtracks().size();
 
240
  size_t number_straight_tracks = evt.straightprtracks().size();
 
241
  for (size_t track_i = 0; track_i < number_helical_tracks; track_i++) {
 
242
    SciFiHelicalPRTrack* helical = evt.helicalprtracks().at(track_i);
 
243
 
 
244
    Kalman::Track data_track = BuildTrack(helical, &_geometry_helper);
 
245
    Kalman::State seed = ComputeSeed(helical, &_geometry_helper, _use_eloss, _seed_value);
 
246
 
 
247
    _helical_track_fitter->SetData(data_track);
 
248
    _helical_track_fitter->SetSeed(seed);
 
249
 
 
250
    _helical_track_fitter->Filter(false);
 
251
    _helical_track_fitter->Smooth(false);
 
252
 
 
253
    SciFiTrack* track = ConvertToSciFiTrack(_helical_track_fitter, &_geometry_helper, helical);
 
254
    calculate_track_rating(track);
 
255
 
 
256
    evt.add_scifitrack(track);
 
257
  }
 
258
  for (size_t track_i = 0; track_i < number_straight_tracks; track_i++) {
 
259
    SciFiStraightPRTrack* straight = evt.straightprtracks().at(track_i);
 
260
    Kalman::Track data_track = BuildTrack(straight, &_geometry_helper);
 
261
 
 
262
    Kalman::State seed = ComputeSeed(straight, &_geometry_helper, _seed_value);
 
263
 
 
264
    _straight_track_fitter->SetData(data_track);
 
265
    _straight_track_fitter->SetSeed(seed);
 
266
 
 
267
    _straight_track_fitter->Filter(false);
 
268
    _straight_track_fitter->Smooth(false);
 
269
 
 
270
    SciFiTrack* track = ConvertToSciFiTrack(_straight_track_fitter, &_geometry_helper, straight);
 
271
    calculate_track_rating(track);
 
272
 
 
273
    evt.add_scifitrack(track);
 
274
  }
517
275
}
518
276
 
 
277
 
519
278
void MapCppTrackerRecon::extrapolate_helical_reference(SciFiEvent& event) const {
520
279
  SciFiHelicalPRTrackPArray helicals = event.helicalprtracks();
521
280
  size_t num_tracks = helicals.size();
528
287
    int tracker = track->get_tracker();
529
288
    ThreeVector reference = _geometry_helper.GetReferencePosition(tracker);
530
289
    CLHEP::HepRotation rotation = _geometry_helper.GetReferenceRotation(tracker);
531
 
    rotation.invert();
532
290
 
533
291
    double r  = track->get_R();
534
292
    double pt = - track->get_charge()*CLHEP::c_light*_geometry_helper.GetFieldValue(tracker)*r;
554
312
  }
555
313
}
556
314
 
 
315
 
557
316
void MapCppTrackerRecon::extrapolate_straight_reference(SciFiEvent& event) const {
558
317
  SciFiStraightPRTrackPArray straights = event.straightprtracks();
559
318
  size_t num_tracks = straights.size();
567
326
    int tracker = track->get_tracker();
568
327
    ThreeVector reference = _geometry_helper.GetReferencePosition(tracker);
569
328
    CLHEP::HepRotation rotation = _geometry_helper.GetReferenceRotation(tracker);
570
 
//    rotation.invert();
571
329
 
572
330
    pos.setX(track->get_x0());
573
331
    pos.setY(track->get_y0());
581
339
    mom.setY(my*default_mom);
582
340
    mom *= rotation;
583
341
    if (tracker == 0) mom *= -1.0;
584
 
//    mom.setZ(sqrt(40000.0 - mom[0]*mom[0] + mom[1]*mom[1]));
585
342
    mom.setZ(default_mom);
586
343
 
587
344
    track->set_reference_position(pos);
589
346
  }
590
347
}
591
348
 
 
349
 
 
350
void MapCppTrackerRecon::set_spacepoint_global_output(const SciFiSpacePointPArray& spoints) const {
 
351
  for (auto sp : spoints) {
 
352
    sp->set_global_position(
 
353
      _geometry_helper.TransformPositionToGlobal(sp->get_position(), sp->get_tracker()));
 
354
  }
 
355
}
 
356
 
 
357
 
 
358
void MapCppTrackerRecon::set_straight_prtrack_global_output(
 
359
  const SciFiStraightPRTrackPArray& trks) const {
 
360
  for (auto trk : trks) {
 
361
    ThreeVector pos = trk->get_spacepoints_pointers()[0]->get_position();
 
362
    std::vector<double> params{ trk->get_x0(), trk->get_mx(), trk->get_y0(), trk->get_my() }; // NOLINT
 
363
    std::vector<double> global_params =
 
364
      _geometry_helper.TransformStraightParamsToGlobal(params, trk->get_tracker());
 
365
 
 
366
    // Update the track
 
367
    trk->set_global_x0(global_params[0]);
 
368
    trk->set_global_mx(global_params[1]);
 
369
    trk->set_global_y0(global_params[2]);
 
370
    trk->set_global_my(global_params[3]);
 
371
  }
 
372
}
 
373
 
 
374
 
592
375
void MapCppTrackerRecon::calculate_track_rating(SciFiTrack* track) const {
593
376
  SciFiBasePRTrack* pr_track = track->pr_track_pointer();
594
377
  int number_spacepoints = pr_track->get_num_points();
665
448
    std::cerr << '\n';
666
449
  }
667
450
  std::cerr << std::endl;
668
 
 
669
 
 
670
 
 
671
 
 
672
 
  /*
673
 
  std::cerr << "Helix Pz Recon:\n";
674
 
  SciFiHelicalPRTrackPArray helices = event.helicalprtracks();
675
 
  for (unsigned int i = 0; i < helices.size(); ++i)
676
 
  {
677
 
    double radius = helices[i]->get_R();
678
 
    double patrec_pt = 1.19*radius;
679
 
    double patrec_pz = patrec_pt / helices[i]->get_dsdz();
680
 
    std::cerr << "Pz = " << patrec_pz << '\n';
681
 
 
682
 
    DoubleArray phis = helices[i]->get_phi();
683
 
    double start_phi = phis[0];
684
 
    double end_phi = phis[ phis.size()-1 ];
685
 
 
686
 
    double delta_phi = end_phi - start_phi;
687
 
 
688
 
    double Z = (1100.0 * 2.0 * CLHEP::pi) / delta_phi;
689
 
    double pz = (1.19 * Z) / (2.0 * CLHEP::pi);
690
 
 
691
 
    std::cerr << "Pz = " << pz << "\n\n";
692
 
 
693
 
    for (unsigned int j = 0; j < phis.size(); ++j)
694
 
    {
695
 
      std::cerr << phis[j] << '\n';
696
 
    }
697
 
  }
698
 
  std::cerr << std::endl;
699
 
  */
700
 
 
701
 
  /*
702
 
  Squeak::mout(Squeak::info) << event.digits().size() << " "
703
 
                              << event.clusters().size() << " "
704
 
                              << event.spacepoints().size() << "; "
705
 
                              << event.straightprtracks().size() << " "
706
 
                              << event.helicalprtracks().size() << "; ";
707
 
  for (size_t track_i = 0; track_i < event.scifitracks().size(); track_i++) {
708
 
    Squeak::mout(Squeak::info) << " Chi2: " << event.scifitracks()[track_i]->f_chi2() << "; "
709
 
                 << " Chi2: " << event.scifitracks()[track_i]->s_chi2() << "; "
710
 
                 << " P-Value: " << event.scifitracks()[track_i]->P_value() << "; ";
711
 
  }
712
 
  Squeak::mout(Squeak::info) << std::endl;
713
 
  */
714
 
}
715
 
 
716
 
void print_helical(SciFiHelicalPRTrack* helical) {
717
 
  SciFiSpacePointPArray spacepoints = helical->get_spacepoints_pointers();
718
 
  size_t numb_spacepoints = spacepoints.size();
719
 
 
720
 
  std::cerr << "HELICAL TRACK:\n";
721
 
 
722
 
  double r = helical->get_R();
723
 
  double sz = helical->get_line_sz_c();
724
 
  double charge = helical->get_charge();
725
 
  std::cerr << "Radius = " << r << " Sz = " << sz << " Charge = " << charge
726
 
            << " Phi_0 = " << sz/r << "\n";
727
 
 
728
 
  for (size_t i = 0; i < numb_spacepoints; ++i) {
729
 
    SciFiSpacePoint *spacepoint = spacepoints[i];
730
 
 
731
 
    std::cerr << spacepoint->get_tracker() << ", "
732
 
              << spacepoint->get_station() << "  -  "
733
 
              << spacepoint->get_position().x() << ", "
734
 
              << spacepoint->get_position().y() << ", "
735
 
              << spacepoint->get_position().z() << "\n";
736
 
  }
737
 
 
738
 
  std::cerr << std::endl;
739
 
}
740
 
 
741
 
void print_straight(SciFiStraightPRTrack* straight) {
742
 
  SciFiSpacePointPArray spacepoints = straight->get_spacepoints_pointers();
743
 
  size_t numb_spacepoints = spacepoints.size();
744
 
 
745
 
  std::cerr << "STRAIGHT TRACK:\n";
746
 
  for (size_t i = 0; i < numb_spacepoints; ++i) {
747
 
    SciFiSpacePoint *spacepoint = spacepoints[i];
748
 
 
749
 
    std::cerr << spacepoint->get_tracker() << ", "
750
 
              << spacepoint->get_station() << "  -  "
751
 
              << spacepoint->get_position().x() << ", "
752
 
              << spacepoint->get_position().y() << ", "
753
 
              << spacepoint->get_position().z() << "\n";
754
 
  }
755
 
 
756
 
  std::cerr << std::endl;
757
 
}
758
 
 
759
 
SciFiClusterPArray find_clusters(SciFiClusterPArray cluster_array, int tracker) {
760
 
  SciFiClusterPArray found_clusters;
761
 
 
762
 
  for (SciFiClusterPArray::iterator cl_it = cluster_array.begin();
763
 
                                                           cl_it != cluster_array.end(); ++cl_it) {
764
 
    if ((*cl_it)->get_tracker() == tracker) {
765
 
      found_clusters.push_back((*cl_it));
766
 
    }
767
 
  }
768
 
 
769
 
  return found_clusters;
770
 
}
771
 
 
772
 
SciFiSpacePointPArray find_spacepoints(SciFiSpacePointPArray spacepoint_array, int tracker) {
773
 
  SciFiSpacePointPArray found_spacepoints;
774
 
 
775
 
  for (SciFiSpacePointPArray::iterator sp_it = spacepoint_array.begin();
776
 
                                                        sp_it != spacepoint_array.end(); ++sp_it) {
777
 
    if ((*sp_it)->get_tracker() == tracker) {
778
 
      found_spacepoints.push_back((*sp_it));
779
 
    }
780
 
  }
781
 
 
782
 
  return found_spacepoints;
783
 
}
784
 
 
785
 
void MapCppTrackerRecon::set_spacepoint_global_output(SciFiSpacePointPArray spoints) const {
786
 
  for (auto sp : spoints) {
787
 
    int tracker_id = sp->get_tracker();
788
 
 
789
 
    ThreeVector reference_position = _geometry_helper.GetReferencePosition(tracker_id);
790
 
    CLHEP::HepRotation reference_rotation = _geometry_helper.GetReferenceRotation(tracker_id);
791
 
 
792
 
    ThreeVector global_position = sp->get_position();
793
 
    global_position *= reference_rotation;
794
 
    global_position += reference_position;
795
 
 
796
 
    sp->set_global_position(global_position);
797
 
  }
798
 
}
799
 
 
800
 
/*
801
 
// Find the plane id of the middle plane for this station
802
 
int plane_id = 0;
803
 
switch ( station_id ) {
804
 
  case 1:
805
 
    plane_id = 2;
806
 
    break;
807
 
  case 2:
808
 
    plane_id = 5;
809
 
    break;
810
 
  case 3:
811
 
    plane_id = 8;
812
 
    break;
813
 
  case 4:
814
 
    plane_id = 11;
815
 
    break;
816
 
  case 5:
817
 
    plane_id = 14;
818
 
    break;
819
 
  default:
820
 
    plane_id = 0;
821
 
    std::cerr << "WARNING: MapCppTrackerRecon: Failed to find spacepoint plane_id";
822
 
    break;
823
 
}
824
 
 
825
 
// Use the geometry helper to find the global position of this plane
826
 
ThreeVector plane_global_position = _geometry_helper.GeometryMap().find(tracker_id)->second.Planes.find(plane_id)->second.GlobalPosition;
827
 
*/
 
451
}
828
452
 
829
453
} // ~namespace MAUS