/ElectroWeakAnalysis/ZEE/macros/PerformAnalysis.cpp
C++ | 628 lines | 507 code | 62 blank | 59 comment | 51 complexity | 3a0d0be621e118257920558a34169e39 MD5 | raw file
1#include "TROOT.h" 2#include "TFile.h" 3#include "TTree.h" 4#include "TH1.h" 5#include "TF1.h" 6#include "TRandom3.h" 7#include "TString.h" 8#include <cmath> 9#include <iostream> 10#include <fstream> 11 12using namespace std; 13 14void PerformAnalysis(TString process, double eventweight, TString datapath)//, ofstream& results) //int main() 15{ 16 cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%\%" << endl; 17 cout << "\% " << "Perform Analysis" << endl; 18 cout << "\% " << process << endl; 19 cout << "\% " << "Event Weight = " << eventweight << endl; 20 cout << "\% " << "Data Path = " << datapath << endl; 21 cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%\%" << endl; 22 23 // Declare electron cut value variables 24 double cMEt = 30.; 25 double cPt = 30.; 26 double cECALiso_EB = 4.2, cECALiso_EE = 3.4; 27 double cHCALiso_EB = 2.0, cHCALiso_EE = 1.3; 28 double cTrackiso_EB = 2.2, cTrackiso_EE = 1.1; 29 double cDeltaEta_EB = 0.0040, cDeltaEta_EE = 0.0066; 30 double cDeltaPhi_EB = 0.025, cDeltaPhi_EE = 0.020; 31 double csIhIh_EB = 0.0099, csIhIh_EE = 0.0280; 32 // Declare neutrino cut value variables 33 double cHCAL = 6.2; 34 double cHCALEt = 12; 35 double cf1x5 = 0.83, cf2x5 = 0.93; 36 int celecmatch = 0; 37 double cnusIhIh = 0.027; 38 39 cout << "Cut Values:" << endl; 40 41 cout << "MEt cut " << cMEt << endl; 42 43 cout << "Electron selection cuts:" << endl; 44 cout << "Pt cut " << cPt << endl; 45 cout << "ECAL Isolation cut (EB) " << cECALiso_EB << endl; 46 cout << "ECAL Isolation cut (EE) " << cECALiso_EE << endl; 47 cout << "HCAL Isolation cut (EB) " << cHCALiso_EB << endl; 48 cout << "HCAL Isolation cut (EE) " << cHCALiso_EE << endl; 49 cout << "Track Isolation cut (EB) " << cTrackiso_EB << endl; 50 cout << "Track Isolation cut (EE) " << cTrackiso_EE << endl; 51 cout << "Delta Eta cut (EB) " << cDeltaEta_EB << endl; 52 cout << "Delta Eta cut (EE) " << cDeltaEta_EE << endl; 53 cout << "Delta Phi cut (EB) " << cDeltaPhi_EB << endl; 54 cout << "Delta Phi cut (EE) " << cDeltaPhi_EE << endl; 55 cout << "Sigma iEta iEta cut (EB) " << csIhIh_EB << endl; 56 cout << "Sigma iEta iEta cut (EE) " << csIhIh_EE << endl; 57 58 cout << "Probe selection cuts:" << endl; 59 cout << "HCAL Energy cut " << cHCAL << endl; 60 cout << "HCAL Transverse Energy cut " << cHCALEt << endl; 61 cout << "Fraction of energy in 1x5 cut " << cf1x5 << endl; 62 cout << "Fraction of energy in 2x5 cut " << cf2x5 << endl; 63 cout << "Require electron match " << celecmatch << endl; 64 cout << "Sigma iEta iEta cut " << cnusIhIh << endl; 65 66 cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%\%" << endl; 67 // Import probe selection efficiency weights 68 double nueff[345]; 69 ifstream weightsin; 70 weightsin.open("EtaWeights.txt", ifstream::in); 71 72 for(int eta=0; eta < 345; ++eta) 73 { 74 double weight; 75 weightsin >> weight; 76 //cout << eta << "\t" << weight << endl; 77 nueff[eta] = weight; 78 } 79 weightsin.close(); 80 cout << "Imported probe selection efficiencies" << endl; 81 82 TString OutFileName = process+".root"; 83 TFile* outfile = TFile::Open(OutFileName, "recreate"); 84 85 cout << "Created output file \"" << OutFileName << "\"" << endl; 86 87 TH1F* h_dataWMEt_pass_EB = new TH1F("dataWMEt_pass_EB","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 88 TH1F* h_dataWMEt_pass_EE = new TH1F("dataWMEt_pass_EE","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 89 TH1F* h_dataWMEt_fail_EB = new TH1F("dataWMEt_fail_EB","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 90 TH1F* h_dataWMEt_fail_EE = new TH1F("dataWMEt_fail_EE","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 91 92 TH1F* h_mcWMEtin_pass_EB = new TH1F("mcWMEtin_pass_EB","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 93 TH1F* h_mcWMEtin_pass_EE = new TH1F("mcWMEtin_pass_EE","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 94 TH1F* h_mcWMEtin_fail_EB = new TH1F("mcWMEtin_fail_EB","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 95 TH1F* h_mcWMEtin_fail_EE = new TH1F("mcWMEtin_fail_EE","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 96 TH1F* h_mcWMEtout_pass_EB = new TH1F("mcWMEtout_pass_EB","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 97 TH1F* h_mcWMEtout_pass_EE = new TH1F("mcWMEtout_pass_EE","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 98 TH1F* h_mcWMEtout_fail_EB = new TH1F("mcWMEtout_fail_EB","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 99 TH1F* h_mcWMEtout_fail_EE = new TH1F("mcWMEtout_fail_EE","W#rightarrow e#nu MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 100 101 TH1F* h_ErsatzMEt_pass_EB = new TH1F("ErsatzMEt_pass_EB","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 102 TH1F* h_ErsatzMEt_pass_EE = new TH1F("ErsatzMEt_pass_EE","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 103 TH1F* h_ErsatzMEt_fail_EB = new TH1F("ErsatzMEt_fail_EB","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 104 TH1F* h_ErsatzMEt_fail_EE = new TH1F("ErsatzMEt_fail_EE","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 105 106 TH1F* h_ErsatzMEt_probept = new TH1F("ErsatzMEt_probept","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 107 TH1F* h_ErsatzMEt_uncorr = new TH1F("ErsatzMEt_uncorr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 108 TH1F* h_ErsatzMEt_fetacorr = new TH1F("ErsatzMEt_fetacorr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 109 TH1F* h_ErsatzMEt_fbremcorr = new TH1F("ErsatzMEt_fbremcorr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 110 111 TH1F* h_ErsatzMEt_pass_EB_peakfit = new TH1F("ErsatzMEt_pass_EB_peakfit","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 112 TH1F* h_ErsatzMEt_pass_EE_peakfit = new TH1F("ErsatzMEt_pass_EE_peakfit","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 113 TH1F* h_WMEt_pass_EB_peakfit = new TH1F("WMEt_pass_EB_peakfit","W MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 114 TH1F* h_WMEt_pass_EE_peakfit = new TH1F("WMEt_pass_EE_peakfit","W MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 115 116 TH1F* h_ErsatzMEt_pass_EB_shifted = new TH1F("ErsatzMEt_pass_EB_shifted","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 117 TH1F* h_ErsatzMEt_pass_EE_shifted = new TH1F("ErsatzMEt_pass_EE_shifted","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 118 TH1F* h_ErsatzMEt_fail_EB_shifted = new TH1F("ErsatzMEt_fail_EB_shifted","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 119 TH1F* h_ErsatzMEt_fail_EE_shifted = new TH1F("ErsatzMEt_fail_EE_shifted","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 120 121 TH1F* h_acceptance_correction_pass_EB = new TH1F("acceptacne_correction_pass_EB", "Acceptance Correction pass EB;#slash{E_{T}};Acceptance Correction", 100, 0., 100.); 122 TH1F* h_acceptance_correction_pass_EE = new TH1F("acceptacne_correction_pass_EE", "Acceptance Correction pass EE;#slash{E_{T}};Acceptance Correction", 100, 0., 100.); 123 TH1F* h_acceptance_correction_fail_EB = new TH1F("acceptacne_correction_fail_EB", "Acceptance Correction fail EB;#slash{E_{T}};Acceptance Correction", 100, 0., 100.); 124 TH1F* h_acceptance_correction_fail_EE = new TH1F("acceptacne_correction_fail_EE", "Acceptance Correction fail EE;#slash{E_{T}};Acceptance Correction", 100, 0., 100.); 125 126 TH1F* h_ErsatzMEt_pass_EB_corr = new TH1F("ErsatzMEt_pass_EB_corr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 127 TH1F* h_ErsatzMEt_pass_EE_corr = new TH1F("ErsatzMEt_pass_EE_corr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 128 TH1F* h_ErsatzMEt_fail_EB_corr = new TH1F("ErsatzMEt_fail_EB_corr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 129 TH1F* h_ErsatzMEt_fail_EE_corr = new TH1F("ErsatzMEt_fail_EE_corr","Ersatz MET;#slash{E}_{T} (GeV);Arbitrary Units", 100, 0., 100.); 130 131 cout << "Declared Histograms" << endl; 132 133 vector<double> ErsatzMEt_pass_EB; 134 vector<double> ErsatzMEt_pass_EE; 135 vector<double> ErsatzMEt_fail_EB; 136 vector<double> ErsatzMEt_fail_EE; 137 138 vector<double> Weight_pass_EB; 139 vector<double> Weight_pass_EE; 140 vector<double> Weight_fail_EB; 141 vector<double> Weight_fail_EE; 142 143 TRandom3 r; 144 145 TString WFileName = datapath+"WenuTrue.root"; 146 TFile *fileW = TFile::Open(WFileName); 147 cout << "Opened W Monte Carlo file" << endl; 148 TTree *t = (TTree*) fileW->Get("analyse/AnalysisData"); 149 cout << "Got W TTree" << endl; 150 151 long nEntries = t->GetEntriesFast(); 152 cout << "Total number of W events = " << nEntries << endl; 153 154 double elec_pt_W[4], elec_eta_W[4]; 155 double elec_trckIso_W[4]/*, elec_EcalIso_W[4], elec_HcalIso_W[4]*/; 156 double elec_sigIhIh_W[4]/*, elec_dPhi_W[4], elec_dEta_W[4]*/; 157 158 double nu_pt_W, nu_eta_W, nu_ECALeta_W, nu_phi_W; 159 160 //double McW_pt, McW_phi; 161 double CaloMEt_W;//, CaloMEt25, CaloMEt30; 162 //double CaloMt[4];// CaloMt25[4], CaloMt30[4]; 163 164// TBranch* bMcW_pt = t->GetBranch("Boson_pt"); 165// bMcW_pt->SetAddress(&McW_pt); 166// TBranch* bMcW_phi = t->GetBranch("Boson_phi"); 167// bMcW_phi->SetAddress(&McW_phi); 168// TBranch* bNSelElecs = t->GetBranch("nSelElecs"); 169// bNSelElecs->SetAddress(&nSelElecs); 170 TBranch* bElec_eta = t->GetBranch("elec_eta"); 171 bElec_eta->SetAddress(&elec_eta_W); 172 TBranch* bElec_pt = t->GetBranch("elec_pt"); 173 bElec_pt->SetAddress(&elec_pt_W); 174 175 TBranch* bTag_sIhIh_W = t->GetBranch("elec_sIhIh"); 176 bTag_sIhIh_W->SetAddress(&elec_sigIhIh_W); 177 //TBranch* bTag_dPhi = t->GetBranch("elec_dPhiIn"); 178 //bTag_dPhi->SetAddress(&elec_dPhi_W); 179 //TBranch* bTag_dEta = t->GetBranch("elec_dEtaIn"); 180 //bTag_dEta->SetAddress(&elec_dEta_W); 181 TBranch* bTag_tIso_W = t->GetBranch("elec_isoTrack"); 182 bTag_tIso_W->SetAddress(&elec_trckIso_W); 183 //TBranch* bTag_eIso = t->GetBranch("elec_isoEcal"); 184 //bTag_eIso->SetAddress(&elec_EcalIso_W); 185 //TBranch* bTag_hIso = t->GetBranch("elec_isoHcal"); 186 //bTag_hIso->SetAddress(&elec_HcalIso_W); 187 188 // TBranch* = t->GetBranch(""); 189// ->SetAddress(&); 190 TBranch* bMcNu_pt = t->GetBranch("McNu_pt"); 191 bMcNu_pt->SetAddress(&nu_pt_W); 192 TBranch* bMcNu_phi = t->GetBranch("McNu_phi"); 193 bMcNu_phi->SetAddress(&nu_phi_W); 194 TBranch* bMcNu_eta = t->GetBranch("McNu_eta"); 195 bMcNu_eta->SetAddress(&nu_eta_W); 196 TBranch* bMcNu_ECALeta = t->GetBranch("McNu_ECALeta"); 197 bMcNu_ECALeta->SetAddress(&nu_ECALeta_W); 198 TBranch* bCalo_MEt = t->GetBranch("caloMEt"); 199 bCalo_MEt->SetAddress(&CaloMEt_W); 200// TBranch* bCalo_MEt25 = t->GetBranch("caloMEt25"); 201// bCalo_MEt25->SetAddress(&CaloMEt25); 202// TBranch* bCalo_MEt30 = t->GetBranch("caloMEt30"); 203// bCalo_MEt30->SetAddress(&CaloMEt30); 204// TBranch* bCalo_Mt = t->GetBranch("caloMt"); 205// bCalo_Mt->SetAddress(&CaloMt); 206 cout << "Set up branches" << endl; 207 208 //long nentries = t->GetEntries(); 209 //int index = 0; 210 int aaa = 0, bbb = 0, ccc = 0, ddd = 0; 211 for(long i = 0; i < nEntries; ++i) 212 { 213 if(i%100000 == 0) cout <<"Analysing event "<< i << endl; 214 //if (i == ChosenEvents[index]) 215 //{ 216 //index++; 217 //bool iIsChosen = (i == ChosenEvents[index]); 218 t->GetEntry(i); 219 if(elec_pt_W[0] > cPt) 220 { 221 bool pass_e_cuts = false; 222 bool pass_trkiso_cut = false; 223 bool inBarrel = false; 224 bool inEndcap = false; 225 if(fabs(elec_eta_W[0]) < 1.4442) 226 { 227 pass_e_cuts = (elec_sigIhIh_W[0] < csIhIh_EB); 228 pass_trkiso_cut = (elec_trckIso_W[0] < cTrackiso_EB); 229 inBarrel = true; 230 }else if(fabs(elec_eta_W[0]) < 2.5) 231 { 232 pass_e_cuts = (elec_sigIhIh_W[0] < csIhIh_EE); 233 pass_trkiso_cut = (elec_trckIso_W[0] < cTrackiso_EE); 234 inEndcap = true; 235 } 236 if(pass_e_cuts) 237 { 238 if(pass_trkiso_cut) 239 { 240 if(inBarrel) 241 { 242 h_dataWMEt_pass_EB->Fill(CaloMEt_W); 243 h_WMEt_pass_EB_peakfit->Fill(CaloMEt_W); 244 aaa++; 245 if(fabs(nu_eta_W) < 2.5) 246 { 247 h_mcWMEtin_pass_EB->Fill(CaloMEt_W); 248 }else{ 249 h_mcWMEtout_pass_EB->Fill(CaloMEt_W); 250 } 251 } 252 if(inEndcap) 253 { 254 h_dataWMEt_pass_EE->Fill(CaloMEt_W); 255 h_WMEt_pass_EE_peakfit->Fill(CaloMEt_W); 256 bbb++; 257 if(fabs(nu_eta_W) < 2.5) 258 { 259 h_mcWMEtin_pass_EE->Fill(CaloMEt_W); 260 }else{ 261 h_mcWMEtout_pass_EE->Fill(CaloMEt_W); 262 } 263 } 264 }else 265 { 266 if(inBarrel) 267 { 268 h_dataWMEt_fail_EB->Fill(CaloMEt_W); 269 ccc++; 270 if(fabs(nu_eta_W) < 2.5) 271 { 272 h_mcWMEtin_fail_EB->Fill(CaloMEt_W); 273 }else{ 274 h_mcWMEtout_fail_EB->Fill(CaloMEt_W); 275 } 276 } 277 if(inEndcap) 278 { 279 h_dataWMEt_fail_EE->Fill(CaloMEt_W); 280 ddd++; 281 if(fabs(nu_eta_W) < 2.5) 282 { 283 h_mcWMEtin_fail_EE->Fill(CaloMEt_W); 284 }else{ 285 h_mcWMEtout_fail_EE->Fill(CaloMEt_W); 286 } 287 } 288 } 289 } 290 } 291 } 292 fileW->Close(); 293 cout << "Closed W Monte Carlo file" << endl; 294 295 cout << "Number of W events passing selection cuts = " << aaa+bbb+ccc+ddd << endl; 296 cout << "Number Pass EB = " << aaa << endl; 297 cout << "Number Pass EE = " << bbb << endl; 298 cout << "Number Fail EB = " << ccc << endl; 299 cout << "Number Fail EE = " << ddd << endl; 300 301 TString ErsatzFileName = datapath+process+".root"; 302 TFile *fileZ = TFile::Open(ErsatzFileName); 303 cout << "Opened Ersatz data file" << endl; 304 t = (TTree*) fileZ->Get("ErsatzMEt/ErsatzMEt"); 305 cout << "Got ersatz TTree" << endl; 306 nEntries = t->GetEntries(); 307 cout << "Total number of ersatz events = " << nEntries << endl; 308 309 int nErNu; 310 double tag_pt[4], tag_eta[4], tag_phi[4], probe_pt[4], probe_eta[4], probe_phi[4]; 311 double ErsatzV1MEt[4], ErsatzV1aMEt[4], ErsatzV1bMEt[4]; 312 double elec_trkIso[4], elec_ECALIso[4], elec_HCALIso[4]; 313 double elec_sigIhIh[4], elec_dPhi[4], elec_dEta[4]; 314 double tag_rescPt[4], mesc[4]; 315 double nu_e1x5[4], nu_e2x5[4], nu_e5x5[4], nu_sigIhIh[4]; 316 double nu_HCALEt[4], nu_HCAL[4], nu_trckIso[4]; 317 double caloMEt; 318 //int nu_elec[4]; 319 320 TBranch* bNum_ErNu = t->GetBranch("nProbes"); 321 bNum_ErNu->SetAddress(&nErNu); 322 //W selected electron properties 323 TBranch* bTag_eta = t->GetBranch("tag_eta"); 324 bTag_eta->SetAddress(&tag_eta); 325 TBranch* bTag_pt = t->GetBranch("tag_pt"); 326 bTag_pt->SetAddress(&tag_pt); 327 TBranch* bTag_phi = t->GetBranch("tag_phi"); 328 bTag_phi->SetAddress(&tag_phi); 329 TBranch* bTag_rescPt = t->GetBranch("tag_rescPt"); 330 bTag_rescPt->SetAddress(&tag_rescPt); 331// TBranch* = t->GetBranch(""); 332// ->SetAddress(&); 333 TBranch* bTag_sIhIh = t->GetBranch("tag_sIhIh"); 334 bTag_sIhIh->SetAddress(&elec_sigIhIh); 335 TBranch* bTag_dPhi = t->GetBranch("tag_dPhiIn"); 336 bTag_dPhi->SetAddress(&elec_dPhi); 337 TBranch* bTag_dEta = t->GetBranch("tag_dEtaIn"); 338 bTag_dEta->SetAddress(&elec_dEta); 339 TBranch* bTag_tIso = t->GetBranch("tag_isoTrack"); 340 bTag_tIso->SetAddress(&elec_trkIso); 341 TBranch* bTag_eIso = t->GetBranch("tag_isoEcal"); 342 bTag_eIso->SetAddress(&elec_ECALIso); 343 TBranch* bTag_hIso = t->GetBranch("tag_isoHcal"); 344 bTag_hIso->SetAddress(&elec_HCALIso); 345 //ersatz neutrino properties 346 TBranch* bProbe_pt = t->GetBranch("probe_pt"); 347 bProbe_pt->SetAddress(&probe_pt); 348 TBranch* bProbe_eta = t->GetBranch("probe_eta"); 349 bProbe_eta->SetAddress(&probe_eta); 350 TBranch* bProbe_phi = t->GetBranch("probe_phi"); 351 bProbe_phi->SetAddress(&probe_phi); 352 //TBranch* bProbe_elec = t->GetBranch("probe_elecMatch"); 353 //bProbe_elec->SetAddress(&nu_elec); 354 TBranch* bProbe_trckIso = t->GetBranch("probe_isoTrack"); 355 bProbe_trckIso->SetAddress(&nu_trckIso); 356 //TBranch* bProbe_ECALIso = t->GetBranch("probe_isoECAL"); 357 //bProbe_ECALIso->SetAddress(&nu_ECALIso); 358 //TBranch* bProbe_HCALIso = t->GetBranch("probe_isoHCAL"); 359 //bProbe_HCALIso->SetAddress(&nu_HCALIso); 360 TBranch* bProbe_sIhIh = t->GetBranch("probe_sIhIh"); 361 bProbe_sIhIh->SetAddress(&nu_sigIhIh); 362 //TBranch* bProbe_DeltaEta = t->GetBranch("probe_DeltaEta"); 363 //bProbe_DeltaEta->SetAddress(&nu_DeltaEta); 364 //TBranch* bProbe_DeltaPhiIso = t->GetBranch("probe_DeltaPhi"); 365 //bProbe_DeltaPhi->SetAddress(&nu_DeltaPhi); 366 TBranch* bProbe_e1x5 = t->GetBranch("probe_e1x5Max"); 367 bProbe_e1x5->SetAddress(&nu_e1x5); 368 TBranch* bProbe_e2x5 = t->GetBranch("probe_e2x5Max"); 369 bProbe_e2x5->SetAddress(&nu_e2x5); 370 TBranch* bProbe_e5x5 = t->GetBranch("probe_e5x5"); 371 bProbe_e5x5->SetAddress(&nu_e5x5); 372 TBranch* bProbe_HCAL = t->GetBranch("probe_HcalE015"); 373 bProbe_HCAL->SetAddress(&nu_HCAL); 374 TBranch* bProbe_HCALEt = t->GetBranch("probe_HcalEt015"); 375 bProbe_HCALEt->SetAddress(&nu_HCALEt); 376 //Ersatz MEt results 377 TBranch* bErsatzV1_MEt = t->GetBranch("ErsatzV1CaloMEt"); 378 bErsatzV1_MEt->SetAddress(&ErsatzV1MEt); 379 TBranch* bErsatzV1a_MEt = t->GetBranch("ErsatzV1aCaloMEt"); 380 bErsatzV1a_MEt->SetAddress(&ErsatzV1aMEt); 381 TBranch* bErsatzV1b_MEt = t->GetBranch("ErsatzV1bCaloMEt"); 382 bErsatzV1b_MEt->SetAddress(&ErsatzV1bMEt); 383 TBranch* bMesc = t->GetBranch("ErsatzV1_Mesc"); 384 bMesc->SetAddress(&mesc); 385 TBranch* bCaloMEt = t->GetBranch("recoCaloMEt"); 386 bCaloMEt->SetAddress(&caloMEt); 387 cout << "Set up Branches" << endl; 388 389 aaa=0, bbb=0, ccc=0, ddd=0; 390 for(int i=0; i < nEntries; ++i) 391 { 392 if(i%100000 == 0) cout <<"Processing event "<< i << endl; 393 t->GetEntry(i); 394 for(int j = 0; j < nErNu; ++j) 395 { 396 bool passEtCut = false; 397 if(tag_rescPt[j] > cPt) passEtCut = true; 398 /* 399 if(process == "Zee" || process == "BCtoE_30to80" || process == "BCtoE_80to170"){ 400 if(tag_rescPt[j] > cPt) passEtCut = true; 401 }else{ 402 if(tag_pt[j] > (91.188/80.398)*cPt) passEtCut = true; 403 } 404 */ 405 if(passEtCut) 406 { 407 if(fabs(mesc[j]-91.1876) < 21.1876) 408 { 409 bool pass_e_cuts = false; 410 bool pass_trkiso_cut = false; 411 bool inBarrel = false; 412 bool inEndcap = false; 413 if(fabs(tag_eta[j])<1.4442) 414 { 415 pass_e_cuts = (elec_ECALIso[j] < cECALiso_EB && elec_HCALIso[j] < cHCALiso_EB 416 && elec_sigIhIh[j] < csIhIh_EB && elec_dPhi[j] < cDeltaPhi_EB 417 && elec_dEta[j] < cDeltaEta_EB); 418 pass_trkiso_cut = (elec_trkIso[j] < cTrackiso_EB); 419 inBarrel = true; 420 }else if(fabs(tag_eta[j] < 2.5)) 421 { 422 pass_e_cuts = (elec_ECALIso[j] < cECALiso_EE && elec_HCALIso[j] < cHCALiso_EE 423 && elec_sigIhIh[j] < csIhIh_EE && elec_dPhi[j] < cDeltaPhi_EE 424 && elec_dEta[j] < cDeltaEta_EE); 425 pass_trkiso_cut = (elec_trkIso[j] < cTrackiso_EE); 426 inEndcap = true; 427 } 428 if(pass_e_cuts) 429 { 430 bool pass_nu_cuts = false; 431 double f1x5 = nu_e1x5[j]/nu_e5x5[j]; 432 double f2x5 = nu_e2x5[j]/nu_e5x5[j]; 433 if(fabs(probe_eta[j]) < 1.4442) 434 { 435 pass_nu_cuts = (nu_HCAL[j] < cHCAL && (f1x5 > cf1x5 || f2x5 > cf2x5) 436 /*&& nu_elec[j] == celecmatch*/); 437 }else if(fabs(probe_eta[j] < 2.5)){ 438 pass_nu_cuts = (nu_HCALEt[j] < cHCALEt && nu_sigIhIh[j] < cnusIhIh 439 /*&& nu_elec[j] == celecmatch*/); 440 } 441 if(pass_nu_cuts) 442 { 443 int EtaInt = int((probe_eta[j] + 3.)/0.01739); 444 double weight = eventweight/nueff[EtaInt]; 445 if(pass_trkiso_cut) 446 { 447 if(inBarrel) 448 { 449 aaa++; 450 ErsatzMEt_pass_EB.push_back(ErsatzV1bMEt[j]); 451 Weight_pass_EB.push_back(weight); 452 h_ErsatzMEt_pass_EB->Fill(ErsatzV1bMEt[j], weight); 453 h_ErsatzMEt_pass_EB_peakfit->Fill(ErsatzV1bMEt[j], weight); 454 h_ErsatzMEt_probept->Fill(probe_pt[j], weight); 455 h_ErsatzMEt_uncorr->Fill(ErsatzV1MEt[j], weight); 456 h_ErsatzMEt_fetacorr->Fill(ErsatzV1aMEt[j], weight); 457 h_ErsatzMEt_fbremcorr->Fill(ErsatzV1bMEt[j], weight); 458 } 459 if(inEndcap) 460 { 461 bbb++; 462 ErsatzMEt_pass_EE.push_back(ErsatzV1bMEt[j]); 463 Weight_pass_EE.push_back(weight); 464 h_ErsatzMEt_pass_EE->Fill(ErsatzV1bMEt[j], weight); 465 h_ErsatzMEt_pass_EE_peakfit->Fill(ErsatzV1bMEt[j], weight); 466 h_ErsatzMEt_probept->Fill(probe_pt[j], weight); 467 h_ErsatzMEt_uncorr->Fill(ErsatzV1MEt[j], weight); 468 h_ErsatzMEt_fetacorr->Fill(ErsatzV1aMEt[j], weight); 469 h_ErsatzMEt_fbremcorr->Fill(ErsatzV1bMEt[j], weight); 470 } 471 }else{ 472 if(inBarrel) 473 { 474 ccc++; 475 ErsatzMEt_fail_EB.push_back(ErsatzV1bMEt[j]); 476 Weight_fail_EB.push_back(weight); 477 h_ErsatzMEt_fail_EB->Fill(ErsatzV1bMEt[j], weight); 478 } 479 if(inEndcap) 480 { 481 ddd++; 482 ErsatzMEt_fail_EE.push_back(ErsatzV1bMEt[j]); 483 Weight_fail_EE.push_back(weight); 484 h_ErsatzMEt_fail_EE->Fill(ErsatzV1bMEt[j], weight); 485 } 486 } 487 } 488 } 489 } 490 } 491 } 492 } 493 fileZ->Close(); 494 cout << "Closed Ersatz data file" << endl; 495 496 cout << "Number of events passing selection cuts = " << aaa+bbb+ccc+ddd << endl; 497 cout << "Number Pass EB = " << aaa << endl; 498 cout << "Number Pass EE = " << bbb << endl; 499 cout << "Number Fail EB = " << ccc << endl; 500 cout << "Number Fail EE = " << ddd << endl; 501 502 cout << "Apply shift correction ..." << endl; 503 504 int maxbin; 505 506 h_WMEt_pass_EB_peakfit->Scale(1./h_WMEt_pass_EB_peakfit->Integral(0,100)); 507 maxbin = h_WMEt_pass_EB_peakfit->GetMaximumBin(); 508 TF1 peakW_EB = TF1("peakW_EB", "gaus", maxbin-4, maxbin+4); 509 h_WMEt_pass_EB_peakfit->Fit("peakW_EB", "MR"); 510 cout << "W MEt maximum bin = " << maxbin << "\tPeak of Gaussian = " << peakW_EB.GetParameter(1) << endl; 511 h_WMEt_pass_EB_peakfit->Draw(); 512 513 h_ErsatzMEt_pass_EB_peakfit->Scale(1./h_ErsatzMEt_pass_EB_peakfit->Integral(0,100)); 514 maxbin = h_ErsatzMEt_pass_EB_peakfit->GetMaximumBin(); 515 TF1 peakZ_EB = TF1("peakZ_EB", "gaus", maxbin-4, maxbin+4); 516 h_ErsatzMEt_pass_EB_peakfit->Fit("peakZ_EB", "MR"); 517 cout << "Ersatz maximum bin = " << maxbin << "\tPeak of Gaussian = " << peakZ_EB.GetParameter(1) << endl; 518 h_ErsatzMEt_pass_EB_peakfit->Draw(); 519 520 double shift_EB = peakW_EB.GetParameter(1) - peakZ_EB.GetParameter(1); 521 cout << "EB Shift = " << shift_EB << endl; 522 523 h_WMEt_pass_EE_peakfit->Scale(1./h_WMEt_pass_EE_peakfit->Integral(0,100)); 524 maxbin = h_WMEt_pass_EE_peakfit->GetMaximumBin(); 525 TF1 peakW_EE = TF1("peakW_EE", "gaus", maxbin-4, maxbin+4); 526 h_WMEt_pass_EE_peakfit->Fit("peakW_EE", "MR"); 527 cout << "W MEt maximum bin = " << maxbin << "\tPeak of Gaussian = " << peakW_EE.GetParameter(1) << endl; 528 h_WMEt_pass_EE_peakfit->Draw(); 529 530 h_ErsatzMEt_pass_EE_peakfit->Scale(1./h_ErsatzMEt_pass_EE_peakfit->Integral(0,100)); 531 maxbin = h_ErsatzMEt_pass_EE_peakfit->GetMaximumBin(); 532 TF1 peakZ_EE = TF1("peakZ_EE", "gaus", maxbin-4, maxbin+4); 533 h_ErsatzMEt_pass_EE_peakfit->Fit("peakZ_EE", "MR"); 534 cout << "Ersatz maximum bin = " << maxbin << "\tPeak of Gaussian = " << peakZ_EE.GetParameter(1) << endl; 535 h_ErsatzMEt_pass_EE_peakfit->Draw(); 536 537 double shift_EE = peakW_EE.GetParameter(1) - peakZ_EE.GetParameter(1); 538 cout << "EE Shift = " << shift_EE << endl; 539 540 for(unsigned int i=0; i < ErsatzMEt_pass_EB.size(); i++) 541 { 542 ErsatzMEt_pass_EB[i] += shift_EB; 543 h_ErsatzMEt_pass_EB_shifted->Fill(ErsatzMEt_pass_EB[i], Weight_pass_EB[i]); 544 } 545 for(unsigned int i=0; i < ErsatzMEt_pass_EE.size(); i++) 546 { 547 ErsatzMEt_pass_EE[i] += shift_EE; 548 h_ErsatzMEt_pass_EE_shifted->Fill(ErsatzMEt_pass_EE[i], Weight_pass_EE[i]); 549 } 550 for(unsigned int i=0; i < ErsatzMEt_fail_EB.size(); i++) 551 { 552 ErsatzMEt_fail_EB[i] += shift_EB; 553 h_ErsatzMEt_fail_EB_shifted->Fill(ErsatzMEt_fail_EB[i], Weight_fail_EB[i]); 554 } 555 for(unsigned int i=0; i < ErsatzMEt_fail_EE.size(); i++) 556 { 557 ErsatzMEt_fail_EE[i] += shift_EE; 558 h_ErsatzMEt_fail_EE_shifted->Fill(ErsatzMEt_fail_EE[i], Weight_fail_EE[i]); 559 } 560 561 cout << "Apply acceptance correction ..." << endl; 562 563 TH1F* h_ones = new TH1F("ones", "Histogram of Ones;;", 100, 0., 100.); 564 for (int i=0; i<100; i++) 565 { 566 h_ones->Fill(i+0.5); 567 } 568 569 h_acceptance_correction_pass_EB->Divide(h_mcWMEtout_pass_EB, h_mcWMEtin_pass_EB); 570 h_acceptance_correction_pass_EB->Add(h_ones); 571 h_ErsatzMEt_pass_EB_corr->Multiply(h_ErsatzMEt_pass_EB_shifted, h_acceptance_correction_pass_EB); 572 573 h_acceptance_correction_pass_EE->Divide(h_mcWMEtout_pass_EE, h_mcWMEtin_pass_EE); 574 h_acceptance_correction_pass_EE->Add(h_ones); 575 h_ErsatzMEt_pass_EE_corr->Multiply(h_ErsatzMEt_pass_EE_shifted, h_acceptance_correction_pass_EE); 576 577 h_acceptance_correction_fail_EB->Divide(h_mcWMEtout_fail_EB, h_mcWMEtin_fail_EB); 578 h_acceptance_correction_fail_EB->Add(h_ones); 579 h_ErsatzMEt_fail_EB_corr->Multiply(h_ErsatzMEt_fail_EB_shifted, h_acceptance_correction_fail_EB); 580 581 h_acceptance_correction_fail_EE->Divide(h_mcWMEtout_fail_EE, h_mcWMEtin_fail_EE); 582 h_acceptance_correction_fail_EE->Add(h_ones); 583 h_ErsatzMEt_fail_EE_corr->Multiply(h_ErsatzMEt_fail_EE_shifted, h_acceptance_correction_fail_EE); 584 585 cout << "Calculating f ..." << endl; 586 587 double N_pass_EB = h_ErsatzMEt_pass_EB_corr->Integral(0,100); 588 double A_EB = h_ErsatzMEt_pass_EB_corr->Integral(int(cMEt)+1,100); 589 double B_EB = h_ErsatzMEt_pass_EB_corr->Integral(0,int(cMEt)); 590 double N_pass_EE = h_ErsatzMEt_pass_EE_corr->Integral(0,100); 591 double A_EE = h_ErsatzMEt_pass_EE_corr->Integral(int(cMEt)+1,100); 592 double B_EE = h_ErsatzMEt_pass_EE_corr->Integral(0,int(cMEt)); 593 double N_fail_EB = h_ErsatzMEt_fail_EB_corr->Integral(0,100); 594 double D_EB = h_ErsatzMEt_fail_EB_corr->Integral(int(cMEt)+1,100); 595 double C_EB = h_ErsatzMEt_fail_EB_corr->Integral(0,int(cMEt)); 596 double N_fail_EE = h_ErsatzMEt_fail_EE_corr->Integral(0,100); 597 double D_EE = h_ErsatzMEt_fail_EE_corr->Integral(int(cMEt)+1,100); 598 double C_EE = h_ErsatzMEt_fail_EE_corr->Integral(0,int(cMEt)); 599 600 double A = A_EB + A_EE; 601 double B = B_EB + B_EE; 602 double C = C_EB + C_EE; 603 double D = D_EB + D_EE; 604 double N_pass = N_pass_EB + N_pass_EE; 605 double N_fail = N_fail_EB + N_fail_EE; 606 //int N = N_pass + N_fail; 607 608 double eff = 1.0*A/(A+B); 609 double efferror = sqrt(eff*(1.-eff)/N_pass); 610 double f = 1.0*A/B; 611 double ferror = efferror/((1.-eff)*(1.-eff)); 612 613 double effprime = 1.0*D/(D+C); 614 double effprimeerror = sqrt(effprime*(1.-effprime)/N_fail); 615 double fprime = 1.0*D/C; 616 double fprimeerror = effprimeerror/((1.-effprime)*(1.-effprime)); 617 618 cout << "f\tferror\teff\tefferror\tA\tB\tN_pass" << endl; 619 cout << f << "\t" << ferror << "\t" << eff << "\t" << efferror << "\t" << A << "\t" << B << "\t" << N_pass << "\n" << endl; 620 621 cout << "fprime\tfprimeerror\teffprime\teffprimeerror\tD\tC\tN_fail" << endl; 622 cout << fprime << "\t" << fprimeerror << "\t" << effprime << "\t" << effprimeerror << "\t" << D << "\t" << C << "\t" << N_fail << "\n" << endl; 623 624// results << process << "\t" << f << "\t" << ferror << "\t" << A << "\t" << B << "\t" << fprime << "\t" << fprimeerror << "\t" << D << "\t" << C << endl; 625 626 outfile->Write(); 627 outfile->Close(); 628}