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/>.
27
#include "src/common_cpp/DataStructure/Hit.hh"
28
#include "src/common_cpp/DataStructure/SciFiEvent.hh"
29
#include "src/common_cpp/DataStructure/SciFiDigit.hh"
30
#include "src/common_cpp/DataStructure/SciFiCluster.hh"
31
#include "src/common_cpp/DataStructure/SciFiSpacePoint.hh"
32
#include "src/common_cpp/DataStructure/SciFiStraightPRTrack.hh"
33
#include "src/common_cpp/DataStructure/SciFiHelicalPRTrack.hh"
34
#include "src/common_cpp/DataStructure/ThreeVector.hh"
35
#include "src/common_cpp/Plotting/SciFi/TrackerDataManager.hh"
36
#include "src/common_cpp/Plotting/SciFi/TrackerDataPlotterXYZ.hh"
40
TrackerDataManager::TrackerDataManager()
41
: _print_tracks(false),
46
TrackerDataManager::~TrackerDataManager() {
50
void TrackerDataManager::process(const Spill *spill) {
51
if (spill != NULL && spill->GetDaqEventType() == "physics_event") {
52
if (_print_tracks || _print_seeds) std::cout << "TrackerDataManager: Opening spill "
53
<< spill->GetSpillNumber() << std::endl;
54
// Loop over recon events
55
for (size_t i = 0; i < spill->GetReconEvents()->size(); ++i) {
56
_t1._spill_num = spill->GetSpillNumber();
57
_t2._spill_num = spill->GetSpillNumber();
58
// Process the SciFi data objects
59
SciFiEvent *evt = (*spill->GetReconEvents())[i]->GetSciFiEvent();
60
process_digits(spill, evt->digits());
61
process_clusters(evt->clusters());
62
process_spoints(evt->spacepoints());
63
process_strks(evt->straightprtracks());
64
process_htrks(evt->helicalprtracks());
65
} // ~ loop over recon events
66
} // ~ check event is a physics event
69
void TrackerDataManager::process_digits(const Spill *spill, const std::vector<SciFiDigit*> digits) {
70
for (size_t i = 0; i < digits.size(); ++i) {
71
SciFiDigit *dig = digits[i];
72
if ( dig->get_tracker() == 0 ) {
74
_t1._num_events = spill->GetReconEvents()->size();
75
} else if ( dig->get_tracker() == 1 ) {
77
_t2._num_events = spill->GetReconEvents()->size();
82
void TrackerDataManager::process_clusters(const std::vector<SciFiCluster*> clusters) {
83
for (size_t i = 0; i < clusters.size(); ++i) {
84
SciFiCluster *clus = clusters[i];
85
if ( clus->get_tracker() == 0 ) {
87
} else if ( clus->get_tracker() == 1 ) {
93
void TrackerDataManager::process_spoints(const std::vector<SciFiSpacePoint*> spoints) {
94
for (size_t i = 0; i < spoints.size(); ++i) {
95
SciFiSpacePoint *sp = spoints[i];
96
ThreeVector pos = sp->get_position();
97
if ( sp->get_tracker() == 0 ) {
99
_t1._spoints_x.push_back(pos.x());
100
_t1._spoints_y.push_back(pos.y());
101
_t1._spoints_z.push_back(pos.z());
102
} else if ( sp->get_tracker() == 1 ) {
104
_t2._spoints_x.push_back(pos.x());
105
_t2._spoints_y.push_back(pos.y());
106
_t2._spoints_z.push_back(pos.z());
111
void TrackerDataManager::process_strks(const std::vector<SciFiStraightPRTrack*> strks) {
112
for (size_t i = 0; i < strks.size(); ++i) {
113
SciFiStraightPRTrack *trk = strks[i];
114
if ( trk->get_tracker() == 0 ) {
115
if ( trk->get_num_points() == 5 ) ++_t1._num_stracks_5pt;
116
if ( trk->get_num_points() == 4 ) ++_t1._num_stracks_4pt;
117
if ( trk->get_num_points() == 3 ) ++_t1._num_stracks_3pt;
118
_t1._trks_str_xz.push_back(make_str_track(trk->get_x0(), trk->get_mx(), _zmin, _zmax));
119
_t1._trks_str_yz.push_back(make_str_track(trk->get_y0(), trk->get_my(), _zmin, _zmax));
120
} else if ( trk->get_tracker() == 1 ) {
121
if ( trk->get_num_points() == 5 ) ++_t2._num_stracks_5pt;
122
if ( trk->get_num_points() == 4 ) ++_t2._num_stracks_4pt;
123
if ( trk->get_num_points() == 3 ) ++_t2._num_stracks_3pt;
124
_t2._trks_str_xz.push_back(make_str_track(trk->get_x0(), trk->get_mx(), _zmin, _zmax));
125
_t2._trks_str_yz.push_back(make_str_track(trk->get_y0(), trk->get_my(), _zmin, _zmax));
130
void TrackerDataManager::process_htrks(const std::vector<SciFiHelicalPRTrack*> htrks) {
131
for (size_t i = 0; i < htrks.size(); ++i) {
132
SciFiHelicalPRTrack *trk = htrks[i];
134
// Print track info to screen
135
if (_print_tracks) print_track_info(trk, i);
137
// Pull out turning angles for each track seed
138
std::vector<double> phi_i;
139
std::vector<double> s_i;
140
for ( size_t i = 0; i < trk->get_phi().size(); ++i ) {
141
phi_i.push_back(trk->get_phi()[i]);
142
s_i.push_back(trk->get_phi()[i]*trk->get_R());
145
// Pull out the z coord for each track seed
146
std::vector<double> z_i;
147
for ( int i = 0; i < (trk->get_spacepoints()->GetLast()+1); ++i ) {
149
static_cast<SciFiSpacePoint*>(trk->get_spacepoints()->At(i))->get_position().z());
152
double x0 = trk->get_circle_x0();
153
double y0 = trk->get_circle_y0();
154
double rad = trk->get_R();
155
double dsdz = trk->get_dsdz();
156
double sz_c = trk->get_line_sz_c();
159
if ( trk->get_tracker() == 0 ) {
160
if ( trk->get_num_points() == 5 ) ++_t1._num_htracks_5pt;
161
if ( trk->get_num_points() == 4 ) ++_t1._num_htracks_4pt;
162
if ( trk->get_num_points() == 3 ) ++_t1._num_htracks_3pt;
163
if ( trk->get_charge() == 1 ) ++_t1._num_pos_tracks;
164
if ( trk->get_charge() == -1 ) ++_t1._num_neg_tracks;
165
_t1._seeds_z.push_back(z_i);
166
_t1._seeds_phi.push_back(phi_i);
167
_t1._seeds_s.push_back(s_i);
168
_t1._trks_xy.push_back(make_circle(x0, y0, rad));
169
// dsdz = - dsdz; // Needed due to the way we plot...
170
_t1._trks_xz.push_back(make_xz(handness, x0, rad, dsdz, sz_c, _zmin, _zmax));
171
_t1._trks_yz.push_back(make_yz(y0, rad, dsdz, sz_c, _zmin, _zmax));
172
_t1._trks_sz.push_back(make_str_track(sz_c, dsdz, _zmin, _zmax));
173
} else if ( trk->get_tracker() == 1 ) {
174
if ( trk->get_num_points() == 5 ) ++_t2._num_htracks_5pt;
175
if ( trk->get_num_points() == 4 ) ++_t2._num_htracks_4pt;
176
if ( trk->get_num_points() == 3 ) ++_t2._num_htracks_3pt;
177
if ( trk->get_charge() == 1 ) ++_t2._num_pos_tracks;
178
if ( trk->get_charge() == -1 ) ++_t2._num_neg_tracks;
179
_t2._seeds_z.push_back(z_i);
180
_t2._seeds_phi.push_back(phi_i);
181
_t2._seeds_s.push_back(s_i);
182
_t2._trks_xy.push_back(make_circle(x0, y0, rad));
183
_t2._trks_xz.push_back(make_xz(handness, x0, rad, dsdz, sz_c, _zmin, _zmax));
184
_t2._trks_yz.push_back(make_yz(y0, rad, dsdz, sz_c, _zmin, _zmax));
185
_t2._trks_sz.push_back(make_str_track(sz_c, dsdz, _zmin, _zmax));
188
// Loop over track seed spacepoints
189
// if (_print_seeds) {
190
// std::cout << "x\ty\tz\ttime\t\tphi\tpx_mc\tpy_mc\tpt_mc\tpz_mc\t";
191
// std::cout << "px_mc1\tpx_mc2\tpx_mc3\tpy_mc1\tpy_mc2\tpy_mc3\tpz_mc1\tpz_mc2\tpz_mc3\n";
194
for ( int j = 0; j < (trk->get_spacepoints()->GetLast()+1); ++j ) {
195
if ( trk->get_tracker() == 0 ) ++_t1._num_seeds;
196
if ( trk->get_tracker() == 1 ) ++_t2._num_seeds;
197
// if (_print_seeds) print_seed_info(trk, j);
199
} // ~loop over helical tracks
202
void TrackerDataManager::clear_spill() {
207
void TrackerDataManager::draw(std::vector<TrackerDataPlotterBase*> plotters) {
208
// Loop over all the plotters and draw
209
TCanvas* lCanvas = NULL;
210
for ( size_t i = 0; i < plotters.size(); ++i ) {
211
TrackerDataPlotterBase * plt = plotters[i];
213
lCanvas = (*plt)(_t1, _t2);
214
if (plt->GetSaveOutput()) {
215
std::string fName = plt->GetOutputName();
216
lCanvas->SaveAs(fName.c_str());
219
std::cerr << "Error: Empty plotter pointer passed\n";
224
TArc TrackerDataManager::make_circle(double x0, double y0, double rad) {
225
TArc arc = TArc(x0, y0, rad);
226
arc.SetFillStyle(0); // 0 - Transparent
227
arc.SetLineColor(kBlue);
231
TF1 TrackerDataManager::make_str_track(double c, double m, double zmin, double zmax) {
232
// Note: in the function expression, x is just the independent variable, which
233
// in this case is the z coordinate in the tracker coordinate system
234
TF1 trk = TF1("trk", "[0]+([1]*x)", zmin, zmax);
235
trk.SetParameters(c, m);
236
trk.SetLineColor(kRed);
240
TF1 TrackerDataManager::make_xz(int handness, double circle_x0, double rad, double dsdz,
241
double sz_c, double zmin, double zmax) {
242
// The x in the cos term is actually representing z (the indep variable)
243
TF1 func = TF1("xz_func", "[0]-[4]*([1]*cos((1/[1])*([2]*x+[3])))", zmin, zmax);
244
func.SetParameter(0, circle_x0);
245
func.SetParameter(1, rad);
246
func.SetParameter(2, dsdz);
247
func.SetParameter(3, sz_c);
248
func.SetParameter(4, handness);
249
func.SetLineColor(kBlue);
253
TF1 TrackerDataManager::make_yz(double circle_y0, double rad, double dsdz,
254
double sz_c, double zmin, double zmax) {
255
// The x in the cos term is actually representing z (the indep variable)
256
TF1 func = TF1("yz_func", "[0]+([1]*sin((1/[1])*([2]*x+[3])))", zmin, zmax);
257
func.SetParameter(0, circle_y0);
258
func.SetParameter(1, rad);
259
func.SetParameter(2, dsdz);
260
func.SetParameter(3, sz_c);
261
func.SetLineColor(kBlue);
265
void TrackerDataManager::print_track_info(const SciFiHelicalPRTrack * const trk, int trk_num) {
266
double rad = trk->get_R();
267
double pt = rad * _RtoPt;
268
double pz = pt / trk->get_dsdz();
269
std::cout << "Tracker " << trk->get_tracker() << + ", ";
270
std::cout << "Track " << trk_num << ", ";
271
std::cout << "Num points " << trk->get_num_points() << ", ";
272
std::cout << "Charge " << trk->get_charge() << ", ";
273
std::cout << "R = " << std::setprecision(4) << trk->get_R() << "mm, ";
274
std::cout << "X0 = " << std::setprecision(4) << trk->get_circle_x0() << "mm, ";
275
std::cout << "Y0 = " << std::setprecision(4) << trk->get_circle_y0() << "mm,\n";
276
std::cout << "dsdz " << std::setprecision(4) << trk->get_dsdz() << ", ";
277
std::cout << "pt = " << std::setprecision(4) << pt << "MeV/c, ";
278
std::cout << "pz = " << std::setprecision(4) << pz << "MeV/c, ";
279
std::cout << "xy_chi2 = " << std::setprecision(4) << trk->get_circle_chisq() << ", ";
280
std::cout << "sz_c = " << std::setprecision(4) << trk->get_line_sz_c() << ", ";
281
std::cout << "sz_chi2 = " << std::setprecision(4) << trk->get_line_sz_chisq() << "\n";
284
// void TrackerDataManager::print_seed_info(const SciFiHelicalPRTrack * const trk, int seed_num) {
285
// SciFiSpacePoint* sp = static_cast<SciFiSpacePoint*>(trk->get_spacepoints()->At(seed_num));
286
// MAUS::ThreeVector pos = sp->get_position();
287
// std::vector<MAUS::SciFiCluster*> chans = sp->get_channels_pointers();
288
// double t = sp->get_time();
290
// double seed_px_mc = 0.0;
291
// double seed_py_mc = 0.0;
292
// double seed_pz_mc = 0.0;
293
// double px_mc1 = 0.0;
294
// double px_mc2 = 0.0;
295
// double px_mc3 = 0.0;
296
// double py_mc1 = 0.0;
297
// double py_mc2 = 0.0;
298
// double py_mc3 = 0.0;
299
// double pz_mc1 = 0.0;
300
// double pz_mc2 = 0.0;
301
// double pz_mc3 = 0.0;
302
// for (size_t iCl = 0; iCl < chans.size(); ++iCl) {
303
// ThreeVector p_mc = chans[iCl]->get_true_momentum();
304
// seed_px_mc += p_mc.x();
305
// seed_py_mc += p_mc.y();
306
// seed_pz_mc += p_mc.z();
308
// px_mc1 = p_mc.x();
309
// py_mc1 = p_mc.y();
310
// pz_mc1 = p_mc.z();
313
// px_mc2 = p_mc.x();
314
// py_mc2 = p_mc.y();
315
// pz_mc2 = p_mc.z();
318
// px_mc3 = p_mc.x();
319
// py_mc3 = p_mc.y();
320
// pz_mc3 = p_mc.z();
323
// seed_px_mc /= chans.size();
324
// seed_py_mc /= chans.size();
325
// seed_pz_mc /= chans.size();
326
// double pt_mc = sqrt(seed_px_mc*seed_px_mc + seed_py_mc*seed_py_mc);
328
// std::cout << std::setprecision(4) << pos.x() << "\t";
329
// std::cout << std::setprecision(4) << pos.y() << "\t";
330
// std::cout << std::setprecision(4) << pos.z() << "\t";
331
// std::cout << std::setprecision(12) << t << "\t";
332
// std::cout << std::setprecision(4) << trk->get_phi()[seed_num] << "\t";
333
// std::cout << std::setprecision(4) << seed_px_mc << "\t";
334
// std::cout << std::setprecision(4) << seed_py_mc << "\t";
335
// std::cout << std::setprecision(4) << pt_mc << "\t";
336
// std::cout << std::setprecision(4) << seed_pz_mc << "\t";
337
// std::cout << std::setprecision(4) << px_mc1 << "\t";
338
// std::cout << std::setprecision(4) << px_mc2 << "\t";
339
// std::cout << std::setprecision(4) << px_mc3 << "\t";
340
// std::cout << std::setprecision(4) << py_mc1 << "\t";
341
// std::cout << std::setprecision(4) << py_mc2 << "\t";
342
// std::cout << std::setprecision(4) << py_mc3 << "\t";
343
// std::cout << std::setprecision(4) << pz_mc1 << "\t";
344
// std::cout << std::setprecision(4) << pz_mc2 << "\t";
345
// std::cout << std::setprecision(4) << pz_mc3 << "\n";
added
removed
src/common_cpp/Plotting/SciFi/TrackerDataManager.cc
