Back to home page

sPhenix code displayed by LXR

 
 

    


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 // Fixed size dimensions of array or collections stored in the TTree if any.
0009 
0010 // Declaration of leaf types
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;   //!pointer to the analyzed TTree or TChain
0041 Int_t          fCurrent; //!current Tree number in a TChain
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 //TString inFileName = "HCalJetPhiShift_negativePion_magnetOn_15k.root";
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) {  // LOOP OVER CALORIMETER LAYERS
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) {  // LOOP OVER CALORIMETER LAYERS
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++) {  // LOOP OVER EVENTS
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     //    float max_tow_phi[ncal];// = {0.,0.,0.};
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) {  // LOOP OVER CALORIMETER LAYERS
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) {  // LOOP OVER CALORIMETER TOWERS
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 //        if (calo_e[i_cal]->at(itow)<=0.) continue;
0178 //        if (i_cal==2 && calo_e[i_cal]->at(itow)<0.08) { continue; }
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 ){  // MATCH PION TO CLOSEST OHCAL TOWER
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       }  // end tower loop
0199       
0200       if (i_cal==2) {
0201         hEMCal3x3->Fill(pt,energyIn3x3);
0202         hEMCal3x3_FINE->Fill(pt,energyIn3x3);
0203       }
0204 
0205 //      cout<<cal_tag[i_cal]<<"\t"<<max_tow_eta<<"\t"<<eta<<"\t"<<delta_phi(max_tow_phi,phi)<<endl;
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       // if (Cut(ientry) < 0) continue;
0221     }  // end calo loop
0222     
0223     for (int i_cal=0; i_cal<ncal; ++i_cal) {  // LOOP OVER CALORIMETER LAYERS
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   }  // end event loop
0232   
0233 //  new TCanvas;
0234 //  hGeantPhi->Draw();
0235 //  for (int i_cal=0; i_cal<ncal; ++i_cal) {  // LOOP OVER CALORIMETER LAYERS
0236 //    hPhi2D[i_cal]->Draw("COLZ");
0237 //    new TCanvas;
0238 //    hDeltaPhi[i_cal]->Draw("COLZ");
0239 //    new TCanvas;
0240 //    hDeltaPhi_pt[i_cal]->Draw("COLZ");
0241 //    new TCanvas;
0242 //    hCaloPhi[i_cal]->Draw();
0243 //    new TCanvas;
0244 //    hTowerEt[i_cal]->Draw();
0245 //    new TCanvas;
0246 //    hTowerEt_pionEt[i_cal]->Draw("COLZ");
0247 //    new TCanvas;
0248 //    hDeltaPhi_eta[i_cal]->Draw("COLZ");
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) {  // LOOP OVER CALORIMETER LAYERS
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) {  // LOOP OVER CALORIMETER LAYERS
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 //  gStyle->SetPalette(kCMYK);
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];//, MIP_sigma[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 //    hEMCal3x3_FINE_py[i-1]->Scale(1./hEMCal3x3_FINE_py[i-1]->Integral());
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");    // LEGEND 0
0376   leg0->SetBorderSize(0);  leg0->SetLineColor(1);  leg0->SetLineStyle(1);
0377   leg0->SetLineWidth(1);
0378 //  leg0->SetFillColorAlpha(0,0.0);
0379 //  leg0->SetFillStyle(1001);
0380   leg0->SetNColumns(2);
0381   leg0->SetHeader("p_{T}^{pion}\t\t\tGaussian Mean","");
0382 //    leg0->AddEntry((TObject*)0,"EA_{Low} \t", "");
0383 //    leg0->AddEntry((TObject*)0,"\tEA_{High}", "");
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