Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:19 

0001 #include "sPhenixStyle.h"
0002 #include "sPhenixStyle.C"
0003 
0004 
0005 void Jet_reso() 
0006 {
0007   SetsPhenixStyle();
0008   TH1::SetDefaultSumw2();
0009   TH2::SetDefaultSumw2();
0010 
0011   TChain * ct = new TChain("T");
0012   //if you want to run one file use this
0013   ct->Add("../macro/output.root");
0014 
0015   //if you want to combine multiple files use this
0016   /*for(int i = 0; i < 100; i++){
0017     ct->Add(Form("../macro/condor/output_%i.root",i));
0018     }*/
0019 
0020   vector<float> *eta = 0;
0021   vector<float> *phi = 0;
0022   vector<float> *pt = 0;
0023   vector<float> *e = 0;
0024   vector<float> *subtracted_et = 0;
0025   vector<float> *truthEta = 0;
0026   vector<float> *truthPhi = 0;
0027   vector<float> *truthPt = 0;
0028   vector<float> *truthE = 0;
0029   vector<float> *subseedEta = 0;
0030   vector<float> *subseedPhi = 0;
0031   vector<float> *subseedPt = 0;
0032   vector<float> *subseedE = 0;
0033   vector<int> *subseedCut = 0;
0034   vector<float> *rawseedEta = 0;
0035   vector<float> *rawseedPhi = 0;
0036   vector<float> *rawseedPt = 0;
0037   vector<float> *rawseedE = 0;
0038   vector<int> *rawseedCut = 0;
0039   float *totalCalo = 0;
0040   int cent;
0041 
0042   ct->SetBranchAddress("eta",&eta);
0043   ct->SetBranchAddress("phi",&phi);
0044   ct->SetBranchAddress("pt",&pt);
0045   ct->SetBranchAddress("e",&e);
0046   ct->SetBranchAddress("subtracted_et",&subtracted_et);
0047   ct->SetBranchAddress("truthEta",&truthEta);
0048   ct->SetBranchAddress("truthPhi",&truthPhi);
0049   ct->SetBranchAddress("truthPt",&truthPt);
0050   ct->SetBranchAddress("truthE",&truthE);
0051   ct->SetBranchAddress("cent", &cent);
0052 
0053   ct->SetBranchAddress("rawseedEta", &rawseedEta);
0054   ct->SetBranchAddress("rawseedPhi", &rawseedPhi);
0055   ct->SetBranchAddress("rawseedPt", &rawseedPt);
0056   ct->SetBranchAddress("rawseedE", &rawseedE);
0057   ct->SetBranchAddress("rawseedCut", &rawseedCut);
0058   ct->SetBranchAddress("subseedEta", &subseedEta);
0059   ct->SetBranchAddress("subseedPhi", &subseedPhi);
0060   ct->SetBranchAddress("subseedPt", &subseedPt);
0061   ct->SetBranchAddress("subseedE", &subseedE);
0062   ct->SetBranchAddress("subseedCut", &subseedCut);
0063 
0064   const Float_t pt_range[] = {10,15,20,25,30,35,40,45,50,60,80};
0065   const int pt_N = sizeof(pt_range)/sizeof(float) - 1;
0066   const int resp_N = 60;
0067   Float_t resp_bins[resp_N+1];
0068   for(int i = 0; i < resp_N+1; i++){
0069     resp_bins[i] = 1.2/resp_N * i;
0070   }
0071   const int eta_N = 40;
0072   Float_t eta_bins[eta_N+1];
0073   for(int i = 0; i < eta_N+1; i++){
0074     eta_bins[i] = -1.1 + 2.2/eta_N * i;
0075   }
0076   const int phi_N = 40;
0077   Float_t phi_bins[phi_N+1];
0078   for(int i = 0; i < phi_N+1; i++){
0079     phi_bins[i] = -TMath::Pi() + 2*TMath::Pi()/phi_N * i;
0080   }
0081   const int subet_N = 400;
0082   Float_t subet_bins[subet_N+1];
0083   for(int i = 0; i < subet_N+1; i++){
0084     subet_bins[i] = i*0.5;
0085   }
0086   const float cent_bins[] = {-1, 0, 10, 30, 50, 80}; //the first bin is for inclusive in centrality/pp events
0087   const int cent_N = sizeof(cent_bins)/sizeof(float) - 1;
0088 
0089   TH3F *h_MC_Reso_pt = new TH3F("h_MC_Reso", "" , pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);;
0090   h_MC_Reso_pt->GetXaxis()->SetTitle("p_{T}^{truth} [GeV]");
0091   h_MC_Reso_pt->GetYaxis()->SetTitle("p_{T}^{reco}/p_{T}^{truth}");
0092 
0093   TH2F *h_pt_reco = new TH2F("h_pt_reco","",subet_N,subet_bins, cent_N, cent_bins);
0094   h_pt_reco->GetXaxis()->SetTitle("p_{T} [GeV]");
0095   TH2F *h_pt_true = new TH2F("h_pt_true","",pt_N,pt_range, cent_N, cent_bins);
0096   h_pt_true->GetXaxis()->SetTitle("p_{T} [GeV]");
0097   TH2F *h_pt_true_matched = new TH2F("h_pt_true_matched","",subet_N,subet_bins, cent_N, cent_bins);
0098   h_pt_true_matched->GetXaxis()->SetTitle("p_{T} [GeV]");
0099   TH3F *h_eta_phi_reco = new TH3F("h_eta_phi_reco","",eta_N, eta_bins, phi_N, phi_bins, cent_N, cent_bins);
0100   h_eta_phi_reco->GetXaxis()->SetTitle("#eta");
0101   h_eta_phi_reco->GetYaxis()->SetTitle("#phi");
0102   TH3F *h_eta_phi_true = new TH3F("h_eta_phi_true","",eta_N, eta_bins, phi_N, phi_bins, cent_N, cent_bins);
0103   TH3F *h_subEt_pt = new TH3F("h_subEt_pt","",pt_N,pt_range, subet_N, subet_bins, cent_N, cent_bins);
0104 
0105 
0106   TGraphErrors *g_jes[cent_N];
0107   TGraphErrors *g_jer[cent_N];
0108   for(int icent = 0; icent < cent_N; icent++){
0109     g_jes[icent] = new TGraphErrors(pt_N);
0110     g_jer[icent] = new TGraphErrors(pt_N);
0111   }
0112 
0113   int nentries = ct->GetEntries();
0114   for(int i = 0; i < nentries; i++){
0115     ct->GetEntry(i);
0116 
0117     int njets = truthPt->size();
0118     int nrecojets = pt->size();
0119 
0120     for(int tj = 0; tj < njets; tj++){
0121       //fill truth hists
0122       h_pt_true->Fill(truthPt->at(tj), cent);
0123       h_pt_true->Fill(truthPt->at(tj), -1);
0124       h_eta_phi_true->Fill(truthEta->at(tj),truthPhi->at(tj), cent);
0125       h_eta_phi_true->Fill(truthEta->at(tj),truthPhi->at(tj), -1);
0126       //do reco to truth jet matching
0127       float matchEta, matchPhi, matchPt, matchE, matchsubtracted_et, dR;
0128       float dRMax = 100;
0129       for(int rj = 0; rj < nrecojets; rj++){
0130     float dEta = truthEta->at(tj) - eta->at(rj);
0131     float dPhi = truthPhi->at(tj) - phi->at(rj);
0132     while(dPhi > TMath::Pi()) dPhi -= 2*TMath::Pi();
0133     while(dPhi < -TMath::Pi()) dPhi += 2*TMath::Pi();
0134     dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0135     if(dR < dRMax){
0136       matchEta = eta->at(rj);
0137       matchPhi = phi->at(rj);
0138       matchPt = pt->at(rj);
0139       matchE = e->at(rj);
0140       matchsubtracted_et = subtracted_et->at(rj);
0141       dRMax = dR;
0142     }
0143       }
0144 
0145       if(dRMax > 0.3) continue;
0146       h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt/truthPt->at(tj), cent);
0147       h_pt_true_matched->Fill(truthPt->at(tj), cent);
0148       h_eta_phi_reco->Fill(matchEta,matchPhi, cent);
0149       h_pt_reco->Fill(matchPt, cent);
0150       h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt/truthPt->at(tj), -1);
0151       h_pt_true_matched->Fill(truthPt->at(tj), -1);
0152       h_eta_phi_reco->Fill(matchEta,matchPhi, -1);
0153       h_pt_reco->Fill(matchPt, -1);
0154       h_subEt_pt->Fill(matchPt, matchsubtracted_et, cent);
0155       h_subEt_pt->Fill(matchPt, matchsubtracted_et, -1);
0156     }
0157   }
0158 
0159   TCanvas *c = new TCanvas("c","c");
0160   //c->Print("fits04.pdf("); //uncomment these to get a pdf with all the fits
0161   TLegend *leg = new TLegend(.25,.2,.6,.5);
0162   leg->SetFillStyle(0);
0163   gStyle->SetOptFit(1);
0164   for(int icent = 0; icent < cent_N; ++icent){
0165     for (int i = 0; i < pt_N; ++i){
0166       TF1 *func = new TF1("func","gaus",0,1.2);
0167       h_MC_Reso_pt->GetXaxis()->SetRange(i+1,i+1);
0168       h_MC_Reso_pt->GetZaxis()->SetRange(icent+1,icent+1);
0169       TH1F *h_temp = (TH1F*) h_MC_Reso_pt->Project3D("y");
0170       h_temp->Fit(func,"","",0,1.2);
0171       h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
0172       func->SetLineColor(kRed);
0173       h_temp->Draw();
0174       leg->AddEntry("",Form("%2.0f%%-%2.0f%%",h_MC_Reso_pt->GetZaxis()->GetBinLowEdge(icent+1),h_MC_Reso_pt->GetZaxis()->GetBinLowEdge(icent+2)),"");
0175       leg->AddEntry("",Form("%2.0f < p_T < %2.0f GeV",h_MC_Reso_pt->GetXaxis()->GetBinLowEdge(i+1),h_MC_Reso_pt->GetXaxis()->GetBinLowEdge(i+2)),"");
0176       leg->Draw("SAME");
0177       /*-------for calculating the JER uncertainty----*/
0178       float dsigma = func -> GetParError(2);
0179       float denergy = func -> GetParError(1);
0180       float sigma = func -> GetParameter(2);
0181       float energy = func -> GetParameter(1);
0182 
0183       float djer = dsigma/energy + sigma*denergy/pow(energy,2);//correct way to calculate jer
0184                                                                //accounting for the fact that the 
0185                                                                //mean response and width are correlated
0186       //c->Print("fits04.pdf");
0187       leg->Clear();
0188       g_jes[icent]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(1));
0189       g_jes[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),func->GetParError(1));
0190       g_jer[icent]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(2)/func->GetParameter(1));
0191       
0192       //g_jer[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),func->GetParError(2));
0193       g_jer[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),djer);
0194     }
0195   }
0196   //c->Print("fits04.pdf)");
0197 
0198 
0199 
0200   TFile *f_out = new TFile("hists.root","RECREATE");
0201   for(int icent = 0; icent < cent_N; ++icent){
0202     g_jes[icent]->Write(Form("g_jes_cent%i",icent));
0203     g_jer[icent]->Write(Form("g_jer_cent%i",icent));
0204   }
0205 
0206   h_pt_true->Write();
0207   h_eta_phi_true->Write();
0208   h_pt_true_matched->Write();
0209   h_eta_phi_reco->Write();
0210   h_pt_reco->Write();
0211   h_subEt_pt->Write();
0212 }
0213
[ 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!