18
#include "Recon/Global/TrackMatching.hh"
20
18
#include <algorithm>
22
#include "Interface/Squeak.hh"
23
#include "Utils/Exception.hh"
20
#include "CLHEP/Units/PhysicalConstants.h"
22
#include "src/common_cpp/Recon/Global/GlobalTools.hh"
23
#include "src/common_cpp/Recon/Global/Particle.hh"
24
#include "src/common_cpp/Utils/Globals.hh"
25
#include "src/common_cpp/Utils/Exception.hh"
27
#include "src/legacy/BeamTools/BTFieldConstructor.hh"
28
#include "src/legacy/Interface/Squeak.hh"
30
#include "src/common_cpp/Recon/Global/TrackMatching.hh"
29
void TrackMatching::FormTracks(MAUS::GlobalEvent* global_event, std::string mapper_name) {
32
throw(Exception(Exception::recoverable,
33
"Trying to import an empty global event.",
34
"MapCppGlobalTrackMatching::TrackMatching"));
37
MAUS::DataStructure::Global::TrackPArray *ImportedTracks = global_event->get_tracks();
38
MAUS::DataStructure::Global::TrackPArray::iterator ImportedTrackIterator;
39
MAUS::DataStructure::Global::Track* ImportedSciFiTrack =
40
new MAUS::DataStructure::Global::Track();
41
for (ImportedTrackIterator = ImportedTracks->begin();
42
ImportedTrackIterator != ImportedTracks->end();
43
++ImportedTrackIterator) {
44
MAUS::DataStructure::Global::Track* ImportedTrack =
45
(*ImportedTrackIterator);
46
if (ImportedTrack->HasDetector(MAUS::DataStructure::Global::kTracker0) ||
47
ImportedTrack->HasDetector(MAUS::DataStructure::Global::kTracker1)) {
48
std::vector<const MAUS::DataStructure::Global::TrackPoint*>
49
tempSciFiTrackPointArray = ImportedTrack->GetTrackPoints();
50
for (unsigned int i = 0; i < tempSciFiTrackPointArray.size(); i++) {
51
MAUS::DataStructure::Global::TrackPoint* tempSciFiTrackPoint =
52
const_cast<MAUS::DataStructure::Global::TrackPoint*>
53
(tempSciFiTrackPointArray[i]);
54
tempSciFiTrackPoint->set_mapper_name(mapper_name);
55
ImportedSciFiTrack->AddTrackPoint(tempSciFiTrackPoint);
58
if (ImportedSciFiTrack->GetTrackPoints().size() == 0) {
59
delete ImportedSciFiTrack;
60
ImportedSciFiTrack = NULL;
64
std::vector<MAUS::DataStructure::Global::SpacePoint*>
65
*GlobalSpacePointArray = global_event->get_space_points();
66
MAUS::DataStructure::Global::TrackPArray TOFTrackArray;
67
MakeTOFTracks(global_event, GlobalSpacePointArray, TOFTrackArray);
68
MAUS::DataStructure::Global::Track* KLTrack =
69
new MAUS::DataStructure::Global::Track();
70
MakeKLTracks(global_event, GlobalSpacePointArray, KLTrack);
71
if (KLTrack->GetTrackPoints().size() == 0) {
76
// Adding global tracks for case where global event contains both SciFi and TOF tracks
77
// (And KL track if applicable)
78
if (ImportedSciFiTrack != NULL && !TOFTrackArray.empty()) {
79
for (unsigned int j = 0; j < TOFTrackArray.size(); j++) {
80
MAUS::DataStructure::Global::Track* GlobalTrack = TOFTrackArray[j]->Clone();
81
GlobalTrack->set_mapper_name(mapper_name);
82
std::vector<const MAUS::DataStructure::Global::TrackPoint*>
83
tempSciFiTrackPointArray = ImportedSciFiTrack->GetTrackPoints();
84
for (unsigned int k = 0; k < tempSciFiTrackPointArray.size(); k++) {
85
MAUS::DataStructure::Global::TrackPoint* tempSciFiTrackPoint =
86
const_cast<MAUS::DataStructure::Global::TrackPoint*>
87
(tempSciFiTrackPointArray[k]);
88
tempSciFiTrackPoint->set_mapper_name(mapper_name);
89
GlobalTrack->AddTrackPoint(tempSciFiTrackPoint);
91
if (KLTrack != NULL) {
92
std::vector<const MAUS::DataStructure::Global::TrackPoint*>
93
tempKLTrackPointArray = KLTrack->GetTrackPoints();
94
for (unsigned int l = 0; l < tempKLTrackPointArray.size(); l++) {
95
MAUS::DataStructure::Global::TrackPoint* tempKLTrackPoint =
96
const_cast<MAUS::DataStructure::Global::TrackPoint*>
97
(tempKLTrackPointArray[l]);
98
tempKLTrackPoint->set_mapper_name(mapper_name);
99
GlobalTrack->AddTrackPoint(tempKLTrackPoint);
102
global_event->add_track_recursive(GlobalTrack);
106
// Adding global tracks for case where global event contains only TOF tracks
107
// (And KL track if applicable)
108
if (ImportedSciFiTrack == NULL && !TOFTrackArray.empty()) {
109
for (unsigned int i = 0; i < TOFTrackArray.size(); i++) {
110
MAUS::DataStructure::Global::Track* GlobalTrack = TOFTrackArray[i]->Clone();
111
GlobalTrack->set_mapper_name(mapper_name);
112
if (KLTrack != NULL) {
113
std::vector<const MAUS::DataStructure::Global::TrackPoint*>
114
tempKLTrackPointArray = KLTrack->GetTrackPoints();
115
for (unsigned int l = 0; l < tempKLTrackPointArray.size(); l++) {
116
MAUS::DataStructure::Global::TrackPoint* tempKLTrackPoint =
117
const_cast<MAUS::DataStructure::Global::TrackPoint*>
118
(tempKLTrackPointArray[l]);
119
tempKLTrackPoint->set_mapper_name(mapper_name);
120
GlobalTrack->AddTrackPoint(tempKLTrackPoint);
123
global_event->add_track_recursive(GlobalTrack);
127
// Adding global tracks for case where global event contains only SciFi tracks
128
// (And KL track if applicable)
129
if (ImportedSciFiTrack != NULL && TOFTrackArray.empty()) {
130
MAUS::DataStructure::Global::Track* GlobalTrack = ImportedSciFiTrack->Clone();
131
GlobalTrack->set_mapper_name(mapper_name);
132
if (KLTrack != NULL) {
133
std::vector<const MAUS::DataStructure::Global::TrackPoint*>
134
tempKLTrackPointArray = KLTrack->GetTrackPoints();
135
for (unsigned int l = 0; l < tempKLTrackPointArray.size(); l++) {
136
MAUS::DataStructure::Global::TrackPoint* tempKLTrackPoint =
137
const_cast<MAUS::DataStructure::Global::TrackPoint*>
138
(tempKLTrackPointArray[l]);
139
tempKLTrackPoint->set_mapper_name(mapper_name);
140
GlobalTrack->AddTrackPoint(tempKLTrackPoint);
143
global_event->add_track_recursive(GlobalTrack);
146
// Adding global tracks for case where global event contains only a KL track
147
if (ImportedSciFiTrack == NULL && TOFTrackArray.empty() && KLTrack != NULL) {
148
MAUS::DataStructure::Global::Track* GlobalTrack = KLTrack->Clone();
149
GlobalTrack->set_mapper_name(mapper_name);
150
global_event->add_track_recursive(GlobalTrack);
154
void TrackMatching::MakeTOFTracks(
155
MAUS::GlobalEvent* global_event,
156
std::vector<MAUS::DataStructure::Global::SpacePoint*>
157
*GlobalSpacePointArray,
158
MAUS::DataStructure::Global::TrackPArray& TOFTrackArray) {
160
std::string local_mapper_name = "GlobalTOFTrack";
162
std::vector<MAUS::DataStructure::Global::TrackPoint*> TOF0tp;
163
std::vector<MAUS::DataStructure::Global::TrackPoint*> TOF1tp;
164
std::vector<MAUS::DataStructure::Global::TrackPoint*> TOF2tp;
166
std::vector<MAUS::DataStructure::Global::TrackPoint*> tempTOF0tp;
167
std::vector<MAUS::DataStructure::Global::TrackPoint*> tempTOF2tp;
168
double TOF01offset = 30;
169
double TOF12offset = 35;
170
double allowance = 8.0;
171
for (unsigned int i = 0; i < GlobalSpacePointArray->size(); ++i) {
172
MAUS::DataStructure::Global::SpacePoint* sp = GlobalSpacePointArray->at(i);
176
MAUS::DataStructure::Global::TrackPoint* tp0;
177
MAUS::DataStructure::Global::TrackPoint* tp1;
178
MAUS::DataStructure::Global::TrackPoint* tp2;
179
if (sp->get_detector() == MAUS::DataStructure::Global::kTOF0) {
180
tp0 = new MAUS::DataStructure::Global::TrackPoint(sp);
181
TOF0tp.push_back(tp0);
182
} else if (sp->get_detector() == MAUS::DataStructure::Global::kTOF1) {
183
tp1 = new MAUS::DataStructure::Global::TrackPoint(sp);
184
TOF1tp.push_back(tp1);
185
} else if (sp->get_detector() == MAUS::DataStructure::Global::kTOF2) {
186
tp2 = new MAUS::DataStructure::Global::TrackPoint(sp);
187
TOF2tp.push_back(tp2);
193
for (unsigned int i = 0; i < TOF1tp.size(); ++i) {
194
MAUS::DataStructure::Global::Track* TOFtrack =
195
new MAUS::DataStructure::Global::Track();
196
TOFtrack->set_mapper_name(local_mapper_name);
197
for (unsigned int j = 0; j < TOF0tp.size(); ++j) {
198
if ((TOF1tp[i]->get_position().T() - TOF0tp[j]->get_position().T()) <=
199
(TOF01offset + allowance) &&
200
(TOF1tp[i]->get_position().T() - TOF0tp[j]->get_position().T()) >=
201
(TOF01offset - allowance)) {
202
tempTOF0tp.push_back(TOF0tp[j]);
205
for (unsigned int k = 0; k < TOF2tp.size(); ++k) {
206
if ((TOF2tp[k]->get_position().T() - TOF1tp[i]->get_position().T()) <=
207
(TOF12offset + allowance) &&
208
(TOF2tp[k]->get_position().T() - TOF1tp[i]->get_position().T()) >=
209
(TOF12offset - allowance)) {
210
tempTOF2tp.push_back(TOF2tp[k]);
213
if (tempTOF0tp.size() < 2 && tempTOF2tp.size() < 2) {
214
TOF1tp[i]->set_mapper_name(local_mapper_name);
215
TOFtrack->AddTrackPoint(TOF1tp[i]);
216
global_event->add_track_point_recursive(TOF1tp[i]);
217
if (tempTOF0tp.size() == 1) {
218
tempTOF0tp[0]->set_mapper_name(local_mapper_name);
219
TOFtrack->AddTrackPoint(tempTOF0tp[0]);
220
global_event->add_track_point_recursive(tempTOF0tp[0]);
222
if (tempTOF2tp.size() == 1) {
223
tempTOF2tp[0]->set_mapper_name(local_mapper_name);
224
TOFtrack->AddTrackPoint(tempTOF2tp[0]);
225
global_event->add_track_point_recursive(tempTOF2tp[0]);
227
TOFTrackArray.push_back(TOFtrack);
229
Squeak::mout(Squeak::debug) << "Global event returned multiple possible"
230
<< " TOF0 and/or TOF2 space points that "
231
<< "could not be separated into tracks."
240
void TrackMatching::MakeKLTracks(
241
MAUS::GlobalEvent* global_event,
242
std::vector<MAUS::DataStructure::Global::SpacePoint*>
243
*GlobalSpacePointArray,
244
MAUS::DataStructure::Global::Track* KLTrack) {
247
std::string local_mapper_name = "GlobalKLTrack";
249
std::vector<MAUS::DataStructure::Global::TrackPoint*> KLtp;
251
for (unsigned int i = 0; i < GlobalSpacePointArray->size(); ++i) {
252
MAUS::DataStructure::Global::SpacePoint* sp = GlobalSpacePointArray->at(i);
256
MAUS::DataStructure::Global::TrackPoint* tp;
257
if (sp->get_detector() == MAUS::DataStructure::Global::kCalorimeter) {
258
tp = new MAUS::DataStructure::Global::TrackPoint(sp);
265
for (unsigned int i = 0; i < KLtp.size(); ++i) {
266
KLTrack->set_mapper_name(local_mapper_name);
267
KLtp[i]->set_mapper_name(local_mapper_name);
268
KLTrack->AddTrackPoint(KLtp[i]);
269
global_event->add_track_point_recursive(KLtp[i]);
36
TrackMatching::TrackMatching(GlobalEvent* global_event, std::string mapper_name,
37
std::string pid_hypothesis_string,
38
std::map<std::string, std::pair<double, double> > matching_tolerances,
39
double max_step_size, bool energy_loss) {
40
_global_event = global_event;
41
_mapper_name = mapper_name;
42
_pid_hypothesis_string = pid_hypothesis_string;
43
_matching_tolerances = matching_tolerances;
44
_max_step_size = max_step_size;
45
_energy_loss = energy_loss;
48
void TrackMatching::USTrack() {
49
// Get all Tracker 0 tracks
50
DataStructure::Global::TrackPArray *scifi_track_array =
51
GetDetectorTrackArray(DataStructure::Global::kTracker0);
52
// Get Spacepoints for TOF0/1 and convert them to TrackPoints
53
std::vector<DataStructure::Global::TrackPoint*> TOF0_tp =
54
GetDetectorTrackPoints(DataStructure::Global::kTOF0, _mapper_name);
55
std::vector<DataStructure::Global::TrackPoint*> TOF1_tp =
56
GetDetectorTrackPoints(DataStructure::Global::kTOF1, _mapper_name);
58
// Load the magnetic field for RK4 propagation
59
BTFieldConstructor* field = Globals::GetMCFieldConstructor();
60
// Iterate over all Tracker0 Tracks (typically 1)
61
for (auto scifi_track_iter = scifi_track_array->begin();
62
scifi_track_iter != scifi_track_array->end();
64
// Pull out the track so we're not making a mess with ressources
65
DataStructure::Global::Track* tracker0_track = (*scifi_track_iter);
66
// Extract four-position and momentum from first track point (i.e. most
68
DataStructure::Global::TrackPoint* first_tracker_tp =
69
GlobalTools::GetNearestZTrackPoint(tracker0_track, 0);
70
TLorentzVector position = first_tracker_tp->get_position();
71
TLorentzVector momentum = first_tracker_tp->get_momentum();
72
delete first_tracker_tp;
73
// Create the list of PIDs for which we want to create hypothesis tracks
74
int charge_hypothesis = tracker0_track->get_charge();
75
std::vector<DataStructure::Global::PID> pids = PIDHypotheses(
76
charge_hypothesis, _pid_hypothesis_string);
77
// Iterate over all possible PIDs and create an hypothesis track for each
78
for (size_t i = 0; i < pids.size(); i++) {
79
DataStructure::Global::Track* hypothesis_track =
80
new DataStructure::Global::Track();
81
hypothesis_track->set_mapper_name("MapCppGlobalTrackMatching-US");
82
hypothesis_track->set_pid(pids[i]);
83
MatchTrackPoint(position, momentum, TOF0_tp, pids[i], field, "TOF0",
85
MatchTrackPoint(position, momentum, TOF1_tp, pids[i], field, "TOF1",
88
// Now we fill the track with trackpoints from the tracker with energy
89
// calculated from p and m, trackpoints are cloned as we want everything
90
// in the hypothesis track to be "fresh"
91
double mass = Particle::GetInstance().GetMass(pids[i]);
92
AddTrackerTrackPoints(tracker0_track, "MapCppGlobalTrackMatching-US",
93
mass, hypothesis_track);
95
_global_event->add_track_recursive(hypothesis_track);
97
_global_event->add_track_recursive(tracker0_track);
101
void TrackMatching::DSTrack() {
102
// Get all Tracker 1 tracks
103
DataStructure::Global::TrackPArray *scifi_track_array =
104
GetDetectorTrackArray(DataStructure::Global::kTracker1);
105
// Get all EMR tracks
106
DataStructure::Global::TrackPArray *emr_track_array =
107
GetDetectorTrackArray(DataStructure::Global::kEMR);
108
// Get Spacepoints for TOF0/1 and convert them to TrackPoints
109
std::vector<DataStructure::Global::TrackPoint*> TOF2_tp =
110
GetDetectorTrackPoints(DataStructure::Global::kTOF2, _mapper_name);
111
std::vector<DataStructure::Global::TrackPoint*> KL_tp =
112
GetDetectorTrackPoints(DataStructure::Global::kCalorimeter, _mapper_name);
113
// Load the magnetic field for RK4 propagation
114
BTFieldConstructor* field = Globals::GetMCFieldConstructor();
115
// Iterate over all Tracker1 Tracks (typically 1)
116
for (auto scifi_track_iter = scifi_track_array->begin();
117
scifi_track_iter != scifi_track_array->end();
118
++scifi_track_iter) {
119
// Pull out the track so we're not making a mess with ressources
120
DataStructure::Global::Track* tracker1_track = (*scifi_track_iter);
121
// Extract four-position and momentum from last track point (i.e. most
123
DataStructure::Global::TrackPoint* last_tracker_tp =
124
GlobalTools::GetNearestZTrackPoint(tracker1_track, 100000);
125
TLorentzVector position = last_tracker_tp->get_position();
126
TLorentzVector momentum = last_tracker_tp->get_momentum();
127
// Create the list of PIDs for which we want to create hypothesis tracks
128
int charge_hypothesis = tracker1_track->get_charge();
129
std::vector<DataStructure::Global::PID> pids = PIDHypotheses(
130
charge_hypothesis, _pid_hypothesis_string);
131
// Iterate over all possible PIDs and create an hypothesis track for each
132
for (size_t i = 0; i < pids.size(); i++) {
133
DataStructure::Global::Track* hypothesis_track =
134
new DataStructure::Global::Track();
135
hypothesis_track->set_mapper_name("MapCppGlobalTrackMatching-DS");
136
hypothesis_track->set_pid(pids[i]);
138
// This time we add in the tracker trackpoints first
139
double mass = Particle::GetInstance().GetMass(pids[i]);
140
AddTrackerTrackPoints(tracker1_track, "MapCppGlobalTrackMatching-DS",
141
mass, hypothesis_track);
143
MatchTrackPoint(position, momentum, TOF2_tp, pids[i], field, "TOF2",
146
MatchTrackPoint(position, momentum, KL_tp, pids[i], field, "KL",
149
MatchEMRTrack(position, momentum, emr_track_array, pids[i], field,
151
_global_event->add_track_recursive(hypothesis_track);
156
void TrackMatching::throughTrack() {
157
// import tracks already created by tracker recon
158
DataStructure::Global::TrackPArray* global_tracks =
159
_global_event->get_tracks();
160
// Typically used for no fields so we can't discriminate by charge hypothesis
161
std::vector<DataStructure::Global::PID> pids = PIDHypotheses(0,
162
_pid_hypothesis_string);
163
// Iterate over all PIDs
164
for (size_t i = 0; i < pids.size(); i++) {
165
DataStructure::Global::TrackPArray* us_tracks =
166
new DataStructure::Global::TrackPArray();
167
DataStructure::Global::TrackPArray* ds_tracks =
168
new DataStructure::Global::TrackPArray();
170
// Loop over all global tracks and sort into US and DS tracks by mapper name
171
// provided they have the correct PID and hits in TOF1/2 to match by dT
172
USDSTracks(global_tracks, pids[i], us_tracks, ds_tracks);
174
// Do we have both US and DS tracks in the event?
175
if ((us_tracks->size() > 0) and
176
(ds_tracks->size() > 0)) {
177
// Iterate over all possible combinations of US and DS tracks
178
for (auto us_track_iter = us_tracks->begin();
179
us_track_iter != us_tracks->end();
181
for (auto ds_track_iter = ds_tracks->begin();
182
ds_track_iter != ds_tracks->end();
184
// Going to assume that if several TrackPoints for the same TOF are in
185
// the track that they have been checked to be sensible (TODO above)
186
DataStructure::Global::TrackPointCPArray us_trackpoints =
187
(*us_track_iter)->GetTrackPoints();
188
DataStructure::Global::TrackPointCPArray ds_trackpoints =
189
(*ds_track_iter)->GetTrackPoints();
190
double emr_range_primary = (*ds_track_iter)->get_emr_range_primary();
191
if ((us_trackpoints.size() > 0) and
192
(ds_trackpoints.size() > 0)) {
193
MatchUSDS(us_trackpoints, ds_trackpoints, pids[i],
202
DataStructure::Global::TrackPArray* TrackMatching::GetDetectorTrackArray(
203
DataStructure::Global::DetectorPoint detector) {
204
DataStructure::Global::TrackPArray *imported_tracks =
205
_global_event->get_tracks();
206
DataStructure::Global::TrackPArray *track_array =
207
new DataStructure::Global::TrackPArray();
208
for (auto imported_track_iter = imported_tracks->begin();
209
imported_track_iter != imported_tracks->end();
210
++imported_track_iter) {
211
DataStructure::Global::Track* imported_track = (*imported_track_iter);
212
// Check that track is from correct detector
213
if (imported_track->HasDetector(detector)) {
214
// Exclude EMR secondary tracks (tracker tracks also have this at 0 by
216
if (imported_track->get_emr_range_secondary() < 0.001) {
217
track_array->push_back(imported_track);
224
std::vector<DataStructure::Global::TrackPoint*>
225
TrackMatching::GetDetectorTrackPoints(
226
DataStructure::Global::DetectorPoint detector, std::string mapper_name) {
227
std::vector<DataStructure::Global::TrackPoint*> track_points;
228
std::vector<DataStructure::Global::SpacePoint*> *global_spacepoint_array =
229
_global_event->get_space_points();
230
for (size_t i = 0; i < global_spacepoint_array->size(); i++) {
231
DataStructure::Global::SpacePoint* space_point =
232
global_spacepoint_array->at(i);
236
if (space_point->get_detector() == detector) {
237
DataStructure::Global::TrackPoint* track_point =
238
new DataStructure::Global::TrackPoint(space_point);
239
std::string blsection_mapper_name = mapper_name;
240
if ((detector == DataStructure::Global::kTOF0) or
241
(detector == DataStructure::Global::kTOF1)) {
242
blsection_mapper_name.append("-US");
244
blsection_mapper_name.append("-DS");
246
track_point->set_mapper_name(blsection_mapper_name);
247
track_points.push_back(track_point);
253
std::vector<DataStructure::Global::PID> TrackMatching::PIDHypotheses(
254
int charge_hypothesis, std::string pid_hypothesis_string) {
255
std::vector<DataStructure::Global::PID> pids;
256
// First check if a specific PID has been forced by the configuration
257
if (pid_hypothesis_string == "kEPlus") {
258
pids.push_back(DataStructure::Global::kEPlus);
259
} else if (pid_hypothesis_string == "kEMinus") {
260
pids.push_back(DataStructure::Global::kEMinus);
261
} else if (pid_hypothesis_string == "kMuPlus") {
262
pids.push_back(DataStructure::Global::kMuPlus);
263
} else if (pid_hypothesis_string == "kMuMinus") {
264
pids.push_back(DataStructure::Global::kMuMinus);
265
} else if (pid_hypothesis_string == "kPiPlus") {
266
pids.push_back(DataStructure::Global::kPiPlus);
267
} else if (pid_hypothesis_string == "kPiMinus") {
268
pids.push_back(DataStructure::Global::kPiMinus);
270
// No specific PID forced, hence we assign by charge hypothesis
271
if (charge_hypothesis != -1) {
272
pids.push_back(DataStructure::Global::kEPlus);
273
pids.push_back(DataStructure::Global::kMuPlus);
274
pids.push_back(DataStructure::Global::kPiPlus);
276
if (charge_hypothesis != 1) {
277
pids.push_back(DataStructure::Global::kEMinus);
278
pids.push_back(DataStructure::Global::kMuMinus);
279
pids.push_back(DataStructure::Global::kPiMinus);
285
void TrackMatching::MatchTrackPoint(
286
const TLorentzVector &position, const TLorentzVector &momentum,
287
const std::vector<DataStructure::Global::TrackPoint*> &trackpoints,
288
DataStructure::Global::PID pid, BTFieldConstructor* field,
289
std::string detector_name, DataStructure::Global::Track* hypothesis_track) {
290
double mass = Particle::GetInstance().GetMass(pid);
291
double energy = ::sqrt(momentum.Rho()*momentum.Rho() + mass*mass);
292
double x_in[] = {0., position.X(), position.Y(), position.Z(),
293
energy, momentum.X(), momentum.Y(), momentum.Z()};
294
if (trackpoints.size() > 0) {
295
double target_z = trackpoints.at(0)->get_position().Z();
297
GlobalTools::propagate(x_in, target_z, field, _max_step_size, pid,
299
for (size_t i = 0; i < trackpoints.size(); i++) {
300
if (GlobalTools::approx(x_in[1], trackpoints.at(i)->get_position().X(),
301
_matching_tolerances.at(detector_name).first) and
302
GlobalTools::approx(x_in[2], trackpoints.at(i)->get_position().Y(),
303
_matching_tolerances.at(detector_name).second)) {
304
hypothesis_track->AddTrackPoint(trackpoints.at(i));
305
Squeak::mout(Squeak::debug) << "TrackMatching: "
306
<< detector_name << " Match" << std::endl;
308
Squeak::mout(Squeak::debug) << "TrackMatching: " << detector_name
309
<< " Mismatch, dx, dy are:\n"
310
<< x_in[1] - trackpoints.at(i)->get_position().X() << " "
311
<< x_in[2] - trackpoints.at(i)->get_position().Y() << std::endl;
314
} catch (Exception exc) {
315
Squeak::mout(Squeak::error) << exc.what() << std::endl;
320
void TrackMatching::MatchTOF0(
321
const TLorentzVector &position, const TLorentzVector &momentum,
322
const std::vector<DataStructure::Global::TrackPoint*> &trackpoints,
323
DataStructure::Global::PID pid,
324
DataStructure::Global::Track* hypothesis_track) {
325
double mass = Particle::GetInstance().GetMass(pid);
326
double energy = ::sqrt(momentum.Rho()*momentum.Rho() + mass*mass);
327
if (trackpoints.size() > 0) {
328
double z_distance = position.Z() - trackpoints.at(0)->get_position().Z();
329
double velocity = (momentum.Z() / energy) * CLHEP::c_light;
330
// Change later to be set by datacards
331
double deltaTMin = (z_distance/velocity) - 2.0;
332
double deltaTMax = (z_distance/velocity) + 2.0;
333
for (size_t i = 0; i < trackpoints.size(); i++) {
334
double deltaT = position.T() - trackpoints.at(i)->get_position().T();
335
if (deltaT > deltaTMin and deltaT < deltaTMax) {
336
hypothesis_track->AddTrackPoint(trackpoints.at(i));
337
Squeak::mout(Squeak::debug) << "TrackMatching: TOF0 Match"
340
Squeak::mout(Squeak::debug) << "TrackMatching: TOF0 Mismatch, "
341
<< "dT is " << deltaT << " when expected between " << deltaTMin
342
<< " and " << deltaTMax << std::endl;
348
void TrackMatching::MatchEMRTrack(
349
const TLorentzVector &position, const TLorentzVector &momentum,
350
DataStructure::Global::TrackPArray* emr_track_array,
351
DataStructure::Global::PID pid, BTFieldConstructor* field,
352
DataStructure::Global::Track* hypothesis_track) {
353
double mass = Particle::GetInstance().GetMass(pid);
354
double energy = ::sqrt(momentum.Rho()*momentum.Rho() + mass*mass);
355
// Here we need to iterate over the EMR tracks first, as they might have
356
// different starting z
357
for (auto emr_track_iter = emr_track_array->begin();
358
emr_track_iter != emr_track_array->end();
360
std::vector<const DataStructure::Global::TrackPoint*> emr_trackpoints =
361
(*emr_track_iter)->GetTrackPoints();
362
double x_in[] = {0., position.X(), position.Y(), position.Z(),
363
energy, momentum.X(), momentum.Y(), momentum.Z()};
364
TLorentzVector first_hit_pos = emr_trackpoints[0]->get_position();
365
TLorentzVector first_hit_pos_err = emr_trackpoints[0]->get_position_error();
366
double target_z = first_hit_pos.Z();
368
GlobalTools::propagate(x_in, target_z, field, _max_step_size, pid,
370
if (GlobalTools::approx(x_in[1], first_hit_pos.X(),
371
first_hit_pos_err.X()*::sqrt(12)) and
372
GlobalTools::approx(x_in[2], first_hit_pos.Y(),
373
first_hit_pos_err.Y()*::sqrt(12))) {
374
Squeak::mout(Squeak::debug) << "TrackMatching: EMR Match" << std::endl;
375
hypothesis_track->set_emr_range_primary(
376
(*emr_track_iter)->get_emr_range_primary());
377
for (size_t i = 0; i < emr_trackpoints.size(); i++) {
378
DataStructure::Global::TrackPoint* emr_tp =
379
emr_trackpoints[i]->Clone();
380
emr_tp->set_mapper_name("MapCppGlobalTrackMatching-DS");
381
TLorentzVector momentum = emr_tp->get_momentum();
382
double energy = ::sqrt(momentum.Rho()*momentum.Rho() + mass*mass);
383
momentum.SetE(energy);
384
emr_tp->set_momentum(momentum);
385
hypothesis_track->AddTrackPoint(emr_tp);
388
Squeak::mout(Squeak::debug)
389
<< "TrackMatching: EMR Mismatch, dx, dy are:\n"
390
<< x_in[1] - first_hit_pos.X() << " "
391
<< x_in[2] - first_hit_pos.Y() << std::endl;
393
} catch (Exception exc) {
394
Squeak::mout(Squeak::error) << exc.what() << std::endl;
399
void TrackMatching::AddTrackerTrackPoints(
400
DataStructure::Global::Track* tracker_track, std::string mapper_name,
401
double mass, DataStructure::Global::Track* hypothesis_track) {
402
std::vector<const DataStructure::Global::TrackPoint*> tracker_trackpoints =
403
tracker_track->GetTrackPoints();
404
std::sort(tracker_trackpoints.begin(), tracker_trackpoints.end(),
405
GlobalTools::TrackPointSort);
406
for (size_t i = 0; i < tracker_trackpoints.size(); i++) {
407
DataStructure::Global::TrackPoint* tracker_tp =
408
tracker_trackpoints[i]->Clone();
409
tracker_tp->set_mapper_name(mapper_name);
410
TLorentzVector momentum = tracker_tp->get_momentum();
411
double energy = ::sqrt(momentum.Rho()*momentum.Rho() + mass*mass);
412
momentum.SetE(energy);
413
tracker_tp->set_momentum(momentum);
414
hypothesis_track->AddTrackPoint(tracker_tp);
418
void TrackMatching::USDSTracks(
419
DataStructure::Global::TrackPArray* global_tracks,
420
DataStructure::Global::PID pid,
421
DataStructure::Global::TrackPArray* us_tracks,
422
DataStructure::Global::TrackPArray* ds_tracks) {
423
for (auto global_track_iter = global_tracks->begin();
424
global_track_iter != global_tracks->end();
425
++global_track_iter) {
426
if (((*global_track_iter)->get_mapper_name() ==
427
"MapCppGlobalTrackMatching-US") and
428
((*global_track_iter)->HasDetector(DataStructure::Global::kTOF1)) and
429
((*global_track_iter)->get_pid() == pid)) {
430
us_tracks->push_back(*global_track_iter);
431
} else if (((*global_track_iter)->get_mapper_name() ==
432
"MapCppGlobalTrackMatching-DS") and
433
((*global_track_iter)->HasDetector(
434
DataStructure::Global::kTOF2)) and
435
((*global_track_iter)->get_pid() == pid)) {
436
ds_tracks->push_back(*global_track_iter);
441
void TrackMatching::MatchUSDS(
442
DataStructure::Global::TrackPointCPArray us_trackpoints,
443
DataStructure::Global::TrackPointCPArray ds_trackpoints,
444
DataStructure::Global::PID pid, double emr_range_primary) {
445
double TOF1_time = TOFTimeFromTrackPoints(us_trackpoints,
446
DataStructure::Global::kTOF1);
447
double TOF2_time = TOFTimeFromTrackPoints(ds_trackpoints,
448
DataStructure::Global::kTOF2);
449
double TOFdT = TOF2_time - TOF1_time;
450
if ((TOFdT > _matching_tolerances.at("TOF12dT").first) and
451
(TOFdT < _matching_tolerances.at("TOF12dT").second)) {
452
DataStructure::Global::Track* through_track =
453
new DataStructure::Global::Track();
454
through_track->set_mapper_name("MapCppGlobalTrackMatching-Through");
455
through_track->set_pid(pid);
456
Squeak::mout(Squeak::debug) << "TrackMatching: US & DS Matched"
458
for (auto trackpoint_iter = us_trackpoints.begin();
459
trackpoint_iter != us_trackpoints.end();
461
through_track->AddTrackPoint(
462
const_cast<DataStructure::Global::TrackPoint*>(*trackpoint_iter));
464
for (auto trackpoint_iter = ds_trackpoints.begin();
465
trackpoint_iter != ds_trackpoints.end();
467
through_track->AddTrackPoint(
468
const_cast<DataStructure::Global::TrackPoint*>(*trackpoint_iter));
470
through_track->set_emr_range_primary(emr_range_primary);
471
_global_event->add_track(through_track);
475
double TrackMatching::TOFTimeFromTrackPoints(
476
DataStructure::Global::TrackPointCPArray trackpoints,
477
DataStructure::Global::DetectorPoint detector) {
478
double TOF_time = 0.0;
479
for (auto trackpoint_iter = trackpoints.begin();
480
trackpoint_iter != trackpoints.end();
482
if ((*trackpoint_iter)->get_detector() == detector) {
483
TOF_time = (*trackpoint_iter)->get_position().T();
273
489
} // ~namespace global
274
490
} // ~namespace recon
275
491
} // ~namespace MAUS
added
removed
src/common_cpp/Recon/Global/TrackMatching.cc
