Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:10 

0001 #include "sPhenixStyle.h"
0002 #include "sPhenixStyle.C"
0003 
0004 bool isIso(int thisRjet, int thisTjet, vector<float> *rJetpT, vector<float> *tJetpT, vector<float> *rJetEta, vector<float> *tJetEta, vector<float> *rJetPhi, vector<float> *tJetPhi, float minDr, TH1F *&rDist, TH1F *&tDist);
0005 // This macro analyzes output from the jet validation sub module of the 
0006 // JS-Jet module. To first order, it simply computes the jet energy scale (JES)
0007 // and jet energy resolution (JER), but also compute the jet energy linearity
0008 // necessary for the numerical inversion which is used to calibrate the MC
0009 // to itself.
0010   
0011 void Jet_reso_Iso(float jetR = 4, float isoParam = 1, int applyCorr = 0, int isLin = 0) 
0012 {
0013   //jetR: jet radius
0014   //
0015   //isoParameter: used to compute distance required to consider a jet "isolated" 
0016   //and is multipled by the jet radius. See minDist variable below. 
0017   //
0018   //applyCorr: determines whether or not to apply calibrations. 
0019   //
0020   //isLin: measure the jet linearity instead of the response
0021   float minDist = .1*(jetR*isoParam);
0022   if(minDist)std::cout << "isolation requirement: dR > " << minDist << " [rad]" << std::endl;
0023   
0024   SetsPhenixStyle();
0025   TH1::SetDefaultSumw2();
0026   TH2::SetDefaultSumw2();
0027   TH3::SetDefaultSumw2();
0028   
0029   TChain * ct = new TChain("T");
0030   //ct->Add("run40_30GeV_10GeV_Spliced_RecoJets_test.root");  
0031   ct->Add("run40_30GeV_10GeV_Spliced_RecoJets.root");
0032   //if you want to combine multiple files use this
0033   //for(int i = 0; i < 15000; i++){
0034   //  ct->Add(Form("../macro/condor/output/jetIso_%d/isoJetCalibOut_%d.root",i,i));
0035   // }
0036   
0037   TH1F *recoDist = new TH1F("recoDist","distance between reco jets",200,0,5);
0038   TH1F *truthDist = new TH1F("truthDist","distance between truth jets",200,0,5);
0039 
0040   TFile *f_out = new TFile(Form("hists_R0%g_dR%g_Corr%d_isLin%d.root",jetR,isoParam,applyCorr,isLin),"RECREATE");
0041   
0042   vector<float> *eta = 0;
0043   vector<float> *phi = 0;
0044   vector<float> *pt = 0;
0045   //vector<float> *e = 0;
0046   //vector<float> *subtracted_et = 0;
0047   vector<float> *truthEta = 0;
0048   vector<float> *truthPhi = 0;
0049   vector<float> *truthPt = 0;
0050   //vector<float> *truthE = 0;
0051 
0052   vector<float> *subseedEta = 0;
0053   vector<float> *subseedPhi = 0;
0054   vector<float> *subseedPt = 0;
0055   vector<float> *subseedE = 0;
0056   vector<int> *subseedCut = 0;
0057   vector<float> *rawseedEta = 0;
0058   vector<float> *rawseedPhi = 0;
0059   vector<float> *rawseedPt = 0;
0060   //vector<float> *rawseedE = 0;
0061   //vector<int> *rawseedCut = 0;
0062   float *totalCalo = 0;
0063   int cent;
0064   Float_t rWeight, tWeight, RawWeight, SubWeight;
0065   
0066 
0067   ct->SetBranchAddress("eta",&eta);
0068   ct->SetBranchAddress("phi",&phi);
0069   ct->SetBranchAddress("pt",&pt);
0070   //ct->SetBranchAddress("e",&e);
0071   //ct->SetBranchAddress("subtracted_et",&subtracted_et);
0072   ct->SetBranchAddress("truthEta",&truthEta);
0073   ct->SetBranchAddress("truthPhi",&truthPhi);
0074   ct->SetBranchAddress("truthPt",&truthPt);
0075   //ct->SetBranchAddress("truthE",&truthE);
0076   //ct->SetBranchAddress("cent", &cent);
0077 
0078   ct->SetBranchAddress("rawseedEta", &rawseedEta);
0079   ct->SetBranchAddress("rawseedPhi", &rawseedPhi);
0080   ct->SetBranchAddress("rawseedPt", &rawseedPt);
0081   //ct->SetBranchAddress("rawseedE", &rawseedE);
0082   //ct->SetBranchAddress("rawseedCut", &rawseedCut);
0083   ct->SetBranchAddress("subseedEta", &subseedEta);
0084   ct->SetBranchAddress("subseedPhi", &subseedPhi);
0085   ct->SetBranchAddress("subseedPt", &subseedPt);
0086   //ct->SetBranchAddress("subseedE", &subseedE);
0087   //ct->SetBranchAddress("subseedCut", &subseedCut);
0088 
0089   ct->SetBranchAddress("tWeight",&tWeight);
0090   ct->SetBranchAddress("rWeight",&rWeight);
0091   ct->SetBranchAddress("RawWeight",&RawWeight);
0092   ct->SetBranchAddress("SubWeight",&SubWeight);
0093   
0094   //const Float_t pt_range[]  = {10,15,20,25,30,35,40,45,50,60,80};
0095   const Float_t pt_range[]  = {15,20,25,30,35,40,45,50,60,80};
0096 
0097   //if(applyCorr)  pt_range[] = {10,15,20,25,30,35,40,45,50,60,80};
0098   //if(isLin)const Float_t pt_range[] = {10,15,20,25,30,35,40,45,50,55,60,65,70,75,80};
0099    
0100   const int pt_N = sizeof(pt_range)/sizeof(float) - 1;
0101   //const int pt_N_truth = 34;//for including the 10-15 region
0102   const int pt_N_truth = 37;
0103   const int resp_N = 100;
0104   Float_t resp_bins[resp_N];
0105   for(int i = 0; i < resp_N + 1; i++){
0106     //resp_bins[i] = 1.2/resp_N * i;
0107     resp_bins[i] = 2./resp_N * i;
0108   }
0109 
0110 
0111   // Float_t pt_range_reco[pt_N_reco];
0112   //const int pt_N_reco = 71;
0113   // for(int i = 0; i < pt_N_reco + 1; i++)
0114   //   {
0115   //     pt_range_reco[i] = 80./pt_N_reco * i;
0116   //   }
0117   vector<Float_t> pt_range_reco1;
0118   Float_t ptBin = 0.;
0119   // while(ptBin < 80)
0120   //   {
0121   //     pt_range_reco1.push_back(ptBin);
0122   //     if(ptBin <= 20)ptBin += 0.1;
0123   //     else if(ptBin >20 && ptBin <= 30) ptBin += 1;
0124   //     else if(ptBin > 30 && ptBin < 40) ptBin +=0.1;
0125   //     else ptBin += 2.5;
0126   //   }
0127   while(ptBin < 79.9)
0128     {
0129       pt_range_reco1.push_back(ptBin);
0130       ptBin += 0.1;
0131     }
0132   
0133   const int pt_N_reco = pt_range_reco1.size();//sizeof(pt_range_reco1)/sizeof(float) - 1;
0134   Float_t pt_range_reco[pt_N_reco];
0135   for(int i = 0; i < pt_N_reco+1; i++)
0136     {
0137       if(i < pt_N_reco)pt_range_reco[i] = pt_range_reco1.at(i);
0138       else pt_range_reco[i] = 79.9;
0139     }
0140   Float_t pt_range_truth[pt_N_truth];
0141   Float_t rebin[pt_N_truth];
0142   for(int i = 0; i < pt_N_truth+1; i++)
0143   {
0144     //pt_range_truth[i] = 10+70./pt_N_truth * i;
0145     //pt_range_truth[i] = 10+2*i;
0146     pt_range_truth[i] = 15+1.5*i;
0147     //if( pt_range_truth[i] < 15) rebin[i] = 1;
0148     if(pt_range_truth[i] >= 15 && pt_range_truth[i] < 30) rebin[i] = 4;//4;
0149     else if(pt_range_truth[i] >= 30 && pt_range_truth[i] < 40) rebin[i] = 1;
0150     else rebin[i] = 10;
0151   }
0152   
0153   
0154   const int eta_N = 40;
0155   Float_t eta_bins[eta_N];
0156   for(int i = 0; i < eta_N+1; i++){
0157     eta_bins[i] = -1.1 + 2.2/eta_N * i;
0158   }
0159   const int phi_N = 40;
0160   Float_t phi_bins[phi_N];
0161   for(int i = 0; i < phi_N+1; i++){
0162     phi_bins[i] = -TMath::Pi() + 2*TMath::Pi()/phi_N * i;
0163   }
0164   const int subet_N = 400;
0165   Float_t subet_bins[subet_N];
0166   for(int i = 0; i < subet_N+1; i++){
0167     subet_bins[i] = i*0.2;
0168   }
0169   const float cent_bins[] = {-1, 0, 10, 30, 50, 80}; //the first bin is for inclusive in centrality/pp events
0170   const int cent_N = 1;//sizeof(cent_bins)/sizeof(float) - 1;
0171 
0172   TH3F *h_MC_Reso_pt = 0;
0173   if(isLin) h_MC_Reso_pt= new TH3F("h_MC_Reso", "" , pt_N_truth, pt_range_truth, pt_N_reco, pt_range_reco, cent_N, cent_bins);
0174   else  h_MC_Reso_pt= new TH3F("h_MC_Reso", "" , pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
0175 
0176   h_MC_Reso_pt -> GetXaxis()->SetTitle("p_{T}^{truth} [GeV]");
0177   if(isLin) h_MC_Reso_pt->GetYaxis()->SetTitle("p_{T}^{reco}");
0178   else h_MC_Reso_pt->GetYaxis()->SetTitle("p_{T}^{reco}/p_{T}^{truth}");
0179   TH2F *h_pt_reco = new TH2F("h_pt_reco","",subet_N,subet_bins, cent_N, cent_bins);
0180   h_pt_reco->GetXaxis()->SetTitle("p_{T} [GeV]");
0181   TH2F *h_pt_true = new TH2F("h_pt_true","",pt_N,pt_range, cent_N, cent_bins);
0182   h_pt_true->GetXaxis()->SetTitle("p_{T} [GeV]");
0183   TH2F *h_pt_true_matched = new TH2F("h_pt_true_matched","",subet_N,subet_bins, cent_N, cent_bins);
0184   h_pt_true_matched->GetXaxis()->SetTitle("p_{T} [GeV]");
0185   TH3F *h_eta_phi_reco = new TH3F("h_eta_phi_reco","",eta_N, eta_bins, phi_N, phi_bins, cent_N, cent_bins);
0186   h_eta_phi_reco->GetXaxis()->SetTitle("#eta");
0187   h_eta_phi_reco->GetYaxis()->SetTitle("#phi");
0188   TH3F *h_eta_phi_true = new TH3F("h_eta_phi_true","",eta_N, eta_bins, phi_N, phi_bins, cent_N, cent_bins);
0189   TH3F *h_subEt_pt = new TH3F("h_subEt_pt","",pt_N,pt_range, subet_N, subet_bins, cent_N, cent_bins);
0190 
0191 
0192   TGraphErrors *g_jes[cent_N];
0193   TGraphErrors *g_jer[cent_N];
0194   for(int icent = 0; icent < cent_N; icent++){
0195     g_jes[icent] = new TGraphErrors(pt_N);
0196     g_jer[icent] = new TGraphErrors(pt_N);
0197   }
0198 
0199   TFile *corrFile;
0200   TF1 *correction = new TF1("jet energy correction","1",0,80);
0201   if(applyCorr)
0202     {
0203       corrFile = new TFile("JES_IsoCorr_NumInv.root","READ");
0204       corrFile -> cd();
0205       correction = (TF1*)corrFile -> Get(Form("corrFit_Iso%g",isoParam));
0206     }
0207 
0208   int nentries = ct->GetEntries();
0209   for(int i = 0; i < nentries; i++){
0210     ct->GetEntry(i);
0211     cent = 101;
0212     int njets = truthPt->size();
0213     int nrecojets = pt->size();
0214     if(applyCorr)
0215       {
0216     for(int rj = 0; rj < nrecojets; rj++)
0217       {
0218         pt -> at(rj) *= correction -> Eval(pt -> at(rj));//correct all the jets in one go
0219       }
0220       }
0221     for(int tj = 0; tj < njets; tj++){
0222       //fill truth hists
0223       //if(truthPt -> at(tj) < 10) continue;
0224       h_pt_true->Fill(truthPt->at(tj), cent);
0225       h_pt_true->Fill(truthPt->at(tj), -1);
0226       h_eta_phi_true->Fill(truthEta->at(tj),truthPhi->at(tj), cent);
0227       h_eta_phi_true->Fill(truthEta->at(tj),truthPhi->at(tj), -1);
0228       //do reco to truth jet matching
0229       float matchEta, matchPhi, matchPt, matchE, matchsubtracted_et, dR;
0230       float dRMax = 100;
0231       int recoMatchJet = -9999;
0232       
0233       for(int rj = 0; rj < nrecojets; rj++){
0234     float dEta = truthEta->at(tj) - eta->at(rj);
0235     float dPhi = truthPhi->at(tj) - phi->at(rj);
0236     while(dPhi > TMath::Pi()) dPhi -= 2*TMath::Pi();
0237     while(dPhi < -TMath::Pi()) dPhi += 2*TMath::Pi();
0238     dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0239     
0240     if(dR < dRMax){
0241       matchEta = eta->at(rj);
0242       matchPhi = phi->at(rj);
0243       matchPt = pt->at(rj);
0244       //matchE = e->at(rj);
0245       //matchsubtracted_et = subtracted_et->at(rj);
0246       dRMax = dR;
0247       recoMatchJet = rj;
0248     }
0249       }
0250       
0251       if(dRMax > 0.1*0.75*jetR) continue;
0252       if(minDist && !isIso(recoMatchJet, tj, pt, truthPt, eta, truthEta, phi, truthPhi, minDist, recoDist, truthDist)) continue;
0253       //if(truthPt->at(tj) >= 30 && rWeight == 1) continue;
0254       if(!isLin)
0255     {//if we're applying the correction, we're looking at the JES, fill with ptReco/ptTruth
0256       
0257       h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt/truthPt->at(tj), cent, rWeight);
0258       h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt/truthPt->at(tj), -1, rWeight);
0259     }
0260       else 
0261     {//We're looking at the linearity, fill with just the matchp
0262       h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt, cent, rWeight);
0263       h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt, -1, rWeight);
0264     }
0265       
0266       h_pt_true_matched->Fill(truthPt->at(tj), cent, tWeight);
0267       h_eta_phi_reco->Fill(matchEta,matchPhi, cent);
0268       h_pt_reco->Fill(matchPt, cent, rWeight);
0269    
0270       h_pt_true_matched->Fill(truthPt->at(tj), -1, tWeight);
0271       h_eta_phi_reco->Fill(matchEta,matchPhi, -1);
0272       h_pt_reco->Fill(matchPt, -1, rWeight);
0273       //h_subEt_pt->Fill(matchPt, matchsubtracted_et, cent);
0274       //h_subEt_pt->Fill(matchPt, matchsubtracted_et, -1);
0275     }
0276   }
0277 
0278   TCanvas *c = new TCanvas("c","c");
0279   //c->Print("fits04.pdf("); //uncomment these to get a pdf with all the fits
0280   TLegend *leg = new TLegend(.6,.9,.9,.9);
0281   leg->SetFillStyle(0);
0282   gStyle->SetOptFit(0);
0283   f_out -> cd();
0284   int stop;//Sets the number of bins to look at, which are different depending on whether or not you're doing the calibration
0285   float fitEnd;//The fit is super generic so it really doesn't care if you're fitting from 0-2 (JES) or 0-80 (linearity), but just to be thorough
0286   float fitStart;
0287   if(isLin) 
0288     {
0289       stop = pt_N_truth;
0290       fitEnd = 80.;
0291       fitStart = .1;
0292     }
0293   else 
0294     {
0295       stop = pt_N;
0296       fitEnd = 2.;
0297       fitStart = 0.1;
0298     }
0299   for(int icent = 0; icent < cent_N; ++icent){
0300     for (int i = 0; i < stop; ++i){
0301       TF1 *func = new TF1("func","gaus",fitStart,fitEnd);
0302       h_MC_Reso_pt->GetXaxis()->SetRange(i+1,i+1);
0303       h_MC_Reso_pt->GetZaxis()->SetRange(icent+1,icent+1);
0304       TH1F *h_temp = (TH1F*) h_MC_Reso_pt->Project3D("y");
0305       if(isLin)h_temp -> Rebin(rebin[i]);
0306       func -> SetParameter(1, h_temp -> GetBinCenter(h_temp -> GetMaximumBin()));
0307 
0308       h_temp -> Fit(func,"","",fitStart,fitEnd);
0309       h_temp -> Write(Form("hFit_cent%d_pt%d_Fit0",icent+1,i));
0310       func -> SetParameter(1, h_temp -> GetBinCenter(h_temp -> GetMaximumBin()));
0311 
0312       h_temp->Fit(func,"","",func->GetParameter(1)-0.75*func->GetParameter(2),func->GetParameter(1)+0.75*func->GetParameter(2));
0313       if(isLin)h_temp -> GetXaxis() -> SetRangeUser(func->GetParameter(1)-4*func->GetParameter(2),func->GetParameter(1)+4*func->GetParameter(2));
0314       func -> SetLineColor(2);
0315       func -> SetRange(func->GetParameter(1)-func->GetParameter(2),func->GetParameter(1)+func->GetParameter(2));
0316       h_temp->Draw();
0317       func -> Draw("same");
0318       gPad -> SetLogy();
0319       func->SetLineColor(2);
0320 
0321       //gPad -> WaitPrimitive();
0322       h_temp -> Write(Form("hFit_cent%d_pt%d_Fit1",icent+1,i));
0323       //leg->AddEntry("",Form("%2.0f%%-%2.0f%%",h_MC_Reso_pt->GetZaxis()->GetBinLowEdge(icent+1),h_MC_Reso_pt->GetZaxis()->GetBinLowEdge(icent+2)),"");
0324       leg->AddEntry("",Form("%g < p_T < %g GeV",h_MC_Reso_pt->GetXaxis()->GetBinLowEdge(i+1),h_MC_Reso_pt->GetXaxis()->GetBinLowEdge(i+2)),"");
0325       leg->Draw("SAME");
0326       c->Print(Form("plots/hFit_cent%d_pt%d_Fit1_Corr%d_isLin%d_Param%g.pdf",icent+1,i,applyCorr,isLin,isoParam));
0327       leg->Clear();
0328       float dsigma = func->GetParError(2);
0329       float denergy = func->GetParError(1);
0330       float sigma = func->GetParameter(2);
0331       float energy = func->GetParameter(1);
0332       
0333       float djer = dsigma/energy + sigma*denergy/pow(energy,2);//correct way to compute error on the resolution.
0334       if(isLin)
0335     {
0336       g_jes[icent]->SetPoint(i,0.5*(pt_range_truth[i]+pt_range_truth[i+1]),func->GetParameter(1));
0337       g_jes[icent]->SetPointError(i,0.5*(pt_range_truth[i+1]-pt_range_truth[i]),func->GetParError(1));
0338       
0339       g_jer[icent]->SetPoint(i,0.5*(pt_range_truth[i]+pt_range_truth[i+1]),func->GetParameter(2)/func->GetParameter(1));
0340       g_jer[icent]->SetPointError(i,0.5*(pt_range_truth[i+1]-pt_range_truth[i]),djer);
0341     }
0342       else
0343     {
0344       g_jes[icent]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(1));
0345       g_jes[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),func->GetParError(1));
0346       
0347       g_jer[icent]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(2)/func->GetParameter(1));
0348       g_jer[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),djer);
0349     }
0350     }
0351   }
0352   //c->Print("fits04.pdf)");
0353 
0354 
0355 
0356   for(int icent = 0; icent < cent_N; ++icent){
0357     g_jes[icent]->Write(Form("g_jes_cent%i",icent));
0358     g_jer[icent]->Write(Form("g_jer_cent%i",icent));
0359   }
0360   
0361   recoDist -> Write();
0362   truthDist -> Write();
0363   
0364   h_pt_true->Write();
0365   h_eta_phi_true->Write();
0366   h_pt_true_matched->Write();
0367   h_eta_phi_reco->Write();
0368   h_pt_reco->Write();
0369   h_subEt_pt->Write();
0370   
0371   f_out -> Close();
0372 }
0373 
0374 bool isIso(int thisRjet, int thisTjet, vector<float> *rJetpT, vector<float> *tJetpT, vector<float> *rJetEta, vector<float> *tJetEta, vector<float> *rJetPhi, vector<float> *tJetPhi, float minDr, TH1F *&rDist, TH1F *&tDist)
0375 {
0376   if(thisRjet < 0) return false;
0377   float thisRJetEta = rJetEta->at(thisRjet);
0378   float thisTJetEta = tJetEta->at(thisTjet);
0379   float thisRJetPhi = rJetPhi->at(thisRjet);
0380   float thisTJetPhi = tJetPhi->at(thisTjet);
0381   float minPt = 3.5;
0382   float dr = -1;
0383   
0384   //uncomment these cout statements if you need to convince yourself this makes sense. 
0385   //std::cout << "On rJet: " << thisRjet << " with size " <<  rJetpT->size() << std::endl;
0386   for(int i = 0; i < rJetpT->size(); i++)
0387     {
0388       if(i == thisRjet) 
0389     {
0390       //std::cout << "On index: " << i << "; thisRJet: " << thisRjet << std::endl;
0391       continue;
0392     }
0393       if(rJetpT->at(i) < minPt)
0394     {
0395       //std::cout << "Discarding reco jet with pt: " << rJetpT->at(i) << std::endl;
0396       continue;
0397     }
0398       float deta = pow(thisRJetEta - rJetEta->at(i),2);
0399       float dphi = thisRJetPhi - rJetPhi->at(i);
0400       while(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
0401       while(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
0402       dphi = pow(dphi,2);
0403 
0404       dr = sqrt(deta + dphi);
0405      
0406       if(dr < minDr) return false;
0407     
0408     }
0409  
0410   if(dr)rDist -> Fill(dr);//dr will stay at -1 if above loop is totally bypassed, which usually happens becuase of the p_T filter. 
0411   dr = -1;
0412   for(int i = 0; i < tJetpT->size(); i++)
0413     {
0414       if(i == thisTjet) continue;
0415       //if(tJetpT < minPt) continue;
0416       
0417       float deta = pow(thisTJetEta - tJetEta->at(i),2);
0418       float dphi = pow(thisTJetPhi - tJetPhi->at(i),2);
0419       while(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
0420       while(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
0421       dphi = pow(dphi,2);
0422 
0423       dr = sqrt(deta + dphi);
0424       
0425       if(dr < minDr) return false;
0426     }     
0427   
0428   if(dr)tDist -> Fill(dr);
0429   
0430   return true;
0431 }
[ 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!