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 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.
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/>.
24
#include "src/common_cpp/Utils/JsonWrapper.hh"
25
#include "src/common_cpp/Utils/CppErrorHandler.hh"
26
#include "src/common_cpp/Utils/Exception.hh"
27
#include "src/legacy/Interface/dataCards.hh"
28
#include "src/common_cpp/API/PyWrapMapBase.hh"
29
#include "src/map/MapCppCuts/MapCppCuts.hh"
33
std::string class_docstring =
34
std::string("MapCppCuts is the mapper for the MAUS Cuts class.\n");
36
std::string birth_docstring =
37
std::string("Checks if the right configuration is passed to the processor.\n");
39
std::string process_docstring =
40
std::string("Set up event(s) with some cut_event data and\n")+
41
std::string("check if mapper sets the correct cut values.\n");
43
std::string death_docstring =
44
std::string("Does nothing.\n");
46
PyMODINIT_FUNC init_MapCppCuts(void) {
47
PyWrapMapBase<MapCppCuts>::PyWrapMapBaseModInit
48
("MapCppCuts", class_docstring, birth_docstring, process_docstring, death_docstring);
51
MapCppCuts::MapCppCuts()
52
: MapBase<Data>("MapCppCuts") {
55
MapCppCuts::~MapCppCuts() {
58
void MapCppCuts::_birth(const std::string& argJsonConfigDocument) {
59
// Called at the beginning of each run.
60
// Check if the JSON document can be parsed, else return error
62
_configJSON = JsonWrapper::StringToJson(argJsonConfigDocument);
65
JsonWrapper::GetProperty(_configJSON, "cuts_parameters", JsonWrapper::objectValue);
68
JsonWrapper::GetProperty(_set_Cut_params, "min_tof", JsonWrapper::realValue).asDouble();
71
JsonWrapper::GetProperty(_set_Cut_params, "max_tof", JsonWrapper::realValue).asDouble();
74
JsonWrapper::GetProperty(_set_Cut_params, "min_mom_US", JsonWrapper::realValue).asDouble();
77
JsonWrapper::GetProperty(_set_Cut_params, "max_mom_US", JsonWrapper::realValue).asDouble();
80
JsonWrapper::GetProperty(_set_Cut_params, "min_mom_DS", JsonWrapper::realValue).asDouble();
83
JsonWrapper::GetProperty(_set_Cut_params, "max_mom_DS", JsonWrapper::realValue).asDouble();
86
JsonWrapper::GetProperty(_set_Cut_params, "min_mom_loss",
87
JsonWrapper::realValue).asDouble();
90
JsonWrapper::GetProperty(_set_Cut_params, "max_mom_loss",
91
JsonWrapper::realValue).asDouble();
94
JsonWrapper::GetProperty(_set_Cut_params, "min_mass",
95
JsonWrapper::realValue).asDouble();
98
JsonWrapper::GetProperty(_set_Cut_params, "max_mass",
99
JsonWrapper::realValue).asDouble();
102
JsonWrapper::GetProperty(_set_Cut_params, "good_pval",
103
JsonWrapper::realValue).asDouble();
106
void MapCppCuts::_process(MAUS::Data *data) const {
109
throw MAUS::Exceptions::Exception(Exceptions::recoverable,
110
"data was NULL", "MapCppCuts::_process");
113
// Get spill, break if there's no DAQ data
114
Spill *spill = data->GetSpill();
115
if (spill->GetDAQData() == NULL)
117
if (spill->GetDaqEventType() != "physics_event")
120
ReconEventPArray *events = spill->GetReconEvents();
121
int nPartEvents = events->size();
122
for (int i = 0; i < nPartEvents; i++) {
124
// Set up empty vector
125
std::vector<bool> all_cut_values(12, false);
128
Get_single_TOF0_hit_cut((*events)[i]);
129
Get_single_TOF1_hit_cut((*events)[i]);
130
Get_TimeOfFlight((*events)[i]);
131
double tof = Get_TimeOfFlight((*events)[i]);
133
GetTOF_hitTimes_cut(tof, _min_tof, _max_tof);
134
all_cut_values[0] = Get_single_TOF0_hit_cut((*events)[i]);
135
all_cut_values[1] = Get_single_TOF1_hit_cut((*events)[i]);
136
all_cut_values[2] = GetTOF_hitTimes_cut(tof, _min_tof, _max_tof);
139
Get_single_track_cut((*events)[i]);
140
Get_station_hits_cut_US((*events)[i]);
141
all_cut_values[3] = Get_single_track_cut((*events)[i]);
142
all_cut_values[4] = Get_station_hits_cut_US((*events)[i]);
144
// Typically we want momentum nearest to absorber (station 1) (US)
145
process_US_mom((*events)[i]);
146
double US_mom = process_US_mom((*events)[i]);
149
Get_US_mom_cut(US_mom, _min_mom_us, _max_mom_us);
150
Get_momentum_loss_cut(tof, _min_mom_loss, _max_mom_loss, US_mom);
151
Get_p_value_cut((*events)[i], _good_pval);
152
Get_mass_cut(tof, _min_mass, _max_mass, US_mom);
153
all_cut_values[5] = Get_US_mom_cut(US_mom, _min_mom_us, _max_mom_us);
154
all_cut_values[6] = Get_momentum_loss_cut(tof, _min_mom_loss,
155
_max_mom_loss, US_mom);
156
all_cut_values[7] = Get_p_value_cut((*events)[i], _good_pval);
157
all_cut_values[8] = Get_mass_cut(tof, _min_mass, _max_mass,
161
Get_good_particle_cut(all_cut_values[0], all_cut_values[1], all_cut_values[2],
162
all_cut_values[3], all_cut_values[4], all_cut_values[5], all_cut_values[6],
163
all_cut_values[7], all_cut_values[8]);
164
all_cut_values[9] = Get_good_particle_cut(all_cut_values[0], all_cut_values[1],
165
all_cut_values[2], all_cut_values[3], all_cut_values[4], all_cut_values[5],
166
all_cut_values[6], all_cut_values[7], all_cut_values[8]);
168
// Downstream cuts 10, 11
169
Get_station_hits_cut_DS((*events)[i]);
170
all_cut_values[10] = Get_station_hits_cut_DS((*events)[i]);
172
// Typically we want momentum nearest to absorber (station 5) (DS)
173
process_DS_mom((*events)[i]);
174
double DS_mom = process_DS_mom((*events)[i]);
175
Get_DS_mom_cut(DS_mom, _min_mom_ds, _max_mom_ds);
176
all_cut_values[11] = Get_DS_mom_cut(DS_mom, _min_mom_ds, _max_mom_ds);
178
// Pass the vector to the event
179
(*events)[i]->GetCutEvent()->SetCutStore(all_cut_values);
183
// GET TIME OF FLIGHT
184
double MapCppCuts::Get_TimeOfFlight(MAUS::ReconEvent* event) const {
186
MAUS::TOFEvent * MyEvent = event->GetTOFEvent();
187
if (MyEvent == NULL) {
191
MAUS::TOFEventSpacePoint SpacePoint = MyEvent->GetTOFEventSpacePoint();
192
if (SpacePoint.GetTOF0SpacePointArray().size() != 1) return 0;
193
if (SpacePoint.GetTOF1SpacePointArray().size() != 1) return 0;
195
MAUS::TOFSpacePoint TOF0_sp = SpacePoint.GetTOF0SpacePointArray()[0];
196
MAUS::TOFSpacePoint TOF1_sp = SpacePoint.GetTOF1SpacePointArray()[0];
197
double TOF0_hitTime = TOF0_sp.GetTime();
198
double TOF1_hitTime = TOF1_sp.GetTime();
199
double dt = TOF1_hitTime - TOF0_hitTime;
205
std::vector<double> MapCppCuts::Process_momentum(MAUS::ReconEvent* event) const {
207
int tracker, station;
208
double TKU_plane1_px, TKU_plane1_py, TKU_plane1_pz, TKU_plane1_p;
209
double TKU_plane2_px, TKU_plane2_py, TKU_plane2_pz, TKU_plane2_p;
210
double TKU_plane3_px, TKU_plane3_py, TKU_plane3_pz, TKU_plane3_p;
211
double TKU_plane4_px, TKU_plane4_py, TKU_plane4_pz, TKU_plane4_p;
212
double TKU_plane5_px, TKU_plane5_py, TKU_plane5_pz, TKU_plane5_p;
214
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
215
std::vector<double> all_station_mom(5, 0);
216
if (MyEvent == NULL) return all_station_mom;
218
MAUS::ThreeVector momentum;
219
std::vector<MAUS::SciFiTrack*> tracks = MyEvent->scifitracks();
220
std::vector<MAUS::SciFiTrack*>::iterator tr_iter;
221
std::vector<MAUS::SciFiTrackPoint*>::iterator tp_iter;
223
// for (tr_iter = tracks.begin(); tr_iter != tracks.end(); tr_iter++) {
224
if (tracks.size() != 0) {
225
tr_iter = tracks.begin();
226
std::vector<MAUS::SciFiTrackPoint*> tr_points = (*tr_iter)->scifitrackpoints();
227
if (tr_points.size() != 0) {
228
// for (tp_iter = tr_points.begin(); tp_iter != tr_points.end(); tp_iter++) {
229
tp_iter = tr_points.begin();
230
MAUS::SciFiTrackPoint* point = (*tp_iter);
231
tracker = point->tracker();
232
station = point->station();
233
momentum = point->mom();
236
TKU_plane1_px = momentum.x();
237
TKU_plane1_py = momentum.y();
238
TKU_plane1_pz = momentum.z();
239
TKU_plane1_p = sqrt(TKU_plane1_px*TKU_plane1_px
240
+ TKU_plane1_py*TKU_plane1_py
241
+ TKU_plane1_pz*TKU_plane1_pz);
242
all_station_mom[0] = TKU_plane1_p;
243
} else if (station == 2) {
244
TKU_plane2_px = momentum.x();
245
TKU_plane2_py = momentum.y();
246
TKU_plane2_pz = momentum.z();
247
TKU_plane2_p = sqrt(TKU_plane2_px*TKU_plane2_px
248
+ TKU_plane2_py*TKU_plane2_py
249
+ TKU_plane2_pz*TKU_plane2_pz);
250
all_station_mom[1] = TKU_plane2_p;
251
} else if (station == 3) {
252
TKU_plane3_px = momentum.x();
253
TKU_plane3_py = momentum.y();
254
TKU_plane3_pz = momentum.z();
255
TKU_plane3_p = sqrt(TKU_plane3_px*TKU_plane3_px
256
+ TKU_plane3_py*TKU_plane3_py
257
+ TKU_plane3_pz*TKU_plane3_pz);
258
all_station_mom[2] = TKU_plane3_p;
259
} else if (station == 4) {
260
TKU_plane4_px = momentum.x();
261
TKU_plane4_py = momentum.y();
262
TKU_plane4_pz = momentum.z();
263
TKU_plane4_p = sqrt(TKU_plane4_px*TKU_plane4_px
264
+ TKU_plane4_py*TKU_plane4_py
265
+ TKU_plane4_pz*TKU_plane4_pz);
266
all_station_mom[3] = TKU_plane4_p;
267
} else if (station == 5) {
268
TKU_plane5_px = momentum.x();
269
TKU_plane5_py = momentum.y();
270
TKU_plane5_pz = momentum.z();
271
TKU_plane5_p = sqrt(TKU_plane5_px*TKU_plane5_px
272
+ TKU_plane5_py*TKU_plane5_py
273
+ TKU_plane5_pz*TKU_plane5_pz);
274
all_station_mom[4] = TKU_plane5_p;
277
} // trackpoint iterator
279
return all_station_mom;
284
double MapCppCuts::process_US_mom(MAUS::ReconEvent* event) const {
286
int tracker, station, plane;
288
double TKU_plane2_px, TKU_plane2_py, TKU_plane2_pz, TKU_plane2_p;
290
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
291
if (MyEvent == NULL) return mom_US;
293
MAUS::ThreeVector momentum;
294
std::vector<MAUS::SciFiTrack*> tracks = MyEvent->scifitracks();
295
std::vector<MAUS::SciFiTrack*>::iterator track_it;
297
if (tracks.size() != 0) {
298
for (track_it = tracks.begin(); track_it != tracks.end(); track_it++) {
299
std::vector<MAUS::SciFiTrackPoint*> tr_point = (*track_it)->scifitrackpoints();
300
std::vector<MAUS::SciFiTrackPoint*>::iterator tp_iter;
301
for (tp_iter = tr_point.begin(); tp_iter != tr_point.end(); tp_iter++) {
302
MAUS::SciFiTrackPoint* point = (*tp_iter);
303
tracker = point->tracker();
304
plane = point->plane();
305
station = point->station();
306
momentum = point->mom();
307
TKU_plane2_px = momentum.x();
308
TKU_plane2_py = momentum.y();
309
TKU_plane2_pz = momentum.z();
310
TKU_plane2_p = sqrt(TKU_plane2_px*TKU_plane2_px
311
+ TKU_plane2_py*TKU_plane2_py
312
+ TKU_plane2_pz*TKU_plane2_pz);
314
if (tracker == 0 and plane == 2 and station == 5) {
315
mom_US = TKU_plane2_p;
324
double MapCppCuts::process_DS_mom(MAUS::ReconEvent* event) const {
326
int tracker, station, plane;
328
double TKD_plane2_px, TKD_plane2_py, TKD_plane2_pz, TKD_plane2_p;
330
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
331
if (MyEvent == NULL) return mom_DS;
333
MAUS::ThreeVector momentum;
334
std::vector<MAUS::SciFiTrack*> tracks = MyEvent->scifitracks();
335
std::vector<MAUS::SciFiTrack*>::iterator track_it;
337
if (tracks.size() != 0) {
338
for (track_it = tracks.begin(); track_it != tracks.end(); track_it++) {
339
std::vector<MAUS::SciFiTrackPoint*> tr_point = (*track_it)->scifitrackpoints();
340
std::vector<MAUS::SciFiTrackPoint*>::iterator tp_iter;
341
for (tp_iter = tr_point.begin(); tp_iter != tr_point.end(); tp_iter++) {
342
MAUS::SciFiTrackPoint* point = (*tp_iter);
343
tracker = point->tracker();
344
plane = point->plane();
345
station = point->station();
346
momentum = point->mom();
347
TKD_plane2_px = momentum.x();
348
TKD_plane2_py = momentum.y();
349
TKD_plane2_pz = momentum.z();
350
TKD_plane2_p = sqrt(TKD_plane2_px*TKD_plane2_px
351
+ TKD_plane2_py*TKD_plane2_py
352
+ TKD_plane2_pz*TKD_plane2_pz);
354
if (tracker == 1 and plane == 2 and station == 5) {
355
mom_DS = TKD_plane2_p;
364
// CUT CHECK FUNCTIONS
365
bool MapCppCuts::Get_single_TOF0_hit_cut(MAUS::ReconEvent* event) const {
367
bool willCut = false;
369
MAUS::TOFEvent * MyEvent = event->GetTOFEvent();
370
if (MyEvent == NULL) {
374
MAUS::TOFEventSpacePoint SpacePoint = MyEvent->GetTOFEventSpacePoint();
375
if (SpacePoint.GetTOF0SpacePointArray().size() == 1) {
381
bool MapCppCuts::Get_single_TOF1_hit_cut(MAUS::ReconEvent* event) const {
383
bool willCut = false;
385
MAUS::TOFEvent * MyEvent = event->GetTOFEvent();
386
if (MyEvent == NULL) {
390
MAUS::TOFEventSpacePoint SpacePoint = MyEvent->GetTOFEventSpacePoint();
391
if (SpacePoint.GetTOF1SpacePointArray().size() == 1) {
397
bool MapCppCuts::GetTOF_hitTimes_cut(double dt, double min_t,
398
double max_t) const {
400
bool willCut = false;
402
if (dt >= min_t && dt <= max_t) {
408
bool MapCppCuts::Get_single_track_cut(MAUS::ReconEvent* event) const {
410
bool willCut = false;
412
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
413
if (MyEvent == NULL) return willCut;
414
std::vector<MAUS::SciFiTrack*> tracks = MyEvent->scifitracks();
415
if (tracks.size() == 1) {
421
bool MapCppCuts::Get_station_hits_cut_DS(MAUS::ReconEvent* event) const {
423
bool willCut = false;
424
int tracker, station_hits;
426
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
427
if (MyEvent == NULL) return willCut;
429
std::vector<MAUS::SciFiHelicalPRTrack*> pr_tracks = MyEvent->helicalprtracks();
430
std::vector<MAUS::SciFiHelicalPRTrack*>::iterator pr_track_it;
432
if (pr_tracks.size() == 1) {
433
for (pr_track_it = pr_tracks.begin(); pr_track_it != pr_tracks.end(); pr_track_it++) {
434
MAUS::SciFiHelicalPRTrack* track = (*pr_track_it);
435
tracker = track->get_tracker();
437
station_hits = track->get_num_points();
438
if (station_hits == 5) {
448
bool MapCppCuts::Get_station_hits_cut_US(MAUS::ReconEvent* event) const {
450
bool willCut = false;
451
int tracker, station_hits;
453
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
454
if (MyEvent == NULL) return willCut;
456
std::vector<MAUS::SciFiHelicalPRTrack*> pr_tracks = MyEvent->helicalprtracks();
457
std::vector<MAUS::SciFiHelicalPRTrack*>::iterator pr_track_it;
459
if (pr_tracks.size() == 1) {
460
for (pr_track_it = pr_tracks.begin(); pr_track_it != pr_tracks.end(); pr_track_it++) {
461
MAUS::SciFiHelicalPRTrack* track = (*pr_track_it);
462
tracker = track->get_tracker();
464
station_hits = track->get_num_points();
465
if (station_hits == 5) {
475
bool MapCppCuts::Get_US_mom_cut(double mom, double min_mom, double max_mom) const {
477
bool willCut = false;
478
if (mom >= min_mom && mom <= max_mom) {
484
bool MapCppCuts::Get_DS_mom_cut(double mom, double min_mom, double max_mom) const {
486
bool willCut = false;
487
if (mom >= min_mom && mom <= max_mom) {
493
bool MapCppCuts::Get_momentum_loss_cut(double dt, double min_mom_loss,
494
double max_mom_loss, double mom) const {
496
bool willCut = false;
497
if (dt == 0) return willCut;
499
double m = 105.6583715;
501
double beta_tof = dt_e/dt;
503
if (beta_tof > 1.0) return willCut;
504
if (beta_tof < 0) return willCut;
506
double min_mom_cut = (mom + min_mom_loss) / m;
507
double max_mom_cut = (mom + max_mom_loss) / m;
508
double gamma_tof = 1.0/(sqrt(1.0 - beta_tof*beta_tof));
509
double beta_gamma_tof = beta_tof*gamma_tof;
511
if (beta_gamma_tof >= min_mom_cut && beta_gamma_tof <= max_mom_cut) {
517
bool MapCppCuts::Get_p_value_cut(MAUS::ReconEvent* event, double good_pval) const {
522
MAUS::SciFiEvent * MyEvent = event->GetSciFiEvent();
523
if (MyEvent == NULL) return willCut;
525
std::vector<MAUS::SciFiTrack*> tracks = MyEvent->scifitracks();
526
std::vector<MAUS::SciFiTrack*>::iterator track_it;
528
if (tracks.size() == 1) {
529
for (track_it = tracks.begin(); track_it != tracks.end(); track_it++) {
530
MAUS::SciFiTrack* track = (*track_it);
531
tku_pval = track->P_value();
532
if (tku_pval > good_pval) {
541
bool MapCppCuts::Get_mass_cut(double dt, double min_mass, double max_mass,
544
bool willCut = false;
545
if (dt == 0) return willCut;
547
double mom_corr = 18.82;
549
double beta_tof = dt_e/dt;
551
if (beta_tof > 1.0) return willCut;
552
if (beta_tof < 0) return willCut;
554
double gamma_tof = 1.0/(sqrt(1.0 - beta_tof*beta_tof));
555
double beta_gamma_tof = beta_tof*gamma_tof;
556
double mass = (mom + mom_corr)/beta_gamma_tof;
558
if (mass >= min_mass && mass <= max_mass) {
564
bool MapCppCuts::Get_good_particle_cut(bool single_tof0, bool single_tof1, bool good_tof,
565
bool single_track, bool station_hits, bool good_mom, bool good_momLoss,
566
bool good_pval, bool good_massval) const {
568
bool willCut = false;
569
if (single_tof0 == 1 && single_tof1 == 1 && good_tof == 1 &&
570
single_track == 1 && station_hits == 1 && good_mom == 1 &&
571
good_momLoss == 1 && good_pval == 1 && good_massval == 1) {
577
void MapCppCuts::_death() {
added
removed
src/map/MapCppCuts/MapCppCuts.cc
