CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TtFullLepKinSolutionProducer.h
Go to the documentation of this file.
1 #ifndef TtFullLepKinSolutionProducer_h
2 #define TtFullLepKinSolutionProducer_h
3 
4 //
5 //
6 #include <memory>
7 #include <string>
8 #include <vector>
9 #include "TLorentzVector.h"
10 #include "DataFormats/Math/interface/deltaR.h"
11 #include "FWCore/MessageLogger/interface/MessageLogger.h"
12 #include "FWCore/Framework/interface/EDProducer.h"
13 #include "FWCore/Framework/interface/Event.h"
14 #include "FWCore/ParameterSet/interface/ParameterSet.h"
15 #include "FWCore/Utilities/interface/InputTag.h"
16 #include "DataFormats/Candidate/interface/LeafCandidate.h"
17 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullLepKinSolver.h"
18 
19 
20 class TtFullLepKinSolutionProducer : public edm::EDProducer {
21 
22  public:
23 
24  explicit TtFullLepKinSolutionProducer(const edm::ParameterSet & iConfig);
25  ~TtFullLepKinSolutionProducer();
26 
27  virtual void beginJob();
28  virtual void produce(edm::Event & evt, const edm::EventSetup & iSetup);
29  virtual void endJob();
30 
31  private:
32 
33  // next methods are avoidable but they make the code legible
34  inline bool PTComp(const reco::Candidate*, const reco::Candidate*) const;
35  inline bool LepDiffCharge(const reco::Candidate* , const reco::Candidate*) const;
36  inline bool HasPositiveCharge(const reco::Candidate*) const;
37 
38  struct Compare{
39  bool operator()(std::pair<double, int> a, std::pair<double, int> b){
40  return a.first > b.first;
41  }
42  };
43 
44  edm::EDGetTokenT<std::vector<pat::Jet> > jetsToken_;
45  edm::EDGetTokenT<std::vector<pat::Electron> > electronsToken_;
46  edm::EDGetTokenT<std::vector<pat::Muon> > muonsToken_;
47  edm::EDGetTokenT<std::vector<pat::MET> > metsToken_;
48 
49  std::string jetCorrLevel_;
50  int maxNJets_, maxNComb_;
51  bool eeChannel_, emuChannel_, mumuChannel_, searchWrongCharge_;
52  double tmassbegin_, tmassend_, tmassstep_;
53  std::vector<double> nupars_;
54 
55  TtFullLepKinSolver* solver;
56 };
57 
58 inline bool TtFullLepKinSolutionProducer::PTComp(const reco::Candidate* l1, const reco::Candidate* l2) const
59 {
60  return (l1->pt() > l2->pt());
61 }
62 
63 inline bool TtFullLepKinSolutionProducer::LepDiffCharge(const reco::Candidate* l1, const reco::Candidate* l2) const
64 {
65  return (l1->charge() != l2->charge());
66 }
67 
68 inline bool TtFullLepKinSolutionProducer::HasPositiveCharge(const reco::Candidate* l) const
69 {
70  return (l->charge() > 0);
71 }
72 
73 TtFullLepKinSolutionProducer::TtFullLepKinSolutionProducer(const edm::ParameterSet & iConfig)
74 {
75  // configurables
76  jetsToken_ = consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("jets"));
77  electronsToken_ = consumes<std::vector<pat::Electron> >(iConfig.getParameter<edm::InputTag>("electrons"));
78  muonsToken_ = consumes<std::vector<pat::Muon> >(iConfig.getParameter<edm::InputTag>("muons"));
79  metsToken_ = consumes<std::vector<pat::MET> >(iConfig.getParameter<edm::InputTag>("mets"));
80  jetCorrLevel_= iConfig.getParameter<std::string> ("jetCorrectionLevel");
81  maxNJets_ = iConfig.getParameter<int> ("maxNJets");
82  maxNComb_ = iConfig.getParameter<int> ("maxNComb");
83  eeChannel_ = iConfig.getParameter<bool>("eeChannel");
84  emuChannel_ = iConfig.getParameter<bool>("emuChannel");
85  mumuChannel_ = iConfig.getParameter<bool>("mumuChannel");
86  searchWrongCharge_ = iConfig.getParameter<bool> ("searchWrongCharge");
87  tmassbegin_ = iConfig.getParameter<double>("tmassbegin");
88  tmassend_ = iConfig.getParameter<double>("tmassend");
89  tmassstep_ = iConfig.getParameter<double>("tmassstep");
90  nupars_ = iConfig.getParameter<std::vector<double> >("neutrino_parameters");
91 
92  // define what will be produced
93  produces<std::vector<std::vector<int> > > (); // vector of the particle inices (b, bbar, e1, e2, mu1, mu2)
94  produces<std::vector<reco::LeafCandidate> >("fullLepNeutrinos");
95  produces<std::vector<reco::LeafCandidate> >("fullLepNeutrinoBars");
96  produces<std::vector<double> >("solWeight"); //weight for a specific kin solution
97  produces<bool>("isWrongCharge"); //true if leptons have the same charge
98 }
99 
100 TtFullLepKinSolutionProducer::~TtFullLepKinSolutionProducer()
101 {
102 }
103 
104 void TtFullLepKinSolutionProducer::beginJob()
105 {
106  solver = new TtFullLepKinSolver(tmassbegin_, tmassend_, tmassstep_, nupars_);
107 }
108 
109 void TtFullLepKinSolutionProducer::endJob()
110 {
111  delete solver;
112 }
113 
114 void TtFullLepKinSolutionProducer::produce(edm::Event & evt, const edm::EventSetup & iSetup)
115 {
116  //create vectors fo runsorted output
117  std::vector<std::vector<int> > idcsV;
118  std::vector<reco::LeafCandidate> nusV;
119  std::vector<reco::LeafCandidate> nuBarsV;
120  std::vector<std::pair<double, int> > weightsV;
121 
122  //create pointer for products
123  std::auto_ptr<std::vector<std::vector<int> > > pIdcs(new std::vector<std::vector<int> >);
124  std::auto_ptr<std::vector<reco::LeafCandidate> > pNus(new std::vector<reco::LeafCandidate>);
125  std::auto_ptr<std::vector<reco::LeafCandidate> > pNuBars(new std::vector<reco::LeafCandidate>);
126  std::auto_ptr<std::vector<double> > pWeight(new std::vector<double>);
127  std::auto_ptr<bool> pWrongCharge(new bool);
128 
129  edm::Handle<std::vector<pat::Jet> > jets;
130  evt.getByToken(jetsToken_, jets);
131  edm::Handle<std::vector<pat::Electron> > electrons;
132  evt.getByToken(electronsToken_, electrons);
133  edm::Handle<std::vector<pat::Muon> > muons;
134  evt.getByToken(muonsToken_, muons);
135  edm::Handle<std::vector<pat::MET> > mets;
136  evt.getByToken(metsToken_, mets);
137 
138  int selMuon1 = -1, selMuon2 = -1;
139  int selElectron1 = -1, selElectron2 = -1;
140  bool ee = false;
141  bool emu = false;
142  bool mumu = false;
143  bool isWrongCharge = false;
144  bool jetsFound = false;
145  bool METFound = false;
146  bool electronsFound = false;
147  bool electronMuonFound = false;
148  bool muonsFound = false;
149 
150  //select Jets (TopJet vector is sorted on ET)
151  if(jets->size()>=2) { jetsFound = true; }
152 
153  //select MET (TopMET vector is sorted on ET)
154  if(mets->size()>=1) { METFound = true; }
155 
156  // If we have electrons and muons available,
157  // build a solutions with electrons and muons.
158  if(muons->size() + electrons->size() >=2) {
159  // select leptons
160  if(electrons->size() == 0) mumu = true;
161  else if(muons->size() == 0) ee = true;
162  else if(electrons->size() == 1) {
163  if(muons->size() == 1) emu = true;
164  else if(PTComp(&(*electrons)[0], &(*muons)[1])) emu = true;
165  else mumu = true;
166  }
167  else if(electrons->size() > 1) {
168  if(PTComp(&(*electrons)[1], &(*muons)[0])) ee = true;
169  else if(muons->size() == 1) emu = true;
170  else if(PTComp(&(*electrons)[0], &(*muons)[1])) emu = true;
171  else mumu = true;
172  }
173  if(ee && eeChannel_) {
174  if(LepDiffCharge(&(*electrons)[0], &(*electrons)[1]) || searchWrongCharge_) {
175  if(HasPositiveCharge(&(*electrons)[0]) || !LepDiffCharge(&(*electrons)[0], &(*electrons)[1])) {
176  selElectron1 = 0;
177  selElectron2 = 1;
178  }
179  else{
180  selElectron1 = 1;
181  selElectron2 = 0;
182  }
183  electronsFound = true;
184  if(!LepDiffCharge(&(*electrons)[0], &(*electrons)[1])) isWrongCharge = true;
185  }
186  }
187  else if(emu && emuChannel_) {
188  if(LepDiffCharge(&(*electrons)[0], &(*muons)[0]) || searchWrongCharge_) {
189  selElectron1 = 0;
190  selMuon1 = 0;
191  electronMuonFound = true;
192  if(!LepDiffCharge(&(*electrons)[0], &(*muons)[0])) isWrongCharge = true;
193  }
194  }
195  else if(mumu && mumuChannel_) {
196  if(LepDiffCharge(&(*muons)[0], &(*muons)[1]) || searchWrongCharge_) {
197  if(HasPositiveCharge(&(*muons)[0]) || !LepDiffCharge(&(*muons)[0], &(*muons)[1])) {
198  selMuon1 = 0;
199  selMuon2 = 1;
200  }
201  else {
202  selMuon1 = 1;
203  selMuon2 = 0;
204  }
205  muonsFound = true;
206  if(!LepDiffCharge(&(*muons)[0], &(*muons)[1])) isWrongCharge = true;
207  }
208  }
209  }
210 
211  *pWrongCharge = isWrongCharge;
212 
213  // Check that the above work makes sense
214  if(int(ee)+int(emu)+int(mumu)>1) {
215  edm::LogWarning("TtFullLepKinSolutionProducer") << "Lepton selection criteria uncorrectly defined";
216  }
217 
218  // Check if the leptons for the required Channel are available
219  bool correctLeptons = ((electronsFound && eeChannel_) || (muonsFound && mumuChannel_) || (electronMuonFound && emuChannel_) );
220  // Check for equally charged leptons if for wrong charge combinations is searched
221  if(isWrongCharge) correctLeptons *= searchWrongCharge_;
222 
223  if(correctLeptons && METFound && jetsFound) {
224 
225  // run over all jets if input parameter maxNJets is -1 or
226  // adapt maxNJets if number of present jets is smaller than selected
227  // number of jets
228  int stop=maxNJets_;
229  if(jets->size()<static_cast<unsigned int>(stop) || stop<0) stop=jets->size();
230 
231  // counter for number of found kinematic solutions
232  int nSol=0;
233 
234  // consider all permutations
235  for (int ib = 0; ib<stop; ib++) {
236  // second loop of the permutations
237  for (int ibbar = 0; ibbar<stop; ibbar++) {
238  // avoid the diagonal: b and bbar must be distinct jets
239  if(ib==ibbar) continue;
240 
241  std::vector<int> idcs;
242 
243  // push back the indices of the jets
244  idcs.push_back(ib);
245  idcs.push_back(ibbar);
246 
247  TLorentzVector LV_l1;
248  TLorentzVector LV_l2;
249  TLorentzVector LV_b = TLorentzVector((*jets)[ib].correctedJet(jetCorrLevel_, "bottom").px(),
250  (*jets)[ib].correctedJet(jetCorrLevel_, "bottom").py(),
251  (*jets)[ib].correctedJet(jetCorrLevel_, "bottom").pz(),
252  (*jets)[ib].correctedJet(jetCorrLevel_, "bottom").energy() );
253  TLorentzVector LV_bbar = TLorentzVector((*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").px(),
254  (*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").py(),
255  (*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").pz(),
256  (*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").energy());
257 
258  double xconstraint = 0, yconstraint = 0;
259 
260  if (ee) {
261  idcs.push_back(selElectron1);
262  LV_l1.SetXYZT((*electrons)[selElectron1].px(),
263  (*electrons)[selElectron1].py(),
264  (*electrons)[selElectron1].pz(),
265  (*electrons)[selElectron1].energy());
266  xconstraint += (*electrons)[selElectron1].px();
267  yconstraint += (*electrons)[selElectron1].py();
268 
269  idcs.push_back(selElectron2);
270  LV_l2.SetXYZT((*electrons)[selElectron2].px(),
271  (*electrons)[selElectron2].py(),
272  (*electrons)[selElectron2].pz(),
273  (*electrons)[selElectron2].energy());
274  xconstraint += (*electrons)[selElectron2].px();
275  yconstraint += (*electrons)[selElectron2].py();
276 
277  idcs.push_back(-1);
278  idcs.push_back(-1);
279  }
280 
281  else if (emu) {
282  if(!isWrongCharge){
283  if(HasPositiveCharge(&(*electrons)[selElectron1])){
284  idcs.push_back(selElectron1);
285  LV_l1.SetXYZT((*electrons)[selElectron1].px(),
286  (*electrons)[selElectron1].py(),
287  (*electrons)[selElectron1].pz(),
288  (*electrons)[selElectron1].energy());
289  xconstraint += (*electrons)[selElectron1].px();
290  yconstraint += (*electrons)[selElectron1].py();
291 
292  idcs.push_back(-1);
293  idcs.push_back(-1);
294 
295  idcs.push_back(selMuon1);
296  LV_l2.SetXYZT((*muons)[selMuon1].px(),
297  (*muons)[selMuon1].py(),
298  (*muons)[selMuon1].pz(),
299  (*muons)[selMuon1].energy());
300  xconstraint += (*muons)[selMuon1].px();
301  yconstraint += (*muons)[selMuon1].py();
302  }
303  else{
304  idcs.push_back(-1);
305 
306  idcs.push_back(selMuon1);
307  LV_l1.SetXYZT((*muons)[selMuon1].px(),
308  (*muons)[selMuon1].py(),
309  (*muons)[selMuon1].pz(),
310  (*muons)[selMuon1].energy());
311  xconstraint += (*muons)[selMuon1].px();
312  yconstraint += (*muons)[selMuon1].py();
313 
314  idcs.push_back(selElectron1);
315  LV_l2.SetXYZT((*electrons)[selElectron1].px(),
316  (*electrons)[selElectron1].py(),
317  (*electrons)[selElectron1].pz(),
318  (*electrons)[selElectron1].energy());
319  xconstraint += (*electrons)[selElectron1].px();
320  yconstraint += (*electrons)[selElectron1].py();
321 
322  idcs.push_back(-1);
323  }
324  }
325  else{ // means "if wrong charge"
326  if(HasPositiveCharge(&(*electrons)[selElectron1])){ // both leps positive
327  idcs.push_back(selElectron1);
328  LV_l1.SetXYZT((*electrons)[selElectron1].px(),
329  (*electrons)[selElectron1].py(),
330  (*electrons)[selElectron1].pz(),
331  (*electrons)[selElectron1].energy());
332  xconstraint += (*electrons)[selElectron1].px();
333  yconstraint += (*electrons)[selElectron1].py();
334 
335  idcs.push_back(-1);
336 
337  idcs.push_back(selMuon1);
338  LV_l2.SetXYZT((*muons)[selMuon1].px(),
339  (*muons)[selMuon1].py(),
340  (*muons)[selMuon1].pz(),
341  (*muons)[selMuon1].energy());
342  xconstraint += (*muons)[selMuon1].px();
343  yconstraint += (*muons)[selMuon1].py();
344 
345  idcs.push_back(-1);
346  }
347  else{ // both leps negative
348  idcs.push_back(-1);
349 
350  idcs.push_back(selElectron1);
351  LV_l2.SetXYZT((*electrons)[selElectron1].px(),
352  (*electrons)[selElectron1].py(),
353  (*electrons)[selElectron1].pz(),
354  (*electrons)[selElectron1].energy());
355  xconstraint += (*electrons)[selElectron1].px();
356  yconstraint += (*electrons)[selElectron1].py();
357 
358  idcs.push_back(-1);
359 
360  idcs.push_back(selMuon1);
361  LV_l1.SetXYZT((*muons)[selMuon1].px(),
362  (*muons)[selMuon1].py(),
363  (*muons)[selMuon1].pz(),
364  (*muons)[selMuon1].energy());
365  xconstraint += (*muons)[selMuon1].px();
366  yconstraint += (*muons)[selMuon1].py();
367  }
368  }
369  }
370 
371  else if (mumu) {
372  idcs.push_back(-1);
373  idcs.push_back(-1);
374 
375  idcs.push_back(selMuon1);
376  LV_l1.SetXYZT((*muons)[selMuon1].px(),
377  (*muons)[selMuon1].py(),
378  (*muons)[selMuon1].pz(),
379  (*muons)[selMuon1].energy());
380  xconstraint += (*muons)[selMuon1].px();
381  yconstraint += (*muons)[selMuon1].py();
382 
383  idcs.push_back(selMuon2);
384  LV_l2.SetXYZT((*muons)[selMuon2].px(),
385  (*muons)[selMuon2].py(),
386  (*muons)[selMuon2].pz(),
387  (*muons)[selMuon2].energy());
388  xconstraint += (*muons)[selMuon2].px();
389  yconstraint += (*muons)[selMuon2].py();
390  }
391 
392  xconstraint += (*jets)[ib].px() + (*jets)[ibbar].px() + (*mets)[0].px();
393  yconstraint += (*jets)[ib].py() + (*jets)[ibbar].py() + (*mets)[0].py();
394 
395  // calculate neutrino momenta and eventweight
396  solver->SetConstraints(xconstraint, yconstraint);
397  TtFullLepKinSolver::NeutrinoSolution nuSol= solver->getNuSolution( LV_l1, LV_l2 , LV_b, LV_bbar);
398 
399  // add solution to the vectors of solutions if solution exists
400  if(nuSol.weight>0){
401  // add the leptons and jets indices to the vector of combinations
402  idcsV.push_back(idcs);
403 
404  // add the neutrinos
405  nusV.push_back(nuSol.neutrino);
406  nuBarsV.push_back(nuSol.neutrinoBar);
407 
408  // add the solution weight
409  weightsV.push_back(std::make_pair(nuSol.weight, nSol));
410 
411  nSol++;
412  }
413  }
414  }
415  }
416 
417  if(weightsV.size()==0){
418  //create dmummy vector
419  std::vector<int> idcs;
420  for(int i=0; i<6; ++i)
421  idcs.push_back(-1);
422 
423  idcsV.push_back(idcs);
424  weightsV.push_back(std::make_pair(-1,0));
425  reco::LeafCandidate nu;
426  nusV.push_back(nu);
427  reco::LeafCandidate nuBar;
428  nuBarsV.push_back(nuBar);
429  }
430 
431  // check if all vectors have correct length
432  int weightL = weightsV.size();
433  int nuL = nusV.size();
434  int nuBarL = nuBarsV.size();
435  int idxL = idcsV.size();
436 
437  if(weightL!=nuL || weightL!=nuBarL || weightL!=idxL){
438  edm::LogWarning("TtFullLepKinSolutionProducer") << "Output vectors are of different length:"
439  << "\n weight: " << weightL
440  << "\n nu: " << nuL
441  << "\n nubar: " << nuBarL
442  << "\n idcs: " << idxL;
443  }
444 
445  // sort vectors by weight in decreasing order
446  if(weightsV.size()>1){
447  sort(weightsV.begin(), weightsV.end(), Compare());
448  }
449 
450  // determine the number of solutions which is written in the event
451  int stop=weightL;
452  if(maxNComb_>0 && maxNComb_<stop) stop=maxNComb_;
453 
454  for(int i=0; i<stop; ++i){
455  pWeight->push_back(weightsV[i].first);
456  pNus ->push_back(nusV[weightsV[i].second]);
457  pNuBars->push_back(nuBarsV[weightsV[i].second]);
458  pIdcs ->push_back(idcsV[weightsV[i].second]);
459  }
460 
461  // put the results in the event
462  evt.put(pIdcs);
463  evt.put(pNus, "fullLepNeutrinos");
464  evt.put(pNuBars, "fullLepNeutrinoBars");
465  evt.put(pWeight, "solWeight");
466  evt.put(pWrongCharge, "isWrongCharge");
467 }
468 
469 #endif
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
TtFullLepKinSolutionProducer::Compare::operator()
bool operator()(std::pair< double, int > a, std::pair< double, int > b)
Definition: TtFullLepKinSolutionProducer.h:39
TtFullLepKinSolutionProducer::muonsToken_
edm::EDGetTokenT< std::vector< pat::Muon > > muonsToken_
Definition: TtFullLepKinSolutionProducer.h:46
i
int i
Definition: DBlmapReader.cc:9
TtFullLepKinSolutionProducer::LepDiffCharge
bool LepDiffCharge(const reco::Candidate *, const reco::Candidate *) const
Definition: TtFullLepKinSolutionProducer.h:63
TtFullLepKinSolver::getNuSolution
NeutrinoSolution getNuSolution(const TLorentzVector &LV_l, const TLorentzVector &LV_l_, const TLorentzVector &LV_b, const TLorentzVector &LV_b_)
Definition: TtFullLepKinSolver.cc:138
cuy.ib
int ib
Definition: cuy.py:660
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
reco::Candidate::pt
virtual double pt() const =0
transverse momentum
TtFullLepKinSolutionProducer::HasPositiveCharge
bool HasPositiveCharge(const reco::Candidate *) const
Definition: TtFullLepKinSolutionProducer.h:68
TtFullLepKinSolver::SetConstraints
void SetConstraints(const double xx=0, const double yy=0)
Definition: TtFullLepKinSolver.cc:131
HI_PhotonSkim_cff.electrons
tuple electrons
Definition: HI_PhotonSkim_cff.py:77
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:250
TtFullLepKinSolutionProducer::PTComp
bool PTComp(const reco::Candidate *, const reco::Candidate *) const
Definition: TtFullLepKinSolutionProducer.h:58
edm::Event::put
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:121
fwrapper::jets
vector< PseudoJet > jets
Definition: fastjetfortran_madfks.cc:42
reco::Candidate::charge
virtual int charge() const =0
electric charge
TtFullLepKinSolutionProducer::metsToken_
edm::EDGetTokenT< std::vector< pat::MET > > metsToken_
Definition: TtFullLepKinSolutionProducer.h:47
TtFullLepKinSolutionProducer::jetsToken_
edm::EDGetTokenT< std::vector< pat::Jet > > jetsToken_
Definition: TtFullLepKinSolutionProducer.h:44
b
double b
Definition: hdecay.h:120
TtFullLepKinSolutionProducer::TtFullLepKinSolutionProducer
TtFullLepKinSolutionProducer(const edm::ParameterSet &iConfig)
Definition: TtFullLepKinSolutionProducer.h:73
TtFullLepKinSolutionProducer::electronsToken_
edm::EDGetTokenT< std::vector< pat::Electron > > electronsToken_
Definition: TtFullLepKinSolutionProducer.h:45
TtFullLepKinSolutionProducer::produce
virtual void produce(edm::Event &evt, const edm::EventSetup &iSetup)
Definition: TtFullLepKinSolutionProducer.h:114
patZpeak.muons
tuple muons
Definition: patZpeak.py:38
a
double a
Definition: hdecay.h:121
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
Generated for CMSSW Reference Manual by doxygen 1.8.5