File indexing completed on 2025-08-05 08:13:16
0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004 #include <TFile.h>
0005 #include <TTree.h>
0006 #include <TChain.h>
0007
0008
0009
0010
0011 Int_t event;
0012 Float_t eta;
0013 Float_t phi;
0014 Float_t e;
0015 Float_t pt;
0016 Float_t vx;
0017 Float_t vy;
0018 Float_t vz;
0019 vector<int> id;
0020
0021 int nTow_in;
0022 vector<float> eta_in;
0023 vector<float> phi_in;
0024 vector<float> e_in;
0025 vector<int> ieta_in;
0026 vector<int> iphi_in;
0027 int nTow_out;
0028 vector<float> eta_out;
0029 vector<float> phi_out;
0030 vector<float> e_out;
0031 vector<int> ieta_out;
0032 vector<int> iphi_out;
0033 int nTow_emc;
0034 vector<float> eta_emc;
0035 vector<float> phi_emc;
0036 vector<float> e_emc;
0037 vector<int> ieta_emc;
0038 vector<int> iphi_emc;
0039
0040 TTree *fChain;
0041 Int_t fCurrent;
0042
0043
0044 const int ncal = 3;
0045 string cal_name[] = { "iHCal", "oHCal", "EMCal" };
0046 string cal_tag[] = { "_in", "_out", "_emc" };
0047 TString inFileName = "HCalJetPhiShift_10k.root";
0048
0049 int n_pt_bins[ncal] = {80,100,100};
0050 double pt_bin_lo[ncal] = {0.,0.,0.};
0051 double pt_bin_hi[ncal] = {4.,50.,20.};
0052 double eta_max[ncal] = {1.1,1.1,1.152};
0053 int n_eta_bins[ncal] = {24,24,96};
0054 int n_phi_bins[ncal] = {64,64,256};
0055
0056 int *nTowers[ncal] = {&nTow_in, &nTow_out, &nTow_emc};
0057 vector<float> *calo_eta[ncal] = {&eta_in, &eta_out, &eta_emc};
0058 vector<float> *calo_phi[ncal] = {&phi_in, &phi_out, &phi_emc};
0059 vector<float> *calo_e[ncal] = {&e_in, &e_out, &e_emc};
0060 vector<int> *calo_ieta[ncal] = {&ieta_in, &ieta_out, &ieta_emc};
0061 vector<int> *calo_iphi[ncal] = {&iphi_in, &iphi_out, &iphi_emc};
0062
0063 double delta_phi(double phi1, double phi2){
0064 double dPhi = phi1 - phi2;
0065 while (dPhi>=M_PI) { dPhi -= 2.*M_PI; }
0066 while (dPhi<=-M_PI) { dPhi += 2.*M_PI; }
0067 return dPhi;
0068 };
0069
0070 double dR(double phi1, double phi2, double eta1, double eta2){
0071 double dphi_sqrd = delta_phi( phi1, phi2)*delta_phi( phi1, phi2);
0072 double deta_sqrd = (eta1-eta2)*(eta1-eta2);
0073 return sqrt(dphi_sqrd+deta_sqrd);
0074 };
0075
0076 double phi_in_range(double phi){
0077 while (phi>=M_PI) { phi -= 2.*M_PI; }
0078 while (phi<=-M_PI) { phi += 2.*M_PI; }
0079 return phi;
0080 };
0081
0082 bool tower_in_3x3(float calo_eta, float calo_phi, float eta, float phi){
0083 double min_eta = eta - 0.024 - 0.012;
0084 double max_eta = eta + 0.024 + 0.012;
0085 if ( calo_eta<min_eta || calo_eta>max_eta ) { return false; }
0086 double min_phi = phi - (M_PI/128.) - (M_PI/256.);
0087 double max_phi = phi + (M_PI/128.) + (M_PI/256.);
0088 if ( calo_phi<min_phi || calo_phi>max_phi ) { return false; }
0089 if ( calo_eta>=min_eta && calo_eta<=max_eta && calo_phi>=min_phi && calo_phi<=max_phi ) { return true; }
0090 return false;
0091 };
0092
0093 TH1D *hGeantPhi = new TH1D("hGeantPhi",";truth #phi",64,-M_PI,M_PI);
0094 TH2D *hEMCal3x3 = new TH2D("hEMCal3x3","pion p_{T};E_{EMCal}^{3x3} about max-E tower [GeV]",30,0.,30.,120,0.,60.);
0095 TH2D *hEMCal3x3_FINE = new TH2D("hEMCal3x3_FINE","pion p_{T};E_{EMCal}^{3x3} about max-E tower [GeV]",30,0.,30.,200.,0.,2.);
0096 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.);
0097 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);
0098 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);
0099
0100 TH3D *hE_weighted_eta_phi[ncal], *h_eta_phi[ncal];
0101 TH2D *hPhi2D[ncal], *hDeltaPhi_pt[ncal], *hDeltaPhi_eta[ncal], *hTowerEt_pionEt[ncal], *hCaloEnergy_pionPt[ncal], *hEnergy_fraction[ncal], *hDeltaPhi_E[ncal], *hDeltaPhi_iPhi[ncal], *hDeltaPhi_fraction[ncal];
0102
0103 TH1D *hDeltaPhi[ncal], *hCaloPhi[ncal], *hTowerEt[ncal];
0104
0105 void ExploreTTrees() {
0106
0107 for (int i_cal=0; i_cal<ncal; ++i_cal) {
0108 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);
0109 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);
0110 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);
0111 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.);
0112 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]);
0113 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.);
0114 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.);
0115 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.);
0116 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.);
0117 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]);
0118 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.);
0119 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);
0120 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);
0121 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]);
0122 }
0123
0124 TFile *f = new TFile(inFileName,"READ");
0125 TTree *fChain = (TTree*)f->Get("T");
0126
0127 fChain->SetBranchAddress("event", &event);
0128 fChain->SetBranchAddress("eta", &eta);
0129 fChain->SetBranchAddress("phi", &phi);
0130 fChain->SetBranchAddress("e", &e);
0131 fChain->SetBranchAddress("pt", &pt);
0132 fChain->SetBranchAddress("vx", &vx);
0133 fChain->SetBranchAddress("vy", &vy);
0134 fChain->SetBranchAddress("vz", &vz);
0135
0136 for (int i_cal=0; i_cal<ncal; ++i_cal) {
0137 fChain->SetBranchAddress(Form("nTow%s",cal_tag[i_cal].c_str()), &nTowers[i_cal]);
0138 fChain->SetBranchAddress(Form("eta%s",cal_tag[i_cal].c_str()), &calo_eta[i_cal]);
0139 fChain->SetBranchAddress(Form("phi%s",cal_tag[i_cal].c_str()), &calo_phi[i_cal]);
0140 fChain->SetBranchAddress(Form("e%s",cal_tag[i_cal].c_str()), &calo_e[i_cal]);
0141 fChain->SetBranchAddress(Form("ieta%s",cal_tag[i_cal].c_str()), &calo_ieta[i_cal]);
0142 fChain->SetBranchAddress(Form("iphi%s",cal_tag[i_cal].c_str()), &calo_iphi[i_cal]);
0143 }
0144
0145 Long64_t nentries = fChain->GetEntriesFast();
0146
0147 for (Long64_t jentry=0; jentry<nentries;jentry++) {
0148 Long64_t ientry = fChain->GetTreeNumber();
0149 if (ientry < 0) break;
0150
0151 fChain->GetEntry(jentry);
0152
0153 hGeantPhi->Fill(phi_in_range(phi));
0154
0155 float total_E[ncal] = {0.,0.,0.};
0156 float max_out_phi;
0157
0158
0159 float closestOuter_tow_dR = 9999.;
0160 float closestOuter_tow_eta;
0161 float closestOuter_tow_phi;
0162
0163 float energyIn3x3 = 0;
0164
0165 for (int i_cal=0; i_cal<ncal; ++i_cal) {
0166 float max_tow_e = 0.;
0167 float max_tow_eta;
0168 int max_tow_iphi;
0169 float max_tow_phi;
0170
0171 for (int itow=0; itow<calo_e[i_cal]->size(); ++itow) {
0172
0173 if (i_cal==2 && tower_in_3x3(calo_eta[i_cal]->at(itow), calo_phi[i_cal]->at(itow), eta, phi) ){
0174 energyIn3x3 += calo_e[i_cal]->at(itow);
0175 }
0176
0177
0178
0179
0180 total_E[i_cal] += calo_e[i_cal]->at(itow);
0181 if (calo_e[i_cal]->at(itow)>max_tow_e){
0182 max_tow_e = calo_e[i_cal]->at(itow);
0183 max_tow_eta = calo_eta[i_cal]->at(itow);
0184 max_tow_phi = calo_phi[i_cal]->at(itow);
0185 max_tow_iphi = calo_iphi[i_cal]->at(itow);
0186 }
0187 if( i_cal==1 ){
0188 float delta_R = dR(phi, calo_phi[i_cal]->at(itow), eta, calo_eta[i_cal]->at(itow));
0189 if (delta_R<closestOuter_tow_dR){
0190 closestOuter_tow_dR = delta_R;
0191 closestOuter_tow_eta = calo_eta[i_cal]->at(itow);
0192 closestOuter_tow_phi = calo_phi[i_cal]->at(itow);
0193 }
0194 }
0195 hCaloEnergy_pionPt[i_cal]->Fill(pt,calo_e[i_cal]->at(itow));
0196 h_eta_phi[i_cal]->Fill(pt,calo_eta[i_cal]->at(itow),phi_in_range(calo_phi[i_cal]->at(itow)));
0197 hE_weighted_eta_phi[i_cal]->Fill(pt,calo_eta[i_cal]->at(itow),phi_in_range(calo_phi[i_cal]->at(itow)),calo_e[i_cal]->at(itow));
0198 }
0199
0200 if (i_cal==2) {
0201 hEMCal3x3->Fill(pt,energyIn3x3);
0202 hEMCal3x3_FINE->Fill(pt,energyIn3x3);
0203 }
0204
0205
0206
0207 if (i_cal==1) {max_out_phi=max_tow_phi;}
0208
0209 double e_fraction = (double) total_E[i_cal]/e;
0210 hEnergy_fraction[i_cal]->Fill(e,e_fraction);
0211 hPhi2D[i_cal]->Fill(phi_in_range(phi),max_tow_phi);
0212 hDeltaPhi[i_cal]->Fill(delta_phi(max_tow_phi,phi));
0213 hDeltaPhi_pt[i_cal]->Fill(pt,delta_phi(max_tow_phi,phi));
0214 hDeltaPhi_E[i_cal]->Fill(delta_phi(max_tow_phi,phi),max_tow_e);
0215 hDeltaPhi_iPhi[i_cal]->Fill(delta_phi(max_tow_phi,phi),max_tow_iphi);
0216 hDeltaPhi_eta[i_cal]->Fill(eta,delta_phi(max_tow_phi,phi));
0217 hCaloPhi[i_cal]->Fill(phi_in_range(max_tow_phi));
0218 hTowerEt[i_cal]->Fill(max_tow_e);
0219 hTowerEt_pionEt[i_cal]->Fill(max_tow_e,e);
0220
0221 }
0222
0223 for (int i_cal=0; i_cal<ncal; ++i_cal) {
0224 hDeltaPhi_fraction[i_cal]->Fill(delta_phi(max_out_phi,phi),total_E[i_cal]/(total_E[0]+total_E[1]+total_E[2]));
0225 }
0226
0227 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));
0228 hDeltaPhi_fraction_oHcalOverHcals->Fill(pt, total_E[1]/(total_E[0]+total_E[1]), delta_phi(closestOuter_tow_phi,phi));
0229
0230 hE_inner_vs_outer->Fill(pt, total_E[0], total_E[1]);
0231 }
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252 string outputroot = (string) inFileName;
0253 string remove_this = ".root";
0254 size_t pos = outputroot.find(remove_this);
0255 if (pos != string::npos) { outputroot.erase(pos, remove_this.length()); }
0256 TString outFileName = outputroot + "_histos.root";
0257
0258 TFile *outFile = new TFile(outFileName,"RECREATE");
0259 hEMCal3x3->Write();
0260 hEMCal3x3_FINE->Write();
0261 hGeantPhi->Write();
0262 hE_inner_vs_outer->Write();
0263 hDeltaPhi_fraction_HcalOverAll->Write();
0264 hDeltaPhi_fraction_oHcalOverHcals->Write();
0265 for (int i_cal=0; i_cal<ncal; ++i_cal) {
0266 hE_weighted_eta_phi[i_cal]->Write();
0267 h_eta_phi[i_cal]->Write();
0268 hPhi2D[i_cal]->Write();
0269 hDeltaPhi_pt[i_cal]->Write();
0270 hDeltaPhi_eta[i_cal]->Write();
0271 hTowerEt_pionEt[i_cal]->Write();
0272 hCaloEnergy_pionPt[i_cal]->Write();
0273 hEnergy_fraction[i_cal]->Write();
0274 hDeltaPhi_fraction[i_cal]->Write();
0275 hDeltaPhi_E[i_cal]->Write();
0276 hDeltaPhi_iPhi[i_cal]->Write();
0277 hDeltaPhi[i_cal]->Write();
0278 hCaloPhi[i_cal]->Write();
0279 hTowerEt[i_cal]->Write();
0280 }
0281
0282 for (int i_cal=0; i_cal<ncal; ++i_cal) {
0283 new TCanvas;
0284 hDeltaPhi_fraction[i_cal]->Draw("COLZ");
0285 }
0286
0287 vector<TH3D *> histos = {hDeltaPhi_fraction_HcalOverAll, hDeltaPhi_fraction_oHcalOverHcals};
0288
0289 TCanvas * can = new TCanvas( "can" , "" ,700 ,500 );
0290 can->SetLogz();
0291 const char us[] = "_";
0292
0293 auto fa1 = new TF1("fa1","0.",0,1);
0294
0295 for (int i=0; i<histos.size(); ++i ) {
0296 TH2D *hTemp = (TH2D*)histos[i]->Project3D("ZY");
0297 hTemp->SetLineColor(kRed);
0298 hTemp->SetMarkerColor(kRed);
0299 hTemp->Draw("COLZ");
0300 hTemp->ProfileX()->Draw("SAME");
0301 can->SaveAs(Form("plots/ExploreTTrees/%s.pdf",histos[i]->GetName()),"PDF");
0302 hTemp->ProfileX()->Draw("");
0303 fa1->Draw("SAME");
0304 can->SaveAs(Form("plots/ExploreTTrees/%s_pfx.pdf",histos[i]->GetName()),"PDF");
0305 for (int j=0; j<3; ++j) {
0306 double ptlo = (double)10.*j;
0307 double pthi = (double)10.*(j+1);
0308 histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
0309 hTemp = (TH2D*)histos[i]->Project3D("ZY");
0310 hTemp->SetLineColor(kRed);
0311 hTemp->SetMarkerColor(kRed);
0312 hTemp->Draw("COLZ");
0313 hTemp->ProfileX()->Draw("SAME");
0314 can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0315 hTemp->ProfileX()->Draw("");
0316 fa1->Draw("SAME");
0317 can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i_pfx.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0318 }
0319 for (int j=1; j<10; ++j) {
0320 double ptlo = (double)1.*j;
0321 double pthi = (double)1.*(j+1);
0322 histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
0323 hTemp = (TH2D*)histos[i]->Project3D("ZY");
0324 hTemp->SetLineColor(kRed);
0325 hTemp->SetMarkerColor(kRed);
0326 hTemp->Draw("COLZ");
0327 hTemp->ProfileX()->Draw("SAME");
0328 can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0329 hTemp->ProfileX()->Draw("");
0330 fa1->Draw("SAME");
0331 can->SaveAs(Form("plots/ExploreTTrees/%s%s%i%s%i_pfx.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0332 }
0333 histos[i]->GetXaxis()->SetRangeUser(0.,30.);
0334 }
0335
0336
0337
0338 hEMCal3x3->Draw("COLZ");
0339 can->SaveAs("plots/ExploreTTrees/hEMCal3x3.pdf","PDF");
0340
0341 hEMCal3x3_FINE->Draw("COLZ");
0342 can->SaveAs("plots/ExploreTTrees/hEMCal3x3_FINE.pdf","PDF");
0343
0344
0345 TF1 *f1 = new TF1("f1","gaus",0.,0.6);
0346
0347 float MIP_value[29];
0348
0349 TH1D *hEMCal3x3_FINE_py[30];
0350 new TCanvas;
0351 for (int i=1; i<30; ++i) {
0352 int ptbin = i+1;
0353 int ptlo = i;
0354 int pthi = i+1;
0355 hEMCal3x3_FINE_py[i-1] = (TH1D*)hEMCal3x3_FINE->ProjectionY(Form("hEMCal3x3_FINE_py__%i_%iGeV",ptlo,pthi),ptbin,ptbin);
0356
0357 hEMCal3x3_FINE_py[i-1]->SetTitle(Form("%i-%i GeV/c pion",ptlo,pthi));
0358 hEMCal3x3_FINE_py[i-1]->SetStats(0);
0359 hEMCal3x3_FINE_py[i-1]->GetYaxis()->SetRangeUser(0.0,80.0);
0360 hEMCal3x3_FINE_py[i-1]->Draw("SAMEPLCPMC");
0361 hEMCal3x3_FINE_py[i-1]->Fit(f1,"","",0.,0.6);
0362 MIP_value[i-1] = f1->GetMaximumX();
0363 }
0364
0365 float mean_value = 0;
0366 for (int i=0; i<29; ++i) {
0367 cout<<i+1<<"-"<<i+2<<"GeV/c \t"<<MIP_value[i]<<endl;
0368 mean_value += MIP_value[i];
0369 }
0370 mean_value/=29.;
0371 cout<<endl<<mean_value<<endl;
0372
0373 TLegend *leg0;
0374
0375 leg0 = new TLegend(0.7, 0.3, 0.98, 0.98,NULL,"brNDC");
0376 leg0->SetBorderSize(0); leg0->SetLineColor(1); leg0->SetLineStyle(1);
0377 leg0->SetLineWidth(1);
0378
0379
0380 leg0->SetNColumns(2);
0381 leg0->SetHeader("p_{T}^{pion}\t\t\tGaussian Mean","");
0382
0383
0384
0385 for (int i=0; i<29; ++i) {
0386 string title = "";
0387 leg0->AddEntry(hEMCal3x3_FINE_py[i],Form("%i-%i GeV/c",i+1,i+2),"lpf");
0388 leg0->AddEntry("",Form("%f",MIP_value[i]),"");
0389 }
0390 leg0->Draw("SAME");
0391
0392 }
0393