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