← Back to branch summary

~mbogomilov/maus/devel3

« back to all changes in this revision

Viewing changes to tests/cpp_unit/Optics/PolynomialOpticsModelTest.cc

  • Committer: Durga Rajaram
  • Date: 2014-01-14 04:39:44 UTC
  • mfrom: (663.38.24 merge)
  • mto: (698.1.1 release)
  • mto: This revision was merged to the branch mainline in revision 693.
  • Revision ID: durga@fnal.gov-20140114043944-qf0i8nksnwu74cll
candidate 0.7.6 - trunk r1026

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/cpp_unit/Optics/PolynomialOpticsModelTest.cc
35
35
#include "CLHEP/Vector/ThreeVector.h"
36
36
 
37
37
#include "BeamTools/BTTracker.hh"
38
 
#include "Interface/Squeal.hh"
 
38
#include "Utils/Exception.hh"
39
39
#include "src/common_cpp/Globals/GlobalsManager.hh"
40
40
#include "src/common_cpp/Optics/CovarianceMatrix.hh"
41
41
#include "src/common_cpp/Optics/PolynomialOpticsModel.hh"
70
70
    (*config)["particle_decay"] = Json::Value(false);
71
71
    simulation->GetPhysicsList()->Setup();
72
72
 
73
 
    (*config)["simulation_reference_particle"] = JsonWrapper::StringToJson(
74
 
      std::string("{\"position\":{\"x\":0.0,\"y\":0.0,\"z\":1000.0},")+
75
 
      std::string("\"momentum\":{\"x\":0.0,\"y\":0.0,\"z\":200.0},")+
76
 
      std::string("\"particle_id\":-13,\"energy\":226.1939223,\"time\":0.0,")+
77
 
      std::string("\"random_seed\":2}")
78
 
    );
 
73
    std::stringstream reference_particle_string;
 
74
    reference_particle_string
 
75
      << std::setprecision(1)
 
76
      << "{\"position\":{\"x\":0.0,\"y\":0.0,\"z\":" << kPrimaryPlane << "},"
 
77
      << "\"momentum\":{\"x\":0.0,\"y\":0.0,\"z\":200.0},"
 
78
      << "\"particle_id\":-13,\"energy\":226.1939223,\"time\":0.0,"
 
79
      << "\"random_seed\":2}";
 
80
    (*config)["simulation_reference_particle"]
 
81
      = JsonWrapper::StringToJson(reference_particle_string.str());
79
82
 
80
83
    Json::Value ellipse(Json::objectValue);
81
84
    ellipse["Emittance_T"] = Json::Value(10.0);
115
118
    virtual_planes_ = new MAUS::VirtualPlaneManager();
116
119
    simulation->SetVirtualPlanes(virtual_planes_);
117
120
    MAUS::VirtualPlane start_plane = MAUS::VirtualPlane::BuildVirtualPlane(
118
 
        CLHEP::HepRotation(), CLHEP::Hep3Vector(0., 0., 1000.), -1, true,
119
 
        1000., BTTracker::z, MAUS::VirtualPlane::ignore, false);
 
121
        CLHEP::HepRotation(), CLHEP::Hep3Vector(0., 0., kPrimaryPlane),
 
122
        -1, true,
 
123
        kPrimaryPlane, BTTracker::z, MAUS::VirtualPlane::ignore, false);
120
124
    virtual_planes_->AddPlane(new MAUS::VirtualPlane(start_plane), NULL);
121
125
 
122
126
    MAUS::VirtualPlane mid_plane = MAUS::VirtualPlane::BuildVirtualPlane(
123
 
        CLHEP::HepRotation(), CLHEP::Hep3Vector(0., 0., 2000.), -1, true,
124
 
        2000., BTTracker::z, MAUS::VirtualPlane::ignore, false);
 
127
        CLHEP::HepRotation(), CLHEP::Hep3Vector(0., 0., kPrimaryPlane+1000.),
 
128
        -1, true,
 
129
        kPrimaryPlane+1000., BTTracker::z, MAUS::VirtualPlane::ignore, false);
125
130
    virtual_planes_->AddPlane(new MAUS::VirtualPlane(mid_plane), NULL);
126
131
 
127
132
    MAUS::VirtualPlane end_plane = MAUS::VirtualPlane::BuildVirtualPlane(
128
 
        CLHEP::HepRotation(), CLHEP::Hep3Vector(0., 0., 3000.), -1, true,
129
 
        3000., BTTracker::z, MAUS::VirtualPlane::ignore, false);
 
133
        CLHEP::HepRotation(), CLHEP::Hep3Vector(0., 0., kPrimaryPlane+2000.),
 
134
        -1, true,
 
135
        kPrimaryPlane+2000., BTTracker::z, MAUS::VirtualPlane::ignore, false);
130
136
    virtual_planes_->AddPlane(new MAUS::VirtualPlane(end_plane), NULL);
131
137
  }
132
138
 
151
157
    }
152
158
  }
153
159
 
 
160
  static const double kPrimaryPlane;
154
161
  static const double kCovariances[36];
155
162
  static const MAUS::CovarianceMatrix kCovarianceMatrix;
156
163
  MAUS::VirtualPlaneManager* virtual_planes_;
162
169
// *************************************************
163
170
// PolynomialOpticsModelTest static const initializations
164
171
// *************************************************
 
172
const double PolynomialOpticsModelTest::kPrimaryPlane = 10;
165
173
const double PolynomialOpticsModelTest::kCovariances[36] = {
166
174
  0., 1., 2., 3., 4., 5.,
167
175
  1., 2., 3., 4., 5., 6.,
180
188
 
181
189
TEST_F(PolynomialOpticsModelTest, Constructor) {
182
190
  const PolynomialOpticsModel optics_model(
183
 
      *MAUS::Globals::GetConfigurationCards());
184
 
}
 
191
      MAUS::Globals::GetConfigurationCards());
 
192
}
 
193
 
 
194
/*
 
195
TEST_F(PolynomialOpticsModelTest, Build) {
 
196
  // Bad build -- incomplete particle tracks
 
197
  Json::Value * config = MAUS::Globals::GetConfigurationCards();
 
198
  (*config)["reference_physics_processes"] = Json::Value("mean_energy_loss");
 
199
  (*config)["physics_processes"] = Json::Value("mean_energy_loss");
 
200
  (*config)["particle_decay"] = Json::Value(true);
 
201
  (*config)["muon_half_life"] = Json::Value(1.0);  // 1 ns -> ~25 cm
 
202
  std::string config_string = JsonWrapper::JsonToString(*config);
 
203
  MAUS::GlobalsManager::DeleteGlobals();
 
204
  MAUS::GlobalsManager::InitialiseGlobals(config_string);
 
205
  PolynomialOpticsModel bad_optics_model_1(config);
 
206
  bool success = false;
 
207
  try {
 
208
    bad_optics_model_1.Build();
 
209
  } catch (MAUS::Exception exc) {
 
210
    success = true;
 
211
  }
 
212
  EXPECT_TRUE(success);
 
213
}
 
214
*/
185
215
 
186
216
TEST_F(PolynomialOpticsModelTest, Accessors) {
187
217
  PolynomialOpticsModel optics_model(
188
 
      *MAUS::Globals::GetConfigurationCards());
189
 
  double first_plane = optics_model.first_plane();
190
 
  ASSERT_DOUBLE_EQ(1000., first_plane);
191
 
 
192
 
  optics_model.set_first_plane(2000.);
193
 
  first_plane = optics_model.first_plane();
194
 
  ASSERT_DOUBLE_EQ(2000., first_plane);
 
218
      MAUS::Globals::GetConfigurationCards());
 
219
  double primary_plane = optics_model.primary_plane();
 
220
  ASSERT_DOUBLE_EQ(kPrimaryPlane, primary_plane);
 
221
 
 
222
  optics_model.set_primary_plane(kPrimaryPlane+1000.);
 
223
  primary_plane = optics_model.primary_plane();
 
224
  ASSERT_DOUBLE_EQ(kPrimaryPlane+1000., primary_plane);
 
225
}
 
226
 
 
227
TEST_F(PolynomialOpticsModelTest, AvailablePositions) {
 
228
  PolynomialOpticsModel optics_model(
 
229
      MAUS::Globals::GetConfigurationCards());
 
230
 
 
231
  // Bad position query before model is built
 
232
  bool success = false;
 
233
  try {
 
234
    optics_model.GetAvailableMapPositions();
 
235
  } catch (MAUS::Exception exc) {
 
236
    success = true;
 
237
  }
 
238
  EXPECT_TRUE(success);
 
239
 
 
240
  // One position per virtual detectors
 
241
  optics_model.Build();
 
242
  const std::vector<int64_t> positions
 
243
    = optics_model.GetAvailableMapPositions();
 
244
  const size_t position_count = positions.size();
 
245
  // start, mid, and end virtual detectors
 
246
  EXPECT_EQ(position_count, static_cast<size_t>(3));
195
247
}
196
248
 
197
249
TEST_F(PolynomialOpticsModelTest, Transport) {
198
250
  PolynomialOpticsModel optics_model(
199
 
      *MAUS::Globals::GetConfigurationCards());
 
251
      MAUS::Globals::GetConfigurationCards());
 
252
 
200
253
  // The configuration specifies a 2m drift between 1m and 3m.
201
254
  optics_model.Build();
202
255
  // Check transport to start plane
203
256
  MAUS::PhaseSpaceVector input_vector(0., 226., 1., 0., 3., 0.);
204
257
  MAUS::PhaseSpaceVector output_vector = optics_model.Transport(input_vector,
205
 
                                                                1000.);
 
258
                                                                kPrimaryPlane);
206
259
  for (int index = 0; index < 6; ++index) {
207
 
    EXPECT_NEAR(input_vector[index], output_vector[index], 1.0e-4);
 
260
    EXPECT_NEAR(input_vector[index], output_vector[index], 5.0e-4);
208
261
  }
209
262
 
210
263
  MAUS::CovarianceMatrix output_errors
211
 
      = optics_model.Transport(kCovarianceMatrix, 1000.);
 
264
      = optics_model.Transport(kCovarianceMatrix, kPrimaryPlane);
212
265
  for (int row = 1; row <= 6; ++row) {
213
266
    for (int column = 1; column <= 6; ++column) {
214
267
      EXPECT_NEAR(kCovarianceMatrix(row, column), output_errors(row, column),
215
 
                  1.0e-4);
 
268
                  5.0e-4);
216
269
    }
217
270
  }
218
271
 
219
272
  // Check transport to end plane
220
273
  MAUS::PhaseSpaceVector expected_vector(7.5466, 226., 1., 0., 3., 0.);
221
 
  output_vector = optics_model.Transport(input_vector, 3000.);
 
274
  output_vector = optics_model.Transport(input_vector, kPrimaryPlane+2000.);
222
275
  for (int index = 0; index < 6; ++index) {
223
276
    EXPECT_NEAR(expected_vector[index], output_vector[index], 5.0e-4);
224
277
  }
225
278
 
226
279
  // Check transport to mid plane
227
280
  expected_vector = MAUS::PhaseSpaceVector(3.7733, 226., 1., 0., 3., 0.);
228
 
  output_vector = optics_model.Transport(input_vector, 2000.);
 
281
  output_vector = optics_model.Transport(input_vector, kPrimaryPlane+1000.);
229
282
  for (int index = 0; index < 6; ++index) {
230
283
    EXPECT_NEAR(expected_vector[index], output_vector[index], 5.0e-4);
231
284
  }
234
287
  // Inverse() function in PolynomialMap is not implemented)
235
288
  bool transport_failed = false;
236
289
  try {
237
 
    output_vector = optics_model.Transport(input_vector, 2000., 3000.);
238
 
  } catch (Squeal squeal) {
 
290
    output_vector = optics_model.Transport(input_vector,
 
291
                                           kPrimaryPlane+1000.,
 
292
                                           kPrimaryPlane+2000.);
 
293
  } catch (MAUS::Exception exc) {
239
294
    transport_failed = true;
240
295
  }
241
296
  EXPECT_TRUE(transport_failed);
276
331
       iter < algorithms.end();
277
332
       ++iter) {
278
333
    (*config)["PolynomialOpticsModel_algorithm"] = Json::Value(*iter);
279
 
    PolynomialOpticsModel optics_model(*config);
 
334
    PolynomialOpticsModel optics_model(config);
280
335
    bool algorithm_failed = false;
281
336
    try {
282
337
      optics_model.Build();
283
 
    } catch (Squeal squeal) {
 
338
    } catch (MAUS::Exception exc) {
284
339
      algorithm_failed = true;
285
340
      std::cout << "DEBUG PolynomialOpticsModelTest_UnsupportedAlgorithms: "
286
341
                << "Algorithm \"" << *iter << "\" failed." << std::endl;