← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

Viewing changes to tests/cpp_unit/Recon/Global/MinuitTrackFitterTest.cc

  • Committer: Durga Rajaram
  • Date: 2014-01-14 07:07:02 UTC
  • mfrom: (659.1.80 relcand)
  • Revision ID: durga@fnal.gov-20140114070702-2l1fuj1w6rraw7xe
Tags: MAUS-v0.7.6
MAUS-v0.7.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/cpp_unit/Recon/Global/MinuitTrackFitterTest.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
/* Author: Peter Lane
 
18
 */
 
19
 
 
20
#include <cmath>
 
21
#include <iostream>
 
22
#include <sstream>
 
23
#include <string>
 
24
 
 
25
#include <cstdlib>
 
26
 
 
27
#include "TLorentzVector.h"
 
28
#include "TMinuit.h"
 
29
#include "TObject.h"
 
30
#include "gtest/gtest.h"
 
31
#include "json/reader.h"
 
32
#include "json/value.h"
 
33
 
 
34
#include "Utils/Exception.hh"
 
35
#include "Utils/JsonWrapper.hh"
 
36
#include "DataStructure/Global/Track.hh"
 
37
#include "DataStructure/Global/TrackPoint.hh"
 
38
#include "src/common_cpp/Optics/OpticsModel.hh"
 
39
#include "src/common_cpp/Optics/CovarianceMatrix.hh"
 
40
#include "src/common_cpp/Optics/PhaseSpaceVector.hh"
 
41
#include "Recon/Global/DataStructureHelper.hh"
 
42
#include "Recon/Global/MinuitTrackFitter.hh"
 
43
#include "Recon/Global/Particle.hh"
 
44
 
 
45
namespace GlobalDS = MAUS::DataStructure::Global;
 
46
namespace Recon = MAUS::recon::global;
 
47
 
 
48
Json::Value SetupConfig(int verbose_level);
 
49
 
 
50
class DummyOpticsModel : public MAUS::OpticsModel {
 
51
 public:
 
52
  explicit DummyOpticsModel(Json::Value const * const configuration)
 
53
      : MAUS::OpticsModel(configuration) { }
 
54
  void Build() { }
 
55
  MAUS::CovarianceMatrix Transport(const MAUS::CovarianceMatrix & covariances,
 
56
                                   const double end_plane) const
 
57
  {return covariances;}
 
58
  MAUS::CovarianceMatrix Transport(const MAUS::CovarianceMatrix & covariances,
 
59
                                   const double start_plane,
 
60
                                   const double end_plane) const
 
61
  {return covariances;}
 
62
  MAUS::PhaseSpaceVector Transport(const MAUS::PhaseSpaceVector & vector,
 
63
                                   const double end_plane) const
 
64
  {return vector;}
 
65
  MAUS::PhaseSpaceVector Transport(const MAUS::PhaseSpaceVector & vector,
 
66
                                   const double start_plane,
 
67
                                   const double end_plane) const
 
68
  {return vector;}
 
69
};
 
70
 
 
71
class MinuitTrackFitterTest : public testing::Test, public TObject {
 
72
 public:
 
73
  MinuitTrackFitterTest() {
 
74
    for (int index = 0; index < 3; ++index) {
 
75
      GlobalDS::TrackPoint * track_point = new GlobalDS::TrackPoint();
 
76
      TLorentzVector four_position(kStartVector[2]+index, kStartVector[4]+index,
 
77
                                    0.0, kStartVector[0]+index);
 
78
      track_point->set_position(four_position);
 
79
      TLorentzVector position_error(1.0, 1.0, 1.0, 1.0);
 
80
      track_point->set_position_error(position_error);
 
81
 
 
82
      TLorentzVector four_momentum(kStartVector[3]+index, kStartVector[5]+index,
 
83
                                    280.7780432, kStartVector[1]+index);
 
84
      track_point->set_momentum(four_momentum);
 
85
      TLorentzVector momentum_error(1.0, 1.0, 1.0, 1.0);
 
86
      track_point->set_momentum_error(momentum_error);
 
87
      detector_events_.push_back(track_point);
 
88
    }
 
89
  }
 
90
 
 
91
  ~MinuitTrackFitterTest() {
 
92
    for (int index = 0; index < 3; ++index) {
 
93
      delete detector_events_[index];
 
94
    }
 
95
  }
 
96
 
 
97
 protected:
 
98
  static const std::string kConfigurationString;
 
99
  static const Json::Value kConfiguration;
 
100
  static const double kMuonMass;
 
101
  static const Double_t kStartVector[6];
 
102
  std::vector<const GlobalDS::TrackPoint *> detector_events_;
 
103
  static const DummyOpticsModel kOpticsModel;
 
104
};
 
105
 
 
106
// *************************************************
 
107
// MinuitTrackFitterTest static const initializations
 
108
// *************************************************
 
109
 
 
110
const std::string MinuitTrackFitterTest::kConfigurationString = "{"
 
111
"\"global_recon_minuit_minimizer\":\"MINIMIZE\",\n"
 
112
"\"global_recon_minuit_max_iterations\":2000,\n"
 
113
"\"global_recon_minuit_max_edm\":1e-9,\n"
 
114
"\"global_recon_minuit_rounds\":7,\n"
 
115
"\"global_recon_minuit_parameters\":[\n"
 
116
  "{\"name\":\"Time\", \"fixed\":false, \"initial_value\":0.,"
 
117
  " \"value_step\":0.1, \"min_value\":-3.0, \"max_value\":3.0},\n"
 
118
  "{\"name\":\"Energy\", \"fixed\":false, \"initial_value\":226.19,\n"
 
119
  " \"value_step\":0.01, \"min_value\":250., \"max_value\":350.},\n"
 
120
  "{\"name\":\"x\", \"fixed\":false, \"initial_value\":0.,\n"
 
121
  " \"value_step\":0.1, \"min_value\":-15., \"max_value\":15.},\n"
 
122
  "{\"name\":\"Px\", \"fixed\":false, \"initial_value\":0.,\n"
 
123
  " \"value_step\":0.1, \"min_value\":-15., \"max_value\":15.},\n"
 
124
  "{\"name\":\"y\", \"fixed\":false, \"initial_value\":0.,"
 
125
  " \"value_step\":0.1, \"min_value\":-15., \"max_value\":15.},\n"
 
126
  "{\"name\":\"Py\", \"fixed\":false, \"initial_value\":0.,"
 
127
  " \"value_step\":0.1, \"min_value\":-15., \"max_value\":15.}"
 
128
"]"
 
129
"}";
 
130
 
 
131
const double MinuitTrackFitterTest::kMuonMass
 
132
  = Recon::Particle::GetInstance().GetMass(GlobalDS::kMuMinus);
 
133
 
 
134
const Json::Value MinuitTrackFitterTest::kConfiguration
 
135
  = JsonWrapper::StringToJson(kConfigurationString);
 
136
 
 
137
const DummyOpticsModel MinuitTrackFitterTest::kOpticsModel
 
138
  = DummyOpticsModel(&kConfiguration);
 
139
 
 
140
const Double_t MinuitTrackFitterTest::kStartVector[6] = {
 
141
  0.0, 300.0, 5.0, 7.0, 9.0, 11.0 };
 
142
 
 
143
// ***********
 
144
// test cases
 
145
// ***********
 
146
 
 
147
TEST_F(MinuitTrackFitterTest, Constructor) {
 
148
  bool success = true;
 
149
  try {
 
150
    Recon::MinuitTrackFitter fitter(&kOpticsModel, 0.);
 
151
  } catch (MAUS::Exception exc) {
 
152
    success = false;
 
153
  }
 
154
  EXPECT_TRUE(success);
 
155
}
 
156
 
 
157
TEST_F(MinuitTrackFitterTest, ScoreTrack) {
 
158
  Double_t score = Recon::MinuitTrackFitter::ScoreTrack(
 
159
                      kStartVector,
 
160
                      kOpticsModel,
 
161
                      kMuonMass,
 
162
                      detector_events_);
 
163
  // chi^2 = 6*0^2 + 6*1^2 + 6*2^2 = 30. (see init. of detector_events_)
 
164
  EXPECT_NEAR(30.0, score, 1e-9);
 
165
}
 
166
 
 
167
TEST_F(MinuitTrackFitterTest, Fit) {
 
168
    Recon::MinuitTrackFitter fitter(&kOpticsModel, 0.);
 
169
 
 
170
    std::vector<const GlobalDS::TrackPoint *> raw_track_points;
 
171
    GlobalDS::Track raw_track;
 
172
    for (int index = 0; index < 3; ++index) {
 
173
      GlobalDS::TrackPoint * track_point = new GlobalDS::TrackPoint();
 
174
      TLorentzVector four_position(kStartVector[2]+index, kStartVector[4]+index,
 
175
                                    0.0, kStartVector[0]+index);
 
176
      track_point->set_position(four_position);
 
177
      TLorentzVector position_error(1.0, 1.0, 1.0, 1.0);
 
178
      track_point->set_position_error(position_error);
 
179
 
 
180
      TLorentzVector four_momentum(kStartVector[3]+index, kStartVector[5]+index,
 
181
                                    280.7780432, kStartVector[1]+index);
 
182
      track_point->set_momentum(four_momentum);
 
183
      TLorentzVector momentum_error(1.0, 1.0, 1.0, 1.0);
 
184
      track_point->set_momentum_error(momentum_error);
 
185
 
 
186
      raw_track_points.push_back(track_point);
 
187
      raw_track.AddTrackPoint(track_point);
 
188
    }
 
189
 
 
190
    GlobalDS::Track fit_track;
 
191
    fitter.Fit(&raw_track, &fit_track, "MinuitTrackFitterTest");
 
192
 
 
193
    // Check fit_track
 
194
    for (int index = 0; index < 3; ++index) {
 
195
      delete raw_track_points[index];
 
196
    }
 
197
 
 
198
    std::vector<const GlobalDS::TrackPoint *> fit_track_points
 
199
      = fit_track.GetTrackPoints();
 
200
    ASSERT_EQ(size_t(4), fit_track_points.size());
 
201
 
 
202
    const GlobalDS::TrackPoint * fit_primary_point = fit_track_points[0];
 
203
    Recon::DataStructureHelper helper = Recon::DataStructureHelper::GetInstance();
 
204
    MAUS::PhaseSpaceVector fit_primary_vector
 
205
      = helper.TrackPoint2PhaseSpaceVector(*fit_primary_point);
 
206
 
 
207
    /* Since the dummy optics model does identity transport, one would expect
 
208
     * the fit to settle on the 2nd detector event given the way the
 
209
     * detector events are defined
 
210
     */
 
211
    for (int index = 0; index < 6; ++index) {
 
212
      EXPECT_NEAR(kStartVector[index]+1., fit_primary_vector[index], 1.e-6);
 
213
    }
 
214
}