Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExploreTTrees.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ExploreTTrees.C
1 #include <iostream>
2 #include <string>
3 #include <vector>
4 #include <TFile.h>
5 #include <TTree.h>
6 #include <TChain.h>
7 
8 // Fixed size dimensions of array or collections stored in the TTree if any.
9 
10 // Declaration of leaf types
11 Int_t event;
12 Float_t eta;
13 Float_t phi;
14 Float_t e;
15 Float_t pt;
16 Float_t vx;
17 Float_t vy;
18 Float_t vz;
19 vector<int> id;
20 
21 int nTow_in;
22 vector<float> eta_in;
23 vector<float> phi_in;
24 vector<float> e_in;
25 vector<int> ieta_in;
26 vector<int> iphi_in;
28 vector<float> eta_out;
29 vector<float> phi_out;
30 vector<float> e_out;
31 vector<int> ieta_out;
32 vector<int> iphi_out;
34 vector<float> eta_emc;
35 vector<float> phi_emc;
36 vector<float> e_emc;
37 vector<int> ieta_emc;
38 vector<int> iphi_emc;
39 
40 TTree *fChain;
41 Int_t fCurrent;
42 
43 
44 const int ncal = 3;
45 string cal_name[] = { "iHCal", "oHCal", "EMCal" };
46 string cal_tag[] = { "_in", "_out", "_emc" };
47 TString inFileName = "HCalJetPhiShift_10k.root";
48 //TString inFileName = "HCalJetPhiShift_negativePion_magnetOn_15k.root";
49 int n_pt_bins[ncal] = {80,100,100};
50 double pt_bin_lo[ncal] = {0.,0.,0.};
51 double pt_bin_hi[ncal] = {4.,50.,20.};
52 double eta_max[ncal] = {1.1,1.1,1.152};
53 int n_eta_bins[ncal] = {24,24,96};
54 int n_phi_bins[ncal] = {64,64,256};
55 
57 vector<float> *calo_eta[ncal] = {&eta_in, &eta_out, &eta_emc};
58 vector<float> *calo_phi[ncal] = {&phi_in, &phi_out, &phi_emc};
59 vector<float> *calo_e[ncal] = {&e_in, &e_out, &e_emc};
60 vector<int> *calo_ieta[ncal] = {&ieta_in, &ieta_out, &ieta_emc};
61 vector<int> *calo_iphi[ncal] = {&iphi_in, &iphi_out, &iphi_emc};
62 
63 double delta_phi(double phi1, double phi2){
64  double dPhi = phi1 - phi2;
65  while (dPhi>=M_PI) { dPhi -= 2.*M_PI; }
66  while (dPhi<=-M_PI) { dPhi += 2.*M_PI; }
67  return dPhi;
68 };
69 
70 double dR(double phi1, double phi2, double eta1, double eta2){
71  double dphi_sqrd = delta_phi( phi1, phi2)*delta_phi( phi1, phi2);
72  double deta_sqrd = (eta1-eta2)*(eta1-eta2);
73  return sqrt(dphi_sqrd+deta_sqrd);
74 };
75 
76 double phi_in_range(double phi){
77  while (phi>=M_PI) { phi -= 2.*M_PI; }
78  while (phi<=-M_PI) { phi += 2.*M_PI; }
79  return phi;
80 };
81 
82 bool tower_in_3x3(float calo_eta, float calo_phi, float eta, float phi){
83  double min_eta = eta - 0.024 - 0.012;
84  double max_eta = eta + 0.024 + 0.012;
85  if ( calo_eta<min_eta || calo_eta>max_eta ) { return false; }
86  double min_phi = phi - (M_PI/128.) - (M_PI/256.);
87  double max_phi = phi + (M_PI/128.) + (M_PI/256.);
88  if ( calo_phi<min_phi || calo_phi>max_phi ) { return false; }
89  if ( calo_eta>=min_eta && calo_eta<=max_eta && calo_phi>=min_phi && calo_phi<=max_phi ) { return true; }
90  return false;
91 };
92 
93 TH1D *hGeantPhi = new TH1D("hGeantPhi",";truth #phi",64,-M_PI,M_PI);
94 TH2D *hEMCal3x3 = new TH2D("hEMCal3x3","pion p_{T};E_{EMCal}^{3x3} about max-E tower [GeV]",30,0.,30.,120,0.,60.);
95 TH2D *hEMCal3x3_FINE = new TH2D("hEMCal3x3_FINE","pion p_{T};E_{EMCal}^{3x3} about max-E tower [GeV]",30,0.,30.,200.,0.,2.);
96 TH3D *hE_inner_vs_outer = new TH3D("hE_inner_vs_outer",";pion p_{T};total E deposited in iHCal;total E deposited in oHCal",60,0.,30.,100, 0., 10.,120,0.,60.);
97 TH3D *hDeltaPhi_fraction_HcalOverAll = new TH3D("hDeltaPhi_fraction_HcalOverAll",";pion p_{T};E_{HCals}/E_{all calos};oHCal #Delta#phi",60,0.,30.,100,0.,1.,160,-0.4,0.4);
98 TH3D *hDeltaPhi_fraction_oHcalOverHcals = new TH3D("hDeltaPhi_fraction_oHcalOverHcals",";pion p_{T};E_{oHCal}/E_{HCals};oHCal #Delta#phi",60,0.,30.,100,0.,1.,160,-0.4,0.4);
99 
102 
104 
106 
107  for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
108  hE_weighted_eta_phi[i_cal] = new TH3D(Form("hE_weighted_eta_phi_%s",cal_name[i_cal].c_str()),";pion p_{T};calo #eta;calo #phi",60,0.,30.,n_eta_bins[i_cal],-eta_max[i_cal],eta_max[i_cal],n_phi_bins[i_cal],-M_PI,M_PI);
109  h_eta_phi[i_cal] = new TH3D(Form("h_eta_phi_%s",cal_name[i_cal].c_str()),";pion p_{T};calo #eta;calo #phi",60,0.,30.,n_eta_bins[i_cal],-eta_max[i_cal],eta_max[i_cal],n_phi_bins[i_cal],-M_PI,M_PI);
110  hPhi2D[i_cal] = new TH2D(Form("hPhi2D_%s",cal_name[i_cal].c_str()),";truth #phi;tower #phi",n_phi_bins[i_cal],-M_PI,M_PI,n_phi_bins[i_cal],-M_PI,M_PI);
111  hDeltaPhi_pt[i_cal] = new TH2D(Form("hDeltaPhi_pt_%s",cal_name[i_cal].c_str()),";pion p_{T};tower #phi - truth #phi",120,0.,60.,320,-4.,4.);
112  hDeltaPhi_E[i_cal] = new TH2D(Form("hDeltaPhi_E_%s",cal_name[i_cal].c_str()),";tower #phi - truth #phi;tower E_{T}",320,-4.,4.,n_pt_bins[i_cal], pt_bin_lo[i_cal], pt_bin_hi[i_cal]);
113  hDeltaPhi_fraction[i_cal] = new TH2D(Form("hDeltaPhi_fraction_%s",cal_name[i_cal].c_str()),";tower #phi - truth #phi;fraction of total calo E in layer",320,-4.,4.,100,0,1.);
114  hDeltaPhi_iPhi[i_cal] = new TH2D(Form("hDeltaPhi_iPhi_%s",cal_name[i_cal].c_str()),";tower #phi - truth #phi;tower iphi",320,-4.,4.,256,0.,256.);
115  hDeltaPhi_eta[i_cal] = new TH2D(Form("hDeltaPhi_eta_%s",cal_name[i_cal].c_str()),";pion #eta;tower #phi - truth #phi",n_phi_bins[i_cal],-eta_max[i_cal],eta_max[i_cal],320,-4.,4.);
116  hTowerEt_pionEt[i_cal] = new TH2D(Form("hTowerEt_pionEt_%s",cal_name[i_cal].c_str()),";tower E_{T};pion E_{T}",n_pt_bins[i_cal], pt_bin_lo[i_cal], pt_bin_hi[i_cal],120,0.,60.);
117  hCaloEnergy_pionPt[i_cal] = new TH2D(Form("hCaloEnergy_pionPt_%s",cal_name[i_cal].c_str()),";pion p_{T};total E deposited in calo",60,0.,30.,n_pt_bins[i_cal], pt_bin_lo[i_cal], pt_bin_hi[i_cal]);
118  hEnergy_fraction[i_cal] = new TH2D(Form("hEnergy_fraction_%s",cal_name[i_cal].c_str()),";pion E;fraction of pion E in calo",120,0.,60.,100,0.,10.);
119  hDeltaPhi[i_cal] = new TH1D(Form("hDeltaPhi_%s",cal_name[i_cal].c_str()),";tower #phi - truth #phi",128,-2.*M_PI,2.*M_PI);
120  hCaloPhi[i_cal] = new TH1D(Form("hCaloPhi_%s",cal_name[i_cal].c_str()),";tower #phi",n_phi_bins[i_cal],-M_PI,M_PI);
121  hTowerEt[i_cal] = new TH1D(Form("hTowerEt_%s",cal_name[i_cal].c_str()),";tower E_{T}",n_pt_bins[i_cal], pt_bin_lo[i_cal], pt_bin_hi[i_cal]);
122  }
123 
124  TFile *f = new TFile(inFileName,"READ");
125  TTree *fChain = (TTree*)f->Get("T");
126 
127  fChain->SetBranchAddress("event", &event);
128  fChain->SetBranchAddress("eta", &eta);
129  fChain->SetBranchAddress("phi", &phi);
130  fChain->SetBranchAddress("e", &e);
131  fChain->SetBranchAddress("pt", &pt);
132  fChain->SetBranchAddress("vx", &vx);
133  fChain->SetBranchAddress("vy", &vy);
134  fChain->SetBranchAddress("vz", &vz);
135 
136  for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
137  fChain->SetBranchAddress(Form("nTow%s",cal_tag[i_cal].c_str()), &nTowers[i_cal]);
138  fChain->SetBranchAddress(Form("eta%s",cal_tag[i_cal].c_str()), &calo_eta[i_cal]);
139  fChain->SetBranchAddress(Form("phi%s",cal_tag[i_cal].c_str()), &calo_phi[i_cal]);
140  fChain->SetBranchAddress(Form("e%s",cal_tag[i_cal].c_str()), &calo_e[i_cal]);
141  fChain->SetBranchAddress(Form("ieta%s",cal_tag[i_cal].c_str()), &calo_ieta[i_cal]);
142  fChain->SetBranchAddress(Form("iphi%s",cal_tag[i_cal].c_str()), &calo_iphi[i_cal]);
143  }
144 
145  Long64_t nentries = fChain->GetEntriesFast();
146 
147  for (Long64_t jentry=0; jentry<nentries;jentry++) { // LOOP OVER EVENTS
148  Long64_t ientry = fChain->GetTreeNumber();
149  if (ientry < 0) break;
150 
151  fChain->GetEntry(jentry);
152 
153  hGeantPhi->Fill(phi_in_range(phi));
154 
155  float total_E[ncal] = {0.,0.,0.};
156  float max_out_phi;
157  // float max_tow_phi[ncal];// = {0.,0.,0.};
158 
159  float closestOuter_tow_dR = 9999.;
160  float closestOuter_tow_eta;
161  float closestOuter_tow_phi;
162 
163  float energyIn3x3 = 0;
164 
165  for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
166  float max_tow_e = 0.;
167  float max_tow_eta;
168  int max_tow_iphi;
169  float max_tow_phi;
170 
171  for (int itow=0; itow<calo_e[i_cal]->size(); ++itow) { // LOOP OVER CALORIMETER TOWERS
172 
173  if (i_cal==2 && tower_in_3x3(calo_eta[i_cal]->at(itow), calo_phi[i_cal]->at(itow), eta, phi) ){
174  energyIn3x3 += calo_e[i_cal]->at(itow);
175  }
176 
177 // if (calo_e[i_cal]->at(itow)<=0.) continue;
178 // if (i_cal==2 && calo_e[i_cal]->at(itow)<0.08) { continue; }
179 
180  total_E[i_cal] += calo_e[i_cal]->at(itow);
181  if (calo_e[i_cal]->at(itow)>max_tow_e){
182  max_tow_e = calo_e[i_cal]->at(itow);
183  max_tow_eta = calo_eta[i_cal]->at(itow);
184  max_tow_phi = calo_phi[i_cal]->at(itow);
185  max_tow_iphi = calo_iphi[i_cal]->at(itow);
186  }
187  if( i_cal==1 ){ // MATCH PION TO CLOSEST OHCAL TOWER
188  float delta_R = dR(phi, calo_phi[i_cal]->at(itow), eta, calo_eta[i_cal]->at(itow));
189  if (delta_R<closestOuter_tow_dR){
190  closestOuter_tow_dR = delta_R;
191  closestOuter_tow_eta = calo_eta[i_cal]->at(itow);
192  closestOuter_tow_phi = calo_phi[i_cal]->at(itow);
193  }
194  }
195  hCaloEnergy_pionPt[i_cal]->Fill(pt,calo_e[i_cal]->at(itow));
196  h_eta_phi[i_cal]->Fill(pt,calo_eta[i_cal]->at(itow),phi_in_range(calo_phi[i_cal]->at(itow)));
198  } // end tower loop
199 
200  if (i_cal==2) {
201  hEMCal3x3->Fill(pt,energyIn3x3);
202  hEMCal3x3_FINE->Fill(pt,energyIn3x3);
203  }
204 
205 // cout<<cal_tag[i_cal]<<"\t"<<max_tow_eta<<"\t"<<eta<<"\t"<<delta_phi(max_tow_phi,phi)<<endl;
206 
207  if (i_cal==1) {max_out_phi=max_tow_phi;}
208 
209  double e_fraction = (double) total_E[i_cal]/e;
210  hEnergy_fraction[i_cal]->Fill(e,e_fraction);
211  hPhi2D[i_cal]->Fill(phi_in_range(phi),max_tow_phi);
212  hDeltaPhi[i_cal]->Fill(delta_phi(max_tow_phi,phi));
213  hDeltaPhi_pt[i_cal]->Fill(pt,delta_phi(max_tow_phi,phi));
214  hDeltaPhi_E[i_cal]->Fill(delta_phi(max_tow_phi,phi),max_tow_e);
215  hDeltaPhi_iPhi[i_cal]->Fill(delta_phi(max_tow_phi,phi),max_tow_iphi);
216  hDeltaPhi_eta[i_cal]->Fill(eta,delta_phi(max_tow_phi,phi));
217  hCaloPhi[i_cal]->Fill(phi_in_range(max_tow_phi));
218  hTowerEt[i_cal]->Fill(max_tow_e);
219  hTowerEt_pionEt[i_cal]->Fill(max_tow_e,e);
220  // if (Cut(ientry) < 0) continue;
221  } // end calo loop
222 
223  for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
224  hDeltaPhi_fraction[i_cal]->Fill(delta_phi(max_out_phi,phi),total_E[i_cal]/(total_E[0]+total_E[1]+total_E[2]));
225  }
226 
227  hDeltaPhi_fraction_HcalOverAll->Fill(pt, (total_E[0]+total_E[1])/(total_E[0]+total_E[1]+total_E[2]), delta_phi(closestOuter_tow_phi,phi));
228  hDeltaPhi_fraction_oHcalOverHcals->Fill(pt, total_E[1]/(total_E[0]+total_E[1]), delta_phi(closestOuter_tow_phi,phi));
229 
230  hE_inner_vs_outer->Fill(pt, total_E[0], total_E[1]);
231  } // end event loop
232 
233 // new TCanvas;
234 // hGeantPhi->Draw();
235 // for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
236 // hPhi2D[i_cal]->Draw("COLZ");
237 // new TCanvas;
238 // hDeltaPhi[i_cal]->Draw("COLZ");
239 // new TCanvas;
240 // hDeltaPhi_pt[i_cal]->Draw("COLZ");
241 // new TCanvas;
242 // hCaloPhi[i_cal]->Draw();
243 // new TCanvas;
244 // hTowerEt[i_cal]->Draw();
245 // new TCanvas;
246 // hTowerEt_pionEt[i_cal]->Draw("COLZ");
247 // new TCanvas;
248 // hDeltaPhi_eta[i_cal]->Draw("COLZ");
249 // }
250 
251 
252  string outputroot = (string) inFileName;
253  string remove_this = ".root";
254  size_t pos = outputroot.find(remove_this);
255  if (pos != string::npos) { outputroot.erase(pos, remove_this.length()); }
256  TString outFileName = outputroot + "_histos.root";
257 
258  TFile *outFile = new TFile(outFileName,"RECREATE");
259  hEMCal3x3->Write();
260  hEMCal3x3_FINE->Write();
261  hGeantPhi->Write();
262  hE_inner_vs_outer->Write();
265  for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
266  hE_weighted_eta_phi[i_cal]->Write();
267  h_eta_phi[i_cal]->Write();
268  hPhi2D[i_cal]->Write();
269  hDeltaPhi_pt[i_cal]->Write();
270  hDeltaPhi_eta[i_cal]->Write();
271  hTowerEt_pionEt[i_cal]->Write();
272  hCaloEnergy_pionPt[i_cal]->Write();
273  hEnergy_fraction[i_cal]->Write();
274  hDeltaPhi_fraction[i_cal]->Write();
275  hDeltaPhi_E[i_cal]->Write();
276  hDeltaPhi_iPhi[i_cal]->Write();
277  hDeltaPhi[i_cal]->Write();
278  hCaloPhi[i_cal]->Write();
279  hTowerEt[i_cal]->Write();
280  }
281 
282  for (int i_cal=0; i_cal<ncal; ++i_cal) { // LOOP OVER CALORIMETER LAYERS
283  new TCanvas;
284  hDeltaPhi_fraction[i_cal]->Draw("COLZ");
285  }
286 
288 
289  TCanvas * can = new TCanvas( "can" , "" ,700 ,500 );
290  can->SetLogz();
291  const char us[] = "_";
292 
293  auto fa1 = new TF1("fa1","0.",0,1);
294 
295  for (int i=0; i<histos.size(); ++i ) {
296  TH2D *hTemp = (TH2D*)histos[i]->Project3D("ZY");
297  hTemp->SetLineColor(kRed);
298  hTemp->SetMarkerColor(kRed);
299  hTemp->Draw("COLZ");
300  hTemp->ProfileX()->Draw("SAME");
301  can->SaveAs(Form("plots/ExploreTTrees/%s.pdf",histos[i]->GetName()),"PDF");
302  hTemp->ProfileX()->Draw("");
303  fa1->Draw("SAME");
304  can->SaveAs(Form("plots/ExploreTTrees/%s_pfx.pdf",histos[i]->GetName()),"PDF");
305  for (int j=0; j<3; ++j) {
306  double ptlo = (double)10.*j;
307  double pthi = (double)10.*(j+1);
308  histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
309  hTemp = (TH2D*)histos[i]->Project3D("ZY");
310  hTemp->SetLineColor(kRed);
311  hTemp->SetMarkerColor(kRed);
312  hTemp->Draw("COLZ");
313  hTemp->ProfileX()->Draw("SAME");
314  can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
315  hTemp->ProfileX()->Draw("");
316  fa1->Draw("SAME");
317  can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i_pfx.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
318  }
319  for (int j=1; j<10; ++j) {
320  double ptlo = (double)1.*j;
321  double pthi = (double)1.*(j+1);
322  histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
323  hTemp = (TH2D*)histos[i]->Project3D("ZY");
324  hTemp->SetLineColor(kRed);
325  hTemp->SetMarkerColor(kRed);
326  hTemp->Draw("COLZ");
327  hTemp->ProfileX()->Draw("SAME");
328  can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
329  hTemp->ProfileX()->Draw("");
330  fa1->Draw("SAME");
331  can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i_pfx.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
332  }
333  histos[i]->GetXaxis()->SetRangeUser(0.,30.);
334  }
335 
336 // gStyle->SetPalette(kCMYK);
337 
338  hEMCal3x3->Draw("COLZ");
339  can->SaveAs("plots/ExploreTTrees/hEMCal3x3.pdf","PDF");
340 
341  hEMCal3x3_FINE->Draw("COLZ");
342  can->SaveAs("plots/ExploreTTrees/hEMCal3x3_FINE.pdf","PDF");
343 
344 
345  TF1 *f1 = new TF1("f1","gaus",0.,0.6);
346 
347  float MIP_value[29];//, MIP_sigma[29];
348 
349  TH1D *hEMCal3x3_FINE_py[30];
350  new TCanvas;
351  for (int i=1; i<30; ++i) {
352  int ptbin = i+1;
353  int ptlo = i;
354  int pthi = i+1;
355  hEMCal3x3_FINE_py[i-1] = (TH1D*)hEMCal3x3_FINE->ProjectionY(Form("hEMCal3x3_FINE_py__%i_%iGeV",ptlo,pthi),ptbin,ptbin);
356 // hEMCal3x3_FINE_py[i-1]->Scale(1./hEMCal3x3_FINE_py[i-1]->Integral());
357  hEMCal3x3_FINE_py[i-1]->SetTitle(Form("%i-%i GeV/c pion",ptlo,pthi));
358  hEMCal3x3_FINE_py[i-1]->SetStats(0);
359  hEMCal3x3_FINE_py[i-1]->GetYaxis()->SetRangeUser(0.0,80.0);
360  hEMCal3x3_FINE_py[i-1]->Draw("SAMEPLCPMC");
361  hEMCal3x3_FINE_py[i-1]->Fit(f1,"","",0.,0.6);
362  MIP_value[i-1] = f1->GetMaximumX();
363  }
364 
365  float mean_value = 0;
366  for (int i=0; i<29; ++i) {
367  cout<<i+1<<"-"<<i+2<<"GeV/c \t"<<MIP_value[i]<<endl;
368  mean_value += MIP_value[i];
369  }
370  mean_value/=29.;
371  cout<<endl<<mean_value<<endl;
372 
373  TLegend *leg0;
374 
375  leg0 = new TLegend(0.7, 0.3, 0.98, 0.98,NULL,"brNDC"); // LEGEND 0
376  leg0->SetBorderSize(0); leg0->SetLineColor(1); leg0->SetLineStyle(1);
377  leg0->SetLineWidth(1);
378 // leg0->SetFillColorAlpha(0,0.0);
379 // leg0->SetFillStyle(1001);
380  leg0->SetNColumns(2);
381  leg0->SetHeader("p_{T}^{pion}\t\t\tGaussian Mean","");
382 // leg0->AddEntry((TObject*)0,"EA_{Low} \t", "");
383 // leg0->AddEntry((TObject*)0,"\tEA_{High}", "");
384 
385  for (int i=0; i<29; ++i) {
386  string title = "";
387  leg0->AddEntry(hEMCal3x3_FINE_py[i],Form("%i-%i GeV/c",i+1,i+2),"lpf");
388  leg0->AddEntry("",Form("%f",MIP_value[i]),"");
389  }
390  leg0->Draw("SAME");
391 
392 }
393