KinFitter/test/testKinFitFWLite.C

0001 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0002 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
0003 #include "PhysicsTools/KinFitter/interface/TKinFitter.h"
0004
0005 #include "TLorentzVector.h"
0006
0007 #include <iostream>
0008
0009 ////////////////////////////////////////////////////////////////////////////////
0010 ////////////////////////////////////////////////////////////////////////////////
0011 // Modified version of the ExampleEtEtaPhi2CFit.C macro from CMS AN 2005/025.
0012 // Jet error parametrization from CMS AN 2005/005.
0013 //
0014 // To run this macro in a root session do:
0015 // root [0] gSystem->Load("libPhysicsToolsKinFitter.so");
0016 // root [1] .x PhysicsTools/KinFitter/test/testKinFitFWLite.C+
0017 ////////////////////////////////////////////////////////////////////////////////
0018 ////////////////////////////////////////////////////////////////////////////////
0019
0020 Double_t ErrEt(Float_t Et, Float_t Eta) {
0021   Double_t InvPerr2, a, b, c;
0022   if(fabs(Eta) < 1.4){
0023     a = 5.6;
0024     b = 1.25;
0025     c = 0.033;
0026   }
0027   else{
0028     a = 4.8;
0029     b = 0.89;
0030     c = 0.043;
0031   }
0032   InvPerr2 = (a * a) + (b * b) * Et + (c * c) * Et * Et;
0033   return InvPerr2;
0034 }
0035
0036 Double_t ErrEta(Float_t Et, Float_t Eta) {
0037   Double_t InvPerr2, a, b, c;
0038   if(fabs(Eta) < 1.4){
0039     a = 1.215;
0040     b = 0.037;
0041     c = 7.941 * 0.0001;
0042   }
0043   else{
0044     a = 1.773;
0045     b = 0.034;
0046     c = 3.56 * 0.0001;
0047   }
0048   InvPerr2 = a/(Et * Et) + b/Et + c;
0049   return InvPerr2;
0050 }
0051
0052 Double_t ErrPhi(Float_t Et, Float_t Eta) {
0053   Double_t InvPerr2, a, b, c;
0054   if(fabs(Eta) < 1.4){
0055     a = 6.65;
0056     b = 0.04;
0057     c = 8.49 * 0.00001;
0058   }
0059   else{
0060     a = 2.908;
0061     b = 0.021;
0062     c = 2.59 * 0.0001;
0063   }
0064   InvPerr2 = a/(Et * Et) + b/Et + c;
0065   return InvPerr2;
0066 }
0067
0068 void print(TKinFitter *fitter)
0069 {
0070   std::cout << "=============================================" << std ::endl;
0071   std::cout << "-> Number of measured Particles  : " << fitter->nbMeasParticles() << std::endl;
0072   std::cout << "-> Number of unmeasured particles: " << fitter->nbUnmeasParticles() << std::endl;
0073   std::cout << "-> Number of constraints         : " << fitter->nbConstraints() << std::endl;
0074   std::cout << "-> Number of degrees of freedom  : " << fitter->getNDF() << std::endl;
0075   std::cout << "-> Number of parameters A        : " << fitter->getNParA() << std::endl;
0076   std::cout << "-> Number of parameters B        : " << fitter->getNParB() << std::endl;
0077   std::cout << "-> Maximum number of iterations  : " << fitter->getMaxNumberIter() << std::endl;
0078   std::cout << "-> Maximum deltaS                : " << fitter->getMaxDeltaS() << std::endl;
0079   std::cout << "-> Maximum F                     : " << fitter->getMaxF() << std::endl;
0080   std::cout << "+++++++++++++++++++++++++++++++++++++++++++++" << std ::endl;
0081   std::cout << "-> Status                        : " << fitter->getStatus() << std::endl;
0082   std::cout << "-> Number of iterations          : " << fitter->getNbIter() << std::endl;
0083   std::cout << "-> S                             : " << fitter->getS() << std::endl;
0084   std::cout << "-> F                             : " << fitter->getF() << std::endl;
0085   std::cout << "=============================================" << std ::endl;
0086 }
0087
0088 void testKinFitFWLite()
0089 {
0090
0091   TLorentzVector v1(-77.92, 16.24, 117.64, 142.87);
0092   TLorentzVector v2( 15.41, 28.78,   6.06,  34.08);
0093   TLorentzVector v3(-99.57, 92.41,   2.54, 136.23);
0094
0095   TMatrixD m1(3,3);
0096   TMatrixD m2(3,3);
0097   TMatrixD m3(3,3);
0098   m1.Zero();
0099   m2.Zero();
0100   m3.Zero();
0101
0102   //In this example the covariant matrix depends on the transverse energy and eta of the jets
0103   m1(0,0) = ErrEt (v1.Et(), v1.Eta()); // et
0104   m1(1,1) = ErrEta(v1.Et(), v1.Eta()); // eta
0105   m1(2,2) = ErrPhi(v1.Et(), v1.Eta()); // phi
0106   m2(0,0) = ErrEt (v2.Et(), v2.Eta()); // et
0107   m2(1,1) = ErrEta(v2.Et(), v2.Eta()); // eta
0108   m2(2,2) = ErrPhi(v2.Et(), v2.Eta()); // phi
0109   m3(0,0) = ErrEt (v3.Et(), v3.Eta()); // et
0110   m3(1,1) = ErrEta(v3.Et(), v3.Eta()); // eta
0111   m3(2,2) = ErrPhi(v3.Et(), v3.Eta()); // phi
0112   TFitParticleEtEtaPhi *jet1 = new TFitParticleEtEtaPhi( "Jet1", "Jet1", &v1, &m1 );
0113   TFitParticleEtEtaPhi *jet2 = new TFitParticleEtEtaPhi( "Jet2", "Jet2", &v2, &m2 );
0114   TFitParticleEtEtaPhi *jet3 = new TFitParticleEtEtaPhi( "Jet3", "Jet3", &v3, &m3 );
0115
0116   //vec1 and vec2 must make a W boson
0117   TFitConstraintM *mCons1 = new TFitConstraintM( "WMassConstraint", "WMass-Constraint", 0, 0 , 80.4);
0118   mCons1->addParticles1( jet1, jet2 );
0119   //vec1 and vec2 and vec3 must make a top quark
0120   TFitConstraintM *mCons2 = new TFitConstraintM( "TopMassConstraint", "TopMass-Constraint", 0, 0 , 175.);
0121   mCons2->addParticles1( jet1, jet2, jet3 );
0122
0123   //Definition of the fitter
0124   //Add three measured particles(jets)
0125   //Add two constraints
0126   TKinFitter* fitter = new TKinFitter("fitter", "fitter");
0127   fitter->addMeasParticle( jet1 );
0128   fitter->addMeasParticle( jet2 );
0129   fitter->addMeasParticle( jet3 );
0130   fitter->addConstraint( mCons1 );
0131   fitter->addConstraint( mCons2 );
0132
0133   //Set convergence criteria
0134   fitter->setMaxNbIter( 30 );
0135   fitter->setMaxDeltaS( 1e-2 );
0136   fitter->setMaxF( 1e-1 );
0137   fitter->setVerbosity(1);
0138
0139   //Perform the fit
0140   std::cout << "Performing kinematic fit..." << std::endl;
0141   print(fitter);
0142   fitter->fit();
0143   std::cout << "Done." << std::endl;
0144   print(fitter);
0145
0146   delete jet1;
0147   delete jet2;
0148   delete jet3;
0149   delete mCons1;
0150   delete mCons2;
0151   delete fitter;
0152
0153 }