CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TtSemiLepHitFitProducer.h
Go to the documentation of this file.
1 #ifndef TtSemiLepHitFitProducer_h
2 #define TtSemiLepHitFitProducer_h
3 
7 
12 
18 
19 template <typename LeptonCollection>
21 
22  public:
23 
26 
27  private:
28  // produce
29  virtual void produce(edm::Event&, const edm::EventSetup&);
30 
34 
36  int maxNJets_;
38  int maxNComb_;
39 
41  double maxEtaMu_;
43  double maxEtaEle_;
45  double maxEtaJet_;
46 
54  bool useBTag_;
55 
57  double mW_;
58  double mTop_;
59 
62 
64  double jes_;
65  double jesB_;
66 
67  struct FitResult {
68  int Status;
69  double Chi2;
70  double Prob;
71  double MT;
72  double SigMT;
79  std::vector<int> JetCombi;
80  bool operator< (const FitResult& rhs) { return Chi2 < rhs.Chi2; };
81  };
82 
84 
91 
96 
98 };
99 
100 template<typename LeptonCollection>
102  jetsToken_ (consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"))),
103  lepsToken_ (consumes<LeptonCollection>(cfg.getParameter<edm::InputTag>("leps"))),
104  metsToken_ (consumes<std::vector<pat::MET> >(cfg.getParameter<edm::InputTag>("mets"))),
105  maxNJets_ (cfg.getParameter<int> ("maxNJets" )),
106  maxNComb_ (cfg.getParameter<int> ("maxNComb" )),
107  bTagAlgo_ (cfg.getParameter<std::string> ("bTagAlgo" )),
108  minBTagValueBJet_ (cfg.getParameter<double> ("minBDiscBJets" )),
109  maxBTagValueNonBJet_ (cfg.getParameter<double> ("maxBDiscLightJets" )),
110  useBTag_ (cfg.getParameter<bool> ("useBTagging" )),
111  mW_ (cfg.getParameter<double> ("mW" )),
112  mTop_ (cfg.getParameter<double> ("mTop" )),
113  jetCorrectionLevel_ (cfg.getParameter<std::string> ("jetCorrectionLevel" )),
114  jes_ (cfg.getParameter<double> ("jes" )),
115  jesB_ (cfg.getParameter<double> ("jesB" )),
116 
117  // The following five initializers read the config parameters for the
118  // ASCII text files which contains the physics object resolutions.
119  hitfitDefault_ (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitDefault"),
120  edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/setting/RunHitFitConfiguration.txt")))),
121  hitfitElectronResolution_(cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitElectronResolution"),
122  edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafElectronResolution.txt")))),
123  hitfitMuonResolution_ (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitMuonResolution"),
124  edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafMuonResolution.txt")))),
125  hitfitUdscJetResolution_ (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitUdscJetResolution"),
126  edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafUdscJetResolution.txt")))),
127  hitfitBJetResolution_ (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitBJetResolution"),
128  edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafBJetResolution.txt")))),
129  hitfitMETResolution_ (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitMETResolution"),
130  edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafKtResolution.txt")))),
131 
132  // The following four initializers instantiate the translator between PAT objects
133  // and HitFit objects using the ASCII text files which contains the resolutions.
134  electronTranslator_(hitfitElectronResolution_.fullPath()),
135  muonTranslator_ (hitfitMuonResolution_.fullPath()),
136  jetTranslator_ (hitfitUdscJetResolution_.fullPath(), hitfitBJetResolution_.fullPath(), jetCorrectionLevel_, jes_, jesB_),
137  metTranslator_ (hitfitMETResolution_.fullPath())
138 
139 {
140  // Create an instance of RunHitFit and initialize it.
146  mW_,
147  mW_,
148  mTop_);
149 
150  maxEtaMu_ = 2.4;
151  maxEtaEle_ = 2.5;
152  maxEtaJet_ = 2.5;
153 
154  edm::LogVerbatim( "TopHitFit" )
155  << "\n"
156  << "+++++++++++ TtSemiLepHitFitProducer ++++++++++++ \n"
157  << " Due to the eta ranges for which resolutions \n"
158  << " are provided in \n"
159  << " TopQuarkAnalysis/TopHitFit/data/resolution/ \n"
160  << " so far, the following cuts are currently \n"
161  << " implemented in the TtSemiLepHitFitProducer: \n"
162  << " |eta(muons )| <= " << maxEtaMu_ << " \n"
163  << " |eta(electrons)| <= " << maxEtaEle_ << " \n"
164  << " |eta(jets )| <= " << maxEtaJet_ << " \n"
165  << "++++++++++++++++++++++++++++++++++++++++++++++++ \n";
166 
167  produces< std::vector<pat::Particle> >("PartonsHadP");
168  produces< std::vector<pat::Particle> >("PartonsHadQ");
169  produces< std::vector<pat::Particle> >("PartonsHadB");
170  produces< std::vector<pat::Particle> >("PartonsLepB");
171  produces< std::vector<pat::Particle> >("Leptons");
172  produces< std::vector<pat::Particle> >("Neutrinos");
173 
174  produces< std::vector<std::vector<int> > >();
175  produces< std::vector<double> >("Chi2");
176  produces< std::vector<double> >("Prob");
177  produces< std::vector<double> >("MT");
178  produces< std::vector<double> >("SigMT");
179  produces< std::vector<int> >("Status");
180  produces< int >("NumberOfConsideredJets");
181 }
182 
183 template<typename LeptonCollection>
185 {
186  delete HitFit;
187 }
188 
189 template<typename LeptonCollection>
191 {
192  std::auto_ptr< std::vector<pat::Particle> > pPartonsHadP( new std::vector<pat::Particle> );
193  std::auto_ptr< std::vector<pat::Particle> > pPartonsHadQ( new std::vector<pat::Particle> );
194  std::auto_ptr< std::vector<pat::Particle> > pPartonsHadB( new std::vector<pat::Particle> );
195  std::auto_ptr< std::vector<pat::Particle> > pPartonsLepB( new std::vector<pat::Particle> );
196  std::auto_ptr< std::vector<pat::Particle> > pLeptons ( new std::vector<pat::Particle> );
197  std::auto_ptr< std::vector<pat::Particle> > pNeutrinos ( new std::vector<pat::Particle> );
198 
199  std::auto_ptr< std::vector<std::vector<int> > > pCombi ( new std::vector<std::vector<int> > );
200  std::auto_ptr< std::vector<double> > pChi2 ( new std::vector<double> );
201  std::auto_ptr< std::vector<double> > pProb ( new std::vector<double> );
202  std::auto_ptr< std::vector<double> > pMT ( new std::vector<double> );
203  std::auto_ptr< std::vector<double> > pSigMT ( new std::vector<double> );
204  std::auto_ptr< std::vector<int> > pStatus( new std::vector<int> );
205  std::auto_ptr< int > pJetsConsidered( new int );
206 
208  evt.getByToken(jetsToken_, jets);
209 
211  evt.getByToken(metsToken_, mets);
212 
214  evt.getByToken(lepsToken_, leps);
215 
216  // -----------------------------------------------------
217  // skip events with no appropriate lepton candidate in
218  // or empty MET or less jets than partons
219  // -----------------------------------------------------
220 
221  const unsigned int nPartons = 4;
222 
223  // Clear the internal state
224  HitFit->clear();
225 
226  // Add lepton into HitFit
227  bool foundLepton = false;
228  if(!leps->empty()) {
229  double maxEtaLep = maxEtaMu_;
230  if( !dynamic_cast<const reco::Muon*>(&((*leps)[0])) ) // assume electron if it is not a muon
231  maxEtaLep = maxEtaEle_;
232  for(unsigned iLep=0; iLep<(*leps).size() && !foundLepton; ++iLep) {
233  if(std::abs((*leps)[iLep].eta()) <= maxEtaLep) {
234  HitFit->AddLepton((*leps)[iLep]);
235  foundLepton = true;
236  }
237  }
238  }
239 
240  // Add jets into HitFit
241  unsigned int nJetsFound = 0;
242  for(unsigned iJet=0; iJet<(*jets).size() && (int)nJetsFound!=maxNJets_; ++iJet) {
243  double jet_eta = (*jets)[iJet].eta();
244  if ( (*jets)[iJet].isCaloJet() ) {
245  jet_eta = ((reco::CaloJet*) ((*jets)[iJet]).originalObject())->detectorP4().eta();
246  }
247  if(std::abs(jet_eta) <= maxEtaJet_) {
248  HitFit->AddJet((*jets)[iJet]);
249  nJetsFound++;
250  }
251  }
252  *pJetsConsidered = nJetsFound;
253 
254  // Add missing transverse energy into HitFit
255  if(!mets->empty())
256  HitFit->SetMet((*mets)[0]);
257 
258  if( !foundLepton || mets->empty() || nJetsFound<nPartons ) {
259  // the kinFit getters return empty objects here
260  pPartonsHadP->push_back( pat::Particle() );
261  pPartonsHadQ->push_back( pat::Particle() );
262  pPartonsHadB->push_back( pat::Particle() );
263  pPartonsLepB->push_back( pat::Particle() );
264  pLeptons ->push_back( pat::Particle() );
265  pNeutrinos ->push_back( pat::Particle() );
266  // indices referring to the jet combination
267  std::vector<int> invalidCombi;
268  for(unsigned int i = 0; i < nPartons; ++i)
269  invalidCombi.push_back( -1 );
270  pCombi->push_back( invalidCombi );
271  // chi2
272  pChi2->push_back( -1. );
273  // chi2 probability
274  pProb->push_back( -1. );
275  // fitted top mass
276  pMT->push_back( -1. );
277  pSigMT->push_back( -1. );
278  // status of the fitter
279  pStatus->push_back( -1 );
280  // feed out all products
281  evt.put(pCombi);
282  evt.put(pPartonsHadP, "PartonsHadP");
283  evt.put(pPartonsHadQ, "PartonsHadQ");
284  evt.put(pPartonsHadB, "PartonsHadB");
285  evt.put(pPartonsLepB, "PartonsLepB");
286  evt.put(pLeptons , "Leptons" );
287  evt.put(pNeutrinos , "Neutrinos" );
288  evt.put(pChi2 , "Chi2" );
289  evt.put(pProb , "Prob" );
290  evt.put(pMT , "MT" );
291  evt.put(pSigMT , "SigMT" );
292  evt.put(pStatus , "Status" );
293  evt.put(pJetsConsidered, "NumberOfConsideredJets");
294  return;
295  }
296 
297  std::list<FitResult> FitResultList;
298 
299  //
300  // BEGIN DECLARATION OF VARIABLES FROM KINEMATIC FIT
301  //
302 
303  // In this part are variables from the
304  // kinematic fit procedure
305 
306  // Number of all permutations of the event
307  size_t nHitFit = 0 ;
308 
309  // Number of jets in the event
310  size_t nHitFitJet = 0 ;
311 
312  // Results of the fit for all jet permutations of the event
313  std::vector<hitfit::Fit_Result> hitFitResult;
314 
315  //
316  // R U N H I T F I T
317  //
318  // Run the kinematic fit and get how many permutations are possible
319  // in the fit
320 
321  nHitFit = HitFit->FitAllPermutation();
322 
323  //
324  // BEGIN PART WHICH EXTRACTS INFORMATION FROM HITFIT
325  //
326 
327  // Get the number of jets
328  nHitFitJet = HitFit->GetUnfittedEvent()[0].njets();
329 
330  // Get the fit results for all permutations
331  hitFitResult = HitFit->GetFitAllPermutation();
332 
333  // Loop over all permutations and extract the information
334  for (size_t fit = 0 ; fit != nHitFit ; ++fit) {
335 
336  // Get the event after the fit
337  hitfit::Lepjets_Event fittedEvent = hitFitResult[fit].ev();
338 
339  /*
340  Get jet permutation according to TQAF convention
341  11 : leptonic b
342  12 : hadronic b
343  13 : hadronic W
344  14 : hadronic W
345  */
346  std::vector<int> hitCombi(4);
347  for (size_t jet = 0 ; jet != nHitFitJet ; ++jet) {
348  int jet_type = fittedEvent.jet(jet).type();
349 
350  switch(jet_type) {
351  case 11: hitCombi[TtSemiLepEvtPartons::LepB ] = jet;
352  break;
353  case 12: hitCombi[TtSemiLepEvtPartons::HadB ] = jet;
354  break;
355  case 13: hitCombi[TtSemiLepEvtPartons::LightQ ] = jet;
356  break;
357  case 14: hitCombi[TtSemiLepEvtPartons::LightQBar] = jet;
358  break;
359  }
360  }
361 
362  // Store the kinematic quantities in the corresponding containers.
363 
364  hitfit::Lepjets_Event_Jet hadP_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQ ]);
365  hitfit::Lepjets_Event_Jet hadQ_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQBar]);
366  hitfit::Lepjets_Event_Jet hadB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::HadB ]);
367  hitfit::Lepjets_Event_Jet lepB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LepB ]);
368  hitfit::Lepjets_Event_Lep lepL_ = fittedEvent.lep(0);
369 
370  /*
372  std::string bTagAlgo_;
374  double minBTagValueBJet_;
376  double maxBTagValueNonBJet_;
378  bool useBTag_;
379  */
380 
381  if ( hitFitResult[fit].chisq() > 0 // only take into account converged fits
382  && (!useBTag_ || ( useBTag_ // use btag information if chosen
383  && jets->at(hitCombi[TtSemiLepEvtPartons::LightQ ]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
384  && jets->at(hitCombi[TtSemiLepEvtPartons::LightQBar]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
385  && jets->at(hitCombi[TtSemiLepEvtPartons::HadB ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
386  && jets->at(hitCombi[TtSemiLepEvtPartons::LepB ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
387  )
388  )
389  ) {
391  result.Status = 0;
392  result.Chi2 = hitFitResult[fit].chisq();
393  result.Prob = exp(-1.0*(hitFitResult[fit].chisq())/2.0);
394  result.MT = hitFitResult[fit].mt();
395  result.SigMT= hitFitResult[fit].sigmt();
397  math::XYZTLorentzVector(hadB_.p().x(), hadB_.p().y(),
398  hadB_.p().z(), hadB_.p().t()), math::XYZPoint()));
400  math::XYZTLorentzVector(hadP_.p().x(), hadP_.p().y(),
401  hadP_.p().z(), hadP_.p().t()), math::XYZPoint()));
403  math::XYZTLorentzVector(hadQ_.p().x(), hadQ_.p().y(),
404  hadQ_.p().z(), hadQ_.p().t()), math::XYZPoint()));
406  math::XYZTLorentzVector(lepB_.p().x(), lepB_.p().y(),
407  lepB_.p().z(), lepB_.p().t()), math::XYZPoint()));
409  math::XYZTLorentzVector(lepL_.p().x(), lepL_.p().y(),
410  lepL_.p().z(), lepL_.p().t()), math::XYZPoint()));
412  math::XYZTLorentzVector(fittedEvent.met().x(), fittedEvent.met().y(),
413  fittedEvent.met().z(), fittedEvent.met().t()), math::XYZPoint()));
414  result.JetCombi = hitCombi;
415 
416  FitResultList.push_back(result);
417  }
418 
419  }
420 
421  // sort results w.r.t. chi2 values
422  FitResultList.sort();
423 
424  // -----------------------------------------------------
425  // feed out result
426  // starting with the JetComb having the smallest chi2
427  // -----------------------------------------------------
428 
429  if( ((unsigned)FitResultList.size())<1 ) { // in case no fit results were stored in the list (all fits aborted)
430  pPartonsHadP->push_back( pat::Particle() );
431  pPartonsHadQ->push_back( pat::Particle() );
432  pPartonsHadB->push_back( pat::Particle() );
433  pPartonsLepB->push_back( pat::Particle() );
434  pLeptons ->push_back( pat::Particle() );
435  pNeutrinos ->push_back( pat::Particle() );
436  // indices referring to the jet combination
437  std::vector<int> invalidCombi;
438  for(unsigned int i = 0; i < nPartons; ++i)
439  invalidCombi.push_back( -1 );
440  pCombi->push_back( invalidCombi );
441  // chi2
442  pChi2->push_back( -1. );
443  // chi2 probability
444  pProb->push_back( -1. );
445  // fitted top mass
446  pMT->push_back( -1. );
447  pSigMT->push_back( -1. );
448  // status of the fitter
449  pStatus->push_back( -1 );
450  }
451  else {
452  unsigned int iComb = 0;
453  for(typename std::list<FitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end(); ++result) {
454  if(maxNComb_ >= 1 && iComb == (unsigned int) maxNComb_) break;
455  iComb++;
456  // partons
457  pPartonsHadP->push_back( result->HadP );
458  pPartonsHadQ->push_back( result->HadQ );
459  pPartonsHadB->push_back( result->HadB );
460  pPartonsLepB->push_back( result->LepB );
461  // lepton
462  pLeptons->push_back( result->LepL );
463  // neutrino
464  pNeutrinos->push_back( result->LepN );
465  // indices referring to the jet combination
466  pCombi->push_back( result->JetCombi );
467  // chi2
468  pChi2->push_back( result->Chi2 );
469  // chi2 probability
470  pProb->push_back( result->Prob );
471  // fitted top mass
472  pMT->push_back( result->MT );
473  pSigMT->push_back( result->SigMT );
474  // status of the fitter
475  pStatus->push_back( result->Status );
476  }
477  }
478  evt.put(pCombi);
479  evt.put(pPartonsHadP, "PartonsHadP");
480  evt.put(pPartonsHadQ, "PartonsHadQ");
481  evt.put(pPartonsHadB, "PartonsHadB");
482  evt.put(pPartonsLepB, "PartonsLepB");
483  evt.put(pLeptons , "Leptons" );
484  evt.put(pNeutrinos , "Neutrinos" );
485  evt.put(pChi2 , "Chi2" );
486  evt.put(pProb , "Prob" );
487  evt.put(pMT , "MT" );
488  evt.put(pSigMT , "SigMT" );
489  evt.put(pStatus , "Status" );
490  evt.put(pJetsConsidered, "NumberOfConsideredJets");
491 }
492 
493 #endif
int maxNComb_
maximal number of combinations to be written to the event
double maxEtaEle_
maximum eta value for electrons, needed to limited range in which resolutions are provided ...
int i
Definition: DBlmapReader.cc:9
double maxBTagValueNonBJet_
max value of bTag for a non-b-jet
static const unsigned int nPartons
hitfit::JetTranslatorBase< pat::Jet > jetTranslator_
tuple cfg
Definition: looper.py:259
Jets made from CaloTowers.
Definition: CaloJet.h:29
Represent a lepton in an instance of Lepjets_Event class. This class hold the following information: ...
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
hitfit::LeptonTranslatorBase< pat::Electron > electronTranslator_
TtSemiLepHitFitProducer(const edm::ParameterSet &)
edm::EDGetTokenT< std::vector< pat::MET > > metsToken_
hitfit::LeptonTranslatorBase< pat::Muon > muonTranslator_
Lepjets_Event_Jet & jet(std::vector< Lepjets_Event_Jet >::size_type i)
Return a reference to jet at index position i.
double maxEtaJet_
maximum eta value for jets, needed to limited range in which resolutions are provided ...
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
int & type()
Return a reference to the type code.
edm::FileInPath hitfitUdscJetResolution_
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:120
A class to hold functions to calculate kinematic quantities of interest in events.
bool useBTag_
switch to tell whether to use b-tagging or not
vector< PseudoJet > jets
Represent a simple event consisting of lepton(s) and jet(s). An instance of this class holds a list o...
Definition: Lepjets_Event.h:66
tuple result
Definition: query.py:137
hitfit::RunHitFit< pat::Electron, pat::Muon, pat::Jet, pat::MET > PatHitFit
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Template class of function object to translate missing transverse energy object to HitFit&#39;s Fourvec o...
A class to represent a jet in an instance of Lepjets_Event class. The class is derived from the Lepje...
int maxNJets_
maximal number of jets (-1 possible to indicate &#39;all&#39;)
Template class of experiment-independent interface to HitFit.
edm::EDGetTokenT< std::vector< pat::Jet > > jetsToken_
Template class of experiment-independent interface to HitFit. This class is intended to be used insid...
Definition: RunHitFit.h:156
hitfit::METTranslatorBase< pat::MET > metTranslator_
Template class of function object to translate jet physics object to HitFit&#39;s Lepjets_Event_Jet objec...
Fourvec & met()
Return a reference to the missing transverse energy.
Lepjets_Event_Lep & lep(std::vector< Lepjets_Event_Lep >::size_type i)
Return a reference to lepton at index position i.
edm::EDGetTokenT< LeptonCollection > lepsToken_
Template class of function object to translate lepton physics object to HitFit&#39;s Lepjets_Event_Lep ob...
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::FileInPath hitfitElectronResolution_
Analysis-level particle class.
Definition: Particle.h:32
std::string jetCorrectionLevel_
jet correction level
Fourvec & p()
Return a reference to the four-momentum.
virtual void produce(edm::Event &, const edm::EventSetup &)
double maxEtaMu_
maximum eta value for muons, needed to limited range in which resolutions are provided ...
std::string fullPath() const
Definition: FileInPath.cc:165
Definition: Chi2.h:17
double minBTagValueBJet_
min value of bTag for a b-jet
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
std::string bTagAlgo_
input tag for b-tagging algorithm
double jes_
jet energy scale