JetSpectrum/offline/analysis_data.C

0001 #include <TChain.h>
0002 #include <TFile.h>
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TF1.h>
0006 #include <TMath.h>
0007 #include <iostream>
0008 #include <vector>
0009 #include <string>
0010 #include <cmath>
0011 #include "unfold_Def.h"
0012
0013 bool check_bad_trigger(int gl1_trigger_vector_scaled[n_hcal_phibin]);
0014 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et);
0015 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi);
0016 void get_jet_filter(std::vector<bool>& jet_filter, std::vector<float>* jet_e, std::vector<float>* jet_pt, std::vector<float>* jet_eta, float zvertex, float jet_radius);
0017 void get_calibjet(std::vector<float>& calibjet_pt, std::vector<float>& calibjet_eta, std::vector<float>& calibjet_phi, std::vector<bool>& calibjet_background_dijet, std::vector<bool>& calibjet_background_frac, std::vector<bool>& jet_filter, std::vector<float>* jet_pt, std::vector<float>* jet_eta, std::vector<float>* jet_phi, std::vector<bool>* jet_background_dijet, std::vector<bool>* jet_background_frac, TF1* f_corr);
0018
0019 ////////////////////////////////////////// Main Function //////////////////////////////////////////
0020 void analysis_data(int runnumber, int nseg, int iseg)  {
0021   /////////////// General Set up ///////////////
0022   TFile *f_out = new TFile(Form("output_data/output_%d_%d_%d.root", runnumber, iseg, iseg+nseg),"RECREATE");
0023
0024   /////////////// Read Files ///////////////
0025   const char* baseDirJet = "/sphenix/tg/tg01/jets/hanpuj/CaloDataAna_skimmed";
0026   TChain chain("ttree");
0027   for (int i = iseg; i < iseg + nseg; ++i) {
0028     chain.Add(Form("%s/Run%d/OutDir%d/output_data.root", baseDirJet, runnumber, i));
0029   }
0030   chain.SetBranchStatus("*", 0);
0031
0032   float zvertex; chain.SetBranchStatus("z_vertex", 1); chain.SetBranchAddress("z_vertex", &zvertex);
0033   int gl1_trigger_vector_scaled[64]; chain.SetBranchStatus("gl1_trigger_vector_scaled", 1); chain.SetBranchAddress("gl1_trigger_vector_scaled", &gl1_trigger_vector_scaled);
0034   std::vector<float>* unsubjet_e = nullptr; chain.SetBranchStatus("unsubjet04_e", 1); chain.SetBranchAddress("unsubjet04_e", &unsubjet_e);
0035   std::vector<float>* unsubjet_pt = nullptr; chain.SetBranchStatus("unsubjet04_pt", 1); chain.SetBranchAddress("unsubjet04_pt", &unsubjet_pt);
0036   std::vector<float>* unsubjet_et = nullptr; chain.SetBranchStatus("unsubjet04_et", 1); chain.SetBranchAddress("unsubjet04_et", &unsubjet_et);
0037   std::vector<float>* unsubjet_eta = nullptr; chain.SetBranchStatus("unsubjet04_eta", 1); chain.SetBranchAddress("unsubjet04_eta", &unsubjet_eta);
0038   std::vector<float>* unsubjet_phi = nullptr; chain.SetBranchStatus("unsubjet04_phi", 1); chain.SetBranchAddress("unsubjet04_phi", &unsubjet_phi);
0039   std::vector<float>* unsubjet_emcal_calo_e = nullptr; chain.SetBranchStatus("unsubjet04_emcal_calo_e", 1); chain.SetBranchAddress("unsubjet04_emcal_calo_e", &unsubjet_emcal_calo_e);
0040   std::vector<float>* unsubjet_ihcal_calo_e = nullptr; chain.SetBranchStatus("unsubjet04_ihcal_calo_e", 1); chain.SetBranchAddress("unsubjet04_ihcal_calo_e", &unsubjet_ihcal_calo_e);
0041   std::vector<float>* unsubjet_ohcal_calo_e = nullptr; chain.SetBranchStatus("unsubjet04_ohcal_calo_e", 1); chain.SetBranchAddress("unsubjet04_ohcal_calo_e", &unsubjet_ohcal_calo_e);
0042
0043   /////////////// JES func ///////////////
0044   TFile *corrFile = new TFile("/sphenix/user/hanpuj/JES_MC_Calibration/offline/JES_Calib_Default.root", "READ");
0045   if (!corrFile) {
0046     std::cout << "Error: cannot open JES_Calib_Default.root" << std::endl;
0047     return;
0048   }
0049   TF1 *f_corr = (TF1*)corrFile->Get("JES_Calib_Default_Func");
0050   if (!f_corr) {
0051     std::cout << "Error: cannot open f_corr" << std::endl;
0052     return;
0053   }
0054
0055   /////////////// Trigger Efficiency ///////////////
0056   TFile *f_trigger = new TFile("output_jettrigeff.root", "READ");
0057   if (!f_trigger) {
0058     std::cout << "Error: cannot open output_jettrigeff.root" << std::endl;
0059     return;
0060   }
0061   TF1* f_turnon = (TF1*)f_trigger->Get("f_fit");
0062   TF1* f_turnon_down = (TF1*)f_trigger->Get("f_fit_shiftdown"); // for trigger efficiency uncertainty
0063   TF1* f_turnon_up = (TF1*)f_trigger->Get("f_fit_shiftup"); // for trigger efficiency uncertainty
0064   if (!f_turnon || !f_turnon_down || !f_turnon_up) {
0065     std::cout << "Error: cannot open trigger efficiency functions" << std::endl;
0066     return;
0067   }
0068
0069   /////////////// Histograms ///////////////
0070   TH1D *h_zvertex = new TH1D("h_zvertex", ";Z-vertex [cm]", 60, -30, 30);
0071   TH1D *h_recojet_pt_record = new TH1D("h_recojet_pt_record", ";p_{T} [GeV]", 1000, 0, 100);
0072   TH1D *h_recojet_pt_record_dijet = new TH1D("h_recojet_pt_record_dijet", ";p_{T} [GeV]", 1000, 0, 100);
0073   TH1D *h_recojet_pt_record_frac = new TH1D("h_recojet_pt_record_frac", ";p_{T} [GeV]", 1000, 0, 100);
0074   TH1D *h_calibjet_pt = new TH1D("h_calibjet_pt", ";p_{T} [GeV]", calibnpt, calibptbins); // raw calib jet pT spectrum
0075   TH1D *h_calibjet_pt_record = new TH1D("h_calibjet_pt_record", ";p_{T} [GeV]", 1000, 0, 100);
0076   TH1D *h_calibjet_pt_dijet_eff = new TH1D("h_calibjet_pt_dijet_eff", ";p_{T} [GeV]", calibnpt, calibptbins); // with trigger efficiency correction applied
0077   TH1D *h_calibjet_pt_dijet_effdown = new TH1D("h_calibjet_pt_dijet_effdown", ";p_{T} [GeV]", calibnpt, calibptbins); // for trigger efficiency uncertainty
0078   TH1D *h_calibjet_pt_dijet_effup = new TH1D("h_calibjet_pt_dijet_effup", ";p_{T} [GeV]", calibnpt, calibptbins); // for trigger efficiency uncertainty
0079   TH1D *h_calibjet_pt_frac_eff = new TH1D("h_calibjet_pt_frac_eff", ";p_{T} [GeV]", calibnpt, calibptbins); // with trigger efficiency correction applied
0080   TH1D *h_calibjet_pt_frac_effdown = new TH1D("h_calibjet_pt_frac_effdown", ";p_{T} [GeV]", calibnpt, calibptbins); // for trigger efficiency uncertainty
0081   TH1D *h_calibjet_pt_frac_effup = new TH1D("h_calibjet_pt_frac_effup", ";p_{T} [GeV]", calibnpt, calibptbins); // for trigger efficiency uncertainty
0082
0083   /////////////// Event Loop ///////////////
0084   std::cout << "Data analysis started." << std::endl;
0085   int n_events = chain.GetEntries();
0086   std::cout << "Total number of events: " << n_events << std::endl;
0087   // Event variables setup.
0088   std::vector<bool>* unsubjet_background_dijet = new std::vector<bool>;
0089   std::vector<bool>* unsubjet_background_frac = new std::vector<bool>;
0090   std::vector<bool> jet_filter;
0091   std::vector<float> calibjet_pt, calibjet_eta, calibjet_phi;
0092   std::vector<bool> calibjet_background_dijet, calibjet_background_frac;
0093   for (int ie = 0; ie < n_events; ++ie) { // event loop start
0094     // Load event.
0095     if (ie % 1000 == 0) {
0096       std::cout << "Processing event " << ie << "..." << std::endl;
0097     }
0098     chain.GetEntry(ie);
0099
0100     // Trigger and Z-vertex cut.
0101     if (check_bad_trigger(gl1_trigger_vector_scaled)) continue;
0102     if (fabs(zvertex) > 30) continue;
0103
0104     // Get leading jet and subleading jet. Do basic jet cuts.
0105     int leadingunsubjet_index = -1;
0106     int subleadingunsubjet_index = -1;
0107     get_leading_subleading_jet(leadingunsubjet_index, subleadingunsubjet_index, unsubjet_pt);
0108     if (leadingunsubjet_index < 0) continue;
0109
0110     // Fill z-vertex histogram.
0111     h_zvertex->Fill(zvertex);
0112
0113     // Beam background cut. Jet level recording.
0114     unsubjet_background_dijet->clear();
0115     unsubjet_background_frac->clear();
0116     // Dijet cut.
0117     bool match_dijet = false;
0118     if (subleadingunsubjet_index >= 0 && unsubjet_e->at(subleadingunsubjet_index) / (float) unsubjet_e->at(leadingunsubjet_index) > 0.3) {
0119       match_dijet = match_leading_subleading_jet(unsubjet_phi->at(leadingunsubjet_index), unsubjet_phi->at(subleadingunsubjet_index));
0120     }
0121     bool isbeambackground_dijet = true;
0122     if (match_dijet) isbeambackground_dijet = false;
0123     for (int ij = 0; ij < unsubjet_e->size(); ++ij) {
0124       if (isbeambackground_dijet) {
0125         unsubjet_background_dijet->push_back(true);
0126       } else {
0127         unsubjet_background_dijet->push_back(false);
0128       }
0129     }
0130     // Fraction cut.
0131     for (int ij = 0; ij < unsubjet_e->size(); ++ij) {
0132       double emfrac = unsubjet_emcal_calo_e->at(ij) / (double)(unsubjet_e->at(ij));
0133       double ihfrac = unsubjet_ihcal_calo_e->at(ij) / (double)(unsubjet_e->at(ij));
0134       double ohfrac = unsubjet_ohcal_calo_e->at(ij) / (double)(unsubjet_e->at(ij));
0135       if (emfrac < 0.1 || emfrac > 0.9 || ohfrac < 0.1 || ohfrac > 0.9 || ihfrac > 0.9) {
0136         unsubjet_background_frac->push_back(true);
0137       } else {
0138         unsubjet_background_frac->push_back(false);
0139       }
0140     }
0141
0142     // Trigger efficiency.
0143     double leadingunsubjet_et = unsubjet_et->at(leadingunsubjet_index);
0144
0145     double jettrigeff = f_turnon->Eval(leadingunsubjet_et);
0146     if (jettrigeff < 0.01) jettrigeff = 0.01;
0147     double jettrig_scale = 1.0 / jettrigeff;
0148
0149     double jettrigeff_down = f_turnon_down->Eval(leadingunsubjet_et);
0150     if (jettrigeff_down < 0.01) jettrigeff_down = 0.01;
0151     double jettrig_scale_down = 1.0 / jettrigeff_down;
0152
0153     double jettrigeff_up = f_turnon_up->Eval(leadingunsubjet_et);
0154     if (jettrigeff_up < 0.01) jettrigeff_up = 0.01;
0155     double jettrig_scale_up = 1.0 / jettrigeff_up;
0156
0157     // Fill histograms.
0158     get_jet_filter(jet_filter, unsubjet_e, unsubjet_pt, unsubjet_eta, zvertex, jet_radius);
0159     for (int ij = 0; ij < jet_filter.size(); ++ij) {
0160       if (jet_filter.at(ij)) continue;
0161       h_recojet_pt_record->Fill(unsubjet_pt->at(ij));
0162       if (!unsubjet_background_dijet->at(ij)) h_recojet_pt_record_dijet->Fill(unsubjet_pt->at(ij));
0163       if (!unsubjet_background_frac->at(ij)) h_recojet_pt_record_frac->Fill(unsubjet_pt->at(ij));
0164     }
0165
0166     get_calibjet(calibjet_pt, calibjet_eta, calibjet_phi, calibjet_background_dijet, calibjet_background_frac, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_dijet, unsubjet_background_frac, f_corr);
0167     for (int ij = 0; ij < calibjet_pt.size(); ++ij) {
0168       h_calibjet_pt->Fill(calibjet_pt[ij]);
0169       h_calibjet_pt_record->Fill(calibjet_pt[ij]);
0170       if (!calibjet_background_dijet[ij]) {
0171         h_calibjet_pt_dijet_eff->Fill(calibjet_pt[ij], jettrig_scale);
0172         h_calibjet_pt_dijet_effdown->Fill(calibjet_pt[ij], jettrig_scale_down);
0173         h_calibjet_pt_dijet_effup->Fill(calibjet_pt[ij], jettrig_scale_up);
0174       }
0175       if (!calibjet_background_frac[ij]) {
0176         h_calibjet_pt_frac_eff->Fill(calibjet_pt[ij], jettrig_scale);
0177         h_calibjet_pt_frac_effdown->Fill(calibjet_pt[ij], jettrig_scale_down);
0178         h_calibjet_pt_frac_effup->Fill(calibjet_pt[ij], jettrig_scale_up);
0179       }
0180     }
0181   } // event loop end
0182
0183   // Write histograms.
0184   std::cout << "Writing histograms..." << std::endl;
0185   f_out->cd();
0186   h_zvertex->Write();
0187   h_recojet_pt_record->Write();
0188   h_recojet_pt_record_dijet->Write();
0189   h_recojet_pt_record_frac->Write();
0190   h_calibjet_pt->Write();
0191   h_calibjet_pt_record->Write();
0192   h_calibjet_pt_dijet_eff->Write();
0193   h_calibjet_pt_dijet_effdown->Write();
0194   h_calibjet_pt_dijet_effup->Write();
0195   h_calibjet_pt_frac_eff->Write();
0196   h_calibjet_pt_frac_effdown->Write();
0197   h_calibjet_pt_frac_effup->Write();
0198   f_out->Close();
0199   std::cout << "All done!" << std::endl;
0200 }
0201
0202 ////////////////////////////////////////// Helper Functions //////////////////////////////////////////
0203 float get_deta(float eta1, float eta2) {
0204   return eta1 - eta2;
0205 }
0206
0207 float get_dphi(float phi1, float phi2) {
0208   float dphi1 = phi1 - phi2;
0209   float dphi2 = phi1 - phi2 + 2*TMath::Pi();
0210   float dphi3 = phi1 - phi2 - 2*TMath::Pi();
0211   if (fabs(dphi1) > fabs(dphi2)) {
0212     dphi1 = dphi2;
0213   }
0214   if (fabs(dphi1) > fabs(dphi3)) {
0215     dphi1 = dphi3;
0216   }
0217   return dphi1;
0218 }
0219
0220 float get_dR(float eta1, float phi1, float eta2, float phi2) {
0221   float deta = get_deta(eta1, eta2);
0222   float dphi = get_dphi(phi1, phi2);
0223   return sqrt(deta*deta + dphi*dphi);
0224 }
0225
0226 int get_hcal_up(int iphi) {
0227   if (iphi == 63) {
0228     return 0;
0229   } else {
0230     return iphi + 1;
0231   }
0232 }
0233
0234 int get_hcal_bottom(int iphi) {
0235   if (iphi == 0) {
0236     return 63;
0237   } else {
0238     return iphi - 1;
0239   }
0240 }
0241
0242 float get_emcal_mineta_zcorrected(float zvertex) {
0243   float minz_EM = -130.23;
0244   float radius_EM = 93.5;
0245   float z = minz_EM - zvertex;
0246   float eta_zcorrected = asinh(z / (float)radius_EM);
0247   return eta_zcorrected;
0248 }
0249
0250 float get_emcal_maxeta_zcorrected(float zvertex) {
0251   float maxz_EM = 130.23;
0252   float radius_EM = 93.5;
0253   float z = maxz_EM - zvertex;
0254   float eta_zcorrected = asinh(z / (float)radius_EM);
0255   return eta_zcorrected;
0256 }
0257
0258 float get_ihcal_mineta_zcorrected(float zvertex) {
0259   float minz_IH = -170.299;
0260   float radius_IH = 127.503;
0261   float z = minz_IH - zvertex;
0262   float eta_zcorrected = asinh(z / (float)radius_IH);
0263   return eta_zcorrected;
0264 }
0265
0266 float get_ihcal_maxeta_zcorrected(float zvertex) {
0267   float maxz_IH = 170.299;
0268   float radius_IH = 127.503;
0269   float z = maxz_IH - zvertex;
0270   float eta_zcorrected = asinh(z / (float)radius_IH);
0271   return eta_zcorrected;
0272 }
0273
0274 float get_ohcal_mineta_zcorrected(float zvertex) {
0275   float minz_OH = -301.683;
0276   float radius_OH = 225.87;
0277   float z = minz_OH - zvertex;
0278   float eta_zcorrected = asinh(z / (float)radius_OH);
0279   return eta_zcorrected;
0280 }
0281
0282 float get_ohcal_maxeta_zcorrected(float zvertex) {
0283   float maxz_OH = 301.683;
0284   float radius_OH = 225.87;
0285   float z = maxz_OH - zvertex;
0286   float eta_zcorrected = asinh(z / (float)radius_OH);
0287   return eta_zcorrected;
0288 }
0289
0290 bool check_bad_jet_eta(float jet_eta, float zertex, float jet_radius) {
0291   float emcal_mineta = get_emcal_mineta_zcorrected(zertex);
0292   float emcal_maxeta = get_emcal_maxeta_zcorrected(zertex);
0293   float ihcal_mineta = get_ihcal_mineta_zcorrected(zertex);
0294   float ihcal_maxeta = get_ihcal_maxeta_zcorrected(zertex);
0295   float ohcal_mineta = get_ohcal_mineta_zcorrected(zertex);
0296   float ohcal_maxeta = get_ohcal_maxeta_zcorrected(zertex);
0297   float minlimit = emcal_mineta;
0298   if (ihcal_mineta > minlimit) minlimit = ihcal_mineta;
0299   if (ohcal_mineta > minlimit) minlimit = ohcal_mineta;
0300   float maxlimit = emcal_maxeta;
0301   if (ihcal_maxeta < maxlimit) maxlimit = ihcal_maxeta;
0302   if (ohcal_maxeta < maxlimit) maxlimit = ohcal_maxeta;
0303   minlimit += jet_radius;
0304   maxlimit -= jet_radius;
0305   return jet_eta < minlimit || jet_eta > maxlimit;
0306 }
0307
0308 ////////////////////////////////////////// Functions //////////////////////////////////////////
0309 bool check_bad_trigger(int gl1_trigger_vector_scaled[64]) {
0310   return gl1_trigger_vector_scaled[18] != 1;
0311 }
0312
0313 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et) {
0314   leadingjet_index = -1;
0315   subleadingjet_index = -1;
0316   float leadingjet_et = -9999;
0317   float subleadingjet_et = -9999;
0318   for (int ij = 0; ij < jet_et->size(); ++ij) {
0319     float jetet = jet_et->at(ij);
0320     if (jetet > leadingjet_et) {
0321       subleadingjet_et = leadingjet_et;
0322       subleadingjet_index = leadingjet_index;
0323       leadingjet_et = jetet;
0324       leadingjet_index = ij;
0325     } else if (jetet > subleadingjet_et) {
0326       subleadingjet_et = jetet;
0327       subleadingjet_index = ij;
0328     }
0329   }
0330 }
0331
0332 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi) {
0333   float dijet_min_phi = 3*TMath::Pi()/4.;
0334   float dphi = get_dphi(leadingjet_phi, subleadingjet_phi);
0335   return dphi > dijet_min_phi;
0336 }
0337
0338 void get_jet_filter(std::vector<bool>& jet_filter, std::vector<float>* jet_e, std::vector<float>* jet_pt, std::vector<float>* jet_eta, float zvertex, float jet_radius) {
0339   jet_filter.clear();
0340   int njet = jet_e->size();
0341   for (int ij = 0; ij < njet; ++ij) {
0342     jet_filter.push_back(jet_e->at(ij) < 0 || check_bad_jet_eta(jet_eta->at(ij), zvertex, jet_radius) || jet_pt->at(ij) < 1);
0343   }
0344 }
0345
0346 void get_calibjet(std::vector<float>& calibjet_pt, std::vector<float>& calibjet_eta, std::vector<float>& calibjet_phi, std::vector<bool>& calibjet_background_dijet, std::vector<bool>& calibjet_background_frac, std::vector<bool>& jet_filter, std::vector<float>* jet_pt, std::vector<float>* jet_eta, std::vector<float>* jet_phi, std::vector<bool>* jet_background_dijet, std::vector<bool>* jet_background_frac, TF1* f_corr) {
0347   calibjet_pt.clear();
0348   calibjet_eta.clear();
0349   calibjet_phi.clear();
0350   calibjet_background_dijet.clear();
0351   calibjet_background_frac.clear();
0352   for (int ij = 0; ij < jet_filter.size(); ++ij) {
0353     if (jet_filter.at(ij)) continue;
0354     double calib_pt = f_corr->Eval(jet_pt->at(ij));
0355     if (calib_pt < calibptbins[0] || calib_pt > calibptbins[calibnpt]) continue;
0356     calibjet_pt.push_back(calib_pt);
0357     calibjet_eta.push_back(jet_eta->at(ij));
0358     calibjet_phi.push_back(jet_phi->at(ij));
0359     calibjet_background_dijet.push_back(jet_background_dijet->at(ij));
0360     calibjet_background_frac.push_back(jet_background_frac->at(ij));
0361   }
0362 }