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.
7 
10 
12 
13 static const unsigned int nPartons=6;
14 
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(const 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  const 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 
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
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
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
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_);
164 
165 }
166 
168 void
170 {
171  printSetup();
172  setupJets();
174 
175  // add measured particles
182 
183  // add constraints
184  for(unsigned int i=0; i<constraints_.size(); i++){
186  }
187 
188  // initialize helper class used to bring the resolutions into covariance matrices
189  if(udscResolutions_->size() && bResolutions_->size())
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
254 
255 
259  }
260  return fitter_->getStatus();
261 }
262 
266 {
267  TtHadEvtSolution fitsol(*asol);
268 
269  std::vector<pat::Jet> jets;
270  jets.resize(6);
271  jets[TtFullHadEvtPartons::LightQ ] = fitsol.getCalHadp();
273  jets[TtFullHadEvtPartons::B ] = fitsol.getCalHadb();
274  jets[TtFullHadEvtPartons::LightP ] = fitsol.getCalHadj();
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_);
286  fitsol.setFitHadk(fittedLightP_);
288  fitsol.setFitHadbbar(fittedBBar_);
289 
290  // store the fit's chi2 probability
291  fitsol.setProbChi2( fitProb() );
292  }
293  return fitsol;
294 }
295 
296 
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  const std::vector<edm::ParameterSet>& udscResolutions, const std::vector<edm::ParameterSet>& bResolutions,
329  const std::vector<double>& jetEnergyResolutionScaleFactors, const std::vector<double>& jetEnergyResolutionEtaBinning,
330  std::string jetCorrectionLevel, int maxNJets, int maxNComb,
331  unsigned int maxNrIter, double maxDeltaS, double maxF, unsigned int jetParam, const 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
357 }
358 
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
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 
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
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
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 
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 
605 TtFullHadKinFitter::KinFit::param(unsigned int configParameter)
606 {
608  switch(configParameter){
612  default:
613  throw cms::Exception("Configuration")
614  << "Chosen jet parametrization is not supported: " << configParameter << "\n";
615  break;
616  }
617  return result;
618 }
619 
621 TtFullHadKinFitter::KinFit::constraint(unsigned configParameter)
622 {
624  switch(configParameter){
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(const 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 }
TAbsFitParticle * lightQ_
double maxF_
maximal deviation for contstraints
float jecFactor(const std::string &level, const std::string &flavor="none", const std::string &set="") const
Definition: Jet.cc:270
int i
Definition: DBlmapReader.cc:9
const std::vector< double > * jetEnergyResolutionEtaBinning_
double maxDeltaS_
maximal chi2 equivalent
TtFullHadKinFitter::Param param(unsigned int configParameter)
pat::Jet getCalHadj() const
void setFitHadj(const pat::Particle &aFitHadj)
CovarianceMatrix * covM_
get object resolutions and put them into a matrix
Param
supported parameterizations
Definition: TopKinFitter.h:22
void setupConstraints()
initialize constraints
TtHadEvtSolution addKinFitInfo(TtHadEvtSolution *asol)
add kin fit information to the old event solution (in for legacy reasons)
TAbsFitParticle * lightP_
Int_t fit()
Definition: TKinFitter.cc:309
TtFullHadKinFitter::Constraint constraint(unsigned int configParameter)
virtual void setIni4Vec(const TLorentzVector *pini)=0
void setFitHadb(const pat::Particle &aFitHadb)
pat::Jet getCalHadq() const
pat::Particle fittedLightPBar_
std::string param(const Param &param) const
convert Param to human readable form
Definition: TopKinFitter.cc:23
pat::Jet getCalHadp() const
void setFitHadp(const pat::Particle &aFitHadp)
virtual void setP4(const LorentzVector &p4)
set 4-momentum
pat::Particle fittedB_
output particles
TMatrixD setupMatrix(const pat::PATObject< T > &object, const TopKinFitter::Param param, const std::string &resolutionProvider="")
return covariance matrix for a PAT object
int fit(const std::vector< pat::Jet > &jets)
kinematic fit interface
void printSetup() const
print fitter setup
TtFullHadKinFitter * fitter
kinematic fit interface
double mW_
W mass value used for constraints.
std::vector< unsigned > constraints_
numbering of different possible kinematic constraints
math::Error< 5 >::type CovarianceMatrix
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:48
void setProbChi2(double c)
std::vector< double > jetEnergyResolutionScaleFactors_
scale factors for the jet energy resolution
std::vector< edm::ParameterSet > udscResolutions_
store the resolutions for the jets
double fitProb() const
return fit probability
Definition: TopKinFitter.h:36
std::vector< edm::ParameterSet > bResolutions_
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_
virtual double energy() const
energy
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&lt;unsigned int&gt;))
~TtFullHadKinFitter()
default destructor
const std::vector< edm::ParameterSet > * bResolutions_
void addConstraint(TAbsFitConstraint *constraint)
Definition: TKinFitter.cc:281
vector< PseudoJet > jets
tuple result
Definition: query.py:137
TAbsFitParticle * b_
input particles
pat::Jet getCalHadb() const
Int_t getStatus()
Definition: TKinFitter.h:41
virtual void setCovMatrix(const TMatrixD *theCovMatrix)
pat::Jet getCalHadk() const
double mW_
W mass value used for constraints.
Definition: TopKinFitter.h:54
void addMeasParticle(TAbsFitParticle *particle)
Definition: TKinFitter.cc:209
std::map< Constraint, TFitConstraintM * > massConstr_
supported constraints
std::vector< TtFullHadKinFitter::Constraint > constraints(const std::vector< unsigned int > &configParameters)
pat::Particle fittedLightP_
TAbsFitParticle * lightQBar_
pat::Particle fittedBBar_
double maxDeltaS_
maximal allowed chi2 (not normalized to degrees of freedom)
Definition: TopKinFitter.h:50
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
virtual double px() const
x coordinate of momentum vector
TAbsFitParticle * bBar_
~KinFit()
default destructor
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
Analysis-level particle class.
Definition: Particle.h:32
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:73
void setFitHadk(const pat::Particle &aFitHadk)
static const unsigned int nPartons
virtual double pz() const
z coordinate of momentum vector
const TLorentzVector * getCurr4Vec()
unsigned int jetParam_
numbering of different possible jet parametrizations
void setFitHadq(const pat::Particle &aFitHadq)
TtFullHadKinFitter()
default constructor
void setupJets()
initialize jet inputs
TKinFitter * fitter_
kinematic fitter
Definition: TopKinFitter.h:46
volatile std::atomic< bool > shutdown_flag false
tuple status
Definition: ntuplemaker.py:245
bool next_combination(BidIt n_begin, BidIt n_end, BidIt r_begin, BidIt r_end)
Definition: combination.h:22
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
pat::Jet getCalHadbbar() const
double mTop_
top mass value used for constraints
Definition: TopKinFitter.h:56
KinFit()
default constructor
virtual double py() const
y coordinate of momentum vector
TAbsFitParticle * lightPBar_
double maxF_
maximal allowed distance from constraints
Definition: TopKinFitter.h:52