~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-27 21:34:57 UTC
  • Revision ID: cheid001-20150427213457-cqvn7u73q1715sfm
style updates and test corrections, still doesn't pass all tests

Show diffs side-by-side

added added

removed removed

Lines of Context:
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
 
#include <stdlib.h>
18
 
 
19
 
#include "src/common_cpp/Recon/Kalman/KalmanSeed.hh"
20
 
 
21
 
#include "gtest/gtest.h"
22
 
 
23
 
/*
24
 
TODO: Update tests according to changes in the KalmanSeed class.
25
 
*/
26
 
 
27
 
namespace MAUS {
28
 
 
29
 
class KalmanSeedTest : public ::testing::Test {
30
 
 protected:
31
 
  KalmanSeedTest()  {}
32
 
  virtual ~KalmanSeedTest() {}
33
 
  virtual void SetUp()    {
34
 
    // Clusters need id,
35
 
    // spacepoints need station number.
36
 
    int id_0     = 1;
37
 
    int id_1     = 2;
38
 
    int id_3     = 3;
39
 
    int id_4     = 4;
40
 
    int id_6     = 5;
41
 
    int id_7     = 6;
42
 
    int id_9     = 7;
43
 
    int id_10    = 8;
44
 
    int id_12    = 9;
45
 
    int id_13    = 10;
46
 
 
47
 
    int station_1 = 1;
48
 
    int station_2 = 2;
49
 
    int station_3 = 3;
50
 
    int station_4 = 4;
51
 
    int station_5 = 5;
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
67
 
 
68
 
    SciFiCluster * c0 = new SciFiCluster();
69
 
    c0->set_id(id_0);
70
 
    SciFiCluster * c1 = new SciFiCluster();
71
 
    c1->set_id(id_1);
72
 
    SciFiCluster * c3 = new SciFiCluster();
73
 
    c3->set_id(id_3);
74
 
    SciFiCluster * c4 = new SciFiCluster();
75
 
    c4->set_id(id_4);
76
 
    SciFiCluster * c6 = new SciFiCluster();
77
 
    c6->set_id(id_6);
78
 
    SciFiCluster * c7 = new SciFiCluster();
79
 
    c7->set_id(id_7);
80
 
    SciFiCluster * c9 = new SciFiCluster();
81
 
    c9->set_id(id_9);
82
 
    SciFiCluster * c10 = new SciFiCluster();
83
 
    c10->set_id(id_10);
84
 
    SciFiCluster * c12 = new SciFiCluster();
85
 
    c12->set_id(id_12);
86
 
    SciFiCluster * c13 = new SciFiCluster();
87
 
    c13->set_id(id_13);
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 );
111
 
    SciFiSpacePoint *sp_1 = new SciFiSpacePoint();
112
 
    sp_1->set_tracker(tracker);
113
 
    sp_1->add_channel(c0);
114
 
    sp_1->add_channel(c1);
115
 
    sp_1->set_position(pos);
116
 
    sp_1->set_station(station_1);
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 );
124
 
    SciFiSpacePoint *sp_2 = new SciFiSpacePoint();
125
 
    sp_2->set_tracker(tracker);
126
 
    sp_2->add_channel(c3);
127
 
    sp_2->add_channel(c4);
128
 
    sp_2->set_position(pos);
129
 
    sp_2->set_station(station_2);
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 );
137
 
    SciFiSpacePoint *sp_3 = new SciFiSpacePoint();
138
 
    sp_3->set_tracker(tracker);
139
 
    sp_3->add_channel(c6);
140
 
    sp_3->add_channel(c7);
141
 
    sp_3->set_position(pos);
142
 
    sp_3->set_station(station_3);
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 );
150
 
    SciFiSpacePoint *sp_4 = new SciFiSpacePoint();
151
 
    sp_4->set_tracker(tracker);
152
 
    sp_4->add_channel(c9);
153
 
    sp_4->add_channel(c10);
154
 
    sp_4->set_position(pos);
155
 
    sp_4->set_station(station_4);
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 );
163
 
    SciFiSpacePoint *sp_5 = new SciFiSpacePoint();
164
 
    sp_5->set_tracker(tracker);
165
 
    sp_5->add_channel(c12);
166
 
    sp_5->add_channel(c13);
167
 
    sp_5->set_position(pos);
168
 
    sp_5->set_station(station_5);
169
 
 
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);
286
 
  }
287
 
  virtual void TearDown() {
288
 
    std::vector<SciFiSpacePoint*>::iterator spoint;
289
 
    std::vector<SciFiCluster*>::iterator cluster;
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);
310
 
  }
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
 
 
331
 
  static const double err = 1.e-6;
332
 
};
333
 
 
334
 
TEST_F(KalmanSeedTest, test_ordering) {
335
 
  MAUS::KalmanSeed seed;
336
 
  // This will load the clusters stored in the _spacepoints
337
 
  seed.RetrieveClusters(_helical_spacepoints);
338
 
  std::vector<SciFiCluster*> seed_clusters = seed.GetClusters();
339
 
 
340
 
  // Now, let's check that they are ordered.
341
 
  int temp_id = -1;
342
 
  for ( size_t i = 0; i < seed_clusters.size(); ++i ) {
343
 
    int id = seed_clusters[i]->get_id();
344
 
    EXPECT_TRUE(temp_id < id);
345
 
    temp_id = id;
346
 
  }
347
 
  // Check order of spacepoints.
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();
351
 
    // Check that they are stored in increasing order.
352
 
    EXPECT_TRUE(temp_z > z_pos);
353
 
    temp_z = z_pos;
354
 
  }
355
 
  // Check if any clusters were lost.
356
 
  EXPECT_EQ(_helical_clusters.size(), seed_clusters.size());
357
 
}
358
 
 
359
 
TEST_F(KalmanSeedTest, test_straight_state_vector) {
360
 
  //
361
 
  // Set up Seed object with spacepoints.
362
 
  //
363
 
  MAUS::KalmanSeed seed;
364
 
  //
365
 
  // Now set up a Straight Pattern Recognition Track
366
 
  //
367
 
  MAUS::SciFiStraightPRTrack straight_track;
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);
379
 
  straight_track.set_my(my);
380
 
  straight_track.set_mx(mx);
381
 
  straight_track.set_tracker(0);
382
 
  straight_track.set_spacepoints_pointers(_straight_spacepoints);
383
 
  // Set up stuff for posterior use.
384
 
  seed.Build(&straight_track);
385
 
  //
386
 
  // Build an initial state vector.
387
 
  //
388
 
  TMatrixD a = seed.initial_state_vector();
389
 
  //
390
 
  // Check the result is the expected.
391
 
  //
392
 
 
393
 
  double x = x_0 + 1100.0*mx;
394
 
  double y = y_0 + 1100.0*my;
395
 
  EXPECT_EQ(a.GetNrows(), 4);
396
 
  EXPECT_NEAR(a(0, 0), x,  err);
397
 
  EXPECT_NEAR(a(1, 0), mx, err);
398
 
  EXPECT_NEAR(a(2, 0), y,  err);
399
 
  EXPECT_NEAR(a(3, 0), my, err);
400
 
}
401
 
 
402
 
TEST_F(KalmanSeedTest, test_helical_state_vector) {
403
 
  int charge = -1.;
404
 
  double c  = CLHEP::c_light;
405
 
  //
406
 
  // Set up Seed object with spacepoints.
407
 
  //
408
 
  MAUS::KalmanSeed seed;
409
 
  seed.SetField(0.004);
410
 
  //
411
 
  // Now set up a Helical Pattern Recognition Track
412
 
  //
413
 
  MAUS::SciFiHelicalPRTrack helical_track;
414
 
  double pt = fabs(charge*c*field*r);
415
 
 
416
 
  helical_track.set_dsdz(dsdz);
417
 
  helical_track.set_charge(charge);
418
 
  helical_track.set_R(r);
419
 
  helical_track.set_tracker(0);
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 );
433
 
  // Set up stuff for posterior use.
434
 
  seed.Build(&helical_track);
435
 
  //
436
 
  // Build an initial state vector.
437
 
  //
438
 
  TMatrixD a = seed.initial_state_vector();
439
 
  //
440
 
  // Check the result is the expected.
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!
450
 
  EXPECT_EQ(a.GetNrows(), 5);
451
 
  EXPECT_NEAR(a(0, 0), x, err);
452
 
  EXPECT_NEAR(a(1, 0), px, 1.0);
453
 
  EXPECT_NEAR(a(2, 0), y, err);
454
 
  EXPECT_NEAR(a(3, 0), py, 1.0);
455
 
  EXPECT_NEAR(a(4, 0), kappa, 2.0);
456
 
}
457
 
 
458
 
} // ~namespace MAUS
 
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
#include <stdlib.h>
 
18
 
 
19
#include "src/common_cpp/Recon/Kalman/KalmanSeed.hh"
 
20
 
 
21
#include "gtest/gtest.h"
 
22
 
 
23
/*
 
24
TODO: Update tests according to changes in the KalmanSeed class.
 
25
*/
 
26
 
 
27
namespace MAUS {
 
28
 
 
29
class KalmanSeedTest : public ::testing::Test {
 
30
 protected:
 
31
  KalmanSeedTest()  {}
 
32
  virtual ~KalmanSeedTest() {}
 
33
  virtual void SetUp()    {
 
34
    // Clusters need id,
 
35
    // spacepoints need station number.
 
36
    int id_0     = 1;
 
37
    int id_1     = 2;
 
38
    int id_3     = 3;
 
39
    int id_4     = 4;
 
40
    int id_6     = 5;
 
41
    int id_7     = 6;
 
42
    int id_9     = 7;
 
43
    int id_10    = 8;
 
44
    int id_12    = 9;
 
45
    int id_13    = 10;
 
46
 
 
47
    int station_1 = 1;
 
48
    int station_2 = 2;
 
49
    int station_3 = 3;
 
50
    int station_4 = 4;
 
51
    int station_5 = 5;
 
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
 
67
 
 
68
    SciFiCluster * c0 = new SciFiCluster();
 
69
    c0->set_id(id_0);
 
70
    SciFiCluster * c1 = new SciFiCluster();
 
71
    c1->set_id(id_1);
 
72
    SciFiCluster * c3 = new SciFiCluster();
 
73
    c3->set_id(id_3);
 
74
    SciFiCluster * c4 = new SciFiCluster();
 
75
    c4->set_id(id_4);
 
76
    SciFiCluster * c6 = new SciFiCluster();
 
77
    c6->set_id(id_6);
 
78
    SciFiCluster * c7 = new SciFiCluster();
 
79
    c7->set_id(id_7);
 
80
    SciFiCluster * c9 = new SciFiCluster();
 
81
    c9->set_id(id_9);
 
82
    SciFiCluster * c10 = new SciFiCluster();
 
83
    c10->set_id(id_10);
 
84
    SciFiCluster * c12 = new SciFiCluster();
 
85
    c12->set_id(id_12);
 
86
    SciFiCluster * c13 = new SciFiCluster();
 
87
    c13->set_id(id_13);
 
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);
 
111
    SciFiSpacePoint *sp_1 = new SciFiSpacePoint();
 
112
    sp_1->set_tracker(tracker);
 
113
    sp_1->add_channel(c0);
 
114
    sp_1->add_channel(c1);
 
115
    sp_1->set_position(pos);
 
116
    sp_1->set_station(station_1);
 
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);
 
124
    SciFiSpacePoint *sp_2 = new SciFiSpacePoint();
 
125
    sp_2->set_tracker(tracker);
 
126
    sp_2->add_channel(c3);
 
127
    sp_2->add_channel(c4);
 
128
    sp_2->set_position(pos);
 
129
    sp_2->set_station(station_2);
 
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);
 
137
    SciFiSpacePoint *sp_3 = new SciFiSpacePoint();
 
138
    sp_3->set_tracker(tracker);
 
139
    sp_3->add_channel(c6);
 
140
    sp_3->add_channel(c7);
 
141
    sp_3->set_position(pos);
 
142
    sp_3->set_station(station_3);
 
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);
 
150
    SciFiSpacePoint *sp_4 = new SciFiSpacePoint();
 
151
    sp_4->set_tracker(tracker);
 
152
    sp_4->add_channel(c9);
 
153
    sp_4->add_channel(c10);
 
154
    sp_4->set_position(pos);
 
155
    sp_4->set_station(station_4);
 
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);
 
163
    SciFiSpacePoint *sp_5 = new SciFiSpacePoint();
 
164
    sp_5->set_tracker(tracker);
 
165
    sp_5->add_channel(c12);
 
166
    sp_5->add_channel(c13);
 
167
    sp_5->set_position(pos);
 
168
    sp_5->set_station(station_5);
 
169
 
 
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);
 
286
  }
 
287
  virtual void TearDown() {
 
288
    std::vector<SciFiSpacePoint*>::iterator spoint;
 
289
    std::vector<SciFiCluster*>::iterator cluster;
 
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);
 
310
  }
 
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
 
 
331
  static const double err = 1.e-6;
 
332
};
 
333
 
 
334
TEST_F(KalmanSeedTest, test_ordering) {
 
335
  MAUS::KalmanSeed seed;
 
336
  // This will load the clusters stored in the _spacepoints
 
337
  seed.RetrieveClusters(_helical_spacepoints);
 
338
  std::vector<SciFiCluster*> seed_clusters = seed.GetClusters();
 
339
 
 
340
  // Now, let's check that they are ordered.
 
341
  int temp_id = -1;
 
342
  for ( size_t i = 0; i < seed_clusters.size(); ++i ) {
 
343
    int id = seed_clusters[i]->get_id();
 
344
    EXPECT_TRUE(temp_id < id);
 
345
    temp_id = id;
 
346
  }
 
347
  // Check order of spacepoints.
 
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();
 
351
    // Check that they are stored in increasing order.
 
352
    EXPECT_TRUE(temp_z > z_pos);
 
353
    temp_z = z_pos;
 
354
  }
 
355
  // Check if any clusters were lost.
 
356
  EXPECT_EQ(_helical_clusters.size(), seed_clusters.size());
 
357
}
 
358
 
 
359
TEST_F(KalmanSeedTest, test_straight_state_vector) {
 
360
  //
 
361
  // Set up Seed object with spacepoints.
 
362
  //
 
363
  MAUS::KalmanSeed seed;
 
364
  //
 
365
  // Now set up a Straight Pattern Recognition Track
 
366
  //
 
367
  MAUS::SciFiStraightPRTrack straight_track;
 
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);
 
379
  straight_track.set_my(my);
 
380
  straight_track.set_mx(mx);
 
381
  straight_track.set_tracker(0);
 
382
  straight_track.set_spacepoints_pointers(_straight_spacepoints);
 
383
  // Set up stuff for posterior use.
 
384
  seed.Build(&straight_track);
 
385
  //
 
386
  // Build an initial state vector.
 
387
  //
 
388
  TMatrixD a = seed.initial_state_vector();
 
389
  //
 
390
  // Check the result is the expected.
 
391
  //
 
392
 
 
393
  double x = x_0 + 1100.0*mx;
 
394
  double y = y_0 + 1100.0*my;
 
395
  EXPECT_EQ(a.GetNrows(), 4);
 
396
  EXPECT_NEAR(a(0, 0), x,  err);
 
397
  EXPECT_NEAR(a(1, 0), mx, err);
 
398
  EXPECT_NEAR(a(2, 0), y,  err);
 
399
  EXPECT_NEAR(a(3, 0), my, err);
 
400
}
 
401
 
 
402
TEST_F(KalmanSeedTest, test_helical_state_vector) {
 
403
  int charge = -1.;
 
404
  double c  = CLHEP::c_light;
 
405
  //
 
406
  // Set up Seed object with spacepoints.
 
407
  //
 
408
  MAUS::KalmanSeed seed;
 
409
  seed.SetField(0.004);
 
410
  //
 
411
  // Now set up a Helical Pattern Recognition Track
 
412
  //
 
413
  MAUS::SciFiHelicalPRTrack helical_track;
 
414
  double pt = fabs(charge*c*field*r);
 
415
 
 
416
  helical_track.set_dsdz(dsdz);
 
417
  helical_track.set_charge(charge);
 
418
  helical_track.set_R(r);
 
419
  helical_track.set_tracker(0);
 
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);
 
433
  // Set up stuff for posterior use.
 
434
  seed.Build(&helical_track);
 
435
  //
 
436
  // Build an initial state vector.
 
437
  //
 
438
  TMatrixD a = seed.initial_state_vector();
 
439
  //
 
440
  // Check the result is the expected.
 
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!
 
450
  EXPECT_EQ(a.GetNrows(), 5);
 
451
  EXPECT_NEAR(a(0, 0), x, err);
 
452
  EXPECT_NEAR(a(1, 0), px, 1.0);
 
453
  EXPECT_NEAR(a(2, 0), y, err);
 
454
  EXPECT_NEAR(a(3, 0), py, 1.0);
 
455
  EXPECT_NEAR(a(4, 0), kappa, 2.0);
 
456
}
 
457
 
 
458
} // ~namespace MAUS