← Back to branch summary

~mbogomilov/maus/devel3

« back to all changes in this revision

Viewing changes to src/common_cpp/Recon/Global/DataStructureHelper.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:
src/common_cpp/Recon/Global/DataStructureHelper.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 <iostream>
 
21
#include <sstream>
 
22
#include <map>
 
23
#include <vector>
 
24
 
 
25
#include "TLorentzVector.h"
 
26
 
 
27
#include "DataStructure/GlobalEvent.hh"
 
28
#include "DataStructure/Primary.hh"
 
29
#include "DataStructure/ThreeVector.hh"
 
30
#include "DataStructure/Global/PrimaryChain.hh"
 
31
#include "DataStructure/Global/Track.hh"
 
32
#include "DataStructure/Global/TrackPoint.hh"
 
33
#include "DataStructure/Global/ReconEnums.hh"
 
34
#include "JsonCppProcessors/GlobalEventProcessor.hh"
 
35
#include "Recon/Global/DataStructureHelper.hh"
 
36
#include "Recon/Global/Detector.hh"
 
37
#include "Recon/Global/Particle.hh"
 
38
#include "src/common_cpp/Simulation/MAUSPrimaryGeneratorAction.hh"
 
39
#include "Utils/Globals.hh"
 
40
#include "Config/MiceModule.hh"
 
41
 
 
42
// legacy classes outside the MAUS namespace
 
43
namespace MAUS {
 
44
namespace recon {
 
45
namespace global {
 
46
 
 
47
using MAUS::PhaseSpaceVector;
 
48
using MAUS::DataStructure::Global::DetectorPoint;
 
49
using MAUS::DataStructure::Global::PID;
 
50
using MAUS::DataStructure::Global::PrimaryChain;
 
51
using MAUS::DataStructure::Global::Track;
 
52
using MAUS::DataStructure::Global::TrackPoint;
 
53
using MAUS::recon::global::Detector;
 
54
using MAUS::recon::global::Particle;
 
55
 
 
56
const DataStructureHelper& DataStructureHelper::GetInstance() {
 
57
  static DataStructureHelper instance;
 
58
  return instance;
 
59
}
 
60
 
 
61
// This differes from the function in MiceModule in that it searches based on
 
62
// just the local name, not the full node path name. The former is just a lazy
 
63
// way of doing a tree traversal without parsing the node path name. This
 
64
// function has no knowledge of where the node is located and will return nodes
 
65
// of the same local name under different mother nodes.
 
66
std::vector<const MiceModule *> DataStructureHelper::FindModulesByName(
 
67
    const MiceModule * module,
 
68
    std::string name) const {
 
69
  std::vector<const MiceModule*> modules;
 
70
 
 
71
  if (module->name() == name)
 
72
    modules.push_back(module);
 
73
 
 
74
  for (int index = 0; index < module->daughters(); ++index) {
 
75
    std::vector<const MiceModule *> tmp
 
76
      = FindModulesByName(module->daughter(index), name);
 
77
    std::vector<const MiceModule *>::iterator daughter_module;
 
78
    for (daughter_module = tmp.begin();
 
79
         daughter_module != tmp.end();
 
80
         ++daughter_module) {
 
81
      modules.push_back(*daughter_module);
 
82
    }
 
83
  }
 
84
 
 
85
  return modules;
 
86
}
 
87
 
 
88
double DataStructureHelper::GetDetectorZPosition(
 
89
    const GlobalDS::DetectorPoint detector_id) const {
 
90
  MiceModule const * const geometry
 
91
    = Globals::GetInstance()->GetReconstructionMiceModules();
 
92
  std::vector<const MiceModule *> modules;
 
93
 
 
94
  std::stringstream detector_name;
 
95
  std::cout << "Root MiceModule " << geometry->fullName() << " has "
 
96
            << geometry->daughters() << " daughters." << std::endl;
 
97
  std::cout << "Selecting on detector " << detector_id << std::endl;
 
98
  switch (detector_id) {
 
99
    case GlobalDS::kTOF0: {
 
100
      detector_name << "TOF0.dat";
 
101
      modules = FindModulesByName(geometry, detector_name.str());
 
102
      break;
 
103
    }
 
104
    case GlobalDS::kTOF1: {
 
105
      detector_name << "TOF1Detector.dat";
 
106
      modules = FindModulesByName(geometry, detector_name.str());
 
107
      break;
 
108
    }
 
109
    case GlobalDS::kTOF2: {
 
110
      detector_name << "TOF2Detector.dat";
 
111
      modules = FindModulesByName(geometry, detector_name.str());
 
112
      break;
 
113
    }
 
114
    case GlobalDS::kTracker0_1:
 
115
    case GlobalDS::kTracker0_2:
 
116
    case GlobalDS::kTracker0_3:
 
117
    case GlobalDS::kTracker0_4:
 
118
    case GlobalDS::kTracker0_5:
 
119
    case GlobalDS::kTracker1_1:
 
120
    case GlobalDS::kTracker1_2:
 
121
    case GlobalDS::kTracker1_3:
 
122
    case GlobalDS::kTracker1_4:
 
123
    case GlobalDS::kTracker1_5: {
 
124
      GlobalDS::DetectorPoint station = GlobalDS::DetectorPoint(
 
125
        detector_id - GlobalDS::kTracker0);
 
126
      if (station > 5) {  // 5 stations per tracker
 
127
        station = GlobalDS::DetectorPoint(detector_id - GlobalDS::kTracker1);
 
128
        detector_name << "Tracker1Station";
 
129
      } else {
 
130
        detector_name << "TrackerStation";
 
131
      }
 
132
      detector_name << station << ".dat";
 
133
      std::vector<const MiceModule *> mothers
 
134
        = FindModulesByName(geometry, detector_name.str());
 
135
      if (mothers.size() == 1) {
 
136
        std::string view_name = "TrackerViewW.dat";
 
137
        if (detector_id == GlobalDS::kTracker0_5) {
 
138
          view_name = "Tracker0Station5ViewW.dat";
 
139
        }
 
140
        modules = FindModulesByName(mothers[0], view_name);
 
141
        detector_name << "/" << view_name;  // for exception message if needed
 
142
      } else if (mothers.size() > 1) {
 
143
        std::stringstream message;
 
144
        message << "Found multiple detector geometry modules named \""
 
145
                << detector_name.str() << "\".";
 
146
        throw(Exception(Exception::recoverable,
 
147
                      message.str(),
 
148
                      "DataStructureHelper::GetDetectorZPosition()"));
 
149
      } else {
 
150
        std::stringstream message;
 
151
        message << "Couldn't find detector geometry module \""
 
152
                << detector_name.str() << "\".";
 
153
        throw(Exception(Exception::recoverable,
 
154
                      message.str(),
 
155
                      "DataStructureHelper::GetDetectorZPosition()"));
 
156
      }
 
157
      break;
 
158
    }
 
159
    default: detector_name << "unknown";
 
160
  }
 
161
 
 
162
  std::cout << "Detector name: " << detector_name.str() << std::endl;
 
163
 
 
164
  if (modules.size() == 0) {
 
165
    std::stringstream message;
 
166
    message << "Couldn't find reconstruction mapping detector \""
 
167
            << detector_name.str() << "\".";
 
168
    throw(Exception(Exception::recoverable,
 
169
                  message.str(),
 
170
                  "DataStructureHelper::GetDetectorZPosition()"));
 
171
  } else if (modules.size() > 1) {
 
172
    std::stringstream message;
 
173
    message << "Found multiple reconstruction mapping detectors named \""
 
174
            << detector_name.str() << "\".";
 
175
    throw(Exception(Exception::recoverable,
 
176
                  message.str(),
 
177
                  "DataStructureHelper::GetDetectorZPosition()"));
 
178
  }
 
179
 
 
180
  std::cout << "Detector z position: " << modules[0]->globalPosition().z() << std::endl;
 
181
 
 
182
  return modules[0]->globalPosition().z();
 
183
}
 
184
 
 
185
void DataStructureHelper::GetDetectorAttributes(
 
186
    const Json::Value& json_document,
 
187
    std::map<DetectorPoint, Detector>& detectors) const {
 
188
  // FIXME(plane1@hawk.iit.edu) Once the detector groups provide this
 
189
  // information this will need to be changed
 
190
  Json::Value detector_attributes_json = JsonWrapper::GetProperty(
 
191
      json_document, "global_recon_detector_attributes",
 
192
      JsonWrapper::arrayValue);
 
193
 
 
194
  // *** Get detector info ***
 
195
  const Json::Value::UInt detector_count = detector_attributes_json.size();
 
196
  for (Json::Value::UInt index = 0; index < detector_count; ++index) {
 
197
    const Json::Value detector_json = detector_attributes_json[index];
 
198
    const Json::Value id_json = JsonWrapper::GetProperty(
 
199
        detector_json, "id", JsonWrapper::intValue);
 
200
    const DetectorPoint id = DetectorPoint(id_json.asInt());
 
201
 
 
202
    const Json::Value uncertainties_json = JsonWrapper::GetProperty(
 
203
        detector_json, "uncertainties", JsonWrapper::arrayValue);
 
204
    const CovarianceMatrix uncertainties
 
205
        = GetJsonCovarianceMatrix(uncertainties_json);
 
206
 
 
207
    const Detector detector(id, uncertainties);
 
208
    detectors.insert(std::pair<DetectorPoint, Detector>(id, detector));
 
209
  }
 
210
}
 
211
 
 
212
PhaseSpaceVector DataStructureHelper::TrackPoint2PhaseSpaceVector(
 
213
    const TrackPoint& track_point) const {
 
214
  TLorentzVector position = track_point.get_position();
 
215
  TLorentzVector momentum = track_point.get_momentum();
 
216
  return MAUS::PhaseSpaceVector(
 
217
    position.T(), momentum.E(),
 
218
    position.X(), momentum.Px(),
 
219
    position.Y(), momentum.Py());
 
220
}
 
221
 
 
222
TrackPoint DataStructureHelper::PhaseSpaceVector2TrackPoint(
 
223
      const PhaseSpaceVector& vector,
 
224
      const double z,
 
225
      const PID particle_id) const {
 
226
    TrackPoint track_point;
 
227
    const TLorentzVector position(vector.x(), vector.y(), z, vector.t());
 
228
    track_point.set_position(position);
 
229
 
 
230
    const double energy = vector.E();
 
231
    const double px = vector.Px();
 
232
    const double py = vector.Py();
 
233
    const double mass = Particle::GetInstance().GetMass(particle_id);
 
234
    double pz = ::sqrt(energy*energy - mass*mass - px*px - py*py);
 
235
    std::cout << "DEBUG DataStructureHelper::PhaseSpaceVector2TrackPoint: "
 
236
              << "mass: " << mass
 
237
              << "\t: " << vector.t() << "\tE: " << energy
 
238
              << "\tx: " << vector.x() << "\tPx: " << px
 
239
              << "\ty: " << vector.y() << "\tPy: " << py
 
240
              << "\tz: " << z << "\tPz: " << pz << std::endl;
 
241
    if (pz != pz) {
 
242
      pz = 0.;
 
243
    }
 
244
 
 
245
    const TLorentzVector momentum(px, py, pz, energy);
 
246
    track_point.set_momentum(momentum);
 
247
 
 
248
    track_point.set_detector(MAUS::DataStructure::Global::kUndefined);
 
249
    return track_point;
 
250
}
 
251
 
 
252
std::vector<PrimaryChain*>* DataStructureHelper::GetPrimaryChains(
 
253
    const Json::Value& recon_event) const {
 
254
  Json::Value global_event_json = JsonWrapper::GetProperty(
 
255
      recon_event, "global_event", JsonWrapper::objectValue);
 
256
  MAUS::GlobalEventProcessor deserializer;
 
257
  GlobalEvent * global_event = deserializer.JsonToCpp(global_event_json);
 
258
 
 
259
  return global_event->get_primary_chains();
 
260
}
 
261
 
 
262
MAUS::Primary DataStructureHelper::PGParticle2Primary(
 
263
    MAUS::MAUSPrimaryGeneratorAction::PGParticle& pgparticle) const {
 
264
  ThreeVector position(pgparticle.x,
 
265
                       pgparticle.y,
 
266
                       pgparticle.z);
 
267
  ThreeVector momentum(pgparticle.px,
 
268
                       pgparticle.py,
 
269
                       pgparticle.pz);
 
270
  Primary primary;
 
271
  primary.SetPosition(position);
 
272
  primary.SetMomentum(momentum);
 
273
  primary.SetTime(pgparticle.time);
 
274
  primary.SetEnergy(pgparticle.energy);
 
275
  primary.SetParticleId(pgparticle.pid);
 
276
  primary.SetRandomSeed(pgparticle.seed);
 
277
 
 
278
  return primary;
 
279
}
 
280
 
 
281
CovarianceMatrix DataStructureHelper::GetJsonCovarianceMatrix(
 
282
    const Json::Value& value) const {
 
283
  if (value.size() < static_cast<Json::Value::UInt>(6)) {
 
284
    throw(Exception(Exception::recoverable,
 
285
                 "Not enough row elements to convert JSON to CovarianceMatrix",
 
286
                 "DataStructureHelper::GetJsonCovarianceMatrix()"));
 
287
  }
 
288
 
 
289
  const size_t size = 6;
 
290
  double matrix_data[size*size];
 
291
  for (size_t row = 0; row < 6; ++row) {
 
292
    const Json::Value row_json = value[row];
 
293
    if (row_json.size() < static_cast<Json::Value::UInt>(6)) {
 
294
      throw(Exception(
 
295
          Exception::recoverable,
 
296
          "Not enough column elements to convert JSON to CovarianceMatrix",
 
297
          "DataStructureHelper::GetJsonCovarianceMatrix()"));
 
298
    }
 
299
    for (size_t column = 0; column < 6; ++column) {
 
300
      const Json::Value element_json = row_json[column];
 
301
      matrix_data[row*size+column] = element_json.asDouble();
 
302
    }
 
303
  }
 
304
 
 
305
  CovarianceMatrix matrix(matrix_data);
 
306
  return matrix;
 
307
}
 
308
 
 
309
}  // namespace global
 
310
}  // namespace recon
 
311
}  // namespace MAUS
 
312