| 29 |
|
#include <TLorentzVector.h> // 4-vector class |
| 30 |
|
#include <TVector3.h> // 3D vector class |
| 31 |
|
|
| 32 |
+ |
|
| 33 |
|
// |
| 34 |
|
// ntuple format headers |
| 35 |
|
// |
| 40 |
|
#include "TJet.hh" |
| 41 |
|
#include "SiMVAElectronSelection.h" |
| 42 |
|
#include "RunLumiRangeMap.h" |
| 43 |
+ |
#include "EfficiencyWeightsInterface.h" |
| 44 |
|
|
| 45 |
|
// |
| 46 |
|
// utility headers |
| 56 |
|
//#define BDTFAIL_THEIR_EVENTS |
| 57 |
|
//#define THEIR_EVENTS |
| 58 |
|
|
| 59 |
+ |
Double_t etaFromY(Double_t pt, Double_t y, Double_t phi, Double_t m) |
| 60 |
+ |
{ |
| 61 |
+ |
Double_t a = (1+exp(2*y))/(exp(2*y)-1); // intermediate term |
| 62 |
+ |
if(a*a<1) { cout << "a too small" << endl; assert(0); } |
| 63 |
+ |
Double_t E = sqrt( a*a*(pt*pt+m*m)/(a*a-1) ); |
| 64 |
+ |
Double_t pz = E*E - pt*pt - m*m; |
| 65 |
+ |
if(pz<0) { cout << "imag. pz" << endl; assert(0); } |
| 66 |
+ |
pz = sqrt(pz); |
| 67 |
+ |
if(y<0) pz *= -1; |
| 68 |
+ |
TLorentzVector v; |
| 69 |
+ |
v.SetPxPyPzE(pt*cos(phi),pt*sin(phi),pz,E); |
| 70 |
+ |
Double_t th = v.Theta(); |
| 71 |
+ |
return -log(tan(th/2)); |
| 72 |
+ |
}; |
| 73 |
+ |
|
| 74 |
+ |
|
| 75 |
+ |
int getChannel(mithep::TGenInfo * ginfo) { |
| 76 |
+ |
|
| 77 |
+ |
int gchannel=-1; |
| 78 |
+ |
|
| 79 |
+ |
if( abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kElectron ) gchannel=0; |
| 80 |
+ |
else if( abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kMuon ) gchannel=1; |
| 81 |
+ |
else if( (abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kMuon) || |
| 82 |
+ |
(abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kElectron) ) gchannel=2; |
| 83 |
+ |
|
| 84 |
+ |
return gchannel; |
| 85 |
+ |
} |
| 86 |
+ |
|
| 87 |
+ |
|
| 88 |
+ |
|
| 89 |
|
//=== MAIN ================================================================================================= |
| 90 |
|
int main(int argc, char** argv) |
| 91 |
|
{ |
| 138 |
|
TFile *file = new TFile(ofname, "RECREATE"); |
| 139 |
|
TH1F * h_evt = new TH1F("hevt", "hevt", 2, 0, 2 ); |
| 140 |
|
TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 1024, 0, 1024 ); |
| 109 |
– |
// TNtuple * passtuple = new TNtuple( "passtuple", "passtuple", "run:evt:lumi:channel:mZ1:mZ2:m4l:pt4l:w" ); |
| 141 |
|
TTree * passtuple = new TTree("passtuple", "passtuple" ); |
| 142 |
< |
unsigned run, evt, lumi, channel; |
| 142 |
> |
unsigned run, evt, lumi, channel, gchannel; |
| 143 |
|
float mZ1, mZ2, m4l, pt4l; |
| 144 |
|
unsigned tZ1, tZ2; |
| 145 |
< |
double w; |
| 145 |
> |
double w, woff, werroff, won, werron; |
| 146 |
> |
if( ctrl.mc ) { |
| 147 |
> |
passtuple->Branch("gchannel", &gchannel); |
| 148 |
> |
passtuple->Branch("woff", &woff); |
| 149 |
> |
passtuple->Branch("werroff", &werroff); |
| 150 |
> |
passtuple->Branch("won", &won); |
| 151 |
> |
passtuple->Branch("werron", &werron); |
| 152 |
> |
} |
| 153 |
> |
passtuple->Branch("channel", &channel); |
| 154 |
|
passtuple->Branch("run", &run); |
| 155 |
|
passtuple->Branch("evt", &evt); |
| 156 |
|
passtuple->Branch("lumi", &lumi); |
| 161 |
|
passtuple->Branch("m4l", &m4l); |
| 162 |
|
passtuple->Branch("pt4l", &pt4l); |
| 163 |
|
passtuple->Branch("w", &w); |
| 164 |
+ |
|
| 165 |
|
initCiCSelection(); |
| 166 |
|
if(ctrl.eleSele=="lik") initLikSelection(); |
| 167 |
|
if(ctrl.eleSele=="bdt") initBDTSelection(); |
| 168 |
|
// if(ctrl.eleSele=="si" ) initSiMVAElectronSelection(); |
| 169 |
|
initRunLumiRangeMap(); |
| 170 |
+ |
initEfficiencyWeights(); |
| 171 |
+ |
|
| 172 |
+ |
|
| 173 |
+ |
|
| 174 |
|
|
| 175 |
|
// |
| 176 |
|
// Loop |
| 184 |
|
|
| 185 |
|
// Data structures to store info from TTrees |
| 186 |
|
mithep::TEventInfo *info = new mithep::TEventInfo(); |
| 187 |
+ |
mithep::TGenInfo *ginfo = new mithep::TGenInfo(); |
| 188 |
|
TClonesArray *electronArr = new TClonesArray("mithep::TElectron"); |
| 189 |
|
TClonesArray *muonArr = new TClonesArray("mithep::TMuon"); |
| 190 |
|
|
| 191 |
|
|
| 147 |
– |
TBranch *infoBr; |
| 148 |
– |
TBranch *electronBr; |
| 149 |
– |
TBranch *muonBr; |
| 150 |
– |
|
| 151 |
– |
// chain->SetBranchStatus("*",0); //disable all branches |
| 152 |
– |
// chain->SetBranchStatus("Info", 1); |
| 153 |
– |
// chain->SetBranchStatus("Muon", 1); |
| 154 |
– |
// chain->SetBranchStatus("Electron", 1); |
| 155 |
– |
|
| 192 |
|
chain->SetBranchAddress("Info", &info); |
| 193 |
|
chain->SetBranchAddress("Electron", &electronArr); |
| 194 |
|
chain->SetBranchAddress("Muon", &muonArr); |
| 195 |
+ |
if( ctrl.mc ) { chain->SetBranchAddress("Gen", &ginfo);} |
| 196 |
|
|
| 197 |
|
|
| 198 |
|
int count=0, pass=0; |
| 199 |
< |
cout << "nEntries: " << chain->GetEntries() << endl; |
| 200 |
< |
cerr << "nEntries: " << chain->GetEntries() << endl; |
| 199 |
> |
float denom[3]={0.,0.,0.}; |
| 200 |
> |
float numer[3]={0.,0.,0.}; |
| 201 |
> |
float numercorr[3]={0.,0.,0.}; |
| 202 |
> |
|
| 203 |
|
UInt_t imax = chain->GetEntries(); |
| 204 |
+ |
cout << "nEntries: " << imax << endl; |
| 205 |
+ |
|
| 206 |
+ |
|
| 207 |
|
for(UInt_t ientry=0; ientry<imax; ientry++) { |
| 208 |
|
chain->GetEntry(ientry); |
| 209 |
|
if(!(ientry%100000)) cerr << "entry: " << ientry << endl; |
| 254 |
|
#endif |
| 255 |
|
|
| 256 |
|
unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple); |
| 257 |
+ |
// if( !evtfail ) { |
| 258 |
+ |
// cout << "run: " << info->runNum |
| 259 |
+ |
// << "\tevt: " << info->evtNum |
| 260 |
+ |
// << "\tlumi: " << info->lumiSec |
| 261 |
+ |
// << "\tfile: " << chain->GetFile()->GetEndpointUrl()->GetFile() |
| 262 |
+ |
// << endl; |
| 263 |
+ |
// } |
| 264 |
|
|
| 265 |
|
#ifdef THEIR_EVENTS |
| 266 |
|
if(evtfail!=0) cout << "Failed!" << endl; |
| 268 |
|
|
| 269 |
|
h_evtfail->Fill( evtfail, eventweight ); |
| 270 |
|
|
| 271 |
+ |
if( ctrl.mc ) { |
| 272 |
+ |
// BR is implicit here ... |
| 273 |
+ |
int gchannel = getChannel(ginfo); |
| 274 |
+ |
if( gchannel >= 0 ) { |
| 275 |
+ |
// denom[gchannel] += eventweight; |
| 276 |
+ |
// if( !evtfail ) numer[gchannel] += eventweight; |
| 277 |
+ |
denom[gchannel] += 1; |
| 278 |
+ |
if( !evtfail ) numer[gchannel] += 1; |
| 279 |
+ |
|
| 280 |
+ |
// get eff weight for current evt |
| 281 |
+ |
if( !evtfail ) { |
| 282 |
+ |
double thewoff, thewon; |
| 283 |
+ |
passtuple->SetBranchAddress("woff", &thewoff); // must put this here ... |
| 284 |
+ |
passtuple->SetBranchAddress("won", &thewon); // must put this here ... |
| 285 |
+ |
// cout << "thewoff: " << thewoff << "\tthewon: " << thewon << endl; |
| 286 |
+ |
int ret = passtuple->GetEntry(passtuple->GetEntries()-1); |
| 287 |
+ |
numercorr[gchannel] += thewoff*thewon; |
| 288 |
+ |
} |
| 289 |
+ |
} |
| 290 |
+ |
} |
| 291 |
+ |
|
| 292 |
|
count++; |
| 293 |
|
if( !evtfail ) pass++; |
| 294 |
|
|
| 329 |
|
// map<TString,TMVA::Reader*>::iterator it; |
| 330 |
|
// for(it=readers.begin(); it!=readers.end(); it++) delete (*it).second; |
| 331 |
|
|
| 332 |
+ |
|
| 333 |
+ |
float acc_4e = numer[0]/denom[0]; |
| 334 |
+ |
float acc_4mu = numer[1]/denom[1]; |
| 335 |
+ |
float acc_2e2mu = numer[2]/denom[2]; |
| 336 |
+ |
float acc_err_4e = (1./denom[0])*sqrt(numer[0]*(1-acc_4e)); |
| 337 |
+ |
float acc_err_4mu = (1./denom[1])*sqrt(numer[1]*(1-acc_4mu)); |
| 338 |
+ |
float acc_err_2e2mu = (1./denom[2])*sqrt(numer[2]*(1-acc_2e2mu)); |
| 339 |
+ |
|
| 340 |
+ |
float acc_c_4e = numercorr[0]/denom[0]; |
| 341 |
+ |
float acc_c_4mu = numercorr[1]/denom[1]; |
| 342 |
+ |
float acc_c_2e2mu = numercorr[2]/denom[2]; |
| 343 |
+ |
float acc_c_err_4e = (1./denom[0])*sqrt(numer[0]*(1-acc_4e)); |
| 344 |
+ |
float acc_c_err_4mu = (1./denom[1])*sqrt(numer[1]*(1-acc_4mu)); |
| 345 |
+ |
float acc_c_err_2e2mu = (1./denom[2])*sqrt(numer[2]*(1-acc_2e2mu)); |
| 346 |
+ |
|
| 347 |
+ |
if( ctrl.mc ) { |
| 348 |
+ |
cerr << "---------------------" << endl; |
| 349 |
+ |
// cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << numer[0] << endl; // numer already weighted ... |
| 350 |
+ |
// cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << numer[1] << endl; |
| 351 |
+ |
// cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << numer[2] << endl; |
| 352 |
+ |
cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << acc_4e << " +/- " << acc_err_4e << endl; |
| 353 |
+ |
cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << acc_4mu << " +/- " << acc_err_4mu << endl; |
| 354 |
+ |
cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_2e2mu << " +/- " << acc_err_2e2mu << endl; |
| 355 |
+ |
cerr << "-------corrected------" << endl; |
| 356 |
+ |
cerr << "acc_c_4e: " << fixed << setw(8) << setprecision(7) << acc_c_4e << " +/- " << acc_c_err_4e << endl; |
| 357 |
+ |
cerr << "acc_c_4mu: " << fixed << setw(8) << setprecision(7) << acc_c_4mu << " +/- " << acc_c_err_4mu << endl; |
| 358 |
+ |
cerr << "acc_c_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_c_2e2mu << " +/- " << acc_c_err_2e2mu << endl; |
| 359 |
+ |
cerr << "---------------------" << endl; |
| 360 |
+ |
} |
| 361 |
+ |
|
| 362 |
|
cerr << "done!" << endl; |
| 363 |
|
} |
| 364 |
|
|
