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