KinFitter/src/TFitConstraintM.cc

0001 // Classname: TFitConstraintM
0002 // Author: Jan E. Sundermann, Verena Klose (TU Dresden)
0003
0004 //________________________________________________________________
0005 //
0006 // TFitConstraintM::
0007 // --------------------
0008 //
0009 // Fit constraint: mass conservation ( m_i - m_j - MassConstraint == 0 )
0010 //
0011
0012 #include <iostream>
0013 #include <iomanip>
0014 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "TLorentzVector.h"
0017 #include "TClass.h"
0018
0019 //----------------
0020 // Constructor --
0021 //----------------
0022 TFitConstraintM::TFitConstraintM() : TAbsFitConstraint(), _ParList1(0), _ParList2(0), _TheMassConstraint(0) {}
0023
0024 TFitConstraintM::TFitConstraintM(std::vector<TAbsFitParticle*>* ParList1,
0025                                  std::vector<TAbsFitParticle*>* ParList2,
0026                                  Double_t Mass)
0027     : TAbsFitConstraint(), _ParList1(0), _ParList2(0) {
0028   // ParList1: Vector containing first list of constrained particles
0029   //           ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0030   // ParList2: Vector containing second list of constrained particles
0031   //           ( sum[ m_i ] - sum[ m_j ]  - MassConstraint == 0 )
0032
0033   if (ParList1) {
0034     _ParList1 = (*ParList1);
0035   }
0036   if (ParList2) {
0037     _ParList2 = (*ParList2);
0038   }
0039   if (Mass >= 0) {
0040     _TheMassConstraint = Mass;
0041   } else if (Mass < 0) {
0042     edm::LogWarning("NegativeMassConstr")
0043         << "Mass constraint in TFitConstraintM cannot be set to a negative value, will be set to 0.";
0044     _TheMassConstraint = 0.;
0045   }
0046 }
0047
0048 TFitConstraintM::TFitConstraintM(const TString& name,
0049                                  const TString& title,
0050                                  std::vector<TAbsFitParticle*>* ParList1,
0051                                  std::vector<TAbsFitParticle*>* ParList2,
0052                                  Double_t Mass)
0053     : TAbsFitConstraint(name, title), _ParList1(0), _ParList2(0) {
0054   // ParList1: Vector containing first list of constrained particles
0055   //           ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0056   // ParList2: Vector containing second list of constrained particles
0057   //           ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0058
0059   if (ParList1) {
0060     _ParList1 = (*ParList1);
0061   }
0062   if (ParList2) {
0063     _ParList2 = (*ParList2);
0064   }
0065   if (Mass >= 0) {
0066     _TheMassConstraint = Mass;
0067   } else if (Mass < 0) {
0068     edm::LogWarning("NegativeMassConstr")
0069         << "Mass constraint in TFitConstraintM cannot be set to a negative value, will be set to 0.";
0070     _TheMassConstraint = 0.;
0071   }
0072 }
0073 void TFitConstraintM::addParticle1(TAbsFitParticle* particle) {
0074   // Add one constrained particle to first list of particles
0075   // ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0076
0077   _ParList1.push_back(particle);
0078 }
0079
0080 void TFitConstraintM::addParticle2(TAbsFitParticle* particle) {
0081   // Add one constrained particle to second list of particles
0082   // ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0083
0084   _ParList2.push_back(particle);
0085 }
0086
0087 void TFitConstraintM::addParticles1(TAbsFitParticle* p1,
0088                                     TAbsFitParticle* p2,
0089                                     TAbsFitParticle* p3,
0090                                     TAbsFitParticle* p4,
0091                                     TAbsFitParticle* p5,
0092                                     TAbsFitParticle* p6,
0093                                     TAbsFitParticle* p7,
0094                                     TAbsFitParticle* p8,
0095                                     TAbsFitParticle* p9,
0096                                     TAbsFitParticle* p10) {
0097   // Add many constrained particle to first list of particles
0098   // ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0099
0100   if (p1)
0101     addParticle1(p1);
0102   if (p2)
0103     addParticle1(p2);
0104   if (p3)
0105     addParticle1(p3);
0106   if (p4)
0107     addParticle1(p4);
0108   if (p5)
0109     addParticle1(p5);
0110   if (p6)
0111     addParticle1(p6);
0112   if (p7)
0113     addParticle1(p7);
0114   if (p8)
0115     addParticle1(p8);
0116   if (p9)
0117     addParticle1(p9);
0118   if (p10)
0119     addParticle1(p10);
0120 }
0121
0122 void TFitConstraintM::addParticles2(TAbsFitParticle* p1,
0123                                     TAbsFitParticle* p2,
0124                                     TAbsFitParticle* p3,
0125                                     TAbsFitParticle* p4,
0126                                     TAbsFitParticle* p5,
0127                                     TAbsFitParticle* p6,
0128                                     TAbsFitParticle* p7,
0129                                     TAbsFitParticle* p8,
0130                                     TAbsFitParticle* p9,
0131                                     TAbsFitParticle* p10) {
0132   // Add many constrained particle to second list of particles
0133   // ( sum[ m_i ] - sum[ m_j ] - MassConstraint == 0 )
0134
0135   if (p1)
0136     addParticle2(p1);
0137   if (p2)
0138     addParticle2(p2);
0139   if (p3)
0140     addParticle2(p3);
0141   if (p4)
0142     addParticle2(p4);
0143   if (p5)
0144     addParticle2(p5);
0145   if (p6)
0146     addParticle2(p6);
0147   if (p7)
0148     addParticle2(p7);
0149   if (p8)
0150     addParticle2(p8);
0151   if (p9)
0152     addParticle2(p9);
0153   if (p10)
0154     addParticle2(p10);
0155 }
0156
0157 //--------------
0158 // Destructor --
0159 //--------------
0160 TFitConstraintM::~TFitConstraintM() {}
0161
0162 //--------------
0163 // Operations --
0164 //--------------
0165 TMatrixD* TFitConstraintM::getDerivative(TAbsFitParticle* particle) {
0166   // returns derivative df/dP with P=(p,E) and f the constraint (f=0).
0167   // The matrix contains one row (df/dp, df/dE).
0168
0169   TMatrixD* DerivativeMatrix = new TMatrixD(1, 4);
0170   (*DerivativeMatrix) *= 0.;
0171
0172   // Pf[4] is the 4-Mom (p,E) of the sum of particles on
0173   // the list particle is part of
0174
0175   Double_t Factor = 0.;
0176   TLorentzVector Pf(0., 0., 0., 0.);
0177
0178   if (OnList(&_ParList1, particle)) {
0179     UInt_t Npart = _ParList1.size();
0180     for (unsigned int i = 0; i < Npart; i++) {
0181       const TLorentzVector* FourVec = (_ParList1[i])->getCurr4Vec();
0182       Pf += (*FourVec);
0183     }
0184     if (Pf.M() == 0.) {
0185       edm::LogInfo("KinFitter") << "Division by zero in " << IsA()->GetName() << " (named " << GetName() << ", titled "
0186                                 << GetTitle() << ") will lead to Inf in derivative matrix for particle "
0187                                 << particle->GetName() << ".";
0188     }
0189     Factor = 1. / Pf.M();
0190   } else if (OnList(&_ParList2, particle)) {
0191     UInt_t Npart = _ParList2.size();
0192     for (unsigned int i = 0; i < Npart; i++) {
0193       const TLorentzVector* FourVec = (_ParList2[i])->getCurr4Vec();
0194       Pf += (*FourVec);
0195     }
0196     if (Pf.M() == 0.) {
0197       edm::LogInfo("KinFitter") << "Division by zero in " << IsA()->GetName() << " (named " << GetName() << ", titled "
0198                                 << GetTitle() << ") will lead to Inf in derivative matrix for particle "
0199                                 << particle->GetName() << ".";
0200     }
0201     Factor = -1. / Pf.M();
0202   } else {
0203     Factor = 0.;
0204   }
0205
0206   (*DerivativeMatrix)(0, 0) = -Pf[0];
0207   (*DerivativeMatrix)(0, 1) = -Pf[1];
0208   (*DerivativeMatrix)(0, 2) = -Pf[2];
0209   (*DerivativeMatrix)(0, 3) = +Pf[3];
0210   (*DerivativeMatrix) *= Factor;
0211
0212   return DerivativeMatrix;
0213 }
0214
0215 Double_t TFitConstraintM::getInitValue() {
0216   // Get initial value of constraint (before the fit)
0217
0218   Double_t InitValue(0);
0219   InitValue = CalcMass(&_ParList1, true) - CalcMass(&_ParList2, true) - _TheMassConstraint;
0220   return InitValue;
0221 }
0222
0223 Double_t TFitConstraintM::getCurrentValue() {
0224   // Get value of constraint after the fit
0225
0226   Double_t CurrentValue(0);
0227   CurrentValue = CalcMass(&_ParList1, false) - CalcMass(&_ParList2, false) - _TheMassConstraint;
0228   return CurrentValue;
0229 }
0230
0231 Bool_t TFitConstraintM::OnList(std::vector<TAbsFitParticle*>* List, TAbsFitParticle* particle) {
0232   // Checks whether list contains given particle
0233
0234   Bool_t ok(false);
0235   UInt_t Npart = List->size();
0236   for (unsigned int i = 0; i < Npart; i++) {
0237     ok = (particle == (*List)[i]);
0238     if (ok)
0239       break;
0240   }
0241   return ok;
0242 }
0243
0244 Double_t TFitConstraintM::CalcMass(std::vector<TAbsFitParticle*>* List, Bool_t IniVal) {
0245   // Calculates initial/current invariant mass of provided list of particles
0246
0247   TLorentzVector P(0., 0., 0., 0.);
0248   UInt_t Npart = List->size();
0249   for (unsigned int i = 0; i < Npart; i++) {
0250     const TLorentzVector* FourVec = nullptr;
0251     if (IniVal)
0252       FourVec = ((*List)[i])->getIni4Vec();
0253     else
0254       FourVec = ((*List)[i])->getCurr4Vec();
0255     P += (*FourVec);
0256   }
0257   return P.M();
0258 }
0259
0260 TString TFitConstraintM::getInfoString() {
0261   // Collect information to be used for printout
0262
0263   std::stringstream info;
0264   info << std::scientific << std::setprecision(6);
0265
0266   info << "__________________________" << std::endl << std::endl;
0267   info << "OBJ: " << IsA()->GetName() << "\t" << GetName() << "\t" << GetTitle() << std::endl;
0268
0269   info << "initial value: " << getInitValue() << std::endl;
0270   info << "current value: " << getCurrentValue() << std::endl;
0271   info << "mass: " << _TheMassConstraint << std::endl;
0272
0273   return info.str();
0274 }
0275
0276 void TFitConstraintM::print() {
0277   // Print constraint contents
0278
0279   edm::LogVerbatim("KinFitter") << this->getInfoString();
0280 }