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 #include "respFitterDef.h"
0004 
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 int getEtaBin(float tJeteta, int nEtaBins, std::vector<float> etaBins);
0011 bool isInRange(float truthJetPt, float mcWeight);
0012 float getLeadJetPt(std::vector<float> *tJetPt);
0013 
0014 
0015 void Jet_Resp_Condor_histMaker(float jetR = 4, int applyCorr = 0, int isLin = 0, const char* fin= "", int segment = 0, int trig = 0) 
0016 {
0017   //jetR: jet radius
0018   //
0019   //applyCorr: determines whether or not to apply calibrations. 
0020   //
0021   //isLin: measure the jet linearity instead of the response
0022   //if variables appear undefined, see respFitterDef.h
0023   vector<float> etaBins;
0024   for(float i = etaStart; i <= etaEnd; i += etaGap)
0025     {
0026       etaBins.push_back(i);
0027     }
0028   const int nEtaBins = etaBins.size() - 1 + extraBins; //see response fitter defs 
0029   std::cout << "nEtaBins " << nEtaBins << std::endl;
0030   float maxDist = 0.1*0.75*jetR;
0031   
0032   SetsPhenixStyle();
0033   TH1::SetDefaultSumw2();
0034   TH2::SetDefaultSumw2();
0035 
0036   TChain * ct = new TChain("T");
0037   ct -> Add(fin);
0038   TH1F *nEvents = new TH1F("nEvents","nEvents",3,0.5,3.5);
0039   //if(trig == 0) nEvents -> SetBinContent(1,ct->GetEntries());
0040   //else if(trig == 10) nEvents -> SetBinContent(2,ct->GetEntries());
0041   // else nEvents -> SetBinContent(3,ct -> GetEntries());
0042   
0043   TH1F *recoDist = new TH1F("recoDist","distance between reco jets",200,0,5);
0044   TH1F *truthDist = new TH1F("truthDist","distance between truth jets",200,0,5);
0045 
0046   TFile *f_out = new TFile(Form("/gpfs/mnt/gpfs02/sphenix/user/ahodges/macros_git/analysis/JS-JetOrig/JetValidation/offline/respFitting/output/hists_R0%g_Corr%d_isLin%d_2D_trig%d_seg%d.root",jetR,applyCorr,isLin,trig,segment),"RECREATE");
0047   
0048   vector<float> *eta = 0;
0049   vector<float> *phi = 0;
0050   vector<float> *pt = 0;
0051   vector<float> *e = 0;
0052   vector<float> *subtracted_et = 0;
0053   vector<float> *truthEta = 0;
0054   vector<float> *truthPhi = 0;
0055   vector<float> *truthPt = 0;
0056   vector<float> *truthE = 0;
0057 
0058   vector<float> *subseedEta = 0;
0059   vector<float> *subseedPhi = 0;
0060   vector<float> *subseedPt = 0;
0061   vector<float> *subseedE = 0;
0062   vector<int> *subseedCut = 0;
0063   vector<float> *rawseedEta = 0;
0064   vector<float> *rawseedPhi = 0;
0065   vector<float> *rawseedPt = 0;
0066   vector<float> *rawseedE = 0;
0067   vector<int> *rawseedCut = 0;
0068   float *totalCalo = 0;
0069   int cent;
0070   Float_t rWeight, tWeight, RawWeight, SubWeight;
0071   float mcWeight;
0072   
0073   ct->SetBranchAddress("eta",&eta);
0074   ct->SetBranchAddress("phi",&phi);
0075   ct->SetBranchAddress("pt",&pt);
0076   ct->SetBranchAddress("e",&e);
0077   ct->SetBranchAddress("subtracted_et",&subtracted_et);
0078   ct->SetBranchAddress("truthEta",&truthEta);
0079   ct->SetBranchAddress("truthPhi",&truthPhi);
0080   ct->SetBranchAddress("truthPt",&truthPt);
0081   ct->SetBranchAddress("truthE",&truthE);
0082   ct->SetBranchAddress("cent", &cent);
0083   ct->SetBranchAddress("mcWeight",&mcWeight);
0084     
0085   ct->SetBranchAddress("rawseedEta", &rawseedEta);
0086   ct->SetBranchAddress("rawseedPhi", &rawseedPhi);
0087   ct->SetBranchAddress("rawseedPt", &rawseedPt);
0088   ct->SetBranchAddress("rawseedE", &rawseedE);
0089   ct->SetBranchAddress("rawseedCut", &rawseedCut);
0090   ct->SetBranchAddress("subseedEta", &subseedEta);
0091   ct->SetBranchAddress("subseedPhi", &subseedPhi);
0092   ct->SetBranchAddress("subseedPt", &subseedPt);
0093   ct->SetBranchAddress("subseedE", &subseedE);
0094   ct->SetBranchAddress("subseedCut", &subseedCut);
0095 
0096   
0097   float resp_bins[resp_N+1];
0098   for(int i = 0; i < resp_N+1; i++){
0099     resp_bins[i] = 2.0/resp_N * i;
0100   }
0101 
0102   Float_t pt_range_reco[pt_N_reco+1];
0103   for(int i = 0; i < pt_N_reco+1; i++)
0104     {
0105       pt_range_reco[i] = 80./pt_N_reco*i;
0106     }
0107   
0108   Float_t pt_range_truth[pt_N_truth];
0109   for(int i = 0; i < pt_N_truth+1; i++)
0110     {
0111       //pt_range_truth[i] = 10+70./pt_N_truth * i;
0112       pt_range_truth[i] = 5+75./pt_N_truth * i;
0113     }
0114 
0115 
0116   
0117   TH2D *h_MC_Reso_pt[nEtaBins];
0118 
0119   if(isLin) 
0120     {
0121       for(int i = 0; i < nEtaBins; i++)
0122         {
0123           h_MC_Reso_pt[i] = new TH2D(Form("h_MC_Reso_eta%d",i), "" , pt_N_truth, pt_range_truth, pt_N_reco, pt_range_reco);
0124           h_MC_Reso_pt[i] -> GetXaxis()->SetTitle("p_{T}^{truth} [GeV]");
0125           h_MC_Reso_pt[i] -> GetYaxis()->SetTitle("p_{T}^{reco}");
0126         }
0127     }
0128   else
0129     {
0130       for(int i = 0; i < nEtaBins; i++) 
0131     {
0132       h_MC_Reso_pt[i]= new TH2D(Form("h_MC_Reso_eta%d",i),Form("h_MC_Reso_eta%d",i) , pt_N, pt_range, resp_N, resp_bins);
0133       h_MC_Reso_pt[i] -> GetXaxis()->SetTitle("p_{T}^{truth} [GeV]");
0134       h_MC_Reso_pt[i] -> GetYaxis()->SetTitle("p_{T}^{reco}/p_{T}^{truth}");
0135     }
0136     }
0137 
0138  
0139   TH1F *h_pt_true_lead = new TH1F("h_pt_true_lead","",pt_N,pt_range);
0140   h_pt_true_lead->GetXaxis()->SetTitle("p_{T} [GeV]");
0141   TH1F *h_pt_true_lead_preCut = new TH1F("h_pt_true_lead_preCut","",50,0,100);
0142   h_pt_true_lead_preCut->GetXaxis()->SetTitle("p_{T} [GeV]");
0143 
0144   
0145   Float_t subet_bins[subet_N+1];
0146   for(int i = 0; i < subet_N+1; i++){
0147     subet_bins[i] = i*0.5;
0148   }
0149 
0150   TH1F *h_pt_true_matched[nEtaBins-1];
0151   for(int i = 0; i < nEtaBins-1; i++)
0152     {
0153       h_pt_true_matched[i] = new TH1F(Form("h_pt_true_matched_eta%d",i),"",subet_N,subet_bins);
0154       h_pt_true_matched[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
0155     }
0156   
0157   
0158   TFile *corrFile = new TFile("JES_IsoCorr_NumInv.root");
0159   TF1 *correction[nEtaBins-1];//last bin doesn't have a unique calibration
0160   TF1 *correctionExtrap[nEtaBins-1];
0161   for(int i = 0; i < nEtaBins-1; i++) correction[i] = new TF1(Form("jes_Corr_eta%d",i),"1",0,80);
0162   if(applyCorr && corrFile)
0163     {
0164       corrFile -> cd();
0165       for(int i = 0; i<nEtaBins-1; i++)correction[i] = (TF1*)corrFile -> Get(Form("corrFit_eta%d",i));
0166       for(int i = 0; i<nEtaBins-1; i++)correctionExtrap[i] = (TF1*)corrFile -> Get(Form("corrFit_eta%d_extrap",i));
0167     }
0168 
0169   int nentries = ct->GetEntries();
0170   std::cout << "about to run over " << nentries << " events" << std::endl;
0171   float centileIterator = 0.1;
0172   for(int i = 0; i < nentries; i++){
0173     ct->GetEntry(i);
0174     if(i == (int)floor(nentries*centileIterator))
0175       {
0176     std::cout << centileIterator*100.  << "% of the way done" << std::endl;
0177     centileIterator += 0.1;
0178       }
0179     if(trig == 0)nEvents -> Fill(1);
0180     else if(trig == 10) nEvents -> Fill(2);
0181     else if(trig == 30) nEvents -> Fill(3);
0182     int njets = truthPt->size();
0183     int nrecojets = pt->size();
0184     for(int tj = 0; tj < njets; tj++){
0185       //fill truth hists
0186       // to avoid statistical artifacts from double counting, the datasets can't overlap
0187       float maxTJetPt = getLeadJetPt(truthPt);
0188       h_pt_true_lead_preCut->Fill(maxTJetPt);
0189 
0190       if(!isInRange(maxTJetPt, mcWeight))
0191     {
0192       if(verbosity)std::cout << "ERRCODE0: max jet pt out of range: " << std::endl;
0193       continue;
0194     }
0195       h_pt_true_lead->Fill(maxTJetPt, mcWeight);
0196       
0197       
0198       
0199       //do reco to truth jet matching
0200       float matchEta, matchPhi, matchPt, matchE, matchsubtracted_et, dR;
0201       float dRMax = 100;
0202       int recoMatchJet = -9999;
0203 
0204       for(int rj = 0; rj < nrecojets; rj++){
0205     if(abs(eta->at(rj)) > 1.0-jetR*0.1)
0206       {
0207         if(verbosity)std::cout << "ERRCODE1: reco jet outside eta range" << std::endl;
0208         continue;
0209       }
0210 
0211     float dEta = truthEta->at(tj) - eta->at(rj);
0212     float dPhi = truthPhi->at(tj) - phi->at(rj);
0213     while(dPhi > TMath::Pi()) dPhi -= 2*TMath::Pi();
0214     while(dPhi < -TMath::Pi()) dPhi += 2*TMath::Pi();
0215     dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0216     
0217     if(dR < dRMax){
0218       matchEta = eta->at(rj);
0219       matchPhi = phi->at(rj);
0220       matchPt = pt->at(rj);
0221       matchE = e->at(rj);
0222       matchsubtracted_et = subtracted_et->at(rj);
0223       dRMax = dR;
0224       recoMatchJet = rj;
0225     }
0226       }
0227       int etaBin = getEtaBin(matchEta, nEtaBins-extraBins, etaBins);
0228       if(etaBin == -1) 
0229         {
0230       if(verbosity)std::cout << "ERRCODE2: match eta: " << matchEta << "; match eta bin: " << etaBin << std::endl;
0231       continue; //outside acceptance
0232     }
0233       if(dRMax > maxDist)
0234     {
0235       if(verbosity)std::cout << "ERRCODE3: drMax > maxDist: " << dRMax << std::endl;
0236       continue;
0237     }
0238      
0239       if(!isLin)
0240     {
0241       float corr = 1;
0242       if(applyCorr)
0243         {
0244           if(correctionExtrap[etaBin] -> Eval(matchPt) > 1)
0245         {
0246           corr = correctionExtrap[etaBin] -> Eval(matchPt);
0247         }
0248           else
0249         {
0250           corr = correction[etaBin] -> Eval(matchPt);
0251         }
0252         }
0253       h_MC_Reso_pt[etaBin]->Fill(truthPt->at(tj),corr*matchPt/truthPt->at(tj), mcWeight);//eta differential JES
0254       
0255       h_MC_Reso_pt[nEtaBins-2]->Fill(truthPt->at(tj),corr*matchPt/truthPt->at(tj), mcWeight);//inclusive JES with eta differential correction
0256 
0257       if(applyCorr)corr = correction[nEtaBins - 2] -> Eval(matchPt);
0258       h_MC_Reso_pt[nEtaBins-1]->Fill(truthPt->at(tj),corr*matchPt/truthPt->at(tj), mcWeight);//inclusive JES with eta inclusive correction
0259       
0260     }
0261       else 
0262     {
0263       //We're looking at the linearity, fill with just the matchpt
0264       h_MC_Reso_pt[etaBin]->Fill(truthPt->at(tj),matchPt, mcWeight);
0265       h_MC_Reso_pt[nEtaBins-2]->Fill(truthPt->at(tj),matchPt, mcWeight);
0266     }
0267       
0268       h_pt_true_matched[etaBin]->Fill(truthPt->at(tj));
0269       h_pt_true_matched[nEtaBins-2]->Fill(truthPt->at(tj));
0270    
0271     }
0272   }
0273   f_out -> cd();
0274   for(int i = 0; i < nEtaBins; i++) 
0275     {
0276       h_MC_Reso_pt[i] -> Write();
0277     }
0278   recoDist -> Write();
0279   truthDist -> Write();
0280   for(int i = 0; i < nEtaBins-1; i++)
0281     {
0282       h_pt_true_matched[i]->Write();
0283     }
0284   h_pt_true_lead -> Write();
0285   h_pt_true_lead_preCut -> Write();
0286   nEvents -> Write();
0287   f_out -> Close();
0288   std::cout << "All done!" << std::endl;
0289 }
0290 
0291 int getEtaBin(float jetEta, int nEtaBins, std::vector<float> etaBins)
0292 {
0293   //std::cout << "jetEta: " << jetEta << std::endl;
0294   int etaBin = -1;
0295   if(jetEta > etaBins.at(nEtaBins)) return -1;
0296   if(jetEta < etaBins.at(0)) return -1;
0297   for(int i = 0; i < nEtaBins; i++)
0298     {
0299       if(jetEta > etaBins.at(i) && jetEta <= etaBins.at(i+1))
0300     {
0301       etaBin = i;
0302       //return etaBin;
0303     }
0304     }
0305   //std::cout << "etaBin: " << etaBin << std::endl;
0306   return etaBin;
0307 }
0308 
0309 bool isInRange(float truthJetPt, float mcWeight)
0310 {
0311   float ptCutOff1 = -1; 
0312   float ptCutOff2 = -1;
0313   
0314   
0315   if(abs(mcWeight/39.06e-3 - 1) < 1e-7)//Minimum bias
0316     {
0317       //std::cout << "MB event found" << std::endl;
0318       ptCutOff1 = 0;
0319       ptCutOff2 = 14;
0320     }
0321   else if(abs(mcWeight/3.210e-6 - 1) < 1e-7)
0322     {
0323       //std::cout << "10GeV event found" << std::endl;
0324       ptCutOff1 = 14;
0325       ptCutOff2 = 37;
0326     }
0327   else if(abs(mcWeight/2.178e-9 - 1) < 1e-7)
0328     {
0329       //std::cout << "30GeV event found" << std::endl;
0330       ptCutOff1 = 37;
0331       ptCutOff2 = 3000;
0332     }
0333 
0334   if(truthJetPt < ptCutOff2 && truthJetPt >= ptCutOff1) return true;
0335   return false;
0336 }
0337  
0338 float getLeadJetPt(std::vector<float> *tJetPt)
0339 {
0340   float maxpt = 0;
0341   for(int i = 0; i < tJetPt->size(); i++)
0342     {
0343       if(tJetPt->at(i) > maxpt) maxpt = tJetPt->at(i);
0344     }
0345   
0346   return maxpt;
0347 }
[ 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!