← Back to branch summary

~vpec/maus/tof_calib_read

« back to all changes in this revision

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

  • Committer: Adam Dobbs
  • Date: 2015-06-26 14:33:50 UTC
  • mfrom: (659.1.113 release-candidate)
  • Revision ID: phuccj@gmail.com-20150626143350-m56pbthi31ahqvxj
Tags: MAUS-v0.9.6
MAUS-v0.9.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerRecon/MapCppTrackerRecon.cc
15
15
 *
16
16
 */
17
17
 
 
18
#include "src/map/MapCppTrackerRecon/MapCppTrackerRecon.hh"
 
19
 
 
20
#include <sstream>
 
21
 
18
22
#include "src/common_cpp/API/PyWrapMapBase.hh"
19
 
#include "src/map/MapCppTrackerRecon/MapCppTrackerRecon.hh"
 
23
#include "src/common_cpp/Recon/Kalman/MAUSTrackWrapper.hh"
20
24
 
21
25
namespace MAUS {
 
26
 
 
27
// void print_state(Kalman::State& state);
 
28
 
 
29
void print_helical(SciFiHelicalPRTrack* helical);
 
30
void print_straight(SciFiStraightPRTrack* straight);
 
31
 
 
32
SciFiClusterPArray find_clusters(SciFiClusterPArray cluster_array, int tracker);
 
33
SciFiSpacePointPArray find_spacepoints(SciFiSpacePointPArray spacepoint_array, int tracker);
 
34
 
 
35
 
22
36
PyMODINIT_FUNC init_MapCppTrackerRecon(void) {
23
37
  PyWrapMapBase<MapCppTrackerRecon>::PyWrapMapBaseModInit
24
38
                                        ("MapCppTrackerRecon", "", "", "", "");
34
48
  if (!Globals::HasInstance()) {
35
49
    GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
36
50
  }
37
 
  Json::Value *json = Globals::GetConfigurationCards();
38
 
  _helical_pr_on  = (*json)["SciFiPRHelicalOn"].asBool();
39
 
  _straight_pr_on = (*json)["SciFiPRStraightOn"].asBool();
40
 
  _kalman_on      = (*json)["SciFiKalmanOn"].asBool();
41
 
  _size_exception = (*json)["SciFiClustExcept"].asInt();
42
 
  _min_npe        = (*json)["SciFiNPECut"].asDouble();
 
51
  Json::Value* json = Globals::GetConfigurationCards();
 
52
  _helical_pr_on   = (*json)["SciFiPRHelicalOn"].asBool();
 
53
  _straight_pr_on  = (*json)["SciFiPRStraightOn"].asBool();
 
54
  _kalman_on       = (*json)["SciFiKalmanOn"].asBool();
 
55
  _patrec_on       = (*json)["SciFiPatRecOn"].asBool();
 
56
  _size_exception  = (*json)["SciFiClustExcept"].asInt();
 
57
  _min_npe         = (*json)["SciFiNPECut"].asDouble();
 
58
  _use_mcs         = (*json)["SciFiKalman_use_MCS"].asBool();
 
59
  _use_eloss       = (*json)["SciFiKalman_use_Eloss"].asBool();
 
60
  _use_patrec_seed = (*json)["SciFiSeedPatRec"].asBool();
 
61
  _correct_pz      = (*json)["SciFiKalmanCorrectPz"].asBool();
43
62
 
44
63
  MiceModule* module = Globals::GetReconstructionMiceModules();
45
64
  std::vector<const MiceModule*> modules =
46
65
    module->findModulesByPropertyString("SensitiveDetector", "SciFi");
47
66
  _geometry_helper = SciFiGeometryHelper(modules);
48
67
  _geometry_helper.Build();
 
68
 
 
69
 
 
70
  _cluster_recon = SciFiClusterRec(_size_exception, _min_npe, _geometry_helper.GeometryMap());
 
71
 
 
72
  _spacepoint_recon = SciFiSpacePointRec();
 
73
 
 
74
  _pattern_recognition.LoadGlobals();
 
75
  _pattern_recognition.set_helical_pr_on(_helical_pr_on);
 
76
  _pattern_recognition.set_straight_pr_on(_straight_pr_on);
 
77
  _pattern_recognition.set_bz_t1(_geometry_helper.GetFieldValue(0));
 
78
  _pattern_recognition.set_bz_t2(_geometry_helper.GetFieldValue(1));
 
79
 
 
80
  if (_use_patrec_seed) {
 
81
    _seed_value = -1.0;
 
82
  } else {
 
83
    _seed_value = (*json)["SciFiSeedCovariance"].asDouble();
 
84
  }
 
85
 
 
86
#ifdef KALMAN_TEST
 
87
  HelicalPropagator* spacepoint_helical_prop = new HelicalPropagator(&_geometry_helper);
 
88
  spacepoint_helical_prop->SetCorrectPz(_correct_pz);
 
89
  spacepoint_helical_prop->SetIncludeMCS(_use_mcs);
 
90
  spacepoint_helical_prop->SetSubtractELoss(_use_eloss);
 
91
  _spacepoint_helical_track_fitter = new Kalman::TrackFit(spacepoint_helical_prop,
 
92
      new SciFiSpacepointMeasurements<5>(0.1));
 
93
 
 
94
  StraightPropagator* spacepoint_straight_prop = new StraightPropagator(&_geometry_helper);
 
95
  spacepoint_straight_prop->SetIncludeMCS(_use_mcs);
 
96
  _spacepoint_straight_track_fitter = new Kalman::TrackFit(spacepoint_straight_prop,
 
97
      new SciFiSpacepointMeasurements<4>(0.1));
 
98
 
 
99
  _spacepoint_recon_plane = 0;
 
100
#else
 
101
  HelicalPropagator* helical_prop = new HelicalPropagator(&_geometry_helper);
 
102
  helical_prop->SetCorrectPz(_correct_pz);
 
103
  helical_prop->SetIncludeMCS(_use_mcs);
 
104
  helical_prop->SetSubtractELoss(_use_eloss);
 
105
  _helical_measurement = new SciFiHelicalMeasurements(&_geometry_helper);
 
106
  _helical_track_fitter = new Kalman::TrackFit(helical_prop, _helical_measurement);
 
107
 
 
108
  StraightPropagator* straight_prop = new StraightPropagator(&_geometry_helper);
 
109
  straight_prop->SetIncludeMCS(_use_mcs);
 
110
  _straight_measurement = new SciFiStraightMeasurements(&_geometry_helper);
 
111
  _straight_track_fitter = new Kalman::TrackFit(straight_prop, _straight_measurement);
 
112
#endif
49
113
}
50
114
 
51
115
void MapCppTrackerRecon::_death() {
54
118
void MapCppTrackerRecon::_process(Data* data) const {
55
119
  Spill& spill = *(data->GetSpill());
56
120
 
57
 
  if ( spill.GetReconEvents() ) {
58
 
    for ( unsigned int k = 0; k < spill.GetReconEvents()->size(); k++ ) {
 
121
  if (spill.GetReconEvents()) {
 
122
    for (unsigned int k = 0; k < spill.GetReconEvents()->size(); k++) {
59
123
      SciFiEvent *event = spill.GetReconEvents()->at(k)->GetSciFiEvent();
 
124
      if (!event) {
 
125
        continue;
 
126
      }
60
127
 
61
128
      // Build Clusters.
62
 
      if ( event->digits().size() ) {
63
 
        cluster_recon(*event);
 
129
      if (event->digits().size()) {
 
130
        _cluster_recon.process(*event);
64
131
      }
65
132
      // Build SpacePoints.
66
 
      if ( event->clusters().size() ) {
67
 
        spacepoint_recon(*event);
 
133
      if (event->clusters().size()) {
 
134
        _spacepoint_recon.process(*event);
68
135
      }
69
136
      // Pattern Recognition.
70
 
      if ( event->spacepoints().size() ) {
71
 
        pattern_recognition(*event);
 
137
      if (_patrec_on && event->spacepoints().size()) {
 
138
        _pattern_recognition.process(*event);
 
139
        extrapolate_helical_reference(*event);
 
140
        extrapolate_straight_reference(*event);
72
141
      }
73
142
      // Kalman Track Fit.
74
 
      if ( _kalman_on ) {
75
 
        if ( event->straightprtracks().size() || event->helicalprtracks().size() ) {
76
 
          track_fit(*event);
77
 
        }
 
143
      if (_kalman_on) {
 
144
        track_fit(*event);
78
145
      }
79
 
      // print_event_info(*event);
 
146
//      print_event_info(*event);
80
147
    }
81
148
  } else {
82
149
    std::cout << "No recon events found\n";
83
150
  }
84
151
}
85
152
 
86
 
void MapCppTrackerRecon::cluster_recon(SciFiEvent &evt) const {
87
 
  SciFiClusterRec clustering(_size_exception, _min_npe, _geometry_helper.GeometryMap());
88
 
  clustering.process(evt);
89
 
}
90
 
 
91
 
void MapCppTrackerRecon::spacepoint_recon(SciFiEvent &evt) const {
92
 
  SciFiSpacePointRec spacepoints;
93
 
  spacepoints.process(evt);
94
 
}
95
 
 
96
 
void MapCppTrackerRecon::pattern_recognition(SciFiEvent &evt) const {
97
 
  PatternRecognition pr1; // Pat rec constructor calls Globals again
98
 
  // We let the Map's helical and straight flags overide the interal pat rec variables for these,
99
 
  // (pulled by pat rec from Globals) as this way we can customise which type runs for
100
 
  // testing purposes.
101
 
  pr1.set_helical_pr_on(_helical_pr_on);
102
 
  pr1.set_straight_pr_on(_straight_pr_on);
103
 
  pr1.process(evt);
104
 
}
105
 
 
106
153
void MapCppTrackerRecon::track_fit(SciFiEvent &evt) const {
107
 
  std::vector<KalmanSeed*> seeds;
108
 
  size_t number_helical_tracks  = evt.helicalprtracks().size();
109
 
  size_t number_straight_tracks = evt.straightprtracks().size();
110
 
 
111
 
  for ( size_t track_i = 0; track_i < number_helical_tracks; track_i++ ) {
112
 
    int tracker = evt.helicalprtracks()[track_i]->get_tracker();
113
 
    double Bz = _geometry_helper.GetFieldValue(tracker);
114
 
    KalmanSeed *seed = new KalmanSeed(_geometry_helper.GeometryMap());
115
 
    seed->SetField(Bz);
116
 
    seed->Build<SciFiHelicalPRTrack>(evt.helicalprtracks()[track_i]);
117
 
    seeds.push_back(seed);
118
 
  }
119
 
 
120
 
  for ( size_t track_i = 0; track_i < number_straight_tracks; track_i++ ) {
121
 
    KalmanSeed *seed = new KalmanSeed(_geometry_helper.GeometryMap());
122
 
    seed->Build<SciFiStraightPRTrack>(evt.straightprtracks()[track_i]);
123
 
    seeds.push_back(seed);
124
 
  }
125
 
 
126
 
  if ( seeds.size() ) {
127
 
    KalmanTrackFit fit;
128
 
    fit.SaveGeometry(_geometry_helper.RefPos(), _geometry_helper.Rot());
129
 
    fit.Process(seeds, evt);
130
 
  }
131
 
}
132
 
 
133
 
// void MapCppTrackerRecon::print_event_info(SciFiEvent &event) {
134
 
//   std::cerr << event.digits().size() << " "
135
 
//                               << event.clusters().size() << " "
136
 
//                               << event.spacepoints().size() << "; "
137
 
//                               << event.straightprtracks().size() << " "
138
 
//                               << event.helicalprtracks().size() << "; ";
139
 
//   for ( size_t track_i = 0; track_i < event.scifitracks().size(); track_i++ ) {
140
 
//     std::cerr << " Chi2: " << event.scifitracks()[track_i]->f_chi2() << "; "
141
 
//                            << " Chi2: " << event.scifitracks()[track_i]->s_chi2() << "; "
142
 
//                            << " P-Value: " << event.scifitracks()[track_i]->P_value() << "; ";
143
 
//   }
144
 
//   std::cerr << std::endl;
145
 
//   /*
146
 
//   Squeak::mout(Squeak::info) << event.digits().size() << " "
147
 
//                               << event.clusters().size() << " "
148
 
//                               << event.spacepoints().size() << "; "
149
 
//                               << event.straightprtracks().size() << " "
150
 
//                               << event.helicalprtracks().size() << "; ";
151
 
//   for ( size_t track_i = 0; track_i < event.scifitracks().size(); track_i++ ) {
152
 
//     Squeak::mout(Squeak::info) << " Chi2: " << event.scifitracks()[track_i]->f_chi2() << "; "
153
 
//                  << " Chi2: " << event.scifitracks()[track_i]->s_chi2() << "; "
154
 
//                  << " P-Value: " << event.scifitracks()[track_i]->P_value() << "; ";
155
 
//   }
156
 
//   Squeak::mout(Squeak::info) << std::endl;
157
 
//   */
158
 
// }
159
 
 
 
154
#ifdef KALMAN_TEST
 
155
  std::cerr << "KALMAN_TEST ACTIVE\n";
 
156
  if (_patrec_on) {
 
157
    std::cerr << "PATREC ACTIVE\n";
 
158
    size_t number_helical_tracks = evt.helicalprtracks().size();
 
159
    size_t number_straight_tracks = evt.straightprtracks().size();
 
160
 
 
161
    for (size_t track_i = 0; track_i < number_helical_tracks; track_i++) {
 
162
      SciFiHelicalPRTrack* helical = evt.helicalprtracks().at(track_i);
 
163
      SciFiSpacePointPArray spacepoints = helical->get_spacepoints_pointers();
 
164
 
 
165
      Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
 
166
                                                                          _spacepoint_recon_plane);
 
167
 
 
168
      Kalman::State seed = ComputeSeed(helical, &_geometry_helper, false, _seed_value);
 
169
 
 
170
      std::cerr << Kalman::print_state(seed, "Helical Seed");
 
171
      seed.GetCovariance().Print();
 
172
 
 
173
      std::cerr << "Helical Track\n";
 
174
      print_helical(helical);
 
175
 
 
176
      _spacepoint_helical_track_fitter->SetSeed(seed);
 
177
      _spacepoint_helical_track_fitter->SetData(data_track);
 
178
 
 
179
      _spacepoint_helical_track_fitter->Filter(false);
 
180
      _spacepoint_helical_track_fitter->Smooth(false);
 
181
 
 
182
      Kalman::Track predicted = _spacepoint_helical_track_fitter->GetPredicted();
 
183
      Kalman::Track filtered = _spacepoint_helical_track_fitter->GetFiltered();
 
184
      Kalman::Track smoothed = _spacepoint_helical_track_fitter->GetSmoothed();
 
185
      std::cerr << "Data Track\n";
 
186
      std::cerr << Kalman::print_track(data_track);
 
187
      std::cerr << "Measured Predicted\n";
 
188
      SciFiSpacepointMeasurements<5>* temp_measurement = new SciFiSpacepointMeasurements<5>(1.0);
 
189
      for (unsigned int k = 0; k < predicted.GetLength(); ++k) {
 
190
        Kalman::State temp_state = temp_measurement->Measure(predicted[k]);
 
191
        std::cerr << Kalman::print_state(temp_state);
 
192
      }
 
193
      std::cerr << "Predicted Track\n";
 
194
      std::cerr << Kalman::print_track(predicted);
 
195
      std::cerr << "Filtered Track\n";
 
196
      std::cerr << Kalman::print_track(filtered);
 
197
      std::cerr << "Smoothed Track\n";
 
198
      std::cerr << Kalman::print_track(smoothed);
 
199
      delete temp_measurement;
 
200
 
 
201
      SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter, &_geometry_helper);
 
202
  //    SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter->GetFiltered(),
 
203
  //                                                             &_geometry_helper);
 
204
      evt.add_scifitrack(track);
 
205
    }
 
206
    for (size_t track_i = 0; track_i < number_straight_tracks; track_i++) {
 
207
      SciFiStraightPRTrack* straight = evt.straightprtracks().at(track_i);
 
208
      SciFiSpacePointPArray spacepoints = straight->get_spacepoints_pointers();
 
209
 
 
210
      Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
 
211
                                                                          _spacepoint_recon_plane);
 
212
      Kalman::State seed = ComputeSeed(straight, &_geometry_helper, _seed_value);
 
213
 
 
214
      std::cerr << Kalman::print_state(seed, "Helical Seed");
 
215
      std::cerr << "Straight Track\n";
 
216
      print_straight(straight);
 
217
 
 
218
      _spacepoint_straight_track_fitter->SetSeed(seed);
 
219
      _spacepoint_straight_track_fitter->SetData(data_track);
 
220
 
 
221
      _spacepoint_straight_track_fitter->Filter(false);
 
222
      _spacepoint_straight_track_fitter->Smooth(false);
 
223
 
 
224
      Kalman::Track predicted = _spacepoint_straight_track_fitter->GetPredicted();
 
225
      Kalman::Track filtered = _spacepoint_straight_track_fitter->GetFiltered();
 
226
      Kalman::Track smoothed = _spacepoint_straight_track_fitter->GetSmoothed();
 
227
 
 
228
      std::cerr << "Data Track\n";
 
229
      std::cerr << Kalman::print_track(data_track);
 
230
      std::cerr << "Predicted Track\n";
 
231
      std::cerr << Kalman::print_track(predicted);
 
232
      std::cerr << "Filtered Track\n";
 
233
      std::cerr << Kalman::print_track(filtered);
 
234
      std::cerr << "Smoothed Track\n";
 
235
      std::cerr << Kalman::print_track(smoothed);
 
236
 
 
237
      SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_straight_track_fitter,
 
238
                                                                                &_geometry_helper);
 
239
      evt.add_scifitrack(track);
 
240
    }
 
241
  } else { // PATREC Not Active
 
242
    std::cerr << "PATREC NOT ACTIVE\n";
 
243
 
 
244
    SciFiSpacePointPArray event_spacepoints = evt.spacepoints();
 
245
    for (int tracker = 0; tracker < 2; ++tracker) {
 
246
 
 
247
      SciFiSpacePointPArray spacepoints = find_spacepoints(event_spacepoints, tracker);
 
248
      std::cerr << " USING " << spacepoints.size() << " SPACEPOINTS from TRACKER : "
 
249
                                                                                << tracker << "\n";
 
250
      if (spacepoints.size() != 5) break;
 
251
 
 
252
      int id = ((4*3 + 1) + _spacepoint_recon_plane) * (tracker == 0 ? -1 : 1);
 
253
      Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
 
254
                                                                          _spacepoint_recon_plane);
 
255
      Kalman::State seed(5);
 
256
      seed.SetId(id);
 
257
 
 
258
      bool seed_found = false;
 
259
      for (SciFiSpacePointPArray::iterator it_s = spacepoints.begin();
 
260
                                                               it_s != spacepoints.end(); ++it_s) {
 
261
        if ((*it_s)->get_station() == 5) {
 
262
 
 
263
          ThreeVector pos = (*it_s)->get_true_position();
 
264
          ThreeVector mom = (*it_s)->get_true_momentum();
 
265
          std::cerr << pos.x() << ", " << pos.y() << ", " << pos.z() << " ;  " << mom.x()
 
266
                                                     << ", " << mom.y() << ", " << mom.z() << "\n";
 
267
 
 
268
          TMatrixD state(5, 1);
 
269
          state(0, 0) = pos.x();
 
270
          state(1, 0) = mom.x();
 
271
          state(2, 0) = pos.y();
 
272
          state(3, 0) = mom.y();
 
273
          state(4, 0) = 1.0 / mom.z();
 
274
 
 
275
          seed.SetVector(state);
 
276
          seed.SetPosition(pos.z());
 
277
          seed_found = true;
 
278
          break;
 
279
        }
 
280
      }
 
281
      if (!seed_found) continue;
 
282
 
 
283
      TMatrixD seed_cov(5, 5);
 
284
      seed_cov(0, 0) = _seed_value;
 
285
      seed_cov(1, 1) = _seed_value;
 
286
      seed_cov(2, 2) = _seed_value;
 
287
      seed_cov(3, 3) = _seed_value;
 
288
      seed_cov(4, 4) = _seed_value;
 
289
      seed.SetCovariance(seed_cov);
 
290
 
 
291
      std::cerr << Kalman::print_state(seed, "Helical Seed");
 
292
 
 
293
      _spacepoint_helical_track_fitter->SetSeed(seed);
 
294
      _spacepoint_helical_track_fitter->SetData(data_track);
 
295
 
 
296
      _spacepoint_helical_track_fitter->Filter(false);
 
297
      _spacepoint_helical_track_fitter->Smooth(false);
 
298
 
 
299
      Kalman::Track predicted = _spacepoint_helical_track_fitter->GetPredicted();
 
300
      Kalman::Track filtered = _spacepoint_helical_track_fitter->GetFiltered();
 
301
      Kalman::Track smoothed = _spacepoint_helical_track_fitter->GetSmoothed();
 
302
      std::cerr << "Data Track\n";
 
303
      std::cerr << Kalman::print_track(data_track);
 
304
      std::cerr << "Measured Predicted\n";
 
305
      SciFiSpacepointMeasurements<5>* temp_measurement = new SciFiSpacepointMeasurements<5>(1.0);
 
306
      for (unsigned int k = 0; k < predicted.GetLength(); ++k) {
 
307
        Kalman::State temp_state = temp_measurement->Measure(predicted[k]);
 
308
        std::cerr << Kalman::print_state(temp_state);
 
309
      }
 
310
      std::cerr << "Predicted Track\n";
 
311
      std::cerr << Kalman::print_track(predicted);
 
312
      std::cerr << "Filtered Track\n";
 
313
      std::cerr << Kalman::print_track(filtered);
 
314
      std::cerr << "Smoothed Track\n";
 
315
      std::cerr << Kalman::print_track(smoothed);
 
316
      delete temp_measurement;
 
317
 
 
318
      SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter, &_geometry_helper);
 
319
  //    SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter->GetFiltered(),
 
320
  //                                                             &_geometry_helper);
 
321
      evt.add_scifitrack(track);
 
322
    }
 
323
  }
 
324
#else
 
325
  if (_patrec_on) {
 
326
    size_t number_helical_tracks = evt.helicalprtracks().size();
 
327
    size_t number_straight_tracks = evt.straightprtracks().size();
 
328
    for (size_t track_i = 0; track_i < number_helical_tracks; track_i++) {
 
329
      SciFiHelicalPRTrack* helical = evt.helicalprtracks().at(track_i);
 
330
 
 
331
      Kalman::Track data_track = BuildTrack(helical, &_geometry_helper);
 
332
      Kalman::State seed = ComputeSeed(helical, &_geometry_helper, _use_eloss, _seed_value);
 
333
 
 
334
      _helical_track_fitter->SetData(data_track);
 
335
      _helical_track_fitter->SetSeed(seed);
 
336
 
 
337
      _helical_track_fitter->Filter(false);
 
338
      _helical_track_fitter->Smooth(false);
 
339
 
 
340
      SciFiTrack* track = ConvertToSciFiTrack(_helical_track_fitter, &_geometry_helper);
 
341
 
 
342
      evt.add_scifitrack(track);
 
343
    }
 
344
    for (size_t track_i = 0; track_i < number_straight_tracks; track_i++) {
 
345
      SciFiStraightPRTrack* straight = evt.straightprtracks().at(track_i);
 
346
      Kalman::Track data_track = BuildTrack(straight, &_geometry_helper);
 
347
 
 
348
      Kalman::State seed = ComputeSeed(straight, &_geometry_helper, _seed_value);
 
349
 
 
350
      _straight_track_fitter->SetData(data_track);
 
351
      _straight_track_fitter->SetSeed(seed);
 
352
 
 
353
      _straight_track_fitter->Filter(false);
 
354
      _straight_track_fitter->Smooth(false);
 
355
 
 
356
      SciFiTrack* track = ConvertToSciFiTrack(_straight_track_fitter, &_geometry_helper);
 
357
 
 
358
      evt.add_scifitrack(track);
 
359
    }
 
360
  } else { // PatRec Not Switched on. => Assume EVERYTHING is a single track.
 
361
    throw Exception(Exception::nonRecoverable,
 
362
        "No pattern recognition, required to seed Kalman",
 
363
        "MapCppTrackerRecon::track_fit()" );
 
364
 
 
365
    /*
 
366
    SciFiClusterPArray event_clusters = evt.clusters();
 
367
    for (int tracker = 0; tracker < 2; ++tracker) {
 
368
 
 
369
      SciFiClusterPArray clusters = find_clusters(event_clusters, tracker);
 
370
      if (clusters.size() != 15) break;
 
371
 
 
372
      Kalman::Track data_track = BuildTrack(clusters, &_geometry_helper);
 
373
      Kalman::State seed(5);
 
374
 
 
375
      bool seed_found = false;
 
376
      for (SciFiClusterPArray::iterator it_c = clusters.begin(); it_c != clusters.end(); ++it_c) {
 
377
        if (((*it_c)->get_station() == 5) && ((*it_c)->get_plane() == 2)) {
 
378
 
 
379
          ThreeVector pos = (*it_c)->get_true_position();
 
380
          ThreeVector mom = (*it_c)->get_true_momentum();
 
381
 
 
382
          TMatrixD state(5, 1);
 
383
          state(0, 0) = pos.x();
 
384
          state(1, 0) = mom.x();
 
385
          state(2, 0) = pos.y();
 
386
          state(3, 0) = mom.y();
 
387
          state(4, 0) = 1.0 / mom.z();
 
388
 
 
389
          seed.SetVector(state);
 
390
          seed.SetPosition(pos.z());
 
391
          seed_found = true;
 
392
          break;
 
393
        }
 
394
      }
 
395
      if (!seed_found) continue;
 
396
 
 
397
      TMatrixD seed_cov(5, 5);
 
398
      seed_cov(0, 0) = _seed_value;
 
399
      seed_cov(1, 1) = _seed_value;
 
400
      seed_cov(2, 2) = _seed_value;
 
401
      seed_cov(3, 3) = _seed_value;
 
402
      seed_cov(4, 4) = _seed_value;
 
403
      seed.SetCovariance(seed_cov);
 
404
 
 
405
      _helical_track_fitter->SetData(data_track);
 
406
      _helical_track_fitter->SetSeed(seed);
 
407
 
 
408
      _helical_track_fitter->Filter(false);
 
409
      _helical_track_fitter->Smooth(false);
 
410
 
 
411
      SciFiTrack* track = ConvertToSciFiTrack(_helical_track_fitter, &_geometry_helper);
 
412
 
 
413
      evt.add_scifitrack(track);
 
414
    }
 
415
    */
 
416
  }
 
417
#endif
 
418
}
 
419
 
 
420
void MapCppTrackerRecon::extrapolate_helical_reference(SciFiEvent& event) const {
 
421
  SciFiHelicalPRTrackPArray helicals = event.helicalprtracks();
 
422
  size_t num_tracks = helicals.size();
 
423
 
 
424
  for (size_t i = 0; i < num_tracks; ++i) {
 
425
    SciFiHelicalPRTrack* track = helicals[i];
 
426
    ThreeVector pos;
 
427
    ThreeVector mom;
 
428
 
 
429
    int tracker = track->get_tracker();
 
430
    ThreeVector reference = _geometry_helper.GetReferencePosition(tracker);
 
431
    CLHEP::HepRotation rotation = _geometry_helper.GetReferenceRotation(tracker);
 
432
    rotation.invert();
 
433
 
 
434
    double r  = track->get_R();
 
435
    double pt = - track->get_charge()*CLHEP::c_light*_geometry_helper.GetFieldValue(tracker)*r;
 
436
    double dsdz = - track->get_dsdz();
 
437
    double x0 = track->get_circle_x0();
 
438
    double y0 = track->get_circle_y0();
 
439
    double s = track->get_line_sz_c();
 
440
    double phi = s / r;
 
441
 
 
442
    pos.setX(x0 + r*cos(phi));
 
443
    pos.setY(y0 + r*sin(phi));
 
444
    pos.setZ(0.0);
 
445
    pos *= rotation;
 
446
    pos += reference;
 
447
 
 
448
    mom.setX(-pt*sin(phi));
 
449
    mom.setY(pt*cos(phi));
 
450
    mom.setZ(-pt/dsdz);
 
451
    mom *= rotation;
 
452
 
 
453
    track->set_reference_position(pos);
 
454
    track->set_reference_momentum(mom);
 
455
  }
 
456
}
 
457
 
 
458
void MapCppTrackerRecon::extrapolate_straight_reference(SciFiEvent& event) const {
 
459
  SciFiStraightPRTrackPArray straights = event.straightprtracks();
 
460
  size_t num_tracks = straights.size();
 
461
 
 
462
  for (size_t i = 0; i < num_tracks; ++i) {
 
463
    SciFiStraightPRTrack* track = straights[i];
 
464
    ThreeVector pos;
 
465
    ThreeVector mom;
 
466
 
 
467
    int tracker = track->get_tracker();
 
468
    ThreeVector reference = _geometry_helper.GetReferencePosition(tracker);
 
469
    CLHEP::HepRotation rotation = _geometry_helper.GetReferenceRotation(tracker);
 
470
    rotation.invert();
 
471
 
 
472
    pos.setX(track->get_x0());
 
473
    pos.setY(track->get_y0());
 
474
    pos.setZ(0.0);
 
475
    pos *= rotation;
 
476
    pos += reference;
 
477
 
 
478
    double mx = track->get_mx();
 
479
    double my = track->get_my();
 
480
    mom.setX(mx*200.0);
 
481
    mom.setY(my*200.0);
 
482
    mom *= rotation;
 
483
    mom.setZ(sqrt(40000.0 - mom[0]*mom[0] + mom[1]*mom[1]));
 
484
 
 
485
    track->set_reference_position(pos);
 
486
    track->set_reference_momentum(mom);
 
487
  }
 
488
}
 
489
 
 
490
 
 
491
void MapCppTrackerRecon::print_event_info(SciFiEvent &event) const {
 
492
  std::cerr << "\n ######  SciFi Event Details  ######\n\n";
 
493
 
 
494
 
 
495
  std::cerr << "Num Digits      : " << event.digits().size() << "\n"
 
496
            << "Num Clusters    : " << event.clusters().size() << "\n"
 
497
            << "Num Spacepoints : " << event.spacepoints().size() << "\n"
 
498
            << "Num Straights   : " << event.straightprtracks().size() << "\n"
 
499
            << "Num Helicals    : " << event.helicalprtracks().size() << "\n\n";
 
500
 
 
501
  SciFiTrackPArray tracks = event.scifitracks();
 
502
  for (unsigned int i = 0; i < tracks.size(); ++i) {
 
503
    SciFiTrackPointPArray trackpoints = tracks[i]->scifitrackpoints();
 
504
    for (unsigned int j = 0; j < trackpoints.size(); ++j) {
 
505
      std::cerr << "Track Point = " << trackpoints[j]->pos()[0] << ", "
 
506
                                    << trackpoints[j]->pos()[1] << ", "
 
507
                                    << trackpoints[j]->pos()[2] << " | "
 
508
                                    << trackpoints[j]->mom()[0] << ", "
 
509
                                    << trackpoints[j]->mom()[1] << ", "
 
510
                                    << trackpoints[j]->mom()[2] << " | "
 
511
//                                  << "DATA = " << (cluster ? cluster->get_alpha() : 0) << '\n
 
512
                                    << "(" << trackpoints[j]->pull() << ", "
 
513
                                    << trackpoints[j]->residual()
 
514
                                    << ", " << trackpoints[j]->smoothed_residual() << ")\n";
 
515
    }
 
516
    std::cerr << "Tracker : " << tracks[i]->tracker()
 
517
              << ", Chi2 = " << tracks[i]->chi2()
 
518
              << ", NDF = " << tracks[i]->ndf()
 
519
              << ", P_Value = " << tracks[i]->P_value()
 
520
              << ", Charge = " << tracks[i]->charge() << '\n';
 
521
 
 
522
    ThreeVector seed_pos = tracks[i]->GetSeedPosition();
 
523
    ThreeVector seed_mom = tracks[i]->GetSeedMomentum();
 
524
    std::vector<double> seed_cov = tracks[i]->GetSeedCovariance();
 
525
 
 
526
    std::cerr << "Seed State = (" << seed_pos[0] << ", " << seed_pos[1] << ", " << seed_pos[2]
 
527
              << ")  (" << seed_mom[0] << ", " << seed_mom[1] << ", " << seed_mom[2] << ")\n";
 
528
 
 
529
    std::cerr << "Seed Covariance = ";
 
530
    for (unsigned int j = 0; j < seed_cov.size(); ++j) {
 
531
      std::cerr << seed_cov[j] << ", ";
 
532
    }
 
533
    std::cerr << '\n';
 
534
  }
 
535
  std::cerr << std::endl;
 
536
 
 
537
 
 
538
 
 
539
 
 
540
  /*
 
541
  std::cerr << "Helix Pz Recon:\n";
 
542
  SciFiHelicalPRTrackPArray helices = event.helicalprtracks();
 
543
  for (unsigned int i = 0; i < helices.size(); ++i)
 
544
  {
 
545
    double radius = helices[i]->get_R();
 
546
    double patrec_pt = 1.19*radius;
 
547
    double patrec_pz = patrec_pt / helices[i]->get_dsdz();
 
548
    std::cerr << "Pz = " << patrec_pz << '\n';
 
549
 
 
550
    DoubleArray phis = helices[i]->get_phi();
 
551
    double start_phi = phis[0];
 
552
    double end_phi = phis[ phis.size()-1 ];
 
553
 
 
554
    double delta_phi = end_phi - start_phi;
 
555
 
 
556
    double Z = (1100.0 * 2.0 * CLHEP::pi) / delta_phi;
 
557
    double pz = (1.19 * Z) / (2.0 * CLHEP::pi);
 
558
 
 
559
    std::cerr << "Pz = " << pz << "\n\n";
 
560
 
 
561
    for (unsigned int j = 0; j < phis.size(); ++j)
 
562
    {
 
563
      std::cerr << phis[j] << '\n';
 
564
    }
 
565
  }
 
566
  std::cerr << std::endl;
 
567
  */
 
568
 
 
569
  /*
 
570
  Squeak::mout(Squeak::info) << event.digits().size() << " "
 
571
                              << event.clusters().size() << " "
 
572
                              << event.spacepoints().size() << "; "
 
573
                              << event.straightprtracks().size() << " "
 
574
                              << event.helicalprtracks().size() << "; ";
 
575
  for (size_t track_i = 0; track_i < event.scifitracks().size(); track_i++) {
 
576
    Squeak::mout(Squeak::info) << " Chi2: " << event.scifitracks()[track_i]->f_chi2() << "; "
 
577
                 << " Chi2: " << event.scifitracks()[track_i]->s_chi2() << "; "
 
578
                 << " P-Value: " << event.scifitracks()[track_i]->P_value() << "; ";
 
579
  }
 
580
  Squeak::mout(Squeak::info) << std::endl;
 
581
  */
 
582
}
 
583
 
 
584
void print_helical(SciFiHelicalPRTrack* helical) {
 
585
  SciFiSpacePointPArray spacepoints = helical->get_spacepoints_pointers();
 
586
  size_t numb_spacepoints = spacepoints.size();
 
587
 
 
588
  std::cerr << "HELICAL TRACK:\n";
 
589
 
 
590
  double r = helical->get_R();
 
591
  double sz = helical->get_line_sz_c();
 
592
  double charge = helical->get_charge();
 
593
  std::cerr << "Radius = " << r << " Sz = " << sz << " Charge = " << charge
 
594
            << " Phi_0 = " << sz/r << "\n";
 
595
 
 
596
  for (size_t i = 0; i < numb_spacepoints; ++i) {
 
597
    SciFiSpacePoint *spacepoint = spacepoints[i];
 
598
 
 
599
    std::cerr << spacepoint->get_tracker() << ", "
 
600
              << spacepoint->get_station() << "  -  "
 
601
              << spacepoint->get_position().x() << ", "
 
602
              << spacepoint->get_position().y() << ", "
 
603
              << spacepoint->get_position().z() << "\n";
 
604
  }
 
605
 
 
606
  std::cerr << std::endl;
 
607
}
 
608
 
 
609
void print_straight(SciFiStraightPRTrack* straight) {
 
610
  SciFiSpacePointPArray spacepoints = straight->get_spacepoints_pointers();
 
611
  size_t numb_spacepoints = spacepoints.size();
 
612
 
 
613
  std::cerr << "STRAIGHT TRACK:\n";
 
614
  for (size_t i = 0; i < numb_spacepoints; ++i) {
 
615
    SciFiSpacePoint *spacepoint = spacepoints[i];
 
616
 
 
617
    std::cerr << spacepoint->get_tracker() << ", "
 
618
              << spacepoint->get_station() << "  -  "
 
619
              << spacepoint->get_position().x() << ", "
 
620
              << spacepoint->get_position().y() << ", "
 
621
              << spacepoint->get_position().z() << "\n";
 
622
  }
 
623
 
 
624
  std::cerr << std::endl;
 
625
}
 
626
 
 
627
SciFiClusterPArray find_clusters(SciFiClusterPArray cluster_array, int tracker) {
 
628
  SciFiClusterPArray found_clusters;
 
629
 
 
630
  for (SciFiClusterPArray::iterator cl_it = cluster_array.begin();
 
631
                                                           cl_it != cluster_array.end(); ++cl_it) {
 
632
    if ((*cl_it)->get_tracker() == tracker) {
 
633
      found_clusters.push_back((*cl_it));
 
634
    }
 
635
  }
 
636
 
 
637
  return found_clusters;
 
638
}
 
639
 
 
640
SciFiSpacePointPArray find_spacepoints(SciFiSpacePointPArray spacepoint_array, int tracker) {
 
641
  SciFiSpacePointPArray found_spacepoints;
 
642
 
 
643
  for (SciFiSpacePointPArray::iterator sp_it = spacepoint_array.begin();
 
644
                                                        sp_it != spacepoint_array.end(); ++sp_it) {
 
645
    if ((*sp_it)->get_tracker() == tracker) {
 
646
      found_spacepoints.push_back((*sp_it));
 
647
    }
 
648
  }
 
649
 
 
650
  return found_spacepoints;
 
651
}
160
652
} // ~namespace MAUS