CMS 3D CMS Logo

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