| 1 |
– |
// system headers |
| 2 |
– |
#include <map> |
| 3 |
– |
#include <utility> |
| 4 |
– |
|
| 5 |
– |
// mit headers |
| 6 |
– |
#include "Vertex.h" |
| 7 |
– |
|
| 8 |
– |
// 4L stuff |
| 9 |
– |
#include "SelectionStatus.h" |
| 10 |
– |
#include "EventData.h" |
| 11 |
– |
#include "SimpleLepton.h" |
| 12 |
– |
#include "EfficiencyWeightsInterface.h" |
| 13 |
– |
#include "ElectronSelection.h" |
| 14 |
– |
#include "MuonSelection.h" |
| 15 |
– |
#include "IsolationSelection.h" |
| 1 |
|
#include "ReferenceSelection.h" |
| 17 |
– |
#include "Selection.h" |
| 18 |
– |
#include "CommonDefs.h" |
| 19 |
– |
#include "SelectionDefs.h" |
| 20 |
– |
#include "FSR.h" |
| 21 |
– |
#include "SelectionFuncs.h" |
| 22 |
– |
#include "ElectronMomentumCorrection.h" |
| 23 |
– |
#include "MuonMomentumCorrection.h" |
| 24 |
– |
|
| 2 |
|
|
| 3 |
|
extern vector<SimpleLepton> failingLeptons; |
| 4 |
|
extern vector<SimpleLepton> passingLeptons; |
| 8 |
|
extern map<TString, map<TString,int>* > counts; |
| 9 |
|
extern ElectronMomentumCorrection electron_momentum_correction; |
| 10 |
|
extern MuCorr *muCorr; |
| 11 |
+ |
extern TrigInfo ti; |
| 12 |
|
|
| 13 |
|
//---------------------------------------------------------------------------------------- |
| 14 |
|
void updateSimpleLepton(SimpleLepton &tmplep); |
| 37 |
– |
//---------------------------------------------------------------------------------------- |
| 38 |
– |
void increment(TString cutName, int z1type=0, int z2type=0) |
| 39 |
– |
//---------------------------------------------------------------------------------------- |
| 40 |
– |
// increment counters that keep track up event counts at each cut |
| 41 |
– |
{ |
| 42 |
– |
assert(counts.find(cutName) != counts.end()); |
| 43 |
– |
(*(counts[cutName]))["all"] ++; |
| 44 |
– |
if(z1type!=0 && z2type!=0) { |
| 45 |
– |
assert(counts.find(cutName) != counts.end()); |
| 46 |
– |
(*(counts[cutName]))[getChannel(z1type,z2type)] ++; |
| 47 |
– |
} |
| 48 |
– |
} |
| 15 |
|
//-------------------------------------------------------------------------------------------------- |
| 16 |
|
EventData apply_HZZ4L_reference_selection(ControlFlags &ctrl, // input control |
| 17 |
|
const EventHeader *info, // input event info |
| 52 |
|
) |
| 53 |
|
//-------------------------------------------------------------------------------------------------- |
| 54 |
|
{ |
| 89 |
– |
increment("start"); |
| 55 |
|
EventData ret; |
| 56 |
|
failingLeptons.clear(); |
| 57 |
|
passingLeptons.clear(); |
| 70 |
|
// correct muon momentum |
| 71 |
|
if(ctrl.correct_muon_momentum) { |
| 72 |
|
assert(muCorr->corr2011 && muCorr->corr2012); |
| 73 |
+ |
if(ctrl.debug) cout << "mu corr: " << endl; |
| 74 |
|
for(int i=0; i<muonArr->GetEntries(); i++) { |
| 75 |
|
const Muon *const_mu = (Muon*)((*muonArr)[i]); |
| 76 |
|
Muon *mu = const_cast<Muon*>(const_mu); |
| 77 |
+ |
float ptBefore = mu->Pt(); |
| 78 |
|
correct_muon_momentum(ctrl,muCorr,mu,info->RunNum()); |
| 79 |
+ |
if(ctrl.debug) cout << " " << setw(12) << ptBefore << " --> " << setw(12) << mu->Pt() << " (" << fabs(ptBefore-mu->Pt())/ptBefore << ")" << endl; |
| 80 |
|
} |
| 81 |
|
} |
| 82 |
|
|
| 94 |
|
// Trigger Selection |
| 95 |
|
//*********************************************************** |
| 96 |
|
if( passes_HLT ) { |
| 97 |
< |
increment("trigger"); |
| 97 |
> |
increment(counts,"trigger"); |
| 98 |
|
} else { |
| 99 |
|
if(ctrl.debug) cout << "fails trigger" << endl; |
| 100 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 110 |
|
if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl; |
| 111 |
|
for(int i=0; i<muonArr->GetEntries(); i++) |
| 112 |
|
{ |
| 113 |
< |
const Muon *mu = (Muon*)((*muonArr)[i]); |
| 114 |
< |
|
| 113 |
> |
const Muon *mu = (Muon*)((*muonArr)[i]); |
| 114 |
> |
|
| 115 |
|
SelectionStatus musel; |
| 148 |
– |
if(ctrl.debug) cout << "musel.status before anything: " << musel.getStatus() << endl; |
| 149 |
– |
|
| 116 |
|
musel |= (*MuonPreSelector)(ctrl,mu,vtx,pfCandidates); |
| 151 |
– |
if(ctrl.debug) cout << "musel.status after presel: " << musel.getStatus() << endl; |
| 152 |
– |
if( !(musel.getStatus() & SelectionStatus::PRESELECTION) ) continue; |
| 153 |
– |
if( ctrl.debug ) cout << endl; |
| 154 |
– |
|
| 155 |
– |
musel |= (*MuonIDSelector)(ctrl,mu,vtx,pfCandidates ); |
| 156 |
– |
if(ctrl.debug) cout << "musel.status after ID: " << musel.getStatus() << endl; |
| 157 |
– |
if( ctrl.debug ) cout << endl; |
| 117 |
|
|
| 118 |
< |
// if we *do* FSR, isolation is calculated only after FSR recovery |
| 119 |
< |
if( !(ctrl.doFSR) ) { |
| 120 |
< |
musel |= (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto); |
| 121 |
< |
if(ctrl.debug) cout << "musel.status after iso: " << musel.getStatus() << endl; |
| 118 |
> |
if(ctrl.fakeScheme == "none") { // regular selection |
| 119 |
> |
if( !musel.passPre() ) continue; |
| 120 |
> |
} else { // fake denominator selection |
| 121 |
> |
SelectionStatus denomSel; |
| 122 |
> |
denomSel |= muonReferencePreSelection(ctrl,mu,vtx,pfCandidates); |
| 123 |
> |
if(ctrl.muSele != "none") { // don't apply *any* preselection, to get a super hi stat fake sample for BDT training |
| 124 |
> |
if( !denomSel.passPre() ) continue; |
| 125 |
> |
} |
| 126 |
|
} |
| 127 |
|
|
| 128 |
< |
if( ctrl.debug ) |
| 129 |
< |
cout << "muon:: pt: " << mu->Pt() << "\teta: " << mu->Eta() << "\tstatus: " << hex << musel.getStatus() << dec << endl; |
| 128 |
> |
if(ctrl.debug) cout << "muon:: pt: " << mu->Pt() << "\teta: " << mu->Eta() << endl; |
| 129 |
> |
|
| 130 |
> |
musel |= (*MuonIDSelector)(ctrl,mu,vtx,pfCandidates ); |
| 131 |
> |
|
| 132 |
> |
// NOTE: if we do FSR this is *changed* later on |
| 133 |
> |
musel |= (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto); |
| 134 |
> |
if(ctrl.debug) cout << "mu status after iso (before fsr) " << musel.getStatus() << endl; |
| 135 |
|
|
| 136 |
|
SimpleLepton tmplep; |
| 137 |
|
tmplep.vec.SetPtEtaPhiM(mu->Pt(), mu->Eta(), mu->Phi(), MUON_MASS); |
| 155 |
|
const Electron *const_ele = (Electron*)((*electronArr)[i]); |
| 156 |
|
Electron *ele = const_cast<Electron*>(const_ele); |
| 157 |
|
|
| 158 |
+ |
if(ctrl.debug) cout << setprecision(8) << "el:: scEt " << setw(12) << ele->SCluster()->Et() << " P: " |
| 159 |
+ |
<< setw(12) << ele->P() << " pT: " << setw(12) << ele->Pt() << " scEta " << setw(12) << ele->SCluster()->Eta() << setprecision(5) << endl; |
| 160 |
+ |
|
| 161 |
|
// evaluate MVA *before* regression correction (but cut on mvaVal *after*) |
| 162 |
|
double mvaVal = getElectronIDMVAval(ctrl, ele, vtx); |
| 163 |
|
|
| 164 |
< |
electron_momentum_correction.correct_electron_momentum(ctrl, ele,puEnergyDensity->At(0)->RhoKt6PFJets(),vtxArr->GetEntries()); |
| 164 |
> |
electron_momentum_correction.correct_electron_momentum(ctrl, ele, info, puEnergyDensity->At(0)->RhoKt6PFJets(), vtxArr->GetEntries()); |
| 165 |
|
|
| 166 |
< |
SelectionStatus elesel; |
| 167 |
< |
if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl; |
| 166 |
> |
if(ctrl.debug) cout << setprecision(8) << " corr el: scEt " << setw(12) << ele->SCluster()->Et() << " P: " |
| 167 |
> |
<< setw(12) << ele->P() << " pT: " << setw(12) << ele->Pt() << " scEta " << setw(12) << ele->SCluster()->Eta() << setprecision(5) << endl; |
| 168 |
|
|
| 169 |
+ |
SelectionStatus elesel; |
| 170 |
|
elesel |= (*ElectronPreSelector)(ctrl,ele,vtx); |
| 199 |
– |
if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl; |
| 200 |
– |
if( !(elesel.getStatus() & SelectionStatus::PRESELECTION) ) continue; |
| 171 |
|
|
| 172 |
< |
// elesel |= (*ElectronIDSelector)(ctrl,ele,vtx); |
| 172 |
> |
if(ctrl.fakeScheme == "none") { // regular selection |
| 173 |
> |
if( !elesel.passPre() ) continue; |
| 174 |
> |
} else { // fake denominator selection |
| 175 |
> |
SelectionStatus denomSel; |
| 176 |
> |
denomSel |= electronReferencePreSelection(ctrl,ele,vtx); |
| 177 |
> |
if(ctrl.eleSele != "none") { // don't apply *any* preselection, to get a super hi stat fake sample for BDT training |
| 178 |
> |
if( !denomSel.passPre() ) continue; |
| 179 |
> |
} |
| 180 |
> |
} |
| 181 |
> |
|
| 182 |
> |
|
| 183 |
|
elesel |= passElectronIDMVA(ctrl, mvaVal, ele); |
| 204 |
– |
if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl; |
| 184 |
|
|
| 185 |
< |
if( !(ctrl.doFSR) ) { |
| 186 |
< |
elesel |= (*ElectronIsoSelector)(ctrl,ele,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto); |
| 187 |
< |
if( ctrl.debug ) cout << "--> status after iso: " << hex << elesel.getStatus() << dec << endl; |
| 209 |
< |
} |
| 185 |
> |
// NOTE: if we do FSR this is *changed* later on |
| 186 |
> |
elesel |= (*ElectronIsoSelector)(ctrl,ele,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto); |
| 187 |
> |
if( ctrl.debug ) cout << "el status after iso (no fsr) " << hex << elesel.getStatus() << dec << endl; |
| 188 |
|
|
| 211 |
– |
if( ctrl.debug ) |
| 212 |
– |
cout << "\tscEt: " << ele->SCluster()->Et() |
| 213 |
– |
<< "\tscEta: " << ele->SCluster()->Eta() |
| 214 |
– |
<< "\tstatus: " << hex << elesel.getStatus() << dec |
| 215 |
– |
<< endl; |
| 216 |
– |
|
| 189 |
|
SimpleLepton tmplep; |
| 190 |
|
tmplep.vec.SetPtEtaPhiM( ele->Pt(), ele->Eta(), ele->Phi(), ELECTRON_MASS ); |
| 191 |
|
tmplep.type = 11; |
| 245 |
|
bool erase_this_electron=false; |
| 246 |
|
for (vector<SimpleLepton>::iterator it2=lepvec.begin(); |
| 247 |
|
it2 != lepvec.end(); it2++ ) { |
| 248 |
< |
if ( it2 == it1 ) continue; |
| 249 |
< |
if ( abs(it2->type) != 13 ) continue; |
| 250 |
< |
// from ANL: |
| 251 |
< |
// ??- if( !it2->status.passPre() || (!it2->isPFMuon && !it2->isGlobalMuon) ) continue; |
| 280 |
< |
if( !it2->status.passPre() ) continue; |
| 281 |
< |
if( !it2->isTight ) continue; // NOTE: isTight is set to it2->isPFMuon || it2->isGlobalMuon |
| 248 |
> |
if ( it2 == it1 ) continue; |
| 249 |
> |
if ( abs(it2->type) != 13 ) continue; |
| 250 |
> |
if( !it2->status.passPre() ) continue; |
| 251 |
> |
if( !it2->isTight ) continue; // NOTE: isTight is set to it2->isPFMuon || it2->isGlobalMuon |
| 252 |
|
|
| 253 |
|
TVector3 mvec = it2->vec.Vect(); |
| 254 |
|
if ( evec.DrEtaPhi(mvec) < 0.05 ) { |
| 265 |
|
lepvec.erase(electrons_to_erase[i]); |
| 266 |
|
} |
| 267 |
|
|
| 268 |
+ |
// fill passing and failing leptons for fakes |
| 269 |
+ |
// (NOTE: isolation hasn't been corrected post-fsr, so this isn't strictly correct) |
| 270 |
+ |
// also note: sip cut not applied |
| 271 |
+ |
for (int i=0; i<lepvec.size(); i++ ) { |
| 272 |
+ |
if( lepvec[i].isLoose ) |
| 273 |
+ |
passingLeptons.push_back(lepvec[i]); |
| 274 |
+ |
else |
| 275 |
+ |
failingLeptons.push_back(lepvec[i]); |
| 276 |
+ |
} |
| 277 |
+ |
|
| 278 |
|
//******************************************************** |
| 279 |
< |
// Step 3: Good Leptons (remove non-id-and-pre'd leptons) |
| 279 |
> |
// Step 3: Good Leptons (remove non-id-and-pre'd leptons, and apply |ip3ds| < 4) |
| 280 |
|
//******************************************************** |
| 281 |
|
vector<double> pt_of_leptons_to_erase; |
| 282 |
|
for (int i=0; i<lepvec.size(); i++ ) { |
| 289 |
|
<< "\tpt: " << lepvec[i].vec.Pt() |
| 290 |
|
<< "\teta: " << lepvec[i].vec.Eta() |
| 291 |
|
<< endl; |
| 292 |
< |
failingLeptons.push_back(lepvec[i]); // these should pass preselection |
| 292 |
> |
// failingLeptons.push_back(lepvec[i]); // these should pass preselection |
| 293 |
|
} else { |
| 294 |
< |
passingLeptons.push_back(lepvec[i]); |
| 294 |
> |
// passingLeptons.push_back(lepvec[i]); |
| 295 |
|
} |
| 296 |
|
if( !already_pushed && fabs(lepvec[i].ip3dSig)>4 ) |
| 297 |
|
pt_of_leptons_to_erase.push_back(lepvec[i].vec.Pt()); |
| 339 |
|
//******************************************************** |
| 340 |
|
// Step 4: Z candidate preselection |
| 341 |
|
//******************************************************** |
| 342 |
< |
vector<pair<int,int> > ZCandidates; |
| 343 |
< |
vector<pair<SimpleLepton,SimpleLepton> > ZCandidatesLeptons; |
| 342 |
> |
vector<pair<int,int> > ZCandidates; // indices in lepvec of the two leptons |
| 343 |
> |
vector<pair<SimpleLepton,SimpleLepton> > ZCandidatesLeptons; // fsr-corrected leptons that form the Z |
| 344 |
|
vector<pair<SelectionStatus,SelectionStatus> > ZCandidatesSelectionStatus; |
| 345 |
|
for(int i = 0; i < lepvec.size(); ++i) { |
| 346 |
|
for(int j = i+1; j < lepvec.size(); ++j) { |
| 355 |
|
if( ctrl.doFSR ) { |
| 356 |
|
if(ctrl.debug) cout << endl << "----------------> FSR ("<<i<<","<<j<<") <----------------------" << endl; |
| 357 |
|
photonsToVeto.clear(); |
| 358 |
< |
TLorentzVector old_lepvec_i = lepvec[i].vec; |
| 359 |
< |
TLorentzVector old_lepvec_j = lepvec[j].vec; |
| 380 |
< |
float old_pt_i = lepvec[i].vec.Pt(); |
| 358 |
> |
// copies of lepvec that have 1) the photons added to the momenta and 2) the fsrRecoveryAttempted flag set |
| 359 |
> |
float old_pt_i = lepvec[i].vec.Pt(); |
| 360 |
|
float old_pt_j = lepvec[j].vec.Pt(); |
| 361 |
|
float old_M = Zvec.M(); |
| 383 |
– |
|
| 384 |
– |
bool recovery_index_i = false; |
| 385 |
– |
bool recovery_index_j = false; |
| 362 |
|
|
| 363 |
+ |
const ChargedParticle *lepton(NULL); |
| 364 |
|
if(ctrl.debug) cout << "i: " << i << endl; |
| 365 |
< |
if( abs(lepvec[i].type) == 13 ) { |
| 366 |
< |
const Muon *mu = (Muon*)((*muonArr)[lepvec[i].index]); |
| 367 |
< |
Muon * newmu = const_cast<Muon *>(mu); |
| 368 |
< |
if( recover_typeI_Photon( ctrl, ret, newmu, i, lepvec_i, pfCandidates, electronArr, &Zvec, photonsToVeto ) ){ |
| 392 |
< |
if(ctrl.debug) cout << " FSR TYPEI :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl; |
| 393 |
< |
recovery_index_i = true; |
| 394 |
< |
} |
| 395 |
< |
if( recover_typeII_Photon( ctrl, ret, newmu, i, lepvec_i, pfCandidates ) ){ |
| 396 |
< |
if(ctrl.debug) cout << " FSR TYPEII :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl; |
| 397 |
< |
recovery_index_i = true; |
| 398 |
< |
} |
| 399 |
< |
} else { |
| 400 |
< |
const Electron *el = (Electron*)((*electronArr)[lepvec[i].index]); |
| 401 |
< |
Electron* newel = const_cast<Electron*>(el); |
| 402 |
< |
if( recover_typeI_Photon( ctrl, ret, newel, i, lepvec_i, pfCandidates, electronArr, &Zvec, photonsToVeto ) ){ |
| 403 |
< |
if(ctrl.debug) cout << " FSR TYPEI :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl; |
| 404 |
< |
recovery_index_i = true; |
| 405 |
< |
} |
| 406 |
< |
} |
| 365 |
> |
if(abs(lepvec[i].type) == 13) lepton = dynamic_cast<const ChargedParticle *>((*muonArr)[lepvec[i].index]); |
| 366 |
> |
else lepton = dynamic_cast<const ChargedParticle *>((*electronArr)[lepvec[i].index]); |
| 367 |
> |
pair<TLorentzVector,int> photon_i = findFsrPhoton(ctrl, ret, lepton, i, lepvec_i, pfCandidates, electronArr, &Zvec); |
| 368 |
> |
if(photon_i.second>=0 && ctrl.debug) cout << " FSR :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl; |
| 369 |
|
|
| 370 |
+ |
lepton = NULL; |
| 371 |
|
if(ctrl.debug) cout << "j: " << j << endl; |
| 372 |
< |
ChargedParticle *cp; |
| 373 |
< |
if( abs(lepvec[j].type) == 13 ) { |
| 374 |
< |
const Muon *mu = (Muon*)((*muonArr)[lepvec[j].index]); |
| 375 |
< |
Muon * newmu = const_cast<Muon *>(mu); |
| 376 |
< |
if( recover_typeI_Photon( ctrl, ret, newmu, j, lepvec_j, pfCandidates, electronArr, &Zvec, photonsToVeto ) ){ |
| 377 |
< |
if(ctrl.debug) cout << " FSR TYPEI :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl; |
| 378 |
< |
recovery_index_j = true; |
| 379 |
< |
} |
| 380 |
< |
if( recover_typeII_Photon( ctrl, ret, newmu, j, lepvec_j, pfCandidates ) ){ |
| 381 |
< |
if(ctrl.debug) cout << " FSR TYPEII :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl; |
| 382 |
< |
recovery_index_j = true; |
| 383 |
< |
} |
| 384 |
< |
} else { |
| 385 |
< |
const Electron *el = (Electron*)((*electronArr)[lepvec[j].index]); |
| 386 |
< |
Electron* newel = const_cast<Electron*>(el); |
| 387 |
< |
if( recover_typeI_Photon( ctrl, ret, newel, j, lepvec_j, pfCandidates, electronArr, &Zvec, photonsToVeto ) ){ |
| 388 |
< |
if(ctrl.debug) cout << " FSR TYPEI :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl; |
| 389 |
< |
recovery_index_j = true; |
| 390 |
< |
} |
| 391 |
< |
} |
| 429 |
< |
|
| 430 |
< |
if(recovery_index_i && recovery_index_j){ |
| 431 |
< |
if(abs((old_lepvec_i+lepvec_j[j].vec).M() - Z_MASS) > abs((lepvec_i[i].vec+old_lepvec_j).M() - Z_MASS)){ |
| 432 |
< |
//remove the photon from lepvec_j |
| 433 |
< |
lepvec_j[j].fsrRecoveryAttempted = false; |
| 434 |
< |
lepvec_j[j].vec = old_lepvec_j; |
| 435 |
< |
} |
| 436 |
< |
else { |
| 437 |
< |
lepvec_i[i].fsrRecoveryAttempted = false; |
| 438 |
< |
lepvec_i[i].vec = old_lepvec_i; |
| 372 |
> |
if( abs(lepvec[j].type) == 13 ) lepton = dynamic_cast<const ChargedParticle *>((*muonArr)[lepvec[j].index]); |
| 373 |
> |
else lepton = dynamic_cast<const ChargedParticle *>((*electronArr)[lepvec[j].index]); |
| 374 |
> |
pair<TLorentzVector,int> photon_j = findFsrPhoton(ctrl, ret, lepton, j, lepvec_j, pfCandidates, electronArr, &Zvec); |
| 375 |
> |
if(photon_j.second>=0 && ctrl.debug) cout << " FSR :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl; |
| 376 |
> |
|
| 377 |
> |
bool useI(false),useJ(false); |
| 378 |
> |
if(photon_i.second >= 0) { // photon for i |
| 379 |
> |
if(photon_j.second < 0) { |
| 380 |
> |
useI = true; |
| 381 |
> |
} else { // both have photons |
| 382 |
> |
TLorentzVector lep_i(lepvec_i[i].vec),lep_j(lepvec_j[j].vec); |
| 383 |
> |
float oldM = (lep_i + lep_j).M(); |
| 384 |
> |
float newM_i = (lep_i + photon_i.first + lep_j).M(); |
| 385 |
> |
float newM_j = (lep_i + lep_j + photon_j.first).M(); |
| 386 |
> |
if(ctrl.debug) cout << " two FSRs " << setw(12) << oldM << setw(12) << newM_i << setw(12) << newM_j << endl; |
| 387 |
> |
if(fabs(newM_i - Z_MASS) < fabs(newM_j - Z_MASS)) { // if mass with photon i is closer to z pole than with photon j |
| 388 |
> |
useI = true; |
| 389 |
> |
} else { |
| 390 |
> |
useJ = true; |
| 391 |
> |
} |
| 392 |
|
} |
| 393 |
+ |
} else if(photon_j.second >= 0) { // photon for j |
| 394 |
+ |
useJ = true; |
| 395 |
|
} |
| 396 |
|
|
| 397 |
< |
// |
| 398 |
< |
// now fix isolation |
| 399 |
< |
// |
| 397 |
> |
assert(!(useI && useJ)); |
| 398 |
> |
if(useI) { |
| 399 |
> |
addPhotonToEventData(ret,photon_i.first); |
| 400 |
> |
lepvec_i[i].vec += photon_i.first; |
| 401 |
> |
lepvec_i[i].fsrRecoveryAttempted = true; |
| 402 |
> |
photonsToVeto.push_back((const PFCandidate*)(pfCandidates->At(photon_i.second))); |
| 403 |
> |
} else if(useJ) { |
| 404 |
> |
addPhotonToEventData(ret,photon_j.first); |
| 405 |
> |
lepvec_j[j].vec += photon_j.first; |
| 406 |
> |
lepvec_j[j].fsrRecoveryAttempted = true; |
| 407 |
> |
photonsToVeto.push_back((const PFCandidate*)(pfCandidates->At(photon_j.second))); |
| 408 |
> |
} |
| 409 |
> |
|
| 410 |
> |
// recompute isolation excluding the fsr photons |
| 411 |
|
if( abs(lepvec[i].type) == 11 ) { |
| 412 |
|
const Electron *el = (Electron*)((*electronArr)[lepvec_i[i].index]); |
| 413 |
< |
lepvec_i[i].status |= |
| 448 |
< |
(*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto); |
| 413 |
> |
lepvec_i[i].status |= (*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto); |
| 414 |
|
|
| 415 |
|
} else { |
| 416 |
|
const Muon *mu = (Muon*)((*muonArr)[lepvec_i[i].index]); |
| 417 |
< |
lepvec_i[i].status |= |
| 453 |
< |
(*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto); |
| 417 |
> |
lepvec_i[i].status |= (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto); |
| 418 |
|
} |
| 455 |
– |
|
| 419 |
|
updateSimpleLepton(lepvec_i[i]); |
| 420 |
|
|
| 421 |
|
if( abs(lepvec[j].type) == 11 ) { |
| 422 |
|
const Electron *el = (Electron*)((*electronArr)[lepvec_j[j].index]); |
| 423 |
< |
lepvec_j[j].status |= |
| 461 |
< |
(*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto); |
| 423 |
> |
lepvec_j[j].status |= (*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto); |
| 424 |
|
} else { |
| 425 |
|
const Muon *mu = (Muon*)((*muonArr)[lepvec_j[j].index]); |
| 426 |
< |
lepvec_j[j].status |= |
| 465 |
< |
(*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto); |
| 426 |
> |
lepvec_j[j].status |= (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto); |
| 427 |
|
} |
| 467 |
– |
|
| 428 |
|
updateSimpleLepton(lepvec_j[j]); |
| 429 |
|
|
| 430 |
|
float new_M = (lepvec_i[i].vec+lepvec_j[j].vec).M(); |
| 431 |
|
float new_pt_i = lepvec_i[i].vec.Pt(); |
| 432 |
|
float new_pt_j = lepvec_j[j].vec.Pt(); |
| 433 |
< |
if( ctrl.debug ) { |
| 433 |
> |
if( ctrl.debug && ret.fsrPhotons.size() > 0) { |
| 434 |
|
cout << " ----> post FSR Z:"; |
| 435 |
< |
cout << " oldM: " << old_M << "\tnewM:" << new_M << endl; |
| 436 |
< |
cout << " old_pt_i: " << old_pt_i << "\tnew_pt_i:" << new_pt_i << endl; |
| 437 |
< |
cout << " old_pt_j: " << old_pt_j << "\tnew_pt_j:" << new_pt_j << endl; |
| 435 |
> |
cout << " oldM: " << setprecision(8) << old_M << "\tnewM:" << setprecision(8) << new_M << endl; |
| 436 |
> |
cout << " old_pt_i: " << setprecision(8) << old_pt_i << "\tnew_pt_i:" << setprecision(8) << new_pt_i << endl; |
| 437 |
> |
cout << " old_pt_j: " << setprecision(8) << old_pt_j << "\tnew_pt_j:" << setprecision(8) << new_pt_j << endl; |
| 438 |
|
} |
| 439 |
|
|
| 440 |
|
} // doFSR |
| 441 |
< |
|
| 441 |
> |
|
| 442 |
> |
// apply isolation criteria with recomputed isolation |
| 443 |
|
if( !(lepvec_i[i].status.loose()) || !(lepvec_j[j].status.loose()) ) { |
| 444 |
|
if(ctrl.debug) { |
| 445 |
|
bitset<16> tmpbits_i(lepvec_i[i].status.getStatus()),tmpbits_j(lepvec_j[j].status.getStatus()); |
| 447 |
|
} |
| 448 |
|
continue; |
| 449 |
|
} |
| 450 |
+ |
|
| 451 |
|
ZCandidates.push_back(pair<int,int> (i,j) ); |
| 452 |
|
ZCandidatesLeptons.push_back(pair<SimpleLepton,SimpleLepton> (lepvec_i[i],lepvec_j[j]) ); |
| 453 |
< |
if( ctrl.debug ) cout << "Z candidate ("<< i << "," << j << ")" |
| 492 |
< |
<< "\tmass: " << (lepvec_i[i].vec+lepvec_j[j].vec).M() << endl; |
| 453 |
> |
if( ctrl.debug ) cout << "Z candidate ("<< i << "," << j << ")" << "\tmass: " << (lepvec_i[i].vec+lepvec_j[j].vec).M() << endl; |
| 454 |
|
} |
| 455 |
|
} |
| 456 |
+ |
|
| 457 |
|
if( ZCandidates.size() > 0 ) { |
| 458 |
|
ret.status.selectionBits.flip(PASS_ZCANDIDATE); |
| 459 |
|
if( ctrl.debug ) cout << "event has >0 Z candidates (" << ZCandidates.size() << " of em)" << endl; |
| 463 |
|
} |
| 464 |
|
|
| 465 |
|
// |
| 466 |
< |
// !!!!!!!!!!!!!! Z1 SELECTED HERE |
| 466 |
> |
// Select Z1 |
| 467 |
|
// |
| 468 |
|
int best_Z1_index; |
| 469 |
|
float best_Z1_mass = 9999.; |
| 470 |
|
TLorentzVector Z1vec; |
| 471 |
|
for( int i=0; i<ZCandidates.size(); i++ ) { |
| 472 |
< |
// TLorentzVector tmpZ1vec = (lepvec[ZCandidates[i].first].vec) + |
| 511 |
< |
// (lepvec[ZCandidates[i].second].vec); |
| 512 |
< |
TLorentzVector tmpZ1vec = (ZCandidatesLeptons[i].first.vec) + |
| 513 |
< |
(ZCandidatesLeptons[i].second.vec); |
| 472 |
> |
TLorentzVector tmpZ1vec = (ZCandidatesLeptons[i].first.vec) + (ZCandidatesLeptons[i].second.vec); |
| 473 |
|
if( fabs(tmpZ1vec.M()-Z_MASS) < fabs(best_Z1_mass-Z_MASS) ) { |
| 474 |
|
best_Z1_index=i; |
| 475 |
|
best_Z1_mass=tmpZ1vec.M(); |
| 476 |
|
Z1vec = tmpZ1vec; |
| 477 |
|
} |
| 478 |
|
} |
| 520 |
– |
// update the leptons, damn FSR ... |
| 521 |
– |
// ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].first]); |
| 522 |
– |
// ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].second]); |
| 479 |
|
ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].first); |
| 480 |
|
ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].second); |
| 481 |
|
|
| 482 |
< |
cout << "best mZ1: " << best_Z1_mass << " from (" << |
| 483 |
< |
ZCandidates[best_Z1_index].first << "," << ZCandidates[best_Z1_index].second << ")" << endl; |
| 482 |
> |
if(ctrl.debug) |
| 483 |
> |
cout << "best mZ1: " << best_Z1_mass << " from (" << ZCandidates[best_Z1_index].first << "," << ZCandidates[best_Z1_index].second << ")" << endl; |
| 484 |
|
|
| 485 |
|
//****************************************************************************** |
| 486 |
|
// Step 6.3 : require Z1 with 40<m<120 |
| 496 |
|
// Step 6.3 : 4 good leptons |
| 497 |
|
//****************************************************************************** |
| 498 |
|
if( lepvec.size() >= 4 ) { |
| 499 |
< |
if( ctrl.debug) cout << "pass4L: " << info->EvtNum() << endl; |
| 499 |
> |
if( ctrl.debug) cout << "four leps" << endl; |
| 500 |
|
ret.status.selectionBits.flip(PASS_4L); |
| 501 |
|
} else { |
| 502 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 503 |
|
return ret; |
| 504 |
|
} |
| 549 |
– |
int nEl=0, nMu=0; |
| 550 |
– |
for( int i=0; i<4; i++ ) { |
| 551 |
– |
if(abs(lepvec[i].type) == 11 ) nEl++; |
| 552 |
– |
if(abs(lepvec[i].type) == 13 ) nMu++; |
| 553 |
– |
} |
| 505 |
|
|
| 506 |
|
//******************************************************** |
| 507 |
|
// Step 5: ZZ candidates |
| 508 |
|
//******************************************************** |
| 509 |
|
int nZZCandidates=0; |
| 510 |
|
vector<pair<int,int> > ZZCandidates; |
| 560 |
– |
int Z2type; |
| 511 |
|
if(ctrl.debug) cout << "looping over " << ZCandidates.size() << " Z candidates: " << endl; |
| 512 |
|
for(int z2index=0; z2index<ZCandidates.size(); ++z2index) { |
| 513 |
|
int z1index = best_Z1_index; |
| 514 |
< |
if ( z2index == z1index ) { if(ctrl.debug) cout << " Z2 candidate failing z2i = z1i" << endl; continue; } |
| 515 |
< |
if( ZCandidates[z1index].first == ZCandidates[z2index].first ) { if(ctrl.debug) cout << " Z2 candidate failing icheck1" << endl; continue; } |
| 516 |
< |
if( ZCandidates[z1index].first == ZCandidates[z2index].second ) { if(ctrl.debug) cout << " Z2 candidate failing icheck2" << endl; continue; } |
| 517 |
< |
if( ZCandidates[z1index].second == ZCandidates[z2index].first ) { if(ctrl.debug) cout << " Z2 candidate failing icheck3" << endl; continue; } |
| 518 |
< |
if( ZCandidates[z1index].second == ZCandidates[z2index].second ) { if(ctrl.debug) cout << " Z2 candidate failing icheck4" << endl; continue; } |
| 514 |
> |
if ( z2index == z1index ) { if(ctrl.debug) cout << " Z2 fails z2i = z1i" << endl; continue; } |
| 515 |
> |
if( ZCandidates[z1index].first == ZCandidates[z2index].first ) { if(ctrl.debug) cout << " Z2 fails icheck1" << endl; continue; } |
| 516 |
> |
if( ZCandidates[z1index].first == ZCandidates[z2index].second ) { if(ctrl.debug) cout << " Z2 fails icheck2" << endl; continue; } |
| 517 |
> |
if( ZCandidates[z1index].second == ZCandidates[z2index].first ) { if(ctrl.debug) cout << " Z2 fails icheck3" << endl; continue; } |
| 518 |
> |
if( ZCandidates[z1index].second == ZCandidates[z2index].second ) { if(ctrl.debug) cout << " Z2 fails icheck4" << endl; continue; } |
| 519 |
|
|
| 520 |
|
ZZCandidates.push_back(pair<int,int> (z1index,z2index)); |
| 571 |
– |
Z2type = abs(ZCandidatesLeptons[z2index].first.type); |
| 521 |
|
} |
| 522 |
< |
if( ZZCandidates.size() > 0 ) { |
| 574 |
< |
if( ctrl.debug) cout << "passZZ: " << info->EvtNum() << endl; |
| 522 |
> |
if(ZZCandidates.size() > 0) { |
| 523 |
|
ret.status.selectionBits.flip(PASS_ZZCANDIDATE); |
| 524 |
< |
if( ZZCandidates.size() > 1 ) Z2type=999; |
| 525 |
< |
if( ctrl.debug ) cout << "nZZcandidates: " << ZZCandidates.size() << endl; |
| 578 |
< |
if( ctrl.debug ) { |
| 579 |
< |
cout << "evt: " << info->EvtNum() << "\tnZZcandidates: " << ZZCandidates.size() |
| 580 |
< |
<< "\tZ1f: " << abs(ret.Z1leptons[0].type) << "\tZ2f: " << Z2type << endl; |
| 581 |
< |
cout << "------------------------------------------------------- ( un-fsrd )" << endl; |
| 582 |
< |
for( int l=0; l<lepvec.size(); l++ ) lepvec[l].print(); |
| 524 |
> |
if(ctrl.debug) { |
| 525 |
> |
cout << ZZCandidates.size() << " zz cands" << endl; |
| 526 |
|
cout << "-------------------------------------------------------" << endl; |
| 527 |
+ |
// for( int l=0; l<lepvec.size(); l++ ) lepvec[l].print(); |
| 528 |
+ |
for(unsigned ican=0; ican<ZCandidatesLeptons.size(); ican++) { |
| 529 |
+ |
ZCandidatesLeptons[ican].first.print(); |
| 530 |
+ |
ZCandidatesLeptons[ican].second.print(); |
| 531 |
+ |
} |
| 532 |
+ |
cout << "-------------------------------------------------------" << endl; |
| 533 |
+ |
if(ctrl.debug) |
| 534 |
+ |
for(unsigned iph=0; iph<ret.fsrPhotons.size(); iph++) |
| 535 |
+ |
cout << " fsr photon " |
| 536 |
+ |
<< setw(12) << ret.fsrPhotons[iph].vec.Pt() |
| 537 |
+ |
<< " eta: " << setw(12) << ret.fsrPhotons[iph].vec.Eta() |
| 538 |
+ |
<< " phi: " << setw(12) << ret.fsrPhotons[iph].vec.Phi() << endl; |
| 539 |
|
} |
| 540 |
|
} else { |
| 541 |
+ |
if(ctrl.debug) cout << "no zz cands" << endl; |
| 542 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 543 |
|
return ret; |
| 544 |
|
} |
| 545 |
|
|
| 546 |
|
// |
| 547 |
< |
// !!!!!!!!!!!!!! Z2 SELECTED HERE |
| 547 |
> |
// Select z2 |
| 548 |
|
// |
| 549 |
|
int best_Z2_index; |
| 550 |
|
float best_Z2_pt = -1.; |
| 557 |
|
best_Z2_pt=Z2.Pt(); |
| 558 |
|
} |
| 559 |
|
} |
| 560 |
< |
// update the leptons, damn FSR ... |
| 605 |
< |
// ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].first]); |
| 606 |
< |
// ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].second]); |
| 607 |
< |
ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].first); |
| 560 |
> |
ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].first); // push back the *fsr-corrected* leptons |
| 561 |
|
ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].second); |
| 562 |
|
|
| 563 |
|
int theZ1type = abs(ZCandidatesLeptons[best_Z1_index].first.type); |
| 566 |
|
if(ctrl.debug) cout << "best mZ2: " << (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M() << " from (" |
| 567 |
|
<< ZCandidates[best_Z2_index].first << "," << ZCandidates[best_Z2_index].second << ")" << endl; |
| 568 |
|
|
| 569 |
< |
increment("sfOsHiPt",theZ1type,theZ2type); |
| 569 |
> |
increment(counts,"sfOsHiPt",theZ1type,theZ2type); |
| 570 |
|
|
| 571 |
|
if(ctrl.debug) cout << "ZZ :: evt: " << info->EvtNum() |
| 572 |
|
<< "\tmZ1: " << (ret.Z1leptons[0].vec+ret.Z1leptons[1].vec).M() |
| 581 |
|
(ZCandidatesLeptons[best_Z2_index].second.vec); |
| 582 |
|
if( Z2vec.M() > 4 && Z2vec.M() < 120 ) { |
| 583 |
|
ret.status.selectionBits.flip(PASS_GOODZ2); |
| 584 |
< |
increment("4<mZ2<120",theZ1type,theZ2type); |
| 584 |
> |
increment(counts,"4<mZ2<120",theZ1type,theZ2type); |
| 585 |
|
} else { |
| 586 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 587 |
|
return ret; |
| 603 |
|
} |
| 604 |
|
if( nlep_above_10 > 1 && nlep_above_20 > 0 ) { |
| 605 |
|
ret.status.selectionBits.flip(PASS_ZZ_20_10); |
| 606 |
< |
increment("pt>20,10",theZ1type,theZ2type); |
| 606 |
> |
increment(counts,"pt>20,10",theZ1type,theZ2type); |
| 607 |
|
if( ctrl.debug ) cout << "passess 20/10 ..." << endl; |
| 608 |
|
} else { |
| 609 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 628 |
|
} |
| 629 |
|
if( !resonance ) { |
| 630 |
|
ret.status.selectionBits.flip(PASS_RESONANCE); |
| 631 |
< |
increment("mll>4",theZ1type,theZ2type); |
| 631 |
> |
increment(counts,"mll>4",theZ1type,theZ2type); |
| 632 |
|
if( ctrl.debug ) cout << "\tpasses resonance killing ... " << endl; |
| 633 |
|
} else { |
| 634 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 655 |
|
if( zzvec.M() > 70. ) { |
| 656 |
|
if(ctrl.debug) cout << "passes mzz > 70, mzz: " << zzvec.M() << endl; |
| 657 |
|
ret.status.selectionBits.flip(PASS_m4l_GT_70); |
| 658 |
< |
increment("m4l>70",theZ1type,theZ2type); |
| 658 |
> |
increment(counts,"m4l>70",theZ1type,theZ2type); |
| 659 |
|
} else { |
| 660 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 661 |
|
return ret; |
| 663 |
|
|
| 664 |
|
if( (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M() > 12 ){ |
| 665 |
|
if(ctrl.debug) cout << "passes mZ2 > 12" << endl; |
| 666 |
< |
increment("mZ2>12",theZ1type,theZ2type); |
| 666 |
> |
increment(counts,"mZ2>12",theZ1type,theZ2type); |
| 667 |
|
} else { |
| 668 |
|
ret.status.setStatus(SelectionStatus::FAIL); |
| 669 |
|
return ret; |
| 671 |
|
|
| 672 |
|
if( zzvec.M() > 100 ) { |
| 673 |
|
if(ctrl.debug) cout << "passes mzz > 100, mzz: " << zzvec.M() << endl; |
| 674 |
< |
increment("m4l>100",theZ1type,theZ2type); |
| 674 |
> |
increment(counts,"m4l>100",theZ1type,theZ2type); |
| 675 |
|
} else { |
| 676 |
< |
cout << "NOTE: failed m4l>100 (" << zzvec.M() << "), but saving to ntuple anyway" << endl; |
| 676 |
> |
// cout << "NOTE: failed m4l>100 (" << zzvec.M() << "), but saving to ntuple anyway" << endl; |
| 677 |
|
} |
| 678 |
|
|
| 679 |
|
//*************************************************************** |
| 686 |
|
(ZCandidatesLeptons[best_Z2_index].second.vec); |
| 687 |
|
TLorentzVector theZZ = theZ1 + theZ2; |
| 688 |
|
|
| 689 |
< |
if(ret.fsrPhotons.size() > 0) { |
| 690 |
< |
if(ctrl.debug) |
| 691 |
< |
cout << " fsr photon! pt: " |
| 692 |
< |
<< setw(12) << ret.fsrPhotons[0].vec.Pt() |
| 693 |
< |
<< " eta: " << setw(12) << ret.fsrPhotons[0].vec.Eta() |
| 694 |
< |
<< " phi: " << setw(12) << ret.fsrPhotons[0].vec.Phi() << endl; |
| 695 |
< |
} |
| 696 |
< |
if(ret.fsrPhotons.size() > 1) { |
| 697 |
< |
if(ctrl.debug) |
| 698 |
< |
cout << " fsr photon! pt: " |
| 699 |
< |
<< setw(12) << ret.fsrPhotons[1].vec.Pt() |
| 700 |
< |
<< " eta: " << setw(12) << ret.fsrPhotons[1].vec.Eta() |
| 701 |
< |
<< " phi: " << setw(12) << ret.fsrPhotons[1].vec.Phi() << endl; |
| 702 |
< |
} |
| 689 |
> |
// if(ret.fsrPhotons.size() > 0) { |
| 690 |
> |
// if(ctrl.debug) |
| 691 |
> |
// cout << " fsr photon! pt: " |
| 692 |
> |
// << setw(12) << ret.fsrPhotons[0].vec.Pt() |
| 693 |
> |
// << " eta: " << setw(12) << ret.fsrPhotons[0].vec.Eta() |
| 694 |
> |
// << " phi: " << setw(12) << ret.fsrPhotons[0].vec.Phi() << endl; |
| 695 |
> |
// } |
| 696 |
> |
// if(ret.fsrPhotons.size() > 1) { |
| 697 |
> |
// if(ctrl.debug) |
| 698 |
> |
// cout << " fsr photon! pt: " |
| 699 |
> |
// << setw(12) << ret.fsrPhotons[1].vec.Pt() |
| 700 |
> |
// << " eta: " << setw(12) << ret.fsrPhotons[1].vec.Eta() |
| 701 |
> |
// << " phi: " << setw(12) << ret.fsrPhotons[1].vec.Phi() << endl; |
| 702 |
> |
// } |
| 703 |
|
if( ctrl.debug ) cout << "run: " << info->RunNum() |
| 704 |
|
<< "\tevt: " << info->EvtNum() |
| 705 |
|
<< "\tZ1channel: " << theZ1type |
