← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

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

merging in changes in merge branch

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerReconTest/MapCppTrackerReconTest.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/MapCppTrackerReconTest/MapCppTrackerReconTest.hh"
 
19
 
 
20
#include <sstream>
 
21
 
 
22
#include "src/common_cpp/API/PyWrapMapBase.hh"
 
23
#include "src/common_cpp/Recon/Kalman/MAUSTrackWrapper.hh"
 
24
 
 
25
namespace MAUS {
 
26
 
 
27
void print_helical(SciFiHelicalPRTrack* helical);
 
28
void print_straight(SciFiStraightPRTrack* straight);
 
29
 
 
30
SciFiClusterPArray find_clusters(SciFiClusterPArray cluster_array, int tracker);
 
31
SciFiSpacePointPArray find_spacepoints(SciFiSpacePointPArray spacepoint_array, int tracker);
 
32
 
 
33
 
 
34
PyMODINIT_FUNC init_MapCppTrackerReconTest(void) {
 
35
  PyWrapMapBase<MapCppTrackerReconTest>::PyWrapMapBaseModInit
 
36
                                        ("MapCppTrackerReconTest", "", "", "", "");
 
37
}
 
38
 
 
39
MapCppTrackerReconTest::MapCppTrackerReconTest() : MapBase<Data>("MapCppTrackerReconTest"),
 
40
  _up_straight_pr_on(true),
 
41
  _down_straight_pr_on(true),
 
42
  _up_helical_pr_on(true),
 
43
  _down_helical_pr_on(true),
 
44
  _spacepoint_helical_track_fitter(NULL),
 
45
  _spacepoint_straight_track_fitter(NULL) {
 
46
}
 
47
 
 
48
MapCppTrackerReconTest::~MapCppTrackerReconTest() {
 
49
}
 
50
 
 
51
void MapCppTrackerReconTest::_birth(const std::string& argJsonConfigDocument) {
 
52
  if (!Globals::HasInstance()) {
 
53
    GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
 
54
  }
 
55
  Json::Value* json = Globals::GetConfigurationCards();
 
56
  int user_up_helical_pr_on    = (*json)["SciFiPRHelicalTkUSOn"].asInt();
 
57
  int user_down_helical_pr_on  = (*json)["SciFiPRHelicalTkDSOn"].asInt();
 
58
  int user_up_straight_pr_on   = (*json)["SciFiPRStraightTkUSOn"].asInt();
 
59
  int user_down_straight_pr_on = (*json)["SciFiPRStraightTkDSOn"].asInt();
 
60
  _kalman_on          = (*json)["SciFiKalmanOn"].asBool();
 
61
  _patrec_on          = (*json)["SciFiPatRecOn"].asBool();
 
62
  _spacepoints_on     = (*json)["SciFiSpacepointReconOn"].asBool();
 
63
  _size_exception     = (*json)["SciFiClustExcept"].asInt();
 
64
  _min_npe            = (*json)["SciFiNPECut"].asDouble();
 
65
  _use_mcs            = (*json)["SciFiKalman_use_MCS"].asBool();
 
66
  _use_eloss          = (*json)["SciFiKalman_use_Eloss"].asBool();
 
67
  _use_patrec_seed    = (*json)["SciFiSeedPatRec"].asBool();
 
68
  _correct_pz         = (*json)["SciFiKalmanCorrectPz"].asBool();
 
69
 
 
70
  MiceModule* module = Globals::GetReconstructionMiceModules();
 
71
  std::vector<const MiceModule*> modules =
 
72
    module->findModulesByPropertyString("SensitiveDetector", "SciFi");
 
73
  _geometry_helper = SciFiGeometryHelper(modules);
 
74
  _geometry_helper.Build();
 
75
 
 
76
  _cluster_recon = SciFiClusterRec(_size_exception, _min_npe, _geometry_helper.GeometryMap());
 
77
 
 
78
  _spacepoint_recon = SciFiSpacePointRec();
 
79
 
 
80
  double up_field = _geometry_helper.GetFieldValue(0);
 
81
  double down_field = _geometry_helper.GetFieldValue(1);
 
82
 
 
83
  if (user_up_helical_pr_on == 2)
 
84
    _up_helical_pr_on = true;
 
85
  else if (user_up_helical_pr_on == 1)
 
86
    _up_helical_pr_on = false;
 
87
  else if (user_up_helical_pr_on == 0 && fabs(up_field) < 0.00001) // 10 Gaus
 
88
    _up_helical_pr_on = false;
 
89
  else if (user_up_helical_pr_on == 0 && fabs(up_field) >= 0.00001)
 
90
    _up_helical_pr_on = true;
 
91
 
 
92
  if (user_down_helical_pr_on == 2)
 
93
    _down_helical_pr_on = true;
 
94
  else if (user_down_helical_pr_on == 1)
 
95
    _down_helical_pr_on = false;
 
96
  else if (user_down_helical_pr_on == 0 && fabs(down_field) < 0.00001)
 
97
    _down_helical_pr_on = false;
 
98
  else if (user_down_helical_pr_on == 0 && fabs(down_field) >= 0.00001)
 
99
    _down_helical_pr_on = true;
 
100
 
 
101
  if (user_up_straight_pr_on == 2)
 
102
    _up_straight_pr_on = true;
 
103
  else if (user_up_straight_pr_on == 1)
 
104
    _up_straight_pr_on = false;
 
105
  else if (user_up_straight_pr_on == 0 && fabs(up_field) < 0.00001)
 
106
    _up_straight_pr_on = true;
 
107
  else if (user_up_straight_pr_on == 0 && fabs(up_field) >= 0.00001)
 
108
    _up_straight_pr_on = false;
 
109
 
 
110
  if (user_down_straight_pr_on == 2)
 
111
    _down_straight_pr_on = true;
 
112
  else if (user_down_straight_pr_on == 1)
 
113
    _down_straight_pr_on = false;
 
114
  else if (user_down_straight_pr_on == 0 && fabs(down_field) < 0.00001)
 
115
    _down_straight_pr_on = true;
 
116
  else if (user_down_straight_pr_on == 0 && fabs(down_field) >= 0.00001)
 
117
    _down_straight_pr_on = false;
 
118
 
 
119
  _pattern_recognition = PatternRecognition();
 
120
  _pattern_recognition.LoadGlobals();
 
121
  _pattern_recognition.set_up_helical_pr_on(_up_helical_pr_on);
 
122
  _pattern_recognition.set_down_helical_pr_on(_down_helical_pr_on);
 
123
  _pattern_recognition.set_up_straight_pr_on(_up_straight_pr_on);
 
124
  _pattern_recognition.set_down_straight_pr_on(_down_straight_pr_on);
 
125
  _pattern_recognition.set_bz_t1(up_field);
 
126
  _pattern_recognition.set_bz_t2(down_field);
 
127
 
 
128
  if (_use_patrec_seed) {
 
129
    _seed_value = -1.0;
 
130
  } else {
 
131
    _seed_value = (*json)["SciFiSeedCovariance"].asDouble();
 
132
  }
 
133
 
 
134
  HelicalPropagator* spacepoint_helical_prop = new HelicalPropagator(&_geometry_helper);
 
135
  spacepoint_helical_prop->SetCorrectPz(_correct_pz);
 
136
  spacepoint_helical_prop->SetIncludeMCS(_use_mcs);
 
137
  spacepoint_helical_prop->SetSubtractELoss(_use_eloss);
 
138
  _spacepoint_helical_track_fitter = new Kalman::TrackFit(spacepoint_helical_prop,
 
139
      new SciFiSpacepointMeasurements<5>(0.186128));
 
140
//      new SciFiSpacepointMeasurements<5>(0.4));
 
141
 
 
142
  StraightPropagator* spacepoint_straight_prop = new StraightPropagator(&_geometry_helper);
 
143
  spacepoint_straight_prop->SetIncludeMCS(_use_mcs);
 
144
  _spacepoint_straight_track_fitter = new Kalman::TrackFit(spacepoint_straight_prop,
 
145
      new SciFiSpacepointMeasurements<4>(0.186128));
 
146
//      new SciFiSpacepointMeasurements<4>(0.4));
 
147
 
 
148
  _spacepoint_recon_plane = 0;
 
149
}
 
150
 
 
151
void MapCppTrackerReconTest::_death() {
 
152
  if (_spacepoint_helical_track_fitter) {
 
153
    delete _spacepoint_helical_track_fitter;
 
154
    _spacepoint_helical_track_fitter = NULL;
 
155
  }
 
156
  if (_spacepoint_straight_track_fitter) {
 
157
    delete _spacepoint_straight_track_fitter;
 
158
    _spacepoint_straight_track_fitter = NULL;
 
159
  }
 
160
}
 
161
 
 
162
void MapCppTrackerReconTest::_process(Data* data) const {
 
163
  Spill& spill = *(data->GetSpill());
 
164
 
 
165
  /* return if not physics spill */
 
166
  if (spill.GetDaqEventType() != "physics_event")
 
167
    return;
 
168
 
 
169
  if (spill.GetReconEvents()) {
 
170
    for (unsigned int k = 0; k < spill.GetReconEvents()->size(); k++) {
 
171
      SciFiEvent *event = spill.GetReconEvents()->at(k)->GetSciFiEvent();
 
172
      if (!event) {
 
173
        continue;
 
174
      }
 
175
 
 
176
      if ( _spacepoints_on ) {
 
177
        // Build Clusters.
 
178
        if (event->digits().size()) {
 
179
          _cluster_recon.process(*event);
 
180
        }
 
181
        // Build SpacePoints.
 
182
        if (event->clusters().size()) {
 
183
          _spacepoint_recon.process(*event);
 
184
          set_spacepoint_global_output(event->spacepoints());
 
185
        }
 
186
      }
 
187
 
 
188
      // Pattern Recognition.
 
189
      if (_patrec_on && event->spacepoints().size()) {
 
190
        _pattern_recognition.process(*event);
 
191
        extrapolate_helical_reference(*event);
 
192
        extrapolate_straight_reference(*event);
 
193
      }
 
194
      // Kalman Track Fit.
 
195
      if (_kalman_on) {
 
196
        track_fit(*event);
 
197
      }
 
198
//      print_event_info(*event);
 
199
    }
 
200
  } else {
 
201
    std::cout << "No recon events found\n";
 
202
  }
 
203
}
 
204
 
 
205
void MapCppTrackerReconTest::track_fit(SciFiEvent &evt) const {
 
206
  std::cerr << "KALMAN_TEST ACTIVE\n";
 
207
  size_t number_helical_tracks = evt.helicalprtracks().size();
 
208
  size_t number_straight_tracks = evt.straightprtracks().size();
 
209
 
 
210
  for (size_t track_i = 0; track_i < number_helical_tracks; track_i++) {
 
211
    SciFiHelicalPRTrack* helical = evt.helicalprtracks().at(track_i);
 
212
    SciFiSpacePointPArray spacepoints = helical->get_spacepoints_pointers();
 
213
 
 
214
    Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
 
215
                                                                        _spacepoint_recon_plane);
 
216
 
 
217
    Kalman::State seed = ComputeSeed(helical, &_geometry_helper, false, _seed_value);
 
218
 
 
219
    std::cerr << Kalman::print_state(seed, "Helical Seed");
 
220
    seed.GetCovariance().Print();
 
221
 
 
222
    std::cerr << "Helical Track\n";
 
223
    print_helical(helical);
 
224
 
 
225
    _spacepoint_helical_track_fitter->SetSeed(seed);
 
226
    _spacepoint_helical_track_fitter->SetData(data_track);
 
227
 
 
228
    _spacepoint_helical_track_fitter->Filter(false);
 
229
    _spacepoint_helical_track_fitter->Smooth(false);
 
230
 
 
231
    Kalman::Track predicted = _spacepoint_helical_track_fitter->GetPredicted();
 
232
    Kalman::Track filtered = _spacepoint_helical_track_fitter->GetFiltered();
 
233
    Kalman::Track smoothed = _spacepoint_helical_track_fitter->GetSmoothed();
 
234
    std::cerr << "Data Track\n";
 
235
    std::cerr << Kalman::print_track(data_track);
 
236
    std::cerr << "Measured Predicted\n";
 
237
    SciFiSpacepointMeasurements<5>* temp_measurement = new SciFiSpacepointMeasurements<5>(1.0);
 
238
    for (unsigned int k = 0; k < predicted.GetLength(); ++k) {
 
239
      Kalman::State temp_state = temp_measurement->Measure(predicted[k]);
 
240
      std::cerr << Kalman::print_state(temp_state);
 
241
    }
 
242
    std::cerr << "Predicted Track\n";
 
243
    std::cerr << Kalman::print_track(predicted);
 
244
    std::cerr << "Filtered Track\n";
 
245
    std::cerr << Kalman::print_track(filtered);
 
246
    std::cerr << "Smoothed Track\n";
 
247
    std::cerr << Kalman::print_track(smoothed);
 
248
    delete temp_measurement;
 
249
 
 
250
    SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_helical_track_fitter,
 
251
                                                                     &_geometry_helper, helical);
 
252
    std::cerr << "Chi Squared : " << track->chi2() << " With " << track->ndf() << std::endl;
 
253
    evt.add_scifitrack(track);
 
254
  }
 
255
  for (size_t track_i = 0; track_i < number_straight_tracks; track_i++) {
 
256
    SciFiStraightPRTrack* straight = evt.straightprtracks().at(track_i);
 
257
    SciFiSpacePointPArray spacepoints = straight->get_spacepoints_pointers();
 
258
 
 
259
    Kalman::Track data_track = BuildSpacepointTrack(spacepoints, &_geometry_helper,
 
260
                                                                        _spacepoint_recon_plane);
 
261
    Kalman::State seed = ComputeSeed(straight, &_geometry_helper, _seed_value);
 
262
 
 
263
    std::cerr << Kalman::print_state(seed, "Straight Seed");
 
264
    std::cerr << "Straight Track\n";
 
265
    print_straight(straight);
 
266
 
 
267
    _spacepoint_straight_track_fitter->SetSeed(seed);
 
268
    _spacepoint_straight_track_fitter->SetData(data_track);
 
269
 
 
270
    _spacepoint_straight_track_fitter->Filter(false);
 
271
    _spacepoint_straight_track_fitter->Smooth(false);
 
272
 
 
273
    Kalman::Track predicted = _spacepoint_straight_track_fitter->GetPredicted();
 
274
    Kalman::Track filtered = _spacepoint_straight_track_fitter->GetFiltered();
 
275
    Kalman::Track smoothed = _spacepoint_straight_track_fitter->GetSmoothed();
 
276
 
 
277
    std::cerr << "Data Track\n";
 
278
    std::cerr << Kalman::print_track(data_track);
 
279
    std::cerr << "Predicted Track\n";
 
280
    std::cerr << Kalman::print_track(predicted);
 
281
    std::cerr << "Filtered Track\n";
 
282
    std::cerr << Kalman::print_track(filtered);
 
283
    std::cerr << "Smoothed Track\n";
 
284
    std::cerr << Kalman::print_track(smoothed);
 
285
 
 
286
    SciFiTrack* track = ConvertToSciFiTrack(_spacepoint_straight_track_fitter,
 
287
                                                                    &_geometry_helper, straight);
 
288
    std::cerr << "Chi Squared : " << track->chi2() << " With " << track->ndf() << std::endl;
 
289
    evt.add_scifitrack(track);
 
290
  }
 
291
}
 
292
 
 
293
void MapCppTrackerReconTest::extrapolate_helical_reference(SciFiEvent& event) const {
 
294
  SciFiHelicalPRTrackPArray helicals = event.helicalprtracks();
 
295
  size_t num_tracks = helicals.size();
 
296
 
 
297
  for (size_t i = 0; i < num_tracks; ++i) {
 
298
    SciFiHelicalPRTrack* track = helicals[i];
 
299
    ThreeVector pos;
 
300
    ThreeVector mom;
 
301
 
 
302
    int tracker = track->get_tracker();
 
303
    ThreeVector reference = _geometry_helper.GetReferencePosition(tracker);
 
304
    CLHEP::HepRotation rotation = _geometry_helper.GetReferenceRotation(tracker);
 
305
 
 
306
    double r  = track->get_R();
 
307
    double pt = - track->get_charge()*CLHEP::c_light*_geometry_helper.GetFieldValue(tracker)*r;
 
308
    double dsdz = - track->get_dsdz();
 
309
    double x0 = track->get_circle_x0();
 
310
    double y0 = track->get_circle_y0();
 
311
    double s = track->get_line_sz_c();
 
312
    double phi = s / r;
 
313
 
 
314
    pos.setX(x0 + r*cos(phi));
 
315
    pos.setY(y0 + r*sin(phi));
 
316
    pos.setZ(0.0);
 
317
    pos *= rotation;
 
318
    pos += reference;
 
319
 
 
320
    mom.setX(-pt*sin(phi));
 
321
    mom.setY(pt*cos(phi));
 
322
    mom.setZ(-pt/dsdz);
 
323
    mom *= rotation;
 
324
 
 
325
    track->set_reference_position(pos);
 
326
    track->set_reference_momentum(mom);
 
327
  }
 
328
}
 
329
 
 
330
void MapCppTrackerReconTest::extrapolate_straight_reference(SciFiEvent& event) const {
 
331
  SciFiStraightPRTrackPArray straights = event.straightprtracks();
 
332
  size_t num_tracks = straights.size();
 
333
 
 
334
  for (size_t i = 0; i < num_tracks; ++i) {
 
335
    SciFiStraightPRTrack* track = straights[i];
 
336
    ThreeVector pos;
 
337
    ThreeVector mom;
 
338
    double default_mom = _geometry_helper.GetDefaultMomentum();
 
339
 
 
340
    int tracker = track->get_tracker();
 
341
    ThreeVector reference = _geometry_helper.GetReferencePosition(tracker);
 
342
    CLHEP::HepRotation rotation = _geometry_helper.GetReferenceRotation(tracker);
 
343
 
 
344
    pos.setX(track->get_x0());
 
345
    pos.setY(track->get_y0());
 
346
    pos.setZ(0.0);
 
347
    pos *= rotation;
 
348
    pos += reference;
 
349
 
 
350
    double mx = track->get_mx();
 
351
    double my = track->get_my();
 
352
    mom.setX(mx*default_mom);
 
353
    mom.setY(my*default_mom);
 
354
    mom *= rotation;
 
355
    if (tracker == 0) mom *= -1.0;
 
356
//    mom.setZ(sqrt(40000.0 - mom[0]*mom[0] + mom[1]*mom[1]));
 
357
    mom.setZ(default_mom);
 
358
 
 
359
    track->set_reference_position(pos);
 
360
    track->set_reference_momentum(mom);
 
361
  }
 
362
}
 
363
 
 
364
 
 
365
void MapCppTrackerReconTest::print_event_info(SciFiEvent &event) const {
 
366
  std::cerr << "\n ######  SciFi Event Details  ######\n\n";
 
367
 
 
368
 
 
369
  std::cerr << "Num Digits      : " << event.digits().size() << "\n"
 
370
            << "Num Clusters    : " << event.clusters().size() << "\n"
 
371
            << "Num Spacepoints : " << event.spacepoints().size() << "\n"
 
372
            << "Num Straights   : " << event.straightprtracks().size() << "\n"
 
373
            << "Num Helicals    : " << event.helicalprtracks().size() << "\n\n";
 
374
 
 
375
  SciFiTrackPArray tracks = event.scifitracks();
 
376
  for (unsigned int i = 0; i < tracks.size(); ++i) {
 
377
    SciFiTrackPointPArray trackpoints = tracks[i]->scifitrackpoints();
 
378
    for (unsigned int j = 0; j < trackpoints.size(); ++j) {
 
379
      std::cerr << "Track Point = " << trackpoints[j]->pos()[0] << ", "
 
380
                                    << trackpoints[j]->pos()[1] << ", "
 
381
                                    << trackpoints[j]->pos()[2] << " | "
 
382
                                    << trackpoints[j]->mom()[0] << ", "
 
383
                                    << trackpoints[j]->mom()[1] << ", "
 
384
                                    << trackpoints[j]->mom()[2] << " | "
 
385
//                                  << "DATA = " << (cluster ? cluster->get_alpha() : 0) << '\n
 
386
                                    << "(" << trackpoints[j]->pull() << ", "
 
387
                                    << trackpoints[j]->residual()
 
388
                                    << ", " << trackpoints[j]->smoothed_residual() << ")\n";
 
389
    }
 
390
    std::cerr << "Tracker : " << tracks[i]->tracker()
 
391
              << ", Chi2 = " << tracks[i]->chi2()
 
392
              << ", NDF = " << tracks[i]->ndf()
 
393
              << ", P_Value = " << tracks[i]->P_value()
 
394
              << ", Charge = " << tracks[i]->charge() << '\n';
 
395
 
 
396
    ThreeVector seed_pos = tracks[i]->GetSeedPosition();
 
397
    ThreeVector seed_mom = tracks[i]->GetSeedMomentum();
 
398
    std::vector<double> seed_cov = tracks[i]->GetSeedCovariance();
 
399
 
 
400
    std::cerr << "Seed State = (" << seed_pos[0] << ", " << seed_pos[1] << ", " << seed_pos[2]
 
401
              << ")  (" << seed_mom[0] << ", " << seed_mom[1] << ", " << seed_mom[2] << ")\n";
 
402
 
 
403
    std::cerr << "Seed Covariance = ";
 
404
    for (unsigned int j = 0; j < seed_cov.size(); ++j) {
 
405
      std::cerr << seed_cov[j] << ", ";
 
406
    }
 
407
    std::cerr << '\n';
 
408
  }
 
409
  std::cerr << std::endl;
 
410
}
 
411
 
 
412
void print_helical(SciFiHelicalPRTrack* helical) {
 
413
  SciFiSpacePointPArray spacepoints = helical->get_spacepoints_pointers();
 
414
  size_t numb_spacepoints = spacepoints.size();
 
415
 
 
416
  std::cerr << "HELICAL TRACK:\n";
 
417
 
 
418
  double r = helical->get_R();
 
419
  double sz = helical->get_line_sz_c();
 
420
  double charge = helical->get_charge();
 
421
  std::cerr << "Radius = " << r << " Sz = " << sz << " Charge = " << charge
 
422
            << " Phi_0 = " << sz/r << "\n";
 
423
 
 
424
  for (size_t i = 0; i < numb_spacepoints; ++i) {
 
425
    SciFiSpacePoint *spacepoint = spacepoints[i];
 
426
 
 
427
    std::cerr << spacepoint->get_tracker() << ", "
 
428
              << spacepoint->get_station() << "  -  "
 
429
              << spacepoint->get_position().x() << ", "
 
430
              << spacepoint->get_position().y() << ", "
 
431
              << spacepoint->get_position().z() << "\n";
 
432
  }
 
433
 
 
434
  std::cerr << std::endl;
 
435
}
 
436
 
 
437
void print_straight(SciFiStraightPRTrack* straight) {
 
438
  SciFiSpacePointPArray spacepoints = straight->get_spacepoints_pointers();
 
439
  size_t numb_spacepoints = spacepoints.size();
 
440
 
 
441
  std::cerr << "STRAIGHT TRACK:\n";
 
442
  for (size_t i = 0; i < numb_spacepoints; ++i) {
 
443
    SciFiSpacePoint *spacepoint = spacepoints[i];
 
444
 
 
445
    std::cerr << spacepoint->get_tracker() << ", "
 
446
              << spacepoint->get_station() << "  -  "
 
447
              << spacepoint->get_position().x() << ", "
 
448
              << spacepoint->get_position().y() << ", "
 
449
              << spacepoint->get_position().z() << "\n";
 
450
  }
 
451
 
 
452
  std::cerr << std::endl;
 
453
}
 
454
 
 
455
SciFiClusterPArray find_clusters(SciFiClusterPArray cluster_array, int tracker) {
 
456
  SciFiClusterPArray found_clusters;
 
457
 
 
458
  for (SciFiClusterPArray::iterator cl_it = cluster_array.begin();
 
459
                                                           cl_it != cluster_array.end(); ++cl_it) {
 
460
    if ((*cl_it)->get_tracker() == tracker) {
 
461
      found_clusters.push_back((*cl_it));
 
462
    }
 
463
  }
 
464
 
 
465
  return found_clusters;
 
466
}
 
467
 
 
468
SciFiSpacePointPArray find_spacepoints(SciFiSpacePointPArray spacepoint_array, int tracker) {
 
469
  SciFiSpacePointPArray found_spacepoints;
 
470
 
 
471
  for (SciFiSpacePointPArray::iterator sp_it = spacepoint_array.begin();
 
472
                                                        sp_it != spacepoint_array.end(); ++sp_it) {
 
473
    if ((*sp_it)->get_tracker() == tracker) {
 
474
      found_spacepoints.push_back((*sp_it));
 
475
    }
 
476
  }
 
477
 
 
478
  return found_spacepoints;
 
479
}
 
480
 
 
481
void MapCppTrackerReconTest::set_spacepoint_global_output(SciFiSpacePointPArray spoints) const {
 
482
  for (auto sp : spoints) {
 
483
    int tracker_id = sp->get_tracker();
 
484
 
 
485
    ThreeVector reference_position = _geometry_helper.GetReferencePosition(tracker_id);
 
486
    CLHEP::HepRotation reference_rotation = _geometry_helper.GetReferenceRotation(tracker_id);
 
487
 
 
488
    ThreeVector global_position = sp->get_position();
 
489
    global_position *= reference_rotation;
 
490
    global_position += reference_position;
 
491
 
 
492
    sp->set_global_position(global_position);
 
493
  }
 
494
}
 
495
 
 
496
} // ~namespace MAUS