Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:23:37 

0001 #include <iostream>
0002 #include "PhysicsTools/KinFitter/interface/TSLToyGen.h"
0003 #include "TMatrixD.h"
0004 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0005 #include "PhysicsTools/KinFitter/interface/TFitConstraintEp.h"
0006 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
0007 #include "TH1.h"
0008 #include "TMath.h"
0009 #include "TRandom.h"
0010 #include "TString.h"
0011 
0012 TSLToyGen::TSLToyGen(const TAbsFitParticle* bReco,
0013                      const TAbsFitParticle* lepton,
0014                      const TAbsFitParticle* X,
0015                      const TAbsFitParticle* neutrino)
0016     : _inimeasParticles(0),
0017       _iniunmeasParticles(0),
0018       _measParticles(0),
0019       _unmeasParticles(0),
0020       _Y4S(0., 0., 0.)
0021 
0022 {
0023   // Clone input particles
0024   _iniBreco = bReco->clone(bReco->GetName() + (TString) "INI");
0025   _breco = bReco->clone(bReco->GetName() + (TString) "SMEAR");
0026   _iniLepton = lepton->clone(lepton->GetName() + (TString) "INI");
0027   _lepton = lepton->clone(lepton->GetName() + (TString) "SMEAR");
0028   _iniX = X->clone(X->GetName() + (TString) "INI");
0029   _X = X->clone(X->GetName() + (TString) "SMEAR");
0030   _iniNeutrino = neutrino->clone(neutrino->GetName() + (TString) "INI");
0031   _neutrino = neutrino->clone(neutrino->GetName() + (TString) "SMEAR");
0032 
0033   _printPartIni = false;
0034   _printConsIni = false;
0035   _printSmearedPartBefore = false;
0036   _printConsBefore = false;
0037   _printConsAfter = false;
0038   _printPartAfter = false;
0039   _withMassConstraint = false;
0040   _withMPDGCons = false;
0041   _doCheckConstraintsTruth = true;
0042 }
0043 
0044 TSLToyGen::~TSLToyGen() {
0045   delete _iniBreco;
0046   delete _iniLepton;
0047   delete _iniX;
0048   delete _iniNeutrino;
0049 
0050   delete _breco;
0051   delete _lepton;
0052   delete _X;
0053   delete _neutrino;
0054 }
0055 
0056 Bool_t TSLToyGen::doToyExperiments(Int_t nbExperiments) {
0057   // define fitter
0058   TKinFitter fitter;
0059 
0060   std::vector<TAbsFitParticle*> ParVec(0);
0061   ParVec.push_back(_breco);
0062   ParVec.push_back(_lepton);
0063   ParVec.push_back(_X);
0064   ParVec.push_back(_neutrino);
0065 
0066   fitter.addMeasParticle(_breco);
0067   _inimeasParticles.push_back(_iniBreco);
0068   _measParticles.push_back(_breco);
0069   fitter.addMeasParticle(_lepton);
0070   _inimeasParticles.push_back(_iniLepton);
0071   _measParticles.push_back(_lepton);
0072   fitter.addMeasParticle(_X);
0073   _inimeasParticles.push_back(_iniX);
0074   _measParticles.push_back(_X);
0075   fitter.addUnmeasParticle(_neutrino);
0076   _iniunmeasParticles.push_back(_iniNeutrino);
0077   _iniunmeasParticles.push_back(_neutrino);
0078 
0079   // Calculate Y4S
0080   _Y4S.SetXYZ(0., 0., 0.);
0081   for (unsigned int p = 0; p < _inimeasParticles.size(); p++) {
0082     _Y4S += _inimeasParticles[p]->getIni4Vec()->Vect();
0083   }
0084   _Y4S += _iniNeutrino->getIni4Vec()->Vect();
0085   //_Y4S.SetXYZ(-0.1212, -0.0033, 5.8784);
0086   Double_t EY4S = TMath::Sqrt(_Y4S.Mag2() + 10.58 * 10.58);
0087   //  std::cout << "_Y4S : " <<_Y4S.x() << " / " << _Y4S.y() << " / " << _Y4S.z() << " / " <<EY4S<< std::endl;
0088 
0089   TFitConstraintEp pXCons("pX", "pX", &ParVec, TFitConstraintEp::pX, _Y4S.x());
0090   TFitConstraintEp pYCons("pY", "pY", &ParVec, TFitConstraintEp::pY, _Y4S.y());
0091   TFitConstraintEp pZCons("pZ", "pZ", &ParVec, TFitConstraintEp::pZ, _Y4S.z());
0092   TFitConstraintEp ECons("E", "E", &ParVec, TFitConstraintEp::E, EY4S);
0093   TFitConstraintM MCons("MassConstraint", "Mass-Constraint", nullptr, nullptr, 0);
0094   MCons.addParticle1(_breco);
0095   MCons.addParticles2(_lepton, _neutrino, _X);
0096   TFitConstraintM MPDGCons("MPDGCons", "MPDGCons", nullptr, nullptr, 5.279);
0097   MPDGCons.addParticles1(_lepton, _neutrino, _X);
0098   //   TFitConstraintE EBCons( "EBXlnuCons", "EBXlnuCons", 0, 0 );
0099   //   EBCons.addParticle1( _breco );
0100   //   EBCons.addParticles2( _lepton, _neutrino, _X );
0101 
0102   fitter.addConstraint(&pXCons);
0103   fitter.addConstraint(&pYCons);
0104   fitter.addConstraint(&pZCons);
0105   fitter.addConstraint(&ECons);
0106   if (_withMassConstraint)
0107     fitter.addConstraint(&MCons);
0108   if (_withMPDGCons)
0109     fitter.addConstraint(&MPDGCons);
0110   // fitter.addConstraint(&EBCons);
0111 
0112   fitter.setMaxNbIter(50);
0113   fitter.setMaxDeltaS(5e-5);
0114   fitter.setMaxF(1e-4);
0115   fitter.setVerbosity(0);
0116 
0117   if (_printPartIni) {
0118     std::cout << std::endl << "----------------------------------" << std::endl;
0119     std::cout << "--- PRINTING INITIAL PARTICLES ---" << std::endl;
0120     std::cout << "----------------------------------" << std::endl;
0121     _iniBreco->print();
0122     _iniLepton->print();
0123     _iniX->print();
0124     _iniNeutrino->print();
0125     std::cout << std::endl << std::endl;
0126   }
0127 
0128   if (_printConsIni) {
0129     std::cout << std::endl << "-------------------------------------------------" << std::endl;
0130     std::cout << "INITIAL CONSTRAINTS " << std::endl;
0131     std::cout << "-------------------------------------------------" << std::endl;
0132     std::cout << "     M: " << MCons.getCurrentValue() << "  MPDG: " << MPDGCons.getCurrentValue()
0133               << "    px: " << pXCons.getCurrentValue() << "    py: " << pYCons.getCurrentValue()
0134               << "    pz: " << pZCons.getCurrentValue() << "     E: " << ECons.getCurrentValue() << std::endl
0135               << std::endl;
0136   }
0137 
0138   // Check initial constraints
0139   if (_doCheckConstraintsTruth) {
0140     if (fitter.getF() > fitter.getMaxF()) {
0141       //std::cout << "Initial constraints are not fulfilled." << std::endl;
0142       return false;
0143     }
0144   }
0145 
0146   // create histograms
0147   createHists();
0148 
0149   // perform pseudo experiments
0150   for (int i = 0; i < nbExperiments; i++) {
0151     smearParticles();
0152 
0153     if (_printSmearedPartBefore) {
0154       std::cout << std::endl << "-------------------------------------------------------" << std::endl;
0155       std::cout << "--- PRINTING SMEARED PARTICLES BEFORE FIT FOR experiment # " << i + 1 << std::endl;
0156       std::cout << "-------------------------------------------------------" << std::endl;
0157       _breco->print();
0158       _lepton->print();
0159       _X->print();
0160       _neutrino->print();
0161     }
0162 
0163     if (_printConsBefore) {
0164       std::cout << std::endl << "-------------------------------------------------" << std::endl;
0165       std::cout << "INITIAL (SMEARED) CONSTRAINTS FOR experiment # " << i + 1 << std::endl;
0166       std::cout << "-------------------------------------------------" << std::endl;
0167       std::cout << "     M: " << MCons.getCurrentValue() << "    px: " << pXCons.getCurrentValue()
0168                 << "    py: " << pYCons.getCurrentValue() << "    pz: " << pZCons.getCurrentValue()
0169                 << "     E: " << ECons.getCurrentValue() << std::endl
0170                 << std::endl;
0171     }
0172 
0173     fitter.fit();
0174 
0175     if (_printConsAfter) {
0176       std::cout << std::endl << "-------------------------------------------------" << std::endl;
0177       std::cout << " CONSTRAINTS AFTER FIT FOR experiment # " << i + 1 << std::endl;
0178       std::cout << "-------------------------------------------------" << std::endl;
0179       std::cout << "     M: " << MCons.getCurrentValue() << "  MPDG: " << MPDGCons.getCurrentValue()
0180                 << "    px: " << pXCons.getCurrentValue() << "    py: " << pYCons.getCurrentValue()
0181                 << "    pz: " << pZCons.getCurrentValue() << "     E: " << ECons.getCurrentValue() << std::endl
0182                 << std::endl;
0183     }
0184 
0185     if (_printPartAfter) {
0186       std::cout << std::endl << "--------------------------------------------------------" << std::endl;
0187       std::cout << "--- PRINTING PARTICLES AFTER FIT FOR experiment # " << i + 1 << std::endl;
0188       std::cout << "--------------------------------------------------------" << std::endl;
0189       _breco->print();
0190       _lepton->print();
0191       _X->print();
0192       _neutrino->print();
0193     }
0194 
0195     _histStatus->Fill(fitter.getStatus());
0196     _histNIter->Fill(fitter.getNbIter());
0197     if (fitter.getStatus() == 0) {
0198       _histPChi2->Fill(TMath::Prob(fitter.getS(), fitter.getNDF()));
0199       _histChi2->Fill(fitter.getS());
0200       fillPull1();
0201       fillPull2();
0202       fillPar();
0203       fillM();
0204     }
0205 
0206     if (i % 176 == 0) {
0207       std::cout << "\r";
0208       std::cout << " ------ " << (Double_t)i / nbExperiments * 100. << " % PROCESSED ------";
0209       std::cout.flush();
0210     }
0211   }
0212 
0213   return true;
0214 }
0215 
0216 void TSLToyGen::smearParticles() {
0217   // Smear measured particles
0218 
0219   for (unsigned int p = 0; p < _measParticles.size(); p++) {
0220     TAbsFitParticle* particle = _measParticles[p];
0221     TAbsFitParticle* iniParticle = _inimeasParticles[p];
0222     TMatrixD parIni(*(iniParticle->getParIni()));
0223     const TMatrixD* covM = iniParticle->getCovMatrix();
0224     for (int m = 0; m < iniParticle->getNPar(); m++) {
0225       parIni(m, 0) += gRandom->Gaus(0., TMath::Sqrt((*covM)(m, m)));
0226     }
0227     TLorentzVector* ini4Vec = iniParticle->calc4Vec(&parIni);
0228     particle->setIni4Vec(ini4Vec);
0229     delete ini4Vec;
0230     TLorentzVector vectrue(*_inimeasParticles[p]->getIni4Vec());
0231     TMatrixD* partrue = _measParticles[p]->transform(vectrue);
0232     //_measParticles[p]->setParIni(partrue);
0233     delete partrue;
0234   }
0235 
0236   // Calculate neutrino
0237   TVector3 nuP3 = _Y4S;
0238   for (unsigned int p = 0; p < _measParticles.size(); p++) {
0239     nuP3 -= _measParticles[p]->getCurr4Vec()->Vect();
0240   }
0241   TLorentzVector ini4VecNeutrino;
0242   ini4VecNeutrino.SetXYZM(nuP3.X(), nuP3.Y(), nuP3.Z(), 0.);
0243   _neutrino->setIni4Vec(&ini4VecNeutrino);
0244 }
0245 
0246 void TSLToyGen::fillPull1() {
0247   Int_t histindex = 0;
0248   for (unsigned int p = 0; p < _measParticles.size(); p++) {
0249     //const TAbsFitParticle* particle = _measParticles[p];
0250     TLorentzVector vectrue(*_inimeasParticles[p]->getIni4Vec());
0251     TMatrixD* partrue = _measParticles[p]->transform(vectrue);
0252     const TMatrixD* parfit = _measParticles[p]->getParCurr();
0253 
0254     TMatrixD parpull(*parfit);
0255     parpull -= (*partrue);
0256     const TMatrixD* covMatrixFit = _measParticles[p]->getCovMatrixFit();
0257     for (int i = 0; i < parpull.GetNrows(); i++) {
0258       ((TH1D*)_histsDiff1[histindex])->Fill(parpull(i, 0));
0259       parpull(i, 0) /= TMath::Sqrt((*covMatrixFit)(i, i));
0260       ((TH1D*)_histsPull1[histindex])->Fill(parpull(i, 0));
0261       ((TH1D*)_histsError1[histindex])->Fill(TMath::Sqrt((*covMatrixFit)(i, i)));
0262       histindex++;
0263     }
0264     delete partrue;
0265   }
0266 }
0267 
0268 void TSLToyGen::fillPull2() {
0269   Int_t histindex = 0;
0270   for (unsigned int p = 0; p < _measParticles.size(); p++) {
0271     const TMatrixD* pull = _measParticles[p]->getPull();
0272     const TMatrixD* VDeltaY = _measParticles[p]->getCovMatrixDeltaY();
0273     TMatrixD pardiff(*(_measParticles[p]->getParCurr()));
0274     pardiff -= *(_measParticles[p]->getParIni());
0275     for (int i = 0; i < pull->GetNrows(); i++) {
0276       ((TH1D*)_histsPull2[histindex])->Fill((*pull)(i, 0));
0277       ((TH1D*)_histsError2[histindex])->Fill(TMath::Sqrt((*VDeltaY)(i, i)));
0278       ((TH1D*)_histsDiff2[histindex])->Fill(pardiff(i, 0));
0279       histindex++;
0280     }
0281   }
0282 }
0283 
0284 void TSLToyGen::fillPar() {
0285   Int_t histindex = 0;
0286   for (unsigned int p = 0; p < _measParticles.size(); p++) {
0287     const TMatrixD* partrue = _inimeasParticles[p]->getParIni();
0288     const TMatrixD* parsmear = _measParticles[p]->getParIni();
0289     const TMatrixD* parfit = _measParticles[p]->getParCurr();
0290     for (int i = 0; i < partrue->GetNrows(); i++) {
0291       ((TH1D*)_histsParTrue[histindex])->Fill((*partrue)(i, 0));
0292       ((TH1D*)_histsParSmear[histindex])->Fill((*parsmear)(i, 0));
0293       ((TH1D*)_histsParFit[histindex])->Fill((*parfit)(i, 0));
0294       histindex++;
0295     }
0296   }
0297 }
0298 
0299 void TSLToyGen::fillM() {
0300   _histMBrecoTrue->Fill(_iniBreco->getIni4Vec()->M());
0301   _histMBrecoSmear->Fill(_breco->getIni4Vec()->M());
0302   _histMBrecoFit->Fill(_breco->getCurr4Vec()->M());
0303 
0304   _histMXTrue->Fill(_iniX->getIni4Vec()->M());
0305   _histMXSmear->Fill(_X->getIni4Vec()->M());
0306   _histMXFit->Fill(_X->getCurr4Vec()->M());
0307 
0308   TLorentzVector xlnutrue = *(_iniLepton->getIni4Vec());
0309   xlnutrue += *(_iniX->getIni4Vec());
0310   xlnutrue += *(_iniNeutrino->getIni4Vec());
0311   _histMXlnuTrue->Fill(xlnutrue.M());
0312 
0313   TLorentzVector xlnusmear = *(_lepton->getIni4Vec());
0314   xlnusmear += *(_X->getIni4Vec());
0315   xlnusmear += *(_neutrino->getIni4Vec());
0316   _histMXlnuSmear->Fill(xlnusmear.M());
0317 
0318   TLorentzVector xlnufit = *(_lepton->getCurr4Vec());
0319   xlnufit += *(_X->getCurr4Vec());
0320   xlnufit += *(_neutrino->getCurr4Vec());
0321   _histMXlnuFit->Fill(xlnufit.M());
0322 }
0323 
0324 void TSLToyGen::createHists() {
0325   _histStatus = new TH1D("hStatus", "Status of the Fit", 16, -1, 15);
0326   _histNIter = new TH1D("hNIter", "Number of iterations", 100, 0, 100);
0327   _histPChi2 = new TH1D("hPChi2", "Chi2 probability", 100, 0., 1.);
0328   _histChi2 = new TH1D("hChi2", "Chi2 ", 200, 0., 20.);
0329   _histMBrecoTrue = new TH1D("histMBrecoTrue", "histMBrecoTrue", 2000, 4., 6.);
0330   _histMBrecoSmear = new TH1D("histMBrecoSmear", "histMBrecoSmear", 2000, 4., 6.);
0331   _histMBrecoFit = new TH1D("histMBrecoFit", "histMBrecoFit", 2000, 4., 6.);
0332   _histMXTrue = new TH1D("histMXTrue", "histMXTrue", 600, 0., 6.);
0333   _histMXSmear = new TH1D("histMXSmear", "histMXSmear", 600, 0., 6.);
0334   _histMXFit = new TH1D("histMXFit", "histMXFit", 600, 0., 6.);
0335   _histMXlnuTrue = new TH1D("histMXlnuTrue", "histMXlnuTrue", 3000, 4., 7.);
0336   _histMXlnuSmear = new TH1D("histMXlnuSmear", "histMXlnuSmear", 500, 3., 8.);
0337   _histMXlnuFit = new TH1D("histMXlnuFit", "histMXlnuFit", 3000, 4., 7.);
0338 
0339   _histsParTrue.Clear();
0340   _histsParSmear.Clear();
0341   _histsParFit.Clear();
0342   _histsPull1.Clear();
0343   _histsError1.Clear();
0344   _histsDiff1.Clear();
0345   _histsPull2.Clear();
0346   _histsError2.Clear();
0347   _histsDiff2.Clear();
0348 
0349   TObjArray histarrays;
0350   histarrays.Add(&_histsParTrue);
0351   histarrays.Add(&_histsParSmear);
0352   histarrays.Add(&_histsParFit);
0353   histarrays.Add(&_histsPull1);
0354   histarrays.Add(&_histsError1);
0355   histarrays.Add(&_histsDiff1);
0356   histarrays.Add(&_histsPull2);
0357   histarrays.Add(&_histsError2);
0358   histarrays.Add(&_histsDiff2);
0359 
0360   TString histnames[] = {
0361       "hParTrue", "hParSmear", "hParFit", "hPull1", "hError1", "hDiff1", "hPull2", "hError2", "hDiff2"};
0362 
0363   TArrayD arrmin(histarrays.GetEntries());
0364   TArrayD arrmax(histarrays.GetEntries());
0365   arrmin[0] = 0.;
0366   arrmax[0] = 2.;  // Parameters
0367   arrmin[1] = 0.;
0368   arrmax[1] = 2.;
0369   arrmin[2] = 0.;
0370   arrmax[2] = 2.;
0371 
0372   arrmin[3] = -3.;
0373   arrmax[3] = 3.;  // Pull1
0374   arrmin[4] = 0.;
0375   arrmax[4] = .2;  //Error1
0376   arrmin[5] = -.5;
0377   arrmax[5] = .5;  //Diff1
0378   arrmin[6] = -3.;
0379   arrmax[6] = 3.;  // Pull2
0380   arrmin[7] = 0.;
0381   arrmax[7] = 0.2;  //Error2
0382   arrmin[8] = -.5;
0383   arrmax[8] = .5;  //Diff2
0384 
0385   for (unsigned int p = 0; p < _measParticles.size(); p++) {
0386     const TAbsFitParticle* particle = _measParticles[p];
0387     //    const TMatrixD* covMatrix = particle->getCovMatrix();
0388 
0389     for (int i = 0; i < particle->getNPar(); i++) {
0390       for (int h = 0; h < histarrays.GetEntries(); h++) {
0391         TString name = histnames[h] + (TString)particle->GetName();
0392         name += i;
0393         if (h < 3) {
0394           const TMatrixD* parfit = _measParticles[p]->getParCurr();
0395           arrmin[h] = (*parfit)(i, 0) * 0.5;
0396           arrmax[h] = (*parfit)(i, 0) * 1.5;
0397         }
0398         TH1D* newhisto = new TH1D(name, name, 100, arrmin[h], arrmax[h]);
0399         ((TObjArray*)histarrays[h])->Add(newhisto);
0400 
0401         //  newhisto->SetCanExtend(TH1::kAllAxes);
0402       }
0403     }
0404   }
0405 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!