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 <TAxis.h>
0007 #include <TChain.h>
0008 
0009 // Fixed size dimensions of array or collections stored in the TTree if any.
0010 
0011 // Declaration of leaf types
0012 Int_t           event;
0013 Float_t         eta;
0014 Float_t         phi;
0015 Float_t         e;
0016 Float_t         pt;
0017 Float_t         vx;
0018 Float_t         vy;
0019 Float_t         vz;
0020 vector<int>     id;
0021 
0022 int nTow_in;
0023 vector<float>   eta_in;
0024 vector<float>   phi_in;
0025 vector<float>   e_in;
0026 vector<int>   ieta_in;
0027 vector<int>   iphi_in;
0028 int nTow_out;
0029 vector<float>   eta_out;
0030 vector<float>   phi_out;
0031 vector<float>   e_out;
0032 vector<int>   ieta_out;
0033 vector<int>   iphi_out;
0034 int nTow_emc;
0035 vector<float>   eta_emc;
0036 vector<float>   phi_emc;
0037 vector<float>   e_emc;
0038 vector<int>   ieta_emc;
0039 vector<int>   iphi_emc;
0040 
0041 TTree          *fChain;   //!pointer to the analyzed TTree or TChain
0042 Int_t          fCurrent; //!current Tree number in a TChain
0043 
0044 
0045 const int ncal = 3;
0046 string cal_name[] = { "iHCal", "oHCal", "EMCal" };
0047 string cal_tag[] = { "_in", "_out", "_emc" };
0048 TString inFileName = "HCalJetPhiShift_10k.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] = {&phi_in, &phi_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 phi_in_range(double phi){
0071   while (phi>=M_PI) { phi -= 2.*M_PI; }
0072   while (phi<=-M_PI) { phi += 2.*M_PI; }
0073   return phi;
0074 };
0075 
0076 bool tower_in_square(float calo_eta, float calo_phi, float eta, float phi, int N){
0077   double min_eta = eta - 0.024*N - 0.012;
0078   double max_eta = eta + 0.024*N + 0.012;
0079   if ( calo_eta<min_eta || calo_eta>max_eta ) { return false; }
0080   double min_phi = phi - (M_PI/128.)*N - (M_PI/256.);
0081   double max_phi = phi + (M_PI/128.)*N + (M_PI/256.);
0082   if ( calo_phi<min_phi || calo_phi>max_phi ) { return false; }
0083   if ( calo_eta>=min_eta && calo_eta<=max_eta && calo_phi>=min_phi && calo_phi<=max_phi ) { return true; }
0084   return false;
0085 };
0086 
0087 bool tower_in_ring(float calo_eta, float calo_phi, float eta, float phi, int N){
0088   double n_rings = (double) N;
0089   double eta_lo = eta - 0.024*n_rings;
0090   double eta_hi = eta + 0.024*n_rings;
0091   double phi_lo = phi_in_range(phi - (M_PI/128.)*n_rings);// - (M_PI/256.);
0092   double phi_hi = phi_in_range(phi + (M_PI/128.)*n_rings);// + (M_PI/256.);
0093 //  if (fabs(delta_phi(calo_phi,phi_lo))>(M_PI/256.) || fabs(delta_phi(calo_phi,phi_hi))>(M_PI/256.)) { return false; }
0094 //  if ( fabs(eta_lo-calo_eta)>0.012 || fabs(eta_hi-calo_eta)>0.012 ) { return false; }
0095   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.)) {
0096     if ( calo_eta>eta_lo-0.012 && calo_eta<eta_hi+0.012 ) { return true; }
0097   }
0098   if ( fabs(eta_lo-calo_eta)<=0.012 || fabs(eta_hi-calo_eta)<=0.012 ) {
0099     if ( calo_phi>phi_lo-(M_PI/256.) && calo_phi<phi_hi+(M_PI/256.) ) { return true; }
0100   }
0101   return false;
0102 };
0103 
0104 const int i_cal = 2; //  THIS IS AN EMCAL STUDY
0105 //const int ref_cal = -1; //  STUDY WRT PION
0106 //const int ref_cal = 0; //  STUDY WRT IHCAL
0107 //const int ref_cal = 1; //  STUDY WRT OHCAL
0108 const int ref_cal = 2; //  STUDY WRT EMCAL
0109 
0110 void EMCalTowerRingStudy() {
0111   
0112   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.);
0113   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.);
0114   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.);
0115 
0116   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.);
0117   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.);
0118 //  hTotalCaloSum_per_nRings->GetXaxis()->SetNdivisions(29);
0119   for (int i=0; i<30; ++i){ hTotalCaloSum_per_nRings->GetXaxis()->SetBinLabel(i+1,Form("%ix%i",2*i+1,2*i+1)); }
0120   hTotalCaloSum_per_nRings->GetXaxis()->SetAlphanumeric();
0121 
0122   TFile *f = new TFile(inFileName,"READ");
0123   TTree *fChain = (TTree*)f->Get("T");
0124   
0125   fChain->SetBranchAddress("event", &event);
0126   fChain->SetBranchAddress("eta", &eta);
0127   fChain->SetBranchAddress("phi", &phi);
0128   fChain->SetBranchAddress("e", &e);
0129   fChain->SetBranchAddress("pt", &pt);
0130   fChain->SetBranchAddress("vx", &vx);
0131   fChain->SetBranchAddress("vy", &vy);
0132   fChain->SetBranchAddress("vz", &vz);
0133   
0134   for (int ical=0; ical<ncal; ++ical) {  // LOOP OVER CALORIMETER LAYERS
0135     fChain->SetBranchAddress(Form("nTow%s",cal_tag[ical].c_str()), &nTowers[ical]);
0136     fChain->SetBranchAddress(Form("eta%s",cal_tag[ical].c_str()), &calo_eta[ical]);
0137     fChain->SetBranchAddress(Form("phi%s",cal_tag[ical].c_str()), &calo_phi[ical]);
0138     fChain->SetBranchAddress(Form("e%s",cal_tag[ical].c_str()), &calo_e[ical]);
0139     fChain->SetBranchAddress(Form("ieta%s",cal_tag[ical].c_str()), &calo_ieta[ical]);
0140     fChain->SetBranchAddress(Form("iphi%s",cal_tag[ical].c_str()), &calo_iphi[ical]);
0141   }
0142   
0143   
0144   Long64_t nentries = fChain->GetEntriesFast();
0145   
0146   for (Long64_t jentry=0; jentry<nentries;jentry++) {  // LOOP OVER EVENTS
0147     Long64_t ientry = fChain->GetTreeNumber();
0148     if (ientry < 0) break;
0149     
0150     fChain->GetEntry(jentry);
0151     
0152     double total_calo_E[ncal] = {0.,0.,0.};
0153     
0154     float ref_e = 0.;
0155     float ref_eta;
0156     float ref_phi;
0157     int   ref_ieta;
0158     int   ref_iphi;
0159     
0160     if (ref_cal==-1) {
0161       ref_e = e;
0162       ref_eta = eta;
0163       ref_phi = phi;
0164 //      ref_phi = phi_in_range(phi+M_PI);
0165     }
0166     else if (ref_cal>=0 && ref_cal<=2){
0167       for (int itow=0; itow<calo_e[ref_cal]->size(); ++itow) {  // LOOP OVER REFERENCE CAL TOWERS
0168         if (calo_e[ref_cal]->at(itow)>ref_e) {
0169           ref_e = calo_e[ref_cal]->at(itow);
0170           ref_eta = calo_eta[ref_cal]->at(itow);
0171           ref_phi = calo_phi[ref_cal]->at(itow);
0172 //          ref_phi = phi_in_range( calo_phi[ref_cal]->at(itow) + M_PI );
0173         }
0174       }
0175     }
0176     else { cerr<<"ope"<<endl; }
0177     
0178     for (int index=0; index<2; ++index) {  // LOOP OVER HCALS
0179       for (int itow=0; itow<calo_e[index]->size(); ++itow) {  // LOOP OVER CAL TOWERS
0180         total_calo_E[index] += calo_e[index]->at(itow);
0181       }
0182     }
0183     
0184     int n_towers_in_square = 0;
0185     double square_e_sum = 0.;
0186     for (int n_rings = 0; n_rings<30; ++n_rings) {
0187       
0188       double ring_e_sum = 0.;
0189       int n_towers_in_ring = 0;
0190       for (int itow=0; itow<calo_e[i_cal]->size(); ++itow) {  // LOOP OVER EMCAL TOWERS
0191 //        if (i_cal==2 && calo_e[i_cal]->at(itow)<0.1) { continue; }
0192         if ( tower_in_ring(calo_eta[i_cal]->at(itow),calo_phi[i_cal]->at(itow),ref_eta,ref_phi,n_rings) ) {
0193           total_calo_E[2] += calo_e[i_cal]->at(itow);
0194           ring_e_sum += calo_e[i_cal]->at(itow);
0195           square_e_sum += calo_e[i_cal]->at(itow);
0196           n_towers_in_square += 1;
0197           n_towers_in_ring += 1;
0198         }
0199       }  // end tower loop
0200       hEsum_per_nRings->Fill(pt,square_e_sum,n_rings);
0201       hEsum_in_ringN->Fill(pt,ring_e_sum,n_rings);
0202       hEMCalE_over_pionE->Fill(pt,ring_e_sum/e,n_rings);
0203       
0204       hTotalCaloSum_per_nRings->Fill(e,square_e_sum+total_calo_E[0]+total_calo_E[1],n_rings);
0205       hTotalCaloEfrac_per_nRings->Fill(e,(square_e_sum+total_calo_E[0]+total_calo_E[1])/e,n_rings);
0206     }
0207     
0208     // if (Cut(ientry) < 0) continue;
0209     
0210   }  // end event loop
0211   
0212   
0213   
0214   string outputroot = (string) inFileName;
0215   string remove_this = ".root";
0216   size_t pos = outputroot.find(remove_this);
0217   if (pos != string::npos) { outputroot.erase(pos, remove_this.length()); }
0218   if (ref_cal>=0 && ref_cal<=2){
0219     outputroot += cal_tag[ref_cal] + "Ref";
0220   }
0221   TString outFileName = outputroot + "_EMCalRings.root";
0222   
0223   TFile *outFile = new TFile(outFileName,"RECREATE");
0224   hEsum_per_nRings->Write();
0225   hEsum_in_ringN->Write();
0226   hTotalCaloSum_per_nRings->Write();
0227   hTotalCaloEfrac_per_nRings->Write();
0228   hEMCalE_over_pionE->Write();
0229   //  for (int ical=0; ical<ncal; ++ical) {  // LOOP OVER CALORIMETER LAYERS
0230   //  }
0231 
0232   vector<TH3D *> histos = {hEsum_per_nRings, hEsum_in_ringN, hTotalCaloSum_per_nRings, hTotalCaloEfrac_per_nRings, hEMCalE_over_pionE};
0233   auto fa1 = new TF1("fa1","0.",0,1);
0234 
0235   TCanvas * can = new TCanvas( "can" , "" ,700 ,500 );
0236   can->SetLogz();
0237   const char us[] = "_";
0238   
0239   for (int i=0; i<histos.size(); ++i ) {
0240     histos[i]->Project3D("YZ")->Draw("COLZ");
0241     can->SaveAs(Form("plots/EMCalTowerRingStudy/%s.pdf",histos[i]->GetName()),"PDF");
0242     for (int j=0; j<3; ++j) {
0243       double ptlo = (double)10.*j;
0244       double pthi = (double)10.*(j+1);
0245       histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
0246       histos[i]->Project3D("YZ")->Draw("COLZ");
0247       can->SaveAs(Form("plots/EMCalTowerRingStudy/%s%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0248     }
0249     if (i==4){
0250       histos[i]->GetXaxis()->SetRangeUser(1.,2.);
0251       TH2D *hDenom2D = (TH2D*)histos[i]->Project3D("YZ");
0252       TH1D *hDenom = (TH1D*)hDenom2D->ProfileX();
0253       hDenom->SetName("hDenom");
0254 
0255       for (int j=1; j<10; ++j) {
0256         double ptlo = (double)1.*j;
0257         double pthi = (double)1.*(j+1);
0258         histos[i]->GetXaxis()->SetRangeUser(ptlo,pthi);
0259         histos[i]->Project3D("YZ")->Draw("COLZ");
0260         TH2D *hTemp = (TH2D*)histos[i]->Project3D("YZ");
0261         hTemp->SetLineColor(kRed);
0262         hTemp->SetMarkerColor(kRed);
0263         TH1D *hNum = (TH1D*)hTemp->ProfileX();
0264         hNum->Draw("SAME");
0265         can->SaveAs(Form("plots/EMCalTowerRingStudy/%s_FINE_%s%i%s%i.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0266         hNum->Divide(hDenom);
0267         hNum->GetYaxis()->SetRangeUser(0.1,2.2);
0268         hNum->Draw();
0269         fa1->Draw("SAME");
0270         can->SaveAs(Form("plots/EMCalTowerRingStudy/%s_FINE_%s%i%s%i__ratioTo1_2.pdf",histos[i]->GetName(),us,(int)ptlo,us,(int)pthi),"PDF");
0271         hNum->Reset();
0272       }
0273     }
0274     histos[i]->GetXaxis()->SetRangeUser(0.,30.);
0275   }
0276 
0277 //  hEsum_per_nRings->Project3D("YZ")->Draw("COLZ");
0278 //  can->SaveAs("plots/EMCalTowerRingStudy/hEsum_per_nRings.pdf","PDF");
0279 //  hEsum_in_ringN->Project3D("YZ")->Draw("COLZ");
0280 //  can->SaveAs("plots/EMCalTowerRingStudy/hEsum_in_ringN.pdf","PDF");
0281 //  hTotalCaloSum_per_nRings->Project3D("YZ")->Draw("COLZ");
0282 //  can->SaveAs("plots/EMCalTowerRingStudy/hTotalCaloSum_per_nRings.pdf","PDF");
0283 //  hTotalCaloEfrac_per_nRings->Project3D("YZ")->Draw("COLZ");
0284 //  can->SaveAs("plots/EMCalTowerRingStudy/hTotalCaloEfrac_per_nRings.pdf","PDF");
0285 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!