Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
track_calorimeter_matching.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file track_calorimeter_matching.C
1 #include <TChain.h>
2 #include <TCanvas.h>
3 #include <TH1.h>
4 #include <TH1F.h>
5 #include <TH2D.h>
6 #include <TF1.h>
7 #include <TFile.h>
8 #include <TGraph.h>
9 #include <TStyle.h>
10 #include <TROOT.h>
11 #include <TLegend.h>
12 #include <TLine.h>
13 #include <TLatex.h>
14 #include <TRandom1.h>
15 #include <TPolyLine.h>
16 #include <iostream>
17 #include <fstream>
18 #include <TMath.h>
19 #include <TLorentzVector.h>
20 #include <TEfficiency.h>
21 #include <TGraphAsymmErrors.h>
22 
23 void track_calorimeter_matching()
24 {
25  int verbose = 2;
26 
27  int n_maps_layer = 3;
28  int n_intt_layer = 4;
29  int n_tpc_layer = 40;
30 
31  int nlayers = n_maps_layer+n_intt_layer+n_tpc_layer; // maximum number of tracking layers for tpc+intt+maps
32  static const int nmissed = 30; // maximum number of missed layers to allow for a track
33 
34  int embed_flag = 1;
35 
36  double quality_cut = 3.0;
37  double dca_cut = 0.1; // cm
38 
39  TH2D * hEoverP_cemc[2];
40  hEoverP_cemc[0] = new TH2D("hEoverP0_cemc","",200, 0, 20, 200, 0, 2.0);
41  hEoverP_cemc[0]->SetMarkerStyle(20);
42  hEoverP_cemc[0]->SetMarkerSize(0.5);
43  hEoverP_cemc[1] = new TH2D("hEoverP1_cemc","",200, 0, 20, 200, 0, 2.0);
44  hEoverP_cemc[1]->SetMarkerStyle(20);
45  hEoverP_cemc[1]->SetMarkerSize(0.5);
46 
47 
48  int n_embed_rtrack = 0;
49  int naccept = 0;
50 
51  // The condor job output files
52  for(int i=0;i<1000;i++)
53  //for(int i=0;i<2;i++)
54  {
55  // Open the evaluator output file
56 
57  TChain* ntp_track = new TChain("ntp_track","reco tracks");
58  TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
59  TChain* ntp_vertex = new TChain("ntp_vertex","events");
60  TChain *ntp_cluster = new TChain("ntp_cluster","clusters");
61 
62  // This include file contains the definitions of the ntuple variables, and the chain definitions
63 #include "ntuple_variables.C"
64 
65  char name[500];
66  sprintf(name,"/sphenix/user/frawley/latest/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
67 
68  cout << "Adding file number " << i << " with name " << name << endl;
69 
70  ntp_vertex->Add(name);
71  ntp_track->Add(name);
72  ntp_gtrack->Add(name);
73 
74  // skip this file if there are no tracks
75  if(!ntp_gtrack->GetEntries())
76  continue;
77 
78  if(verbose> 0)
79  cout << " ntp_vertex entries: " << ntp_vertex->GetEntries()
80  << " ntp_gtrack entries: " << ntp_gtrack->GetEntries()
81  << " ntp_track entries: " << ntp_track->GetEntries()
82  << endl;
83 
84  //==================================
85  // Begin event loop
86  //==================================
87 
88  // These keep track of the starting position in the track ntuples for each event
89  int nr = 0;
90  //int ng = 0;
91 
92  for(int iev=0;iev<ntp_vertex->GetEntries();iev++)
93  {
94  if(verbose) cout << " iev = " << iev << " nr " << nr << endl;
95 
96  int recoget = ntp_vertex->GetEntry(iev);
97  if(!recoget)
98  {
99  cout << "Failed to get ntp_vertex entry " << iev << endl;
100  exit(1);
101  }
102  if(iev%1 == 0)
103  if(verbose)
104  cout << "Get event " << iev << " with vertex x " << evx
105  << " vertex y " << evy << " vertex z " << evz
106  << " ngtracks " << ngtracks << " ntracks " << ntracks
107  << endl;
108 
109  // loop over tracks
110 
111  for(int ir=nr;ir<nr+ntracks;ir++)
112  {
113  int recoget = ntp_track->GetEntry(ir);
114  if(!recoget)
115  {
116  if(verbose) cout << "Failed to get ntp_track entry " << ir << endl;
117  break;
118  }
119  if(verbose > 0) cout << " ir = " << ir << " rtrackid " << rtrackid << " revent " << revent
120  << " rquality " << rquality << " rgtrackid = " << rgtrackid << " rgnhits " << rgnhits
121  << " rnhits " << rnhits << " rnfromtruth " << rnfromtruth << endl;
122 
123  // Make some overall reconstructed track cuts
124  if(rgnhits < nlayers-nmissed)
125  {
126  //skip tracks that do not pass through enough layers in truth, same cut made on truth tracks earlier
127  cout << " --- skip because rgnhits too small " << endl;
128  continue;
129  }
130 
131  if(rquality > quality_cut)
132  {
133  if(verbose > 0) cout << " --- failed quality cut - rejected " << endl;
134  continue;
135  }
136 
137  if(fabs(rvz - evz) > 0.1)
138  {
139  if(verbose) cout << "skip because failed z vertex cut with rvz = " << rvz << " evz = " << evz << endl;
140  continue;
141  }
142 
143  if( isnan(rdca2d) )
144  {
145  if(verbose) cout << "skip because dca2d is nan" << endl;
146  continue;
147  }
148 
149  if(rdca2d > 0.1)
150  {
151  if(verbose) cout << "skip because failed dca2d cut " << endl;
152  continue;
153  }
154 
155  //double rgpT = sqrt(rgpx*rgpx+rgpy*rgpy);
156  //double rpT = sqrt(rpx*rpx+rpy*rpy);
157 
158  if(rgembed == embed_flag)
159  {
160  n_embed_rtrack++;
161  naccept++;
162 
163  double reta = asinh(rpz/sqrt(rpx*rpx+rpy*rpy));
164  double rmom = sqrt(rpx*rpx+rpy*rpy+rpz*rpz);
165  double EoverP_cemc = cemc_e / rmom;
166  cout << " reta " << reta << " rmom " << rmom << " EoverP " << EoverP_cemc << endl;
167 
168  if(rgflavor ==11)
169  hEoverP_cemc[0]->Fill(rmom, EoverP_cemc);
170  else
171  hEoverP_cemc[1]->Fill(rmom, EoverP_cemc);
172 
173  }
174 
175  } // end loop over reco'd tracks
176  if(verbose > 0) cout << " Done with tracks loop: n_embed_rtrack = " << n_embed_rtrack << " naccept = " << naccept << endl;
177  //====================================================
178 
179  // set the starting positions in the track ntuples for the next event
180  nr += ntracks;
181  //ng += ngtracks;
182  } // end loop over events
183 
184 
185  delete ntp_gtrack;
186  delete ntp_track;
187  delete ntp_cluster;
188  delete ntp_vertex;
189 
190  } // end loop over files
191 
192  TCanvas *cep_cemc = new TCanvas("cep_cemc","",5,5,800,500);
193  cep_cemc->Divide(2,1);
194 
195  cep_cemc->cd(1);
196  hEoverP_cemc[0]->Draw();
197 
198  cep_cemc->cd(2);
199  hEoverP_cemc[1]->Draw();
200 
201  TCanvas *cproj_cemc = new TCanvas("cproj_cemc","",20,20,1400,800);
202  cproj_cemc->Divide(3,2);
203 
204  double dp = 2.0;
205  int np = 12.0 / dp;
206  for(int ip=0;ip<np;ip++)
207  {
208  cproj_cemc->cd(ip+1);
209 
210  double plo = ip * dp;
211  double phi = ip * dp + dp;
212  int binlo = hEoverP_cemc[0]->GetXaxis()->FindBin(plo);
213  int binhi = hEoverP_cemc[0]->GetXaxis()->FindBin(phi);
214 
215  char label[500];
216  sprintf(label,"p = %.1f - %.1f GeV/c",plo, phi);
217  TH1D *h = new TH1D("h",label, 200, 0 , 2.0);
218  hEoverP_cemc[0]->ProjectionY("h",binlo,binhi);
219  h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
220  h->GetXaxis()->SetTitle("E/p for CEMC");
221  h->GetXaxis()->SetTitleOffset(1.0);
222 
223  h->DrawCopy();
224 
225  delete h;
226  }
227 
228 
229 }