Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EMCalTowerRingStudy.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EMCalTowerRingStudy.C
1 #include <iostream>
2 #include <string>
3 #include <vector>
4 #include <TFile.h>
5 #include <TTree.h>
6 #include <TAxis.h>
7 #include <TChain.h>
8 
9 // Fixed size dimensions of array or collections stored in the TTree if any.
10 
11 // Declaration of leaf types
12 Int_t event;
13 Float_t eta;
14 Float_t phi;
15 Float_t e;
16 Float_t pt;
17 Float_t vx;
18 Float_t vy;
19 Float_t vz;
20 vector<int> id;
21 
22 int nTow_in;
23 vector<float> eta_in;
24 vector<float> phi_in;
25 vector<float> e_in;
26 vector<int> ieta_in;
27 vector<int> iphi_in;
28 int nTow_out;
29 vector<float> eta_out;
30 vector<float> phi_out;
31 vector<float> e_out;
32 vector<int> ieta_out;
33 vector<int> iphi_out;
34 int nTow_emc;
35 vector<float> eta_emc;
36 vector<float> phi_emc;
37 vector<float> e_emc;
38 vector<int> ieta_emc;
39 vector<int> iphi_emc;
40 
41 TTree *fChain;
42 Int_t fCurrent;
43 
44 
45 const int ncal = 3;
46 string cal_name[] = { "iHCal", "oHCal", "EMCal" };
47 string cal_tag[] = { "_in", "_out", "_emc" };
48 TString inFileName = "HCalJetPhiShift_10k.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 
56 int *nTowers[ncal] = {&nTow_in, &nTow_out, &nTow_emc};
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] = {&phi_in, &phi_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 phi_in_range(double phi){
71  while (phi>=M_PI) { phi -= 2.*M_PI; }
72  while (phi<=-M_PI) { phi += 2.*M_PI; }
73  return phi;
74 };
75 
76 bool tower_in_square(float calo_eta, float calo_phi, float eta, float phi, int N){
77  double min_eta = eta - 0.024*N - 0.012;
78  double max_eta = eta + 0.024*N + 0.012;
79  if ( calo_eta<min_eta || calo_eta>max_eta ) { return false; }
80  double min_phi = phi - (M_PI/128.)*N - (M_PI/256.);
81  double max_phi = phi + (M_PI/128.)*N + (M_PI/256.);
82  if ( calo_phi<min_phi || calo_phi>max_phi ) { return false; }
83  if ( calo_eta>=min_eta && calo_eta<=max_eta && calo_phi>=min_phi && calo_phi<=max_phi ) { return true; }
84  return false;
85 };
86 
87 bool tower_in_ring(float calo_eta, float calo_phi, float eta, float phi, int N){
88  double n_rings = (double) N;
89  double eta_lo = eta - 0.024*n_rings;
90  double eta_hi = eta + 0.024*n_rings;
91  double phi_lo = phi_in_range(phi - (M_PI/128.)*n_rings);// - (M_PI/256.);
92  double phi_hi = phi_in_range(phi + (M_PI/128.)*n_rings);// + (M_PI/256.);
93 // if (fabs(delta_phi(calo_phi,phi_lo))>(M_PI/256.) || fabs(delta_phi(calo_phi,phi_hi))>(M_PI/256.)) { return false; }
94 // if ( fabs(eta_lo-calo_eta)>0.012 || fabs(eta_hi-calo_eta)>0.012 ) { return false; }
95  if (fabs(delta_phi(phi_in_range(calo_phi),phi_in_range(phi_lo)))<(M_PI/256.) || fabs(delta_phi(phi_in_range(calo_phi),phi_in_range(phi_hi)))<(M_PI/256.)) {
96  if ( calo_eta>eta_lo-0.012 && calo_eta<eta_hi+0.012 ) { return true; }
97  }
98  if ( fabs(eta_lo-calo_eta)<=0.012 || fabs(eta_hi-calo_eta)<=0.012 ) {
99  if ( calo_phi>phi_lo-(M_PI/256.) && calo_phi<phi_hi+(M_PI/256.) ) { return true; }
100  }
101  return false;
102 };
103 
104 const int i_cal = 2; // THIS IS AN EMCAL STUDY
105 //const int ref_cal = -1; // STUDY WRT PION
106 //const int ref_cal = 0; // STUDY WRT IHCAL
107 //const int ref_cal = 1; // STUDY WRT OHCAL
108 const int ref_cal = 2; // STUDY WRT EMCAL
109 
110 void EMCalTowerRingStudy() {
111 
112  TH3D *hEsum_per_nRings = new TH3D("hEsum_per_nRings",";pion p_{T} [GeV/c];E sum in tower square [GeV];n_{rings}",30,0.,30.,100,0.,100.,30,0.,30.);
113  TH3D *hEsum_in_ringN = new TH3D("hEsum_in_ringN",";pion p_{T} [GeV/c];E sum in tower ring [GeV];n_{rings}",30,0.,30.,200,0.,40.,30,0.,30.);
114  TH3D *hEMCalE_over_pionE = new TH3D("hEMCalE_over_pionE",";pion p_{T} [GeV/c];E_{EMCal}/E_{pion};n_{rings}",30,0.,30.,115,0.,1.15,30,0.,30.);
115 
116  TH3D *hTotalCaloSum_per_nRings = new TH3D("hTotalCaloSum_per_nRings",";pion p_{T} [GeV/c];E sum in all calos [GeV];",30,0.,30.,200,0.,100.,30,0.,30.);
117  TH3D *hTotalCaloEfrac_per_nRings = new TH3D("hTotalCaloEfrac_per_nRings",";pion p_{T} [GeV/c];total calo E / pion E;",30,0.,30.,70,0.,7.,30,0.,30.);
118 // hTotalCaloSum_per_nRings->GetXaxis()->SetNdivisions(29);
119  for (int i=0; i<30; ++i){ hTotalCaloSum_per_nRings->GetXaxis()->SetBinLabel(i+1,Form("%ix%i",2*i+1,2*i+1)); }
120  hTotalCaloSum_per_nRings->GetXaxis()->SetAlphanumeric();
121 
122  TFile *f = new TFile(inFileName,"READ");
123  TTree *fChain = (TTree*)f->Get("T");
124 
125  fChain->SetBranchAddress("event", &event);
126  fChain->SetBranchAddress("eta", &eta);
127  fChain->SetBranchAddress("phi", &phi);
128  fChain->SetBranchAddress("e", &e);
129  fChain->SetBranchAddress("pt", &pt);
130  fChain->SetBranchAddress("vx", &vx);
131  fChain->SetBranchAddress("vy", &vy);
132  fChain->SetBranchAddress("vz", &vz);
133 
134  for (int ical=0; ical<ncal; ++ical) { // LOOP OVER CALORIMETER LAYERS
135  fChain->SetBranchAddress(Form("nTow%s",cal_tag[ical].c_str()), &nTowers[ical]);
136  fChain->SetBranchAddress(Form("eta%s",cal_tag[ical].c_str()), &calo_eta[ical]);
137  fChain->SetBranchAddress(Form("phi%s",cal_tag[ical].c_str()), &calo_phi[ical]);
138  fChain->SetBranchAddress(Form("e%s",cal_tag[ical].c_str()), &calo_e[ical]);
139  fChain->SetBranchAddress(Form("ieta%s",cal_tag[ical].c_str()), &calo_ieta[ical]);
140  fChain->SetBranchAddress(Form("iphi%s",cal_tag[ical].c_str()), &calo_iphi[ical]);
141  }
142 
143 
144  Long64_t nentries = fChain->GetEntriesFast();
145 
146  for (Long64_t jentry=0; jentry<nentries;jentry++) { // LOOP OVER EVENTS
147  Long64_t ientry = fChain->GetTreeNumber();
148  if (ientry < 0) break;
149 
150  fChain->GetEntry(jentry);
151 
152  double total_calo_E[ncal] = {0.,0.,0.};
153 
154  float ref_e = 0.;
155  float ref_eta;
156  float ref_phi;
157  int ref_ieta;
158  int ref_iphi;
159 
160  if (ref_cal==-1) {
161  ref_e = e;
162  ref_eta = eta;
163  ref_phi = phi;
164 // ref_phi = phi_in_range(phi+M_PI);
165  }
166  else if (ref_cal>=0 && ref_cal<=2){
167  for (int itow=0; itow<calo_e[ref_cal]->size(); ++itow) { // LOOP OVER REFERENCE CAL TOWERS
168  if (calo_e[ref_cal]->at(itow)>ref_e) {
169  ref_e = calo_e[ref_cal]->at(itow);
170  ref_eta = calo_eta[ref_cal]->at(itow);
171  ref_phi = calo_phi[ref_cal]->at(itow);
172 // ref_phi = phi_in_range( calo_phi[ref_cal]->at(itow) + M_PI );
173  }
174  }
175  }
176  else { cerr<<"ope"<<endl; }
177 
178  for (int index=0; index<2; ++index) { // LOOP OVER HCALS
179  for (int itow=0; itow<calo_e[index]->size(); ++itow) { // LOOP OVER CAL TOWERS
180  total_calo_E[index] += calo_e[index]->at(itow);
181  }
182  }
183 
184  int n_towers_in_square = 0;
185  double square_e_sum = 0.;
186  for (int n_rings = 0; n_rings<30; ++n_rings) {
187 
188  double ring_e_sum = 0.;
189  int n_towers_in_ring = 0;
190  for (int itow=0; itow<calo_e[i_cal]->size(); ++itow) { // LOOP OVER EMCAL TOWERS
191 // if (i_cal==2 && calo_e[i_cal]->at(itow)<0.1) { continue; }
192  if ( tower_in_ring(calo_eta[i_cal]->at(itow),calo_phi[i_cal]->at(itow),ref_eta,ref_phi,n_rings) ) {
193  total_calo_E[2] += calo_e[i_cal]->at(itow);
194  ring_e_sum += calo_e[i_cal]->at(itow);
195  square_e_sum += calo_e[i_cal]->at(itow);
196  n_towers_in_square += 1;
197  n_towers_in_ring += 1;
198  }
199  } // end tower loop
200  hEsum_per_nRings->Fill(pt,square_e_sum,n_rings);
201  hEsum_in_ringN->Fill(pt,ring_e_sum,n_rings);
202  hEMCalE_over_pionE->Fill(pt,ring_e_sum/e,n_rings);
203 
204  hTotalCaloSum_per_nRings->Fill(e,square_e_sum+total_calo_E[0]+total_calo_E[1],n_rings);
205  hTotalCaloEfrac_per_nRings->Fill(e,(square_e_sum+total_calo_E[0]+total_calo_E[1])/e,n_rings);
206  }
207 
208  // if (Cut(ientry) < 0) continue;
209 
210  } // end event loop
211 
212 
213 
214  string outputroot = (string) inFileName;
215  string remove_this = ".root";
216  size_t pos = outputroot.find(remove_this);
217  if (pos != string::npos) { outputroot.erase(pos, remove_this.length()); }
218  if (ref_cal>=0 && ref_cal<=2){
219  outputroot += cal_tag[ref_cal] + "Ref";
220  }
221  TString outFileName = outputroot + "_EMCalRings.root";
222 
223  TFile *outFile = new TFile(outFileName,"RECREATE");
224  hEsum_per_nRings->Write();
225  hEsum_in_ringN->Write();
226  hTotalCaloSum_per_nRings->Write();
227  hTotalCaloEfrac_per_nRings->Write();
228  hEMCalE_over_pionE->Write();
229  // for (int ical=0; ical<ncal; ++ical) { // LOOP OVER CALORIMETER LAYERS
230  // }
231 
232  vector<TH3D *> histos = {hEsum_per_nRings, hEsum_in_ringN, hTotalCaloSum_per_nRings, hTotalCaloEfrac_per_nRings, hEMCalE_over_pionE};
233  auto fa1 = new TF1("fa1","0.",0,1);
234 
235  TCanvas * can = new TCanvas( "can" , "" ,700 ,500 );
236  can->SetLogz();
237  const char us[] = "_";
238 
239  for (int i=0; i<histos.size(); ++i ) {
240  histos[i]->Project3D("YZ")->Draw("COLZ");
241  can->SaveAs(Form("plots/EMCalTowerRingStudy/%s.pdf",histos[i]->GetName()),"PDF");
242  for (int j=0; j<3; ++j) {
243  double ptlo = (double)10.*j;
244  double pthi = (double)10.*(j+1);
245  histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
246  histos[i]->Project3D("YZ")->Draw("COLZ");
247  can->SaveAs(Form("plots/EMCalTowerRingStudy/%s%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
248  }
249  if (i==4){
250  histos[i]->GetXaxis()->SetRangeUser(1.,2.);
251  TH2D *hDenom2D = (TH2D*)histos[i]->Project3D("YZ");
252  TH1D *hDenom = (TH1D*)hDenom2D->ProfileX();
253  hDenom->SetName("hDenom");
254 
255  for (int j=1; j<10; ++j) {
256  double ptlo = (double)1.*j;
257  double pthi = (double)1.*(j+1);
258  histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
259  histos[i]->Project3D("YZ")->Draw("COLZ");
260  TH2D *hTemp = (TH2D*)histos[i]->Project3D("YZ");
261  hTemp->SetLineColor(kRed);
262  hTemp->SetMarkerColor(kRed);
263  TH1D *hNum = (TH1D*)hTemp->ProfileX();
264  hNum->Draw("SAME");
265  can->SaveAs(Form("plots/EMCalTowerRingStudy/%s_FINE_%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
266  hNum->Divide(hDenom);
267  hNum->GetYaxis()->SetRangeUser(0.1,2.2);
268  hNum->Draw();
269  fa1->Draw("SAME");
270  can->SaveAs(Form("plots/EMCalTowerRingStudy/%s_FINE_%s%i%s%i__ratioTo1_2.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
271  hNum->Reset();
272  }
273  }
274  histos[i]->GetXaxis()->SetRangeUser(0.,30.);
275  }
276 
277 // hEsum_per_nRings->Project3D("YZ")->Draw("COLZ");
278 // can->SaveAs("plots/EMCalTowerRingStudy/hEsum_per_nRings.pdf","PDF");
279 // hEsum_in_ringN->Project3D("YZ")->Draw("COLZ");
280 // can->SaveAs("plots/EMCalTowerRingStudy/hEsum_in_ringN.pdf","PDF");
281 // hTotalCaloSum_per_nRings->Project3D("YZ")->Draw("COLZ");
282 // can->SaveAs("plots/EMCalTowerRingStudy/hTotalCaloSum_per_nRings.pdf","PDF");
283 // hTotalCaloEfrac_per_nRings->Project3D("YZ")->Draw("COLZ");
284 // can->SaveAs("plots/EMCalTowerRingStudy/hTotalCaloEfrac_per_nRings.pdf","PDF");
285 }