~maus-scifi/maus/tracker_devel

« back to all changes in this revision

Viewing changes to tests/cpp_unit/Recon/Kalman/KalmanSeedTest.cc

  • Committer: cheid001
  • Date: 2015-04-01 11:46:50 UTC
  • mfrom: (1144.1.2 tracker_devel)
  • Revision ID: cheid001-20150401114650-w005e6s21mcnstve
merging updates

Show diffs side-by-side

added added

removed removed

Lines of Context:
49
49
    int station_3 = 3;
50
50
    int station_4 = 4;
51
51
    int station_5 = 5;
52
 
    int tracker = 1;
 
52
    int tracker = 0;
 
53
 
 
54
    charge = -1;
 
55
    field = 0.004;
 
56
 
 
57
    z_positions[0] = 0.0;
 
58
    z_positions[1] = 350.0;
 
59
    z_positions[2] = 650.0;
 
60
    z_positions[3] = 900.0;
 
61
    z_positions[4] = 1100.0;
 
62
 
 
63
    ThreeVector pos;
 
64
    double x, y, z;
 
65
 
 
66
    // SET UP STRSIGHT TRACK
53
67
 
54
68
    SciFiCluster * c0 = new SciFiCluster();
55
69
    c0->set_id(id_0);
71
85
    c12->set_id(id_12);
72
86
    SciFiCluster * c13 = new SciFiCluster();
73
87
    c13->set_id(id_13);
74
 
    _clusters.push_back(c0);
75
 
    _clusters.push_back(c1);
76
 
    _clusters.push_back(c3);
77
 
    _clusters.push_back(c4);
78
 
    _clusters.push_back(c6);
79
 
    _clusters.push_back(c7);
80
 
    _clusters.push_back(c9);
81
 
    _clusters.push_back(c10);
82
 
    _clusters.push_back(c12);
83
 
    _clusters.push_back(c13);
84
 
 
85
 
    x = 7.0;
86
 
    y = 8.0;
87
 
    z = 9.0;
88
 
    ThreeVector pos(x, y, z);
89
 
 
 
88
    _straight_clusters.push_back(c0);
 
89
    _straight_clusters.push_back(c1);
 
90
    _straight_clusters.push_back(c3);
 
91
    _straight_clusters.push_back(c4);
 
92
    _straight_clusters.push_back(c6);
 
93
    _straight_clusters.push_back(c7);
 
94
    _straight_clusters.push_back(c9);
 
95
    _straight_clusters.push_back(c10);
 
96
    _straight_clusters.push_back(c12);
 
97
    _straight_clusters.push_back(c13);
 
98
 
 
99
    // Straight tracks
 
100
    x_0 = 10.0;
 
101
    y_0 = 0.0;
 
102
    mx = 20.0 / 1100.0;
 
103
    my = 0.0;
 
104
 
 
105
    x = x_0;
 
106
    y = 0.0;
 
107
    z = z_positions[0];
 
108
    pos.setX( x );
 
109
    pos.setY( y );
 
110
    pos.setZ( z );
90
111
    SciFiSpacePoint *sp_1 = new SciFiSpacePoint();
91
112
    sp_1->set_tracker(tracker);
92
113
    sp_1->add_channel(c0);
94
115
    sp_1->set_position(pos);
95
116
    sp_1->set_station(station_1);
96
117
 
 
118
    x = x_0 + mx * z_positions[1];
 
119
    y = 0.0;
 
120
    z = z_positions[1];
 
121
    pos.setX( x );
 
122
    pos.setY( y );
 
123
    pos.setZ( z );
97
124
    SciFiSpacePoint *sp_2 = new SciFiSpacePoint();
98
125
    sp_2->set_tracker(tracker);
99
126
    sp_2->add_channel(c3);
101
128
    sp_2->set_position(pos);
102
129
    sp_2->set_station(station_2);
103
130
 
 
131
    x = x_0 + mx * z_positions[2];
 
132
    y = 0.0;
 
133
    z = z_positions[2];
 
134
    pos.setX( x );
 
135
    pos.setY( y );
 
136
    pos.setZ( z );
104
137
    SciFiSpacePoint *sp_3 = new SciFiSpacePoint();
105
138
    sp_3->set_tracker(tracker);
106
139
    sp_3->add_channel(c6);
108
141
    sp_3->set_position(pos);
109
142
    sp_3->set_station(station_3);
110
143
 
 
144
    x = x_0 + mx * z_positions[3];
 
145
    y = 0.0;
 
146
    z = z_positions[3];
 
147
    pos.setX( x );
 
148
    pos.setY( y );
 
149
    pos.setZ( z );
111
150
    SciFiSpacePoint *sp_4 = new SciFiSpacePoint();
112
151
    sp_4->set_tracker(tracker);
113
152
    sp_4->add_channel(c9);
115
154
    sp_4->set_position(pos);
116
155
    sp_4->set_station(station_4);
117
156
 
 
157
    x = x_0 + mx * z_positions[4];
 
158
    y = 0.0;
 
159
    z = z_positions[4];
 
160
    pos.setX( x );
 
161
    pos.setY( y );
 
162
    pos.setZ( z );
118
163
    SciFiSpacePoint *sp_5 = new SciFiSpacePoint();
119
164
    sp_5->set_tracker(tracker);
120
165
    sp_5->add_channel(c12);
122
167
    sp_5->set_position(pos);
123
168
    sp_5->set_station(station_5);
124
169
 
125
 
    _spacepoints.push_back(sp_1);
126
 
    _spacepoints.push_back(sp_2);
127
 
    _spacepoints.push_back(sp_3);
128
 
    _spacepoints.push_back(sp_4);
129
 
    _spacepoints.push_back(sp_5);
 
170
    _straight_spacepoints.push_back(sp_1);
 
171
    _straight_spacepoints.push_back(sp_2);
 
172
    _straight_spacepoints.push_back(sp_3);
 
173
    _straight_spacepoints.push_back(sp_4);
 
174
    _straight_spacepoints.push_back(sp_5);
 
175
 
 
176
    // SET UP HELICAL TRACKS
 
177
 
 
178
    c0 = new SciFiCluster();
 
179
    c0->set_id(id_0);
 
180
    c1 = new SciFiCluster();
 
181
    c1->set_id(id_1);
 
182
    c3 = new SciFiCluster();
 
183
    c3->set_id(id_3);
 
184
    c4 = new SciFiCluster();
 
185
    c4->set_id(id_4);
 
186
    c6 = new SciFiCluster();
 
187
    c6->set_id(id_6);
 
188
    c7 = new SciFiCluster();
 
189
    c7->set_id(id_7);
 
190
    c9 = new SciFiCluster();
 
191
    c9->set_id(id_9);
 
192
    c10 = new SciFiCluster();
 
193
    c10->set_id(id_10);
 
194
    c12 = new SciFiCluster();
 
195
    c12->set_id(id_12);
 
196
    c13 = new SciFiCluster();
 
197
    c13->set_id(id_13);
 
198
    _helical_clusters.push_back(c0);
 
199
    _helical_clusters.push_back(c1);
 
200
    _helical_clusters.push_back(c3);
 
201
    _helical_clusters.push_back(c4);
 
202
    _helical_clusters.push_back(c6);
 
203
    _helical_clusters.push_back(c7);
 
204
    _helical_clusters.push_back(c9);
 
205
    _helical_clusters.push_back(c10);
 
206
    _helical_clusters.push_back(c12);
 
207
    _helical_clusters.push_back(c13);
 
208
 
 
209
    // Helix Parameters:
 
210
    x_c = 0.0;
 
211
    y_c = 0.0;
 
212
    r = 10.0;
 
213
    s_0 = 0.0;
 
214
    dsdz = 0.005;
 
215
 
 
216
    x = x_c + r*cos( ( s_0 - charge*z_positions[0]*dsdz ) / r );
 
217
    y = y_c + r*sin( ( s_0 - charge*z_positions[0]*dsdz ) / r );
 
218
    z = z_positions[0];
 
219
    pos.setX( x );
 
220
    pos.setY( y );
 
221
    pos.setZ( z );
 
222
    sp_1 = new SciFiSpacePoint();
 
223
    sp_1->set_tracker(tracker);
 
224
    sp_1->add_channel(c0);
 
225
    sp_1->add_channel(c1);
 
226
    sp_1->set_position(pos);
 
227
    sp_1->set_station(station_1);
 
228
 
 
229
    x = x_c + r*cos( ( s_0 - charge*z_positions[1]*dsdz ) / r );
 
230
    y = y_c + r*sin( ( s_0 - charge*z_positions[1]*dsdz ) / r );
 
231
    z = z_positions[1];
 
232
    pos.setX( x );
 
233
    pos.setY( y );
 
234
    pos.setZ( z );
 
235
    sp_2 = new SciFiSpacePoint();
 
236
    sp_2->set_tracker(tracker);
 
237
    sp_2->add_channel(c3);
 
238
    sp_2->add_channel(c4);
 
239
    sp_2->set_position(pos);
 
240
    sp_2->set_station(station_2);
 
241
 
 
242
    x = x_c + r*cos( ( s_0 - charge*z_positions[2]*dsdz ) / r );
 
243
    y = y_c + r*sin( ( s_0 - charge*z_positions[2]*dsdz ) / r );
 
244
    z = z_positions[2];
 
245
    pos.setX( x );
 
246
    pos.setY( y );
 
247
    pos.setZ( z );
 
248
    sp_3 = new SciFiSpacePoint();
 
249
    sp_3->set_tracker(tracker);
 
250
    sp_3->add_channel(c6);
 
251
    sp_3->add_channel(c7);
 
252
    sp_3->set_position(pos);
 
253
    sp_3->set_station(station_3);
 
254
 
 
255
    x = x_c + r*cos( ( s_0 - charge*z_positions[3]*dsdz ) / r );
 
256
    y = y_c + r*sin( ( s_0 - charge*z_positions[3]*dsdz ) / r );
 
257
    z = z_positions[3];
 
258
    pos.setX( x );
 
259
    pos.setY( y );
 
260
    pos.setZ( z );
 
261
    sp_4 = new SciFiSpacePoint();
 
262
    sp_4->set_tracker(tracker);
 
263
    sp_4->add_channel(c9);
 
264
    sp_4->add_channel(c10);
 
265
    sp_4->set_position(pos);
 
266
    sp_4->set_station(station_4);
 
267
 
 
268
    x = x_c + r*cos( ( s_0 - charge*z_positions[4]*dsdz ) / r );
 
269
    y = y_c + r*sin( ( s_0 - charge*z_positions[4]*dsdz ) / r );
 
270
    z = z_positions[4];
 
271
    pos.setX( x );
 
272
    pos.setY( y );
 
273
    pos.setZ( z );
 
274
    sp_5 = new SciFiSpacePoint();
 
275
    sp_5->set_tracker(tracker);
 
276
    sp_5->add_channel(c12);
 
277
    sp_5->add_channel(c13);
 
278
    sp_5->set_position(pos);
 
279
    sp_5->set_station(station_5);
 
280
 
 
281
    _helical_spacepoints.push_back(sp_1);
 
282
    _helical_spacepoints.push_back(sp_2);
 
283
    _helical_spacepoints.push_back(sp_3);
 
284
    _helical_spacepoints.push_back(sp_4);
 
285
    _helical_spacepoints.push_back(sp_5);
130
286
  }
131
287
  virtual void TearDown() {
132
 
    // delete spacepoints ------------------------
133
288
    std::vector<SciFiSpacePoint*>::iterator spoint;
134
 
    for (spoint = _spacepoints.begin(); spoint!= _spacepoints.end(); ++spoint) {
135
 
      delete (*spoint);
136
 
    }
137
 
    _spacepoints.resize(0);
138
 
    // delete clusters ------------------------
139
289
    std::vector<SciFiCluster*>::iterator cluster;
140
 
    for (cluster = _clusters.begin(); cluster!= _clusters.end(); ++cluster) {
141
 
      delete (*cluster);
142
 
    }
143
 
    _clusters.resize(0);
 
290
    // delete spacepoints ------------------------
 
291
    for (spoint = _helical_spacepoints.begin(); spoint!= _helical_spacepoints.end(); ++spoint) {
 
292
      delete (*spoint);
 
293
    }
 
294
    _helical_spacepoints.resize(0);
 
295
    // delete clusters ------------------------
 
296
    for (cluster = _helical_clusters.begin(); cluster!= _helical_clusters.end(); ++cluster) {
 
297
      delete (*cluster);
 
298
    }
 
299
    // delete spacepoints ------------------------
 
300
    _helical_clusters.resize(0);
 
301
    for (spoint = _straight_spacepoints.begin(); spoint!= _straight_spacepoints.end(); ++spoint) {
 
302
      delete (*spoint);
 
303
    }
 
304
    _straight_spacepoints.resize(0);
 
305
    // delete clusters ------------------------
 
306
    for (cluster = _straight_clusters.begin(); cluster!= _straight_clusters.end(); ++cluster) {
 
307
      delete (*cluster);
 
308
    }
 
309
    _straight_clusters.resize(0);
144
310
  }
145
 
  double x;
146
 
  double y;
147
 
  double z;
 
311
 
 
312
  double z_positions[5];
 
313
  double charge;
 
314
  double field;
 
315
  // Straight track example
 
316
  double x_0;
 
317
  double y_0;
 
318
  double mx;
 
319
  double my;
 
320
  std::vector<SciFiCluster*>    _straight_clusters;
 
321
  std::vector<SciFiSpacePoint*> _straight_spacepoints;
 
322
  // Helical track example
 
323
  double dsdz;
 
324
  double x_c;
 
325
  double y_c;
 
326
  double r;
 
327
  double s_0;
 
328
  std::vector<SciFiCluster*>    _helical_clusters;
 
329
  std::vector<SciFiSpacePoint*> _helical_spacepoints;
 
330
 
148
331
  static const double err = 1.e-6;
149
 
  std::vector<SciFiSpacePoint*> _spacepoints;
150
 
  std::vector<SciFiCluster*>    _clusters;
151
332
};
152
333
 
153
334
TEST_F(KalmanSeedTest, test_ordering) {
154
335
  MAUS::KalmanSeed seed;
155
336
  // This will load the clusters stored in the _spacepoints
156
 
  seed.RetrieveClusters(_spacepoints);
 
337
  seed.RetrieveClusters(_helical_spacepoints);
157
338
  std::vector<SciFiCluster*> seed_clusters = seed.GetClusters();
158
339
 
159
340
  // Now, let's check that they are ordered.
160
 
  int temp = -1;
 
341
  int temp_id = -1;
161
342
  for ( size_t i = 0; i < seed_clusters.size(); ++i ) {
162
343
    int id = seed_clusters[i]->get_id();
163
 
    EXPECT_TRUE(temp < id);
164
 
    temp = id;
 
344
    EXPECT_TRUE(temp_id < id);
 
345
    temp_id = id;
165
346
  }
166
347
  // Check order of spacepoints.
167
 
  temp = -1;
168
 
  for ( size_t j = 0; j < _spacepoints.size(); ++j ) {
169
 
    int station_number = _spacepoints[j]->get_station();
 
348
  double temp_z = 100000.0;
 
349
  for ( size_t j = 0; j < _helical_spacepoints.size(); ++j ) {
 
350
    double z_pos = _helical_spacepoints[j]->get_position().z();
170
351
    // Check that they are stored in increasing order.
171
 
    EXPECT_TRUE(temp < station_number);
172
 
    temp = station_number;
 
352
    EXPECT_TRUE(temp_z > z_pos);
 
353
    temp_z = z_pos;
173
354
  }
174
355
  // Check if any clusters were lost.
175
 
  EXPECT_EQ(_clusters.size(), seed_clusters.size());
 
356
  EXPECT_EQ(_helical_clusters.size(), seed_clusters.size());
176
357
}
177
358
 
178
359
TEST_F(KalmanSeedTest, test_straight_state_vector) {
184
365
  // Now set up a Straight Pattern Recognition Track
185
366
  //
186
367
  MAUS::SciFiStraightPRTrack straight_track;
187
 
  double mx = 1.;
188
 
  double my = 2.;
 
368
  // Make a fake covariance matrix.
 
369
  // Should be makde by pattern recognition
 
370
  std::vector<double> seed_covariance( 16, 0.0 );
 
371
  seed_covariance[0] = 1000.0;
 
372
  seed_covariance[5] = 1000.0;
 
373
  seed_covariance[10] = 1000.0;
 
374
  seed_covariance[15] = 1000.0;
 
375
  straight_track.set_covariance( seed_covariance );
 
376
 
 
377
  straight_track.set_x0(x_0);
 
378
  straight_track.set_y0(y_0);
189
379
  straight_track.set_my(my);
190
380
  straight_track.set_mx(mx);
191
381
  straight_track.set_tracker(0);
192
 
  straight_track.set_spacepoints_pointers(_spacepoints);
 
382
  straight_track.set_spacepoints_pointers(_straight_spacepoints);
193
383
  // Set up stuff for posterior use.
194
384
  seed.Build(&straight_track);
195
385
  //
199
389
  //
200
390
  // Check the result is the expected.
201
391
  //
 
392
 
 
393
  double x = x_0 + 1100.0*mx;
 
394
  double y = y_0 + 1100.0*my;
202
395
  EXPECT_EQ(a.GetNrows(), 4);
203
396
  EXPECT_NEAR(a(0, 0), x,  err);
204
397
  EXPECT_NEAR(a(1, 0), mx, err);
207
400
}
208
401
 
209
402
TEST_F(KalmanSeedTest, test_helical_state_vector) {
 
403
  int charge = -1.;
 
404
  double c  = CLHEP::c_light;
210
405
  //
211
406
  // Set up Seed object with spacepoints.
212
407
  //
213
408
  MAUS::KalmanSeed seed;
214
 
  seed.SetField(-0.004);
 
409
  seed.SetField(0.004);
215
410
  //
216
411
  // Now set up a Helical Pattern Recognition Track
217
412
  //
218
413
  MAUS::SciFiHelicalPRTrack helical_track;
219
 
  double r     = 3.; // mm
220
 
  double tan_lambda = 1./200.;
221
 
  double dsdz  = 1./tan_lambda;
222
 
  // Is this Pi really necessary?
223
 
  double PI = acos(-1.);
224
 
  double phi_0 = 0.;
225
 
  int charge = -1.;
226
 
  double field = -0.004;
227
 
  double c  = CLHEP::c_light;
228
 
  double pt = charge*c*fabs(field)*r;
229
 
  double pz = pt*tan_lambda;
 
414
  double pt = fabs(charge*c*field*r);
230
415
 
231
 
  DoubleArray phi;
232
 
  phi.push_back(phi_0);
 
416
  helical_track.set_dsdz(dsdz);
233
417
  helical_track.set_charge(charge);
234
 
  helical_track.set_phi(phi);
 
418
  helical_track.set_R(r);
235
419
  helical_track.set_tracker(0);
236
 
  helical_track.set_R(r);
237
 
  helical_track.set_dsdz(dsdz);
238
 
  helical_track.set_phi0(phi_0);
239
 
  helical_track.set_spacepoints_pointers(_spacepoints);
 
420
  helical_track.set_circle_x0( x_c );
 
421
  helical_track.set_circle_y0( y_c );
 
422
  helical_track.set_line_sz_c( s_0 );
 
423
  helical_track.set_spacepoints_pointers(_helical_spacepoints);
 
424
  // Make a fake covariance matrix.
 
425
  // Should be makde by pattern recognition
 
426
  std::vector<double> seed_covariance( 25, 0.0 );
 
427
  seed_covariance[0] = 1000.0;
 
428
  seed_covariance[6] = 1000.0;
 
429
  seed_covariance[12] = 1000.0;
 
430
  seed_covariance[18] = 1000.0;
 
431
  seed_covariance[24] = 1000.0;
 
432
  helical_track.set_covariance( seed_covariance );
240
433
  // Set up stuff for posterior use.
241
434
  seed.Build(&helical_track);
242
435
  //
246
439
  //
247
440
  // Check the result is the expected.
248
441
  //
 
442
  double phi = ( s_0 - charge*z_positions[4]*dsdz ) / r;
 
443
  double x = x_c + r*cos( phi );
 
444
  double y = y_c + r*sin( phi );
 
445
  double px = pt*cos( phi + TMath::PiOver2() );
 
446
  double py = pt*sin( phi + TMath::PiOver2() );
 
447
  double kappa = ( dsdz*charge ) / pt;
 
448
 
 
449
  // Allow less stringent coonstraints on momenutm. Should be tighter!
249
450
  EXPECT_EQ(a.GetNrows(), 5);
250
451
  EXPECT_NEAR(a(0, 0), x, err);
251
 
  // allow 2 MeV error due to the removal of PR bias.
252
 
  EXPECT_NEAR(a(1, 0), pt*cos(phi_0+PI/2.), 2);
 
452
  EXPECT_NEAR(a(1, 0), px, 1.0);
253
453
  EXPECT_NEAR(a(2, 0), y, err);
254
 
  EXPECT_NEAR(a(3, 0), pt*sin(phi_0+PI/2.), 2);
 
454
  EXPECT_NEAR(a(3, 0), py, 1.0);
 
455
  EXPECT_NEAR(a(4, 0), kappa, 2.0);
255
456
}
256
457
 
257
458
} // ~namespace MAUS