← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

Viewing changes to src/common_cpp/Plotting/SciFi/TrackerDataAnalyserMomentum.cc

  • Committer: Durga Rajaram
  • Date: 2014-07-16 15:13:05 UTC
  • mfrom: (659.1.92 cand)
  • Revision ID: durga@fnal.gov-20140716151305-q27rv1y9p03v9lks
Tags: MAUS-v0.9, MAUS-v0.9.0
MAUS-v0.9.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Plotting/SciFi/TrackerDataAnalyserMomentum.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
 
 
18
// C++ headers
 
19
#include <cmath>
 
20
#include <map>
 
21
#include <sstream>
 
22
 
 
23
// ROOT headers
 
24
#include "TAxis.h"
 
25
#include "TVirtualPad.h"
 
26
#include "TStyle.h"
 
27
#include "TF1.h"
 
28
#include "TH1.h"
 
29
#include "TMath.h"
 
30
#include "Rtypes.h"
 
31
#include "TRefArray.h"
 
32
#include "TRef.h"
 
33
 
 
34
// MAUS headers
 
35
#include "src/common_cpp/Plotting/SciFi/TrackerDataAnalyserMomentum.hh"
 
36
#include "src/common_cpp/DataStructure/Spill.hh"
 
37
#include "src/common_cpp/DataStructure/Data.hh"
 
38
#include "src/common_cpp/DataStructure/MCEvent.hh"
 
39
#include "src/common_cpp/DataStructure/ReconEvent.hh"
 
40
#include "src/common_cpp/DataStructure/SciFiEvent.hh"
 
41
#include "src/common_cpp/DataStructure/SciFiSpacePoint.hh"
 
42
#include "src/common_cpp/DataStructure/SciFiStraightPRTrack.hh"
 
43
#include "src/common_cpp/DataStructure/SciFiHelicalPRTrack.hh"
 
44
#include "src/common_cpp/DataStructure/ThreeVector.hh"
 
45
#include "src/common_cpp/Recon/SciFi/SimpleLine.hh"
 
46
#include "src/common_cpp/Recon/SciFi/SimpleCircle.hh"
 
47
 
 
48
namespace MAUS {
 
49
 
 
50
TrackerDataAnalyserMomentum::TrackerDataAnalyserMomentum()
 
51
  : _spill_num(0),
 
52
    _tracker_num(0),
 
53
    _charge(0),
 
54
    _num_points(0),
 
55
    _n_bad_tracks(0),
 
56
    _mc_track_id(0),
 
57
    _mc_pid(0),
 
58
    _n_matched(0),
 
59
    _n_mismatched(0),
 
60
    _n_missed(0),
 
61
    _pt_rec(0.0),
 
62
    _pz_rec(0.0),
 
63
    _pt_mc(0.0),
 
64
    _pz_mc(0.0),
 
65
    _out_file(NULL),
 
66
    _tree(NULL),
 
67
    _t1_pt_res(NULL),
 
68
    _t1_pz_res(NULL),
 
69
    _t1_pz_res_log(NULL),
 
70
    _t2_pt_res(NULL),
 
71
    _t2_pz_res(NULL),
 
72
    _t2_pz_res_log(NULL),
 
73
    _t1_pt_res_pt(NULL),
 
74
    _t1_pz_res_pt(NULL),
 
75
    _t2_pt_res_pt(NULL),
 
76
    _t2_pz_res_pt(NULL),
 
77
    _t1_pz_resol(NULL),
 
78
    _t2_pz_resol(NULL),
 
79
    _cResiduals(NULL),
 
80
    _cGraphs(NULL),
 
81
    _cResolutions(NULL),
 
82
    _cutPzRec(0.0) {
 
83
  // Do nothing
 
84
}
 
85
 
 
86
TrackerDataAnalyserMomentum::~TrackerDataAnalyserMomentum() {
 
87
  if (_out_file) delete _out_file;
 
88
  if (_tree) delete _tree;
 
89
  if (_t1_pt_res) delete _t1_pt_res;
 
90
  if (_t1_pz_res) delete _t1_pz_res;
 
91
  if (_t1_pz_res_log) delete _t1_pz_res_log;
 
92
  if (_t2_pt_res) delete _t2_pt_res;
 
93
  if (_t2_pz_res) delete _t2_pz_res;
 
94
  if (_t2_pz_res_log) delete _t2_pz_res_log;
 
95
  if (_t1_pt_res_pt) delete _t1_pz_res_pt;
 
96
  if (_t1_pz_res_pt) delete _t1_pt_res_pt;
 
97
  if (_t2_pt_res_pt) delete _t2_pz_res_pt;
 
98
  if (_t2_pz_res_pt) delete _t2_pt_res_pt;
 
99
  if (_t1_pz_resol) delete _t1_pz_resol;
 
100
  if (_t2_pz_resol) delete _t2_pz_resol;
 
101
  if (_cResiduals) delete _cResiduals;
 
102
  if (_cGraphs) delete _cGraphs;
 
103
  if (_cResolutions) delete _cResolutions;
 
104
}
 
105
 
 
106
void TrackerDataAnalyserMomentum::setUp() {
 
107
  // Set up the output ROOT file
 
108
  _out_file = new TFile("momentum_analysis_output.root", "recreate");
 
109
 
 
110
  // Set up the output tree
 
111
  _tree = new TTree("tree", "Momentum Residual Data");
 
112
  _tree->Branch("spill_num", &_spill_num, "spill_num/I");
 
113
  _tree->Branch("tracker_num", &_tracker_num, "tracker_num/I");
 
114
  _tree->Branch("charge", &_charge, "charge/I");
 
115
  _tree->Branch("num_points", &_num_points, "num_points/I");
 
116
  _tree->Branch("n_bad_tracks", &_n_bad_tracks, "n_bad_tracks/I");
 
117
  _tree->Branch("mc_track_id", &_mc_track_id, "mc_track_id/I");
 
118
  _tree->Branch("mc_pid", &_mc_pid, "mc_pid/I");
 
119
  _tree->Branch("n_matched", &_n_matched, "n_matched/I");
 
120
  _tree->Branch("n_mismatched", &_n_mismatched, "n_mismatched/I");
 
121
  _tree->Branch("n_missed", &_n_missed, "n_missed/I");
 
122
  _tree->Branch("pt_rec", &_pt_rec, "pt_rec/D");
 
123
  _tree->Branch("pz_rec", &_pz_rec, "pz_rec/D");
 
124
  _tree->Branch("pt_mc", &_pt_mc, "pt_mc/D");
 
125
  _tree->Branch("pz_mc", &_pz_mc, "pz_mc/D");
 
126
 
 
127
  // Set up the residual histograms
 
128
  int res_n_bins = 100;
 
129
 
 
130
  _t1_pt_res = new TH1D("t1_pt_res", "T1 p_{t} Residual", res_n_bins, -5, 5);
 
131
  _t1_pt_res->GetXaxis()->SetTitle("p_{T}^{MC} - p_{T} (MeV/c)");
 
132
 
 
133
  _t1_pz_res = new TH1D("t1_pz_res", "T1 p_{z} Residual", res_n_bins, -30, 30);
 
134
  _t1_pz_res->GetXaxis()->SetTitle("p_{z}^{MC} - p_{z} (MeV/c)");
 
135
 
 
136
  _t1_pz_res_log = new TH1D("t1_pz_res_log", "T1 p_{z} Residual", res_n_bins, -500, 500);
 
137
  _t1_pz_res_log->GetXaxis()->SetTitle("p_{z}^{MC} - p_{z} (MeV/c)");
 
138
 
 
139
  _t2_pt_res = new TH1D("t2_pt_res", "T2 p_{t} Residual", res_n_bins, -5, 5);
 
140
  _t2_pt_res->GetXaxis()->SetTitle("p_{T}^{MC} - p_{T} (MeV/c)");
 
141
 
 
142
  _t2_pz_res = new TH1D("t2_pz_res", "T2 p_{z} Residual", res_n_bins, -30, 30);
 
143
  _t2_pz_res->GetXaxis()->SetTitle("p_{z}^{MC} - p_{z} (MeV/c)");
 
144
 
 
145
  _t2_pz_res_log = new TH1D("t2_pz_res_log", "T2 p_{z} Residual", res_n_bins, -500, 500);
 
146
  _t2_pz_res_log->GetXaxis()->SetTitle("p_{z}^{MC} - p_{z} (MeV/c)");
 
147
 
 
148
  // Set up the residual graphs
 
149
  _t1_pz_res_pt = new TGraph();
 
150
  _t2_pz_res_pt = new TGraph();
 
151
 
 
152
  // Set global styles
 
153
  gStyle->SetOptStat(111111);
 
154
  gStyle->SetLabelSize(0.05, "XYZ");
 
155
  gStyle->SetTitleOffset(1.15, "X");
 
156
  gStyle->SetStatW(0.3);
 
157
  gStyle->SetStatH(0.2);
 
158
}
 
159
 
 
160
void TrackerDataAnalyserMomentum::accumulate(Spill* spill) {
 
161
  if (spill != NULL && spill->GetDaqEventType() == "physics_event") {
 
162
    _spill_num = spill->GetSpillNumber();
 
163
 
 
164
    // Loop over MC events in the spill
 
165
    for (size_t iMevt = 0; iMevt < spill->GetMCEvents()->size(); ++iMevt) {
 
166
      MCEvent *mc_evt = (*spill->GetMCEvents())[iMevt];
 
167
      SciFiEvent *evt = (*spill->GetReconEvents())[iMevt]->GetSciFiEvent();
 
168
      calc_pat_rec_efficiency(mc_evt, evt);
 
169
    }
 
170
  } else {
 
171
    std::cout << "Not a usable spill" << std::endl;
 
172
  }
 
173
}
 
174
 
 
175
void TrackerDataAnalyserMomentum::calc_pat_rec_efficiency(MCEvent *mc_evt, SciFiEvent* evt) {
 
176
  std::vector<SciFiHelicalPRTrack*> htrks = evt->helicalprtracks();
 
177
 
 
178
  // Create the lookup bridge between MC and Recon
 
179
  SciFiLookup lkup;
 
180
  lkup.make_hits_map(mc_evt);
 
181
 
 
182
  // Reset tracks counters
 
183
  _n_bad_tracks = 0;
 
184
 
 
185
  // Loop over helical pattern recognition tracks in event
 
186
  for (size_t iTrk = 0; iTrk < htrks.size(); ++iTrk) {
 
187
    SciFiHelicalPRTrack* trk = htrks[iTrk];
 
188
    _pt_rec = 0.0;
 
189
    _pz_rec = 0.0;
 
190
    if ((trk->get_R() != 0.0) & (trk->get_dsdz() != 0.0)) {
 
191
      _num_points = trk->get_num_points();
 
192
      std::vector< std::vector<int> > spoint_mc_tids;  // Vector of MC track ids for each spoint
 
193
 
 
194
      // Calc recon momentum
 
195
      _pt_rec = 1.2 * trk->get_R();
 
196
      _pz_rec = _pt_rec / trk->get_dsdz();
 
197
 
 
198
      // Calc MC momentum
 
199
      // Loop over seed spacepoints associated with the track
 
200
      std::vector<SciFiSpacePoint*> spnts = trk->get_spacepoints_pointers();
 
201
      for ( size_t k = 0; k < spnts.size(); ++k ) {
 
202
        std::map<int, int> track_id_counter;  // Map of track id to freq for each spoint
 
203
        std::vector<SciFiCluster*> clusters = spnts[k]->get_channels_pointers();
 
204
 
 
205
        // Loop over clusters
 
206
        for ( size_t l = 0; l < clusters.size(); ++l ) {
 
207
          std::vector<SciFiDigit*> digits = clusters[l]->get_digits_pointers();
 
208
 
 
209
          // Loop over digits
 
210
          for ( size_t m = 0; m < digits.size(); ++m ) {
 
211
            // Perform the digits to hits lookup
 
212
            std::vector<SciFiHit*> hits;
 
213
            if (!lkup.get_hits(digits[m], hits)) {
 
214
              std::cerr << "Lookup failed\n";
 
215
              continue;
 
216
            }
 
217
 
 
218
            // Loop over MC hits
 
219
            for ( size_t n = 0; n < hits.size(); ++n ) {
 
220
              int track_id = hits[n]->GetTrackId();
 
221
              if (track_id_counter.count(track_id))
 
222
                track_id_counter[track_id] = track_id_counter[track_id] + 1;
 
223
              else
 
224
                track_id_counter[track_id] = 1;
 
225
            } // ~Loop over MC hits
 
226
          } // ~Loop over digits
 
227
        } // ~// Loop over clusters
 
228
 
 
229
        // Find the track_ids for this spacepoint
 
230
        std::map<int, int>::iterator mit1;
 
231
        std::vector<int> mc_tracks_ids;
 
232
        for ( mit1 = track_id_counter.begin(); mit1 != track_id_counter.end(); ++mit1 ) {
 
233
          mc_tracks_ids.push_back(mit1->first);
 
234
        }
 
235
        spoint_mc_tids.push_back(mc_tracks_ids);
 
236
      } // ~Loop over seed spacepoints
 
237
 
 
238
      // Is there a track id associated with 3 or more spoints?
 
239
      int track_id = 0;
 
240
      int counter = 0;
 
241
      bool success = find_mc_tid(spoint_mc_tids, track_id, counter);
 
242
      // If we have not found a common track id, abort for this track
 
243
      if (!success) {
 
244
        std::cerr << "Malformed track, skipping\n";
 
245
        ++_n_bad_tracks;
 
246
        break;
 
247
      }
 
248
      // If we have found a common track id amoung the spoints, fill the tree
 
249
      _mc_track_id = track_id;
 
250
      _n_matched = counter;
 
251
      _n_mismatched = spoint_mc_tids.size() - counter;
 
252
      _n_missed = 5 - counter; // TODO: improve
 
253
 
 
254
      // Calc the MC track momentum using hits only with this track id
 
255
      _pt_mc = 0.0;
 
256
      _pz_mc = 0.0;
 
257
      int counter2 = 0;
 
258
      for ( size_t k = 0; k < spnts.size(); ++k ) {
 
259
        ThreeVector mom_mc;
 
260
        bool success = find_mc_spoint_momentum(track_id, spnts[k], lkup, mom_mc);
 
261
        if (!success) continue;
 
262
        _pt_mc += sqrt(mom_mc.x()*mom_mc.x() + mom_mc.y()*mom_mc.y());
 
263
        _pz_mc += mom_mc.z();
 
264
        ++counter2;
 
265
      }
 
266
      _pt_mc /= counter2;
 
267
      _pz_mc /= counter2;
 
268
 
 
269
      // Fill the vectors and tree with the extracted recon and MC momenta
 
270
      if (trk->get_tracker() == 0) {
 
271
        _t1_pt_res->Fill(_pt_mc - _pt_rec);
 
272
        _t1_pz_res->Fill(_pz_mc + trk->get_charge()*_pz_rec);
 
273
        _t1_pz_res_log->Fill(_pz_mc + trk->get_charge()*_pz_rec);
 
274
        _vec_t1_pt_mc.push_back(_pt_mc);
 
275
        _vec_t1_pt_res.push_back(_pt_mc - _pt_rec);
 
276
        _vec_t1_pz.push_back(_pz_rec);
 
277
        _vec_t1_pz_res.push_back(_pz_mc + trk->get_charge()*_pz_rec);
 
278
      } else if (trk->get_tracker() == 1) {
 
279
        _t2_pt_res->Fill(_pt_mc - _pt_rec);
 
280
        _t2_pz_res->Fill(_pz_mc - trk->get_charge()*_pz_rec);
 
281
        _t2_pz_res_log->Fill(_pz_mc - trk->get_charge()*_pz_rec);
 
282
        _vec_t2_ptMc.push_back(_pt_mc);
 
283
        _vec_t2_pt_res.push_back(_pt_mc - _pt_rec);
 
284
        _vec_t2_pz.push_back(_pz_rec);
 
285
        _vec_t2_pz_res.push_back(_pz_mc - trk->get_charge()*_pz_rec);
 
286
      }
 
287
      _tracker_num = trk->get_tracker();
 
288
      _charge = trk->get_charge();
 
289
      _tree->Fill();
 
290
    } else {
 
291
      std::cout << "Bad track, skipping" << std::endl;
 
292
    }
 
293
  }
 
294
}
 
295
 
 
296
bool TrackerDataAnalyserMomentum::find_mc_spoint_momentum(const int track_id,
 
297
                const SciFiSpacePoint* sp, SciFiLookup &lkup, ThreeVector &mom) {
 
298
 
 
299
  std::vector<SciFiCluster*> clusters = sp->get_channels_pointers();
 
300
  std::vector<SciFiCluster*>::iterator clus_it;
 
301
  for (clus_it = clusters.begin(); clus_it != clusters.end(); ++clus_it) {
 
302
    std::vector<SciFiDigit*> digits = (*clus_it)->get_digits_pointers();
 
303
    std::vector<SciFiDigit*>::iterator dig_it;
 
304
    for (dig_it = digits.begin(); dig_it != digits.end(); ++dig_it) {
 
305
      std::vector<SciFiHit*> hits;
 
306
      if (!lkup.get_hits((*dig_it), hits)) {
 
307
        std::cerr << "Lookup failed\n";
 
308
        continue;
 
309
      }
 
310
      std::vector<SciFiHit*>::iterator hit_it;
 
311
      for (hit_it = hits.begin(); hit_it != hits.end(); ++hit_it) {
 
312
        if (track_id == (*hit_it)->GetTrackId()) {
 
313
          mom = (*hit_it)->GetMomentum();
 
314
          return true;
 
315
        }
 
316
      }
 
317
    }
 
318
  }
 
319
  return false;
 
320
};
 
321
 
 
322
bool TrackerDataAnalyserMomentum::find_mc_tid(const std::vector< std::vector<int> > &spoint_mc_tids,
 
323
                                              int &tid, int &counter) {
 
324
 
 
325
  std::map<int, int> track_id_counter;  // Map of track id to freq for each spoint
 
326
 
 
327
  for ( size_t i = 0; i < spoint_mc_tids.size(); ++i ) {
 
328
    for ( size_t j = 0; j < spoint_mc_tids[i].size(); ++j ) {
 
329
      int current_tid = spoint_mc_tids[i][j];
 
330
      if (track_id_counter.count(current_tid))
 
331
        track_id_counter[current_tid] = track_id_counter[current_tid] + 1;
 
332
      else
 
333
        track_id_counter[current_tid] = 1;
 
334
    }
 
335
  }
 
336
 
 
337
  std::map<int, int>::iterator mit1;
 
338
  tid = 0;
 
339
  for ( mit1 = track_id_counter.begin(); mit1 != track_id_counter.end(); ++mit1 ) {
 
340
    if ( mit1->second > 2 && tid == 0 ) {
 
341
      tid = mit1->first;
 
342
      counter = mit1->second;
 
343
    } else if ( mit1->second > 2 && tid != 0 ) {
 
344
      tid = -1;  // A malformed track, cannot asscoiate with an mc track id
 
345
      counter = 0;
 
346
      return false;
 
347
    }
 
348
  }
 
349
  return true;
 
350
};
 
351
 
 
352
void TrackerDataAnalyserMomentum::make_all() {
 
353
  make_residual_histograms();
 
354
  make_residual_graphs();
 
355
  make_pz_resolutions();
 
356
  make_resolution_graphs();
 
357
}
 
358
 
 
359
void TrackerDataAnalyserMomentum::make_pz_resolutions() {
 
360
  int nPoints = 9;                        // The number of MC momentum intervals used in the plots
 
361
  std::vector<TCut> vCuts(nPoints);        // The cuts defining the pt_mc intervals
 
362
  std::vector<double> vPtMc(nPoints);        // The centre of the pt_mc intervals
 
363
  std::vector<double> vPtMcErr(nPoints);     // The width of the intervals
 
364
  std::vector<double> vPzRes_t1(nPoints);    // The resulatant pz resolutions
 
365
  std::vector<double> vPzResErr_t1(nPoints); // The errors assocaited with the pz res
 
366
  std::vector<double> vPzRes_t2(nPoints);    // The resulatant pz resolutions
 
367
  std::vector<double> vPzResErr_t2(nPoints); // The errors assocaited with the pz res
 
368
 
 
369
  // Cuts in pt_mc
 
370
  vCuts[0] = "pt_mc>=0&&pt_mc<10";
 
371
  vCuts[1] = "pt_mc>=10&&pt_mc<20";
 
372
  vCuts[2] = "pt_mc>=20&&pt_mc<30";
 
373
  vCuts[3] = "pt_mc>=30&&pt_mc<40";
 
374
  vCuts[4] = "pt_mc>=40&&pt_mc<50";
 
375
  vCuts[5] = "pt_mc>=50&&pt_mc<60";
 
376
  vCuts[6] = "pt_mc>=60&&pt_mc<70";
 
377
  vCuts[7] = "pt_mc>=70&&pt_mc<80";
 
378
  vCuts[8] = "pt_mc>=80&&pt_mc<90";
 
379
  // vCuts[9] = "pt_mc>=90&&pt_mc<100";
 
380
 
 
381
  // The central MC momentum
 
382
  vPtMc[0] = 5.0;
 
383
  vPtMc[1] = 15.0;
 
384
  vPtMc[2] = 25.0;
 
385
  vPtMc[3] = 35.0;
 
386
  vPtMc[4] = 45.0;
 
387
  vPtMc[5] = 55.0;
 
388
  vPtMc[6] = 65.0;
 
389
  vPtMc[7] = 75.0;
 
390
  vPtMc[8] = 85.0;
 
391
  // vPtMc[9] = 95.0;
 
392
 
 
393
  // The error associated with the MC momentum (just the interval half width)
 
394
  vPtMcErr[0] = 5.0;
 
395
  vPtMcErr[1] = 5.0;
 
396
  vPtMcErr[2] = 5.0;
 
397
  vPtMcErr[3] = 5.0;
 
398
  vPtMcErr[4] = 5.0;
 
399
  vPtMcErr[5] = 5.0;
 
400
  vPtMcErr[6] = 5.0;
 
401
  vPtMcErr[7] = 5.0;
 
402
  vPtMcErr[8] = 5.0;
 
403
  // vPtMcErr[9] = 5.0;
 
404
 
 
405
  // Cuts for to select each tracker
 
406
  TCut cutT1 = "tracker_num==0";
 
407
  TCut cutT2 = "tracker_num==1";
 
408
 
 
409
  // Form the reconstructed pz TCut
 
410
  TCut tcut_pzrec = formTCut("TMath::Abs(pz_rec)", "<", _cutPzRec);
 
411
 
 
412
  // Loop over the mometum intervals and calculate the resolution for each
 
413
  for (int i = 0; i < nPoints; ++i) {
 
414
    TCut input_cut = vCuts[i]&&cutT1&&tcut_pzrec;
 
415
    calc_pz_resolution(0, input_cut, vPzRes_t1[i], vPzResErr_t1[i]);
 
416
  }
 
417
  for (int i = 0; i < nPoints; ++i) {
 
418
    TCut input_cut = vCuts[i]&&cutT2&&tcut_pzrec;
 
419
    calc_pz_resolution(1, input_cut, vPzRes_t2[i], vPzResErr_t2[i]);
 
420
  }
 
421
 
 
422
  // Create the resultant resolution plots
 
423
  _t1_pz_resol = new TGraphErrors(nPoints, &vPtMc[0], &vPzRes_t1[0],
 
424
                                  &vPtMcErr[0], &vPzResErr_t1[0]);
 
425
  _t2_pz_resol = new TGraphErrors(nPoints, &vPtMc[0], &vPzRes_t2[0],
 
426
                                  &vPtMcErr[0], &vPzResErr_t2[0]);
 
427
}
 
428
 
 
429
void TrackerDataAnalyserMomentum::make_residual_histograms() {
 
430
  _cResiduals = new TCanvas("cResiduals", "Mometum Residuals");
 
431
  _cResiduals->Divide(3, 2);
 
432
  _cResiduals->cd(1);
 
433
  _t1_pt_res->UseCurrentStyle();
 
434
  _t1_pt_res->Draw();
 
435
  _cResiduals->cd(2);
 
436
  _t1_pz_res->UseCurrentStyle();
 
437
  _t1_pz_res->Draw();
 
438
  TVirtualPad* pad = _cResiduals->cd(3);
 
439
  pad->SetLogy();
 
440
  _t1_pz_res_log->UseCurrentStyle();
 
441
  _t1_pz_res_log->Draw();
 
442
  _cResiduals->cd(4);
 
443
  _t2_pt_res->UseCurrentStyle();
 
444
  _t2_pt_res->Draw();
 
445
  _cResiduals->cd(5);
 
446
  _t2_pz_res->UseCurrentStyle();
 
447
  _t2_pz_res->Draw();
 
448
  pad = _cResiduals->cd(6);
 
449
  pad->SetLogy();
 
450
  _t2_pz_res_log->UseCurrentStyle();
 
451
  _t2_pz_res_log->Draw();
 
452
}
 
453
 
 
454
void TrackerDataAnalyserMomentum::make_residual_graphs() {
 
455
  _cGraphs = new TCanvas("cGraphs", "Mometum Graphs");
 
456
  _cGraphs->Divide(2, 2);
 
457
 
 
458
  _cGraphs->cd(1);
 
459
  _t1_pt_res_pt = new TGraph(_vec_t1_pt_mc.size(), &_vec_t1_pt_mc[0], &_vec_t1_pt_res[0]);
 
460
  _t1_pt_res_pt->SetName("t1_pt_res_pt");
 
461
  _t1_pt_res_pt->SetTitle("T1 p_{t} res vs. p_{t}^{MC}");
 
462
  _t1_pt_res_pt->GetXaxis()->SetTitle("p_{t}^{MC} (MeV/c)");
 
463
  _t1_pt_res_pt->GetYaxis()->SetTitle("p_{t}^{MC} - p_{t} (MeV/c)");
 
464
  _t1_pt_res_pt->GetYaxis()->SetTitleOffset(1.4);
 
465
  _t1_pt_res_pt->Draw("AP");
 
466
 
 
467
  _cGraphs->cd(2);
 
468
  _t1_pz_res_pt = new TGraph(_vec_t1_pt_mc.size(), &_vec_t1_pt_mc[0], &_vec_t1_pz_res[0]);
 
469
  _t1_pz_res_pt->SetName("t1_pz_res_pt");
 
470
  _t1_pz_res_pt->SetTitle("T1 p_{z} res vs. p_{t}^{MC}");
 
471
  _t1_pz_res_pt->GetXaxis()->SetTitle("p_{t}^{MC} (MeV/c)");
 
472
  _t1_pz_res_pt->GetYaxis()->SetTitle("p_{z}^{MC} - p_{z} (MeV/c)");
 
473
  _t1_pz_res_pt->GetYaxis()->SetTitleOffset(1.4);
 
474
  _t1_pz_res_pt->GetYaxis()->SetRangeUser(-200, 200);
 
475
  _t1_pz_res_pt->Draw("AP");
 
476
 
 
477
  _cGraphs->cd(3);
 
478
  _t2_pt_res_pt = new TGraph(_vec_t2_ptMc.size(), &_vec_t2_ptMc[0], &_vec_t2_pt_res[0]);
 
479
  _t2_pt_res_pt->SetName("t2_pt_res_pt");
 
480
  _t2_pt_res_pt->SetTitle("T2 p_{t} res vs. p_{t}^{MC}");
 
481
  _t2_pt_res_pt->GetXaxis()->SetTitle("p_{t}^{MC} (MeV/c)");
 
482
  _t2_pt_res_pt->GetYaxis()->SetTitle("p_{t}^{MC} - p_{t} (MeV/c)");
 
483
  _t2_pt_res_pt->GetYaxis()->SetTitleOffset(1.4);
 
484
  _t2_pt_res_pt->Draw("AP");
 
485
 
 
486
  _cGraphs->cd(4);
 
487
  _t2_pz_res_pt = new TGraph(_vec_t2_ptMc.size(), &_vec_t2_ptMc[0], &_vec_t2_pz_res[0]);
 
488
  _t2_pz_res_pt->SetName("t2_pz_res_pt");
 
489
  _t2_pz_res_pt->SetTitle("T2 p_{z} res vs. p_{t}^{MC}");
 
490
  _t2_pz_res_pt->GetXaxis()->SetTitle("p_{t}^{MC} (MeV/c)");
 
491
  _t2_pz_res_pt->GetYaxis()->SetTitle("p_{z}^{MC} - p_{z} (MeV/c)");
 
492
  _t2_pz_res_pt->GetYaxis()->SetTitleOffset(1.4);
 
493
  _t2_pz_res_pt->GetYaxis()->SetRangeUser(-200, 200);
 
494
  _t2_pz_res_pt->Draw("AP");
 
495
}
 
496
 
 
497
void TrackerDataAnalyserMomentum::make_resolution_graphs() {
 
498
  _cResolutions = new TCanvas("cResolutions", "Mometum Resolution Graphs");
 
499
  _cResolutions->Divide(2);
 
500
 
 
501
  _cResolutions->cd(1);
 
502
  _t1_pz_resol->SetName("t1_pz_resol");
 
503
  _t1_pz_resol->SetTitle("T1 p_{z} Resolution");
 
504
  _t1_pz_resol->GetXaxis()->SetTitle("p_{t}^{MC} (MeV/c)");
 
505
  _t1_pz_resol->GetYaxis()->SetTitle("p_{z} Resolution (MeV/c)");
 
506
  _t1_pz_resol->Draw("AP");
 
507
 
 
508
  _cResolutions->cd(2);
 
509
  _t2_pz_resol->SetName("t2_pz_resol");
 
510
  _t2_pz_resol->SetTitle("T2 p_{z} Resolution");
 
511
  _t2_pz_resol->GetXaxis()->SetTitle("p_{t}^{MC} (MeV/c)");
 
512
  _t2_pz_resol->GetYaxis()->SetTitle("p_{z} Resolution (MeV/c)");
 
513
  _t2_pz_resol->Draw("AP");
 
514
}
 
515
 
 
516
void TrackerDataAnalyserMomentum::calc_pz_resolution(const int trker, const TCut cut,
 
517
                                                     double &res, double &res_err) {
 
518
  if (_tree) {
 
519
    TCanvas lCanvas;
 
520
    // Create a histogram of the pz residual for the pt_mc interval defined by the input cut
 
521
    if (trker == 0) {
 
522
      _tree->Draw("pz_mc+charge*pz_rec>>h1", cut);
 
523
    } else if (trker == 1) {
 
524
      _tree->Draw("pz_mc-charge*pz_rec>>h1", cut);
 
525
    }
 
526
    // Pull the histogram from the current pad
 
527
    TH1 *h1 = reinterpret_cast<TH1*>(gPad->GetPrimitive("h1"));
 
528
    // Fit a gaussian to the histogram
 
529
    h1->Fit("gaus");
 
530
    TF1 *fit1 = h1->GetFunction("gaus");
 
531
    // Extract the sigma of the gaussian fit (equivalent to the pz resolution for this interval)
 
532
    res = fit1->GetParameter(2);
 
533
    res_err = fit1->GetParError(2);
 
534
    // Scale gets messed up unless histogram is deleted each time
 
535
    delete h1;
 
536
  } else {
 
537
    std::cerr << "Tree pointer invalid" << std::endl;
 
538
  }
 
539
}
 
540
 
 
541
bool TrackerDataAnalyserMomentum::save_graphics(std::string save_type) {
 
542
  if ( (save_type == "eps") || (save_type == "pdf") || (save_type == "png") ) {
 
543
    if (_cResiduals) {
 
544
      std::string save_name = "mom_res_histos.";
 
545
      save_name += save_type;
 
546
      _cResiduals->Print(save_name.c_str());
 
547
    }
 
548
    if (_cGraphs) {
 
549
      std::string save_name = "mom_residual_graphs.";
 
550
      save_name += save_type;
 
551
      _cGraphs->Print(save_name.c_str());
 
552
    }
 
553
    if (_cResolutions) {
 
554
      std::string save_name = "mom_resolution_graphs.";
 
555
      save_name += save_type;
 
556
      _cResolutions->Print(save_name.c_str());
 
557
    }
 
558
  } else {
 
559
    return false;
 
560
  }
 
561
  return true;
 
562
}
 
563
 
 
564
void TrackerDataAnalyserMomentum::save_root() {
 
565
  if (_out_file) {
 
566
    _out_file->cd();
 
567
    if (_tree) _tree->Write();
 
568
 
 
569
    if (_t1_pt_res) _t1_pt_res->Write();
 
570
    if (_t1_pz_res) _t1_pz_res->Write();
 
571
    if (_t1_pz_res_log) _t1_pz_res_log->Write();
 
572
    if (_t1_pt_res_pt) _t1_pt_res_pt->Write();
 
573
    if (_t1_pz_res_pt) _t1_pz_res_pt->Write();
 
574
 
 
575
    if (_t2_pt_res) _t2_pt_res->Write();
 
576
    if (_t2_pz_res) _t2_pz_res->Write();
 
577
    if (_t2_pz_res_log) _t2_pz_res_log->Write();
 
578
    if (_t2_pt_res_pt) _t2_pt_res_pt->Write();
 
579
    if (_t2_pz_res_pt) _t2_pz_res_pt->Write();
 
580
 
 
581
    if (_t1_pz_resol) _t1_pz_resol->Write();
 
582
    if (_t2_pz_resol) _t2_pz_resol->Write();
 
583
 
 
584
    if (_cResiduals) _cResiduals->Write();
 
585
    if (_cGraphs) _cGraphs->Write();
 
586
    if (_cResolutions) _cResolutions->Write();
 
587
 
 
588
    _out_file->Write();
 
589
    _out_file->Close();
 
590
  } else {
 
591
    std::cerr << "Invalid ROOT file pointer provided" << std::endl;
 
592
  }
 
593
}
 
594
 
 
595
TCut TrackerDataAnalyserMomentum::formTCut(const std::string &var, const std::string &op,
 
596
                                           double value) {
 
597
  TCut cut = "";
 
598
  std::stringstream ss1;
 
599
  TString s1;
 
600
  ss1 << var << op << value;
 
601
  ss1 >> s1;
 
602
  cut = s1;
 
603
  return cut;
 
604
}
 
605
 
 
606
} // ~namespace MAUS
 
607