babymaker  e95a6a9342d4604277fe7cc6149b6b5b24447d89
bmaker_full.cc
Go to the documentation of this file.
1 
4 // System include files
5 #include <cmath>
6 #include <memory>
7 #include <iostream>
8 #include <stdlib.h>
9 #include <iomanip> // put_time
10 #include <utility> // std::pair
11 
12 // User include files
16 
17 // FW include files
18 #include "FWCore/Framework/interface/Frameworkfwd.h"
19 #include "FWCore/Framework/interface/EDAnalyzer.h"
20 #include "FWCore/Framework/interface/Event.h"
21 #include "FWCore/Framework/interface/MakerMacros.h"
22 #include "FWCore/ParameterSet/interface/ParameterSet.h"
23 
24 // FW physics include files
25 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
26 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h"
27 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
28 #include "DataFormats/TrackReco/interface/Track.h"
29 
30 #include "DataFormats/PatCandidates/interface/PFIsolation.h"
31 
32 // ROOT include files
33 #include "TFile.h"
34 #include "TROOT.h"
35 
36 using namespace std;
37 using namespace utilities;
38 
39 
41 void bmaker_full::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
42  nevents++;
43  isData = iEvent.isRealData();
44  baby.Clear();
45 
47  baby.run() = iEvent.id().run();
48  baby.event() = iEvent.id().event();
49  baby.lumiblock() = iEvent.luminosityBlock();
50  // We are applying the golden JSON with lumisToProcess in bmaker_full_cfg.py
51  if(isData){
52  if (debug) cout<<"INFO: Checking JSON..."<<endl;
53  baby.json4p0() = baby.run() <= 275125;
54  baby.json12p9() = baby.run() <= 276811;
55  } else {
56  baby.json4p0() = baby.json12p9() = true;
57  }
58  baby.type() = event_tools::type(outname.Data());
59 
61  if (debug) cout<<"INFO: Processing trigger info..."<<endl;
62  bool triggerFired;
63  edm::Handle<edm::TriggerResults> triggerBits;
64  iEvent.getByToken(tok_trigResults_hlt_,triggerBits);
65  edm::Handle<pat::PackedTriggerPrescales> triggerPrescales;
66  iEvent.getByToken(tok_patTrig_,triggerPrescales);
67 
68  if(triggerBits.isValid() && triggerPrescales.isValid()){
69  const edm::TriggerNames &names = iEvent.triggerNames(*triggerBits);
70  // Not saving data events that don't have triggers we care about
71  triggerFired = writeTriggers(names, triggerBits, triggerPrescales);
72  }
73  else
74  triggerFired = true;
75  if(!triggerFired && isData) {
76  reportTime(iEvent);
77  return;
78  }
79 
81  if (debug) cout<<"INFO: Writing vertices..."<<endl;
82  edm::Handle<reco::VertexCollection> vtx;
83  iEvent.getByToken(tok_primVertex_, vtx);
84  edm::Handle<std::vector< PileupSummaryInfo > > pu_info;
85  if(!isData) {
86  iEvent.getByToken(tok_addPileup_, pu_info);
87  if(!pu_info.isValid()) iEvent.getByToken(tok_slimAddPileup_, pu_info);
88  }
89  writeVertices(vtx, pu_info);
90 
92  if (debug) cout<<"INFO: Writing leptons..."<<endl;
93  // pfcands, to be used in calculating isolation
94  edm::Handle<pat::PackedCandidateCollection> pfcands;
95  iEvent.getByToken(tok_packPFCands_, pfcands);
96 
97  // Event energy density, to be used in calculating isolation and JECs
98  edm::Handle<double> rhoEventCentral_h;
99  iEvent.getByToken(tok_rhoFastJet_centralNeutral_, rhoEventCentral_h);
100  double rhoEventCentral = static_cast<double>(*rhoEventCentral_h);
101 
102  vCands sig_leps, veto_leps, sig_mus, veto_mus, sig_els, veto_els;
103  vCands all_mus, all_els;
104  edm::Handle<pat::MuonCollection> allmuons;
105  iEvent.getByToken(tok_muons_, allmuons);
106  sig_mus = writeMuons(allmuons, pfcands, vtx, veto_mus, all_mus, rhoEventCentral);
107  edm::Handle<pat::ElectronCollection> allelectrons;
108  iEvent.getByToken(tok_electrons_, allelectrons);
109 
110  sig_els = writeElectrons(allelectrons, pfcands, vtx, veto_els, all_els, rhoEventCentral);
111 
112  writeDiLep(sig_mus, sig_els, veto_mus, veto_els);
113 
114  // Putting muons and electrons together
115  sig_leps = sig_mus;
116  sig_leps.insert(sig_leps.end(), sig_els.begin(), sig_els.end());
117  writeLeptons(sig_leps);
118  // if (baby.nleps()<1) return;
119  veto_leps = veto_mus;
120  veto_leps.insert(veto_leps.end(), veto_els.begin(), veto_els.end());
121 
122  // bool keep_event = false;
123  // if (baby.nleps()==2&&((baby.elel_m()>80&&baby.elel_m()<100)||(baby.mumu_m()>80&&baby.mumu_m()<100)))
124  // keep_event = true;
125 
127  if (debug) cout<<"INFO: Writing photons..."<<endl;
128  vCands photons;
129  edm::Handle<double> rhoEvent_h;
130  iEvent.getByToken(tok_rhoFastJet_all_, rhoEvent_h);
131  edm::Handle<reco::BeamSpot> beamspot;
132  iEvent.getByToken(tok_offBeamSpot_, beamspot);
133  edm::Handle<pat::PhotonCollection> allphotons;
134  iEvent.getByToken(tok_photons_, allphotons);
135  edm::Handle<vector<reco::Conversion> > conversions;
136  iEvent.getByToken(tok_reducedEgamma_conver_, conversions);
137  photons = writePhotons(allphotons, allelectrons, conversions, beamspot, *rhoEvent_h);
138 
140  if (debug) cout<<"INFO: Applying JECs..."<<endl;
141  edm::Handle<pat::JetCollection> alljets;
142  iEvent.getByToken(tok_jets_, alljets);
143  edm::Handle<edm::View<reco::GenJet> > genjets;
144  iEvent.getByToken(tok_genJets_, genjets) ;
145  jetTool->getJetCorrections(genjets, alljets, *rhoEvent_h);
146 
148  if (debug) cout<<"INFO: Writing MET..."<<endl;
149  edm::Handle<pat::METCollection> mets;
150  edm::Handle<pat::METCollection> mets_nohf;
151  edm::Handle<pat::METCollection> mets_puppi;
152 
153  iEvent.getByToken(tok_met_noHF_, mets_nohf);
154  if (!isData) {
155  iEvent.getByToken(tok_met_, mets); // using MuEGClean for default, for now
156  } else {
157  iEvent.getByToken(tok_met_Puppi_, mets_puppi);
158  iEvent.getByToken(tok_met_, mets); //The collection called "slimmedMETs" is corrected for muons but not EG
159  }
160 
161  //Saving these lines here in case we decide to switch to a different default
162  // edm::Handle<pat::METCollection> mets_muegclean;
163 
164 
165  writeMET(mets, mets_nohf, mets_puppi);
166 
168  if (debug) cout<<"INFO: Calculating track veto..."<<endl;
169  vCands tks,ra4tks;
170  if (eventTool->hasGoodPV(vtx)){
171  // RA2/b track veto
172  tks = lepTool->getIsoTracks(pfcands, baby.met(), baby.met_phi());
173  baby.ntks() = tks.size();
174 
175  // RA4 track veto
176  writeTks(pfcands,vtx,rhoEventCentral);
177  }
178 
180  if (debug) cout<<"INFO: Writing jets..."<<endl;
181  vector<LVector> all_baby_jets;
182  vector<unsigned> all_baby_jets_idx;
183  vector<double> jetsMuonEnergyFrac;
184  vector<vector<LVector> > sys_jets;
185  all_baby_jets = writeJets(alljets, all_baby_jets_idx, genjets, sig_leps, veto_leps, photons, tks, sys_jets, jetsMuonEnergyFrac);
186  // if (baby.njets()>=4) keep_event = true;
187  // if (!keep_event && isData) return;
188  if (addBTagWeights) writeBTagWeights(alljets, all_baby_jets, all_baby_jets_idx);
189  writeHiggVars(all_baby_jets, baby.jets_csv(), baby.jets_h1(), baby.jets_h2(),
190  baby.jets_islep(), baby.nbl(), baby.nbm(), baby.nbt(),
191  baby.hig_am(), baby.hig_dm(), baby.hig_drmax(), baby.hig_bin(), baby.mct());
192  writeHiggVars(all_baby_jets, baby.jets_csvd(), baby.jets_h1d(), baby.jets_h2d(),
193  baby.jets_islep(), baby.nbdl(), baby.nbdm(), baby.nbdt(),
194  baby.higd_am(), baby.higd_dm(), baby.higd_drmax(), baby.higd_bin(), baby.mctd());
195  writeBBVars(all_baby_jets, sig_leps);
196  writeFatJets();
197 /* DAK8
198  // DeepAk8
199  fatjetNN_->readEvent(iEvent, iSetup);
200  decorrNN_->readEvent(iEvent, iSetup);
201  edm::Handle<edm::View<pat::Jet>> ak8jets;
202  iEvent.getByToken(tok_deepJetToken_, ak8jets);
203  writeAk8Jets(ak8jets);
204 */
205  if (doSystematics) {
206  for (unsigned isys(0); isys<kSysLast; isys++){
207  bool cluster_leps = false;
208  baby.sys_mj14_nolep().push_back(jetTool->getSysMJ(1.4, sys_jets[isys], baby.jets_islep(), jetTool->JetPtCut, cluster_leps));
209  cluster_leps = true;
210  baby.sys_mj12().push_back(jetTool->getSysMJ(1.2, sys_jets[isys], baby.jets_islep(), jetTool->JetPtCut, cluster_leps));
211  baby.sys_mj14().push_back(jetTool->getSysMJ(1.4, sys_jets[isys], baby.jets_islep(), jetTool->JetPtCut, cluster_leps));
212 
213  int dummy_int(0); vector<bool> dummy_v; float dummy_flt(0);
214  baby.sys_higd_am().resize(kSysLast, 0);
215  baby.sys_higd_dm().resize(kSysLast, 0);
216  baby.sys_higd_drmax().resize(kSysLast, 0);
217  writeHiggVars(sys_jets[isys], baby.jets_csvd(), dummy_v, dummy_v,
218  baby.jets_islep(), dummy_int, dummy_int, dummy_int,
219  baby.sys_higd_am()[isys], baby.sys_higd_dm()[isys],
220  baby.sys_higd_drmax()[isys], dummy_int, dummy_flt, doSystematics);
221  }
222  }
223 
225  // It requires previous storing of MET
226  if (debug) cout<<"INFO: Calculating mT..."<<endl;
227  if(sig_leps.size()>0){
228  float wx = baby.met()*cos(baby.met_phi()) + sig_leps[0]->px();
229  float wy = baby.met()*sin(baby.met_phi()) + sig_leps[0]->py();
230  float wphi = atan2(wy, wx);
231  baby.dphi_wlep() = deltaPhi(wphi, sig_leps[0]->phi());
232 
233  baby.mt() = getMT(baby.met(), baby.met_phi(), sig_leps[0]->pt(), sig_leps[0]->phi());
234  baby.mt_nohf() = getMT(baby.met_nohf(), baby.met_nohf_phi(), sig_leps[0]->pt(), sig_leps[0]->phi());
235  }
236  // calculating met with systematics and corresponding mT
237  if (doSystematics) {
238  baby.sys_met().resize(kSysLast,-999.);
239  baby.sys_mt().resize(kSysLast,-999.);
240  for (unsigned isys(0); isys<kSysLast; isys++) {
241  float sys_met_phi(0.);
242  jetTool->getMETWithJEC(mets, baby.sys_met()[isys], sys_met_phi, isys);
243  if (sig_leps.size()>0)
244  baby.sys_mt()[isys] = getMT(baby.sys_met()[isys], sys_met_phi, sig_leps[0]->pt(), sig_leps[0]->phi());
245  }
246  }
247 
249  if (debug) cout<<"INFO: Writing filters..."<<endl;
250  edm::Handle<edm::TriggerResults> filterBits;
251  iEvent.getByToken(tok_trigResults_pat_,filterBits);
252  //Giovanni filters in reminiaod only available in "PAT"
253  //Try reco if fails..
254  if(!filterBits.isValid() && isData)
255  iEvent.getByToken(tok_trigResults_reco_,filterBits);
256  if(filterBits.isValid()){
257  const edm::TriggerNames &fnames = iEvent.triggerNames(*filterBits);
258  //this method uses baby.pass_jets(), so call only after writeJets()!!
259  writeFilters(fnames, filterBits, vtx, jetsMuonEnergyFrac);
260  }
261 
262 
264  if (debug) cout<<"INFO: Writing HLT objects..."<<endl;
265  edm::Handle<pat::TriggerObjectStandAloneCollection> triggerObjects;
266  iEvent.getByToken(tok_selectedPatTrig_,triggerObjects);
267  // Requires having called writeMuons and writeElectrons for truth-matching
268  if(triggerBits.isValid() && triggerPrescales.isValid()){
269  const edm::TriggerNames &names = iEvent.triggerNames(*triggerBits);
270  writeHLTObjects(names, triggerObjects, all_mus, all_els, photons, iEvent);
271  }
272 
274  if (!isData) {
275  if (debug) cout<<"INFO: Writing MC particles..."<<endl;
276  edm::Handle<reco::GenParticleCollection> genParticles;
277  iEvent.getByToken(tok_pruneGenPart_, genParticles);
278  writeMC(genParticles, all_mus, all_els, photons);
279  writeIFSR(genParticles, all_baby_jets);
280  }
281 
283  if (debug) cout<<"INFO: Calculating MET rebalancing..."<<endl;
284  baby.jetmismeas() = false;
285  if(doMetRebalancing && sig_leps.size()==1) {
286  double temp_met = baby.met();
287  double temp_met_phi = baby.met_phi();
288  rebalancedMET(temp_met, temp_met_phi);
289  baby.met_rebal() = temp_met;
290  baby.mt_rebal() = getMT(temp_met, temp_met_phi, sig_leps[0]->pt(), sig_leps[0]->phi());
291  if(baby.met_rebal()/baby.met()<0.2 && baby.mt_rebal()<150) baby.jetmismeas()=true;
292  } else {
293  // use default values for events that do not have exactly one lepton
294  baby.mt_rebal() = baby.mt();
295  baby.met_rebal() = baby.met();
296  }
297 
299  if (debug) cout<<"INFO: Retrieving hard scatter info..."<<endl;
300  edm::Handle<LHEEventProduct> lhe_info;
301  baby.stitch() = baby.stitch_ht() = baby.stitch_met() = true;
302  if (outname.Contains("SMS-") && outname.Contains("PUSpring16Fast")) {
303  baby.mgluino() = mprod_;
304  baby.mlsp() = mlsp_;
305  } else if (!isData) {
306  iEvent.getByToken(tok_extLHEProducer_, lhe_info);
307  if(!lhe_info.isValid()) iEvent.getByToken(tok_source_, lhe_info);
308  if(lhe_info.isValid()) writeGenInfo(lhe_info);
309  if((outname.Contains("TTJets_SingleLeptFromT_Tune") ||
310  outname.Contains("TTJets_SingleLeptFromTbar_Tune") ||
311  outname.Contains("TTJets_DiLept_Tune"))
312  && outname.Contains("madgraphMLM")) {
313  if (baby.ht_isr_me()>600) {
314  baby.stitch() = false;
315  baby.stitch_ht() = false;
316  }
317  if (baby.met_tru()>150) {
318  baby.stitch_met() = false;
319  }
320  }
321  if((outname.Contains("DYJetsToLL_M-50_TuneCUETP8M") && baby.ht_isr_me()>70)
322  || (outname.Contains("WJetsToLNu_TuneCUETP8M1") && baby.ht_isr_me()>70)){
323  baby.stitch() = false;
324  baby.stitch_ht() = false;
325  baby.stitch_met() = false;
326  }
327  if(outname.Contains("TTJets_Tune") && outname.Contains("madgraphMLM")){
328  if(baby.ntruleps()!=0) {
329  baby.stitch()=false;
330  baby.stitch_met()=false;
331  }
332  if(baby.ht_isr_me()>600) baby.stitch_ht()=false;
333  }
334  if(outname.Contains("DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX")){
335  if (baby.ptll_me()>50) {
336  baby.stitch()=false;
337  }
338  }
339  } // if it is not data
340 
342  // w_btag calculated in writeJets
343  // w_toppt calculated in writeMC
344  edm::Handle<GenEventInfoProduct> gen_event_info;
345  iEvent.getByToken(tok_generator_, gen_event_info);
346  //assuming baby.nisr() is filled
347  writeWeights(sig_leps, gen_event_info, lhe_info);
348 
350  baby.Fill();
351 
352  reportTime(iEvent);
353 
354 }
355 
356 
357 /*
358  ______ _ _ _ _
359  | ___ \ | | (_) | (_)
360  | |_/ /_ __ __ _ _ __ ___| |__ __ ___ __ _| |_ _ _ __ __ _
361  | ___ \ '__/ _` | '_ \ / __| '_ \ \ \ /\ / / '__| | __| | '_ \ / _` |
362  | |_/ / | | (_| | | | | (__| | | | \ V V /| | | | |_| | | | | (_| |
363  \____/|_| \__,_|_| |_|\___|_| |_| \_/\_/ |_| |_|\__|_|_| |_|\__, |
364  __/ |
365  |___/
366 */
367 
368 // Requires having called jetTool->getJetCorrections(alljets, rhoEvent_) beforehand
369 void bmaker_full::writeMET(edm::Handle<pat::METCollection> mets, edm::Handle<pat::METCollection> mets_nohf, edm::Handle<pat::METCollection> mets_puppi){
370  jetTool->getMETWithJEC(mets, baby.met(), baby.met_phi(), kSysLast);
371  jetTool->getMETRaw(mets, baby.met_raw(), baby.met_raw_phi());
372 
373  //Avoid overwriting systematics, don't call getMETwithJEC
374  // jetTool->getMETWithJEC(mets_uncorr, baby.met_uncorr(), baby.met_phi_uncorr(), kSysLast);
375  // jetTool->getMETWithJEC(mets_egclean, baby.met_egclean(), baby.met_phi_egclean(), kSysLast);
376  // jetTool->getMETWithJEC(mets_muclean, baby.met_muclean(), baby.met_phi_muclean(), kSysLast);
377 
378  if (isData) {
379  baby.met_puppi() = mets_puppi->at(0).pt();
380  baby.met_phi_puppi() = mets_puppi->at(0).phi();
381  }
382 
383  baby.met_mini() = mets->at(0).pt();
384  baby.met_mini_phi() = mets->at(0).phi();
385  baby.met_calo() = mets->at(0).caloMETPt();
386  baby.met_calo_phi() = mets->at(0).caloMETPhi();
387  if(mets_nohf.isValid()){
388  baby.met_nohf() = mets_nohf->at(0).pt();
389  baby.met_nohf_phi() = mets_nohf->at(0).phi();
390  }
391  if(!isData){
392  baby.met_tru() = mets->at(0).genMET()->pt();
393  baby.met_tru_phi() = mets->at(0).genMET()->phi();
394  }
395  if(isnan(baby.met())) {
396  cout<<"MET is NaN. Perhaps you are running on old miniAOD with a new release. Setting MET to zero."<<endl;
397  baby.met() = 0;
398  }
399 }
400 
401 // Requires having called jetTool->getJetCorrections(alljets, rhoEvent_) beforehand
402 vector<LVector> bmaker_full::writeJets(edm::Handle<pat::JetCollection> alljets,
403  vector<unsigned> &all_baby_jets_idx,
404  edm::Handle<edm::View <reco::GenJet> > genjets,
405  vCands &sig_leps, vCands &veto_leps, vCands &photons, vCands &tks,
406  vector<vector<LVector> > &sys_jets,
407  vector<double> &jetsMuonEnergyFrac){
408  vector<LVector> all_baby_jets;
409  vCands jets_ra2;
410  LVector jetsys_p4, jetsys_nob_p4, jetsys_nobd_p4;
411  baby.njets() = 0; baby.nbl() = 0; baby.nbm() = 0; baby.nbt() = 0;
412  baby.nbdl() = 0; baby.nbdm() = 0; baby.nbdt() = 0;
413  baby.nbdfl() = 0; baby.nbdfm() = 0; baby.nbdft() = 0;
414  baby.ht() = 0.; baby.st() = 0.; baby.ht_hlt() = 0.;
415  baby.njets_ra2() = 0; baby.njets_clean() = 0; baby.nbm_ra2() = 0; baby.ht_ra2() = 0.; baby.ht_clean() = 0.;
416  baby.pass_jets() = true; baby.pass_jets_nohf() = true; baby.pass_jets_tight() = true;
417  baby.pass_jets_ra2() = true; baby.pass_jets_tight_ra2() = true; baby.pass_fsjets()=true;
418  if (doSystematics) {
419  baby.sys_njets().resize(kSysLast, 0); baby.sys_nbm().resize(kSysLast, 0);
420  baby.sys_nbdl().resize(kSysLast, 0); baby.sys_nbdm().resize(kSysLast, 0); baby.sys_nbdt().resize(kSysLast, 0);
421  baby.sys_pass().resize(kSysLast, true); baby.sys_ht().resize(kSysLast, 0.);
422  baby.sys_st().resize(kSysLast, 0);
423  sys_jets.resize(kSysLast, vector<LVector>());
424  }
425  float mht_px(0.), mht_py(0.), mht_clean_px(0.), mht_clean_py(0.);
426  for (size_t ijet(0); ijet < alljets->size(); ijet++) {
427  const pat::Jet &jet = (*alljets)[ijet];
428  if(fabs(jet.eta()) > 5) continue;
429 
430  LVector jetp4(jetTool->corrJet[ijet]);
431  float csv(jet.bDiscriminator("pfCombinedInclusiveSecondaryVertexV2BJetTags"));
432  float csvd(-999.), csvdf(-999.);
433  TString cmssw_rel = getenv("CMSSW_BASE");
434  if (cmssw_rel.Contains("CMSSW_8"))
435  csvd = jet.bDiscriminator("deepFlavourJetTags:probb")+jet.bDiscriminator("deepFlavourJetTags:probbb");
436  else {
437  csvd = jet.bDiscriminator("pfDeepCSVJetTags:probb")+jet.bDiscriminator("pfDeepCSVJetTags:probbb");
438  csvdf = jet.bDiscriminator("pfDeepFlavourJetTags:probb")+jet.bDiscriminator("pfDeepFlavourJetTags:problepb")
439  +jet.bDiscriminator("pfDeepFlavourJetTags:probbb");
440  }
441 
442  bool isLep = jetTool->leptonInJet(jet, sig_leps);
443  bool looseID = jetTool->idJet(jet, jetTool->kLoose);
444  bool tightID = jetTool->idJet(jet, jetTool->kTight);
445  bool goodPtEta = jetp4.pt() > jetTool->JetPtCut && fabs(jet.eta()) <= jetTool->JetEtaCut;
446  if(isFastSim){
447  if(jetp4.pt() > 20. && fabs(jet.eta()) < 2.5 && !jetTool->matchesGenJet(jet,genjets) && jet.chargedHadronEnergyFraction() < 0.1) baby.pass_fsjets()=false;
448  }
449 
450  // RA4 Jet Quality filters
451  //--------------------------------
452  if(jetp4.pt() > jetTool->JetPtCut && !isLep) {
453  if(goodPtEta && !looseID) baby.pass_jets_nohf() = false;
454  if(!looseID) baby.pass_jets() = false;
455  if(!tightID) baby.pass_jets_tight() = false;
456  }
457 
458  // RA4 Jet Variables
459  //----------------------------
460  if((looseID && goodPtEta) || isLep) {
461  all_baby_jets.push_back(jetp4);
462  all_baby_jets_idx.push_back(ijet);
463  baby.jets_pt().push_back(jetp4.pt());
464  baby.jets_eta().push_back(jet.eta());
465  baby.jets_phi().push_back(jet.phi());
466  baby.jets_m().push_back(jetp4.mass());
467  baby.jets_islep().push_back(isLep);
468  jetsMuonEnergyFrac.push_back(jet.muonEnergyFraction());
469  baby.jets_gen_pt().push_back(isData ? 0. : jetTool->genJetPt[ijet]);
470  baby.jets_pflavor().push_back(jet.partonFlavour());
471  baby.jets_hflavor().push_back(jet.hadronFlavour());
472  baby.jets_ntrub().push_back(0); //filled in writeMC
473  baby.jets_ntruc().push_back(0); //filled in writeMC
474  baby.jets_gs_index().push_back(-1); //filled in writeMC
475  baby.jets_csv().push_back(csv);
476  baby.jets_csvd().push_back(csvd);
477  baby.jets_csvdf().push_back(csvdf);
478 
479  if(!isLep && goodPtEta){
480  jetsys_p4 += jet.p4();
481  baby.njets()++;
482  baby.ht() += jetp4.pt();
483  baby.st() += jetp4.pt();
484  if(csv > jetTool->CSVLoose) baby.nbl()++;
485  if(csv > jetTool->CSVMedium) baby.nbm()++;
486  else jetsys_nob_p4 += jet.p4();
487  if(csv > jetTool->CSVTight) baby.nbt()++;
488 
489  if(csvd > jetTool->DeepCSVLoose) baby.nbdl()++;
490  if(csvd > jetTool->DeepCSVMedium) baby.nbdm()++;
491  else jetsys_nobd_p4 += jet.p4();
492  if(csvd > jetTool->DeepCSVTight) baby.nbdt()++;
493 
494  if(csvdf > jetTool->DeepFlavourLoose) baby.nbdfl()++;
495  if(csvdf > jetTool->DeepFlavourMedium) baby.nbdfm()++;
496  if(csvdf > jetTool->DeepFlavourTight) baby.nbdft()++;
497  }
498  }
499 
500  // HLT HT definition
501  //----------------------------
502  if(jetp4.pt() > jetTool->JetHLTPtCut && fabs(jet.eta()) <= jetTool->JetHLTEtaCut) baby.ht_hlt() += jetp4.pt();
503 
504  // RA2 Jet Variables
505  //----------------------------
506  bool isLep_ra2 = jetTool->jetMatched(jet, veto_leps); // Uses RA2/b's loose matching, dpt/pt < 100%, dR < 0.4
507  bool isPhoton = jetTool->jetMatched(jet, photons); // Uses RA2/b's loose matching, dpt/pt < 100%, dR < 0.4
508  bool isIsoTrack = jetTool->jetMatched(jet, tks); // Uses RA2/b's loose matching, dpt/pt < 100%, dR < 0.4
509  bool applyId_ra2 = !isLep_ra2 && !isPhoton && !isIsoTrack; // Only check ID if jet not matched
510  bool goodJet_ra2 = (looseID || !applyId_ra2);
511  bool tightJet_ra2 = (tightID || !applyId_ra2);
512 
513  if(jetp4.pt() > jetTool->JetPtCut) { // Check jet ID on 30 GeV jets
514  if(!goodJet_ra2) baby.pass_jets_ra2() = false;
515  if(!tightJet_ra2) baby.pass_jets_tight_ra2() = false;
516  }
517 
518  if(goodPtEta && goodJet_ra2) {
519  baby.njets_ra2()++;
520  baby.ht_ra2() += jetp4.pt();
521  jets_ra2.push_back(dynamic_cast<const reco::Candidate *>(&jet));
522  if(csv > jetTool->CSVMedium) baby.nbm_ra2()++;
523  if(!isLep_ra2 && !isPhoton) {
524  baby.njets_clean()++;
525  baby.ht_clean() += jetp4.pt();
526  }
527  }
528  if(goodJet_ra2 && jetp4.pt() > jetTool->JetPtCut && fabs(jet.eta()) <= jetTool->JetMHTEtaCut){
529  mht_px -= jet.px();
530  mht_py -= jet.py();
531  if(!isLep_ra2 && !isPhoton) {
532  mht_clean_px -= jet.px();
533  mht_clean_py -= jet.py();
534  }
535  }
536 
537  // Systematic variations
538  //----------------------------
539  if (doSystematics){
540  for (unsigned isys(0); isys<kSysLast; isys++){
541  jetp4 = jetTool->corrJet[ijet];
542  if (isys == kSysJER) jetp4 *= (1 + jetTool->jerUnc[ijet]);
543  else if (isys == kSysJECUp) jetp4 *= (1 + jetTool->jecUnc[ijet]);
544  else if (isys == kSysJECDn) jetp4 *= (1 - jetTool->jecUnc[ijet]);
545  //for now store pass_jets in the baby.sys_pass
546  //variable will be updated with rest of filters in writeFilters()
547  if(jetp4.pt() > jetTool->JetPtCut && !isLep && !looseID) baby.sys_pass()[isys] = false;
548 
549  //cutting on unaltered goodPtEtaLoose since indices must match baby.jets_islep()
550  if((looseID && goodPtEta) || isLep) sys_jets[isys].push_back(jetp4);
551  //now cutting on altered variable for calculation of njets, ht, st
552  if(looseID && jetp4.pt() > jetTool->JetPtCut && fabs(jet.eta()) <= jetTool->JetEtaCut) {
553  if(!isLep){
554  baby.sys_njets()[isys]++;
555  baby.sys_ht()[isys] += jetp4.pt();
556  baby.sys_st()[isys] += jetp4.pt();
557  if(csv > jetTool->CSVMedium) baby.sys_nbm()[isys]++;
558  if(csvd > jetTool->DeepCSVLoose) baby.sys_nbdl()[isys]++;
559  if(csvd > jetTool->DeepCSVMedium) baby.sys_nbdm()[isys]++;
560  if(csvd > jetTool->DeepCSVTight) baby.sys_nbdt()[isys]++;
561  }
562  }
563  } // loop over systematics
564  } // do systematics
565  } // Loop over jets
566 
567  for (size_t ilep(0); ilep<sig_leps.size(); ilep++){
568  baby.st() += sig_leps[ilep]->pt();
569  if (doSystematics){
570  for (unsigned isys(0); isys<kSysLast; isys++){
571  baby.sys_st()[isys] += sig_leps[ilep]->pt();
572  }
573  }
574  }
575  if(!isData) baby.ht_tru() = jetTool->trueHT(genjets);
576 
577  baby.mht() = hypot(mht_px, mht_py);
578  baby.mht_phi() = atan2(mht_py, mht_px);
579  baby.mht_clean() = hypot(mht_clean_px, mht_clean_py);
580  baby.mht_clean_phi() = atan2(mht_clean_py, mht_clean_px);
581  baby.low_dphi() = jetTool->isLowDphi(jets_ra2, baby.mht_phi(), baby.dphi1(), baby.dphi2(), baby.dphi3(), baby.dphi4());
582 
583  // ISR system for Z->ll and tt->llbb configurations
584  baby.jetsys_pt() = jetsys_p4.pt();
585  baby.jetsys_eta() = jetsys_p4.eta();
586  baby.jetsys_phi() = jetsys_p4.phi();
587  baby.jetsys_m() = jetsys_p4.mass();
588  baby.jetsys_nob_pt() = jetsys_nob_p4.pt();
589  baby.jetsys_nob_eta() = jetsys_nob_p4.eta();
590  baby.jetsys_nob_phi() = jetsys_nob_p4.phi();
591  baby.jetsys_nob_m() = jetsys_nob_p4.mass();
592  baby.jetsys_nobd_pt() = jetsys_nobd_p4.pt();
593  baby.jetsys_nobd_eta() = jetsys_nobd_p4.eta();
594  baby.jetsys_nobd_phi() = jetsys_nobd_p4.phi();
595  baby.jetsys_nobd_m() = jetsys_nobd_p4.mass();
596 
597  if(isFastSim) baby.pass_fsmet() = eventTool->passFSMET(alljets, genjets);
598  else baby.pass_fsmet() = true;
599 
600  return all_baby_jets;
601 } // writeJets
602 
603 void bmaker_full::writeBTagWeights(edm::Handle<pat::JetCollection> alljets,
604  std::vector<reco::Candidate::LorentzVector> &all_baby_jets,
605  std::vector<unsigned> &all_baby_jet_idx){
606 
607  baby.w_btag() = baby.w_btag_loose() = baby.w_btag_tight() = baby.w_bhig() = 1.;
608  baby.w_btag_deep() = baby.w_btag_loose_deep() = baby.w_btag_tight_deep() = baby.w_bhig_deep() = 1.;
609 
610  baby.w_btag_proc() = baby.w_btag_loose_proc() = baby.w_btag_tight_proc() = baby.w_bhig_proc() = 1.;
611  baby.w_btag_deep_proc() = baby.w_btag_loose_deep_proc() = baby.w_btag_tight_deep_proc() = baby.w_bhig_deep_proc() = 1.;
612  if (doSystematics){
613  baby.sys_bctag().resize(2, 1.); baby.sys_udsgtag().resize(2, 1.);
614  baby.sys_bctag_deep().resize(2, 1.); baby.sys_udsgtag_deep().resize(2, 1.);
615  baby.sys_bctag_loose().resize(2, 1.); baby.sys_udsgtag_loose().resize(2, 1.);
616  baby.sys_bctag_loose_deep().resize(2, 1.); baby.sys_udsgtag_loose_deep().resize(2, 1.);
617  baby.sys_bctag_tight().resize(2, 1.); baby.sys_udsgtag_tight().resize(2, 1.);
618  baby.sys_bctag_tight_deep().resize(2, 1.); baby.sys_udsgtag_tight_deep().resize(2, 1.);
619  baby.sys_bchig().resize(2, 1.); baby.sys_udsghig().resize(2, 1.);
620  baby.sys_bchig_deep().resize(2, 1.); baby.sys_udsghig_deep().resize(2, 1.);
621 
622  baby.sys_bctag_proc().resize(2, 1.); baby.sys_udsgtag_proc().resize(2, 1.);
623  baby.sys_bctag_deep_proc().resize(2, 1.); baby.sys_udsgtag_deep_proc().resize(2, 1.);
624  baby.sys_bctag_loose_proc().resize(2, 1.); baby.sys_udsgtag_loose_proc().resize(2, 1.);
625  baby.sys_bctag_loose_deep_proc().resize(2, 1.); baby.sys_udsgtag_loose_deep_proc().resize(2, 1.);
626  baby.sys_bctag_tight_proc().resize(2, 1.); baby.sys_udsgtag_tight_proc().resize(2, 1.);
627  baby.sys_bctag_tight_deep_proc().resize(2, 1.); baby.sys_udsgtag_tight_deep_proc().resize(2, 1.);
628  baby.sys_bchig_proc().resize(2, 1.); baby.sys_udsghig_proc().resize(2, 1.);
629  baby.sys_bchig_deep_proc().resize(2, 1.); baby.sys_udsghig_deep_proc().resize(2, 1.);
630  if (isFastSim) {
631  baby.sys_fs_bctag().resize(2, 1.); baby.sys_fs_udsgtag().resize(2, 1.);
632  baby.sys_fs_bctag_deep().resize(2, 1.); baby.sys_fs_udsgtag_deep().resize(2, 1.);
633  baby.sys_fs_bchig().resize(2, 1.); baby.sys_fs_udsghig().resize(2, 1.);
634  baby.sys_fs_bchig_deep().resize(2, 1.); baby.sys_fs_udsghig_deep().resize(2, 1.);
635  }
636  }
637 
638  for (size_t ijet(0); ijet < all_baby_jet_idx.size(); ijet++) {
639  const pat::Jet &jet = (*alljets)[all_baby_jet_idx[ijet]];
640  LVector &jetp4 = all_baby_jets[ijet];
641  const vector<BTagEntry::OperatingPoint> all_ops = {BTagEntry::OP_LOOSE, BTagEntry::OP_MEDIUM, BTagEntry::OP_TIGHT};
642  if(!baby.jets_islep()[ijet]){
643  const string ctr("central"), vup("up"), vdn("down");
644  //central weight for fastsim taken into account together with the fullsim inside jetTool->jetBTagWeight()
645  baby.w_btag() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, false);
646  baby.w_btag_loose() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, ctr, false);
647  baby.w_btag_tight() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, ctr, false);
648  baby.w_bhig() *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, false);
649  // with deepCSV
650  baby.w_btag_deep() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, true);
651  baby.w_btag_loose_deep() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, ctr, true);
652  baby.w_btag_tight_deep() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, ctr, true);
653  baby.w_bhig_deep() *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, true);
654 
655  // MC Efficiencies taken by process
656  baby.w_btag_proc() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, false, true);
657  baby.w_btag_loose_proc() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, ctr, false, true);
658  baby.w_btag_tight_proc() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, ctr, false, true);
659  baby.w_bhig_proc() *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, false, true);
660  // with deepCSV
661  baby.w_btag_deep_proc() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, true, true);
662  baby.w_btag_loose_deep_proc() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, ctr, true, true);
663  baby.w_btag_tight_deep_proc() *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, ctr, true, true);
664  baby.w_bhig_deep_proc() *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, true, true);
665  if (doSystematics){
666  // now, vary only the full sim scale factor, regardless of whether we run on Fast or Full sim
667  // this is necessary since uncertainties for FastSim and FullSim uncorrelated
668  baby.sys_bctag()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vup, ctr, false);
669  baby.sys_bctag()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vdn, ctr, false);
670  baby.sys_udsgtag()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vup, false);
671  baby.sys_udsgtag()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vdn, false);
672  baby.sys_bctag_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vup, ctr, true);
673  baby.sys_bctag_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vdn, ctr, true);
674  baby.sys_udsgtag_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vup, true);
675  baby.sys_udsgtag_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vdn, true);
676  baby.sys_bctag_loose()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vup, ctr, false);
677  baby.sys_bctag_loose()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vdn, ctr, false);
678  baby.sys_udsgtag_loose()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vup, false);
679  baby.sys_udsgtag_loose()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vdn, false);
680  baby.sys_bctag_loose_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vup, ctr, true);
681  baby.sys_bctag_loose_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vdn, ctr, true);
682  baby.sys_udsgtag_loose_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vup, true);
683  baby.sys_udsgtag_loose_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vdn, true);
684  baby.sys_bctag_tight()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vup, ctr, false);
685  baby.sys_bctag_tight()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vdn, ctr, false);
686  baby.sys_udsgtag_tight()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vup, false);
687  baby.sys_udsgtag_tight()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vdn, false);
688  baby.sys_bctag_tight_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vup, ctr, true);
689  baby.sys_bctag_tight_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vdn, ctr, true);
690  baby.sys_udsgtag_tight_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vup, true);
691  baby.sys_udsgtag_tight_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vdn, true);
692  baby.sys_bchig()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vup, ctr, false);
693  baby.sys_bchig()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vdn, ctr, false);
694  baby.sys_udsghig()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vup, false);
695  baby.sys_udsghig()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vdn, false);
696  baby.sys_bchig_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vup, ctr, true);
697  baby.sys_bchig_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vdn, ctr, true);
698  baby.sys_udsghig_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vup, true);
699  baby.sys_udsghig_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vdn, true);
700 
701  baby.sys_bctag_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vup, ctr, false, true);
702  baby.sys_bctag_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vdn, ctr, false, true);
703  baby.sys_udsgtag_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vup, false, true);
704  baby.sys_udsgtag_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vdn, false, true);
705  baby.sys_bctag_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vup, ctr, true, true);
706  baby.sys_bctag_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, vdn, ctr, true, true);
707  baby.sys_udsgtag_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vup, true, true);
708  baby.sys_udsgtag_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, vdn, true, true);
709  baby.sys_bctag_loose_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vup, ctr, false, true);
710  baby.sys_bctag_loose_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vdn, ctr, false, true);
711  baby.sys_udsgtag_loose_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vup, false, true);
712  baby.sys_udsgtag_loose_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vdn, false, true);
713  baby.sys_bctag_loose_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vup, ctr, true, true);
714  baby.sys_bctag_loose_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, vdn, ctr, true, true);
715  baby.sys_udsgtag_loose_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vup, true, true);
716  baby.sys_udsgtag_loose_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_LOOSE, ctr, vdn, true, true);
717  baby.sys_bctag_tight_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vup, ctr, false, true);
718  baby.sys_bctag_tight_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vdn, ctr, false, true);
719  baby.sys_udsgtag_tight_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vup, false, true);
720  baby.sys_udsgtag_tight_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vdn, false, true);
721  baby.sys_bctag_tight_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vup, ctr, true, true);
722  baby.sys_bctag_tight_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, vdn, ctr, true, true);
723  baby.sys_udsgtag_tight_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vup, true, true);
724  baby.sys_udsgtag_tight_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_TIGHT, ctr, vdn, true, true);
725  baby.sys_bchig_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vup, ctr, false, true);
726  baby.sys_bchig_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vdn, ctr, false, true);
727  baby.sys_udsghig_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vup, false, true);
728  baby.sys_udsghig_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vdn, false, true);
729  baby.sys_bchig_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vup, ctr, true, true);
730  baby.sys_bchig_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, vdn, ctr, true, true);
731  baby.sys_udsghig_deep_proc()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vup, true, true);
732  baby.sys_udsghig_deep_proc()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, vdn, true, true);
733  if (isFastSim) {
734  // now we vary only the FastSim SF
735  baby.sys_fs_bctag()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, vup, ctr, false);
736  baby.sys_fs_bctag()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, vdn, ctr, false);
737  baby.sys_fs_udsgtag()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, ctr, vup, false);
738  baby.sys_fs_udsgtag()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, ctr, vdn, false);
739  baby.sys_fs_bctag_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, vup, ctr, true);
740  baby.sys_fs_bctag_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, vdn, ctr, true);
741  baby.sys_fs_udsgtag_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, ctr, vup, true);
742  baby.sys_fs_udsgtag_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, BTagEntry::OP_MEDIUM, ctr, ctr, ctr, vdn, true);
743  baby.sys_fs_bchig()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, vup, ctr, false);
744  baby.sys_fs_bchig()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, vdn, ctr, false);
745  baby.sys_fs_udsghig()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, ctr, vup, false);
746  baby.sys_fs_udsghig()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, ctr, vdn, false);
747  baby.sys_fs_bchig_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, vup, ctr, true);
748  baby.sys_fs_bchig_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, vdn, ctr, true);
749  baby.sys_fs_udsghig_deep()[0] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, ctr, vup, true);
750  baby.sys_fs_udsghig_deep()[1] *= jetTool->jetBTagWeight(jet, jetp4, all_ops, ctr, ctr, ctr, vdn, true);
751  }
752  }
753  }
754  }
755 }
756 
757 void bmaker_full::writeHiggVars(vector<LVector> &baby_jets_p4, vector<float> &baby_jets_csv,
758  vector<bool> &baby_jets_h1, vector<bool> &baby_jets_h2,
759  vector<bool> &baby_jets_islep, int &baby_nbl, int &baby_nbm, int &baby_nbt,
760  float &baby_hig_am, float &baby_hig_dm, float &baby_hig_drmax,
761  int &baby_hig_bin, float &baby_mct, bool isSystemtic){
762  vector<int> hi_csv(5,-1); // Indices of the 5 jets with highest CSV
763  for (size_t ijet(0); ijet < baby_jets_csv.size(); ijet++) {
764  if (!isSystemtic) {
765  baby_jets_h1.push_back(false);
766  baby_jets_h2.push_back(false);
767  }
768 
769  //matters for WH to not work with lepton jets
770  // check pT in case running with systematics, where jets collection may contain jets below the threshold
771  if (baby_jets_islep[ijet] || baby_jets_p4[ijet].pt()<=jetTool->JetPtCut || fabs(baby_jets_p4[ijet].eta())>jetTool->JetEtaCut) continue;
772 
773  // Finding the N jets with highest CSV values
774  float csv = baby_jets_csv[ijet];
775  for(size_t ind(0); ind<hi_csv.size(); ind++){
776  float icsv(-99999);
777  if(hi_csv[ind]>=0) icsv = baby_jets_csv[hi_csv[ind]];
778  if(hi_csv[ind]==-1 || csv > icsv){
779  for(size_t ind2(hi_csv.size()-1); ind2>=ind+1; ind2--) hi_csv[ind2] = hi_csv[ind2-1];
780  hi_csv[ind] = ijet;
781  break;
782  }
783  } // Loop over highest CSV jets
784  } // Loop over jets
785 
787  if(hi_csv[3]>=0){
788  // hig_p4 has the p4 of the jet if row==col, and if not the sum of the p4 for the row-th and col-th jets
789  unsigned hi_csv_size = 4;
790  if (hi_csv[4]>=0) hi_csv_size = 5;
791  vector<vector<LVector> > hig_p4;
792  for(size_t row=0; row<hi_csv_size; row++){
793  hig_p4.push_back(vector<LVector>());
794  for(size_t col=0; col<=row; col++){
795  LVector jetp4(baby_jets_p4[hi_csv[row]]);
796  hig_p4.back().push_back(jetp4);
797  if(row!=col) hig_p4.back().back() += hig_p4[col][col];
798  } // Loop over columns in hig_p4
799  } // Loop over rows in hig_p4
800  // Loop over all possible Higgs combination in the first nCSVs jets of highest CSV
801  size_t nCSVs = hi_csv.size();
802  nCSVs = 4;
803  vector<int> hig_ind(4,0);
804  float minDm(9999.);
805  for(size_t ind0=0; ind0<nCSVs; ind0++){
806  for(size_t ind1=0; ind1<ind0; ind1++){
807  for(size_t ind2=ind0+1; ind2<nCSVs; ind2++){
808  if(ind2==ind1) continue;
809  for(size_t ind3=0; ind3<ind2; ind3++){
810  if(ind3==ind0 || ind3==ind1) continue;
811  float thisDm = fabs(hig_p4[ind0][ind1].mass() - hig_p4[ind2][ind3].mass());
812  if(thisDm < minDm) {
813  hig_ind[0] = ind0; hig_ind[1] = ind1; hig_ind[2] = ind2; hig_ind[3] = ind3;
814  minDm = thisDm;
815  }
816  } // ind3
817  } // ind2
818  } // ind1
819  } // ind0
820 
821  if (!isSystemtic) {
822  baby_jets_h1[hi_csv[hig_ind[0]]] = true; baby_jets_h1[hi_csv[hig_ind[1]]] = true;
823  baby_jets_h2[hi_csv[hig_ind[2]]] = true; baby_jets_h2[hi_csv[hig_ind[3]]] = true;
824  }
825 
826  LVector hig1 = hig_p4[hig_ind[0]][hig_ind[1]], hig2 = hig_p4[hig_ind[2]][hig_ind[3]];
827  baby_hig_dm = fabs(hig1.mass() - hig2.mass());
828  baby_hig_am = (hig1.mass() + hig2.mass())/2.;
829  baby_hig_drmax = max(deltaR(baby_jets_p4[hi_csv[hig_ind[0]]], baby_jets_p4[hi_csv[hig_ind[1]]]),
830  deltaR(baby_jets_p4[hi_csv[hig_ind[2]]], baby_jets_p4[hi_csv[hig_ind[3]]]));
831 
832  // Setting up the ABCD bin:
833  // 2 -> SIG, 1 -> SB, 0 -> in between, not used
834  if (!isSystemtic) {
835  if(baby_hig_dm<=40) {
836  if (baby_hig_am>100 && baby_hig_am<140) baby_hig_bin = 2;
837  else if (baby_hig_am<=200) baby_hig_bin = 1;
838  else baby_hig_bin = 0;
839  } else baby_hig_bin = 0;
840  // 20 -> 2b, 30 -> 3b, 40 -> 4b
841  if(baby_nbt>=2) {
842  baby_hig_bin += 20;
843  if(baby_nbm>=3) {
844  baby_hig_bin += 10;
845  if(baby_nbl>=4) baby_hig_bin += 10;
846  }
847  }
848  }
849  } // if njets >= 4
850  else if (hi_csv[1]>=0 && !isSystemtic){ // Checking that at least 2 jets exist (important for WH analysis)
851  baby_jets_h1[hi_csv[0]] = true;
852  baby_jets_h1[hi_csv[1]] = true;
853  }
854 
855  if (hi_csv[1]>=0 && !isSystemtic){
856  LVector pb1 = baby_jets_p4[hi_csv[0]], pb2 = baby_jets_p4[hi_csv[1]];
857  baby_mct = sqrt(2*pb1.pt()*pb2.pt() * (1+cos(deltaPhi(pb1.phi(), pb2.phi()))) );
858  }
859 
860  return;
861 }
862 
863 void bmaker_full::writeBBVars(std::vector<reco::Candidate::LorentzVector> &all_baby_jets, vCands &sig_leps){
864  // write deltaR between csvm jets
865  vector<size_t> branks = jet_met_tools::getBRanking(all_baby_jets, baby.jets_csv(), baby.jets_islep());
866  baby.bb_highcsv_idx() = -1;
867  jetTool->fillDeltaRbb(baby.dr_bb(),baby.bb_pt(),baby.bb_m(),baby.bb_jet_idx1(),baby.bb_jet_idx2(), baby.bb_gs_idx(), baby.bb_gs_flavor(),all_baby_jets, baby.jets_csv(), baby.jets_islep(),baby.jets_pt(),branks,baby.bb_highcsv_idx());
868 
869  vector<size_t> deepbranks = jet_met_tools::getBRanking(all_baby_jets, baby.jets_csvd(), baby.jets_islep());
870  baby.bb_highcsv_idx_deep() = -1;
871  jetTool->fillDeltaRbb(baby.dr_bb_deep(),baby.bb_pt_deep(),baby.bb_m_deep(),baby.bb_jet_idx1_deep(),baby.bb_jet_idx2_deep(), baby.bb_gs_idx_deep(), baby.bb_gs_flavor_deep(),all_baby_jets, baby.jets_csvd(), baby.jets_islep(),baby.jets_pt(),deepbranks,baby.bb_highcsv_idx_deep(),true);
872 
873  if(baby.nbm() >= 2 && sig_leps.size()>0){
874  const auto &jet1 = all_baby_jets.at(branks.at(0));
875  const auto &jet2 = all_baby_jets.at(branks.at(1));
876  double px = baby.met()*cos(baby.met_phi())+sig_leps.at(0)->px();
877  double py = baby.met()*sin(baby.met_phi())+sig_leps.at(0)->py();
878  baby.mt2() = getMT2(jet1.mass(), jet1.pt(), jet1.phi(),
879  jet2.mass(), jet2.pt(), jet2.phi(),
880  hypot(px, py), atan2(py, px));
881  baby.mt2_0mass() = getMT2(0., jet1.pt(), jet1.phi(),
882  0., jet2.pt(), jet2.phi(),
883  hypot(px, py), atan2(py, px));
884  }else{
885  baby.mt2() = -1.;
886  }
887 }
888 
890  vector<float> fdummy;
891  vector<int> idummy;
892  int intdummy;
893  bool cluster_leps = false;
894  jetTool->clusterFatJets(intdummy, baby.mj14_nolep(), fdummy, fdummy, fdummy,
895  fdummy, idummy, baby.jets_fjet14_nolep_index(),
896  baby, 1.4, jetTool->JetPtCut, cluster_leps);
897 
898  cluster_leps = true;
899  fdummy.clear(); idummy.clear();
900  // jetTool->clusterFatJets(intdummy, baby.mj08(), fdummy, fdummy, fdummy,
901  // fdummy, idummy, baby.jets_fjet08_index(),
902  // baby, 0.8, jetTool->JetPtCut, cluster_leps);
903  jetTool->clusterFatJets(baby.nfjets08(), baby.mj08(),
904  baby.fjets08_pt(), baby.fjets08_eta(),
905  baby.fjets08_phi(), baby.fjets08_m(),
906  baby.fjets08_nconst(), baby.jets_fjet08_index(),
907  baby, 0.8, jetTool->JetPtCut, cluster_leps);
908 
909  fdummy.clear(); idummy.clear();
910  jetTool->clusterFatJets(intdummy, baby.mj12(), fdummy, fdummy, fdummy,
911  fdummy, idummy, baby.jets_fjet12_index(),
912  baby, 1.2, jetTool->JetPtCut, cluster_leps);
913 
914  fdummy.clear(); idummy.clear();
915  jetTool->clusterFatJets(baby.nfjets14(), baby.mj14(),
916  baby.fjets14_pt(), baby.fjets14_eta(),
917  baby.fjets14_phi(), baby.fjets14_m(),
918  baby.fjets14_nconst(), baby.jets_fjet14_index(),
919  baby, 1.4, jetTool->JetPtCut, cluster_leps);
920 
921 }
922 /* DAK8
923 void bmaker_full::writeAk8Jets(edm::Handle<edm::View<pat::Jet>> &ak8jets){
924 
925  baby.nak8jets() = 0;
926  for (unsigned idx=0; idx<ak8jets->size(); ++idx){
927  const auto &jet = ak8jets->at(idx);
928  JetHelper jet_helper(&jet);
929  bool goodPtEta = jet.pt() > jetTool->JetPtCut && fabs(jet.eta()) <= jetTool->JetEtaCut;
930  if(goodPtEta) {
931  baby.nak8jets()++;
932  baby.ak8jets_pt().push_back(jet.pt());
933  baby.ak8jets_eta().push_back(jet.eta());
934  baby.ak8jets_phi().push_back(jet.phi());
935  baby.ak8jets_m().push_back(jet.mass());
936  // jet score prediction
937  // Nominal
938  const auto& nnpreds = fatjetNN_->predict(jet_helper);
939  FatJetNNHelper nn(nnpreds);
940  // binarized score (normalized by qcd)
941  baby.ak8jets_nom_bin_top().push_back(nn.get_binarized_score_top());
942  baby.ak8jets_nom_bin_w().push_back(nn.get_binarized_score_w());
943  baby.ak8jets_nom_bin_z().push_back(nn.get_binarized_score_z());
944  baby.ak8jets_nom_bin_zbb().push_back(nn.get_binarized_score_zbb());
945  baby.ak8jets_nom_bin_hbb().push_back(nn.get_binarized_score_hbb());
946  baby.ak8jets_nom_bin_h4q().push_back(nn.get_binarized_score_h4q());
947  // raw score
948  baby.ak8jets_nom_raw_top().push_back(nn.get_raw_score_top());
949  baby.ak8jets_nom_raw_w().push_back(nn.get_raw_score_w());
950  baby.ak8jets_nom_raw_z().push_back(nn.get_raw_score_z());
951  baby.ak8jets_nom_raw_zbb().push_back(nn.get_raw_score_zbb());
952  baby.ak8jets_nom_raw_hbb().push_back(nn.get_raw_score_hbb());
953  baby.ak8jets_nom_raw_h4q().push_back(nn.get_raw_score_h4q());
954  baby.ak8jets_nom_raw_qcd().push_back(nn.get_raw_score_qcd());
955  // Decorrelated
956  const auto& mdpreds = decorrNN_->predict(jet_helper);
957  FatJetNNHelper md(mdpreds);
958  // binarized score (normalized by qcd)
959  baby.ak8jets_decor_bin_top().push_back(md.get_binarized_score_top());
960  baby.ak8jets_decor_bin_w().push_back(md.get_binarized_score_w());
961  baby.ak8jets_decor_bin_z().push_back(md.get_binarized_score_z());
962  baby.ak8jets_decor_bin_zbb().push_back(md.get_binarized_score_zbb());
963  baby.ak8jets_decor_bin_hbb().push_back(md.get_binarized_score_hbb());
964  baby.ak8jets_decor_bin_h4q().push_back(md.get_binarized_score_h4q());
965  // raw score
966  baby.ak8jets_decor_raw_top().push_back(md.get_raw_score_top());
967  baby.ak8jets_decor_raw_w().push_back(md.get_raw_score_w());
968  baby.ak8jets_decor_raw_z().push_back(md.get_raw_score_z());
969  baby.ak8jets_decor_raw_zbb().push_back(md.get_raw_score_zbb());
970  baby.ak8jets_decor_raw_hbb().push_back(md.get_raw_score_hbb());
971  baby.ak8jets_decor_raw_h4q().push_back(md.get_raw_score_h4q());
972  baby.ak8jets_decor_raw_qcd().push_back(md.get_raw_score_qcd());
973  // flavor tagger
974  baby.ak8jets_decor_flav_bb().push_back(md.get_flavor_score_bb()); // H->bb + Z->bb + gluon->bb
975  baby.ak8jets_decor_flav_cc().push_back(md.get_flavor_score_cc()); // H->cc + Z->cc + gluon->cc
976  baby.ak8jets_decor_flav_bb_noglu().push_back(md.get_flavor_score_bb_no_gluon()); // H->bb + Z->bb
977  baby.ak8jets_decor_flav_cc_noglu().push_back(md.get_flavor_score_cc_no_gluon()); // H->cc + Z->cc
978  }
979  }
980 }
981 */
982 vCands bmaker_full::writeMuons(edm::Handle<pat::MuonCollection> muons,
983  edm::Handle<pat::PackedCandidateCollection> pfcands,
984  edm::Handle<reco::VertexCollection> vtx,
985  vCands &veto_mus, vCands &all_mus, double rhoEventCentral){
986  vCands sig_mus;
987  veto_mus.clear(); all_mus.clear();
988  baby.nmus() = 0; baby.nvmus() = 0;
989 
990  set<unsigned> badmu_idx, badmu_dupl_idx;
991  if (isData && !outname.Contains("Run2017") && !outname.Contains("Run2018")) {
992  badmu_idx = lepTool->badGlobalMuonSelector(vtx, muons, false);
993  badmu_dupl_idx = lepTool->badGlobalMuonSelector(vtx, muons, true);
994  }
995 
996  for (unsigned ilep(0); ilep < muons->size(); ilep++) {
997  const pat::Muon &lep = (*muons)[ilep];
998  //Save muons that were demoted from isPF in 8_0_26_patch1 in order to debug Giovanni badMuon flags
999  bool demoted(false);
1000  //userInt has old value of isPFMuon()
1001 
1002  if(!outname.Contains("Run2017") && !outname.Contains("Run2018"))
1003  if(isData && !lep.isPFMuon() && lep.userInt("muonsCleaned:oldPF")) demoted = true;
1004 
1005  bool isBadMu(false), isBadDuplMu(false);
1006  if (badmu_idx.find(ilep)!=badmu_idx.end()) isBadMu = true;
1007  if (badmu_dupl_idx.find(ilep)!=badmu_dupl_idx.end()) isBadDuplMu = true;
1008 
1009  bool isBadTrackerMuon(false);
1010  if (lep.pt()>10 && lep.isPFMuon() && lep.isTrackerMuon() && lep.numberOfMatchedStations()>0) {
1011  if (lep.innerTrack().isNonnull()) {
1012  bool goodQualityTrack = (lep.innerTrack()->numberOfValidHits()>=10) || (lep.innerTrack()->numberOfValidHits()>=7 && lep.innerTrack()->numberOfLostHits()==0);
1013  if (!goodQualityTrack) isBadTrackerMuon = true;
1014  } else{
1015  isBadTrackerMuon = true;
1016  }
1017  }
1018 
1019  // Storing leptons that pass all veto cuts except for iso
1020  bool save_mu = lepTool->isVetoMuon(lep, vtx, -99.) || isBadMu || isBadDuplMu || isBadTrackerMuon || demoted;
1021  if(!save_mu) continue;
1022 
1023 // double lep_iso(lepTool->getMinIsolation(dynamic_cast<const reco::Candidate *>(&lep), rhoEventCentral));
1024 // double lep_iso(lepTool->getMinIsolation(lepiso, rhoEventCentral));
1025 // const pat::PFIsolation lepiso = lep.miniPFIsolation();
1026 // double lep_iso = lepiso.chargedHadronIso() + lepiso.neutralHadronIso() + lepiso.photonIso();
1027 
1028  double lep_iso(lepTool->getPFIsolation(pfcands, dynamic_cast<const reco::Candidate *>(&lep), 0.05, 0.2, 10., rhoEventCentral, false));
1029  double lep_reliso(lepTool->getRelIsolation(lep, rhoEventCentral));
1030  double dz(0.), d0(0.);
1031  lepTool->vertexMuon(lep, vtx, dz, d0); // Calculating dz and d0
1032 
1033  baby.mus_bad().push_back(isBadMu);
1034  baby.mus_bad_dupl().push_back(isBadDuplMu);
1035  baby.mus_bad_trkmu().push_back(isBadTrackerMuon);
1036  baby.mus_demoted().push_back(demoted);
1037  baby.mus_pt().push_back(lep.pt());
1038  baby.mus_eta().push_back(lep.eta());
1039  baby.mus_phi().push_back(lep.phi());
1040  baby.mus_dz().push_back(dz);
1041  baby.mus_d0().push_back(d0);
1042  baby.mus_charge().push_back(lep.charge());
1043  baby.mus_sigid().push_back(lepTool->idMuon(lep, vtx, lepTool->kMedium));
1044  baby.mus_tight().push_back(lepTool->idMuon(lep, vtx, lepTool->kTight));
1045  baby.mus_miniso().push_back(lep_iso);
1046  baby.mus_reliso().push_back(lep_reliso);
1047  baby.mus_tm().push_back(false); // Filled in writeMC
1048  baby.mus_inz().push_back(false); // Filled in writeDiLep
1049  baby.mus_inzv().push_back(false); // Filled in writeDiLep
1050  baby.mus_vvvl().push_back(false); // Filled in writeHLTObjects
1051  baby.mus_isomu18().push_back(false); // Filled in writeHLTObjects
1052  baby.mus_mu50().push_back(false); // Filled in writeHLTObjects
1053  baby.mus_mu8().push_back(false); // Filled in writeHLTObjects
1054  if (lep.track().isNonnull()){
1055  baby.mus_trk_quality().push_back(lep.innerTrack()->quality(reco::TrackBase::highPurity));
1056  baby.mus_pterr().push_back(lep.innerTrack()->ptError());
1057  baby.mus_trk_nholes_in().push_back(lep.innerTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
1058  baby.mus_trk_nholes_out().push_back(lep.innerTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS));
1059  baby.mus_trk_algo().push_back(lep.innerTrack()->algo());
1060  } else {
1061  baby.mus_trk_quality().push_back(-99999);
1062  baby.mus_pterr().push_back(-99999);
1063  baby.mus_trk_nholes_in().push_back(-99999);
1064  baby.mus_trk_nholes_out().push_back(-99999);
1065  baby.mus_trk_algo().push_back(-99999);
1066  }
1067  baby.mus_em_e().push_back(lep.calEnergy().em);
1068  baby.mus_had_e().push_back(lep.calEnergy().had);
1069  all_mus.push_back(dynamic_cast<const reco::Candidate *>(&lep)); // For truth-matching in writeMC
1070 
1071  if(lepTool->isVetoMuon(lep, vtx, lep_iso)) {
1072  baby.nvmus()++;
1073  veto_mus.push_back(dynamic_cast<const reco::Candidate *>(&lep));
1074  }
1075  if(lepTool->isSignalMuon(lep, vtx, lep_iso)) {
1076  baby.nmus()++;
1077  sig_mus.push_back(dynamic_cast<const reco::Candidate *>(&lep));
1078  baby.mus_sig().push_back(true);
1079  } else baby.mus_sig().push_back(false);
1080  } // Loop over muons
1081 
1082  return sig_mus;
1083 }
1084 
1085 
1086 vCands bmaker_full::writeElectrons(edm::Handle<pat::ElectronCollection> electrons,
1087  edm::Handle<pat::PackedCandidateCollection> pfcands,
1088  edm::Handle<reco::VertexCollection> vtx,
1089  vCands &veto_els, vCands &all_els, double rhoEventCentral){
1090  vCands sig_els;
1091  veto_els.clear(); all_els.clear();
1092  baby.nels() = 0; baby.nvels() = 0;
1093  for (size_t ilep(0); ilep < electrons->size(); ilep++) {
1094  const pat::Electron &lep = (*electrons)[ilep];
1095  if(!lepTool->isVetoElectron(lep, vtx, -99.)) continue; // Storing leptons that pass all veto cuts except for iso
1096 
1097 // double lep_iso(lepTool->getMinIsolation(dynamic_cast<const reco::Candidate *>(&lep), rhoEventCentral));
1098 // double lep_iso(lepTool->getMinIsolation(lepiso, rhoEventCentral));
1099 // const PFIsolation lepiso = lep.miniPFIsolation();
1100 // const double lep_iso = lepiso.chargedHadronIso() + lepiso.neutralHadronIso() + lepiso.photonIso();
1101 
1102  double lep_iso(lepTool->getPFIsolation(pfcands, dynamic_cast<const reco::Candidate *>(&lep), 0.05, 0.2, 10., rhoEventCentral, false));
1103  double lep_reliso(lepTool->getRelIsolation(lep, rhoEventCentral));
1104  double dz(0.), d0(0.);
1105  lepTool->vertexElectron(lep, vtx, dz, d0); // Calculating dz and d0
1106 
1107  baby.els_pt().push_back(lep.pt());
1108 
1109  baby.els_scpt().push_back(lep.superCluster()->energy()*sin(lep.superClusterPosition().theta()));
1110  baby.els_sceta().push_back(lep.superCluster()->eta());
1111  baby.els_eta().push_back(lep.eta());
1112  baby.els_phi().push_back(lep.phi());
1113  baby.els_dz().push_back(dz);
1114  baby.els_d0().push_back(d0);
1115  baby.els_ip3d().push_back(lep.ip3d());
1116  baby.els_charge().push_back(lep.charge());
1117  baby.els_sigid().push_back(lepTool->idElectron(lep, vtx, lepTool->kMedium));
1118  baby.els_ispf().push_back(lep.numberOfSourceCandidatePtrs()==2 && abs(lep.sourceCandidatePtr(1)->pdgId())==11);
1119  baby.els_tight().push_back(lepTool->idElectron(lep, vtx, lepTool->kTight));
1120  baby.els_miniso().push_back(lep_iso);
1121  baby.els_reliso().push_back(lep_reliso);
1122  baby.els_tm().push_back(false); // Filled in writeMC
1123  baby.els_inz().push_back(false); // Filled in writeDiLep
1124  baby.els_inzv().push_back(false); // Filled in writeDiLep
1125  baby.els_vvvl().push_back(false); // Filled in writeHLTObjects
1126  baby.els_ele23().push_back(false); // Filled in writeHLTObjects
1127  baby.els_ele105().push_back(false); // Filled in writeHLTObjects
1128  baby.els_ele8().push_back(false); // Filled in writeHLTObjects
1129  baby.els_hovere().push_back(lep.hadronicOverEm());
1130  baby.els_eoverp().push_back(lep.eSuperClusterOverP());
1131  baby.els_em_e().push_back(lep.ecalEnergy());
1132  baby.els_deta_sctrk().push_back(lep.deltaEtaSuperClusterTrackAtVtx());
1133  baby.els_dphi_sctrk().push_back(lep.deltaPhiSuperClusterTrackAtVtx());
1134  if (lep.gsfTrack().isNonnull()){
1135  baby.els_trk_pt().push_back(lep.gsfTrack()->pt());
1136  baby.els_trk_pterr().push_back(lep.gsfTrack()->ptError());
1137  baby.els_trk_nholes().push_back(lep.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
1138  } else {
1139  baby.els_trk_pt().push_back(-99999);
1140  baby.els_trk_pterr().push_back(-99999);
1141  baby.els_trk_nholes().push_back(-99999);
1142  }
1143  all_els.push_back(dynamic_cast<const reco::Candidate *>(&lep)); // For truth-matching in writeMC
1144 
1145  if(lepTool->isVetoElectron(lep, vtx, lep_iso)){
1146  baby.nvels()++;
1147  veto_els.push_back(dynamic_cast<const reco::Candidate *>(&lep));
1148  }
1149  if(lepTool->isSignalElectron(lep, vtx, lep_iso)) {
1150  baby.nels()++;
1151  sig_els.push_back(dynamic_cast<const reco::Candidate *>(&lep));
1152  baby.els_sig().push_back(true);
1153  } else baby.els_sig().push_back(false);
1154  } // Loop over electrons
1155 
1156  return sig_els;
1157 }
1158 
1160  baby.nleps() = baby.nmus() + baby.nels();
1161  baby.nvleps() = baby.nvmus() + baby.nvels();
1162  sort(leptons.begin(), leptons.end(), greaterPt);
1163  for(size_t ind(0); ind < leptons.size(); ind++) {
1164  baby.leps_pt().push_back(leptons[ind]->pt());
1165  baby.leps_eta().push_back(leptons[ind]->eta());
1166  baby.leps_phi().push_back(leptons[ind]->phi());
1167  baby.leps_id().push_back(leptons[ind]->pdgId());
1168  } // Loop over leptons
1169 }
1170 
1171 void bmaker_full::writeTks(edm::Handle<pat::PackedCandidateCollection> pfcands,edm::Handle<reco::VertexCollection> vtx, double rhoEventCentral ){
1172  if(baby.leps_id().size()>0){
1173  vector<float> isos;
1174  vector<float> relisos;
1175  vCands ra4tks;
1176  ra4tks = lepTool->getRA4IsoTracks(pfcands, baby.met(), baby.met_phi(),rhoEventCentral,isos,relisos,baby.leps_id().at(0));
1177 
1178  int nveto=0;
1179 
1180  for(unsigned i=0;i<ra4tks.size();i++){
1181  baby.tks_pt().push_back(ra4tks.at(i)->pt());
1182  baby.tks_eta().push_back(ra4tks.at(i)->eta());
1183  baby.tks_phi().push_back(ra4tks.at(i)->phi());
1184  baby.tks_d0().push_back( sqrt(pow(ra4tks.at(i)->vx()-vtx->at(0).x(),2) + pow(vtx->at(0).y()-ra4tks.at(i)->vy(),2)));
1185  baby.tks_dz().push_back(ra4tks.at(i)->vz()-vtx->at(0).z());
1186  baby.tks_pdg().push_back(ra4tks.at(i)->pdgId());
1187  baby.tks_miniso().push_back(isos.at(i));
1188  baby.tks_reliso().push_back(relisos.at(i));
1189  baby.tks_tm().push_back(false); //filled in writeMC
1190 
1191  baby.tks_mt2().push_back(getMT2(baby.leps_pt().at(0),baby.leps_phi().at(0),baby.tks_pt().back(),baby.tks_phi().back(),baby.met(),baby.met_phi()));
1192  baby.tks_mt().push_back(getMT(baby.tks_pt().back(),baby.tks_phi().back(),baby.met(),baby.met_phi()));
1193 
1194  if(fabs(baby.tks_pdg().back())==211 && baby.tks_pt().back()>15. && baby.tks_miniso().back()<0.1 && baby.tks_mt2().back()<60 && baby.tks_dz().back()<0.07 && baby.tks_d0().back()<0.05 ) nveto++;
1195  else if (fabs(baby.tks_pdg().back())==13 && baby.tks_pt().back()>10. && baby.tks_miniso().back()<0.2 && baby.tks_mt2().back()<80 && baby.tks_dz().back()<0.07 && baby.tks_d0().back()<0.05) nveto++;
1196  else if (fabs(baby.tks_pdg().back())==11 && baby.tks_pt().back()>10. && baby.tks_miniso().back()<0.2 && baby.tks_mt2().back()<80 && baby.tks_dz().back()<0.07 && baby.tks_d0().back()<0.05) nveto++;
1197 
1198  }
1199 
1200  baby.nveto()=nveto;
1201 
1202  }
1203 } // if goodPV
1204 
1205 
1206 void bmaker_full::writeDiLep(vCands &sig_mus, vCands &sig_els, vCands &veto_mus, vCands &veto_els){
1219  setElMuMass(sig_els, sig_mus, &baby_base::elmu_m, &baby_base::elmu_pt1, &baby_base::elmu_pt2, &baby_base::elmu_pt,
1222  // setElMuMass(veto_els, veto_mus, &baby_base::elmuv_m, &baby_base::elmuv_pt1, &baby_base::elmuv_pt2, &baby_base::elmuv_pt,
1223  // &baby_base::elmuv_eta, &baby_base::elmuv_phi);
1224 }
1225 
1227  baby_float ll_pt, baby_float ll_eta, baby_float ll_phi, baby_vfloat l_pt, baby_vbool l_inz,
1228  baby_float ll_w){
1229  for(size_t lep1(0); lep1 < leptons.size(); lep1++){
1230  for(size_t lep2(lep1+1); lep2 < leptons.size(); lep2++){
1231  if(leptons[lep1]->charge()*leptons[lep2]->charge()<0){
1232  LVector z_p4(leptons[lep1]->p4());
1233  z_p4 += leptons[lep2]->p4();
1234  (baby.*ll_m)() = z_p4.mass();
1235  (baby.*ll_pt)() = z_p4.pt();
1236  (baby.*ll_eta)() = z_p4.eta();
1237  (baby.*ll_phi)() = z_p4.phi();
1238  float pt1(leptons[lep1]->pt()), pt2(leptons[lep2]->pt());
1239  (baby.*ll_pt1)() = max(pt1, pt2);
1240  (baby.*ll_pt2)() = min(pt1, pt2);
1241  for(size_t ilep(0); ilep < (baby.*l_pt)().size(); ilep++){
1242  if(fabs(pt1 - (baby.*l_pt)()[ilep]) < 1e-7) (baby.*l_inz)()[ilep] = true;
1243  if(fabs(pt2 - (baby.*l_pt)()[ilep]) < 1e-7) (baby.*l_inz)()[ilep] = true;
1244  }
1245  (baby.*ll_w)() = lepton_tools::getScaleFactor({leptons[lep1], leptons[lep2]}).first;
1246  return; // We only set it with the first good ll combination
1247  }
1248  } // Loop over lep2
1249  } // Loop over lep1
1250 }
1251 
1252 void bmaker_full::setElMuMass(vCands leptons1, vCands leptons2, baby_float ll_m, baby_float ll_pt1, baby_float ll_pt2,
1253  baby_float ll_pt, baby_float ll_eta, baby_float ll_phi,
1254  baby_float ll_w){
1255  for(size_t lep1(0); lep1 < leptons1.size(); lep1++){
1256  for(size_t lep2(0); lep2 < leptons2.size(); lep2++){
1257  if(leptons1[lep1]->charge()*leptons2[lep2]->charge()<0){
1258  LVector z_p4(leptons1[lep1]->p4());
1259  z_p4 += leptons2[lep2]->p4();
1260  (baby.*ll_m)() = z_p4.mass();
1261  (baby.*ll_pt)() = z_p4.pt();
1262  (baby.*ll_eta)() = z_p4.eta();
1263  (baby.*ll_phi)() = z_p4.phi();
1264  float pt1(leptons1[lep1]->pt()), pt2(leptons2[lep2]->pt());
1265  (baby.*ll_pt1)() = pt1;
1266  (baby.*ll_pt2)() = pt2;
1267  (baby.*ll_w)() = lepton_tools::getScaleFactor({leptons1[lep1], leptons2[lep2]}).first;
1268  return; // We only set it with the first good ll combination
1269  }
1270  } // Loop over lep2
1271  } // Loop over lep1
1272 }
1273 
1274 
1275 vCands bmaker_full::writePhotons(edm::Handle<pat::PhotonCollection> allphotons,
1276  edm::Handle<std::vector<pat::Electron> > &electrons,
1277  edm::Handle<reco::ConversionCollection> &conversions,
1278  edm::Handle<reco::BeamSpot> &beamspot, double rho){
1279  vCands photons;
1280  baby.nph() = 0;
1281 
1282  for (size_t ind(0); ind < allphotons->size(); ind++) {
1283  const pat::Photon &photon = (*allphotons)[ind];
1284 
1285  if(photon.pt() < 50) continue;
1286  if(!photonTool->idPhoton(photon, electrons, conversions, beamspot, rho)) continue;
1287  if(photon.pt() > photonTool->PhotonPtCut) baby.nph()++;
1288  baby.ph_pt().push_back(photon.pt());
1289  baby.ph_eta().push_back(photon.eta());
1290  baby.ph_phi().push_back(photon.phi());
1291  baby.ph_tm().push_back(false); // Filled in writeMC
1292  baby.ph_ph90().push_back(false); // Filled in writeHLTObjects
1293 
1294  photons.push_back(dynamic_cast<const reco::Candidate *>(&photon));
1295  } // Loop over allphotons
1296 
1297  return photons;
1298 } // writePhotons
1299 
1300 bool bmaker_full::writeTriggers(const edm::TriggerNames &names,
1301  edm::Handle<edm::TriggerResults> triggerBits,
1302  edm::Handle<pat::PackedTriggerPrescales> triggerPrescales){
1303  bool keep(false);
1304  baby.trig().resize(trig_name.size(), false);
1305  baby.trig_prescale().resize(trig_name.size(), -1.);
1306  for (size_t itrig(0); itrig < triggerBits->size(); itrig++) {
1307  for(size_t itn(0); itn < trig_name.size(); itn++){
1308  if(names.triggerName(itrig).find(trig_name[itn])!=string::npos){
1309  baby.trig()[itn] = triggerBits->accept(itrig);
1310  baby.trig_prescale()[itn] = triggerPrescales->getPrescaleForIndex(itrig);
1311  if(baby.trig()[itn]) keep = true;
1312  }
1313  }
1314  }
1315  // OR-ing triggers used in RA4
1316  vector<TString> trigs_ra4({"HLT_PFHT500_PFMET100_PFMHT100_IDTight_v","HLT_Mu15_IsoVVVL_PFHT450_v",
1317  "HLT_Ele15_IsoVVVL_PFHT450_v","HLT_PFMET120_PFMHT120_IDTight_v",
1318  "HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_v","HLT_IsoMu24_v","HLT_Mu50_v",
1319  "HLT_Ele115_CaloIdVT_GsfTrkIdT_v","HLT_Ele27_WPTight_Gsf_v",
1320  "HLT_Ele28_eta2p1_WPTight_Gsf_HT150_v"});
1321  vector<TString> trigs_met({"HLT_PFHT500_PFMET100_PFMHT100_IDTight_v","HLT_PFMET120_PFMHT120_IDTight_v",
1322  "HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_v"});
1323  vector<TString> trigs_vvvl({"HLT_Mu15_IsoVVVL_PFHT350_v", "HLT_Mu15_IsoVVVL_PFHT400_v", "HLT_Mu50_IsoVVVL_PFHT400_v",
1324  "HLT_Ele15_IsoVVVL_PFHT350_v", "HLT_Ele15_IsoVVVL_PFHT400_v", "HLT_Ele50_IsoVVVL_PFHT400_v"});
1325  vector<TString> trigs_lep({"HLT_IsoMu24_v","HLT_Mu50_v","HLT_Ele115_CaloIdVT_GsfTrkIdT_v","HLT_Ele28_eta2p1_WPTight_Gsf_HT150_v"});
1326  baby.trig_ra4() = false;
1327  baby.trig_met() = false;
1328  baby.trig_vvvl() = false;
1329  baby.trig_lep() = false;
1330  for(size_t itn(0); itn < trig_name.size(); itn++){
1331  for(size_t ira4(0); ira4 < trigs_ra4.size(); ira4++)
1332  if(trig_name[itn].Contains(trigs_ra4[ira4])) baby.trig_ra4() = (baby.trig_ra4() || baby.trig()[itn]);
1333  for(size_t imet(0); imet < trigs_met.size(); imet++)
1334  if(trig_name[itn].Contains(trigs_met[imet])) baby.trig_met() = (baby.trig_met() || baby.trig()[itn]);
1335  for(size_t ivvvl(0); ivvvl < trigs_vvvl.size(); ivvvl++)
1336  if(trig_name[itn].Contains(trigs_vvvl[ivvvl])) baby.trig_vvvl() = (baby.trig_vvvl() || baby.trig()[itn]);
1337  for(size_t ilep(0); ilep < trigs_lep.size(); ilep++)
1338  if(trig_name[itn].Contains(trigs_lep[ilep])) baby.trig_lep() = (baby.trig_lep() || baby.trig()[itn]);
1339  } // Loop over trigger names
1340  return keep;
1341 }
1342 
1343 void bmaker_full::writeHLTObjects(const edm::TriggerNames &names,
1344  edm::Handle<pat::TriggerObjectStandAloneCollection> triggerObjects,
1345  vCands &all_mus, vCands &all_els, vCands &photons, const edm::Event& iEvent){
1346  baby.nmus_vvvl() = 0; baby.nmus_isomu18() = 0;
1347  baby.nels_vvvl() = 0; baby.nels_ele23() = 0;
1348  const float relptThreshold(1), drThreshold(0.3);
1349  edm::Handle<edm::TriggerResults> triggerBits;
1350  iEvent.getByToken(tok_trigResults_hlt_,triggerBits);
1351  for (pat::TriggerObjectStandAlone obj : *triggerObjects) {
1352  obj.unpackFilterLabels(iEvent, *triggerBits);
1353  obj.unpackPathNames(names);
1354  TString name(obj.collection());
1355  float objpt(obj.pt());
1356  if(name=="hltPFMETProducer::HLT") baby.onmet() = objpt;
1357  if(name=="hltPFHT::HLT") baby.onht() = objpt; // There's 2 of these, and we want the last one
1358 
1359  if(name=="hltL3MuonCandidates::HLT"){
1360  bool vvvl(obj.hasFilterLabel("hltL3MuVVVLIsoFIlter"));
1361  bool isomu18(obj.hasFilterLabel("hltL3crIsoL1sMu16L1f0L2f10QL3f18QL3trkIsoFiltered0p09") );
1362  bool mu50(obj.hasFilterLabel("hltL3fL1sMu16orMu25L1f0L2f10QL3Filtered50Q"));
1363  bool mu8(obj.hasFilterLabel("hltDoubleMu8Mass8L3Filtered"));
1364  if(vvvl) {
1365  baby.nmus_vvvl()++;
1366  baby.mus_vvvl_pt().push_back(objpt);
1367  baby.mus_vvvl_eta().push_back(obj.eta());
1368  baby.mus_vvvl_phi().push_back(obj.phi());
1369  }
1370  if(isomu18) baby.nmus_isomu18()++;
1371  if(vvvl && baby.onmu_vvvl()<objpt) baby.onmu_vvvl() = objpt;
1372  if(isomu18 && baby.onmu_isomu18()<objpt) baby.onmu_isomu18() = objpt;
1373  if(mu50 && baby.onmu_mu50()<objpt) baby.onmu_mu50() = objpt;
1374  if(mu8 && baby.onmu_mu8()<objpt) baby.onmu_mu8() = objpt;
1375  double mindr(999.);
1376  int minind(-1);
1377  if(vvvl || isomu18 || mu50 || mu8){
1378  for(size_t ind(0); ind < all_mus.size(); ind++) {
1379  double dr(deltaR(obj, *(all_mus[ind])));
1380  //double drelpt(fabs((all_mus[ind]->pt() - objpt)/objpt));
1381  //if(dr > drThreshold || drelpt > relptThreshold) continue;
1382  if(dr > drThreshold) continue;
1383  if(dr < mindr){
1384  mindr = dr;
1385  minind = ind;
1386  }
1387  } // Loop over reco muons
1388  if(minind>=0){
1389  baby.mus_vvvl()[minind] = vvvl;
1390  baby.mus_isomu18()[minind] = isomu18;
1391  baby.mus_mu50()[minind] = mu50;
1392  baby.mus_mu8()[minind] = mu8;
1393  }
1394  } // At least one match
1395  }
1396  if(name=="hltEgammaCandidates::HLT" || name=="hltDoubleEle8HLTPixelMatchElectronProducer::HLT"){
1397  bool vvvl(obj.hasFilterLabel("hltEle15VVVLGsfTrackIsoFilter"));
1398  bool ele23(obj.hasFilterLabel("hltEle23WPLooseGsfTrackIsoFilter") );
1399  bool ele105(obj.hasFilterLabel("hltEle105CaloIdVTGsfTrkIdTGsfDphiFilter"));
1400  bool ele8(obj.hasFilterLabel("hltDoubleEle8Mass8Filter"));
1401  if(vvvl) {
1402  baby.nels_vvvl()++;
1403  baby.els_vvvl_pt().push_back(objpt);
1404  baby.els_vvvl_eta().push_back(obj.eta());
1405  baby.els_vvvl_phi().push_back(obj.phi());
1406  }
1407  if(ele23) baby.nels_ele23()++;
1408  if(vvvl && baby.onel_vvvl()<objpt) baby.onel_vvvl() = objpt;
1409  if(ele23 && baby.onel_ele23()<objpt) baby.onel_ele23() = objpt;
1410  if(ele105 && baby.onel_ele105()<objpt) baby.onel_ele105() = objpt;
1411  if(ele8 && baby.onel_ele8()<objpt) baby.onel_ele8() = objpt;
1412  if(vvvl || ele23 || ele105 || ele8){
1413  double mindr(999.);
1414  int minind(-1);
1415  for(size_t ind(0); ind < all_els.size(); ind++) {
1416  double dr(deltaR(obj, *(all_els[ind])));
1417  double drelpt(fabs((all_els[ind]->pt() - objpt)/objpt));
1418  if(dr > drThreshold || drelpt > relptThreshold) continue;
1419  if(dr < mindr){
1420  mindr = dr;
1421  minind = ind;
1422  }
1423  } // Loop over reco elecrons
1424  if(minind>=0){
1425  baby.els_vvvl()[minind] = vvvl;
1426  baby.els_ele23()[minind] = ele23;
1427  baby.els_ele105()[minind] = ele105;
1428  baby.els_ele8()[minind] = ele8;
1429  }
1430  } // At least one electron match
1431  bool ph90(obj.hasFilterLabel("hltEG90L1SingleEG40HEFilter"));
1432  if(ph90 && baby.onph_ph90()<objpt) baby.onph_ph90() = objpt;
1433  if(ph90){
1434  double mindr(999.);
1435  int minind(-1);
1436  for(size_t ind(0); ind < photons.size(); ind++) {
1437  double dr(deltaR(obj, *(photons[ind])));
1438  double drelpt(fabs((photons[ind]->pt() - objpt)/objpt));
1439  if(dr > drThreshold || drelpt > relptThreshold) continue;
1440  if(dr < mindr){
1441  mindr = dr;
1442  minind = ind;
1443  }
1444  } // Loop over reco photons
1445  if(minind>=0){
1446  baby.ph_ph90()[minind] = ph90;
1447  }
1448  } // At least one photon match
1449  }
1450  } // Loop over trigger objects
1451 }
1452 
1453 // From https://twiki.cern.ch/twiki/bin/viewauth/CMS/MissingETOptionalFiltersRun2
1454 // Last updated on Jan 29, 2018
1455 void bmaker_full::writeFilters(const edm::TriggerNames &fnames,
1456  edm::Handle<edm::TriggerResults> filterBits,
1457  edm::Handle<reco::VertexCollection> vtx,
1458  vector<double> jetsMuonEnergyFrac){
1459  baby.pass_goodv() = true; baby.pass_cschalo() = true; baby.pass_eebadsc() = true;
1460  baby.pass_hbhe() = true; baby.pass_hbheiso() = true; baby.pass_ecaldeadcell() = true;
1461  baby.pass_badpfmu() = true; baby.pass_badchhad() = true; baby.pass_badcalib() = true;
1462  baby.pass_ra2_badmu() = true;
1463  for (size_t i(0); i < filterBits->size(); ++i) {
1464  string name = fnames.triggerName(i);
1465  bool pass = static_cast<bool>(filterBits->accept(i));
1466  if (name=="Flag_goodVertices") baby.pass_goodv() = pass;
1467  else if (name=="Flag_globalTightHalo2016Filter") baby.pass_cschalo() = pass;
1468  else if (name=="Flag_HBHENoiseFilter") baby.pass_hbhe() = pass;
1469  else if (name=="Flag_HBHENoiseIsoFilter") baby.pass_hbheiso() = pass;
1470  else if (name=="Flag_EcalDeadCellTriggerPrimitiveFilter") baby.pass_ecaldeadcell() = pass;
1471  //These filter bad events and must be inverted
1472  else if (name=="Flag_BadPFMuonFilter") baby.pass_badpfmu() = pass;
1473  else if (name=="BadChargedCandidateFilter") baby.pass_badchhad() = !pass;
1474  // else if (name=="Flag_eeBadScFilter") baby.pass_eebadsc() = pass; // Not recommended
1475  else if (name=="Flag_ecalBadCalibFilter") baby.pass_badcalib() = pass;
1476  }
1477 
1478  baby.pass() = baby.pass_goodv() && baby.pass_ecaldeadcell() &&
1479  baby.pass_hbhe() && baby.pass_hbheiso() &&
1480  baby.pass_badpfmu() && baby.pass_badchhad() && baby.pass_badcalib() &&
1481  baby.pass_jets() && baby.pass_fsmet() && baby.pass_fsjets();
1482 
1483  baby.pass_ra2() = baby.pass_goodv() && baby.pass_ecaldeadcell() &&
1484  baby.pass_hbhe() && baby.pass_hbheiso() &&
1485  baby.pass_badpfmu() && baby.pass_badchhad() && baby.pass_badcalib() &&
1486  baby.pass_jets_ra2() && baby.pass_fsmet() && baby.pass_fsjets();
1487 
1488  baby.pass_nohf() = baby.pass_goodv() && baby.pass_ecaldeadcell() &&
1489  baby.pass_hbhe() && baby.pass_hbheiso() &&
1490  baby.pass_badpfmu() && baby.pass_badchhad() && baby.pass_badcalib() &&
1491  baby.pass_jets_nohf() && baby.pass_fsmet() && baby.pass_fsjets();
1492 
1493 
1494  //Some filters not recommended for MC
1495  //https://twiki.cern.ch/twiki/bin/view/CMS/MissingETOptionalFiltersRun2#Moriond_2017
1496 
1497  if(isData){ baby.pass() = baby.pass() && baby.pass_eebadsc();
1498  baby.pass_ra2() = baby.pass_ra2() && baby.pass_eebadsc();
1499  baby.pass_nohf() = baby.pass_nohf() && baby.pass_eebadsc();
1500  }
1501 
1502  // Suggested only for Data and Fullsim
1503  if(!isFastSim){ baby.pass() = baby.pass() && baby.pass_cschalo();
1504  baby.pass_ra2() = baby.pass_ra2() && baby.pass_cschalo();
1505  baby.pass_nohf() = baby.pass_nohf() && baby.pass_cschalo();
1506  }
1507 
1508 
1509  for (size_t ijet(0); ijet < baby.jets_pt().size(); ijet++){
1510  if (abs(baby.jets_eta()[ijet])>2.4 || baby.jets_pt()[ijet]<200.) continue;
1511  if (baby.jets_islep()[ijet]) continue;
1512  if (jetsMuonEnergyFrac[ijet]<=0.5) continue;
1513  if (abs(dPhi(baby.jets_phi()[ijet],baby.met_phi()))<(3.14159-0.4)) continue;
1514  baby.pass_ra2_badmu() = false;
1515  break;
1516  }
1517 
1518  if (doSystematics){
1519  for (unsigned isys(0); isys<kSysLast; isys++){
1520  // sys_pass_jets already stored in the value of this variable in the baby
1521  baby.sys_pass()[isys] = baby.sys_pass()[isys] && baby.pass_goodv() && baby.pass_eebadsc() && baby.pass_cschalo() && baby.pass_hbhe() && baby.pass_hbheiso() && baby.pass_ecaldeadcell();
1522  }
1523  }
1524 }
1525 
1526 void bmaker_full::writeVertices(edm::Handle<reco::VertexCollection> vtx,
1527  edm::Handle<std::vector< PileupSummaryInfo > > pu_info){
1528  baby.npv() = vtx->size();
1529  if(pu_info.isValid()){
1530  for(size_t bc(0); bc<pu_info->size(); ++bc){
1531  if(pu_info->at(bc).getBunchCrossing()==0){
1532  baby.ntrupv() = pu_info->at(bc).getPU_NumInteractions();
1533  baby.ntrupv_mean() = pu_info->at(bc).getTrueNumInteractions();
1534  break;
1535  }
1536  } // Loop over true PVs
1537  } // if pu_info is valid
1538 }
1539 
1540 void bmaker_full::writeGenInfo(edm::Handle<LHEEventProduct> lhe_info){
1541  baby.nisr_me()=0; baby.ht_isr_me()=0.; baby.ptll_me()=0.;
1542  float px1(0.), px2(0.), py1(0.), py2(0.);
1543  unsigned nme_leps(0);
1544  for ( unsigned int icount = 0 ; icount < (unsigned int)lhe_info->hepeup().NUP; icount++ ) {
1545  unsigned int pdgid = abs(lhe_info->hepeup().IDUP[icount]);
1546  int status = lhe_info->hepeup().ISTUP[icount];
1547  int mom1_idx = lhe_info->hepeup().MOTHUP[icount].first;
1548  int mom2_idx = lhe_info->hepeup().MOTHUP[icount].second;
1549  int mom1id = mom1_idx==0 ? 0 : abs(lhe_info->hepeup().IDUP[mom1_idx-1]);
1550  int mom2id = mom2_idx==0 ? 0 : abs(lhe_info->hepeup().IDUP[mom2_idx-1]);
1551  float px = (lhe_info->hepeup().PUP[icount])[0];
1552  float py = (lhe_info->hepeup().PUP[icount])[1];
1553  float pt = sqrt(px*px+py*py);
1554 
1555 
1556  if(status==1 && (pdgid<6 || pdgid==21) && mom1id!=6 && mom2id!=6 && mom1id!=24 && mom2id!=24
1557  && mom1id!=23 && mom2id!=23 && mom1id!=25 && mom2id!=25) {
1558  baby.nisr_me()++;
1559  baby.ht_isr_me() += pt;
1560  }
1561 
1562  if(status==1 && (pdgid==11 || pdgid==13 || pdgid==15)) {
1563  if (nme_leps==0){
1564  px1 = px; py1 = py;
1565  } else if (nme_leps==1){
1566  px2 = px; py2 = py;
1567  } else cout<<"Found more than two leptons in ME"<<endl;
1568  nme_leps++;
1569  }
1570  } // Loop over generator particles
1571  if (nme_leps==2) {
1572  baby.ptll_me() = sqrt(pow(px1+px2,2)+pow(py1+py2,2));
1573  }
1574 
1575  if (outname.Contains("SMS-") && outname.Contains("PUSpring16Fast")){ //Get mgluino and mlsp
1576  typedef std::vector<std::string>::const_iterator comments_const_iterator;
1577 
1578  comments_const_iterator c_begin = lhe_info->comments_begin();
1579  comments_const_iterator c_end = lhe_info->comments_end();
1580 
1581  TString model_params;
1582  for(comments_const_iterator cit=c_begin; cit!=c_end; ++cit) {
1583  size_t found = (*cit).find("model");
1584  if(found != std::string::npos) {
1585  // std::cout <<"BABYMAKER: "<< *cit <<"end"<< std::endl;
1586  model_params = *cit;
1587  }
1588  }
1589 
1590  mcTool->getMassPoints(model_params,baby.mgluino(),baby.mlsp());
1591  }
1592 } // writeGenInfo
1593 
1594 void bmaker_full::writeIFSR(edm::Handle<reco::GenParticleCollection> genParticles,
1595  vector<reco::Candidate::LorentzVector> &all_baby_jets){
1596  bool verbose = false;
1597  baby.jets_isisr().resize(all_baby_jets.size(), false);
1598  int nisr(0);
1599  for (size_t ijet(0); ijet<all_baby_jets.size(); ijet++){
1600  if (baby.jets_islep()[ijet] || all_baby_jets[ijet].pt()<=jetTool->JetPtCut) continue;
1601 
1602  bool matched=false;
1603  for (size_t imc(0); imc < genParticles->size(); imc++) {
1604  if (matched) break;
1605  const reco::GenParticle &mc = (*genParticles)[imc];
1606  if (mc.status()!=23 || abs(mc.pdgId())>5) continue;
1607  int momid = abs(mc.mother()->pdgId());
1608  if(!(momid==6 || momid==23 || momid==24 || momid==25 || momid>1e6)) continue;
1609  //check against daughter in case of hard initial splitting
1610  for (size_t idau(0); idau < mc.numberOfDaughters(); idau++) {
1611  float dR = deltaR(all_baby_jets[ijet], mc.daughter(idau)->p4());
1612  if(dR<0.3){
1613  if (verbose) cout<<"Jet: ("<<all_baby_jets[ijet].pt()<<", "<<all_baby_jets[ijet].eta()<<", "<<all_baby_jets[ijet].phi()
1614  <<"), MC: ("<<mc.daughter(idau)->pt()<<", "<<mc.daughter(idau)->eta()<<", "<<mc.daughter(idau)->phi()<<"), ID "<<mc.daughter(idau)->pdgId()<<". dR "<<dR <<endl;
1615  matched = true;
1616  break;
1617  }
1618  }
1619  } // Loop over MC particles
1620  if(!matched) {
1621  nisr++;
1622  baby.jets_isisr()[ijet] = true;
1623  }
1624  } // Loop over jets
1625  baby.nisr() = nisr;
1626 }
1627 
1628 void bmaker_full::writeMC(edm::Handle<reco::GenParticleCollection> genParticles,
1629  vCands &all_mus, vCands &all_els, vCands &photons){
1630  LVector isr_p4;
1631  float metw_tru_x(0.), metw_tru_y(0.);
1632  float lep_tru_pt(0.), lep_tru_phi(0.);
1633  baby.ntruleps()=0; baby.ntrumus()=0; baby.ntruels()=0; baby.ntrutaush()=0; baby.ntrutausl()=0;
1634  baby.nleps_tm()=0;
1635  baby.fromGS()=false;
1636  baby.m_tt()=0; baby.top_pt()=0; baby.antitop_pt()=0;
1637  vector<float> top_pt;
1638  int topIndex=-1;
1639  int antitopIndex=-1;
1640  const size_t bsmid(1000000);
1641 
1643  vector<pair<int, const reco::GenParticle *> > indices;
1644  int Nsaved=0;
1645  for (size_t imc(0); imc < genParticles->size(); imc++) {
1646  const reco::GenParticle &mc = (*genParticles)[imc];
1647  size_t id = abs(mc.pdgId());
1648  bool isLast = mcTool->isLast(mc, id);
1649  bool isFirst = true;
1650  const reco::GenParticle *mcMom = static_cast<const reco::GenParticle *>(mc.mother());
1651  if(mcMom){
1652  if(mcMom->pdgId() == mc.pdgId()) isFirst = false;
1653  }
1654 
1655  const reco::GenParticle *mom = nullptr;
1656  size_t momid = abs(mcTool->mom(mc, mom));
1657  bool isTop(id==6);
1658  bool isNewPhysics(id>=bsmid);
1659  bool isGluino(id==1000021);
1660  bool isNeutralino(id==1000023 || id==1000025);
1661  bool isZ(id==23);
1662  bool isW(id==24);
1663  bool isH(id==25);
1664  bool bTopOrBSM(id==5 && (momid==6 || momid>=bsmid));
1665  bool nuFromZ((id==12 || id==14 || id==16) && momid==23);
1666  bool eFromTopZ(id==11 && (momid==24 || momid==23));
1667  bool muFromTopZ(id==13 && (momid==24 || momid==23));
1668  bool tauFromTopZ(id==15 && (momid==24 || momid==23));
1669  bool fromTau(mcTool->fromTau(mc));
1670  bool fromWOrWTau(mcTool->fromWOrWTau(mc));
1671  bool chgPionFromTau(id==211 && momid==15 && fromWOrWTau);
1672  bool eFromTopZtau(id==11 && (momid==24 || momid==23 || fromTau ));
1673  bool muFromTopZtau(id==13 && (momid==24 || momid==23 || fromTau ));
1674  bool fromHiggs(momid==25);
1675  bool fromTop(momid==6);
1676  bool mg_me = mc.status()==23;
1677 
1678  //Identify if c or b is from gluon splitting. Save first pair after splitting (arbitrary choice)
1679  bool from_gs = false;
1680  if(isFirst && (id==4 || id==5) && mcTool->isFromGSP(dynamic_cast<const reco::Candidate*>(&mc))) from_gs = true;
1681  if(from_gs && !baby.fromGS()) baby.fromGS() = true;
1682 
1683  //isLast is not safe, since 21 -> -4,4,21 fails isLast but yields 4s that pass isFromGSP
1684  if(id==21){
1685  for(size_t idau(0); idau < mc.numberOfDaughters(); idau++) {
1686  int dauid(abs(mc.daughter(idau)->pdgId()));
1687  if(dauid == 4 || dauid==5) {
1688  if(mcTool->isFromGSP(mc.daughter(idau))){
1689  from_gs=true;
1690  break;
1691  }
1692  }
1693  }
1694  }
1695 
1696 
1698  if((isLast && (isTop || isNewPhysics || isZ || isH))
1699  || (isFirst && (fromHiggs || isW || bTopOrBSM || eFromTopZ || muFromTopZ
1700  || tauFromTopZ || nuFromZ || fromWOrWTau || fromTop))
1701  || mg_me || from_gs){
1702 
1703  //Find jet associated with interesting particles
1704  //Match using anti-kt, R=0.4 metric between mc vector and jets
1705  //Require dR inside 0.4 as well (This represents a jet being closed and removed from list)
1706  double mind=999.;
1707  int minjetidx=-1;
1708  for(size_t ind(0); ind < baby.jets_pt().size(); ind++) {
1709  double dr = dR(baby.jets_phi()[ind],mc.phi(),baby.jets_eta()[ind],mc.eta());
1710  double dij = pow(max(baby.jets_pt()[ind],static_cast<float>(mc.pt())),-2)*dr*pow(0.4,-2);
1711  if(dij<mind && dr < 0.4){
1712  mind=dij;
1713  minjetidx=ind;
1714  }
1715  }
1716 
1717  baby.mc_status().push_back(mc.status());
1718  baby.mc_id().push_back(mc.pdgId());
1719  baby.mc_pt().push_back(mc.pt());
1720  baby.mc_eta().push_back(mc.eta());
1721  baby.mc_phi().push_back(mc.phi());
1722  baby.mc_mass().push_back(mc.mass());
1723  baby.mc_mom().push_back(mcTool->mom(mc,mom));
1724  baby.mc_gs().push_back(from_gs);
1725  baby.mc_gs_dau_jetmatch().push_back(-1); //Filled after loop
1726  baby.mc_gs_dau_dr().push_back(-1.); //Filled after loop
1727  baby.mc_jetidx().push_back(minjetidx);
1728  baby.mc_num_dau().push_back(mc.numberOfDaughters());
1729 
1730  // Saving mother index
1731  baby.mc_momidx().push_back(mcTool->getMomIndex(mc,indices));
1732  indices.push_back(pair<int, const reco::GenParticle *>(Nsaved, &mc));
1733  Nsaved++;
1734 
1735 
1736 
1737  //Store some convenient info in jets collection
1738  if(minjetidx>=0){
1739  if(id==5) baby.jets_ntrub()[minjetidx]++;
1740  if(id==4) baby.jets_ntruc()[minjetidx]++;
1741  if(from_gs && (id==4||id==5)) baby.jets_gs_index()[minjetidx]= baby.mc_momidx().back();
1742  }
1743 
1744 
1745  }
1746 
1747 
1748  if(isLast){
1749  if(isTop) mc.pdgId()>0 ? topIndex=imc : antitopIndex=imc;
1750 
1752  if((isTop && (outname.Contains("TTJets") || outname.Contains("TT_") || outname.Contains("TTTo")))
1753  || (isGluino && (outname.Contains("SMS") || outname.Contains("RPV")))
1754  || (isZ && outname.Contains("DY"))) isr_p4 -= mc.p4();
1755 
1756  if(isNeutralino && outname.Contains("TChiHH")) isr_p4 -= mc.p4();
1757 
1758  if(isTop && (outname.Contains("TTJets") || outname.Contains("TT_") || outname.Contains("TTTo"))){
1759  top_pt.push_back(mc.pt());
1760  }
1761 
1763  if(muFromTopZ) baby.ntrumus()++;
1764  if(eFromTopZ) baby.ntruels()++;
1765  if(tauFromTopZ){
1766  const reco::GenParticle *tauDaughter(0);
1767  if(mcTool->decaysTo(mc, 11, tauDaughter) || mcTool->decaysTo(mc, 13, tauDaughter)){
1768  baby.ntrutausl()++;
1769  } else baby.ntrutaush()++;
1770  }
1771  }
1772 
1773  //Finding lost leptons
1774  const float relptThreshold(0.3), drThreshold(0.1);
1775  float relptThres(2.), drThres(0.15);
1776  //if(mcTool->fromWOrWTau(mc) && (laste||lastmu||lasttau|| (mcTool->fromTau(mc) && abs(mc.pdgId())==211))){
1777  if(eFromTopZtau|| muFromTopZtau || chgPionFromTau){
1778  double mindr(999.);
1779  int minind(-1);
1780 
1781  for(size_t ind(0); ind < baby.leps_pt().size(); ind++) {
1782  double dr = dR(baby.leps_phi()[ind],mc.phi(),baby.leps_eta()[ind],mc.eta());
1783  // cout<<"mc lep eta phi = "<<mc.eta()<<" "<<mc.phi()<<", reco lep eta phi "<<baby.leps_eta()[ind]<<" "<<baby.leps_phi()[ind]<<", dr = "<<dr<<endl;
1784  double drelpt(fabs((baby.leps_pt()[ind] - mc.pt())/mc.pt()));
1785  if(dr > drThres || drelpt > relptThres) continue;
1786  if(dr < mindr){
1787  mindr = dr;
1788  minind = ind;
1789  }
1790  }
1791 
1792 
1793  if(minind<0){ //Lepton is lost
1794  //Try to match to veto track collection
1795  double min_dr(999.);
1796  int minindex(-1);
1797  for(size_t ind(0); ind < baby.tks_pt().size(); ind++) {
1798  double dr = dR(baby.tks_phi()[ind],mc.phi(),baby.tks_eta()[ind],mc.eta());
1799  double drelpt(fabs((baby.tks_pt()[ind] - mc.pt())/mc.pt()));
1800  if(dr > drThreshold || drelpt > relptThres) continue; //use loose pt comparison in case lost lepton is very poorly measured
1801  if(dr < min_dr){
1802  min_dr = dr;
1803  minindex = ind;
1804  }
1805  } // Loop over tks
1806  if(minindex >= 0) {
1807  baby.tks_tm()[minindex] = true;
1808  }
1809 
1810 
1811  }
1812  }
1814  //const float relptThreshold(0.3), drThreshold(0.1); needed earlier
1815  if(id==11 && fromWOrWTau){
1816  double mindr(999.);
1817  int minind(-1);
1818  for(size_t ind(0); ind < all_els.size(); ind++) {
1819  double dr(deltaR(mc, *(all_els[ind])));
1820  double drelpt(fabs((all_els[ind]->pt() - mc.pt())/mc.pt()));
1821  if(dr > drThreshold || drelpt > relptThreshold) continue;
1822  if(dr < mindr){
1823  mindr = dr;
1824  minind = ind;
1825  }
1826  } // Loop over all_els
1827  if(minind >= 0) {
1828  baby.els_tm()[minind] = true;
1829  if(baby.els_sig()[minind]) baby.nleps_tm()++;
1830  }
1831  if(lep_tru_pt < mc.pt()){
1832  lep_tru_pt = mc.pt();
1833  lep_tru_phi = mc.phi();
1834  } // Lepton pt to find mt_tru
1835  } // If it is an electron
1836  if(id==13 && fromWOrWTau){
1837  double mindr(999.);
1838  int minind(-1);
1839  for(size_t ind(0); ind < all_mus.size(); ind++) {
1840  double dr(deltaR(mc, *(all_mus[ind])));
1841  double drelpt(fabs((all_mus[ind]->pt() - mc.pt())/mc.pt()));
1842  if(dr > drThreshold || drelpt > relptThreshold) continue;
1843  if(dr < mindr){
1844  mindr = dr;
1845  minind = ind;
1846  }
1847  } // Loop over all_mus
1848  if(minind >= 0) {
1849  baby.mus_tm()[minind] = true;
1850  if(baby.mus_sig()[minind]) baby.nleps_tm()++;
1851  }
1852  if(lep_tru_pt < mc.pt()){
1853  lep_tru_pt = mc.pt();
1854  lep_tru_phi = mc.phi();
1855  } // Lepton pt to find mt_tru
1856  } // If it is a muon
1857  if(id==22){
1858  double mindr(999.);
1859  int minind(-1);
1860  for(size_t ind(0); ind < photons.size(); ind++) {
1861  double dr(deltaR(mc, *(photons[ind])));
1862  double drelpt(fabs((photons[ind]->pt() - mc.pt())/mc.pt()));
1863  if(dr > drThreshold || drelpt > relptThreshold) continue;
1864  if(dr < mindr){
1865  mindr = dr;
1866  minind = ind;
1867  }
1868  } // Loop over photons
1869  if(minind >= 0) baby.ph_tm()[minind] = true;
1870  } // If it is a photon
1871 
1873  if((id==12 || id==14 || id==16 || id==18 || id==1000012 || id==1000014 || id==1000016
1874  || id==1000022 || id==1000023 || id==1000025 || id==1000035 || id==1000039) &&
1875  id != momid){ // neutrinos decay to themselves
1876  if(fromWOrWTau) {
1877  metw_tru_x += mc.px();
1878  metw_tru_y += mc.py();
1879  }
1880  } // If undetected neutral particle
1881 
1882  } // Loop over genParticles
1883  // calculate invariant mass of ttbar pair
1884  if(topIndex>=0 && antitopIndex>=0) {
1885  reco::Candidate::LorentzVector topP4 = genParticles->at(topIndex).p4();
1886  reco::Candidate::LorentzVector antitopP4 = genParticles->at(antitopIndex).p4();
1887  reco::Candidate::LorentzVector ttbarP4 = topP4+antitopP4;
1888  baby.m_tt()=ttbarP4.mass();
1889  baby.top_pt()=topP4.pt();
1890  baby.antitop_pt()=antitopP4.pt();
1891  }
1892 
1893  //Loop over gluon splitting gluons to mark truth info in bb pairs
1894  for(size_t ind(0); ind < baby.mc_pt().size(); ind++) {
1895  if(baby.mc_gs()[ind] && baby.mc_id()[ind]==21){
1896 
1897  //Found gluon that splits to bb or cc. Now find daughter quarks, and store their index
1898  vector<int> q_pair_idx;
1899  for(size_t indd(0); indd < baby.mc_pt().size(); indd++) {
1900  if(static_cast<size_t>(baby.mc_momidx()[indd])==ind) q_pair_idx.push_back(indd);
1901  }
1902  if(q_pair_idx.size()!=2) continue;
1903 
1904  //Store info about daughters in mc collection, associated with gluon
1905  //store dR of gluon daughters
1906  baby.mc_gs_dau_dr()[ind]= dR(baby.mc_phi()[q_pair_idx[0]], baby.mc_phi()[q_pair_idx[1]],baby.mc_eta()[q_pair_idx[0]], baby.mc_eta()[q_pair_idx[1]]);
1907 
1908 
1909  // mc_gs_dau_jetmatch convention
1910  // -1: default, remains this way if number of daughters is not 2
1911  // 0: both quarks fail jet matching
1912  // 1: 1 quark is matched to a jet, 1 quark is lost
1913  // 2: both quarks are matched to the same jet
1914  // 3: both quarks are matched to different jets
1915 
1916  if(baby.mc_jetidx()[q_pair_idx[0]]<0 && baby.mc_jetidx()[q_pair_idx[1]]<0) baby.mc_gs_dau_jetmatch()[ind]=0;
1917  else if(baby.mc_jetidx()[q_pair_idx[0]]<0 || baby.mc_jetidx()[q_pair_idx[1]]<0) baby.mc_gs_dau_jetmatch()[ind]=1;
1918  else if(baby.mc_jetidx()[q_pair_idx[0]] == baby.mc_jetidx()[q_pair_idx[1]] ) baby.mc_gs_dau_jetmatch()[ind]=2;
1919  else baby.mc_gs_dau_jetmatch()[ind]=3;
1920 
1921 
1922 
1923  //loop over bb pairs, to fill truth info
1924  for(size_t indbb(0);indbb<baby.dr_bb().size();indbb++){
1925  //if quark1 matches jet1 and q2 matches jet2, or vice versa, fill bb_idx and bb_gs_flavor
1926  if(((baby.mc_jetidx()[q_pair_idx[0]]==baby.bb_jet_idx1()[indbb]) && (baby.mc_jetidx()[q_pair_idx[1]]==baby.bb_jet_idx2()[indbb])) ||((baby.mc_jetidx()[q_pair_idx[0]]==baby.bb_jet_idx2()[indbb]) && (baby.mc_jetidx()[q_pair_idx[1]]==baby.bb_jet_idx1()[indbb]))){
1927  baby.bb_gs_idx()[indbb]=ind; //using gluon index
1928  baby.bb_gs_flavor()[indbb]=baby.mc_id()[q_pair_idx[0]];
1929  }
1930  }
1931 
1932  //loop over deep bb pairs, to fill truth info
1933  for(size_t indbb(0);indbb<baby.dr_bb_deep().size();indbb++){
1934  //if quark1 matches jet1 and q2 matches jet2, or vice versa, fill bb_idx and bb_gs_flavor
1935  if(((baby.mc_jetidx()[q_pair_idx[0]]==baby.bb_jet_idx1_deep()[indbb]) && (baby.mc_jetidx()[q_pair_idx[1]]==baby.bb_jet_idx2_deep()[indbb])) ||((baby.mc_jetidx()[q_pair_idx[0]]==baby.bb_jet_idx2_deep()[indbb]) && (baby.mc_jetidx()[q_pair_idx[1]]==baby.bb_jet_idx1_deep()[indbb]))){
1936  baby.bb_gs_idx_deep()[indbb]=ind; //using gluon index
1937  baby.bb_gs_flavor_deep()[indbb]=baby.mc_id()[q_pair_idx[0]];
1938  }
1939  }
1940  }
1941  }
1942 
1943  baby.ntruleps() = baby.ntrumus()+baby.ntruels()+baby.ntrutaush()+baby.ntrutausl();
1944  baby.isr_tru_pt() = isr_p4.pt();
1945  baby.isr_tru_eta() = isr_p4.eta();
1946  baby.isr_tru_phi() = isr_p4.phi();
1947 
1948  if((outname.Contains("TTJets") || outname.Contains("TT_") || outname.Contains("TTTo")) && top_pt.size() == 2) baby.w_toppt() = weightTool->topPtWeight(top_pt.at(0),top_pt.at(1));
1949  else baby.w_toppt() = 1.;
1950 
1951  baby.met_tru_nuw() = hypot(metw_tru_x, metw_tru_y);
1952  baby.met_tru_nuw_phi() = atan2(metw_tru_y, metw_tru_x);
1953 
1954  baby.mt_tru() = getMT(baby.met_tru(), baby.met_tru_phi(), lep_tru_pt, lep_tru_phi);
1955  baby.mt_tru_nuw() = getMT(baby.met_tru_nuw(), baby.met_tru_nuw_phi(), lep_tru_pt, lep_tru_phi);
1956 } // writeMC
1957 
1958 // Finds the jet that minimizes the MET when a variation is performed
1959 void bmaker_full::rebalancedMET( double& minMET, double& minMETPhi)
1960 {
1961  for(unsigned int iJet=0; iJet<baby.jets_pt().size(); iJet++) {
1962  // calculate best rescaling factor for this jet
1963  double rescalingFactor=calculateRescalingFactor(iJet);
1964  double newMETPhi=0;
1965  double newMET=calculateRebalancedMET(iJet, rescalingFactor, newMETPhi);
1966  if(newMET<minMET) {
1967  minMET=newMET;
1968  minMETPhi=newMETPhi;
1969  }
1970  }
1971 }
1972 
1973 // calculate a rebalancing of the jet momentum that minimizes MET
1974 double bmaker_full::calculateRescalingFactor(unsigned int jetIdx)
1975 {
1976 
1977  // don't allow jet pt to be scaled by more than this factor
1978  const double scaleCutoff=1;
1979 
1980  TVector3 jet, metVector;
1981  jet.SetPtEtaPhi(baby.jets_pt().at(jetIdx), baby.jets_eta().at(jetIdx), baby.jets_phi().at(jetIdx));
1982  metVector.SetPtEtaPhi(baby.met(), 0, baby.met_phi());
1983 
1984  double denominator = -jet.Px()*jet.Px()-jet.Py()*jet.Py();
1985  double numerator = jet.Px()*metVector.Px()+jet.Py()*metVector.Py();
1986 
1987  double rescalingFactor=1e6;
1988  if(denominator!=0) rescalingFactor = numerator/denominator;
1989  if(fabs(rescalingFactor)>scaleCutoff) rescalingFactor=scaleCutoff*rescalingFactor/fabs(rescalingFactor);
1990  // the resolution tail is on the _low_ side, not the high side
1991  // so we always need to subtract pT
1992  if(rescalingFactor>0) rescalingFactor=0;
1993 
1994  return rescalingFactor;
1995 }
1996 
1997 double bmaker_full::calculateRebalancedMET(unsigned int jetIdx, double mu, double& METPhi)
1998 {
1999  TVector3 jet, metVector;
2000  jet.SetPtEtaPhi(baby.jets_pt().at(jetIdx), baby.jets_eta().at(jetIdx), baby.jets_phi().at(jetIdx));
2001  metVector.SetPtEtaPhi(baby.met(), 0, baby.met_phi());
2002 
2003  double sumPx = metVector.Px()+mu*jet.Px();
2004  double sumPy = metVector.Py()+mu*jet.Py();
2005 
2006  METPhi=atan(sumPy/sumPx);
2007 
2008  return sqrt(sumPx*sumPx+sumPy*sumPy);
2009 }
2010 
2011 void bmaker_full::writeWeights(const vCands &sig_leps, edm::Handle<GenEventInfoProduct> &gen_event_info,
2012  edm::Handle<LHEEventProduct> &lhe_info){
2013  if (debug) cout<<"INFO: Filling weights..."<<endl;
2014 
2015  // Initializing weights
2016  if(isData) {
2017  baby.eff_trig() = baby.w_btag() = baby.w_btag_loose() = baby.w_btag_tight() = baby.w_pu() = baby.w_pu() = baby.w_lep() = baby.w_fs_lep() = baby.w_toppt() = 1.;
2018  baby.eff_jetid() = baby.w_lumi() = baby.weight() = baby.weight_rpv() = 1.;
2019  return;
2020  }
2021 
2022  // Luminosity weight
2023  const float luminosity = 1000., fneg(xsec::fractionNegWeights(outname));
2024  baby.w_lumi() = xsec*luminosity / (static_cast<double>(nevents_sample)*(1-2*fneg));
2025  if (gen_event_info->weight() < 0) baby.w_lumi() *= -1.;
2026 
2027  // Pile-up weight
2028  baby.w_pu() = weightTool->pileupWeight(baby.ntrupv_mean(),0);
2029  baby.sys_pu().resize(2, 1.);
2030  baby.sys_pu()[0] = weightTool->pileupWeight(baby.ntrupv_mean(), 1);
2031  baby.sys_pu()[1] = weightTool->pileupWeight(baby.ntrupv_mean(), -1);
2032 
2033  // Lepton SFs
2034  pair<double, double> sf_err = lepTool->getScaleFactor(sig_leps);
2035  double sf = sf_err.first;
2036  double unc = sf_err.second;
2037  baby.w_lep() = sf;
2038  baby.sys_lep().resize(2, 1.);
2039  baby.sys_lep().at(0) = sf+unc;
2040  baby.sys_lep().at(1) = sf-unc;
2041 
2042  // Lepton SFs in FastSim
2043  baby.sys_fs_lep().resize(2, 1.);
2044  if(isFastSim){
2045  pair<double, double> sf_err_fs = lepTool->getScaleFactorFs(sig_leps);
2046  double sf_fs = sf_err_fs.first;
2047  double unc_fs = sf_err_fs.second;
2048  baby.w_fs_lep() = sf_fs;
2049  baby.sys_fs_lep()[0] = sf_fs+unc_fs;
2050  baby.sys_fs_lep()[1] = sf_fs-unc_fs;
2051  } else baby.w_fs_lep() = 1.;
2052 
2053  // VVVL trigger efficiency
2054  baby.eff_trig() = weightTool->triggerEfficiency(baby.nmus(), baby.nels(), baby.met(), baby.sys_trig());
2055 
2056  // In FastSim the JetID is broken, so we just apply 0.99 +- 0.01
2057  if(isFastSim) baby.eff_jetid() = 0.99;
2058  else baby.eff_jetid() = 1.;
2060  // w_btag calculated in writeJets
2061  // w_toppt and sys_isr calculated in writeMC
2062  baby.weight() = baby.w_lumi() * baby.w_lep() * baby.w_fs_lep() * baby.w_btag()
2063  * baby.eff_jetid();
2064  baby.weight_rpv() = baby.w_lumi() * baby.w_lep() * baby.w_fs_lep() * baby.w_btag()
2065  * baby.w_pu() * baby.eff_jetid();
2066 
2068  if(lhe_info.isValid()) weightTool->getTheoryWeights(lhe_info);
2069  // Renormalization and Factorization scales
2070  baby.sys_mur().push_back(weightTool->theoryWeight(weight_tools::muRup));
2071  baby.sys_mur().push_back(weightTool->theoryWeight(weight_tools::muRdown));
2072  baby.sys_muf().push_back(weightTool->theoryWeight(weight_tools::muFup));
2073  baby.sys_muf().push_back(weightTool->theoryWeight(weight_tools::muFdown));
2074  baby.sys_murf().push_back(weightTool->theoryWeight(weight_tools::muRup_muFup));
2075  baby.sys_murf().push_back(weightTool->theoryWeight(weight_tools::muRdown_muFdown));
2076  // PDF variations
2077  weightTool->getPDFWeights(baby.sys_pdf(), baby.w_pdf());
2078 
2080  baby.w_isr() = 1.; baby.sys_isr().resize(2,1.);
2081  if (baby.type()/1000==1 || (baby.type()/1000>=100 && !outname.Contains("TChi"))){
2082  const float isr_norm_tt = 1.117;
2083  float isr_wgt = -999.;
2084  if (baby.nisr()==0) isr_wgt = 1.;
2085  else if (baby.nisr()==1) isr_wgt = 0.882;
2086  else if (baby.nisr()==2) isr_wgt = 0.792;
2087  else if (baby.nisr()==3) isr_wgt = 0.702;
2088  else if (baby.nisr()==4) isr_wgt = 0.648;
2089  else if (baby.nisr()==5) isr_wgt = 0.601;
2090  else if (baby.nisr()>=6) isr_wgt = 0.515;
2091  baby.w_isr() = isr_wgt*isr_norm_tt;
2092  //assign relative unc = 50% of the deviation from flat
2093  float absolute_unc = (1-isr_wgt)/2.;
2094  baby.sys_isr()[0] = isr_norm_tt*(isr_wgt+absolute_unc);
2095  baby.sys_isr()[1] = isr_norm_tt*(isr_wgt-absolute_unc);
2096  } else if (outname.Contains("TChi")) {
2097  float isr_wgt = 1.;
2098  if (baby.isr_tru_pt()<=50) isr_wgt = 1.;
2099  else if (baby.isr_tru_pt()<=100) isr_wgt = 1.052;
2100  else if (baby.isr_tru_pt()<=150) isr_wgt = 1.179;
2101  else if (baby.isr_tru_pt()<=200) isr_wgt = 1.150;
2102  else if (baby.isr_tru_pt()<=300) isr_wgt = 1.057;
2103  else if (baby.isr_tru_pt()<=400) isr_wgt = 1.000;
2104  else if (baby.isr_tru_pt()<=600) isr_wgt = 0.912;
2105  else isr_wgt = 0.783;
2106  baby.w_isr() = isr_wgt;
2107  //assign relative unc = 100% of the deviation from flat
2108  if (isr_wgt>1) baby.sys_isr()[0] = 1+2*(isr_wgt-1);
2109  else baby.sys_isr()[0] = 1-2*(1-isr_wgt);
2110  baby.sys_isr()[1] = 1.;
2111  }
2112 }
2113 
2114 /*
2115  _____ _ _
2116  / __ \ | | | |
2117  | / \/ ___ _ __ ___| |_ _ __ _ _ ___| |_ ___ _ __ ___
2118  | | / _ \| '_ \/ __| __| '__| | | |/ __| __/ _ \| '__/ __|
2119  | \__/\ (_) | | | \__ \ |_| | | |_| | (__| || (_) | | \__ \
2120  \____/\___/|_| |_|___/\__|_| \__,_|\___|\__\___/|_| |___/
2121 */
2122 //Constructor (this line searchable)
2123 bmaker_full::bmaker_full(const edm::ParameterSet& iConfig):
2124  outname(TString(iConfig.getParameter<string>("outputFile"))),
2125  inputfiles(iConfig.getParameter<vector<string> >("inputFiles")),
2126  jsonfile(iConfig.getParameter<string>("json")),
2127  condor_subtime(iConfig.getParameter<string>("condor_subtime")),
2128  jec_label(iConfig.getParameter<string>("jec")),
2129  met_label(iConfig.getParameter<edm::InputTag>("met")),
2130  met_nohf_label(iConfig.getParameter<edm::InputTag>("met_nohf")),
2131  jets_label(iConfig.getParameter<edm::InputTag>("jets")),
2132  nevents_sample(iConfig.getParameter<unsigned int>("nEventsSample")),
2133  nevents(0),
2134  doMetRebalancing(iConfig.getParameter<bool>("doMetRebalancing")),
2135  addBTagWeights(iConfig.getParameter<bool>("addBTagWeights")),
2136  isFastSim(iConfig.getParameter<bool>("isFastSim")),
2137  doSystematics(iConfig.getParameter<bool>("doSystematics")),
2138  debug(iConfig.getParameter<bool>("debugMode")),
2139 
2140  // Initialize tokens
2141  tok_trigResults_hlt_(consumes<edm::TriggerResults>(edm::InputTag("TriggerResults","","HLT"))),
2142  tok_patTrig_(consumes<pat::PackedTriggerPrescales>(edm::InputTag("patTrigger"))),
2143  tok_primVertex_(consumes<reco::VertexCollection>(edm::InputTag("offlineSlimmedPrimaryVertices"))),
2144  tok_addPileup_(consumes<std::vector< PileupSummaryInfo > >(edm::InputTag("addPileupInfo"))),
2145  tok_slimAddPileup_(consumes<std::vector< PileupSummaryInfo > >(edm::InputTag("slimmedAddPileupInfo"))),
2146  tok_packPFCands_(consumes<pat::PackedCandidateCollection>(edm::InputTag("packedPFCandidates"))),
2147  tok_rhoFastJet_centralNeutral_(consumes<double>(edm::InputTag("fixedGridRhoFastjetCentralNeutral"))),
2148  tok_muons_(consumes<pat::MuonCollection>(edm::InputTag("slimmedMuons"))),
2149  tok_electrons_(consumes<pat::ElectronCollection>(edm::InputTag("slimmedElectrons"))),
2150  tok_rhoFastJet_all_(consumes<double>(edm::InputTag("fixedGridRhoFastjetAll"))),
2151  tok_offBeamSpot_(consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
2152  tok_photons_(consumes<pat::PhotonCollection>(edm::InputTag("slimmedPhotons"))),
2153  tok_reducedEgamma_conver_(consumes<vector<reco::Conversion> >(edm::InputTag("reducedEgamma","reducedConversions"))),
2154  tok_jets_(consumes<pat::JetCollection>(jets_label)),
2155  tok_genJets_(consumes<edm::View<reco::GenJet> >(edm::InputTag("slimmedGenJets"))),
2156  tok_met_(consumes<pat::METCollection>(met_label)),
2157  tok_met_noHF_(consumes<pat::METCollection>(met_nohf_label)),
2158  tok_met_Puppi_(consumes<pat::METCollection>(edm::InputTag("slimmedMETsPuppi"))),
2159  tok_HBHENoiseFilter_(consumes<bool>(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"))),
2160  tok_HBHEIsoNoiseFilter_(consumes<bool>(edm::InputTag("HBHENoiseFilterResultProducer","HBHEIsoNoiseFilterResult"))),
2161  tok_trigResults_reco_(consumes<edm::TriggerResults>(edm::InputTag("TriggerResults","","RECO"))),
2162  tok_trigResults_pat_(consumes<edm::TriggerResults>(edm::InputTag("TriggerResults","","PAT"))),
2163  tok_selectedPatTrig_(consumes<pat::TriggerObjectStandAloneCollection>(edm::InputTag("slimmedPatTrigger"))),
2164  tok_pruneGenPart_(consumes<reco::GenParticleCollection>(edm::InputTag("prunedGenParticles"))),
2165  tok_extLHEProducer_(consumes<LHEEventProduct>(edm::InputTag("externalLHEProducer"))),
2166  tok_source_(consumes<LHEEventProduct>(edm::InputTag("source"))),
2167  tok_generator_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
2168  tok_genlumiheader_(consumes<GenLumiInfoHeader,edm::InLumi>(edm::InputTag("generator")))
2169 /* DAK8
2170  tok_deepJetToken_(consumes<edm::View<pat::Jet> >(edm::InputTag("slimmedJetsAK8")))
2171 */
2172 {
2173  time(&startTime);
2174 
2175 
2176  lepTool = new lepton_tools();
2178  photonTool = new photon_tools();
2179  mcTool = new mc_tools();
2180  weightTool = new weight_tools();
2181  eventTool = new event_tools(outname);
2182 /*
2183  string deepJetDataPath = "NNKit/data/ak8";
2184  deepJetR = 0.8;
2185  auto cc = consumesCollector();
2186  fatjetNN_ = std::make_unique<FatJetNN>(iConfig, cc, deepJetR);
2187  // load json for input variable transformation
2188  fatjetNN_->load_json(edm::FileInPath(deepJetDataPath+"/full/preprocessing.json").fullPath());
2189  // load DNN model and parameter files
2190  fatjetNN_->load_model(edm::FileInPath(deepJetDataPath+"/full/resnet-symbol.json").fullPath(),
2191  edm::FileInPath(deepJetDataPath+"/full/resnet.params").fullPath());
2192 
2193  decorrNN_ = std::make_unique<FatJetNN>(iConfig, cc, deepJetR);
2194  // load json for input variable transformation
2195  decorrNN_->load_json(edm::FileInPath(deepJetDataPath+"/decorrelated/preprocessing.json").fullPath());
2196  // load DNN model and parameter files
2197  decorrNN_->load_model(edm::FileInPath(deepJetDataPath+"/decorrelated/resnet-symbol.json").fullPath(),
2198  edm::FileInPath(deepJetDataPath+"/decorrelated/resnet.params").fullPath());
2199 */
2200  outfile = new TFile(outname, "recreate");
2201  outfile->cd();
2202  baby.tree_.SetDirectory(outfile);
2203 
2205  // if(xsec<=0) {
2206  // cout<<"BABYMAKER: Cross section not found, aborting"<<endl<<endl;
2207  // exit(1);
2208  // }
2209 
2210  trig_name = vector<TString>();
2211  if(outname.Contains("Run201")){ // Would like to define isData, but need iEvent?
2212  if(outname.Contains("Run2017") || outname.Contains("Run2018")){
2213  trig_name.push_back("HLT_PFHT500_PFMET100_PFMHT100_IDTight_v"); // 0
2214  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT450_PFMET50_v"); // 1
2215  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT600_v"); // 2
2216  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT450_v"); // 3
2217  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT450_CaloBTagCSV_4p5_v"); // 4
2218  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT450_PFMET50_v"); // 5
2219  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT600_v"); // 6
2220  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT450_v"); // 7
2221  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT450_CaloBTagCSV_4p5_v"); // 8
2222  trig_name.push_back("HLT_PFMET120_PFMHT120_IDTight_v"); // 9
2223  trig_name.push_back("HLT_PFMET120_PFMHT120_IDTight_PFHT60_v"); // 10
2224 
2225  trig_name.push_back("HLT_PFMET120_PFMHT120_IDTight_HFCleaned_v"); // 11
2226  trig_name.push_back("HLT_PFMET120_PFMHT120_IDTight_PFHT60_HFCleaned_v"); // 12
2227  trig_name.push_back("HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_PFHT60_v"); // 13
2228  trig_name.push_back("HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_HFCleaned_v"); // 14
2229  trig_name.push_back("HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_v"); // 15
2230  trig_name.push_back("HLT_PFHT700_PFMET85_PFMHT85_IDTight_v"); // 16
2231  trig_name.push_back("HLT_PFHT800_PFMET75_PFMHT75_IDTight_v"); // 17
2232  trig_name.push_back("HLT_PFHT1050_v"); // 18
2233  trig_name.push_back("HLT_IsoMu24_v"); // 19
2234  trig_name.push_back("HLT_IsoMu27_v"); // 20
2235 
2236  trig_name.push_back("HLT_Mu50_v"); // 21
2237  trig_name.push_back("HLT_Mu50_IsoVVVL_PFHT450_v"); // 22
2238  trig_name.push_back("HLT_Ele35_WPTight_Gsf_v"); // 23
2239  trig_name.push_back("HLT_Ele115_CaloIdVT_GsfTrkIdT_v"); // 24
2240  trig_name.push_back("HLT_Ele300_CaloIdVT_GsfTrkIdT_v"); // 25
2241  trig_name.push_back("HLT_Ele27_WPTight_Gsf_v"); // 26
2242  trig_name.push_back("HLT_Ele28_eta2p1_WPTight_Gsf_HT150_v"); // 27
2243  trig_name.push_back("HLT_Ele30_eta2p1_WPTight_Gsf_CentralPFJet35_EleCleaned_v"); // 28
2244  trig_name.push_back("HLT_Ele38_WPTight_Gsf_v"); // 29
2245  trig_name.push_back("HLT_Ele50_IsoVVVL_PFHT450_v"); // 30
2246 
2247  trig_name.push_back("HLT_Photon200_v"); // 31
2248  trig_name.push_back("HLT_Photon300_NoHE_v"); // 32
2249  trig_name.push_back("HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass3p8_v"); // 33
2250  trig_name.push_back("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v"); // 34
2251  trig_name.push_back("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_v"); // 35
2252  trig_name.push_back("HLT_DoubleEle33_CaloIdL_MW_v"); // 36
2253  trig_name.push_back("HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_DZ_PFHT350_v"); // 37
2254  trig_name.push_back("HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT350_v"); // 38
2255  trig_name.push_back("HLT_PFHT180_v"); // 39
2256  trig_name.push_back("HLT_PFHT250_v"); // 40
2257 
2258  trig_name.push_back("HLT_PFHT350_v"); // 41
2259  trig_name.push_back("HLT_PFHT510_v"); // 42
2260  trig_name.push_back("HLT_PFHT680_v"); // 43
2261  trig_name.push_back("HLT_PFHT890_v"); // 44
2262  trig_name.push_back("HLT_PFJet40_v"); // 45
2263  trig_name.push_back("HLT_PFJet140_v"); // 46
2264  trig_name.push_back("HLT_PFJet260_v"); // 47
2265  trig_name.push_back("HLT_PFJet500_v"); // 48
2266  trig_name.push_back("HLT_AK8PFJet500_v"); // 49
2267  trig_name.push_back("HLT_AK8PFJet360_TrimMass30_v"); // 50
2268 
2269  trig_name.push_back("HLT_PFMET250_HBHECleaned_v"); // 51
2270  } else if(outname.Contains("Run2016")){
2271  trig_name.push_back("HLT_PFHT300_PFMET100_v"); // 0
2272  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT350_PFMET50_v"); // 1
2273  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT600_v"); // 2
2274  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT400_v"); // 3
2275  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT350_v"); // 4
2276  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT350_PFMET50_v"); // 5
2277  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT600_v"); // 6
2278  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT400_v"); // 7
2279  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT350_v"); // 8
2280  trig_name.push_back("HLT_DoubleMu8_Mass8_PFHT300_v"); // 9
2281  trig_name.push_back("HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT300_v"); // 10
2282 
2283  trig_name.push_back("HLT_PFHT475_v"); // 11
2284  trig_name.push_back("HLT_PFHT800_v"); // 12
2285  trig_name.push_back("HLT_PFMET100_PFMHT100_IDTight_v"); // 13
2286  trig_name.push_back("HLT_PFMET110_PFMHT110_IDTight_v"); // 14
2287  trig_name.push_back("HLT_PFMETNoMu110_PFMHTNoMu110_IDTight_v"); // 15
2288  trig_name.push_back("HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_v"); // 16
2289  trig_name.push_back("HLT_Mu45_eta2p1_v"); // 17
2290  trig_name.push_back("HLT_IsoMu18_v"); // 18
2291  trig_name.push_back("HLT_IsoMu24_v"); // 19
2292  trig_name.push_back("HLT_IsoMu27_v"); // 20
2293 
2294  trig_name.push_back("HLT_Mu50_v"); // 21
2295  trig_name.push_back("HLT_Ele27_eta2p1_WPLoose_Gsf_v"); // 22
2296  trig_name.push_back("HLT_Ele25_eta2p1_WPTight_Gsf_v"); // 23
2297  trig_name.push_back("HLT_Ele105_CaloIdVT_GsfTrkIdT_v"); // 24
2298  trig_name.push_back("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v"); // 25
2299  trig_name.push_back("HLT_Photon175_v"); // 26
2300  trig_name.push_back("HLT_Photon90_CaloIdL_PFHT500_v"); // 27
2301  trig_name.push_back("HLT_PFMET90_PFMHT90_IDTight_v"); // 28
2302  trig_name.push_back("HLT_Ele23_WPLoose_Gsf_v"); // 29
2303  trig_name.push_back("HLT_PFMET120_PFMHT120_IDTight_v"); // 30
2304 
2305  trig_name.push_back("HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_v"); // 31
2306  trig_name.push_back("HLT_IsoMu22_v"); // 32
2307  trig_name.push_back("HLT_PFMETNoMu100_PFMHTNoMu100_IDTight_v"); // 33
2308  trig_name.push_back("HLT_Mu50_IsoVVVL_PFHT400_v"); // 34
2309  trig_name.push_back("HLT_Mu15_IsoVVVL_BTagCSV_p067_PFHT400_v"); // 35
2310  trig_name.push_back("HLT_Ele50_IsoVVVL_PFHT400_v"); // 36
2311  trig_name.push_back("HLT_Ele15_IsoVVVL_BTagCSV_p067_PFHT400_v"); // 37
2312  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT400_PFMET50_v"); // 38
2313  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT400_PFMET50_v"); // 39
2314  trig_name.push_back("HLT_Ele27_WPTight_Gsf_v"); // 40
2315 
2316  trig_name.push_back("HLT_Ele115_CaloIdVT_GsfTrkIdT_v"); // 41
2317  trig_name.push_back("HLT_IsoMu22_eta2p1_v"); // 42
2318  trig_name.push_back("HLT_PFHT300_PFMET110_v"); // 43
2319  trig_name.push_back("HLT_PFHT200_DiPFJetAve90_PFAlphaT0p63_v"); // 44
2320  trig_name.push_back("HLT_PFHT250_DiPFJetAve90_PFAlphaT0p58_v"); // 45
2321  trig_name.push_back("HLT_PFHT300_DiPFJetAve90_PFAlphaT0p54_v"); // 46
2322  trig_name.push_back("HLT_PFHT200_v"); // 47
2323  trig_name.push_back("HLT_PFHT250_v"); // 48
2324  trig_name.push_back("HLT_PFHT300_v"); // 49
2325  trig_name.push_back("HLT_PFHT350_v"); // 50
2326 
2327  trig_name.push_back("HLT_PFHT400_v"); // 51
2328  trig_name.push_back("HLT_PFHT600_v"); // 52
2329  trig_name.push_back("HLT_PFHT650_v"); // 53
2330  trig_name.push_back("HLT_PFHT900_v"); // 54
2331  trig_name.push_back("HLT_IsoTkMu24_v"); // 55
2332  trig_name.push_back("HLT_PFJet450_v"); // 56
2333  trig_name.push_back("HLT_AK8PFJet450_v"); // 57
2334  } else {
2335  trig_name.push_back("HLT_PFHT350_PFMET100_"); // 0
2336  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT350_PFMET50_v"); // 1
2337  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT600_v"); // 2
2338  trig_name.push_back("HLT_Mu15_IsoVVVL_BTagCSV0p72_PFHT400_v"); // 3
2339  trig_name.push_back("HLT_Mu15_IsoVVVL_PFHT350_v"); // 4
2340  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT350_PFMET50_v"); // 5
2341  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT600_v"); // 6
2342  trig_name.push_back("HLT_Ele15_IsoVVVL_BTagCSV0p72_PFHT400_v"); // 7
2343  trig_name.push_back("HLT_Ele15_IsoVVVL_PFHT350_v"); // 8
2344  trig_name.push_back("HLT_DoubleMu8_Mass8_PFHT300_v"); // 9
2345  trig_name.push_back("HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT300_v"); // 10
2346  trig_name.push_back("HLT_PFHT475_v"); // 11
2347  trig_name.push_back("HLT_PFHT800_v"); // 12
2348  trig_name.push_back("HLT_PFMET120_JetIdCleaned_Mu5_v"); // 13
2349  trig_name.push_back("HLT_PFMET170_JetIdCleaned_v"); // 14
2350  trig_name.push_back("HLT_DoubleIsoMu17_eta2p1_v"); // 15
2351  trig_name.push_back("HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_v"); // 16
2352  trig_name.push_back("HLT_IsoMu20_v"); // 17
2353  trig_name.push_back("HLT_IsoMu18_v"); // 18
2354  trig_name.push_back("HLT_IsoMu24_eta2p1_v"); // 19
2355  trig_name.push_back("HLT_IsoMu27_v"); // 20
2356  trig_name.push_back("HLT_Mu50_v"); // 21
2357  trig_name.push_back("HLT_Ele27_eta2p1_WPLoose_Gsf_v"); // 22
2358  trig_name.push_back("HLT_Ele23_WPLoose_Gsf_v"); // 23
2359  trig_name.push_back("HLT_Ele105_CaloIdVT_GsfTrkIdT_v"); // 24
2360  trig_name.push_back("HLT_DoubleEle33_CaloIdL_GsfTrkIdVL_v"); // 25
2361  trig_name.push_back("HLT_DoubleEle24_22_eta2p1_WPLoose_Gsf_v"); // 26
2362  trig_name.push_back("HLT_Photon90_CaloIdL_PFHT500_v"); // 27
2363  trig_name.push_back("HLT_PFMETNoMu90_JetIdCleaned_PFMHTNoMu90_IDTight_v"); // 28
2364  }
2365  } else {
2366  trig_name.push_back("HLT_PFMET90_PFMHT90_IDTight_v7"); // 0
2367  trig_name.push_back("HLT_PFMET100_PFMHT100_IDTight_v7"); // 1
2368  trig_name.push_back("HLT_PFMET110_PFMHT110_IDTight_v7"); // 2
2369  trig_name.push_back("HLT_PFMET120_PFMHT120_IDTight_v7"); // 3
2370  trig_name.push_back("HLT_PFMETNoMu90_PFMHTNoMu90_IDTight_v7"); // 4
2371  trig_name.push_back("HLT_PFMETNoMu100_PFMHTNoMu100_IDTight_v7"); // 5
2372  trig_name.push_back("HLT_PFMETNoMu110_PFMHTNoMu110_IDTight_v7"); // 6
2373  trig_name.push_back("HLT_PFMETNoMu120_PFMHTNoMu120_IDTight_v7"); // 7
2374  }
2375 
2376 }
2377 
2378 
2380  outfile->cd();
2381  baby.tree_.SetDirectory(outfile);
2382  baby.Write();
2383 
2384  string commit_s = execute("git rev-parse HEAD");
2385  while(!commit_s.empty() && commit_s.at(commit_s.length()-1) == '\n') commit_s.erase(commit_s.length()-1);
2386  TString commit = commit_s;
2387  TString type = baby.BabyType();
2388  TString root_version = gROOT->GetVersion();
2389  TString root_tutorial_dir = gROOT->GetTutorialsDir();
2390  TString user(getenv("ORIGIN_USER"));
2391  if (user=="") user = getenv("USER");
2392  TString cmssw(getenv("CMSSW_BASE"));
2393  time_t curTime;
2394  time(&curTime);
2395  char time_c[100];
2396  struct tm * timeinfo = localtime(&curTime);
2397  strftime(time_c,100,"%Y-%m-%d %H:%M:%S",timeinfo);
2398  TString date(time_c);
2399  int seconds(floor(difftime(curTime,startTime)+0.5));
2400 
2401  vector<TString> sys_names;
2402  sys_names.resize(kSysLast,"");
2403  sys_names[kSysJER] = "jer";
2404  sys_names[kSysJECUp] = "jec_up";
2405  sys_names[kSysJECDn] = "jec_dn";
2406 
2407  TTree treeglobal("treeglobal", "treeglobal");
2408  treeglobal.Branch("nev_sample", &nevents_sample);
2409  treeglobal.Branch("nev_file", &nevents);
2410  treeglobal.Branch("runtime_seconds", &seconds);
2411  treeglobal.Branch("git_commit", &commit);
2412  // treeglobal.Branch("model", &model);
2413  treeglobal.Branch("baby_type", &type);
2414  treeglobal.Branch("root_version", &root_version);
2415  treeglobal.Branch("root_tutorial_dir", &root_tutorial_dir);
2416  treeglobal.Branch("trig_names", &trig_name);
2417  treeglobal.Branch("sys_names", &sys_names);
2418  treeglobal.Branch("xsec", &xsec);
2419  treeglobal.Branch("user", &user);
2420  treeglobal.Branch("cmssw", &cmssw);
2421  treeglobal.Branch("jec", &jec_label);
2422  treeglobal.Branch("json", &jsonfile);
2423  treeglobal.Branch("date", &date);
2424  treeglobal.Branch("inputfiles", &inputfiles);
2425  treeglobal.Branch("condor_subtime", &condor_subtime);
2426  treeglobal.Fill();
2427  treeglobal.SetDirectory(outfile);
2428  treeglobal.Write();
2429 
2430  outfile->Close();
2431 
2432  int minutes((seconds/60)%60), hours(seconds/3600);
2433  TString runtime("");
2434  if(hours<10) runtime += "0";
2435  runtime += hours; runtime += ":";
2436  if(minutes<10) runtime += "0";
2437  runtime += minutes; runtime += ":";
2438  if((seconds%60)<10) runtime += "0";
2439  runtime += seconds%60;
2440  float hertz(nevents); hertz /= seconds;
2441  cout<<endl<<"BABYMAKER: Written "<<nevents<<" events in "<<outname<<". It took "<<seconds<<" seconds to run ("<<runtime<<"), "
2442  <<roundNumber(hertz,1)<<" Hz, "<<roundNumber(1000,2,hertz)<<" ms per event"<<endl<<endl;
2443  cout<<"BABYMAKER: *********** List of input files ***********"<<endl;
2444  for(size_t ifile(0); ifile < inputfiles.size(); ifile++)
2445  cout<<"BABYMAKER: "<<inputfiles[ifile].c_str()<<endl;
2446  cout<<endl;
2447 
2448  delete outfile;
2449 
2450  delete lepTool;
2451  delete photonTool;
2452  delete jetTool;
2453  delete mcTool;
2454  delete weightTool;
2455 }
2456 
2457 void bmaker_full::reportTime(const edm::Event& iEvent){
2458  // Time reporting
2459  if(nevents==1) {
2460  time_t curTime;
2461  time(&curTime);
2462  cout<<endl<<"BABYMAKER: Took "<<roundNumber(difftime(curTime,startTime),1)<<" seconds for set up and run first event"
2463  <<endl<<endl;
2464  time(&startTime);
2465  }
2466  if(debug || (nevents<100&&nevents%10==0) || (nevents<1000&&nevents%100==0)
2467  || (nevents<10000&&nevents%1000==0) || nevents%10000==0) {
2468  time_t curTime;
2469  time(&curTime);
2470  float seconds(difftime(curTime,startTime));
2471  cout<<"BABYMAKER: Run "<<iEvent.id().run()<<", Event "<< setw(8)<<iEvent.id().event()
2472  <<", LumiSection "<< setw(5)<< iEvent.luminosityBlock()
2473  <<". Ran "<<setw(7)<<nevents<<" events in "<<setw(7)<<seconds<<" seconds -> "
2474  <<setw(5)<<roundNumber(nevents-1,1,seconds)<<" Hz, "
2475  <<setw(5)<<roundNumber(seconds*1000,2,nevents-1)<<" ms per event"<<endl;
2476  }
2477 }
2478 
2479 // ------------ method called once each job just before starting event loop ------------
2481 }
2482 
2483 // ------------ method called once each job just after ending the event loop ------------
2485 }
2486 
2487 // ------------ method called when starting to processes a run ------------
2488 /*
2489  void
2490  bmaker_full::beginRun(edm::Run const&, edm::EventSetup const&)
2491  {
2492  }
2493 */
2494 
2495 // ------------ method called when ending the processing of a run ------------
2496 /*
2497  void
2498  bmaker_full::endRun(edm::Run const&, edm::EventSetup const&)
2499  {
2500  }
2501 */
2502 
2503 // ------------ method called when starting to processes a luminosity block ------------
2504 
2505 void bmaker_full::beginLuminosityBlock(edm::LuminosityBlock const& iLumi, edm::EventSetup const& iEvent){
2506  if (outname.Contains("PUSpring16Fast") && outname.Contains("SMS-")){
2507  edm::Handle<GenLumiInfoHeader> gen_header;
2508 // iLumi.getByToken(tok_genlumiheader_, gen_header);
2509  string model = gen_header->configDescription();
2510  mcTool->getMassPoints(model, mprod_, mlsp_);
2511  }
2512 }
2513 
2514 
2515 // ------------ method called when ending the processing of a luminosity block ------------
2516 /*
2517  void
2518  bmaker_full::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
2519  {
2520  }
2521 */
2522 
2523 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
2524 void
2525 bmaker_full::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
2526  //The following says we do not know what parameters are allowed so do no validation
2527  // Please change this to state exactly what you do use, even if it is no parameters
2528  edm::ParameterSetDescription desc;
2529  desc.setUnknown();
2530  descriptions.addDefault(desc);
2531 }
2532 
2533 //define this as a plug-in
float const & elel_pt2() const
Definition: baby_base.cc:4208
float fractionNegWeights(const TString &file)
float const & mumuv_pt1() const
Definition: baby_base.cc:5253
TTree tree_
Definition: baby_base.hh:765
float const & mumu_pt1() const
Definition: baby_base.cc:5176
float const & elmu_w() const
Definition: baby_base.cc:4373
std::vector< float > const & mus_pt() const
Definition: baby_base.cc:6936
float const & elelv_w() const
Definition: baby_base.cc:4296
float const & mumuv_m() const
Definition: baby_base.cc:5220
TString outname
Definition: bmaker_full.hh:156
float const & elel_phi() const
Definition: baby_base.cc:4175
std::vector< bool > const & mus_inz() const
Definition: baby_base.cc:6221
time_t startTime
Definition: bmaker_full.hh:70
std::vector< bool > const & els_inz() const
Definition: baby_base.cc:6089
float crossSection(const TString &file)
void writeWeights(const vCands &sig_leps, edm::Handle< GenEventInfoProduct > &gen_event_info, edm::Handle< LHEEventProduct > &lhe_info)
bool writeTriggers(const edm::TriggerNames &names, edm::Handle< edm::TriggerResults > triggerBits, edm::Handle< pat::PackedTriggerPrescales > triggerPrescales)
void Write()
Definition: baby_base.cc:3282
void writeHLTObjects(const edm::TriggerNames &names, edm::Handle< pat::TriggerObjectStandAloneCollection > triggerObjects, vCands &all_mus, vCands &all_els, vCands &photons, const edm::Event &)
static std::vector< size_t > getBRanking(const std::vector< LVector > &momentum, const std::vector< float > &csv, const std::vector< bool > &is_lep)
float const & elmu_pt2() const
Definition: baby_base.cc:4362
void setDiLepMass(vCands leptons, baby_float ll_m, baby_float ll_pt1, baby_float ll_pt2, baby_float ll_pt, baby_float ll_eta, baby_float ll_phi, baby_vfloat l_pt, baby_vbool l_inz, baby_float ll_w)
float const & elel_m() const
Definition: baby_base.cc:4164
float const & elelv_pt() const
Definition: baby_base.cc:4263
bmaker_full(const edm::ParameterSet &)
virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) override
float const & mumu_eta() const
Definition: baby_base.cc:5132
tuple condor_subtime
Setting global tag From https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookJetEnergyCorrections#J...
std::vector< const reco::Candidate * > vCands
Definition: utilities.hh:15
STL namespace.
list reco
Definition: get_flist.py:50
weight_tools * weightTool
Definition: bmaker_full.hh:78
float const & elelv_pt2() const
Definition: baby_base.cc:4285
std::vector< bool > const & mus_inzv() const
Definition: baby_base.cc:6232
lepton_tools * lepTool
Definition: bmaker_full.hh:73
float const & mumuv_w() const
Definition: baby_base.cc:5275
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: bmaker_full.cc:41
isData
If you only want to re-correct and get the proper uncertainties, no reclustering. ...
float const & mumuv_pt() const
Definition: baby_base.cc:5242
float dR(float eta1, float eta2, float phi1, float phi2)
Definition: utilities.cc:474
baby_full baby
Definition: bmaker_full.hh:68
float const & elmu_m() const
Definition: baby_base.cc:4318
void writeDiLep(vCands &sig_mus, vCands &sig_els, vCands &veto_mus, vCands &veto_els)
float const & mumuv_phi() const
Definition: baby_base.cc:5231
void writeMET(edm::Handle< pat::METCollection > mets, edm::Handle< pat::METCollection > mets_nohf, edm::Handle< pat::METCollection > mets_puppi)
Definition: bmaker_full.cc:369
float getMT(float pt1, float phi1, float pt2, float phi2)
Definition: utilities.cc:44
float const & elelv_pt1() const
Definition: baby_base.cc:4274
vCands writeElectrons(edm::Handle< pat::ElectronCollection > electrons, edm::Handle< pat::PackedCandidateCollection > pfcands, edm::Handle< reco::VertexCollection > vtx, vCands &veto_els, vCands &all_els, double rhoEventCentral)
void rebalancedMET(double &MET, double &METPhi)
TString jec_label
Definition: bmaker_full.hh:160
tuple col
Definition: sub_cond.py:196
float const & mumu_m() const
Definition: baby_base.cc:5143
float const & elmu_pt1() const
Definition: baby_base.cc:4351
void reportTime(const edm::Event &iEvent)
float const & elelv_m() const
Definition: baby_base.cc:4241
float const & mumu_phi() const
Definition: baby_base.cc:5154
float &(baby_base::* baby_float)()
Definition: bmaker_full.hh:57
TString condor_subtime
Definition: bmaker_full.hh:159
float const & elel_pt() const
Definition: baby_base.cc:4186
bool greaterPt(const reco::Candidate *a, const reco::Candidate *b)
Definition: utilities.cc:36
void writeVertices(edm::Handle< reco::VertexCollection > vtx, edm::Handle< std::vector< PileupSummaryInfo > > pu_info)
unsigned int nevents
Definition: bmaker_full.hh:167
void writeTks(edm::Handle< pat::PackedCandidateCollection > pfcands, edm::Handle< reco::VertexCollection > vtx, double rhoEventCentral)
std::vector< float > const & els_pt() const
Definition: baby_base.cc:6529
void writeBBVars(std::vector< reco::Candidate::LorentzVector > &all_baby_jets, vCands &sig_leps)
Definition: bmaker_full.cc:863
bool Contains(const std::string &text, const std::string &pattern)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
TString roundNumber(double num, int decimals, double denom=1.)
Definition: utilities.cc:478
mc_tools * mcTool
Definition: bmaker_full.hh:77
float const & elel_pt1() const
Definition: baby_base.cc:4197
bool doSystematics
Definition: bmaker_full.hh:173
void getMassPoints(TString mpoints, int &mgluino, int &mlsp)
Definition: mc_tools.cc:131
float const & mumuv_pt2() const
Definition: baby_base.cc:5264
float const & mumu_pt2() const
Definition: baby_base.cc:5187
virtual std::string BabyType() const
Definition: baby_full.cc:31
void setElMuMass(vCands leptons1, vCands leptons2, baby_float ll_m, baby_float ll_pt1, baby_float ll_pt2, baby_float ll_pt, baby_float ll_eta, baby_float ll_phi, baby_float ll_w)
std::vector< bool > &(baby_base::* baby_vbool)()
Definition: bmaker_full.hh:59
tuple ind
Definition: resubmit.py:140
reco::Candidate::LorentzVector LVector
Definition: utilities.hh:16
static int type(const std::string &name)
Definition: event_tools.cc:104
float const & elmu_eta() const
Definition: baby_base.cc:4307
void writeIFSR(edm::Handle< reco::GenParticleCollection > genParticles, std::vector< reco::Candidate::LorentzVector > &jets)
void writeBTagWeights(edm::Handle< pat::JetCollection > alljets, std::vector< reco::Candidate::LorentzVector > &all_baby_jets, std::vector< unsigned > &all_baby_jet_idx)
Definition: bmaker_full.cc:603
float const & mumuv_eta() const
Definition: baby_base.cc:5209
void writeGenInfo(edm::Handle< LHEEventProduct > lhe_info)
virtual void endJob() override
float const & elelv_phi() const
Definition: baby_base.cc:4252
float dPhi(float phi1, float phi2)
Definition: utilities.cc:24
std::string execute(const std::string &cmd)
float const & elmu_phi() const
Definition: baby_base.cc:4329
float const & mumu_pt() const
Definition: baby_base.cc:5165
unsigned int nevents_sample
Definition: bmaker_full.hh:166
float const & elel_w() const
Definition: baby_base.cc:4219
vCands writePhotons(edm::Handle< pat::PhotonCollection > allphotons, edm::Handle< std::vector< pat::Electron > > &electrons, edm::Handle< reco::ConversionCollection > &conversions, edm::Handle< reco::BeamSpot > &beamspot, double rho)
float getMT2(float pt1, float phi1, float pt2, float phi2, float met, float met_phi)
Definition: utilities.cc:54
event_tools * eventTool
Definition: bmaker_full.hh:76
void writeFatJets()
Definition: bmaker_full.cc:889
float const & elel_eta() const
Definition: baby_base.cc:4153
static std::pair< double, double > getScaleFactor(const reco::Candidate &cand)
TFile * outfile
Definition: bmaker_full.hh:67
TString jsonfile
Definition: bmaker_full.hh:158
DEFINE_FWK_MODULE(bmaker_full)
std::vector< TString > trig_name
Definition: bmaker_full.hh:153
photon_tools * photonTool
Definition: bmaker_full.hh:75
virtual void beginJob() override
void writeHiggVars(std::vector< LVector > &baby_jets_p4, std::vector< float > &baby_jets_csv, std::vector< bool > &baby_jets_h1, std::vector< bool > &baby_jets_h2, std::vector< bool > &baby_jets_islep, int &baby_nbl, int &baby_nbm, int &baby_nbt, float &baby_hig_am, float &baby_hig_dm, float &baby_hig_drmax, int &baby_hig_bin, float &baby_mct, bool isSystemtatic=false)
Definition: bmaker_full.cc:757
std::vector< float > &(baby_base::* baby_vfloat)()
Definition: bmaker_full.hh:58
vCands writeMuons(edm::Handle< pat::MuonCollection > muons, edm::Handle< pat::PackedCandidateCollection > pfcands, edm::Handle< reco::VertexCollection > vtx, vCands &veto_mus, vCands &all_mus, double rhoEventCentral)
Definition: bmaker_full.cc:982
std::vector< reco::Candidate::LorentzVector > writeJets(edm::Handle< pat::JetCollection > alljets, std::vector< unsigned > &all_baby_jets_idx, edm::Handle< edm::View< reco::GenJet > > genjets, vCands &sig_leps, vCands &veto_leps, vCands &photons, vCands &tks, std::vector< std::vector< LVector > > &sysjets, std::vector< double > &jetsMuonEnergyFrac)
Definition: bmaker_full.cc:402
float const & elelv_eta() const
Definition: baby_base.cc:4230
void writeLeptons(vCands &leptons)
void writeMC(edm::Handle< reco::GenParticleCollection > genParticles, vCands &all_mus, vCands &all_els, vCands &photons)
void writeFilters(const edm::TriggerNames &fnames, edm::Handle< edm::TriggerResults > filterBits, edm::Handle< reco::VertexCollection > vtx, std::vector< double > jetsMuonEnergyFrac)
double calculateRebalancedMET(unsigned int jetIdx, double mu, double &METPhi)
jet_met_tools * jetTool
Definition: bmaker_full.hh:74
std::vector< bool > const & els_inzv() const
Definition: baby_base.cc:6100
std::vector< std::string > inputfiles
Definition: bmaker_full.hh:157
float const & elmu_pt() const
Definition: baby_base.cc:4340
double calculateRescalingFactor(unsigned int jetIdx)
float const & mumu_w() const
Definition: baby_base.cc:5198