File indexing completed on 2023-10-25 09:57:58
0001
0002
0003
0004
0005
0006
0007
0008
0009
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
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
0029
0030
0031
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
0055
0056
0057
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
0075
0076
0077 _ParList1.push_back(particle);
0078 }
0079
0080 void TFitConstraintM::addParticle2(TAbsFitParticle* particle) {
0081
0082
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
0098
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
0133
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
0159
0160 TFitConstraintM::~TFitConstraintM() {}
0161
0162
0163
0164
0165 TMatrixD* TFitConstraintM::getDerivative(TAbsFitParticle* particle) {
0166
0167
0168
0169 TMatrixD* DerivativeMatrix = new TMatrixD(1, 4);
0170 (*DerivativeMatrix) *= 0.;
0171
0172
0173
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
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
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
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
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
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
0278
0279 edm::LogVerbatim("KinFitter") << this->getInfoString();
0280 }