← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

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

merging in changes in merge branch

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerVirtualsDigitization/MapCppTrackerVirtualsDigitization.cc
23
23
#include "src/common_cpp/DataStructure/ReconEvent.hh"
24
24
#include "src/common_cpp/JsonCppProcessors/SpillProcessor.hh"
25
25
 
26
 
#define QUANTISE_ALPHA
27
 
 
28
26
namespace MAUS {
29
27
  PyMODINIT_FUNC init_MapCppTrackerVirtualsDigitization(void) {
30
28
    PyWrapMapBase<MapCppTrackerVirtualsDigitization>::PyWrapMapBaseModInit
44
42
    if (!Globals::HasInstance()) {
45
43
      GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
46
44
    }
 
45
    Json::Value* json = Globals::GetConfigurationCards();
 
46
    _smear_value = (*json)["SciFiTestVirtualSmear"].asDouble();
 
47
    std::string test_method = (*json)["SciFiTestVirtualMethod"].asString();
 
48
 
 
49
    if ( test_method == "virtual" ) {
 
50
      _test_method = 2;
 
51
    } else if ( test_method == "straight" ) {
 
52
      _test_method = 0;
 
53
    } else if ( test_method == "helix" ) {
 
54
      _test_method = 1;
 
55
    } else {
 
56
      throw MAUS::Exception(Exception::nonRecoverable, "Unknown Test Method for virtual digitiser",
 
57
                                                 "MapCppTrackerVirtualsDigitization::_birth(...)");
 
58
    }
 
59
 
 
60
    Json::Value& straight_desc = (*json)["SciFiTestVirtualTracksStraight"];
 
61
    _straight_rms_pos = straight_desc["rms_position"].asDouble();
 
62
    _straight_rms_ang = straight_desc["rms_angle"].asDouble();
 
63
 
 
64
    Json::Value& helix_desc = (*json)["SciFiTestVirtualTracksHelix"];
 
65
    _helix_rms_pos = helix_desc["rms_position"].asDouble();
 
66
    _helix_rms_pt = helix_desc["rms_pt"].asDouble();
 
67
    _helix_pz = helix_desc["pz"].asDouble();
 
68
 
 
69
 
47
70
    static MiceModule* mice_modules = Globals::GetMonteCarloMiceModules();
48
71
    std::vector<const MiceModule*> modules =
49
72
                        mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi");
54
77
    _default_npe = 10.0;
55
78
    _assignment_tolerance = 1.0e-6;
56
79
 
57
 
    _make_digits = false;
58
 
    _make_clusters = true;
59
 
    _make_perfect_clusters = false;
60
 
    _make_spacepoints = false;
61
 
    _make_perfect_spacepoints = false;
 
80
    _make_digits = (*json)["SciFiTestVirtualMakeDigits"].asBool();
 
81
    _make_clusters = (*json)["SciFiTestVirtualMakeClusters"].asBool();
 
82
    _make_spacepoints = (*json)["SciFiTestVirtualMakeSpacepoints"].asBool();
62
83
 
63
84
    _spacepoint_plane_num = 0;
64
85
 
65
 
    _helix_measure = new SciFiHelicalMeasurements(&_geometry_helper);
 
86
//    _helix_measure = new SciFiHelicalMeasurements(&_geometry_helper);
66
87
  }
67
88
 
68
89
 
69
90
  void MapCppTrackerVirtualsDigitization::_death() {
70
 
    delete _helix_measure;
 
91
//    delete _helix_measure;
71
92
  }
72
93
 
73
94
  void MapCppTrackerVirtualsDigitization::_process(MAUS::Data* data) const {
86
107
      ReconEventPArray* revts = new ReconEventPArray();
87
108
      spill.SetReconEvents(revts);
88
109
    }
89
 
    // Construct digits from virtual plane hits for each MC event
 
110
 
90
111
    for (size_t event_i = 0; event_i < spill.GetMCEventSize(); event_i++) {
91
112
      MCEvent *mc_evt = spill.GetMCEvents()->at(event_i);
92
113
      SciFiEvent *sf_evt = new SciFiEvent();
93
114
 
94
 
      if (_make_digits) {
95
 
        SciFiDigitPArray digits;
96
 
        construct_digits(mc_evt->GetVirtualHits(), spill.GetSpillNumber(), event_i, digits);
97
 
        sf_evt->set_digits(digits);
98
 
      }
99
 
 
100
 
      if (_make_clusters) {
101
 
        SciFiClusterPArray clusters;
102
 
        construct_clusters(mc_evt->GetVirtualHits(), spill.GetSpillNumber(), event_i, clusters);
103
 
        sf_evt->set_clusters(clusters);
104
 
      }
105
 
 
106
 
      if (_make_spacepoints) {
107
 
        SciFiSpacePointPArray spacepoints;
108
 
        construct_spacepoints(mc_evt->GetVirtualHits(), spill.GetSpillNumber(), event_i,
109
 
                                                                                      spacepoints);
110
 
        sf_evt->set_spacepoints(spacepoints);
111
 
      }
112
 
 
113
 
      if (_make_perfect_spacepoints) {
114
 
        SciFiSpacePointPArray spacepoints;
115
 
        construct_perfect_spacepoints(spill.GetSpillNumber(), event_i, spacepoints);
116
 
        sf_evt->set_spacepoints(spacepoints);
 
115
      switch (_test_method) {
 
116
        case 0 :
 
117
          _process_straight(sf_evt, spill.GetSpillNumber(), event_i);
 
118
          break;
 
119
        case 1 :
 
120
          _process_helix(sf_evt, spill.GetSpillNumber(), event_i);
 
121
          break;
 
122
        case 2 :
 
123
        // Construct digits from virtual plane hits for each MC event
 
124
        if (_make_digits) {
 
125
          SciFiDigitPArray digits;
 
126
          construct_virtual_digits(mc_evt->GetVirtualHits(),
 
127
                                                          spill.GetSpillNumber(), event_i, digits);
 
128
          sf_evt->set_digits(digits);
 
129
        }
 
130
 
 
131
        if (_make_clusters) {
 
132
          SciFiClusterPArray clusters;
 
133
          construct_virtual_clusters(mc_evt->GetVirtualHits(),
 
134
                                                        spill.GetSpillNumber(), event_i, clusters);
 
135
          sf_evt->set_clusters(clusters);
 
136
        }
 
137
 
 
138
        if (_make_spacepoints) {
 
139
          SciFiSpacePointPArray spacepoints;
 
140
          construct_virtual_spacepoints(mc_evt->GetVirtualHits(),
 
141
                                                     spill.GetSpillNumber(), event_i, spacepoints);
 
142
          sf_evt->set_spacepoints(spacepoints);
 
143
        }
 
144
        break;
 
145
        default:
 
146
          break;
117
147
      }
118
148
 
119
149
      // If there is already a Recon event associated with this MC event, add the SciFiEvent to it,
129
159
    }
130
160
  }
131
161
 
132
 
  void MapCppTrackerVirtualsDigitization::construct_digits(VirtualHitArray* hits, int spill_num,
133
 
                                                   int event_num, SciFiDigitPArray& digits) const {
 
162
 
 
163
  void MapCppTrackerVirtualsDigitization::_process_straight(MAUS::SciFiEvent* sf_event,
 
164
                                                              int spill_num, int event_num) const {
 
165
    double x0 = CLHEP::RandGauss::shoot(0.0, _straight_rms_pos);
 
166
    double y0 = CLHEP::RandGauss::shoot(0.0, _straight_rms_pos);
 
167
    double mx = CLHEP::RandGauss::shoot(0.0, _straight_rms_ang);
 
168
    double my = CLHEP::RandGauss::shoot(0.0, _straight_rms_ang);
 
169
 
 
170
    double time = 0.0;
 
171
 
 
172
    SciFiDigitPArray digits;
 
173
    SciFiClusterPArray clusters;
 
174
    SciFiSpacePointPArray spacepoints;
 
175
 
 
176
    int id;
 
177
    double Z;
 
178
    int tracker;
 
179
    int station;
 
180
    int plane;
 
181
    double alpha;
 
182
    int channelNumber;
 
183
 
 
184
    for (SciFiTrackerMap::const_iterator gmIt = _geometry_map.begin();
 
185
                                                           gmIt != _geometry_map.end(); ++gmIt) {
 
186
      const SciFiPlaneMap& planes = gmIt->second.Planes;
 
187
      tracker = gmIt->first;
 
188
      if (_make_clusters) {
 
189
        for (SciFiPlaneMap::const_iterator plIt = planes.begin(); plIt != planes.end(); ++plIt) {
 
190
          id = plIt->first;
 
191
          station = ((abs(id)-1) / 3) + 1;
 
192
          plane = (abs(id) - 1) % 3;
 
193
 
 
194
          SciFiPlaneGeometry geo = plIt->second;
 
195
 
 
196
          Z = geo.Position.z();
 
197
          ThreeVector position(x0 + Z*mx, y0 + Z*my, Z);
 
198
 
 
199
          alpha = _compute_alpha(position, geo);
 
200
          channelNumber = floor(0.5 + ((7.0 * geo.CentralFibre) - (2.0 * alpha / geo.Pitch)) / 7.0);
 
201
 
 
202
          if (_make_digits) {
 
203
            SciFiDigit* a_digit = new SciFiDigit(spill_num, event_num,
 
204
                                       tracker, station, plane, channelNumber, _default_npe, time);
 
205
            digits.push_back(a_digit);
 
206
          }
 
207
 
 
208
          SciFiCluster* a_cluster = new SciFiCluster();
 
209
          a_cluster->set_plane(plane);
 
210
          a_cluster->set_station(station);
 
211
          a_cluster->set_tracker(tracker);
 
212
          a_cluster->set_channel(channelNumber);
 
213
          a_cluster->set_spill(spill_num);
 
214
          a_cluster->set_event(event_num);
 
215
          a_cluster->set_npe(10);
 
216
          a_cluster->set_position(position);
 
217
          a_cluster->set_direction(geo.Direction);
 
218
          a_cluster->set_alpha(alpha);
 
219
          a_cluster->set_id(id);
 
220
          a_cluster->set_used(false);
 
221
          clusters.push_back(a_cluster);
 
222
        }
 
223
      }
 
224
 
 
225
      if (_make_spacepoints) {
 
226
        for (SciFiPlaneMap::const_iterator plIt = planes.begin(); plIt != planes.end(); ++plIt) {
 
227
          id = plIt->first;
 
228
          station = ((abs(id)-1) / 3) + 1;
 
229
          plane = (abs(id) - 1) % 3;
 
230
          if ( plane != 1 ) continue;
 
231
 
 
232
          SciFiPlaneGeometry geo = plIt->second;
 
233
          Z = geo.Position.z();
 
234
          ThreeVector position(x0 + Z*mx, y0 + Z*my, Z);
 
235
 
 
236
          double smear_x = CLHEP::RandGauss::shoot(0.0, _smear_value);
 
237
          double smear_y = CLHEP::RandGauss::shoot(0.0, _smear_value);
 
238
          position.SetX(position.x() + smear_x);
 
239
          position.SetY(position.y() + smear_y);
 
240
 
 
241
          SciFiSpacePoint* spoint = new SciFiSpacePoint();
 
242
          spoint->set_tracker(tracker);
 
243
          spoint->set_station(station);
 
244
          spoint->set_npe(10);
 
245
          spoint->set_position(position);
 
246
 
 
247
          if (_make_clusters) {
 
248
            for ( int plane = 0; plane < 3; ++plane ) {
 
249
              for ( SciFiClusterPArray::iterator clIt = clusters.begin();
 
250
                                                                 clIt != clusters.end(); ++clIt ) {
 
251
                if ( (*clIt)->get_tracker() == tracker &&
 
252
                     (*clIt)->get_station() == station &&
 
253
                     (*clIt)->get_plane() == plane ) {
 
254
                  spoint->add_channel((*clIt));
 
255
                }
 
256
              }
 
257
            }
 
258
          }
 
259
          spacepoints.push_back(spoint);
 
260
        }
 
261
      }
 
262
    }
 
263
    sf_event->set_digits(digits);
 
264
    sf_event->set_clusters(clusters);
 
265
    sf_event->set_spacepoints(spacepoints);
 
266
  }
 
267
 
 
268
 
 
269
  void MapCppTrackerVirtualsDigitization::_process_helix(MAUS::SciFiEvent* event,
 
270
                                                              int spill_num, int event_num) const {
 
271
    // Not Yet Implemented
 
272
  }
 
273
 
 
274
 
 
275
  void MapCppTrackerVirtualsDigitization::construct_virtual_digits(VirtualHitArray* hits,
 
276
                                    int spill_num, int event_num, SciFiDigitPArray& digits) const {
134
277
    for (unsigned int hit_i = 0; hit_i < hits->size(); hit_i++) {
135
278
      bool found = false;
136
279
      VirtualHit* a_hit    = &hits->at(hit_i);
155
298
            int id = plIt->first;
156
299
            station = ((abs(id)-1) / 3) + 1;
157
300
            plane = (abs(id) - 1) % 3;
158
 
 
159
 
            alpha = _compute_alpha(a_hit, geo, tracker);
 
301
            position.setZ(geo.Position.z());
 
302
            alpha = _compute_alpha(position, geo);
160
303
            channelNumber = floor(0.5 + ((7.0 * geo.CentralFibre) +
161
304
                                                                 (2.0 * alpha / geo.Pitch)) / 7.0);
162
305
 
173
316
  }
174
317
 
175
318
 
176
 
  void MapCppTrackerVirtualsDigitization::construct_clusters(VirtualHitArray* hits, int spill_num,
177
 
                                               int event_num, SciFiClusterPArray& clusters) const {
 
319
  void MapCppTrackerVirtualsDigitization::construct_virtual_clusters(VirtualHitArray* hits,
 
320
                                int spill_num, int event_num, SciFiClusterPArray& clusters) const {
178
321
    for (unsigned int hit_i = 0; hit_i < hits->size(); hit_i++) {
179
322
      bool found = false;
180
323
      VirtualHit* a_hit    = &hits->at(hit_i);
191
334
        tracker = gmIt->first;
192
335
        ThreeVector tracker_position = gmIt->second.Position;
193
336
        HepRotation tracker_rotation = gmIt->second.Rotation;
194
 
        HepRotation rotation = tracker_rotation.inverse();
195
337
 
196
338
        for (SciFiPlaneMap::const_iterator plIt = planes.begin(); plIt != planes.end(); ++plIt) {
197
339
          SciFiPlaneGeometry geo = plIt->second;
200
342
            int id = plIt->first;
201
343
            station = ((abs(id)-1) / 3) + 1;
202
344
            plane = (abs(id) - 1) % 3;
203
 
            alpha = _compute_alpha(a_hit, geo, tracker);
204
 
 
205
 
            ThreeVector momentum = a_hit->GetMomentum();
206
 
 
207
 
            position *= rotation;
 
345
 
 
346
            position *= tracker_rotation;
208
347
            position.setZ(geo.Position.z());
209
348
 
210
 
            momentum *= rotation;
 
349
            alpha = _compute_alpha(position, geo);
211
350
 
212
 
            int channelNumber = floor(0.5 + ((7.0 * geo.CentralFibre) +
 
351
            int channelNumber = floor(0.5 + ((7.0 * geo.CentralFibre) -
213
352
                                                                 (2.0 * alpha / geo.Pitch)) / 7.0);
214
353
 
215
 
            // Change alpha to be quantised by the fibres!
216
 
#ifdef QUANTISE_ALPHA
217
 
            double new_alpha = 3.5 * geo.Pitch *
218
 
                                             static_cast<double>(channelNumber - geo.CentralFibre);
219
 
#else
220
 
            double new_alpha = alpha;
221
 
#endif
222
 
 
223
354
            SciFiCluster* a_cluster = new SciFiCluster();
224
355
 
225
356
            a_cluster->set_plane(plane);
231
362
            a_cluster->set_npe(10);
232
363
            a_cluster->set_position(position);
233
364
            a_cluster->set_direction(geo.Direction);
234
 
            a_cluster->set_alpha(new_alpha);
 
365
            a_cluster->set_alpha(alpha);
235
366
            a_cluster->set_id(id);
236
367
            a_cluster->set_used(false);
237
368
 
238
 
//            delete a_cluster;
239
369
            clusters.push_back(a_cluster);
240
370
 
241
 
            TMatrixD vector(5, 1);
242
 
            vector(0, 0) = position.x();
243
 
            vector(1, 0) = momentum.x();
244
 
            vector(2, 0) = position.y();
245
 
            vector(3, 0) = momentum.y();
246
 
            vector(4, 0) = 1.0 / momentum.z();
247
 
 
248
 
            TMatrix covariance(5, 5);
249
 
 
250
 
            Kalman::State temp_state(vector, covariance, position.z());
251
 
            temp_state.SetId(id);
252
 
            Kalman::State temp_measurement = _helix_measure->Measure(temp_state);
253
 
 
254
371
            found = true;
255
372
            break;
256
373
          }
260
377
    }
261
378
  }
262
379
 
263
 
  void MapCppTrackerVirtualsDigitization::construct_spacepoints(VirtualHitArray* hits,
 
380
 
 
381
  void MapCppTrackerVirtualsDigitization::construct_virtual_spacepoints(VirtualHitArray* hits,
264
382
                          int spill_num, int event_num, SciFiSpacePointPArray& spacepoints) const {
265
383
    for (unsigned int hit_i = 0; hit_i < hits->size(); hit_i++) {
266
384
      bool found = false;
267
385
      VirtualHit* a_hit    = &hits->at(hit_i);
268
386
      ThreeVector position = a_hit->GetPosition();
269
 
      ThreeVector momentum = a_hit->GetMomentum();
270
387
 
271
388
      // Pull tracker/station/plane information from geometry
272
389
      int tracker;
278
395
        tracker = gmIt->first;
279
396
        ThreeVector tracker_position = gmIt->second.Position;
280
397
        HepRotation tracker_rotation = gmIt->second.Rotation;
281
 
        HepRotation rotation = tracker_rotation.inverse();
282
398
 
283
399
        for (SciFiPlaneMap::const_iterator plIt = planes.begin(); plIt != planes.end(); ++plIt) {
284
400
          SciFiPlaneGeometry geo = plIt->second;
286
402
          station = ((abs(id)-1) / 3) + 1;
287
403
          plane = (abs(id) - 1) % 3;
288
404
 
289
 
          if (plane != _spacepoint_plane_num) continue;
290
 
 
291
405
          if (fabs(geo.GlobalPosition.z() - position.z()) < _assignment_tolerance) {
292
406
 
293
 
            position *= rotation;
 
407
            position *= tracker_rotation;
294
408
            position.setZ(geo.Position.z());
295
409
 
296
 
            momentum *= rotation;
 
410
            double smear_x = CLHEP::RandGauss::shoot(0.0, _smear_value);
 
411
            double smear_y = CLHEP::RandGauss::shoot(0.0, _smear_value);
 
412
            position.SetX(position.x() + smear_x);
 
413
            position.SetY(position.y() + smear_y);
297
414
 
298
415
            SciFiSpacePoint* spoint = new SciFiSpacePoint();
299
416
 
313
430
    }
314
431
  }
315
432
 
316
 
  void MapCppTrackerVirtualsDigitization::construct_perfect_spacepoints(int spill_num,
317
 
                                         int event_num, SciFiSpacePointPArray& spacepoints) const {
318
 
    // TODO : Maybe later...
319
 
  }
320
 
 
321
 
  double MapCppTrackerVirtualsDigitization::_compute_alpha(VirtualHit* ahit,
322
 
                                                SciFiPlaneGeometry& geometry, int tracker) const {
323
 
    HepRotation rot = _geometry_helper.GetReferenceRotation(tracker);
324
 
    ThreeVector globalPos = ahit->GetPosition();
 
433
 
 
434
  double MapCppTrackerVirtualsDigitization::_compute_alpha(ThreeVector& position,
 
435
                                                SciFiPlaneGeometry& geometry) const {
 
436
 
 
437
    ThreeVector vec = geometry.Direction.Orthogonal().Unit();
 
438
    double alpha = position.Dot(vec);
 
439
    int channelNumber = floor(0.5 + ((7.0 * geometry.CentralFibre) -
 
440
                                                            (2.0 * alpha / geometry.Pitch)) / 7.0);
 
441
 
 
442
    double new_alpha;
 
443
    if (_smear_value < 0.0) {
 
444
      new_alpha = 3.5 * geometry.Pitch * static_cast<double>(geometry.CentralFibre - channelNumber);
 
445
    } else {
 
446
      double smear_amount = CLHEP::RandGauss::shoot(0.0, _smear_value);
 
447
      new_alpha = alpha + smear_amount;
 
448
    }
 
449
 
 
450
    return new_alpha;
 
451
 
 
452
//    HepRotation rot = _geometry_helper.GetReferenceRotation(tracker);
 
453
//    ThreeVector globalPos = ahit->GetPosition();
 
454
//    ThreeVector pos = position - geometry.GlobalPosition;
325
455
//    ThreeVector pos = globalPos - geometry.GlobalPosition;
326
 
    ThreeVector pos = globalPos;
327
 
 
328
 
    pos *= rot;
329
 
    ThreeVector vec = geometry.Direction.Orthogonal().Unit();
330
 
 
331
 
    double alpha = pos.Dot(vec);
332
 
    return alpha;
 
456
//    ThreeVector pos = globalPos;
 
457
//
 
458
//    pos *= rot;
 
459
//    ThreeVector vec = geometry.Direction.Orthogonal().Unit();
 
460
//
 
461
//    double alpha = pos.Dot(vec);
 
462
//    double alpha = position.Dot(vec);
 
463
//    return alpha;
333
464
    // The old method alpha = fibre num - central fibre
334
465
//    int fibreNumber = (2.0 * dist / geometry.Pitch) + (7.0*geometry.CentralFibre);
335
466
//    int alpha = floor(fibreNumber / 7.0);