00001 #include <iostream>
00002 #include "PhysicsTools/KinFitter/interface/TSLToyGen.h"
00003 #include "TMatrixD.h"
00004 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
00005 #include "PhysicsTools/KinFitter/interface/TFitConstraintEp.h"
00006 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
00007 #include "TH1.h"
00008 #include "TMath.h"
00009 #include "TRandom.h"
00010 #include "TString.h"
00011
00012
00013
00014 TSLToyGen::TSLToyGen( const TAbsFitParticle* bReco, const TAbsFitParticle* lepton,
00015 const TAbsFitParticle* X, const TAbsFitParticle* neutrino):
00016 _inimeasParticles(0),
00017 _iniunmeasParticles(0),
00018 _measParticles(0),
00019 _unmeasParticles(0),
00020 _Y4S(0., 0., 0.)
00021
00022 {
00023
00024
00025 _iniBreco = bReco->clone( bReco->GetName() + (TString) "INI" );
00026 _breco = bReco->clone( bReco->GetName() + (TString) "SMEAR" );
00027 _iniLepton = lepton->clone( lepton->GetName() + (TString) "INI" );
00028 _lepton = lepton->clone( lepton->GetName() + (TString) "SMEAR" );
00029 _iniX = X->clone( X->GetName() + (TString) "INI" );
00030 _X = X->clone( X->GetName() + (TString) "SMEAR" );
00031 _iniNeutrino = neutrino->clone( neutrino->GetName() +(TString) "INI" );
00032 _neutrino = neutrino->clone( neutrino->GetName() + (TString) "SMEAR" );
00033
00034 _printPartIni = false;
00035 _printConsIni = false;
00036 _printSmearedPartBefore = false;
00037 _printConsBefore = false;
00038 _printConsAfter = false;
00039 _printPartAfter = false;
00040 _withMassConstraint = false;
00041 _withMPDGCons = false;
00042 _doCheckConstraintsTruth = true;
00043 }
00044
00045 TSLToyGen::~TSLToyGen() {
00046
00047 delete _iniBreco;
00048 delete _iniLepton;
00049 delete _iniX;
00050 delete _iniNeutrino;
00051
00052 delete _breco;
00053 delete _lepton;
00054 delete _X;
00055 delete _neutrino;
00056
00057 }
00058
00059 Bool_t TSLToyGen::doToyExperiments( Int_t nbExperiments ) {
00060
00061
00062 TKinFitter fitter;
00063
00064 std::vector<TAbsFitParticle*> ParVec(0);
00065 ParVec.push_back(_breco);
00066 ParVec.push_back(_lepton);
00067 ParVec.push_back(_X);
00068 ParVec.push_back(_neutrino);
00069
00070 fitter.addMeasParticle(_breco);
00071 _inimeasParticles.push_back( _iniBreco );
00072 _measParticles.push_back( _breco );
00073 fitter.addMeasParticle(_lepton);
00074 _inimeasParticles.push_back( _iniLepton );
00075 _measParticles.push_back( _lepton );
00076 fitter.addMeasParticle(_X);
00077 _inimeasParticles.push_back( _iniX );
00078 _measParticles.push_back( _X );
00079 fitter.addUnmeasParticle(_neutrino);
00080 _iniunmeasParticles.push_back( _iniNeutrino );
00081 _iniunmeasParticles.push_back( _neutrino );
00082
00083
00084 _Y4S.SetXYZ(0., 0., 0.);
00085 for (unsigned int p = 0; p < _inimeasParticles.size(); p++) {
00086 _Y4S += _inimeasParticles[p]->getIni4Vec()->Vect();
00087 }
00088 _Y4S += _iniNeutrino->getIni4Vec()->Vect();
00089
00090 Double_t EY4S = TMath::Sqrt( _Y4S.Mag2() + 10.58*10.58 );
00091
00092
00093 TFitConstraintEp pXCons( "pX", "pX", &ParVec, TFitConstraintEp::pX, _Y4S.x() );
00094 TFitConstraintEp pYCons( "pY", "pY", &ParVec, TFitConstraintEp::pY, _Y4S.y() );
00095 TFitConstraintEp pZCons( "pZ", "pZ", &ParVec, TFitConstraintEp::pZ, _Y4S.z() );
00096 TFitConstraintEp ECons( "E", "E", &ParVec, TFitConstraintEp::E, EY4S );
00097 TFitConstraintM MCons( "MassConstraint", "Mass-Constraint", 0, 0 ,0);
00098 MCons.addParticle1( _breco );
00099 MCons.addParticles2( _lepton, _neutrino, _X );
00100 TFitConstraintM MPDGCons( "MPDGCons", "MPDGCons", 0, 0 , 5.279 );
00101 MPDGCons.addParticles1( _lepton, _neutrino, _X );
00102
00103
00104
00105
00106 fitter.addConstraint(&pXCons);
00107 fitter.addConstraint(&pYCons);
00108 fitter.addConstraint(&pZCons);
00109 fitter.addConstraint(&ECons);
00110 if (_withMassConstraint)
00111 fitter.addConstraint(&MCons);
00112 if(_withMPDGCons)
00113 fitter.addConstraint(&MPDGCons);
00114
00115
00116 fitter.setMaxNbIter( 50 );
00117 fitter.setMaxDeltaS( 5e-5 );
00118 fitter.setMaxF( 1e-4 );
00119 fitter.setVerbosity(0);
00120
00121 if( _printPartIni ) {
00122 std::cout << std::endl
00123 << "----------------------------------" << std::endl;
00124 std::cout << "--- PRINTING INITIAL PARTICLES ---" << std::endl;
00125 std::cout << "----------------------------------" << std::endl ;
00126 _iniBreco->print();
00127 _iniLepton->print();
00128 _iniX->print();
00129 _iniNeutrino->print();
00130 std::cout << std::endl << std::endl;
00131 }
00132
00133 if( _printConsIni ) {
00134 std::cout << std::endl
00135 << "-------------------------------------------------" << std::endl;
00136 std::cout << "INITIAL CONSTRAINTS " << std::endl ;
00137 std::cout << "-------------------------------------------------" << std::endl;
00138 std::cout << " M: " << MCons.getCurrentValue()
00139 << " MPDG: " << MPDGCons.getCurrentValue()
00140 << " px: " << pXCons.getCurrentValue()
00141 << " py: " << pYCons.getCurrentValue()
00142 << " pz: " << pZCons.getCurrentValue()
00143 << " E: " << ECons.getCurrentValue() << std::endl << std::endl;
00144 }
00145
00146
00147 if ( _doCheckConstraintsTruth ) {
00148 if (fitter.getF() > fitter.getMaxF()) {
00149
00150 return false;
00151 }
00152 }
00153
00154
00155 createHists();
00156
00157
00158 for (int i = 0; i < nbExperiments; i++) {
00159
00160 smearParticles();
00161
00162 if( _printSmearedPartBefore ) {
00163 std::cout << std::endl
00164 << "-------------------------------------------------------" << std::endl ;
00165 std::cout << "--- PRINTING SMEARED PARTICLES BEFORE FIT FOR experiment # " <<i+1 << std::endl;
00166 std::cout << "-------------------------------------------------------" << std::endl;
00167 _breco->print();
00168 _lepton->print();
00169 _X->print();
00170 _neutrino->print();
00171 }
00172
00173
00174 if( _printConsBefore ) {
00175 std::cout << std::endl
00176 << "-------------------------------------------------" << std::endl;
00177 std::cout << "INITIAL (SMEARED) CONSTRAINTS FOR experiment # "<< i+1 << std::endl ;
00178 std::cout << "-------------------------------------------------" << std::endl;
00179 std::cout << " M: " << MCons.getCurrentValue()
00180 << " px: " << pXCons.getCurrentValue()
00181 << " py: " << pYCons.getCurrentValue()
00182 << " pz: " << pZCons.getCurrentValue()
00183 << " E: " << ECons.getCurrentValue() << std::endl << std::endl;
00184 }
00185
00186 fitter.fit();
00187
00188 if( _printConsAfter) {
00189 std::cout << std::endl
00190 << "-------------------------------------------------" << std::endl;
00191 std::cout << " CONSTRAINTS AFTER FIT FOR experiment # "<< i+1 << std::endl ;
00192 std::cout << "-------------------------------------------------" << std::endl;
00193 std::cout << " M: " << MCons.getCurrentValue()
00194 << " MPDG: " << MPDGCons.getCurrentValue()
00195 << " px: " << pXCons.getCurrentValue()
00196 << " py: " << pYCons.getCurrentValue()
00197 << " pz: " << pZCons.getCurrentValue()
00198 << " E: " << ECons.getCurrentValue() << std::endl << std::endl;
00199 }
00200
00201 if( _printPartAfter ) {
00202 std::cout << std::endl
00203 << "--------------------------------------------------------" << std::endl ;
00204 std::cout << "--- PRINTING PARTICLES AFTER FIT FOR experiment # "<< i+1 << std::endl;
00205 std::cout << "--------------------------------------------------------" << std::endl;
00206 _breco->print();
00207 _lepton->print();
00208 _X->print();
00209 _neutrino->print();
00210 }
00211
00212 _histStatus->Fill( fitter.getStatus() );
00213 _histNIter->Fill( fitter.getNbIter() );
00214 if ( fitter.getStatus() == 0 ) {
00215 _histPChi2->Fill( TMath::Prob( fitter.getS(), fitter.getNDF() ) );
00216 _histChi2->Fill( fitter.getS());
00217 fillPull1();
00218 fillPull2();
00219 fillPar();
00220 fillM();
00221 }
00222
00223 if (i % 176 == 0) {
00224 std::cout << "\r";
00225 std::cout <<" ------ "<< (Double_t) i/nbExperiments*100. << " % PROCESSED ------";
00226 std::cout.flush();
00227 }
00228 }
00229
00230 return true;
00231
00232 }
00233
00234 void TSLToyGen::smearParticles() {
00235
00236
00237 for (unsigned int p = 0; p < _measParticles.size(); p++) {
00238
00239 TAbsFitParticle* particle = _measParticles[p];
00240 TAbsFitParticle* iniParticle = _inimeasParticles[p];
00241 TMatrixD parIni( *(iniParticle->getParIni()) );
00242 const TMatrixD* covM = iniParticle->getCovMatrix();
00243 for (int m = 0; m < iniParticle->getNPar(); m++) {
00244 parIni(m, 0) += gRandom->Gaus(0., TMath::Sqrt( (*covM)(m,m) ) );
00245 }
00246 TLorentzVector* ini4Vec = iniParticle->calc4Vec( &parIni );
00247 particle->setIni4Vec( ini4Vec );
00248 delete ini4Vec;
00249 TLorentzVector vectrue( *_inimeasParticles[p]->getIni4Vec() );
00250 TMatrixD* partrue = _measParticles[p]->transform( vectrue );
00251
00252 delete partrue;
00253 }
00254
00255
00256 TVector3 nuP3 = _Y4S;
00257 for (unsigned int p = 0; p < _measParticles.size(); p++) {
00258 nuP3 -= _measParticles[p]->getCurr4Vec()->Vect();
00259 }
00260 TLorentzVector ini4VecNeutrino ;
00261 ini4VecNeutrino.SetXYZM( nuP3.X(), nuP3.Y(), nuP3.Z(), 0. );
00262 _neutrino->setIni4Vec( &ini4VecNeutrino );
00263
00264 }
00265
00266 void TSLToyGen::fillPull1() {
00267
00268 Int_t histindex = 0;
00269 for (unsigned int p = 0; p < _measParticles.size(); p++) {
00270
00271
00272 TLorentzVector vectrue( *_inimeasParticles[p]->getIni4Vec() );
00273 TMatrixD* partrue = _measParticles[p]->transform( vectrue );
00274 const TMatrixD* parfit = _measParticles[p]->getParCurr();
00275
00276 TMatrixD parpull( *parfit );
00277 parpull -= (*partrue);
00278 const TMatrixD* covMatrixFit = _measParticles[p]->getCovMatrixFit();
00279 for (int i = 0; i < parpull.GetNrows(); i++) {
00280
00281 ((TH1D*) _histsDiff1[histindex])->Fill( parpull(i, 0) );
00282 parpull(i, 0) /= TMath::Sqrt( (*covMatrixFit)(i, i) );
00283 ((TH1D*) _histsPull1[histindex])->Fill( parpull(i, 0) );
00284 ((TH1D*) _histsError1[histindex])->Fill( TMath::Sqrt( (*covMatrixFit)(i, i) ) );
00285 histindex++;
00286
00287 }
00288 delete partrue;
00289
00290 }
00291
00292 }
00293
00294 void TSLToyGen::fillPull2() {
00295
00296 Int_t histindex = 0;
00297 for (unsigned int p = 0; p < _measParticles.size(); p++) {
00298
00299 const TMatrixD* pull = _measParticles[p]->getPull();
00300 const TMatrixD* VDeltaY = _measParticles[p]->getCovMatrixDeltaY();
00301 TMatrixD pardiff( *(_measParticles[p]->getParCurr()) );
00302 pardiff -= *(_measParticles[p]->getParIni());
00303 for (int i = 0; i < pull->GetNrows(); i++) {
00304 ( (TH1D*) _histsPull2[histindex])->Fill( (*pull)(i, 0) );
00305 ( (TH1D*) _histsError2[histindex])->Fill( TMath::Sqrt( (*VDeltaY)(i, i) ) );
00306 ( (TH1D*) _histsDiff2[histindex])->Fill( pardiff(i, 0) );
00307 histindex++;
00308 }
00309 }
00310
00311 }
00312
00313 void TSLToyGen::fillPar() {
00314
00315 Int_t histindex = 0;
00316 for (unsigned int p = 0; p < _measParticles.size(); p++) {
00317
00318 const TMatrixD* partrue = _inimeasParticles[p]->getParIni();
00319 const TMatrixD* parsmear = _measParticles[p]->getParIni();
00320 const TMatrixD* parfit = _measParticles[p]->getParCurr();
00321 for (int i = 0; i < partrue->GetNrows(); i++) {
00322 ( (TH1D*) _histsParTrue[histindex])->Fill( (*partrue)(i, 0) );
00323 ( (TH1D*) _histsParSmear[histindex])->Fill( (*parsmear)(i, 0) );
00324 ( (TH1D*) _histsParFit[histindex])->Fill( (*parfit)(i, 0) );
00325 histindex++;
00326 }
00327 }
00328
00329 }
00330
00331 void TSLToyGen::fillM() {
00332
00333 _histMBrecoTrue->Fill( _iniBreco->getIni4Vec()->M() );
00334 _histMBrecoSmear->Fill( _breco->getIni4Vec()->M() );
00335 _histMBrecoFit->Fill( _breco->getCurr4Vec()->M() );
00336
00337 _histMXTrue->Fill( _iniX->getIni4Vec()->M() );
00338 _histMXSmear->Fill( _X->getIni4Vec()->M() );
00339 _histMXFit->Fill( _X->getCurr4Vec()->M() );
00340
00341 TLorentzVector xlnutrue = *(_iniLepton->getIni4Vec());
00342 xlnutrue += *(_iniX->getIni4Vec());
00343 xlnutrue += *(_iniNeutrino->getIni4Vec());
00344 _histMXlnuTrue->Fill( xlnutrue.M() );
00345
00346 TLorentzVector xlnusmear = *(_lepton->getIni4Vec());
00347 xlnusmear += *(_X->getIni4Vec());
00348 xlnusmear += *(_neutrino->getIni4Vec());
00349 _histMXlnuSmear->Fill( xlnusmear.M() );
00350
00351 TLorentzVector xlnufit = *(_lepton->getCurr4Vec());
00352 xlnufit += *(_X->getCurr4Vec());
00353 xlnufit += *(_neutrino->getCurr4Vec());
00354 _histMXlnuFit->Fill( xlnufit.M() );
00355
00356 }
00357
00358 void TSLToyGen::createHists() {
00359
00360 _histStatus = new TH1D( "hStatus", "Status of the Fit", 16, -1, 15);
00361 _histNIter = new TH1D( "hNIter", "Number of iterations", 100, 0, 100);
00362 _histPChi2 = new TH1D( "hPChi2", "Chi2 probability", 100, 0., 1.);
00363 _histChi2 = new TH1D( "hChi2", "Chi2 ", 200, 0., 20.);
00364 _histMBrecoTrue = new TH1D( "histMBrecoTrue", "histMBrecoTrue", 2000, 4., 6.);
00365 _histMBrecoSmear = new TH1D( "histMBrecoSmear", "histMBrecoSmear", 2000, 4., 6.);
00366 _histMBrecoFit = new TH1D( "histMBrecoFit", "histMBrecoFit", 2000, 4., 6.);
00367 _histMXTrue = new TH1D( "histMXTrue", "histMXTrue", 600, 0., 6.);
00368 _histMXSmear = new TH1D( "histMXSmear", "histMXSmear", 600, 0., 6.);
00369 _histMXFit = new TH1D( "histMXFit", "histMXFit", 600, 0., 6.);
00370 _histMXlnuTrue = new TH1D( "histMXlnuTrue", "histMXlnuTrue", 3000, 4., 7.);
00371 _histMXlnuSmear = new TH1D( "histMXlnuSmear", "histMXlnuSmear", 500, 3., 8.);
00372 _histMXlnuFit = new TH1D( "histMXlnuFit", "histMXlnuFit", 3000, 4., 7.);
00373
00374 _histsParTrue.Clear();
00375 _histsParSmear.Clear();
00376 _histsParFit.Clear();
00377 _histsPull1.Clear();
00378 _histsError1.Clear();
00379 _histsDiff1.Clear();
00380 _histsPull2.Clear();
00381 _histsError2.Clear();
00382 _histsDiff2.Clear();
00383
00384 TObjArray histarrays;
00385 histarrays.Add( &_histsParTrue );
00386 histarrays.Add( &_histsParSmear );
00387 histarrays.Add( &_histsParFit );
00388 histarrays.Add( &_histsPull1 );
00389 histarrays.Add( &_histsError1 );
00390 histarrays.Add( &_histsDiff1 );
00391 histarrays.Add( &_histsPull2 );
00392 histarrays.Add( &_histsError2 );
00393 histarrays.Add( &_histsDiff2 );
00394
00395
00396 TString histnames[] = {"hParTrue", "hParSmear", "hParFit", "hPull1", "hError1", "hDiff1", "hPull2", "hError2", "hDiff2" };
00397
00398 TArrayD arrmin( histarrays.GetEntries() );
00399 TArrayD arrmax( histarrays.GetEntries() );
00400 arrmin[0] = 0.; arrmax[0] = 2.;
00401 arrmin[1] = 0.; arrmax[1] = 2.;
00402 arrmin[2] = 0.; arrmax[2] = 2.;
00403
00404 arrmin[3] = -3.; arrmax[3] = 3.;
00405 arrmin[4] = 0.; arrmax[4] = .2;
00406 arrmin[5] = -.5; arrmax[5] = .5;
00407 arrmin[6] = -3.; arrmax[6] = 3.;
00408 arrmin[7] = 0.; arrmax[7] = 0.2;
00409 arrmin[8] = -.5; arrmax[8] = .5;
00410
00411 for (unsigned int p = 0; p < _measParticles.size(); p++) {
00412
00413 const TAbsFitParticle* particle = _measParticles[p];
00414
00415
00416 for (int i = 0; i < particle->getNPar(); i++) {
00417 for (int h = 0; h < histarrays.GetEntries(); h++ ) {
00418 TString name = histnames[h] + (TString) particle->GetName();
00419 name += i;
00420 if ( h < 3) {
00421 const TMatrixD* parfit = _measParticles[p]->getParCurr();
00422 arrmin[h] = (*parfit)(i,0)*0.5;
00423 arrmax[h] = (*parfit)(i,0)*1.5;
00424 }
00425 TH1D* newhisto = new TH1D( name, name, 100, arrmin[h], arrmax[h]) ;
00426 ((TObjArray*) histarrays[h])->Add( newhisto );
00427
00428
00429 }
00430 }
00431 }
00432 }
1.7.3