← Back to branch summary

~jan.greis/maus/1811

« back to all changes in this revision

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

Merging start

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/cpp_unit/Recon/Kalman/KalmanSeedTest.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
 
#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 = 1;
53
 
 
54
 
    SciFiCluster * c0 = new SciFiCluster();
55
 
    c0->set_id(id_0);
56
 
    SciFiCluster * c1 = new SciFiCluster();
57
 
    c1->set_id(id_1);
58
 
    SciFiCluster * c3 = new SciFiCluster();
59
 
    c3->set_id(id_3);
60
 
    SciFiCluster * c4 = new SciFiCluster();
61
 
    c4->set_id(id_4);
62
 
    SciFiCluster * c6 = new SciFiCluster();
63
 
    c6->set_id(id_6);
64
 
    SciFiCluster * c7 = new SciFiCluster();
65
 
    c7->set_id(id_7);
66
 
    SciFiCluster * c9 = new SciFiCluster();
67
 
    c9->set_id(id_9);
68
 
    SciFiCluster * c10 = new SciFiCluster();
69
 
    c10->set_id(id_10);
70
 
    SciFiCluster * c12 = new SciFiCluster();
71
 
    c12->set_id(id_12);
72
 
    SciFiCluster * c13 = new SciFiCluster();
73
 
    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
 
 
90
 
    SciFiSpacePoint *sp_1 = new SciFiSpacePoint();
91
 
    sp_1->set_tracker(tracker);
92
 
    sp_1->add_channel(c0);
93
 
    sp_1->add_channel(c1);
94
 
    sp_1->set_position(pos);
95
 
    sp_1->set_station(station_1);
96
 
 
97
 
    SciFiSpacePoint *sp_2 = new SciFiSpacePoint();
98
 
    sp_2->set_tracker(tracker);
99
 
    sp_2->add_channel(c3);
100
 
    sp_2->add_channel(c4);
101
 
    sp_2->set_position(pos);
102
 
    sp_2->set_station(station_2);
103
 
 
104
 
    SciFiSpacePoint *sp_3 = new SciFiSpacePoint();
105
 
    sp_3->set_tracker(tracker);
106
 
    sp_3->add_channel(c6);
107
 
    sp_3->add_channel(c7);
108
 
    sp_3->set_position(pos);
109
 
    sp_3->set_station(station_3);
110
 
 
111
 
    SciFiSpacePoint *sp_4 = new SciFiSpacePoint();
112
 
    sp_4->set_tracker(tracker);
113
 
    sp_4->add_channel(c9);
114
 
    sp_4->add_channel(c10);
115
 
    sp_4->set_position(pos);
116
 
    sp_4->set_station(station_4);
117
 
 
118
 
    SciFiSpacePoint *sp_5 = new SciFiSpacePoint();
119
 
    sp_5->set_tracker(tracker);
120
 
    sp_5->add_channel(c12);
121
 
    sp_5->add_channel(c13);
122
 
    sp_5->set_position(pos);
123
 
    sp_5->set_station(station_5);
124
 
 
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);
130
 
  }
131
 
  virtual void TearDown() {
132
 
    // delete spacepoints ------------------------
133
 
    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
 
    std::vector<SciFiCluster*>::iterator cluster;
140
 
    for (cluster = _clusters.begin(); cluster!= _clusters.end(); ++cluster) {
141
 
      delete (*cluster);
142
 
    }
143
 
    _clusters.resize(0);
144
 
  }
145
 
  double x;
146
 
  double y;
147
 
  double z;
148
 
  static const double err = 1.e-6;
149
 
  std::vector<SciFiSpacePoint*> _spacepoints;
150
 
  std::vector<SciFiCluster*>    _clusters;
151
 
};
152
 
 
153
 
TEST_F(KalmanSeedTest, test_ordering) {
154
 
  MAUS::KalmanSeed seed;
155
 
  // This will load the clusters stored in the _spacepoints
156
 
  seed.RetrieveClusters(_spacepoints);
157
 
  std::vector<SciFiCluster*> seed_clusters = seed.GetClusters();
158
 
 
159
 
  // Now, let's check that they are ordered.
160
 
  int temp = -1;
161
 
  for ( size_t i = 0; i < seed_clusters.size(); ++i ) {
162
 
    int id = seed_clusters[i]->get_id();
163
 
    EXPECT_TRUE(temp < id);
164
 
    temp = id;
165
 
  }
166
 
  // 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();
170
 
    // Check that they are stored in increasing order.
171
 
    EXPECT_TRUE(temp < station_number);
172
 
    temp = station_number;
173
 
  }
174
 
  // Check if any clusters were lost.
175
 
  EXPECT_EQ(_clusters.size(), seed_clusters.size());
176
 
}
177
 
 
178
 
TEST_F(KalmanSeedTest, test_straight_state_vector) {
179
 
  //
180
 
  // Set up Seed object with spacepoints.
181
 
  //
182
 
  MAUS::KalmanSeed seed;
183
 
  //
184
 
  // Now set up a Straight Pattern Recognition Track
185
 
  //
186
 
  MAUS::SciFiStraightPRTrack straight_track;
187
 
  double mx = 1.;
188
 
  double my = 2.;
189
 
  straight_track.set_my(my);
190
 
  straight_track.set_mx(mx);
191
 
  straight_track.set_tracker(0);
192
 
  straight_track.set_spacepoints_pointers(_spacepoints);
193
 
  // Set up stuff for posterior use.
194
 
  seed.Build(&straight_track);
195
 
  //
196
 
  // Build an initial state vector.
197
 
  //
198
 
  TMatrixD a = seed.initial_state_vector();
199
 
  //
200
 
  // Check the result is the expected.
201
 
  //
202
 
  EXPECT_EQ(a.GetNrows(), 4);
203
 
  EXPECT_NEAR(a(0, 0), x,  err);
204
 
  EXPECT_NEAR(a(1, 0), mx, err);
205
 
  EXPECT_NEAR(a(2, 0), y,  err);
206
 
  EXPECT_NEAR(a(3, 0), my, err);
207
 
}
208
 
 
209
 
TEST_F(KalmanSeedTest, test_helical_state_vector) {
210
 
  //
211
 
  // Set up Seed object with spacepoints.
212
 
  //
213
 
  MAUS::KalmanSeed seed;
214
 
  seed.SetField(-0.004);
215
 
  //
216
 
  // Now set up a Helical Pattern Recognition Track
217
 
  //
218
 
  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;
230
 
 
231
 
  DoubleArray phi;
232
 
  phi.push_back(phi_0);
233
 
  helical_track.set_charge(charge);
234
 
  helical_track.set_phi(phi);
235
 
  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);
240
 
  // Set up stuff for posterior use.
241
 
  seed.Build(&helical_track);
242
 
  //
243
 
  // Build an initial state vector.
244
 
  //
245
 
  TMatrixD a = seed.initial_state_vector();
246
 
  //
247
 
  // Check the result is the expected.
248
 
  //
249
 
  EXPECT_EQ(a.GetNrows(), 5);
250
 
  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);
253
 
  EXPECT_NEAR(a(2, 0), y, err);
254
 
  EXPECT_NEAR(a(3, 0), pt*sin(phi_0+PI/2.), 2);
255
 
}
256
 
 
257
 
} // ~namespace MAUS