1
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
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.
8
* MAUS is distributedTrackReconstructor 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.
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/>.
20
#include "src/map/MapCppGlobalRawTracks/MapCppGlobalRawTracks.hh"
30
#include "TLorentzVector.h"
32
#include "CLHEP/Random/RandGauss.h"
33
#include "CLHEP/Units/PhysicalConstants.h"
34
#include "json/json.h"
37
#include "Config/MiceModule.hh"
38
#include "Utils/Exception.hh"
41
#include "Converter/DataConverters/JsonCppSpillConverter.hh"
42
#include "Converter/DataConverters/CppJsonSpillConverter.hh"
43
#include "DataStructure/Data.hh"
44
#include "DataStructure/GlobalEvent.hh"
45
#include "DataStructure/MCEvent.hh"
46
#include "DataStructure/ReconEvent.hh"
47
#include "DataStructure/SciFiTrackPoint.hh"
48
#include "DataStructure/TOFEventSpacePoint.hh"
49
#include "DataStructure/ThreeVector.hh"
50
#include "DataStructure/Global/ReconEnums.hh"
51
#include "DataStructure/Global/BasePoint.hh"
52
#include "DataStructure/Global/SpacePoint.hh"
53
#include "DataStructure/Global/Track.hh"
54
#include "DataStructure/Global/TrackPoint.hh"
55
#include "src/common_cpp/Optics/CovarianceMatrix.hh"
56
#include "Recon/Global/DataStructureHelper.hh"
57
#include "Recon/Global/Detector.hh"
58
#include "Recon/Global/Particle.hh"
59
#include "src/common_cpp/Simulation/MAUSGeant4Manager.hh"
60
#include "Utils/Globals.hh"
61
#include "Utils/JsonWrapper.hh"
62
#include "Utils/CppErrorHandler.hh"
66
namespace GlobalDS = ::MAUS::DataStructure::Global;
67
using MAUS::recon::global::DataStructureHelper;
68
using MAUS::recon::global::Detector;
69
using MAUS::recon::global::DetectorMap;
70
using MAUS::recon::global::Particle;
72
MapCppGlobalRawTracks::MapCppGlobalRawTracks() {
75
MapCppGlobalRawTracks::~MapCppGlobalRawTracks() {
78
bool MapCppGlobalRawTracks::birth(std::string configuration_string) {
79
// parse the JSON document.
81
const Json::Value configuration
82
= JsonWrapper::StringToJson(configuration_string);
84
DataStructureHelper::GetInstance().GetDetectorAttributes(
85
configuration, detectors_);
86
} catch(Exception& exception) {
87
MAUS::CppErrorHandler::getInstance()->HandleExceptionNoJson(
88
exception, MapCppGlobalRawTracks::kClassname);
90
} catch(std::exception& exc) {
91
MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(
92
exc, MapCppGlobalRawTracks::kClassname);
96
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
97
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
98
= simulator->GetReferenceParticle();
99
switch (reference_pgparticle.pid) {
100
case MAUS::DataStructure::Global::kEPlus:
101
case MAUS::DataStructure::Global::kMuPlus:
102
case MAUS::DataStructure::Global::kPiPlus:
105
case MAUS::DataStructure::Global::kEMinus:
106
case MAUS::DataStructure::Global::kMuMinus:
107
case MAUS::DataStructure::Global::kPiMinus:
111
throw(Exception(Exception::nonRecoverable,
112
"Reference particle is not a pion+/-, muon+/-, or e+/-.",
113
"MapCppGlobalRawTracks::birth()"));
116
return true; // Sucessful parsing
119
std::string MapCppGlobalRawTracks::process(std::string run_data_string) {
120
// parse the JSON document.
121
Json::Value run_data_json
122
= Json::Value(JsonWrapper::StringToJson(run_data_string));
123
if (run_data_json.isNull() || run_data_json.empty()) {
124
return std::string("{\"errors\":{\"bad_json_document\":"
125
"\"Failed to parse input document\"}}");
128
JsonCppSpillConverter deserialize;
129
MAUS::Data * run_data = deserialize(&run_data_json);
131
return std::string("{\"errors\":{\"failed_json_cpp_conversion\":"
132
"\"Failed to convert Json to C++ Data object\"}}");
135
const MAUS::Spill * spill = run_data->GetSpill();
136
MAUS::ReconEventPArray * recon_events = spill->GetReconEvents();
138
return run_data_string;
141
MAUS::ReconEventPArray::const_iterator recon_event;
142
MAUS::GlobalEvent * global_event;
143
for (recon_event = recon_events->begin();
144
recon_event != recon_events->end();
146
global_event = new GlobalEvent();
147
// Load the ReconEvent, and import it into the GlobalEvent
148
AssembleRawTracks(*recon_event, global_event);
151
// Serialize the Spill for passing on to the next map in the workflow
152
CppJsonSpillConverter serialize;
153
Json::FastWriter writer;
154
std::string output = writer.write(*serialize(run_data));
156
std::cout << "DEBUG MapCppGlobalRawTracks::process(): "
157
<< "Output: " << std::endl
158
<< output << std::endl;
165
bool MapCppGlobalRawTracks::death() {
166
return true; // successful
169
/** AssembleRawTracks - load data provided by the DAQ system or via digitized MC
172
* 1) Deserialize TOF and SciFi events
173
* 2) Assemble all raw track permutations
174
* 3) Prune any raw tracks that are unlickly or impossible
175
* 4) Serialize raw tracks
177
void MapCppGlobalRawTracks::AssembleRawTracks(
178
MAUS::ReconEvent * recon_event,
179
MAUS::GlobalEvent * global_event) {
181
// Load TOF and SciFi tracks from the appropriate recon event trees
182
GlobalDS::TrackPArray tracks;
184
std::cout << "DEBUG MapCppGlobalRawTracks::AssembleRawTracks(): "
185
<< "Loading TOF track..." << std::endl;
187
LoadTOFTrack(recon_event, tracks);
188
// size_t tof_track_count = tracks.size();
190
std::cout << "DEBUG MapCppGlobalRawTracks::AssembleRawTracks(): "
191
<< "Loaded " << tof_track_count << " TOF tracks." << std::endl;
194
GlobalDS::TrackPArray sci_fi_tracks;
195
std::cout << "DEBUG MapCppGlobalRawTracks::AssembleRawTracks(): "
196
<< "Loading SciFi track..." << std::endl;
198
LoadSciFiTracks(recon_event, tracks);
200
size_t sci_fi_track_count = tracks.size() - tof_track_count;
201
std::cout << "DEBUG MapCppGlobalRawTracks::AssembleRawTracks(): "
202
<< "Loaded " << sci_fi_track_count << " SciFi tracks."
205
GlobalDS::TrackPArray::iterator track;
206
for (track = tracks.begin(); track != tracks.end(); ++track) {
207
if ((*track)->GetTrackPoints().size() >= 2) {
208
(*track)->SortTrackPointsByZ();
209
(*track)->set_mapper_name(kClassname);
210
global_event->add_track_recursive(*track);
212
std::cout << "DEBUG MapCppGlobalRawTracks::AssembleRawTracks(): "
213
<< "Skipping track with only " << (*track)->GetTrackPoints().size() << " points."
218
recon_event->SetGlobalEvent(global_event);
223
void MapCppGlobalRawTracks::LoadTOFTrack(
224
MAUS::ReconEvent const * const recon_event,
225
GlobalDS::TrackPArray & tof_tracks) {
226
const int particle_event_number = recon_event->GetPartEventNumber();
227
const TOFEvent * tof_event = recon_event->GetTOFEvent();
228
const TOFEventSpacePoint space_point_events
229
= tof_event->GetTOFEventSpacePoint();
230
const std::vector<TOFSpacePoint> tof0_space_points
231
= space_point_events.GetTOF0SpacePointArray();
232
const std::vector<TOFSpacePoint> tof1_space_points
233
= space_point_events.GetTOF1SpacePointArray();
234
const std::vector<TOFSpacePoint> tof2_space_points
235
= space_point_events.GetTOF2SpacePointArray();
237
// FIXME: Each space point in a particular detector is a separate track branch
239
std::cout << "DEBUG LoadTOFTrack: space points in TOF0: "
240
<< tof0_space_points.size() << "\tTOF1: "
241
<< tof1_space_points.size() << "\tTOF2: "
242
<< tof2_space_points.size() << std::endl;
243
if (tof0_space_points.size() != tof1_space_points.size()) {
244
std::cout << "DEBUG LoadTOFTrack: track has different number of "
245
<< "space points in TOF0 (" << tof0_space_points.size() << ")"
246
<< " and TOF1 (" << tof1_space_points.size() << ")...skipping."
249
} else if (tof0_space_points.size() == 0) {
250
std::cout << "DEBUG LoadTOFTrack: track has no space points...skipping."
255
GlobalDS::Track * track = new GlobalDS::Track();
257
TLorentzVector last_position[3];
258
TLorentzVector deltas[3];
260
GlobalDS::TrackPointPArray track_points;
262
// FIXME(Lane) Replace hard coded slab to coordinate conversion once
263
// Mark Reiner's code has been incorporated into TOF space point code.
264
MAUS::recon::global::DetectorMap::const_iterator tof0_mapping
265
= detectors_.find(GlobalDS::kTOF0);
266
if (tof0_mapping == detectors_.end()) {
267
throw(Exception(Exception::nonRecoverable,
268
"Unable to find info for detector TOF0",
269
"MapCppGlobalRawTracks::LoadTOFTrack()"));
271
const Detector& tof0 = tof0_mapping->second;
272
std::vector<TOFSpacePoint>::const_iterator tof0_space_point;
273
for (tof0_space_point = tof0_space_points.begin();
274
tof0_space_point != tof0_space_points.end();
275
++tof0_space_point) {
276
GlobalDS::TrackPoint * track_point = new GlobalDS::TrackPoint();
277
track_point->set_particle_event(particle_event_number);
278
PopulateTOFTrackPoint(tof0, tof0_space_point, 40., 10, track_point);
279
track_points.push_back(track_point);
281
last_position[0] = track_point->get_position();
284
MAUS::recon::global::DetectorMap::const_iterator tof1_mapping
285
= detectors_.find(GlobalDS::kTOF1);
286
if (tof1_mapping == detectors_.end()) {
287
throw(Exception(Exception::nonRecoverable,
288
"Unable to find info for detector TOF1",
289
"MapCppGlobalRawTracks::LoadTOFTrack()"));
291
const Detector& tof1 = tof1_mapping->second;
292
std::vector<TOFSpacePoint>::const_iterator tof1_space_point;
293
for (tof1_space_point = tof1_space_points.begin();
294
tof1_space_point != tof1_space_points.end();
295
++tof1_space_point) {
296
GlobalDS::TrackPoint * track_point = new GlobalDS::TrackPoint();
297
track_point->set_particle_event(particle_event_number);
298
PopulateTOFTrackPoint(tof1, tof1_space_point, 60., 7, track_point);
299
track_points.push_back(track_point);
301
last_position[1] = track_point->get_position();
303
deltas[0] = last_position[1] - last_position[0];
304
deltas[1] = deltas[0];
306
MAUS::recon::global::DetectorMap::const_iterator tof2_mapping
307
= detectors_.find(GlobalDS::kTOF2);
308
if (tof2_mapping == detectors_.end()) {
309
throw(Exception(Exception::nonRecoverable,
310
"Unable to find info for detector TOF2",
311
"MapCppGlobalRawTracks::LoadTOFTrack()"));
313
const Detector& tof2 = tof2_mapping->second;
314
std::vector<TOFSpacePoint>::const_iterator tof2_space_point;
315
for (tof2_space_point = tof2_space_points.begin();
316
tof2_space_point != tof2_space_points.end();
317
++tof2_space_point) {
318
GlobalDS::TrackPoint * track_point = new GlobalDS::TrackPoint();
319
track_point->set_particle_event(particle_event_number);
320
PopulateTOFTrackPoint(tof2, tof2_space_point, 60., 10, track_point);
321
track_points.push_back(track_point);
323
last_position[2] = track_point->get_position();
325
deltas[2] = last_position[2] - last_position[1];
327
// Approximate momenta by using tof0/tof1 and tof1/tof2 position deltas
328
TLorentzVector momenta[3];
329
for (size_t index = 0; index < 3; ++index) {
330
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
331
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
332
= simulator->GetReferenceParticle();
333
momenta[index] = TLorentzVector(reference_pgparticle.px,
334
reference_pgparticle.py,
335
reference_pgparticle.pz,
336
reference_pgparticle.energy);
339
// Set each track point's 4-momentum and add to the track
340
std::vector<GlobalDS::TrackPoint *>::iterator track_point;
341
size_t point_index = 0;
342
for (track_point = track_points.begin();
343
track_point != track_points.end();
345
(*track_point)->set_momentum(momenta[point_index++]);
347
track->AddTrackPoint(*track_point);
350
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
351
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
352
= simulator->GetReferenceParticle();
353
const GlobalDS::PID particle_id = GlobalDS::PID(reference_pgparticle.pid);
354
track->set_pid(particle_id);
356
tof_tracks.push_back(track);
359
/** PopulateTOFTrackPoint
361
void MAUS::MapCppGlobalRawTracks::PopulateTOFTrackPoint(
362
const Detector & detector,
363
const std::vector<TOFSpacePoint>::const_iterator & tof_space_point,
364
const double slab_width,
365
const size_t number_of_slabs,
366
GlobalDS::TrackPoint * track_point) {
367
DataStructureHelper helper = DataStructureHelper::GetInstance();
368
const double z = helper.GetDetectorZPosition(detector.id());
369
std::cerr << "DEBUG MapCppGlobalRawTracks::PopulateTOFTrackPoint: "
371
<< "\tPhys. Event #: " << tof_space_point->GetPhysEventNumber()
373
std::cerr << "DEBUG MapCppGlobalRawTracks::PopulateTOFTrackPoint: " << std::endl
374
<< "\tID: " << detector.id() << std::endl
375
<< "\tz-position: " << z << std::endl
376
<< "\tUncertainties: " << detector.uncertainties() << std::endl;
377
double max_xy = slab_width * ::floor(number_of_slabs / 2.0);
378
if (number_of_slabs % 2 == 0) {
379
max_xy -= slab_width / 2.0;
381
std::cerr << "DEBUG MapCppGlobalRawTracks::PopulateTOFTrackPoint: " << std::endl
382
<< "\t# Slabs: " << number_of_slabs << std::endl
383
<< "\tSlab Width: " << slab_width << std::endl
384
<< "\tMax X/Y Value: " << max_xy << std::endl;
386
GlobalDS::SpacePoint * space_point = new GlobalDS::SpacePoint();
388
space_point->set_detector(detector.id());
390
// FIXME(Lane) not sure what to put here
391
space_point->set_geometry_path("");
393
// GetSlaby() gets the number of the slab that is *oriented* in the y
394
// direction For TOF0 this is the slab that gives an approximate x coordinate)
395
// and vice versa. For TOF1 and TOF2 it yields the y coordinate.
396
const double x = slab_width * tof_space_point->GetSlaby() - max_xy;
397
const double y = slab_width * tof_space_point->GetSlabx() - max_xy;
398
const double t = tof_space_point->GetTime();
399
TLorentzVector position(x, y, z, t);
400
space_point->set_position(position);
401
static_cast<GlobalDS::BasePoint*>(track_point)->operator=(*space_point);
402
track_point->set_space_point(space_point);
403
track_point->set_mapper_name(kClassname);
405
track_point->set_position(position);
406
CovarianceMatrix covariances = detector.uncertainties();
407
TLorentzVector position_errors(::sqrt(covariances(3, 3)),
408
::sqrt(covariances(5, 5)),
410
::sqrt(covariances(1, 1)));
411
track_point->set_position_error(position_errors);
413
TLorentzVector momentum_errors(::sqrt(covariances(4, 4)),
414
::sqrt(covariances(6, 6)),
416
::sqrt(covariances(2, 2)));
417
track_point->set_momentum_error(momentum_errors);
419
track_point->set_charge(tof_space_point->GetChargeProduct());
424
double MapCppGlobalRawTracks::FindEnergy(const double mass,
425
const double delta_z,
426
const double actual_delta_t) const {
427
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
428
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
429
= simulator->GetReferenceParticle();
430
const double reference_energy = reference_pgparticle.energy;
431
const double slab_length = 25.; // mm
432
const double drift_length = delta_z - 2 * slab_length;
433
const double beta_0 = delta_z / actual_delta_t / ::CLHEP::c_light;
434
const double gamma_0 = 1. / std::sqrt(1 - beta_0*beta_0);
435
const double E_0 = gamma_0 * mass;
439
double E_max = reference_energy + 50.; // allow for momentum spread
444
bool energy_too_low = false;
446
for (size_t iteration = 0; iteration < 10; ++iteration) {
447
double slab_energy_loss = 0.;
448
for (size_t index = 0; index < 4; ++index) {
449
beta[index] = ::sqrt(1 - ::pow(mass / (E - slab_energy_loss), 2));
450
if (beta[index] == 0.) {
451
energy_too_low = true;
453
slab_energy_loss += TOFSlabEnergyLoss(beta[index], mass);
454
if ((slab_energy_loss > E) || (slab_energy_loss != slab_energy_loss)) {
455
energy_too_low = true;
458
delta_t[0] = slab_length / ((beta[0] + beta[1]) / 2) / ::CLHEP::c_light;
459
delta_t[1] = drift_length / beta[1] / ::CLHEP::c_light;
460
delta_t[2] = slab_length / ((beta[1] + beta[2]) / 2) / ::CLHEP::c_light;
461
double total_delta_t = 0.;
462
for (size_t index = 0; index < 3; ++index) {
463
total_delta_t += delta_t[index];
465
const double residual_delta_t = actual_delta_t - total_delta_t;
466
std::cout << "DEBUG FindEnergy: Actual dt = " << actual_delta_t
467
<< "\tdt = " << total_delta_t << std::endl;
468
std::cout << "DEBUG FindEnergy: Min E = " << E_min
469
<< "\tMax E = " << E_max << std::endl;
470
if (energy_too_low || (total_delta_t != total_delta_t)
471
|| (residual_delta_t < 0.)) {
475
energy_too_low = false;
476
} else if (residual_delta_t > 0.) {
477
// energy is too high
481
std::cout << "DEBUG FindEnergy: E[" << iteration << "] = " << E << std::endl;
487
/** TOFSlabEnergyLoss
489
double MapCppGlobalRawTracks::TOFSlabEnergyLoss(const double beta,
490
const double mass) const {
491
const size_t slab_length = 25; // mm steps
492
double E_0 = mass / ::sqrt(1 - beta*beta);
493
double beta_i = beta;
495
for (size_t step = 0; step < slab_length; ++step) {
496
const double dEdx = TOFMeanStoppingPower(beta_i, mass);
502
beta_i = ::sqrt(1 - ::pow(mass / E_i, 2));
503
std::cout << "DEBUG TOFSlabEnergyLoss: dE/dx[" << step << "]: " << dEdx
505
<< "\t mass: " << mass << "\tbeta_i: " << beta_i << std::endl;
510
/** TOFMeanStoppingPower
512
double MapCppGlobalRawTracks::TOFMeanStoppingPower(const double beta,
513
const double mass) const {
514
const double gamma = 1. / ::sqrt(1-beta*beta);
515
const double beta2 = beta*beta;
516
const double gamma2 = gamma*gamma;
517
// values are for polystyrene
518
const double density = 1.060; // g cm^-3
519
const double me = 0.510998918; // electron mass (MeV/c^2)
520
const double K = 0.307075; // MeV g^-2 cm^2 mol
521
const double Tmax = 2 * me * beta2 * gamma2
522
/ (1 + 2 * gamma * me/mass + ::pow(me/mass, 2));
523
const double Z_over_A = 0.53768;
524
const double I = 6.87e-5; // MeV (68.7 eV)
526
std::cout << "DEBUG TOFMeanStoppingPower: beta: " << beta << "\tgamma: "
527
<< gamma << "\tTmax: " << Tmax << std::endl;
528
double density_effect = 0.;
529
const double x0 = 0.1647;
530
const double x1 = 2.5031;
531
const double Cbar = 3.2999;
532
const double a = 0.1645;
533
const double k = 3.2224;
534
const double x = beta * gamma;
536
density_effect = 2 * ::log(10) * x - Cbar;
537
} else if ((x0 <= x) && (x < x1)) {
538
density_effect = 2 * ::log(10) * x - Cbar + a * ::pow(x1-x, k);
541
const double log_operand = 2 * me * beta2 * gamma2 * Tmax / (I*I);
542
std::cout << "DEBUG TOFMeanStoppingPower: ln operand: " << log_operand
545
// Bethe-Bloch equation times the density of polystyrene in MeV / mm
546
const double stopping_power
547
= K * Z_over_A / beta2
548
* ( 0.5 * ::log(log_operand) - beta2 - density_effect / 2)
550
std::cout << "DEBUG TOFMeanStoppingPower: Stopping Power: " << stopping_power
553
return stopping_power;
557
* Assumes tracks have already been added to each recon_event.
559
void MapCppGlobalRawTracks::LoadSciFiTracks(
560
MAUS::ReconEvent const * const recon_event,
561
GlobalDS::TrackPArray & tracks) {
562
if (tracks.size() == 0) {
563
// Add a new track if none previously existed
564
GlobalDS::Track * track = new GlobalDS::Track();
565
tracks.push_back(track);
567
const size_t num_global_tracks = tracks.size();
569
const int particle_event_number = recon_event->GetPartEventNumber();
571
SciFiEvent const * const scifi_event = recon_event->GetSciFiEvent();
572
SciFiTrackPArray scifi_tracks = scifi_event->scifitracks();
574
// Duplicate each global track for each additional SciFi track. When done
575
// the order will be A1,B1,C1,...,A2,B2,C2,...; where A,B,C represent unique
576
// global tracks and 1,2,3 index the SciFi track.
577
const size_t num_scifi_tracks = scifi_tracks.size();
578
tracks.resize(num_global_tracks * num_scifi_tracks);
579
GlobalDS::TrackPArray::iterator track;
580
for (size_t copy_index = 1; copy_index < num_scifi_tracks; ++copy_index) {
581
track = tracks.begin();
582
for (size_t trk_index = 0; trk_index < num_global_tracks; ++trk_index) {
583
tracks.push_back(new GlobalDS::Track(**track));
587
track = tracks.begin();
589
std::cout << "DEBUG MapCppGlobalRawTracks::LoadSciFiTrack: " << std::endl
590
<< "\t# Previous Global Tracks: " << num_global_tracks << std::endl
591
<< "\t# SciFi Tracks: " << num_scifi_tracks << std::endl
592
<< "\tNew # Global Tracks: " << tracks.size() << std::endl;
594
SciFiTrackPArray::const_iterator scifi_track = scifi_tracks.begin();
595
for (; scifi_track != scifi_tracks.end(); ++scifi_track) {
596
SciFiTrackPointPArray scifi_track_points
597
= (*scifi_track)->scifitrackpoints();
598
SciFiTrackPointPArray::const_iterator scifi_track_point
599
= scifi_track_points.begin();
600
GlobalDS::TrackPArray::iterator track_begin = track;
601
GlobalDS::TrackPArray::iterator track_end = track + num_global_tracks;
602
for (; scifi_track_point != scifi_track_points.end(); ++scifi_track_point) {
603
const int tracker = (*scifi_track_point)->tracker();
604
const int station = (*scifi_track_point)->station();
605
std::cout << "DEBUG MapCppGlobalRawTracks::LoadSciFiTrack: " << std::endl
606
<< "\tTracker: " << tracker << "\tStation: " << station
607
<< "\tPz: " << (*scifi_track_point)->pz() << std::endl;
608
const GlobalDS::DetectorPoint detector_id = GlobalDS::DetectorPoint(
609
GlobalDS::kTracker0 + 6 * tracker + station);
610
MAUS::recon::global::DetectorMap::const_iterator detector_mapping
611
= detectors_.find(detector_id);
612
if (detector_mapping == detectors_.end()) {
613
std::stringstream message;
614
message << "Couldn't find configuration for detector Tracker "
615
<< tracker << " Station " << station
616
<< "(id=" << detector_id << ")";
617
throw(Exception(Exception::nonRecoverable,
619
"MapCppGlobalRawTracks::LoadSciFiTrack()"));
621
const Detector& detector = detector_mapping->second;
622
GlobalDS::TrackPoint * track_point = new GlobalDS::TrackPoint();
623
track_point->set_particle_event(particle_event_number);
624
PopulateSciFiTrackPoint(detector, scifi_track_point, track_point);
626
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
627
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
628
= simulator->GetReferenceParticle();
629
const GlobalDS::PID particle_id = GlobalDS::PID(reference_pgparticle.pid);
631
// For each SciFi track, add it's track points to the associated set of
632
// global track copies
633
for (; track != track_end; ++track) {
634
if (track == tracks.end()) {
635
throw(Exception(Exception::nonRecoverable,
636
"insufficient number of global tracks to accomodate"
637
"the number of SciFi tracks",
638
"MapCppGlobalRawTracks::LoadSciFiTracks()"));
640
// FIXME(Lane) this will have to change once PID functionality exists
641
(*track)->set_pid(particle_id);
643
std::cout << "DEBUG MapCppGlobalRawTracks::LoadSciFiTracks(): "
644
<< "Adding track point...";
645
(*track)->AddTrackPoint(track_point);
646
std::cout << "DONE!" << std::endl;
647
} // end for loop over global tracks
648
track = track_begin; // reset for next SciFi track point
649
} // end for loop over SciFi track points
650
track += num_global_tracks; // go to the next set of global tracks
651
} // end for loop over SciFi tracks
654
/** PopulateSciFiTrackPoint
656
void MapCppGlobalRawTracks::PopulateSciFiTrackPoint(
657
const MAUS::recon::global::Detector & detector,
658
const SciFiTrackPointPArray::const_iterator & scifi_track_point,
659
GlobalDS::TrackPoint * track_point) {
660
GlobalDS::SpacePoint * space_point = new GlobalDS::SpacePoint();
662
space_point->set_detector(detector.id());
664
// FIXME(Lane) not sure what to put here
665
space_point->set_geometry_path("");
667
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
668
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
669
= simulator->GetReferenceParticle();
671
DataStructureHelper helper = DataStructureHelper::GetInstance();
672
const double x = (detector.id() <= GlobalDS::kTracker0_5)?
673
-(*scifi_track_point)->x():(*scifi_track_point)->x();
674
const double y = (*scifi_track_point)->y();
675
const double z = helper.GetDetectorZPosition(detector.id());
676
std::cout << "DEBUG MapCppGlobalRawTracks::PopulateSciFiTrackPoint: "
677
<< "z position of detector with ID " << detector.id() << " is "
679
// No timestamp in the SciFiEvent/Track/TrackPoint data structure
680
const double time = 0.;
681
TLorentzVector four_position(x, y, z, time);
682
space_point->set_position(four_position);
684
static_cast<GlobalDS::BasePoint*>(track_point)->operator=(*space_point);
685
track_point->set_space_point(space_point);
687
track_point->set_mapper_name(kClassname);
689
const double Px = (*scifi_track_point)->px();
690
const double Py = (*scifi_track_point)->py();
691
const double Pz = (*scifi_track_point)->pz();
692
const double momentum = ::sqrt(Px*Px + Py*Py + Pz*Pz);
693
const GlobalDS::PID particle_id = GlobalDS::PID(reference_pgparticle.pid);
694
const double beta = Beta(particle_id, momentum);
695
const double energy = momentum / beta;
696
TLorentzVector four_momentum(Px, Py, Pz, energy);
697
track_point->set_momentum(four_momentum);
698
std::cout << "DEBUG MapCppGlobalRawTracks::PopulateSciFiTrackPoint(): "
699
<< "\tSciFi Space Point:" << std::endl
700
<< std::setprecision(10)
701
<< "\tTime: " << time << "\tEnergy: " << energy << std::endl
702
<< "\t4-Position: (" << four_position.X() << ", "
703
<< four_position.Y() << ", " << four_position.Z() << ", "
704
<< four_position.T() << ")" << std::endl
705
<< "\tMomentum: (" << four_momentum.Px() << ", "
706
<< four_momentum.Py() << ", " << four_momentum.Pz() << ", "
707
<< four_momentum.Pt() << ")" << std::endl;
708
CovarianceMatrix covariances = detector.uncertainties();
709
TLorentzVector position_errors(::sqrt(covariances(3, 3)),
710
::sqrt(covariances(5, 5)),
712
::sqrt(covariances(1, 1)));
713
track_point->set_position_error(position_errors);
715
TLorentzVector momentum_errors(::sqrt(covariances(4, 4)),
716
::sqrt(covariances(6, 6)),
718
::sqrt(covariances(2, 2)));
719
track_point->set_momentum_error(momentum_errors);
721
// track_point->set_charge((*scifi_track_point)->get_npe());
722
std::cout << "DEBUG MapCppGlobalRawTracks::PopulateSciFiTrackPoint(): END!"
726
/* Take an educated guess at the particle ID based on the axial velocity
727
* of the particle (beta)
729
GlobalDS::PID MapCppGlobalRawTracks::IdentifyParticle(const double beta) {
730
MAUSGeant4Manager * const simulator = MAUSGeant4Manager::GetInstance();
731
MAUSPrimaryGeneratorAction::PGParticle reference_pgparticle
732
= simulator->GetReferenceParticle();
733
const double reference_momentum = reference_pgparticle.pz; // good enough
734
const double beta_pi = Beta(GlobalDS::kPiPlus, reference_momentum);
735
const double beta_mu = Beta(GlobalDS::kMuPlus, reference_momentum);
736
const double beta_e = Beta(GlobalDS::kEPlus, reference_momentum);
737
std::cout << "DEBUG MapCppGlobalRawTracks::IdentifyParticle(): " << std::endl
738
<< "\tbeta = " << beta << "\tbeta_pi = " << beta_pi << std::endl
739
<< "\tbeta_mu = " << beta_mu << "\tbeta_e = " << beta_e << std::endl;
741
const double beta_fit[3] = {
742
std::abs(1. - beta / beta_pi),
743
std::abs(1. - beta / beta_mu),
744
std::abs(1. - beta / beta_e)
746
std::cout << "DEBUG MapCppGlobalRawTracks::IdentifyParticle(): "
747
<< "beta fits {" << std::endl << std::setprecision(4)
748
<< "\t" << 1. - beta / beta_pi
749
<< "\t" << beta_fit[1]
750
<< "\t" << beta_fit[2]
751
<< std::endl << "}" << std::endl;
752
size_t min_index = 0;
753
for (size_t index = 1; index < 3; ++index) {
754
if (beta_fit[index] < beta_fit[min_index]) {
759
GlobalDS::PID particle_id = GlobalDS::kNoPID;
761
case 0: particle_id = GlobalDS::PID(GlobalDS::kPiPlus * beam_polarity_);
763
case 1: particle_id = GlobalDS::PID(GlobalDS::kMuPlus * beam_polarity_);
765
case 2: particle_id = GlobalDS::PID(GlobalDS::kEPlus * beam_polarity_);
772
double MapCppGlobalRawTracks::Beta(GlobalDS::PID pid, const double momentum) {
773
const double mass = Particle::GetInstance().GetMass(pid);
774
const double one_over_beta_squared
775
= 1+ (mass*mass) / (momentum*momentum);
776
return 1. / std::sqrt(one_over_beta_squared);
779
const std::string MapCppGlobalRawTracks::kClassname
780
= "MapCppGlobalRawTracks";
added
removed
src/map/MapCppGlobalRawTracks/MapCppGlobalRawTracks.cc
