← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

Viewing changes to src/map/MapCppTrackerIntegratedTrackFit/MapCppTrackerIntegratedTrackFit.cc

  • Committer: Christopher Hunt
  • Date: 2018-07-04 14:12:12 UTC
  • Revision ID: christopher.hunt08@imperial.ac.uk-20180704141212-orzzl40mi53xc6kh
Added files for field integrated kalman filter. Now to test it...

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerIntegratedTrackFit/MapCppTrackerIntegratedTrackFit.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
 
 
18
#include "src/map/MapCppTrackerIntegratedTrackFit/MapCppTrackerIntegratedTrackFit.hh"
 
19
 
 
20
#include <sstream>
 
21
 
 
22
#include "src/common_cpp/API/PyWrapMapBase.hh"
 
23
#include "src/common_cpp/Recon/Kalman/MAUSTrackWrapper.hh"
 
24
 
 
25
namespace MAUS {
 
26
 
 
27
PyMODINIT_FUNC init_MapCppTrackerIntegratedTrackFit(void) {
 
28
  PyWrapMapBase<MapCppTrackerIntegratedTrackFit>::PyWrapMapBaseModInit
 
29
                                        ("MapCppTrackerIntegratedTrackFit", "", "", "", "");
 
30
}
 
31
 
 
32
 
 
33
MapCppTrackerIntegratedTrackFit::MapCppTrackerIntegratedTrackFit() : MapBase<Data>("MapCppTrackerIntegratedTrackFit"),
 
34
                                           _kalman_on(true),
 
35
                                           _helical_track_fitter(NULL),
 
36
                                           _straight_track_fitter(NULL) {
 
37
}
 
38
 
 
39
 
 
40
MapCppTrackerIntegratedTrackFit::~MapCppTrackerIntegratedTrackFit() {
 
41
}
 
42
 
 
43
 
 
44
void MapCppTrackerIntegratedTrackFit::_birth(const std::string& argJsonConfigDocument) {
 
45
  // Pull out the global settings
 
46
  if (!Globals::HasInstance()) {
 
47
    GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
 
48
  }
 
49
  Json::Value* json = Globals::GetConfigurationCards();
 
50
  _kalman_on          = (*json)["SciFiKalmanOn"].asBool();
 
51
  _use_patrec_seed    = (*json)["SciFiSeedPatRec"].asBool();
 
52
  // Values used to set the track rating:
 
53
  _excellent_num_trackpoints     = (*json)["SciFiExcellentNumTrackpoints"].asInt();
 
54
  _good_num_trackpoints          = (*json)["SciFiGoodNumTrackpoints"].asInt();
 
55
  _poor_num_trackpoints          = (*json)["SciFiPoorNumTrackpoints"].asInt();
 
56
  _excellent_p_value             = (*json)["SciFiExcellentPValue"].asDouble();
 
57
  _good_p_value                  = (*json)["SciFiGoodPValue"].asDouble();
 
58
  _poor_p_value                  = (*json)["SciFiPoorPValue"].asDouble();
 
59
  _excellent_num_spacepoints  = (*json)["SciFiExcellentNumSpacepoints"].asInt();
 
60
  _good_num_spacepoints       = (*json)["SciFiGoodNumSpacepoints"].asInt();
 
61
  _poor_num_spacepoints       = (*json)["SciFiPoorNumSpacepoints"].asInt();
 
62
 
 
63
  SciFiTrackerMap& geo_map = Globals::GetSciFiGeometryHelper()->GeometryMap();
 
64
 
 
65
  if (_use_patrec_seed) {
 
66
    _seed_value = -1.0;
 
67
  } else {
 
68
    _seed_value = (*json)["SciFiSeedCovariance"].asDouble();
 
69
  }
 
70
 
 
71
  // Set up final track fit (Kalman filter)
 
72
  FieldIntegratedPropagator* helical_prop = new FieldIntegratedPropagator(Globals::GetSciFiGeometryHelper());
 
73
  _helical_track_fitter = new Kalman::TrackFit(helical_prop);
 
74
 
 
75
  StraightPropagator* straight_prop = new StraightPropagator(Globals::GetSciFiGeometryHelper());
 
76
  straight_prop->SetIncludeMCS(_use_mcs);
 
77
  _straight_track_fitter = new Kalman::TrackFit(straight_prop);
 
78
 
 
79
  // Each measurement plane has a unique alignment and rotation => they all need their own
 
80
  //  measurement object.
 
81
  for (SciFiTrackerMap::iterator track_it = geo_map.begin();
 
82
                                                           track_it != geo_map.end(); ++track_it) {
 
83
    int tracker_const = (track_it->first == 0 ? -1 : 1);
 
84
    for (SciFiPlaneMap::iterator plane_it = track_it->second.Planes.begin();
 
85
                                           plane_it != track_it->second.Planes.end(); ++plane_it) {
 
86
 
 
87
      int id = plane_it->first * tracker_const;
 
88
      _helical_track_fitter->AddMeasurement(id,
 
89
                                             new MAUS::SciFiHelicalMeasurements(plane_it->second));
 
90
      _straight_track_fitter->AddMeasurement(id,
 
91
                                            new MAUS::SciFiStraightMeasurements(plane_it->second));
 
92
    }
 
93
  }
 
94
}
 
95
 
 
96
 
 
97
void MapCppTrackerIntegratedTrackFit::_death() {
 
98
  if (_helical_track_fitter) {
 
99
    delete _helical_track_fitter;
 
100
    _helical_track_fitter = NULL;
 
101
  }
 
102
  if (_straight_track_fitter) {
 
103
    delete _straight_track_fitter;
 
104
    _straight_track_fitter = NULL;
 
105
  }
 
106
}
 
107
 
 
108
 
 
109
void MapCppTrackerIntegratedTrackFit::_process(Data* data) const {
 
110
  Spill& spill = *(data->GetSpill());
 
111
 
 
112
  /* return if not physics spill */
 
113
  if (spill.GetDaqEventType() != "physics_event")
 
114
    return;
 
115
 
 
116
  if (spill.GetReconEvents()) {
 
117
    for (unsigned int k = 0; k < spill.GetReconEvents()->size(); k++) {
 
118
      SciFiEvent *event = spill.GetReconEvents()->at(k)->GetSciFiEvent();
 
119
      if (!event) {
 
120
        continue;
 
121
      }
 
122
      _set_field_values(event);
 
123
 
 
124
      // Kalman Track Fit.
 
125
      if (_kalman_on) {
 
126
        event->clear_scifitracks();
 
127
        SciFiSeedPArray seeds = event->scifiseeds();
 
128
        for (SciFiSeedPArray::iterator it = seeds.begin(); it != seeds.end(); ++it) {
 
129
          if ((*it)->getTracker() == 1) continue;
 
130
          try {
 
131
            switch ((*it)->getAlgorithm()) {
 
132
              case 0 :
 
133
                event->add_scifitrack(track_fit_straight((*it)));
 
134
                break;
 
135
              case 1 :
 
136
                event->add_scifitrack(track_fit_helix((*it)));
 
137
                break;
 
138
              default:
 
139
                break;
 
140
            }
 
141
          } catch (Exceptions::Exception& e) {
 
142
            event->add_scifitrack(make_empty((*it)));
 
143
            std::cerr << "Track Fit Failed: " << e.what();
 
144
          }
 
145
        }
 
146
      }
 
147
    }
 
148
  } else {
 
149
    std::cout << "No recon events found\n";
 
150
  }
 
151
}
 
152
 
 
153
SciFiTrack* MapCppTrackerIntegratedTrackFit::make_empty(SciFiSeed* seed) const {
 
154
  SciFiTrack* empty_track = new SciFiTrack();
 
155
  empty_track->set_scifi_seed(seed);
 
156
  empty_track->set_tracker(seed->getTracker());
 
157
  empty_track->set_chi2(-1.0);
 
158
  empty_track->set_P_value(-1.0);
 
159
 
 
160
  return empty_track;
 
161
}
 
162
 
 
163
void MapCppTrackerIntegratedTrackFit::_set_field_values(SciFiEvent* event) const {
 
164
  event->set_mean_field_up(Globals::GetSciFiGeometryHelper()->GetFieldValue(0));
 
165
  event->set_mean_field_down(Globals::GetSciFiGeometryHelper()->GetFieldValue(1));
 
166
 
 
167
  event->set_variance_field_up(Globals::GetSciFiGeometryHelper()->GetFieldVariance(0));
 
168
  event->set_variance_field_down(Globals::GetSciFiGeometryHelper()->GetFieldVariance(1));
 
169
 
 
170
  event->set_range_field_up(Globals::GetSciFiGeometryHelper()->GetFieldRange(0));
 
171
  event->set_range_field_down(Globals::GetSciFiGeometryHelper()->GetFieldRange(1));
 
172
}
 
173
 
 
174
 
 
175
SciFiTrack* MapCppTrackerIntegratedTrackFit::track_fit_helix(SciFiSeed* seed) const {
 
176
 
 
177
  SciFiHelicalPRTrack* helical = static_cast<SciFiHelicalPRTrack*>(seed->getPRTrackTobject());
 
178
 
 
179
  Kalman::Track data_track = BuildTrack(helical, Globals::GetSciFiGeometryHelper(), 5);
 
180
  Kalman::State kalman_seed(seed->getStateVector(), seed->getCovariance());
 
181
 
 
182
  std::cerr << "SEED\n";
 
183
  seed->getStateVector().Print();
 
184
  seed->getCovariance().Print();
 
185
 
 
186
  _helical_track_fitter->SetTrack(data_track);
 
187
  _helical_track_fitter->SetSeed(kalman_seed);
 
188
 
 
189
  _helical_track_fitter->Filter(false);
 
190
  _helical_track_fitter->Smooth(false);
 
191
 
 
192
  std::cerr << "SMOOTHED\n";
 
193
  Kalman::Track new_track = _helical_track_fitter->GetTrack();
 
194
  new_track[14].GetSmoothed().GetVector().Print();
 
195
 
 
196
  SciFiTrack* track = ConvertToSciFiTrack(_helical_track_fitter,
 
197
                                                       Globals::GetSciFiGeometryHelper(), helical);
 
198
  calculate_track_rating(track);
 
199
  track->set_scifi_seed_tobject(seed);
 
200
 
 
201
  ThreeVector seed_pos = track->GetSeedPosition();
 
202
  ThreeVector seed_mom = track->GetSeedMomentum();
 
203
 
 
204
  return track;
 
205
}
 
206
 
 
207
 
 
208
SciFiTrack* MapCppTrackerIntegratedTrackFit::track_fit_straight(SciFiSeed* seed) const {
 
209
 
 
210
  SciFiStraightPRTrack* straight = static_cast<SciFiStraightPRTrack*>(seed->getPRTrackTobject());
 
211
 
 
212
  Kalman::Track data_track = BuildTrack(straight, Globals::GetSciFiGeometryHelper(), 4);
 
213
  Kalman::State kalman_seed(seed->getStateVector(), seed->getCovariance());
 
214
 
 
215
  _straight_track_fitter->SetTrack(data_track);
 
216
  _straight_track_fitter->SetSeed(kalman_seed);
 
217
 
 
218
  _straight_track_fitter->Filter(false);
 
219
  _straight_track_fitter->Smooth(false);
 
220
 
 
221
  SciFiTrack* track = ConvertToSciFiTrack(_straight_track_fitter,
 
222
                                                      Globals::GetSciFiGeometryHelper(), straight);
 
223
  calculate_track_rating(track);
 
224
  track->set_scifi_seed_tobject(seed);
 
225
 
 
226
  return track;
 
227
}
 
228
 
 
229
void MapCppTrackerIntegratedTrackFit::calculate_track_rating(SciFiTrack* track) const {
 
230
  SciFiBasePRTrack* pr_track = track->pr_track_pointer();
 
231
  int number_spacepoints = pr_track->get_num_points();
 
232
  int number_trackpoints = track->GetNumberDataPoints();
 
233
 
 
234
  bool excel_numtp = (number_trackpoints >= _excellent_num_trackpoints);
 
235
  bool good_numtp = (number_trackpoints >= _good_num_trackpoints);
 
236
  bool poor_numtp = (number_trackpoints >= _poor_num_trackpoints);
 
237
  bool excel_numsp = (number_spacepoints >= _excellent_num_spacepoints);
 
238
  bool good_numsp = (number_spacepoints >= _good_num_spacepoints);
 
239
  bool poor_numsp = (number_spacepoints >= _poor_num_spacepoints);
 
240
  bool excel_pval = (track->P_value() >= _excellent_p_value);
 
241
  bool good_pval = (track->P_value() >= _good_p_value);
 
242
  bool poor_pval = (track->P_value() >= _poor_p_value);
 
243
 
 
244
  int rating = 0;
 
245
 
 
246
  if (excel_numtp && excel_numsp && excel_pval) {
 
247
    rating = 1;
 
248
  } else if (good_numtp && good_numsp && good_pval) {
 
249
    rating = 2;
 
250
  } else if (poor_numtp && poor_numsp && poor_pval) {
 
251
    rating = 3;
 
252
  } else {
 
253
    rating = 5;
 
254
  }
 
255
 
 
256
  track->SetRating(rating);
 
257
}
 
258
 
 
259
} // ~namespace MAUS