← Back to branch summary

~chris-rogers/maus/emr_mc_digitization

« back to all changes in this revision

Viewing changes to src/common_cpp/Optics/PolynomialOpticsModel.cc

  • Committer: Chris Rogers
  • Date: 2014-04-16 11:48:45 UTC
  • mfrom: (707 merge)
  • mto: This revision was merged to the branch mainline in revision 711.
  • Revision ID: chris.rogers@stfc.ac.uk-20140416114845-h3u3q7pdcxkxvovs
Update to trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Optics/PolynomialOpticsModel.cc
20
20
 
21
21
#include "Optics/PolynomialOpticsModel.hh"
22
22
 
 
23
#include <algorithm>
23
24
#include <cfloat>
24
25
#include <cmath>
25
26
#include <map>
26
27
 
27
28
#include <iomanip>
 
29
#include <list>
28
30
#include <sstream>
29
31
#include <vector>
30
32
 
31
 
#include "Interface/Squeal.hh"
 
33
#include "Utils/Exception.hh"
32
34
#include "Maths/PolynomialMap.hh"
 
35
#include "src/common_cpp/DataStructure/Primary.hh"
 
36
#include "src/common_cpp/DataStructure/ThreeVector.hh"
 
37
#include "src/common_cpp/JsonCppProcessors/PrimaryProcessor.hh"
 
38
#include "src/common_cpp/Optics/PhaseSpaceVector.hh"
33
39
#include "src/common_cpp/Optics/PolynomialTransferMap.hh"
34
 
#include "Recon/Global/TrackPoint.hh"
 
40
#include "Recon/Global/DataStructureHelper.hh"
 
41
#include "Recon/Global/ParticleOpticalVector.hh"
35
42
#include "Simulation/MAUSGeant4Manager.hh"
36
43
#include "Simulation/MAUSPhysicsList.hh"
37
44
 
38
45
namespace MAUS {
39
46
 
40
 
using recon::global::TrackPoint;
 
47
using MAUS::PhaseSpaceVector;
 
48
using recon::global::ParticleOpticalVector;
41
49
 
42
 
PolynomialOpticsModel::PolynomialOpticsModel(const Json::Value & configuration)
43
 
      : TransferMapOpticsModel(configuration), algorithm_(kNone) {
 
50
PolynomialOpticsModel::PolynomialOpticsModel(
 
51
      Json::Value const * const configuration)
 
52
      : TransferMapOpticsModel(configuration),
 
53
        algorithm_(kNone) {
44
54
  // Determine which fitting algorithm to use
45
55
  SetupAlgorithm();
46
56
 
55
65
 
56
66
void PolynomialOpticsModel::Build() {
57
67
  // Create some test hits at the desired First plane
58
 
  const std::vector<TrackPoint> first_plane_hits = BuildFirstPlaneHits();
59
 
 
60
 
  // Iterate through each First plane hit
 
68
  const std::vector<PhaseSpaceVector> primary_vectors = PrimaryVectors();
 
69
std::cout << "DEBUG PolynomialOpticsModel::Build: "
 
70
          << "# primaries: " << primary_vectors.size() << std::endl;
 
71
  Json::Value primaries_json;
 
72
  for (std::vector<PhaseSpaceVector>::const_iterator primary_vector
 
73
        = primary_vectors.begin();
 
74
       primary_vector < primary_vectors.end();
 
75
       ++primary_vector) {
 
76
    // generate a Primary object
 
77
    Primary primary(reference_primary_);
 
78
    ThreeVector position(primary_vector->x(), primary_vector->y(),
 
79
                         reference_primary_.GetPosition().z());
 
80
    ThreeVector momentum(primary_vector->Px(), primary_vector->Py(),
 
81
                         reference_primary_.GetMomentum().z());
 
82
    primary.SetPosition(position);
 
83
    primary.SetMomentum(momentum);
 
84
    primary.SetTime(primary_vector->t());
 
85
    primary.SetEnergy(primary_vector->E());
 
86
 
 
87
    // serialize primary to JSON
 
88
    PrimaryProcessor serializer;
 
89
    Json::Value * primary_json = serializer.CppToJson(primary, "");
 
90
    Json::Value object_value;
 
91
    object_value["primary"] = *primary_json;
 
92
    primaries_json.append(object_value);
 
93
  }
 
94
 
61
95
  MAUSGeant4Manager * simulator = MAUSGeant4Manager::GetInstance();
 
96
 
 
97
  // Force setting of stochastics
62
98
  simulator->GetPhysicsList()->BeginOfRunAction();
63
99
 
64
 
  std::map<int, std::vector<TrackPoint> > station_hits_map;
65
 
  std::vector<TrackPoint>::const_iterator first_plane_hit;
66
 
  for (first_plane_hit = first_plane_hits.begin();
67
 
       first_plane_hit < first_plane_hits.end();
68
 
       ++first_plane_hit) {
69
 
    // Simulate the current particle (First plane hit) through MICE.
70
 
    simulator->RunParticle(
71
 
      recon::global::PrimaryGeneratorParticle(*first_plane_hit));
72
 
 
73
 
    // Identify the hits by station and add them to the mappings from stations
74
 
    // to the hits they recorded.
75
 
    MapStationsToHits(station_hits_map);
76
 
  }
77
 
 
78
 
  // Iterate through each station
79
 
  std::map<int, std::vector<TrackPoint> >::iterator station_hits;
 
100
  // Simulate on the primaries, generating virtual detector tracks for each
 
101
  const Json::Value virtual_tracks
 
102
      = MAUSGeant4Manager::GetInstance()->RunManyParticles(primaries_json);
 
103
std::cout << "DEBUG PolynomialOpticsModel::Build: "
 
104
          << "# virtual tracks: " << virtual_tracks.size() << std::endl;
 
105
  if (virtual_tracks.size() == 0) {
 
106
    throw(Exception(Exception::nonRecoverable,
 
107
                 "No events were generated during simulation.",
 
108
                 "MAUS::TransferMapOpticsModel::Build()"));
 
109
  }
 
110
 
 
111
  // Map stations to hits in each virtual track
 
112
  std::map<int64_t, std::vector<PhaseSpaceVector> > station_hits_map;
 
113
size_t count = 0;
 
114
  for (Json::Value::const_iterator virtual_track = virtual_tracks.begin();
 
115
       virtual_track != virtual_tracks.end();
 
116
       ++virtual_track) {
 
117
    MapStationsToHits(station_hits_map, *virtual_track);
 
118
    ++count;
 
119
  }
 
120
std::cout << "DEBUG PolynomialOpticsModel::Build:" << std::endl
 
121
          << "\t# virtual tracks mapped: " << count << std::endl
 
122
          << "\t# station Zs: " << station_hits_map.size() << std::endl;
 
123
 
 
124
  // Calculate transfer maps from the primary plane to each station plane
 
125
  std::map<int64_t, std::vector<PhaseSpaceVector> >::iterator station_hits;
80
126
  for (station_hits = station_hits_map.begin();
81
127
       station_hits != station_hits_map.end();
82
128
       ++station_hits) {
83
 
    // find the average z coordinate for the station
84
 
    std::vector<TrackPoint>::iterator station_hit;
85
 
 
86
 
    double station_plane = station_hits->second.begin()->z();
87
 
 
88
 
    // Generate a transfer map between the First plane and the current station
89
 
    // and map the station ID to the transfer map
90
 
    transfer_maps_[station_plane]
91
 
      = CalculateTransferMap(first_plane_hits, station_hits->second);
92
 
  }
 
129
    // calculate transfer map and index it by the station z-plane
 
130
std::cout << "DEBUG PolynomialOpticsModel::Build: "
 
131
          << "# virtual track hits for z = " << station_hits->first
 
132
          << ": " << station_hits->second.size() << std::endl;
 
133
    transfer_maps_[station_hits->first]
 
134
      = CalculateTransferMap(primary_vectors, station_hits->second);
 
135
  }
 
136
 
 
137
  built_ = true;
 
138
}
 
139
 
 
140
const std::vector<int64_t> PolynomialOpticsModel::GetAvailableMapPositions()
 
141
    const {
 
142
  if (!built_) {
 
143
    throw(Exception(Exception::nonRecoverable,
 
144
                  "No transfer maps available since the optics model has not "
 
145
                  "been built yet. Call Build() first.",
 
146
                  "MAUS::PolynomialOpticsModel::GetAvailableMapPositions()"));
 
147
  }
 
148
 
 
149
  std::vector<int64_t> positions;
 
150
  std::map<int64_t, const TransferMap *>::const_iterator maps;
 
151
  // insertion sort the map keys (z-positions)
 
152
  for (maps = transfer_maps_.begin(); maps != transfer_maps_.end(); ++maps) {
 
153
    positions.push_back(maps->first);
 
154
  }
 
155
  std::sort(positions.begin(), positions.end());  // just in case
 
156
 
 
157
  return positions;
93
158
}
94
159
 
95
160
void PolynomialOpticsModel::SetupAlgorithm() {
134
199
 * necessary in order to solve the least squares problem which involves
135
200
 * the calculation of a Moore-Penrose psuedo inverse.
136
201
 */
137
 
const std::vector<TrackPoint> PolynomialOpticsModel::BuildFirstPlaneHits() {
 
202
const std::vector<PhaseSpaceVector> PolynomialOpticsModel::PrimaryVectors() {
138
203
  size_t num_poly_coefficients
139
204
    = PolynomialMap::NumberOfPolynomialCoefficients(6, polynomial_order_);
140
 
 
141
 
    std::vector<TrackPoint> first_plane_hits;
 
205
  std::vector<PhaseSpaceVector> primaries;
 
206
 
 
207
  // The (0,0,0,0,0,0) case produces a polynomial vector (1,0,0,...,0)
 
208
  // primaries.push_back(reference_trajectory_);
 
209
 
 
210
  std::cerr << "Primaries:" << std::endl;
142
211
  for (size_t i = 0; i < 6; ++i) {
143
212
    for (size_t j = i; j < 6; ++j) {
144
 
      TrackPoint first_plane_hit;
 
213
      PhaseSpaceVector primary;
145
214
 
146
215
      for (size_t k = 0; k < i; ++k) {
147
 
        first_plane_hit[k] = 1.;
 
216
        primary[k] = deltas_[k];  // non-zero, lower triangular elements
148
217
      }
149
 
      double delta = deltas_[j];
150
 
      first_plane_hit[j] = delta;
 
218
      primary[j] = deltas_[j];  // diagonal element
151
219
 
152
 
      first_plane_hits.push_back(TrackPoint(first_plane_hit + reference_particle_,
153
 
                                 reference_particle_.z(),
154
 
                                 reference_particle_.particle_id()));
 
220
      // std::cerr << (primary + reference_trajectory_) << std::endl;
 
221
      std::cerr << primary << std::endl;
 
222
      primaries.push_back(primary + reference_trajectory_);
155
223
    }
156
224
  }
157
 
  size_t base_block_length = first_plane_hits.size();
 
225
  size_t base_block_length = primaries.size();
158
226
 
159
227
  // adjust for the case where the base block size is greater than N(n, v)
160
228
  while (base_block_length > num_poly_coefficients) {
161
 
    first_plane_hits.pop_back();
162
 
    base_block_length = first_plane_hits.size();
 
229
    primaries.pop_back();
 
230
    base_block_length = primaries.size();
163
231
  }
164
232
 
165
 
  int summand;
166
233
  for (size_t row = base_block_length; row < num_poly_coefficients; ++row) {
167
 
    summand = row / base_block_length;
168
 
    TrackPoint first_plane_hit
169
 
      = TrackPoint(first_plane_hits[row % base_block_length] + summand);
 
234
    PhaseSpaceVector deltas = deltas_ * (row / base_block_length);
 
235
    PhaseSpaceVector primary
 
236
      = primaries[row % base_block_length] - reference_trajectory_ + deltas;
170
237
 
171
 
    first_plane_hits.push_back(TrackPoint(first_plane_hit + reference_particle_,
172
 
                                reference_particle_.z(),
173
 
                                reference_particle_.particle_id()));
 
238
    // std::cerr << (primary + reference_trajectory_) << std::endl;
 
239
    std::cerr << primary << std::endl;
 
240
    primaries.push_back(primary + reference_trajectory_);
174
241
  }
175
 
  return first_plane_hits;
 
242
  return primaries;
176
243
}
177
244
 
178
245
/* Calculate a transfer matrix from an equal number of inputs and output.
179
246
 */
180
247
const TransferMap * PolynomialOpticsModel::CalculateTransferMap(
181
 
    const std::vector<recon::global::TrackPoint> & start_plane_hits,
182
 
    const std::vector<recon::global::TrackPoint> & station_hits)
 
248
    const std::vector<PhaseSpaceVector> & start_plane_hits,
 
249
    const std::vector<PhaseSpaceVector> & station_hits)
183
250
    const {
 
251
  #if 0
 
252
  MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
 
253
  MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
 
254
    = simulator->GetReferenceParticle();
 
255
  const double t0 = reference_pgparticle.time;
 
256
  const double E0 = reference_pgparticle.energy;
 
257
  const double P0 = reference_pgparticle.pz;
 
258
  #endif
184
259
 
185
260
  if (start_plane_hits.size() != station_hits.size()) {
186
 
    throw(Squeal(Squeal::nonRecoverable,
187
 
                  "The number of start plane hits is not the same as the number "
188
 
                  "of hits per station.",
 
261
    std::stringstream message;
 
262
    message << "The number of start plane hits (" << start_plane_hits.size()
 
263
            << ") is not the same as the number of hits per station ("
 
264
            << station_hits.size() << ").";
 
265
    throw(Exception(Exception::nonRecoverable,
 
266
                  message.str(),
189
267
                  "PolynomialOpticsModel::CalculateTransferMap()"));
190
268
  }
191
269
 
192
270
  std::vector< std::vector<double> > points;
193
271
  for (size_t pt_index = 0; pt_index < start_plane_hits.size(); ++pt_index) {
 
272
    /*
 
273
    ParticleOpticalVector start_plane_point(start_plane_hits[pt_index],
 
274
                                            t0, E0, P0);
 
275
    */
 
276
    PhaseSpaceVector start_plane_point(
 
277
      start_plane_hits[pt_index] - reference_trajectory_);
 
278
 
194
279
    std::vector<double> point;
195
280
    for (size_t coord_index = 0; coord_index < 6; ++coord_index) {
196
 
      point.push_back(start_plane_hits[pt_index][coord_index]);
 
281
      #if 1
 
282
        point.push_back(start_plane_point[coord_index]);
 
283
      #else
 
284
        point.push_back(start_plane_hits[pt_index][coord_index]);
 
285
      #endif
197
286
    }
198
287
    points.push_back(point);
199
288
  }
200
289
 
201
290
  std::vector< std::vector<double> > values;
202
291
  for (size_t val_index = 0; val_index < station_hits.size(); ++val_index) {
 
292
    /*
 
293
    ParticleOpticalVector station_value(station_hits[val_index], t0, E0, P0);
 
294
    */
 
295
    PhaseSpaceVector station_value(
 
296
      station_hits[val_index] - reference_trajectory_);
 
297
 
203
298
    std::vector<double> value;
204
299
    for (size_t coord_index = 0; coord_index < 6; ++coord_index) {
205
 
      value.push_back(station_hits[val_index][coord_index]);
 
300
      #if 1
 
301
        value.push_back(station_value[coord_index]);
 
302
      #else
 
303
        value.push_back(station_hits[val_index][coord_index]);
 
304
      #endif
206
305
    }
207
306
    values.push_back(value);
208
307
  }
211
310
 
212
311
  switch (algorithm_) {
213
312
    case kNone:
214
 
      throw(Squeal(Squeal::nonRecoverable,
 
313
      throw(Exception(Exception::nonRecoverable,
215
314
                    "No fitting algorithm specified in configuration.",
216
315
                    "PolynomialOpticsModel::CalculateTransferMap()"));
217
316
    case kLeastSquares:
218
317
      // Fit to first order and then fit to higher orders with the
219
318
      // first order map as a constraint
220
319
      polynomial_map = PolynomialMap::PolynomialLeastSquaresFit(
 
320
          points, values, polynomial_order_, weights_);
 
321
      /*
221
322
          points, values, 1, weights_);
222
323
      if (polynomial_order_ > 1) {
223
324
        PolynomialMap * linear_polynomial_map = polynomial_map;
225
326
            points, values, polynomial_order_,
226
327
            linear_polynomial_map->GetCoefficientsAsVector(), weights_);
227
328
      }
 
329
      */
228
330
      break;
229
331
    case kConstrainedLeastSquares:
230
332
      // constrained least squares
231
333
      // ConstrainedLeastSquaresFit(...);
232
 
      throw(Squeal(Squeal::nonRecoverable,
 
334
      throw(Exception(Exception::nonRecoverable,
233
335
                    "Constrained Polynomial fitting algorithm "
234
336
                    "is not yet implemented.",
235
337
                    "PolynomialOpticsModel::CalculateTransferMap()"));
236
338
    case kConstrainedChiSquared:
237
339
      // constrained chi squared least squares
238
340
      // Chi2ConstrainedLeastSquaresFit(...);
239
 
      throw(Squeal(Squeal::nonRecoverable,
 
341
      throw(Exception(Exception::nonRecoverable,
240
342
                    "Constrained Chi Squared fitting algorithm "
241
343
                    "is not yet implemented.",
242
344
                    "PolynomialOpticsModel::CalculateTransferMap()"));
243
345
    case kSweepingChiSquared:
244
346
      // sweeping chi squared least squares
245
347
      // Chi2SweepingLeastSquaresFit(...);
246
 
      throw(Squeal(Squeal::nonRecoverable,
 
348
      throw(Exception(Exception::nonRecoverable,
247
349
                    "Sweeping Chi Squared fitting algorithm "
248
350
                    "is not yet implemented.",
249
351
                    "PolynomialOpticsModel::CalculateTransferMap()"));
250
352
    case kSweepingChiSquaredWithVariableWalls:
251
353
      // sweeping chi squared with variable walls
252
354
      // Chi2SweepingLeastSquaresFitVariableWalls(...);
253
 
      throw(Squeal(Squeal::nonRecoverable,
 
355
      throw(Exception(Exception::nonRecoverable,
254
356
                    "Sweeping Chi Squared Variable Walls fitting algorithm "
255
357
                    "is not yet implemented.",
256
358
                    "PolynomialOpticsModel::CalculateTransferMap()"));
257
359
    default:
258
 
      throw(Squeal(Squeal::nonRecoverable,
 
360
      throw(Exception(Exception::nonRecoverable,
259
361
                    "Unrecognized fitting algorithm in configuration.",
260
362
                    "PolynomialOpticsModel::CalculateTransferMap()"));
261
363
  }
262
364
 
263
365
  TransferMap * transfer_map = new PolynomialTransferMap(
264
 
    *polynomial_map, reference_particle_);
 
366
    *polynomial_map, reference_trajectory_);
265
367
  delete polynomial_map;
266
368
 
267
369
  return transfer_map;