Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakeNewMatcherPlots.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MakeNewMatcherPlots.cxx
1 // ----------------------------------------------------------------------------
2 // 'MakeNewMatcherPlots.cxx'
3 // Derek Anderson
4 // 08.24.2023
5 //
6 // Short macro to make plots from the
7 // output of the ClusMatchTree module.
8 // ----------------------------------------------------------------------------
9 
10 #include <string>
11 #include <vector>
12 #include <utility>
13 #include <iostream>
14 #include "TH1.h"
15 #include "TH2.h"
16 #include "TFile.h"
17 #include "TTree.h"
18 #include "TError.h"
19 
20 using namespace std;
21 
22 
23 
24 void MakeNewMatcherPlots() {
25 
26  // io parameters
27  const string sOutput("newMatcherPlots_oneMatchPerParticle_oddFrac02120.pt10n1evt500pim.d18m1y2024.root");
28  const string sInTrue("input/merged/sPhenixG4_oneMatchPerParticle_newMatcher.pt10num1evt500pim.d4m1y2024.root");
29  const string sInReco("input/merged/sPhenixG4_oneMatchPerParticle_newMatcher.pt10num1evt500pim.d4m1y2024.root");
30  const string sTreeTrue("T");
31  const string sTreeReco("T");
32 
33  // weird track parameters
34  const pair<float, float> oddPtFrac = {0.2, 1.2};
35 
36  // lower verbosity
37  gErrorIgnoreLevel = kError;
38  cout << "\n Staring new matcher plot script..." << endl;
39 
40  // open files
41  TFile* fOutput = new TFile(sOutput.data(), "recreate");
42  TFile* fInTrue = new TFile(sInTrue.data(), "read");
43  TFile* fInReco = new TFile(sInReco.data(), "read");
44  if (!fOutput || !fInTrue || !fInReco) {
45  cerr << "PANIC: couldn't open a file!\n"
46  << " fOutput = " << fOutput << "fInTrue = " << fInTrue << ", fInReco = " << fInReco << "\n"
47  << endl;
48  return;
49  }
50  cout << " Opened files." << endl;
51 
52  // grab input tree
53  TTree* tInTrue = (TTree*) fInTrue -> Get(sTreeTrue.data());
54  TTree* tInReco = (TTree*) fInReco -> Get(sTreeReco.data());
55  if (!tInTrue || !tInReco) {
56  cerr << "PANIC: couldn't grab input tree!\n"
57  << " sTreeTrue = " << sTreeTrue << ", tInTrue = " << tInTrue << "\n"
58  << " sTreeReco = " << sTreeReco << ", tInReco = " << tInReco << "\n"
59  << endl;
60  return;
61  }
62  cout << " Grabbed input trees." << endl;
63 
64  // declare input leaves (ignore cluster branches for now)
65  Int_t tru_event;
66  Int_t tru_nphg4_part;
67  Float_t tru_centrality;
68  Int_t tru_ntrackmatches;
69  Int_t tru_nphg4;
70  Int_t tru_nsvtx;
71  Int_t tru_trackid;
72  Bool_t tru_is_G4track;
73  Bool_t tru_is_Svtrack;
74  Bool_t tru_is_matched;
75  Float_t tru_trkpt;
76  Float_t tru_trketa;
77  Float_t tru_trkphi;
78  Int_t tru_nclus;
79  Int_t tru_nclustpc;
80  Int_t tru_nclusmvtx;
81  Int_t tru_nclusintt;
82  Float_t tru_matchrat;
83  Float_t tru_matchrat_intt;
84  Float_t tru_matchrat_mvtx;
85  Float_t tru_matchrat_tpc;
86 
87  Int_t rec_event;
88  Int_t rec_nphg4_part;
89  Float_t rec_centrality;
90  Int_t rec_ntrackmatches;
91  Int_t rec_nphg4;
92  Int_t rec_nsvtx;
93  Int_t rec_trackid;
94  Bool_t rec_is_G4track;
95  Bool_t rec_is_Svtrack;
96  Bool_t rec_is_matched;
97  Float_t rec_trkpt;
98  Float_t rec_trketa;
99  Float_t rec_trkphi;
100  Int_t rec_nclus;
101  Int_t rec_nclustpc;
102  Int_t rec_nclusmvtx;
103  Int_t rec_nclusintt;
104  Float_t rec_matchrat;
105  Float_t rec_matchrat_intt;
106  Float_t rec_matchrat_mvtx;
107  Float_t rec_matchrat_tpc;
108 
109  // set branch addresses (ignore cluster branches for now)
110  tInTrue -> SetBranchAddress("event", &tru_event);
111  tInTrue -> SetBranchAddress("nphg4_part", &tru_nphg4_part);
112  tInTrue -> SetBranchAddress("centrality", &tru_centrality);
113  tInTrue -> SetBranchAddress("ntrackmatches", &tru_ntrackmatches);
114  tInTrue -> SetBranchAddress("nphg4", &tru_nphg4);
115  tInTrue -> SetBranchAddress("nsvtx", &tru_nsvtx);
116  tInTrue -> SetBranchAddress("trackid", &tru_trackid);
117  tInTrue -> SetBranchAddress("is_G4track", &tru_is_G4track);
118  tInTrue -> SetBranchAddress("is_Svtrack", &tru_is_Svtrack);
119  tInTrue -> SetBranchAddress("is_matched", &tru_is_matched);
120  tInTrue -> SetBranchAddress("trkpt", &tru_trkpt);
121  tInTrue -> SetBranchAddress("trketa", &tru_trketa);
122  tInTrue -> SetBranchAddress("trkphi", &tru_trkphi);
123  tInTrue -> SetBranchAddress("nclus", &tru_nclus);
124  tInTrue -> SetBranchAddress("nclustpc", &tru_nclustpc);
125  tInTrue -> SetBranchAddress("nclusmvtx", &tru_nclusmvtx);
126  tInTrue -> SetBranchAddress("nclusintt", &tru_nclusintt);
127  tInTrue -> SetBranchAddress("matchrat", &tru_matchrat);
128  tInTrue -> SetBranchAddress("matchrat_intt", &tru_matchrat_intt);
129  tInTrue -> SetBranchAddress("matchrat_mvtx", &tru_matchrat_mvtx);
130  tInTrue -> SetBranchAddress("matchrat_tpc", &tru_matchrat_tpc);
131 
132  tInReco -> SetBranchAddress("event", &rec_event);
133  tInReco -> SetBranchAddress("nphg4_part", &rec_nphg4_part);
134  tInReco -> SetBranchAddress("centrality", &rec_centrality);
135  tInReco -> SetBranchAddress("ntrackmatches", &rec_ntrackmatches);
136  tInReco -> SetBranchAddress("nphg4", &rec_nphg4);
137  tInReco -> SetBranchAddress("nsvtx", &rec_nsvtx);
138  tInReco -> SetBranchAddress("trackid", &rec_trackid);
139  tInReco -> SetBranchAddress("is_G4track", &rec_is_G4track);
140  tInReco -> SetBranchAddress("is_Svtrack", &rec_is_Svtrack);
141  tInReco -> SetBranchAddress("is_matched", &rec_is_matched);
142  tInReco -> SetBranchAddress("trkpt", &rec_trkpt);
143  tInReco -> SetBranchAddress("trketa", &rec_trketa);
144  tInReco -> SetBranchAddress("trkphi", &rec_trkphi);
145  tInReco -> SetBranchAddress("nclus", &rec_nclus);
146  tInReco -> SetBranchAddress("nclustpc", &rec_nclustpc);
147  tInReco -> SetBranchAddress("nclusmvtx", &rec_nclusmvtx);
148  tInReco -> SetBranchAddress("nclusintt", &rec_nclusintt);
149  tInReco -> SetBranchAddress("matchrat", &rec_matchrat);
150  tInReco -> SetBranchAddress("matchrat_intt", &rec_matchrat_intt);
151  tInReco -> SetBranchAddress("matchrat_mvtx", &rec_matchrat_mvtx);
152  tInReco -> SetBranchAddress("matchrat_tpc", &rec_matchrat_tpc);
153  cout << " Set input branches." << endl;
154 
155  // histogram indices
156  enum Var {
157  NTot,
158  NIntt,
159  NMvtx,
160  NTpc,
161  RTot,
162  RIntt,
163  RMvtx,
164  RTpc,
165  Phi,
166  Eta,
167  Pt,
168  Frac,
169  Qual,
170  PtErr,
171  EtaErr,
172  PhiErr,
173  PtRes,
174  EtaRes,
175  PhiRes
176  };
177  enum Type {
178  Truth,
179  Track,
180  Weird,
181  Normal
182  };
183  enum Comp {
184  VsTruthPt,
185  VsNumTpc
186  };
187 
188  // output histogram base names
189  const vector<vector<string>> vecHistBase = {
190  {"hTruthNumTot", "hTrackNumTot", "hWeirdNumTot", "hNormNumTot"},
191  {"hTruthNumIntt", "hTrackNumIntt", "hWeirdNumIntt", "hNormNumIntt"},
192  {"hTruthNumMvtx", "hTrackNumMvtx", "hWeirdNumMvtx", "hNormNumMvtx"},
193  {"hTruthNumTpc", "hTrackNumTpc", "hWeirdNumTpc", "hNormNumTpc"},
194  {"hTruthRatTot", "hTrackRatTot", "hWeirdRatTot", "hNormRatTot"},
195  {"hTruthRatIntt", "hTrackRatIntt", "hWeirdRatIntt", "hNormRatIntt"},
196  {"hTruthRatMvtx", "hTrackRatMvtx", "hWeirdRatMvtx", "hNormRatMvtx"},
197  {"hTruthRatTpc", "hTrackRatTpc", "hWeirdRatTpc", "hNormRatTpc"},
198  {"hTruthPhi", "hTrackPhi", "hWeirdPhi", "hNormPhi"},
199  {"hTruthEta", "hTrackEta", "hWeirdEta", "hNormEta"},
200  {"hTruthPt", "hTrackPt", "hWeirdPt", "hNormPt"},
201  {"hTruthFrac", "hTrackFrac", "hWeirdFrac", "hNormFrac"},
202  {"hTruthQual", "hTrackQual", "hWeirdQual", "hNormQual"},
203  {"hTruthPtErr", "hTrackPtErr", "hWeirdPtErr", "hNormPtErr"},
204  {"hTruthEtaErr", "hTrackEtaErr", "hWeirdEtaErr", "hNormEtaErr"},
205  {"hTruthPhiErr", "hTrackPhiErr", "hWeirdPhiErr", "hNormPhiErr"},
206  {"hTruthPtRes", "hTrackPtRes", "hWeirdPtRes", "hNormPtRes"},
207  {"hTruthEtaRes", "hTrackEtaRes", "hWeirdEtaRes", "hNormEtaRes"},
208  {"hTruthPhiRes", "hTrackPhiRes", "hWeirdPhiRes", "hNormPhiRes"}
209  };
210  const size_t nHistRows = vecHistBase.size();
211  const size_t nHistTypes = vecHistBase[0].size();
212 
213  // 2D histogram name modifiers
214  const vector<string> vecVsHistModifiers = {
215  "VsTruthPt",
216  "VsNumTpc"
217  };
218  const size_t nVsHistMods = vecVsHistModifiers.size();
219 
220  // axis titles
221  const string sCount("counts");
222  const vector<string> vecBaseAxisVars = {
223  "N^{tot} = N_{hit}^{mvtx} + N_{hit}^{intt} + N_{clust}^{tpc}",
224  "N_{hit}^{intt}",
225  "N_{hit}^{mvtx}",
226  "N_{clust}^{tpc}",
227  "N_{reco}^{tot} / N_{true}^{tot}",
228  "N_{reco}^{intt} / N_{true}^{intt}",
229  "N_{reco}^{mvtx} / N_{true}^{mvtx}",
230  "N_{reco}^{tpc} / N_{true}^{tpc}",
231  "#varphi",
232  "#eta",
233  "p_{T} [GeV/c]",
234  "p_{T}^{reco} / p_{T}^{true}"
235  };
236  const vector<string> vecVsAxisVars = {
237  "p_{T}^{true} [GeV/c]",
238  "N_{clust}^{tpc}"
239  };
240 
241  // output histogram no. of bins
242  const uint32_t nNumBins = 101;
243  const uint32_t nRatBins = 120;
244  const uint32_t nEtaBins = 80;
245  const uint32_t nPhiBins = 360;
246  const uint32_t nPtBins = 202;
247  const uint32_t nFracBins = 220;
248  const uint32_t nQualBins = ;
249  const uint32_t nResBins = ;
250 
251  // output histogram bin ranges
252  const pair<float, float> xNumBins = {-0.5, 100.5};
253  const pair<float, float> xRatBins = {-0.5, 5.5};
254  const pair<float, float> xEtaBins = {-2., 2.};
255  const pair<float, float> xPhiBins = {-3.15, 3.15};
256  const pair<float, float> xPtBins = {-1., 100.};
257  const pair<float, float> xFracBins = {-0.5, 10.5};
258  const pair<float, float> xQualBins = {};
259  const pair<float, float> xResBins = {};
260 
261  // output histogram base binning definitions
262  vector<tuple<uint32_t, pair<float, float>>> vecBaseHistBins = {
263  make_tuple(nNumBins, xNumBins),
264  make_tuple(nNumBins, xNumBins),
265  make_tuple(nNumBins, xNumBins),
266  make_tuple(nNumBins, xNumBins),
267  make_tuple(nRatBins, xRatBins),
268  make_tuple(nRatBins, xRatBins),
269  make_tuple(nRatBins, xRatBins),
270  make_tuple(nRatBins, xRatBins),
271  make_tuple(nPhiBins, xPhiBins),
272  make_tuple(nEtaBins, xEtaBins),
273  make_tuple(nPtBins, xPtBins),
274  make_tuple(nFracBins, xFracBins),
275  make_tuple(nResBins, xResBins)
276  };
277 
278  // output 2D histogram x-axis binning
279  vector<tuple<uint32_t, pair<float, float>>> vecVsHistBins = {
280  make_tuple(nPtBins, xPtBins),
281  make_tuple(nNumBins, xNumBins)
282  };
283 
284  // make 1D histograms
285  vector<vector<TH1D*>> vecHist1D(nHistRows);
286  for (size_t iHistRow = 0; iHistRow < nHistRows; iHistRow++) {
287  for (const string sHistBase : vecHistBase[iHistRow]) {
288  vecHist1D[iHistRow].push_back(
289  new TH1D(
290  sHistBase.data(),
291  "",
292  get<0>(vecBaseHistBins[iHistRow]),
293  get<1>(vecBaseHistBins[iHistRow]).first,
294  get<1>(vecBaseHistBins[iHistRow]).second
295  )
296  );
297  } // end type loop
298  } // end row loop
299 
300  // make 2D histograms
301  vector<vector<vector<TH2D*>>> vecHist2D(nHistRows, vector<vector<TH2D*>>(nHistTypes));
302  for (size_t iHistRow = 0; iHistRow < nHistRows; iHistRow++) {
303  for (size_t iHistType = 0; iHistType < nHistTypes; iHistType++) {
304  for (size_t iVsHistMod = 0; iVsHistMod < nVsHistMods; iVsHistMod++) {
305  const string sHistName2D = vecHistBase[iHistRow][iHistType] + vecVsHistModifiers[iVsHistMod];
306  vecHist2D[iHistRow][iHistType].push_back(
307  new TH2D(
308  sHistName2D.data(),
309  "",
310  get<0>(vecVsHistBins[iVsHistMod]),
311  get<1>(vecVsHistBins[iVsHistMod]).first,
312  get<1>(vecVsHistBins[iVsHistMod]).second,
313  get<0>(vecBaseHistBins[iHistRow]),
314  get<1>(vecBaseHistBins[iHistRow]).first,
315  get<1>(vecBaseHistBins[iHistRow]).second
316  )
317  );
318  } // end modifier loop
319  } // end type loop
320  } // end row loop
321 
322  // set errors
323  for (auto histRow1D : vecHist1D) {
324  for (auto hist1D : histRow1D) {
325  hist1D -> Sumw2();
326  }
327  }
328 
329  for (auto histRow2D : vecHist2D) {
330  for (auto histType2D : histRow2D) {
331  for (auto hist2D: histType2D) {
332  hist2D -> Sumw2();
333  }
334  }
335  }
336 
337  // grab no. of entries
338  const int64_t nTrueEntries = tInTrue -> GetEntries();
339  const int64_t nRecoEntries = tInReco -> GetEntries();
340  cout << " Beginning truth particle loop: " << nTrueEntries << " to process" << endl;
341 
342  // loop over truth particles
343  int64_t nTrueBytes = 0;
344  for (int64_t iTrueEntry = 0; iTrueEntry < nTrueEntries; iTrueEntry++) {
345 
346  // grab truth particle entry
347  const int64_t trueBytes = tInTrue -> GetEntry(iTrueEntry);
348  if (trueBytes < 0) {
349  cerr << "PANIC: issue with entry " << iTrueEntry << "! Aborting loop!\n" << endl;
350  break;
351  } else {
352  nTrueBytes += trueBytes;
353  }
354 
355  const int64_t iTrueProg = iTrueEntry + 1;
356  if (iTrueProg == nTrueEntries) {
357  cout << " Processing entry " << iTrueProg << "/" << nTrueEntries << "..." << endl;
358  } else {
359  cout << " Processing entry " << iTrueProg << "/" << nTrueEntries << "...\r" << flush;
360  }
361 
362  // select truth particles
363  if (!tru_is_G4track) continue;
364 
365  // fill truth 1D histograms
366  vecHist1D[Var::NTot][Type::Truth] -> Fill(tru_nclus);
367  vecHist1D[Var::NIntt][Type::Truth] -> Fill(tru_nclusintt);
368  vecHist1D[Var::NMvtx][Type::Truth] -> Fill(tru_nclusmvtx);
369  vecHist1D[Var::NTpc][Type::Truth] -> Fill(tru_nclustpc);
370  vecHist1D[Var::RTot][Type::Truth] -> Fill(1.);
371  vecHist1D[Var::RIntt][Type::Truth] -> Fill(1.);
372  vecHist1D[Var::RMvtx][Type::Truth] -> Fill(1.);
373  vecHist1D[Var::RTpc][Type::Truth] -> Fill(1.);
374  vecHist1D[Var::Phi][Type::Truth] -> Fill(tru_trkphi);
375  vecHist1D[Var::Eta][Type::Truth] -> Fill(tru_trketa);
376  vecHist1D[Var::Pt][Type::Truth] -> Fill(tru_trkpt);
377  vecHist1D[Var::Frac][Type::Truth] -> Fill(1.);
378 
379  // fill truth 2D histograms
380  vecHist2D[Var::NTot][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_nclus);
381  vecHist2D[Var::NIntt][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_nclusintt);
382  vecHist2D[Var::NMvtx][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_nclusmvtx);
383  vecHist2D[Var::NTpc][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_nclustpc);
384  vecHist2D[Var::RTot][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, 1.);
385  vecHist2D[Var::RIntt][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, 1.);
386  vecHist2D[Var::RMvtx][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, 1.);
387  vecHist2D[Var::RTpc][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, 1.);
388  vecHist2D[Var::Phi][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_trkphi);
389  vecHist2D[Var::Eta][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_trketa);
390  vecHist2D[Var::Pt][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, tru_trkpt);
391  vecHist2D[Var::Frac][Type::Truth][Comp::VsTruthPt] -> Fill(tru_trkpt, 1.);
392 
393  vecHist2D[Var::NTot][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_nclus);
394  vecHist2D[Var::NIntt][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_nclusintt);
395  vecHist2D[Var::NMvtx][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_nclusmvtx);
396  vecHist2D[Var::NTpc][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_nclustpc);
397  vecHist2D[Var::RTot][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, 1.);
398  vecHist2D[Var::RIntt][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, 1.);
399  vecHist2D[Var::RMvtx][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, 1.);
400  vecHist2D[Var::RTpc][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, 1.);
401  vecHist2D[Var::Phi][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_trkphi);
402  vecHist2D[Var::Eta][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_trketa);
403  vecHist2D[Var::Pt][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, tru_trkpt);
404  vecHist2D[Var::Frac][Type::Truth][Comp::VsNumTpc] -> Fill(tru_nclustpc, 1.);
405  } // end truth particle loop
406 
407  // announce next entry loop
408  cout << " Finished truth particle loop.\n"
409  << " Beginning reconstructed track loop: " << nRecoEntries << " to process"
410  << endl;
411 
412  // loop over reco tracks
413  // TODO identify matched truth particles
414  int64_t nRecoBytes = 0;
415  for (int64_t iRecoEntry = 0; iRecoEntry < nRecoEntries; iRecoEntry++) {
416 
417  // grab reco track entry
418  const int64_t recoBytes = tInReco -> GetEntry(iRecoEntry);
419  if (recoBytes < 0) {
420  cerr << "PANIC: issue with entry " << iRecoEntry << "! Aborting loop!\n" << endl;
421  break;
422  } else {
423  nRecoBytes += recoBytes;
424  }
425 
426  const int64_t iRecoProg = iRecoEntry + 1;
427  if (iRecoProg == nRecoEntries) {
428  cout << " Processing entry " << iRecoProg << "/" << nRecoEntries << "..." << endl;
429  } else {
430  cout << " Processing entry " << iRecoProg << "/" << nRecoEntries << "...\r" << flush;
431  }
432 
433  // select only tracks matched to truth particle
434  if (!rec_is_matched || !rec_is_Svtrack) continue;
435 
436  // fill all matched reco 1D histograms
437  // FIXME actually calculate pt fraction
438  vecHist1D[Var::NTot][Type::Track] -> Fill(rec_nclus);
439  vecHist1D[Var::NIntt][Type::Track] -> Fill(rec_nclusintt);
440  vecHist1D[Var::NMvtx][Type::Track] -> Fill(rec_nclusmvtx);
441  vecHist1D[Var::NTpc][Type::Track] -> Fill(rec_nclustpc);
442  vecHist1D[Var::RTot][Type::Track] -> Fill(rec_matchrat);
443  vecHist1D[Var::RIntt][Type::Track] -> Fill(rec_matchrat_intt);
444  vecHist1D[Var::RMvtx][Type::Track] -> Fill(rec_matchrat_mvtx);
445  vecHist1D[Var::RTpc][Type::Track] -> Fill(rec_matchrat_tpc);
446  vecHist1D[Var::Phi][Type::Track] -> Fill(rec_trkphi);
447  vecHist1D[Var::Eta][Type::Track] -> Fill(rec_trketa);
448  vecHist1D[Var::Pt][Type::Track] -> Fill(rec_trkpt);
449  vecHist1D[Var::Frac][Type::Track] -> Fill(1.);
450 
451  // fill all matched reco 2D histograms
452  // FIXME use actual truth pt & ntpc of matched particle
453  vecHist2D[Var::NTot][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclus);
454  vecHist2D[Var::NIntt][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclusintt);
455  vecHist2D[Var::NMvtx][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclusmvtx);
456  vecHist2D[Var::NTpc][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclustpc);
457  vecHist2D[Var::RTot][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat);
458  vecHist2D[Var::RIntt][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_intt);
459  vecHist2D[Var::RMvtx][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_mvtx);
460  vecHist2D[Var::RTpc][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_tpc);
461  vecHist2D[Var::Phi][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trkphi);
462  vecHist2D[Var::Eta][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trketa);
463  vecHist2D[Var::Pt][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trkpt);
464  vecHist2D[Var::Frac][Type::Track][Comp::VsTruthPt] -> Fill(rec_trkpt, 1.);
465 
466  vecHist2D[Var::NTot][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclus);
467  vecHist2D[Var::NIntt][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclusintt);
468  vecHist2D[Var::NMvtx][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclusmvtx);
469  vecHist2D[Var::NTpc][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclustpc);
470  vecHist2D[Var::RTot][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat);
471  vecHist2D[Var::RIntt][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_intt);
472  vecHist2D[Var::RMvtx][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_mvtx);
473  vecHist2D[Var::RTpc][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_tpc);
474  vecHist2D[Var::Phi][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trkphi);
475  vecHist2D[Var::Eta][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trketa);
476  vecHist2D[Var::Pt][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trkpt);
477  vecHist2D[Var::Frac][Type::Track][Comp::VsNumTpc] -> Fill(rec_nclustpc, 1.);
478 
479  // fill weird and normal matched reco 1D histograms
480  // FIXME actually check if is a weird track
481  const bool isWeirdTrack = true;
482  if (isWeirdTrack) {
483  vecHist1D[Var::NTot][Type::Weird] -> Fill(rec_nclus);
484  vecHist1D[Var::NIntt][Type::Weird] -> Fill(rec_nclusintt);
485  vecHist1D[Var::NMvtx][Type::Weird] -> Fill(rec_nclusmvtx);
486  vecHist1D[Var::NTpc][Type::Weird] -> Fill(rec_nclustpc);
487  vecHist1D[Var::RTot][Type::Weird] -> Fill(rec_matchrat);
488  vecHist1D[Var::RIntt][Type::Weird] -> Fill(rec_matchrat_intt);
489  vecHist1D[Var::RMvtx][Type::Weird] -> Fill(rec_matchrat_mvtx);
490  vecHist1D[Var::RTpc][Type::Weird] -> Fill(rec_matchrat_tpc);
491  vecHist1D[Var::Phi][Type::Weird] -> Fill(rec_trkphi);
492  vecHist1D[Var::Eta][Type::Weird] -> Fill(rec_trketa);
493  vecHist1D[Var::Pt][Type::Weird] -> Fill(rec_trkpt);
494  vecHist1D[Var::Frac][Type::Weird] -> Fill(1.);
495 
496  vecHist2D[Var::NTot][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclus);
497  vecHist2D[Var::NIntt][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclusintt);
498  vecHist2D[Var::NMvtx][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclusmvtx);
499  vecHist2D[Var::NTpc][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclustpc);
500  vecHist2D[Var::RTot][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat);
501  vecHist2D[Var::RIntt][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_intt);
502  vecHist2D[Var::RMvtx][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_mvtx);
503  vecHist2D[Var::RTpc][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_tpc);
504  vecHist2D[Var::Phi][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trkphi);
505  vecHist2D[Var::Eta][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trketa);
506  vecHist2D[Var::Pt][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trkpt);
507  vecHist2D[Var::Frac][Type::Weird][Comp::VsTruthPt] -> Fill(rec_trkpt, 1.);
508 
509  vecHist2D[Var::NTot][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclus);
510  vecHist2D[Var::NIntt][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclusintt);
511  vecHist2D[Var::NMvtx][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclusmvtx);
512  vecHist2D[Var::NTpc][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclustpc);
513  vecHist2D[Var::RTot][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat);
514  vecHist2D[Var::RIntt][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_intt);
515  vecHist2D[Var::RMvtx][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_mvtx);
516  vecHist2D[Var::RTpc][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_tpc);
517  vecHist2D[Var::Phi][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trkphi);
518  vecHist2D[Var::Eta][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trketa);
519  vecHist2D[Var::Pt][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trkpt);
520  vecHist2D[Var::Frac][Type::Weird][Comp::VsNumTpc] -> Fill(rec_nclustpc, 1.);
521  } else {
522  vecHist1D[Var::NTot][Type::Normal] -> Fill(rec_nclus);
523  vecHist1D[Var::NIntt][Type::Normal] -> Fill(rec_nclusintt);
524  vecHist1D[Var::NMvtx][Type::Normal] -> Fill(rec_nclusmvtx);
525  vecHist1D[Var::NTpc][Type::Normal] -> Fill(rec_nclustpc);
526  vecHist1D[Var::RTot][Type::Normal] -> Fill(rec_matchrat);
527  vecHist1D[Var::RIntt][Type::Normal] -> Fill(rec_matchrat_intt);
528  vecHist1D[Var::RMvtx][Type::Normal] -> Fill(rec_matchrat_mvtx);
529  vecHist1D[Var::RTpc][Type::Normal] -> Fill(rec_matchrat_tpc);
530  vecHist1D[Var::Phi][Type::Normal] -> Fill(rec_trkphi);
531  vecHist1D[Var::Eta][Type::Normal] -> Fill(rec_trketa);
532  vecHist1D[Var::Pt][Type::Normal] -> Fill(rec_trkpt);
533  vecHist1D[Var::Frac][Type::Normal] -> Fill(1.);
534 
535  vecHist2D[Var::NTot][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclus);
536  vecHist2D[Var::NIntt][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclusintt);
537  vecHist2D[Var::NMvtx][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclusmvtx);
538  vecHist2D[Var::NTpc][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_nclustpc);
539  vecHist2D[Var::RTot][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat);
540  vecHist2D[Var::RIntt][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_intt);
541  vecHist2D[Var::RMvtx][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_mvtx);
542  vecHist2D[Var::RTpc][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_matchrat_tpc);
543  vecHist2D[Var::Phi][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trkphi);
544  vecHist2D[Var::Eta][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trketa);
545  vecHist2D[Var::Pt][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, rec_trkpt);
546  vecHist2D[Var::Frac][Type::Normal][Comp::VsTruthPt] -> Fill(rec_trkpt, 1.);
547 
548  vecHist2D[Var::NTot][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclus);
549  vecHist2D[Var::NIntt][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclusintt);
550  vecHist2D[Var::NMvtx][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclusmvtx);
551  vecHist2D[Var::NTpc][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_nclustpc);
552  vecHist2D[Var::RTot][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat);
553  vecHist2D[Var::RIntt][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_intt);
554  vecHist2D[Var::RMvtx][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_mvtx);
555  vecHist2D[Var::RTpc][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_matchrat_tpc);
556  vecHist2D[Var::Phi][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trkphi);
557  vecHist2D[Var::Eta][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trketa);
558  vecHist2D[Var::Pt][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, rec_trkpt);
559  vecHist2D[Var::Frac][Type::Normal][Comp::VsNumTpc] -> Fill(rec_nclustpc, 1.);
560  }
561  } // end reco track loop
562  cout << " Finished reconstructed track loop." << endl;
563 
564  // set axis titles
565  size_t iVar = 0;
566  for (auto histRow1D : vecHist1D) {
567  for (auto hist1D : histRow1D) {
568  hist1D -> GetXaxis() -> SetTitle(vecBaseAxisVars.at(iVar).data());
569  hist1D -> GetYaxis() -> SetTitle(sCount.data());
570  }
571  ++iVar;
572  }
573  iVar = 0;
574 
575  size_t iComp = 0;
576  for (auto histRow2D : vecHist2D) {
577  for (auto histType2D : histRow2D) {
578  iComp = 0;
579  for (auto hist2D : histType2D) {
580  hist2D -> GetXaxis() -> SetTitle(vecVsAxisVars.at(iComp).data());
581  hist2D -> GetYaxis() -> SetTitle(vecBaseAxisVars.at(iVar).data());
582  hist2D -> GetZaxis() -> SetTitle(sCount.data());
583  ++iComp;
584  }
585  }
586  ++iVar;
587  }
588  cout << " Set axis titles." << endl;
589 
590  // save histograms
591  fOutput -> cd();
592  for (auto histRow1D : vecHist1D) {
593  for (auto hist1D : histRow1D) {
594  hist1D -> Write();
595  }
596  }
597  for (auto histRow2D : vecHist2D) {
598  for (auto histType2D : histRow2D) {
599  for (auto hist2D: histType2D) {
600  hist2D -> Write();
601  }
602  }
603  }
604  cout << " Saved histograms." << endl;
605 
606  // close files
607  fOutput -> cd();
608  fOutput -> Close();
609  fInTrue -> cd();
610  fInTrue -> Close();
611  fInReco -> cd();
612  fInReco -> Close();
613 
614  // exit
615  cout << " Finished new matcher plot script!\n" << endl;
616  return;
617 
618 }
619 
620 // end ------------------------------------------------------------------------