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