CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TtFullHadKinFitter.cc
Go to the documentation of this file.
1 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
2 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
3 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
4 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
5 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
6 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
7 
8 #include "AnalysisDataFormats/TopObjects/interface/TtFullHadEvtPartons.h"
9 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullHadKinFitter.h"
10 
11 #include "FWCore/MessageLogger/interface/MessageLogger.h"
12 
13 static const unsigned int nPartons=6;
14 
16 TtFullHadKinFitter::TtFullHadKinFitter():
17  TopKinFitter(),
18  b_(0), bBar_(0), lightQ_(0), lightQBar_(0), lightP_(0), lightPBar_(0),
19  udscResolutions_(0), bResolutions_(0),
20  jetEnergyResolutionScaleFactors_(0), jetEnergyResolutionEtaBinning_(0),
21  jetParam_(kEMom)
22 {
23  setupFitter();
24 }
25 
27 std::vector<TtFullHadKinFitter::Constraint>
28 TtFullHadKinFitter::intToConstraint(std::vector<unsigned int> constraints)
29 {
30  std::vector<TtFullHadKinFitter::Constraint> cConstraints;
31  cConstraints.resize(constraints.size());
32  for(unsigned int i=0;i<constraints.size();++i)
33  {
34  cConstraints[i]=(Constraint)constraints[i];
35  }
36  return cConstraints;
37 }
38 
40 TtFullHadKinFitter::TtFullHadKinFitter(int jetParam, int maxNrIter, double maxDeltaS, double maxF,
41  const std::vector<unsigned int> constraints, double mW, double mTop,
42  const std::vector<edm::ParameterSet>* udscResolutions,
43  const std::vector<edm::ParameterSet>* bResolutions,
44  const std::vector<double>* jetEnergyResolutionScaleFactors,
45  const std::vector<double>* jetEnergyResolutionEtaBinning):
46  TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
47  b_(0), bBar_(0), lightQ_(0), lightQBar_(0), lightP_(0), lightPBar_(0),
48  udscResolutions_(udscResolutions), bResolutions_(bResolutions),
49  jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
50  jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
51  jetParam_((Param)jetParam), constraints_(intToConstraint(constraints))
52 {
53  setupFitter();
54 }
55 
57 TtFullHadKinFitter::TtFullHadKinFitter(Param jetParam, int maxNrIter, double maxDeltaS, double maxF,
58  std::vector<Constraint> constraints, double mW, double mTop,
59  const std::vector<edm::ParameterSet>* udscResolutions,
60  const std::vector<edm::ParameterSet>* bResolutions,
61  const std::vector<double>* jetEnergyResolutionScaleFactors,
62  const std::vector<double>* jetEnergyResolutionEtaBinning):
63  TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
64  b_(0), bBar_(0), lightQ_(0), lightQBar_(0), lightP_(0), lightPBar_(0),
65  udscResolutions_(udscResolutions), bResolutions_(bResolutions),
66  jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
67  jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
68  jetParam_(jetParam), constraints_(constraints)
69 {
70  setupFitter();
71 }
72 
74 TtFullHadKinFitter::~TtFullHadKinFitter()
75 {
76  delete b_;
77  delete bBar_;
78  delete lightQ_;
79  delete lightQBar_;
80  delete lightP_;
81  delete lightPBar_;
82  delete covM_;
83  for(std::map<Constraint, TFitConstraintM*>::iterator it = massConstr_.begin(); it != massConstr_.end(); ++it)
84  delete it->second;
85 }
86 
88 void
89 TtFullHadKinFitter::printSetup() const
90 {
91  std::stringstream constr;
92  for(unsigned int i=0; i<constraints_.size(); ++i){
93  switch(constraints_[i]){
94  case kWPlusMass : constr << " * W+-mass (" << mW_ << " GeV) \n"; break;
95  case kWMinusMass : constr << " * W--mass (" << mW_ << " GeV) \n"; break;
96  case kTopMass : constr << " * t-mass (" << mTop_ << " GeV) \n"; break;
97  case kTopBarMass : constr << " * tBar-mass (" << mTop_ << " GeV) \n"; break;
98  case kEqualTopMasses : constr << " * equal t-masses \n"; break;
99  }
100  }
101  edm::LogVerbatim( "TtFullHadKinFitter" )
102  << "\n"
103  << "+++++++++++ TtFullHadKinFitter Setup ++++++++++++ \n"
104  << " Parametrization: \n"
105  << " * jet : " << param(jetParam_) << "\n"
106  << " Constraints: \n"
107  << constr.str()
108  << " Max(No iterations): " << maxNrIter_ << "\n"
109  << " Max(deltaS) : " << maxDeltaS_ << "\n"
110  << " Max(F) : " << maxF_ << "\n"
111  << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
112 }
113 
115 void
116 TtFullHadKinFitter::setupJets()
117 {
118  TMatrixD empty3x3(3,3);
119  TMatrixD empty4x4(4,4);
120  switch(jetParam_){ // setup jets according to parameterization
121  case kEMom :
122  b_ = new TFitParticleEMomDev ("Jet1", "Jet1", 0, &empty4x4);
123  bBar_ = new TFitParticleEMomDev ("Jet2", "Jet2", 0, &empty4x4);
124  lightQ_ = new TFitParticleEMomDev ("Jet3", "Jet3", 0, &empty4x4);
125  lightQBar_= new TFitParticleEMomDev ("Jet4", "Jet4", 0, &empty4x4);
126  lightP_ = new TFitParticleEMomDev ("Jet5", "Jet5", 0, &empty4x4);
127  lightPBar_= new TFitParticleEMomDev ("Jet6", "Jet6", 0, &empty4x4);
128  break;
129  case kEtEtaPhi :
130  b_ = new TFitParticleEtEtaPhi ("Jet1", "Jet1", 0, &empty3x3);
131  bBar_ = new TFitParticleEtEtaPhi ("Jet2", "Jet2", 0, &empty3x3);
132  lightQ_ = new TFitParticleEtEtaPhi ("Jet3", "Jet3", 0, &empty3x3);
133  lightQBar_= new TFitParticleEtEtaPhi ("Jet4", "Jet4", 0, &empty3x3);
134  lightP_ = new TFitParticleEtEtaPhi ("Jet5", "Jet5", 0, &empty3x3);
135  lightPBar_= new TFitParticleEtEtaPhi ("Jet6", "Jet6", 0, &empty3x3);
136  break;
137  case kEtThetaPhi :
138  b_ = new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3x3);
139  bBar_ = new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3x3);
140  lightQ_ = new TFitParticleEtThetaPhi("Jet3", "Jet3", 0, &empty3x3);
141  lightQBar_= new TFitParticleEtThetaPhi("Jet4", "Jet4", 0, &empty3x3);
142  lightP_ = new TFitParticleEtThetaPhi("Jet5", "Jet5", 0, &empty3x3);
143  lightPBar_= new TFitParticleEtThetaPhi("Jet6", "Jet6", 0, &empty3x3);
144  break;
145  }
146 }
147 
149 void
150 TtFullHadKinFitter::setupConstraints()
151 {
152  massConstr_[kWPlusMass ] = new TFitConstraintM("WPlusMass" , "WPlusMass" , 0, 0, mW_ );
153  massConstr_[kWMinusMass ] = new TFitConstraintM("WMinusMass" , "WMinusMass" , 0, 0, mW_ );
154  massConstr_[kTopMass ] = new TFitConstraintM("TopMass" , "TopMass" , 0, 0, mTop_ );
155  massConstr_[kTopBarMass ] = new TFitConstraintM("TopBarMass" , "TopBarMass" , 0, 0, mTop_ );
156  massConstr_[kEqualTopMasses] = new TFitConstraintM("EqualTopMasses", "EqualTopMasses" , 0, 0, 0 );
157 
158  massConstr_[kWPlusMass ]->addParticles1(lightQ_, lightQBar_);
159  massConstr_[kWMinusMass ]->addParticles1(lightP_, lightPBar_);
160  massConstr_[kTopMass ]->addParticles1(b_, lightQ_, lightQBar_);
161  massConstr_[kTopBarMass ]->addParticles1(bBar_, lightP_, lightPBar_);
162  massConstr_[kEqualTopMasses]->addParticles1(b_, lightQ_, lightQBar_);
163  massConstr_[kEqualTopMasses]->addParticles2(bBar_, lightP_, lightPBar_);
164 
165 }
166 
168 void
169 TtFullHadKinFitter::setupFitter()
170 {
171  printSetup();
172  setupJets();
173  setupConstraints();
174 
175  // add measured particles
176  fitter_->addMeasParticle(b_);
177  fitter_->addMeasParticle(bBar_);
178  fitter_->addMeasParticle(lightQ_);
179  fitter_->addMeasParticle(lightQBar_);
180  fitter_->addMeasParticle(lightP_);
181  fitter_->addMeasParticle(lightPBar_);
182 
183  // add constraints
184  for(unsigned int i=0; i<constraints_.size(); i++){
185  fitter_->addConstraint(massConstr_[constraints_[i]]);
186  }
187 
188  // initialize helper class used to bring the resolutions into covariance matrices
189  if(udscResolutions_->size() && bResolutions_->size())
190  covM_ = new CovarianceMatrix(*udscResolutions_, *bResolutions_,
191  *jetEnergyResolutionScaleFactors_,
192  *jetEnergyResolutionEtaBinning_);
193  else
194  covM_ = new CovarianceMatrix();
195 }
196 
198 int
199 TtFullHadKinFitter::fit(const std::vector<pat::Jet>& jets)
200 {
201  if( jets.size()<6 ){
202  throw edm::Exception( edm::errors::Configuration, "Cannot run the TtFullHadKinFitter with less than 6 jets" );
203  }
204 
205  // get jets in right order
206  const pat::Jet& b = jets[TtFullHadEvtPartons::B ];
207  const pat::Jet& bBar = jets[TtFullHadEvtPartons::BBar ];
208  const pat::Jet& lightQ = jets[TtFullHadEvtPartons::LightQ ];
209  const pat::Jet& lightQBar = jets[TtFullHadEvtPartons::LightQBar];
210  const pat::Jet& lightP = jets[TtFullHadEvtPartons::LightP ];
211  const pat::Jet& lightPBar = jets[TtFullHadEvtPartons::LightPBar];
212 
213  // initialize particles
214  const TLorentzVector p4B( b.px(), b.py(), b.pz(), b.energy() );
215  const TLorentzVector p4BBar( bBar.px(), bBar.py(), bBar.pz(), bBar.energy() );
216  const TLorentzVector p4LightQ( lightQ.px(), lightQ.py(), lightQ.pz(), lightQ.energy() );
217  const TLorentzVector p4LightQBar( lightQBar.px(), lightQBar.py(), lightQBar.pz(), lightQBar.energy() );
218  const TLorentzVector p4LightP( lightP.px(), lightP.py(), lightP.pz(), lightP.energy() );
219  const TLorentzVector p4LightPBar( lightPBar.px(), lightPBar.py(), lightPBar.pz(), lightPBar.energy() );
220 
221  // initialize covariance matrices
222  TMatrixD m1 = covM_->setupMatrix(lightQ, jetParam_);
223  TMatrixD m2 = covM_->setupMatrix(lightQBar, jetParam_);
224  TMatrixD m3 = covM_->setupMatrix(b, jetParam_, "bjets");
225  TMatrixD m4 = covM_->setupMatrix(lightP, jetParam_);
226  TMatrixD m5 = covM_->setupMatrix(lightPBar, jetParam_);
227  TMatrixD m6 = covM_->setupMatrix(bBar , jetParam_, "bjets");
228 
229  // set the kinematics of the objects to be fitted
230  b_ ->setIni4Vec(&p4B );
231  bBar_ ->setIni4Vec(&p4BBar );
232  lightQ_ ->setIni4Vec(&p4LightQ );
233  lightQBar_->setIni4Vec(&p4LightQBar);
234  lightP_ ->setIni4Vec(&p4LightP );
235  lightPBar_->setIni4Vec(&p4LightPBar);
236 
237  // initialize covariance matrices
238  lightQ_ ->setCovMatrix( &m1);
239  lightQBar_->setCovMatrix( &m2);
240  b_ ->setCovMatrix( &m3);
241  lightP_ ->setCovMatrix( &m4);
242  lightPBar_->setCovMatrix( &m5);
243  bBar_ ->setCovMatrix( &m6);
244 
245  // perform the fit!
246  fitter_->fit();
247 
248  // add fitted information to the solution
249  if( fitter_->getStatus()==0 ){
250  // read back jet kinematics
251  fittedB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(b_->getCurr4Vec()->X(), b_->getCurr4Vec()->Y(), b_->getCurr4Vec()->Z(), b_->getCurr4Vec()->E()), math::XYZPoint()));
252  fittedLightQ_ = pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightQ_->getCurr4Vec()->X(), lightQ_->getCurr4Vec()->Y(), lightQ_->getCurr4Vec()->Z(), lightQ_->getCurr4Vec()->E()), math::XYZPoint()));
253  fittedLightQBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightQBar_->getCurr4Vec()->X(), lightQBar_->getCurr4Vec()->Y(), lightQBar_->getCurr4Vec()->Z(), lightQBar_->getCurr4Vec()->E()), math::XYZPoint()));
254 
255 
256  fittedBBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(bBar_->getCurr4Vec()->X(), bBar_->getCurr4Vec()->Y(), bBar_->getCurr4Vec()->Z(), bBar_->getCurr4Vec()->E()), math::XYZPoint()));
257  fittedLightP_ = pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightP_->getCurr4Vec()->X(), lightP_->getCurr4Vec()->Y(), lightP_->getCurr4Vec()->Z(), lightP_->getCurr4Vec()->E()), math::XYZPoint()));
258  fittedLightPBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightPBar_->getCurr4Vec()->X(), lightPBar_->getCurr4Vec()->Y(), lightPBar_->getCurr4Vec()->Z(), lightPBar_->getCurr4Vec()->E()), math::XYZPoint()));
259  }
260  return fitter_->getStatus();
261 }
262 
264 TtHadEvtSolution
265 TtFullHadKinFitter::addKinFitInfo(TtHadEvtSolution * asol)
266 {
267  TtHadEvtSolution fitsol(*asol);
268 
269  std::vector<pat::Jet> jets;
270  jets.resize(6);
271  jets[TtFullHadEvtPartons::LightQ ] = fitsol.getCalHadp();
272  jets[TtFullHadEvtPartons::LightQBar] = fitsol.getCalHadq();
273  jets[TtFullHadEvtPartons::B ] = fitsol.getCalHadb();
274  jets[TtFullHadEvtPartons::LightP ] = fitsol.getCalHadj();
275  jets[TtFullHadEvtPartons::LightPBar] = fitsol.getCalHadk();
276  jets[TtFullHadEvtPartons::BBar ] = fitsol.getCalHadbbar();
277 
278  fit( jets);
279 
280  // add fitted information to the solution
281  if (fitter_->getStatus() == 0) {
282  // finally fill the fitted particles
283  fitsol.setFitHadb(fittedB_);
284  fitsol.setFitHadp(fittedLightQ_);
285  fitsol.setFitHadq(fittedLightQBar_);
286  fitsol.setFitHadk(fittedLightP_);
287  fitsol.setFitHadj(fittedLightPBar_);
288  fitsol.setFitHadbbar(fittedBBar_);
289 
290  // store the fit's chi2 probability
291  fitsol.setProbChi2( fitProb() );
292  }
293  return fitsol;
294 }
295 
296 
298 TtFullHadKinFitter::KinFit::KinFit() :
299  useBTagging_(true),
300  bTags_(2),
301  bTagAlgo_("trackCountingHighPurBJetTags"),
302  minBTagValueBJet_(3.41),
303  maxBTagValueNonBJet_(3.41),
304  udscResolutions_(std::vector<edm::ParameterSet>(0)),
305  bResolutions_(std::vector<edm::ParameterSet>(0)),
306  jetEnergyResolutionScaleFactors_(0),
307  jetEnergyResolutionEtaBinning_(0),
308  jetCorrectionLevel_("L3Absolute"),
309  maxNJets_(-1),
310  maxNComb_(1),
311  maxNrIter_(500),
312  maxDeltaS_(5e-5),
313  maxF_(0.0001),
314  jetParam_(1),
315  mW_(80.4),
316  mTop_(173.),
317  useOnlyMatch_(false),
318  match_(std::vector<int>(0)),
319  invalidMatch_(false)
320 {
321  constraints_.push_back(1);
322  constraints_.push_back(2);
323  constraints_.push_back(5);
324 }
325 
327 TtFullHadKinFitter::KinFit::KinFit(bool useBTagging, unsigned int bTags, std::string bTagAlgo, double minBTagValueBJet, double maxBTagValueNonBJet,
328  std::vector<edm::ParameterSet> udscResolutions, std::vector<edm::ParameterSet> bResolutions,
329  std::vector<double> jetEnergyResolutionScaleFactors, std::vector<double> jetEnergyResolutionEtaBinning,
330  std::string jetCorrectionLevel, int maxNJets, int maxNComb,
331  unsigned int maxNrIter, double maxDeltaS, double maxF, unsigned int jetParam, std::vector<unsigned> constraints, double mW, double mTop) :
332  useBTagging_(useBTagging),
333  bTags_(bTags),
334  bTagAlgo_(bTagAlgo),
335  minBTagValueBJet_(minBTagValueBJet),
336  maxBTagValueNonBJet_(maxBTagValueNonBJet),
337  udscResolutions_(udscResolutions),
338  bResolutions_(bResolutions),
339  jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
340  jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
341  jetCorrectionLevel_(jetCorrectionLevel),
342  maxNJets_(maxNJets),
343  maxNComb_(maxNComb),
344  maxNrIter_(maxNrIter),
345  maxDeltaS_(maxDeltaS),
346  maxF_(maxF),
347  jetParam_(jetParam),
348  constraints_(constraints),
349  mW_(mW),
350  mTop_(mTop),
351  useOnlyMatch_(false),
352  invalidMatch_(false)
353 {
354  // define kinematic fit interface
355  fitter = new TtFullHadKinFitter(param(jetParam_), maxNrIter_, maxDeltaS_, maxF_, TtFullHadKinFitter::KinFit::constraints(constraints_), mW_, mTop_,
356  &udscResolutions_, &bResolutions_, &jetEnergyResolutionScaleFactors_, &jetEnergyResolutionEtaBinning_);
357 }
358 
360 TtFullHadKinFitter::KinFit::~KinFit()
361 {
362  delete fitter;
363 }
364 
365 bool
366 TtFullHadKinFitter::KinFit::doBTagging(const std::vector<pat::Jet>& jets, const unsigned int& bJetCounter, std::vector<int>& combi){
367 
368  if( !useBTagging_ ) {
369  return true;
370  }
371  if( bTags_ == 2 &&
372  jets[combi[TtFullHadEvtPartons::B ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
373  jets[combi[TtFullHadEvtPartons::BBar ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
374  jets[combi[TtFullHadEvtPartons::LightQ ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
375  jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
376  jets[combi[TtFullHadEvtPartons::LightP ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
377  jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ ) {
378  return true;
379  }
380  else if( bTags_ == 1 ){
381  if( bJetCounter == 1 &&
382  (jets[combi[TtFullHadEvtPartons::B ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
383  jets[combi[TtFullHadEvtPartons::BBar ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
384  jets[combi[TtFullHadEvtPartons::LightQ ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
385  jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
386  jets[combi[TtFullHadEvtPartons::LightP ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
387  jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ ) {
388  return true;
389  }
390  else if( bJetCounter > 1 &&
391  jets[combi[TtFullHadEvtPartons::B ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
392  jets[combi[TtFullHadEvtPartons::BBar ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
393  jets[combi[TtFullHadEvtPartons::LightQ ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
394  jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
395  jets[combi[TtFullHadEvtPartons::LightP ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
396  jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ ) {
397  return true;
398  }
399  }
400  else if( bTags_ == 0 ){
401  if( bJetCounter == 0){
402  return true;
403  }
404  else if( bJetCounter == 1 &&
405  (jets[combi[TtFullHadEvtPartons::B ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
406  jets[combi[TtFullHadEvtPartons::BBar ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
407  jets[combi[TtFullHadEvtPartons::LightQ ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
408  jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
409  jets[combi[TtFullHadEvtPartons::LightP ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
410  jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ ) {
411  return true;
412  }
413  else if( bJetCounter > 1 &&
414  jets[combi[TtFullHadEvtPartons::B ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
415  jets[combi[TtFullHadEvtPartons::BBar ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
416  jets[combi[TtFullHadEvtPartons::LightQ ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
417  jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
418  jets[combi[TtFullHadEvtPartons::LightP ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
419  jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ ) {
420  return true;
421  }
422  }
423  else if( bTags_ > 2 ){
424  throw cms::Exception("Configuration")
425  << "Wrong number of bTags (" << bTags_ << " bTags not supported)!\n";
426  return true;
427  }
428  return false;
429 }
430 
432 pat::Jet
433 TtFullHadKinFitter::KinFit::corJet(const pat::Jet& jet, const std::string& quarkType)
434 {
435  // jetCorrectionLevel was not configured
436  if(jetCorrectionLevel_.empty())
437  throw cms::Exception("Configuration")
438  << "Unconfigured jetCorrectionLevel. Please use an appropriate, non-empty string.\n";
439 
440  // quarkType is unknown
441  if( !(quarkType=="wMix" ||
442  quarkType=="uds" ||
443  quarkType=="charm" ||
444  quarkType=="bottom") )
445  throw cms::Exception("Configuration")
446  << quarkType << " is unknown as a quarkType for the jetCorrectionLevel.\n";
447 
448  float jecFactor = 1.;
449  if(quarkType=="wMix") jecFactor = 0.75 * jet.jecFactor(jetCorrectionLevel_, "uds") + 0.25 * jet.jecFactor(jetCorrectionLevel_, "charm");
450  else jecFactor = jet.jecFactor(jetCorrectionLevel_, quarkType);
451 
452  pat::Jet ret = jet;
453  ret.setP4(ret.p4()*jecFactor);
454  return ret;
455 }
456 
457 std::list<TtFullHadKinFitter::KinFitResult>
458 TtFullHadKinFitter::KinFit::fit(const std::vector<pat::Jet>& jets){
459 
460  std::list<TtFullHadKinFitter::KinFitResult> fitResults;
461 
468  if( jets.size()<nPartons || invalidMatch_ ) {
469  // indices referring to the jet combination
470  std::vector<int> invalidCombi;
471  for(unsigned int i = 0; i < nPartons; ++i) invalidCombi.push_back( -1 );
472 
473  KinFitResult result;
474  // status of the fitter
475  result.Status = -1;
476  // chi2
477  result.Chi2 = -1.;
478  // chi2 probability
479  result.Prob = -1.;
480  // the kinFit getters return empty objects here
481  result.B = fitter->fittedB();
482  result.BBar = fitter->fittedBBar();
483  result.LightQ = fitter->fittedLightQ();
484  result.LightQBar= fitter->fittedLightQBar();
485  result.LightP = fitter->fittedLightP();
486  result.LightPBar= fitter->fittedLightPBar();
487  result.JetCombi = invalidCombi;
488  // push back fit result
489  fitResults.push_back( result );
490  return fitResults;
491  }
492 
498  std::vector<int> jetIndices;
499  if(!useOnlyMatch_) {
500  for(unsigned int idx=0; idx<jets.size(); ++idx){
501  if(maxNJets_>=(int)nPartons && maxNJets_==(int)idx) break;
502  jetIndices.push_back(idx);
503  }
504  }
505 
506  std::vector<int> combi;
507  for(unsigned int idx=0; idx<nPartons; ++idx) {
508  useOnlyMatch_?combi.push_back(match_[idx]):combi.push_back(idx);
509  }
510 
511 
512  unsigned int bJetCounter = 0;
513  for(std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet < jets.end(); ++jet){
514  if(jet->bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) ++bJetCounter;
515  }
516 
517  do{
518  for(int cnt=0; cnt<TMath::Factorial(combi.size()); ++cnt){
519  // take into account indistinguishability of the two jets from the two W decays,
520  // and the two decay branches, this reduces the combinatorics by a factor of 2*2*2
521  if( ((combi[TtFullHadEvtPartons::LightQ] < combi[TtFullHadEvtPartons::LightQBar] &&
522  combi[TtFullHadEvtPartons::LightP] < combi[TtFullHadEvtPartons::LightPBar] &&
523  combi[TtFullHadEvtPartons::B] < combi[TtFullHadEvtPartons::BBar] ) ||
524  useOnlyMatch_) && doBTagging(jets, bJetCounter, combi) ) {
525 
526  std::vector<pat::Jet> jetCombi;
527  jetCombi.resize(nPartons);
528  jetCombi[TtFullHadEvtPartons::LightQ ] = corJet(jets[combi[TtFullHadEvtPartons::LightQ ]], "wMix");
529  jetCombi[TtFullHadEvtPartons::LightQBar] = corJet(jets[combi[TtFullHadEvtPartons::LightQBar]], "wMix");
530  jetCombi[TtFullHadEvtPartons::B ] = corJet(jets[combi[TtFullHadEvtPartons::B ]], "bottom");
531  jetCombi[TtFullHadEvtPartons::BBar ] = corJet(jets[combi[TtFullHadEvtPartons::BBar ]], "bottom");
532  jetCombi[TtFullHadEvtPartons::LightP ] = corJet(jets[combi[TtFullHadEvtPartons::LightP ]], "wMix");
533  jetCombi[TtFullHadEvtPartons::LightPBar] = corJet(jets[combi[TtFullHadEvtPartons::LightPBar]], "wMix");
534 
535  // do the kinematic fit
536  int status = fitter->fit(jetCombi);
537 
538  if( status == 0 ) {
539  // fill struct KinFitResults if converged
540  TtFullHadKinFitter::KinFitResult result;
541  result.Status = status;
542  result.Chi2 = fitter->fitS();
543  result.Prob = fitter->fitProb();
544  result.B = fitter->fittedB();
545  result.BBar = fitter->fittedBBar();
546  result.LightQ = fitter->fittedLightQ();
547  result.LightQBar= fitter->fittedLightQBar();
548  result.LightP = fitter->fittedLightP();
549  result.LightPBar= fitter->fittedLightPBar();
550  result.JetCombi = combi;
551  // push back fit result
552  fitResults.push_back( result );
553  }
554  }
555  // don't go through combinatorics if useOnlyMatch was chosen
556  if(useOnlyMatch_){
557  break;
558  }
559  // next permutation
560  std::next_permutation( combi.begin(), combi.end() );
561  }
562  // don't go through combinatorics if useOnlyMatch was chosen
563  if(useOnlyMatch_){
564  break;
565  }
566  }
567  while( stdcomb::next_combination( jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end() ) );
568 
569 
570  // sort results w.r.t. chi2 values
571  fitResults.sort();
572 
578  if( (unsigned)fitResults.size() < 1 ) {
579  // in case no fit results were stored in the list (i.e. when all fits were aborted)
580 
581  KinFitResult result;
582  // status of the fitter
583  result.Status = -1;
584  // chi2
585  result.Chi2 = -1.;
586  // chi2 probability
587  result.Prob = -1.;
588  // the kinFit getters return empty objects here
589  result.B = fitter->fittedB();
590  result.BBar = fitter->fittedBBar();
591  result.LightQ = fitter->fittedLightQ();
592  result.LightQBar= fitter->fittedLightQBar();
593  result.LightP = fitter->fittedLightP();
594  result.LightPBar= fitter->fittedLightPBar();
595  // indices referring to the jet combination
596  std::vector<int> invalidCombi(nPartons, -1);
597  result.JetCombi = invalidCombi;
598  // push back fit result
599  fitResults.push_back( result );
600  }
601  return fitResults;
602 }
603 
604 TtFullHadKinFitter::Param
605 TtFullHadKinFitter::KinFit::param(unsigned int configParameter)
606 {
607  TtFullHadKinFitter::Param result;
608  switch(configParameter){
609  case TtFullHadKinFitter::kEMom : result=TtFullHadKinFitter::kEMom; break;
610  case TtFullHadKinFitter::kEtEtaPhi : result=TtFullHadKinFitter::kEtEtaPhi; break;
611  case TtFullHadKinFitter::kEtThetaPhi : result=TtFullHadKinFitter::kEtThetaPhi; break;
612  default:
613  throw cms::Exception("Configuration")
614  << "Chosen jet parametrization is not supported: " << configParameter << "\n";
615  break;
616  }
617  return result;
618 }
619 
620 TtFullHadKinFitter::Constraint
621 TtFullHadKinFitter::KinFit::constraint(unsigned configParameter)
622 {
623  TtFullHadKinFitter::Constraint result;
624  switch(configParameter){
625  case TtFullHadKinFitter::kWPlusMass : result=TtFullHadKinFitter::kWPlusMass; break;
626  case TtFullHadKinFitter::kWMinusMass : result=TtFullHadKinFitter::kWMinusMass; break;
627  case TtFullHadKinFitter::kTopMass : result=TtFullHadKinFitter::kTopMass; break;
628  case TtFullHadKinFitter::kTopBarMass : result=TtFullHadKinFitter::kTopBarMass; break;
629  case TtFullHadKinFitter::kEqualTopMasses : result=TtFullHadKinFitter::kEqualTopMasses; break;
630  default:
631  throw cms::Exception("Configuration")
632  << "Chosen fit constraint is not supported: " << configParameter << "\n";
633  break;
634  }
635  return result;
636 }
637 
638 std::vector<TtFullHadKinFitter::Constraint>
639 TtFullHadKinFitter::KinFit::constraints(std::vector<unsigned>& configParameters)
640 {
641  std::vector<TtFullHadKinFitter::Constraint> result;
642  for(unsigned i=0; i<configParameters.size(); ++i){
643  result.push_back(constraint(configParameters[i]));
644  }
645  return result;
646 }
TtFullHadKinFitter::lightQ_
TAbsFitParticle * lightQ_
Definition: TtFullHadKinFitter.h:88
TtFullHadKinFitter::KinFit::maxF_
double maxF_
maximal deviation for contstraints
Definition: TtFullHadKinFitter.h:241
pat::Jet::jecFactor
float jecFactor(const std::string &level, const std::string &flavor="none", const std::string &set="") const
Definition: Jet.cc:239
i
int i
Definition: DBlmapReader.cc:9
TtFullHadKinFitter::jetEnergyResolutionEtaBinning_
const std::vector< double > * jetEnergyResolutionEtaBinning_
Definition: TtFullHadKinFitter.h:97
TtFullHadKinFitter::KinFit::maxDeltaS_
double maxDeltaS_
maximal chi2 equivalent
Definition: TtFullHadKinFitter.h:239
TtFullHadKinFitter::KinFit::param
TtFullHadKinFitter::Param param(unsigned int configParameter)
Definition: TtFullHadKinFitter.cc:605
TtHadEvtSolution::getCalHadj
pat::Jet getCalHadj() const
Definition: TtHadEvtSolution.h:82
TtHadEvtSolution::setFitHadj
void setFitHadj(const pat::Particle &aFitHadj)
Definition: TtHadEvtSolution.h:178
TtFullHadKinFitter::covM_
CovarianceMatrix * covM_
get object resolutions and put them into a matrix
Definition: TtFullHadKinFitter.h:113
TopKinFitter::Param
Param
supported parameterizations
Definition: TopKinFitter.h:22
TtFullHadKinFitter::setupConstraints
void setupConstraints()
initialize constraints
Definition: TtFullHadKinFitter.cc:150
TtFullHadKinFitter::addKinFitInfo
TtHadEvtSolution addKinFitInfo(TtHadEvtSolution *asol)
add kin fit information to the old event solution (in for legacy reasons)
Definition: TtFullHadKinFitter.cc:265
TtFullHadKinFitter::lightP_
TAbsFitParticle * lightP_
Definition: TtFullHadKinFitter.h:90
TKinFitter::fit
Int_t fit()
Definition: TKinFitter.cc:309
TtFullHadKinFitter::KinFit::constraint
TtFullHadKinFitter::Constraint constraint(unsigned int configParameter)
Definition: TtFullHadKinFitter.cc:621
TAbsFitParticle::setIni4Vec
virtual void setIni4Vec(const TLorentzVector *pini)=0
run_regression.ret
int ret
Definition: run_regression.py:388
TtHadEvtSolution::setFitHadb
void setFitHadb(const pat::Particle &aFitHadb)
Definition: TtHadEvtSolution.h:180
TtHadEvtSolution::getCalHadq
pat::Jet getCalHadq() const
Definition: TtHadEvtSolution.h:81
TtFullHadKinFitter::fittedLightPBar_
pat::Particle fittedLightPBar_
Definition: TtFullHadKinFitter.h:106
TopKinFitter::param
std::string param(const Param &param) const
convert Param to human readable form
Definition: TopKinFitter.cc:23
TtHadEvtSolution::getCalHadp
pat::Jet getCalHadp() const
Definition: TtHadEvtSolution.h:80
TtHadEvtSolution::setFitHadp
void setFitHadp(const pat::Particle &aFitHadp)
Definition: TtHadEvtSolution.h:176
reco::LeafCandidate::setP4
virtual void setP4(const LorentzVector &p4)
set 4-momentum
Definition: LeafCandidate.h:149
TtFullHadKinFitter::fittedB_
pat::Particle fittedB_
output particles
Definition: TtFullHadKinFitter.h:101
CovarianceMatrix::setupMatrix
TMatrixD setupMatrix(const pat::PATObject< T > &object, const TopKinFitter::Param param, const std::string &resolutionProvider="")
return covariance matrix for a PAT object
Definition: CovarianceMatrix.h:81
TtFullHadKinFitter::fit
int fit(const std::vector< pat::Jet > &jets)
kinematic fit interface
Definition: TtFullHadKinFitter.cc:199
TtFullHadKinFitter::printSetup
void printSetup() const
print fitter setup
Definition: TtFullHadKinFitter.cc:89
TtFullHadKinFitter::KinFit::fitter
TtFullHadKinFitter * fitter
kinematic fit interface
Definition: TtFullHadKinFitter.h:258
TtFullHadKinFitter::KinFit::mW_
double mW_
W mass value used for constraints.
Definition: TtFullHadKinFitter.h:247
TtFullHadKinFitter::KinFit::constraints_
std::vector< unsigned > constraints_
numbering of different possible kinematic constraints
Definition: TtFullHadKinFitter.h:245
TtFullHadKinFitter::fittedLightQBar_
pat::Particle fittedLightQBar_
Definition: TtFullHadKinFitter.h:104
TtFullHadKinFitter::jetEnergyResolutionScaleFactors_
const std::vector< double > * jetEnergyResolutionScaleFactors_
scale factors for the jet energy resolution
Definition: TtFullHadKinFitter.h:96
TtFullHadKinFitter::KinFit::fit
std::list< TtFullHadKinFitter::KinFitResult > fit(const std::vector< pat::Jet > &jets)
do the fitting and return fit result
Definition: TtFullHadKinFitter.cc:458
TopKinFitter::maxNrIter_
int maxNrIter_
maximal allowed number of iterations to be used for the fit
Definition: TopKinFitter.h:48
TtHadEvtSolution::setProbChi2
void setProbChi2(double c)
Definition: TtHadEvtSolution.h:206
TtFullHadKinFitter::KinFit::jetEnergyResolutionScaleFactors_
std::vector< double > jetEnergyResolutionScaleFactors_
scale factors for the jet energy resolution
Definition: TtFullHadKinFitter.h:228
TtFullHadKinFitter::KinFit::udscResolutions_
std::vector< edm::ParameterSet > udscResolutions_
store the resolutions for the jets
Definition: TtFullHadKinFitter.h:226
TopKinFitter::fitProb
double fitProb() const
return fit probability
Definition: TopKinFitter.h:36
TtFullHadKinFitter::KinFit::bResolutions_
std::vector< edm::ParameterSet > bResolutions_
Definition: TtFullHadKinFitter.h:226
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
TtFullHadKinFitter::KinFit::doBTagging
bool doBTagging(const std::vector< pat::Jet > &jets, const unsigned int &bJetCounter, std::vector< int > &combi)
Definition: TtFullHadKinFitter.cc:366
TtFullHadKinFitter::KinFit::corJet
pat::Jet corJet(const pat::Jet &jet, const std::string &quarkType)
helper function to construct the proper corrected jet for its corresponding quarkType ...
Definition: TtFullHadKinFitter.cc:433
TtFullHadKinFitter::KinFit::jetEnergyResolutionEtaBinning_
std::vector< double > jetEnergyResolutionEtaBinning_
Definition: TtFullHadKinFitter.h:229
reco::LeafCandidate::energy
virtual double energy() const
energy
Definition: LeafCandidate.h:119
TtFullHadKinFitter::Constraint
Constraint
supported constraints
Definition: TtFullHadKinFitter.h:33
TtFullHadKinFitter::KinFit::maxNrIter_
unsigned int maxNrIter_
maximal number of iterations to be performed for the fit
Definition: TtFullHadKinFitter.h:237
TtFullHadKinFitter::~TtFullHadKinFitter
~TtFullHadKinFitter()
default destructor
Definition: TtFullHadKinFitter.cc:74
TtFullHadKinFitter::intToConstraint
std::vector< TtFullHadKinFitter::Constraint > intToConstraint(std::vector< unsigned int > constraints)
used to convert vector of int&#39;s to vector of constraints (just used in TtFullHadKinFitter(int, int, double, double, std::vector&lt;unsigned int&gt;))
Definition: TtFullHadKinFitter.cc:28
TtFullHadKinFitter::bResolutions_
const std::vector< edm::ParameterSet > * bResolutions_
Definition: TtFullHadKinFitter.h:94
TKinFitter::addConstraint
void addConstraint(TAbsFitConstraint *constraint)
Definition: TKinFitter.cc:281
fwrapper::jets
vector< PseudoJet > jets
Definition: fastjetfortran_madfks.cc:43
query.result
tuple result
Definition: query.py:137
TtFullHadKinFitter::b_
TAbsFitParticle * b_
input particles
Definition: TtFullHadKinFitter.h:86
TtHadEvtSolution::getCalHadb
pat::Jet getCalHadb() const
Definition: TtHadEvtSolution.h:78
TKinFitter::getStatus
Int_t getStatus()
Definition: TKinFitter.h:41
TAbsFitParticle::setCovMatrix
virtual void setCovMatrix(const TMatrixD *theCovMatrix)
Definition: TAbsFitParticle.cc:117
TtHadEvtSolution::getCalHadk
pat::Jet getCalHadk() const
Definition: TtHadEvtSolution.h:83
TopKinFitter::mW_
double mW_
W mass value used for constraints.
Definition: TopKinFitter.h:54
TKinFitter::addMeasParticle
void addMeasParticle(TAbsFitParticle *particle)
Definition: TKinFitter.cc:209
TtFullHadKinFitter::massConstr_
std::map< Constraint, TFitConstraintM * > massConstr_
supported constraints
Definition: TtFullHadKinFitter.h:99
TtFullHadKinFitter::fittedLightP_
pat::Particle fittedLightP_
Definition: TtFullHadKinFitter.h:105
TtFullHadKinFitter::lightQBar_
TAbsFitParticle * lightQBar_
Definition: TtFullHadKinFitter.h:89
TtFullHadKinFitter::fittedBBar_
pat::Particle fittedBBar_
Definition: TtFullHadKinFitter.h:102
TopKinFitter::maxDeltaS_
double maxDeltaS_
maximal allowed chi2 (not normalized to degrees of freedom)
Definition: TopKinFitter.h:50
TtFullHadKinFitter::KinFit::mTop_
double mTop_
top mass value used for constraints
Definition: TtFullHadKinFitter.h:249
TtFullHadKinFitter::constraints_
std::vector< Constraint > constraints_
vector of constraints to be used
Definition: TtFullHadKinFitter.h:110
TtFullHadKinFitter::jetParam_
Param jetParam_
jet parametrization
Definition: TtFullHadKinFitter.h:108
TtFullHadKinFitter::fittedLightQ_
pat::Particle fittedLightQ_
Definition: TtFullHadKinFitter.h:103
TtFullHadKinFitter::setupFitter
void setupFitter()
setup fitter
Definition: TtFullHadKinFitter.cc:169
reco::LeafCandidate::px
virtual double px() const
x coordinate of momentum vector
Definition: LeafCandidate.h:131
TtFullHadKinFitter::bBar_
TAbsFitParticle * bBar_
Definition: TtFullHadKinFitter.h:87
TtFullHadKinFitter::KinFit::~KinFit
~KinFit()
default destructor
Definition: TtFullHadKinFitter.cc:360
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
pat::Particle
Analysis-level particle class.
Definition: Particle.h:34
TtHadEvtSolution::setFitHadbbar
void setFitHadbbar(const pat::Particle &aFitHadbbar)
Definition: TtHadEvtSolution.h:181
b
double b
Definition: hdecay.h:120
TtFullHadKinFitter::udscResolutions_
const std::vector< edm::ParameterSet > * udscResolutions_
resolutions
Definition: TtFullHadKinFitter.h:93
pat::Jet
Analysis-level calorimeter jet class.
Definition: Jet.h:71
TtHadEvtSolution::setFitHadk
void setFitHadk(const pat::Particle &aFitHadk)
Definition: TtHadEvtSolution.h:179
nPartons
static const unsigned int nPartons
Definition: TtFullHadKinFitter.cc:13
reco::LeafCandidate::pz
virtual double pz() const
z coordinate of momentum vector
Definition: LeafCandidate.h:135
TAbsFitParticle::getCurr4Vec
const TLorentzVector * getCurr4Vec()
Definition: TAbsFitParticle.h:33
TtFullHadKinFitter::KinFit::jetParam_
unsigned int jetParam_
numbering of different possible jet parametrizations
Definition: TtFullHadKinFitter.h:243
TtHadEvtSolution::setFitHadq
void setFitHadq(const pat::Particle &aFitHadq)
Definition: TtHadEvtSolution.h:177
TtFullHadKinFitter::TtFullHadKinFitter
TtFullHadKinFitter()
default constructor
Definition: TtFullHadKinFitter.cc:16
TtFullHadKinFitter::setupJets
void setupJets()
initialize jet inputs
Definition: TtFullHadKinFitter.cc:116
TopKinFitter::fitter_
TKinFitter * fitter_
kinematic fitter
Definition: TopKinFitter.h:46
TtFullHadKinFitter::KinFit::constraints
std::vector< TtFullHadKinFitter::Constraint > constraints(std::vector< unsigned int > &configParameters)
Definition: TtFullHadKinFitter.cc:639
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
stdcomb::next_combination
bool next_combination(BidIt n_begin, BidIt n_end, BidIt r_begin, BidIt r_end)
Definition: combination.h:22
reco::LeafCandidate::p4
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
Definition: LeafCandidate.h:108
TtHadEvtSolution::getCalHadbbar
pat::Jet getCalHadbbar() const
Definition: TtHadEvtSolution.h:79
TopKinFitter::mTop_
double mTop_
top mass value used for constraints
Definition: TopKinFitter.h:56
TtFullHadKinFitter::KinFit::KinFit
KinFit()
default constructor
Definition: TtFullHadKinFitter.cc:298
reco::LeafCandidate::py
virtual double py() const
y coordinate of momentum vector
Definition: LeafCandidate.h:133
TtFullHadKinFitter::lightPBar_
TAbsFitParticle * lightPBar_
Definition: TtFullHadKinFitter.h:91
TopKinFitter::maxF_
double maxF_
maximal allowed distance from constraints
Definition: TopKinFitter.h:52
Generated for CMSSW Reference Manual by doxygen 1.8.5