test/unitTests/main.cc

0001 #include <iostream>
0002 #include <execinfo.h>
0003 #include <signal.h>
0004 #include <stdlib.h>
0005
0006 #include <AbsFinitePlane.h>
0007 #include <AbsFitterInfo.h>
0008 #include <AbsMeasurement.h>
0009 #include <AbsTrackRep.h>
0010 #include <ConstField.h>
0011 #include <DetPlane.h>
0012 #include <Exception.h>
0013 #include <FieldManager.h>
0014 #include <KalmanFittedStateOnPlane.h>
0015 #include <KalmanFitterInfo.h>
0016 #include <MeasuredStateOnPlane.h>
0017 #include <MeasurementOnPlane.h>
0018 #include <PlanarMeasurement.h>
0019 #include <ProlateSpacepointMeasurement.h>
0020 #include <RectangularFinitePlane.h>
0021 #include <ReferenceStateOnPlane.h>
0022 #include <SharedPlanePtr.h>
0023 #include <SpacepointMeasurement.h>
0024 #include <StateOnPlane.h>
0025 #include <Tools.h>
0026 #include <TrackCand.h>
0027 #include <TrackCandHit.h>
0028 #include <Track.h>
0029 #include <TrackPoint.h>
0030 #include <WireMeasurement.h>
0031 #include <WirePointMeasurement.h>
0032
0033 #include <MaterialEffects.h>
0034 #include <RKTools.h>
0035 #include <RKTrackRep.h>
0036 #include <StepLimits.h>
0037 #include <TGeoMaterialInterface.h>
0038
0039 #include <TApplication.h>
0040 #include <TCanvas.h>
0041 #include <TDatabasePDG.h>
0042 #include <TEveManager.h>
0043 #include <TGeoManager.h>
0044 #include <TH1D.h>
0045 #include <TH2D.h>
0046 #include <TRandom.h>
0047 #include <TStyle.h>
0048 #include <TVector3.h>
0049
0050 #include <vector>
0051 #include <map>
0052 #include <stdio.h>
0053 #include <stdlib.h>
0054
0055 #include <TROOT.h>
0056 #include <TFile.h>
0057 #include <TTree.h>
0058 #include "TDatabasePDG.h"
0059 #include <TMath.h>
0060 #include <TString.h>
0061
0062 //#include <callgrind/callgrind.h>
0063
0064
0065 enum e_testStatus {
0066   kPassed,
0067   kFailed,
0068   kException
0069 };
0070
0071 constexpr bool verbose = false;
0072
0073 void handler(int sig) {
0074   void *array[10];
0075   size_t size;
0076
0077   // get void*'s for all entries on the stack
0078   size = backtrace(array, 10);
0079
0080   // print out all the frames to stderr
0081   fprintf(stderr, "Error: signal %d:\n", sig);
0082   backtrace_symbols_fd(array, size, 2);
0083   exit(1);
0084 }
0085
0086 int randomPdg() {
0087   int pdg;
0088
0089   switch(int(gRandom->Uniform(8))) {
0090   case 1:
0091     pdg = -11; break;
0092   case 2:
0093     pdg = 11; break;
0094   case 3:
0095     pdg = 13; break;
0096   case 4:
0097     pdg = -13; break;
0098   case 5:
0099     pdg = 211; break;
0100   case 6:
0101     pdg = -211; break;
0102   case 7:
0103     pdg = 2212; break;
0104   default:
0105     pdg = 211;
0106   }
0107
0108   return pdg;
0109 }
0110
0111
0112 int randomSign() {
0113   if (gRandom->Uniform(1) > 0.5)
0114     return 1;
0115   return -1;
0116 }
0117
0118
0119 e_testStatus compareMatrices(const TMatrixTBase<double>& A, const TMatrixTBase<double>& B, double maxRelErr) {
0120   e_testStatus retVal = kPassed;
0121
0122   // search max abs value
0123   double max(0);
0124   for (int i=0; i<A.GetNrows(); ++i) {
0125     for (int j=0; j<A.GetNcols(); ++j) {
0126       if (fabs(A(i,j)) > max)
0127         max = fabs(A(i,j));
0128     }
0129   }
0130
0131   double maxAbsErr = maxRelErr*max;
0132
0133   for (int i=0; i<A.GetNrows(); ++i) {
0134     for (int j=0; j<A.GetNcols(); ++j) {
0135       double absErr = A(i,j) - B(i,j);
0136       if ( fabs(absErr) > maxAbsErr ) {
0137         double relErr = A(i,j)/B(i,j) - 1;
0138         if ( fabs(relErr) > maxRelErr ) {
0139           if (verbose)  {
0140             std::cout << "compareMatrices: A("<<i<<","<<j<<") = " << A(i,j) << "  B("<<i<<","<<j<<") = " << B(i,j) << "     absErr = " << absErr << "    relErr = " << relErr << "\n";
0141           }
0142           retVal = kFailed;
0143         }
0144       }
0145     }
0146   }
0147   return retVal;
0148 }
0149
0150 e_testStatus isCovMatrix(TMatrixTBase<double>& cov) {
0151
0152   if (!(cov.IsSymmetric())) {
0153     if (verbose) {
0154       std::cout << "isCovMatrix: not symmetric\n";
0155     }
0156     return kFailed;
0157   }
0158
0159   for (int i=0; i<cov.GetNrows(); ++i) {
0160     for (int j=0; j<cov.GetNcols(); ++j) {
0161        if (std::isnan(cov(i,j))) {
0162          if (verbose) {
0163            std::cout << "isCovMatrix: element isnan\n";
0164          }
0165          return kFailed;
0166        }
0167        if (i==j && cov(i,j) < 0) {
0168          if (verbose) {
0169            std::cout << "isCovMatrix: negative diagonal element\n";
0170          }
0171          return kFailed;
0172        }
0173     }
0174   }
0175
0176   return kPassed;
0177 }
0178
0179
0180
0181 e_testStatus checkSetGetPosMom(bool writeHisto = false) {
0182
0183   if (writeHisto)
0184     return kPassed;
0185
0186   double epsilonLen = 1.E-10;
0187   double epsilonMom = 1.E-10;
0188
0189   int pdg = randomPdg();
0190   genfit::AbsTrackRep* rep;
0191   rep = new genfit::RKTrackRep(pdg);
0192
0193   //TVector3 pos(0,0,0);
0194   TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
0195   TVector3 mom(0,0.5,gRandom->Gaus(0,0.3));
0196   mom.SetMag(0.5);
0197   mom *= randomSign();
0198
0199
0200   genfit::StateOnPlane state(rep);
0201   rep->setPosMom(state, pos, mom);
0202
0203   // check if we can set another position in the same plane
0204   if (randomSign() == 1) {
0205     genfit::SharedPlanePtr plane = state.getPlane();
0206     const TVector3& u = plane->getU();
0207     const TVector3& v = plane->getV();
0208
0209     // random position on plane
0210     pos += gRandom->Gaus() * u;
0211     pos += gRandom->Gaus() * v;
0212
0213     // new random momentum
0214     mom.SetXYZ(0,0.5,gRandom->Gaus(0,0.3));
0215     mom.SetMag(0.5);
0216     mom *= randomSign();
0217
0218     rep->setPosMom(state, pos, mom);
0219
0220     // check if plane has changed
0221     if (state.getPlane() != plane) {
0222         if (verbose) {
0223           std::cout << "plane has changed unexpectedly! \n";
0224         }
0225       delete rep;
0226       return kFailed;
0227     }
0228   }
0229
0230
0231   // compare
0232   if ((pos - rep->getPos(state)).Mag() > epsilonLen ||
0233       (mom - rep->getMom(state)).Mag() > epsilonMom) {
0234
0235     if (verbose) {
0236       state.Print();
0237       std::cout << "pos difference = " << (pos - rep->getPos(state)).Mag() << "\n";
0238       std::cout << "mom difference = " << (mom - rep->getMom(state)).Mag() << "\n";
0239
0240       std::cout << std::endl;
0241     }
0242     delete rep;
0243     return kFailed;
0244   }
0245
0246   delete rep;
0247   return kPassed;
0248
0249 }
0250
0251
0252 e_testStatus compareForthBackExtrapolation(bool writeHisto = false) {
0253   static std::map<int, std::vector<TH2D*> > histoMap;
0254
0255   static const int nPDGs(5);
0256   int pdgs[nPDGs]={0, 211,13,11,2212};
0257   static bool fill(true);
0258   if (fill) {
0259     for (int i = 0; i<nPDGs; ++i) {
0260       int pdg = pdgs[i];
0261
0262       std::string Result;//string which will contain the result
0263       std::stringstream convert; // stringstream used for the conversion
0264       convert << pdg;//add the value of Number to the characters in the stream
0265       Result = convert.str();//set Result to the content of the stream
0266
0267       histoMap[pdg].push_back(new TH2D((std::string("deviationRel_")+Result).c_str(), "log(betaGamma) vs relative deviation", 100000, -1.e-2, 1.e-2, 50, -4, 8));
0268       histoMap[pdg].push_back(new TH2D((std::string("deviationAbs_")+Result).c_str(), "log(betaGamma) vs absolute deviation; deviation (keV)", 100000, -90.0, 10.0, 50, -4, 8));
0269       histoMap[pdg].push_back(new TH2D((std::string("ExtrapLen_")+Result).c_str(), "delta ExtrapLen vs relative deviation", 50000, -5.e-2, 5.e-2, 400, -0.1, 0.1));
0270     }
0271     fill = false;
0272   }
0273
0274
0275   if (writeHisto) {
0276     TFile outfile("deviation.root", "recreate");
0277     outfile.cd();
0278
0279     for (int i = 0; i<nPDGs; ++i) {
0280       int pdg = pdgs[i];
0281       histoMap[pdg][0]->Write();
0282       histoMap[pdg][1]->Write();
0283       histoMap[pdg][2]->Write();
0284     }
0285
0286     outfile.Close();
0287
0288     return kPassed;
0289   }
0290
0291
0292   double epsilonLen = 5.E-5; // 0.5 mu
0293   double epsilonMom = 1.E-6; // 1 keV
0294
0295   int pdg = randomPdg();
0296   //pdg = 211;
0297   genfit::AbsTrackRep* rep;
0298   rep = new genfit::RKTrackRep(pdg);
0299
0300   //rep->setDebugLvl(1);
0301   //genfit::MaterialEffects::getInstance()->setDebugLvl(1);
0302
0303
0304
0305   TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
0306   double mass = part->Mass(); // GeV
0307
0308   //TVector3 pos(0,0,0);
0309   TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
0310   TVector3 mom(0,0.5,gRandom->Gaus(0,0.3));
0311   mom.SetMag( exp(gRandom->Uniform(-4, 8)) * mass );
0312   mom *= randomSign();
0313
0314
0315   double betaGamma = log(mom.Mag()/mass);
0316
0317   //mom.Print();
0318
0319   genfit::StateOnPlane state(rep);
0320   rep->setPosMom(state, pos, mom);
0321
0322   genfit::SharedPlanePtr origPlane = state.getPlane();
0323   genfit::SharedPlanePtr plane(new genfit::DetPlane(TVector3(0,randomSign()*10,0), TVector3(0,randomSign()*1,0)));
0324
0325   genfit::StateOnPlane origState(state);
0326
0327   TVector3 mom2;
0328   double momLoss1, momLoss2;
0329
0330   // forth
0331   double extrapLen(0);
0332   try {
0333     extrapLen = rep->extrapolateToPlane(state, plane);
0334
0335     mom2 = state.getMom();
0336     momLoss1 = mom.Mag()-mom2.Mag();
0337
0338     //mom2.Print();
0339     //std::cout << "MomLoss = " << momLoss1 << "\n";
0340   }
0341   catch (genfit::Exception& e) {
0342     if (verbose) {
0343       std::cerr << "Exception in forth Extrapolation. PDG = " << pdg << "; mom: \n";
0344       mom.Print();
0345
0346       std::cerr << e.what();
0347     }
0348     delete rep;
0349     return kException;
0350   }
0351
0352
0353
0354   // back
0355   double backExtrapLen(0);
0356   try {
0357     backExtrapLen = rep->extrapolateToPlane(state, origPlane);
0358
0359     momLoss2 = mom2.Mag()-state.getMom().Mag();
0360
0361     //state.getMom().Print();
0362     //std::cout << "MomLoss = " << momLoss2 << "\n";
0363
0364     double deviation = 1. + momLoss1/momLoss2;
0365     histoMap[abs(pdg)][0]->Fill(deviation, betaGamma);
0366     histoMap[abs(pdg)][1]->Fill((mom.Mag() - state.getMom().Mag())*1e6, betaGamma);
0367     histoMap[abs(pdg)][2]->Fill(deviation, extrapLen+backExtrapLen);
0368
0369     histoMap[0][0]->Fill(deviation, betaGamma);
0370     histoMap[0][1]->Fill((mom.Mag() - state.getMom().Mag())*1e6, betaGamma);
0371     histoMap[0][2]->Fill(deviation, extrapLen+backExtrapLen);
0372
0373     //std::cout << "deviation = " << deviation << "\n";
0374   }
0375   catch (genfit::Exception& e) {
0376     if (verbose) {
0377       std::cerr << "Exception in back Extrapolation. PDG = " << pdg << "; mom:  \n";
0378       mom.Print();
0379       std::cerr << "mom2:  \n";
0380       mom2.Print();
0381     }
0382     std::cerr << e.what();
0383
0384     delete rep;
0385     return kException;
0386   }
0387
0388   // compare
0389   if ((rep->getPos(origState) - rep->getPos(state)).Mag() > epsilonLen ||
0390       (rep->getMom(origState) - rep->getMom(state)).Mag() > epsilonMom ||
0391       fabs(extrapLen + backExtrapLen) > epsilonLen) {
0392
0393     if (verbose) {
0394         origState.Print();
0395         state.Print();
0396
0397         std::cout << "pos difference = " << (rep->getPos(origState) - rep->getPos(state)).Mag() << "\n";
0398         std::cout << "mom difference = " << (rep->getMom(origState) - rep->getMom(state)).Mag() << "\n";
0399         std::cout << "len difference = " << extrapLen + backExtrapLen << "\n";
0400
0401         std::cout << std::endl;
0402     }
0403     delete rep;
0404     return kFailed;
0405   }
0406
0407   delete rep;
0408   return kPassed;
0409
0410 }
0411
0412
0413 e_testStatus compareForthBackJacNoise(bool writeHisto = false) {
0414
0415   if (writeHisto)
0416     return kPassed;
0417
0418   e_testStatus retVal(kPassed);
0419
0420   bool fx( randomSign() > 0);
0421   genfit::MaterialEffects::getInstance()->setNoEffects(!fx);
0422
0423   double deltaJac = 0.005; // relative
0424   double deltaNoise = 0.00001;
0425   double deltaState = 3.E-6; // absolute
0426
0427   if (fx) {
0428     deltaJac = 0.1; // relative
0429     deltaNoise = 0.1;
0430     deltaState = 5.E-4; // absolute
0431   }
0432
0433   int pdg = randomPdg();
0434   genfit::AbsTrackRep* rep;
0435   rep = new genfit::RKTrackRep(pdg);
0436
0437   //TVector3 pos(0,0,0);
0438   //TVector3 mom(0,1,2);
0439   TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
0440   TVector3 mom(0, 0.5, gRandom->Gaus(0, 1));
0441   mom *= randomSign();
0442   mom.SetMag(gRandom->Uniform(2)+0.3);
0443   //mom.SetMag(3);
0444
0445   TMatrixD jac_f, jac_fi, jac_b, jac_bi;
0446   TMatrixDSym noise_f, noise_fi, noise_b, noise_bi;
0447   TVectorD c_f, c_fi, c_b, c_bi;
0448   TVectorD state_man, stateOrig_man;
0449
0450   // original state and plane
0451   genfit::MeasuredStateOnPlane state(rep);
0452   rep->setPosMom(state, pos, mom);
0453
0454   static const double smear = 0.2;
0455   TVector3 normal(mom);
0456   normal.SetMag(1);
0457   normal.SetXYZ(gRandom->Gaus(normal.X(), smear),
0458       gRandom->Gaus(normal.Y(), smear),
0459       gRandom->Gaus(normal.Z(), smear));
0460   genfit::DetPlane* origPlanePtr = new genfit::DetPlane (pos, normal);
0461   //genfit::DetPlane* origPlanePtr = new genfit::DetPlane (pos, TVector3(1,0,0), TVector3(0,0,1));
0462   double rotAngleOrig = gRandom->Uniform(2.*TMath::Pi());
0463   origPlanePtr->rotate(rotAngleOrig);
0464   genfit::SharedPlanePtr origPlane(origPlanePtr);
0465   //genfit::SharedPlanePtr origPlane = state.getPlane();
0466   rep->extrapolateToPlane(state, origPlane);
0467
0468   const genfit::StateOnPlane origState(state);
0469
0470
0471   // dest plane
0472   normal = mom;
0473   normal.SetMag(1);
0474   normal.SetXYZ(gRandom->Gaus(normal.X(), smear),
0475       gRandom->Gaus(normal.Y(), smear),
0476       gRandom->Gaus(normal.Z(), smear));
0477   TVector3 dest(mom);
0478   dest.SetMag(10);
0479   genfit::DetPlane* planePtr = new genfit::DetPlane (dest, normal);
0480   //genfit::DetPlane* planePtr = new genfit::DetPlane (dest, TVector3(1,0,0), TVector3(0,0,1));
0481   double rotAngle = gRandom->Uniform(2.*TMath::Pi());
0482   planePtr->rotate(rotAngle);
0483   genfit::SharedPlanePtr plane(planePtr);
0484
0485  /* genfit::DetPlane* planePtr = new genfit::DetPlane (*origPlane);
0486   planePtr->setO(TVector3(0,randomSign()*10,0));
0487   //planePtr->rotate(rotAngle);
0488   genfit::SharedPlanePtr plane(planePtr);
0489 */
0490
0491   // numerical calculation
0492   TMatrixD jac_f_num;
0493   rep->calcJacobianNumerically(origState, plane, jac_f_num);
0494
0495
0496   // forth
0497   genfit::StateOnPlane extrapolatedState;
0498   try {
0499     //std::cout << "DO FORTH EXTRAPOLATION \n";
0500     rep->extrapolateToPlane(state, plane);
0501     //std::cout << "GET INFO FOR FORTH EXTRAPOLATION \n";
0502     extrapolatedState = state;
0503     rep->getForwardJacobianAndNoise(jac_f, noise_f, c_f);
0504     rep->getBackwardJacobianAndNoise(jac_fi, noise_fi, c_fi);
0505   }
0506   catch (genfit::Exception& e) {
0507     std::cerr << e.what();
0508
0509     delete rep;
0510     genfit::MaterialEffects::getInstance()->setNoEffects(false);
0511     return kException;
0512   }
0513
0514   // back
0515   try {
0516     //std::cout << "DO BACK EXTRAPOLATION \n";
0517     rep->extrapolateToPlane(state, origPlane);
0518     //std::cout << "GET INFO FOR BACK EXTRAPOLATION \n";
0519     rep->getForwardJacobianAndNoise(jac_b, noise_b, c_b);
0520     rep->getBackwardJacobianAndNoise(jac_bi, noise_bi, c_bi);
0521   }
0522   catch (genfit::Exception& e) {
0523     std::cerr << e.what();
0524
0525     delete rep;
0526     genfit::MaterialEffects::getInstance()->setNoEffects(false);
0527     return kException;
0528   }
0529
0530
0531   // compare
0532   if (!((origState.getState() - state.getState()).Abs()  < deltaState) ) {
0533     if (verbose) {
0534       std::cout << "(origState.getState() - state.getState()) ";
0535       (origState.getState() - state.getState()).Print();
0536     }
0537
0538     retVal = kFailed;
0539   }
0540
0541   // check c
0542   if (!((jac_f * origState.getState() + c_f  -  extrapolatedState.getState()).Abs() < deltaState) ||
0543       !((jac_bi * origState.getState() + c_bi  -  extrapolatedState.getState()).Abs() < deltaState) ||
0544       !((jac_b * extrapolatedState.getState() + c_b  -  origState.getState()).Abs() < deltaState) ||
0545       !((jac_fi * extrapolatedState.getState() + c_fi  -  origState.getState()).Abs() < deltaState)   ) {
0546
0547     if (verbose) {
0548       std::cout << "(jac_f * origState.getState() + c_f  -  extrapolatedState.getState()) ";
0549       (jac_f * origState.getState() + c_f - extrapolatedState.getState()).Print();
0550       std::cout << "(jac_bi * origState.getState() + c_bi  -  extrapolatedState.getState()) ";
0551       (jac_bi * origState.getState() + c_bi - extrapolatedState.getState()).Print();
0552       std::cout << "(jac_b * extrapolatedState.getState() + c_b  -  origState.getState()) ";
0553       (jac_b * extrapolatedState.getState() + c_b - origState.getState()).Print();
0554       std::cout << "(jac_fi * extrapolatedState.getState() + c_fi  -  origState.getState()) ";
0555       (jac_fi * extrapolatedState.getState() + c_fi - origState.getState()).Print();
0556     }
0557     retVal = kFailed;
0558   }
0559
0560   if (isCovMatrix(state.getCov()) == kFailed) {
0561     retVal = kFailed;
0562   }
0563
0564
0565   // compare
0566   if (compareMatrices(jac_f, jac_bi, deltaJac) == kFailed) {
0567     if (verbose) {
0568       std::cout << "jac_f = ";
0569       jac_f.Print();
0570       std::cout << "jac_bi = ";
0571       jac_bi.Print();
0572       std::cout << std::endl;
0573     }
0574     retVal = kFailed;
0575   }
0576
0577   // compare
0578   if (compareMatrices(jac_b, jac_fi, deltaJac) == kFailed) {
0579     if (verbose) {
0580       std::cout << "jac_b = ";
0581       jac_b.Print();
0582       std::cout << "jac_fi = ";
0583       jac_fi.Print();
0584       std::cout << std::endl;
0585     }
0586     retVal = kFailed;
0587   }
0588
0589   // compare
0590   if (compareMatrices(noise_f, noise_bi, deltaNoise) == kFailed) {
0591     if (verbose) {
0592       std::cout << "noise_f = ";
0593       noise_f.Print();
0594       std::cout << "noise_bi = ";
0595       noise_bi.Print();
0596       std::cout << std::endl;
0597     }
0598     retVal = kFailed;
0599   }
0600
0601   // compare
0602   if (compareMatrices(noise_b, noise_fi, deltaNoise) == kFailed) {
0603     if (verbose) {
0604       std::cout << "noise_b = ";
0605       noise_b.Print();
0606       std::cout << "noise_fi = ";
0607       noise_fi.Print();
0608       std::cout << std::endl;
0609     }
0610     retVal = kFailed;
0611   }
0612
0613
0614   if (!fx) {
0615     // compare
0616     if (compareMatrices(jac_f, jac_f_num, deltaJac) == kFailed) {
0617       if (verbose) {
0618         std::cout << "jac_f = ";
0619         jac_f.Print();
0620         std::cout << "jac_f_num = ";
0621         jac_f_num.Print();
0622         std::cout << std::endl;
0623       }
0624       retVal = kFailed;
0625     }
0626   }
0627
0628   delete rep;
0629   genfit::MaterialEffects::getInstance()->setNoEffects(false);
0630
0631   return retVal;
0632 }
0633
0634
0635 e_testStatus checkStopAtBoundary(bool writeHisto = false) {
0636
0637   if (writeHisto)
0638     return kPassed;
0639
0640   double epsilonLen = 1.E-4; // 1 mu
0641   //double epsilonMom = 1.E-5; // 10 keV
0642
0643   double matRadius(1.);
0644
0645   int pdg = randomPdg();
0646   genfit::AbsTrackRep* rep;
0647   rep = new genfit::RKTrackRep(pdg);
0648
0649   //TVector3 pos(0,0,0);
0650   TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
0651   TVector3 mom(0,0.5,gRandom->Gaus(0,0.1));
0652   mom *= randomSign();
0653
0654
0655   genfit::StateOnPlane state(rep);
0656   rep->setPosMom(state, pos, mom);
0657
0658   genfit::SharedPlanePtr origPlane = state.getPlane();
0659   genfit::SharedPlanePtr plane(new genfit::DetPlane(TVector3(0,randomSign()*10,0), TVector3(0,randomSign()*1,0)));
0660
0661   genfit::StateOnPlane origState(state);
0662
0663   // forth
0664   try {
0665     rep->extrapolateToPlane(state, plane, true);
0666   }
0667   catch (genfit::Exception& e) {
0668     std::cerr << e.what();
0669
0670     delete rep;
0671     return kException;
0672   }
0673
0674
0675   // compare
0676   if (fabs(rep->getPos(state).Perp() - matRadius) > epsilonLen) {
0677       if (verbose) {
0678         origState.Print();
0679         state.Print();
0680
0681         std::cerr << "radius difference = " << rep->getPos(state).Perp() - matRadius << "\n";
0682
0683         std::cerr << std::endl;
0684       }
0685       delete rep;
0686       return kFailed;
0687     }
0688
0689     delete rep;
0690     return kPassed;
0691
0692 }
0693
0694
0695 e_testStatus checkErrorPropagation(bool writeHisto = false) {
0696
0697   if (writeHisto)
0698     return kPassed;
0699
0700   //double epsilonLen = 1.E-4; // 1 mu
0701   //double epsilonMom = 1.E-5; // 10 keV
0702
0703   int pdg = randomPdg();
0704   genfit::AbsTrackRep* rep;
0705   rep = new genfit::RKTrackRep(pdg);
0706
0707   //TVector3 pos(0,0,0);
0708   TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
0709   TVector3 mom(0,0.5,gRandom->Gaus(0,0.1));
0710   mom *= randomSign();
0711
0712
0713   genfit::MeasuredStateOnPlane state(rep);
0714   rep->setPosMom(state, pos, mom);
0715
0716   genfit::SharedPlanePtr origPlane = state.getPlane();
0717   genfit::SharedPlanePtr plane(new genfit::DetPlane(TVector3(0,randomSign()*50,0), TVector3(0,randomSign()*1,0)));
0718
0719   genfit::MeasuredStateOnPlane origState(state);
0720
0721   // forth
0722   try {
0723     rep->extrapolateToPlane(state, plane);
0724   }
0725   catch (genfit::Exception& e) {
0726     std::cerr << e.what();
0727
0728     delete rep;
0729     return kException;
0730   }
0731
0732
0733   // check
0734   if (isCovMatrix(state.getCov()) == kFailed) {
0735     if (verbose) {
0736       origState.Print();
0737       state.Print();
0738     }
0739     delete rep;
0740     return kFailed;
0741   }
0742
0743   delete rep;
0744   return kPassed;
0745
0746 }
0747
0748
0749 e_testStatus checkExtrapolateToLine(bool writeHisto = false) {
0750
0751   if (writeHisto)
0752     return kPassed;
0753
0754   double epsilonLen = 1.E-4; // 1 mu
0755   double epsilonMom = 1.E-5; // 10 keV
0756
0757   int pdg = randomPdg();
0758   genfit::AbsTrackRep* rep;
0759   rep = new genfit::RKTrackRep(pdg);
0760
0761   //TVector3 pos(0,0,0);
0762   TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
0763   TVector3 mom(0,0.5,0.);
0764   mom *= randomSign();
0765
0766
0767   genfit::StateOnPlane state(rep);
0768   rep->setPosMom(state, pos, mom);
0769
0770   genfit::SharedPlanePtr origPlane = state.getPlane();
0771   genfit::StateOnPlane origState(state);
0772
0773   TVector3 linePoint(gRandom->Gaus(),randomSign()*10+gRandom->Gaus(),gRandom->Gaus());
0774   TVector3 lineDirection(gRandom->Gaus(),gRandom->Gaus(),randomSign()*10+gRandom->Gaus());
0775
0776   // forth
0777   try {
0778     rep->extrapolateToLine(state, linePoint, lineDirection, false);
0779   }
0780   catch (genfit::Exception& e) {
0781     std::cerr << e.what();
0782
0783     delete rep;
0784     return kException;
0785   }
0786
0787
0788   // compare
0789   if (fabs(state.getPlane()->distance(linePoint)) > epsilonLen ||
0790       fabs(state.getPlane()->distance(linePoint+lineDirection)) > epsilonLen ||
0791       (rep->getMom(state).Unit() * state.getPlane()->getNormal()) > epsilonMom) {
0792
0793       if (verbose) {
0794         origState.Print();
0795         state.Print();
0796
0797         std::cout << "distance of linePoint to plane = " << state.getPlane()->distance(linePoint) << "\n";
0798         std::cout << "distance of linePoint+lineDirection to plane = "
0799                   << state.getPlane()->distance(linePoint + lineDirection) << "\n";
0800         std::cout << "direction * plane normal = " << rep->getMom(state).Unit() * state.getPlane()->getNormal() << "\n";
0801       }
0802       delete rep;
0803       return kFailed;
0804     }
0805
0806     delete rep;
0807     return kPassed;
0808
0809 }
0810
0811
0812 e_testStatus checkExtrapolateToPoint(bool writeHisto = false) {
0813
0814   if (writeHisto)
0815     return kPassed;
0816
0817   double epsilonLen = 1.E-4; // 1 mu
0818   double epsilonMom = 1.E-5; // 10 keV
0819
0820   int pdg = randomPdg();
0821   genfit::AbsTrackRep* rep;
0822   rep = new genfit::RKTrackRep(pdg);
0823
0824   //TVector3 pos(0,0,0);
0825   TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
0826   TVector3 mom(0,0.5,gRandom->Gaus(0,0.1));
0827   mom *= randomSign();
0828
0829
0830   genfit::StateOnPlane state(rep);
0831   rep->setPosMom(state, pos, mom);
0832
0833   genfit::SharedPlanePtr origPlane = state.getPlane();
0834   genfit::StateOnPlane origState(state);
0835
0836   TVector3 point(gRandom->Gaus(),randomSign()*10+gRandom->Gaus(),gRandom->Gaus());
0837
0838   // forth
0839   try {
0840     rep->extrapolateToPoint(state, point, false);
0841   }
0842   catch (genfit::Exception& e) {
0843     std::cerr << e.what();
0844
0845     delete rep;
0846     return kException;
0847   }
0848
0849
0850   // compare
0851   if (fabs(state.getPlane()->distance(point)) > epsilonLen ||
0852       fabs((rep->getMom(state).Unit() * state.getPlane()->getNormal())) - 1 > epsilonMom) {
0853       if (verbose) {
0854         origState.Print();
0855         state.Print();
0856
0857         std::cout << "distance of point to plane = " << state.getPlane()->distance(point) << "\n";
0858         std::cout << "direction * plane normal = " << rep->getMom(state).Unit() * state.getPlane()->getNormal() << "\n";
0859       }
0860       delete rep;
0861       return kFailed;
0862     }
0863
0864     delete rep;
0865     return kPassed;
0866
0867 }
0868
0869
0870 e_testStatus checkExtrapolateToCylinder(bool writeHisto = false) {
0871
0872   if (writeHisto)
0873     return kPassed;
0874
0875   double epsilonLen = 1.E-4; // 1 mu
0876   //double epsilonMom = 1.E-5; // 10 keV
0877
0878   int pdg = randomPdg();
0879   genfit::AbsTrackRep* rep;
0880   rep = new genfit::RKTrackRep(pdg);
0881
0882   //TVector3 pos(0,0,0);
0883   TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
0884   TVector3 mom(0,0.5,0.);
0885   mom *= randomSign();
0886
0887
0888   genfit::StateOnPlane state(rep);
0889   rep->setPosMom(state, pos, mom);
0890
0891   genfit::SharedPlanePtr origPlane = state.getPlane();
0892   genfit::StateOnPlane origState(state);
0893
0894   const TVector3 linePoint(gRandom->Gaus(0,5), gRandom->Gaus(0,5), gRandom->Gaus(0,5));
0895   const TVector3 lineDirection(gRandom->Gaus(),gRandom->Gaus(),2+gRandom->Gaus());
0896   const double radius = gRandom->Uniform(10);
0897
0898   // forth
0899   try {
0900     rep->extrapolateToCylinder(state, radius, linePoint, lineDirection, false);
0901   }
0902   catch (genfit::Exception& e) {
0903     std::cerr << e.what();
0904
0905     delete rep;
0906
0907     return kException;
0908   }
0909
0910   TVector3 pocaOnLine(lineDirection);
0911   double t = 1./(pocaOnLine.Mag2()) * ((rep->getPos(state)*pocaOnLine) - (linePoint*pocaOnLine));
0912   pocaOnLine *= t;
0913   pocaOnLine += linePoint;
0914
0915   TVector3 radiusVec = rep->getPos(state) - pocaOnLine;
0916
0917   // compare
0918   if (fabs(state.getPlane()->getNormal()*radiusVec.Unit())-1 > epsilonLen ||
0919       fabs(lineDirection*radiusVec) > epsilonLen ||
0920       fabs(radiusVec.Mag()-radius) > epsilonLen) {
0921       if (verbose) {
0922         origState.Print();
0923         state.Print();
0924
0925         std::cout << "lineDirection*radiusVec = " << lineDirection * radiusVec << "\n";
0926         std::cout << "radiusVec.Mag()-radius = " << radiusVec.Mag() - radius << "\n";
0927       }
0928       delete rep;
0929       return kFailed;
0930     }
0931
0932     delete rep;
0933     return kPassed;
0934
0935 }
0936
0937
0938 e_testStatus checkExtrapolateToSphere(bool writeHisto = false) {
0939
0940   if (writeHisto)
0941     return kPassed;
0942
0943   double epsilonLen = 1.E-4; // 1 mu
0944   //double epsilonMom = 1.E-5; // 10 keV
0945
0946   int pdg = randomPdg();
0947   genfit::AbsTrackRep* rep;
0948   rep = new genfit::RKTrackRep(pdg);
0949
0950   //TVector3 pos(0,0,0);
0951   TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
0952   TVector3 mom(0,0.5,0.);
0953   mom *= randomSign();
0954
0955
0956   genfit::StateOnPlane state(rep);
0957   rep->setPosMom(state, pos, mom);
0958
0959   genfit::SharedPlanePtr origPlane = state.getPlane();
0960   genfit::StateOnPlane origState(state);
0961
0962   const TVector3 centerPoint(gRandom->Gaus(0,10), gRandom->Gaus(0,10), gRandom->Gaus(0,10));
0963   const double radius = gRandom->Uniform(10);
0964
0965   // forth
0966   try {
0967     rep->extrapolateToSphere(state, radius, centerPoint, false);
0968   }
0969   catch (genfit::Exception& e) {
0970     std::cerr << e.what();
0971
0972     delete rep;
0973
0974     return kException;
0975   }
0976
0977
0978   TVector3 radiusVec = rep->getPos(state) - centerPoint;
0979
0980   // compare
0981   if (fabs(state.getPlane()->getNormal()*radiusVec.Unit())-1 > epsilonLen ||
0982       fabs(radiusVec.Mag()-radius) > epsilonLen) {
0983       if (verbose) {
0984         origState.Print();
0985         state.Print();
0986
0987         std::cout << "state.getPlane()->getNormal()*radiusVec = " << state.getPlane()->getNormal() * radiusVec << "\n";
0988         std::cout << "radiusVec.Mag()-radius = " << radiusVec.Mag() - radius << "\n";
0989       }
0990       delete rep;
0991       return kFailed;
0992     }
0993
0994     delete rep;
0995     return kPassed;
0996
0997 }
0998
0999
1000 e_testStatus checkExtrapolateBy(bool writeHisto = false) {
1001
1002   if (writeHisto)
1003     return kPassed;
1004
1005   double epsilonLen = 1.E-3; // 10 mu
1006
1007   int pdg = randomPdg();
1008   genfit::AbsTrackRep* rep;
1009   rep = new genfit::RKTrackRep(pdg);
1010
1011   //TVector3 pos(0,0,0);
1012   TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
1013   TVector3 mom(0,0.5,0.);
1014   mom *= randomSign();
1015
1016
1017   genfit::StateOnPlane state(rep);
1018   rep->setPosMom(state, pos, mom);
1019
1020   TVector3 posOrig(state.getPos());
1021
1022   genfit::SharedPlanePtr origPlane = state.getPlane();
1023   genfit::StateOnPlane origState(state);
1024
1025   double step = gRandom->Uniform(-15.,15.);
1026   double extrapolatedLen(0);
1027
1028   // forth
1029   try {
1030     extrapolatedLen = rep->extrapolateBy(state, step, false);
1031   }
1032   catch (genfit::Exception& e) {
1033     return kException;
1034   }
1035
1036   TVector3 posExt(state.getPos());
1037
1038
1039
1040
1041   // compare
1042   if (fabs(extrapolatedLen-step) > epsilonLen ||
1043       (posOrig - posExt).Mag() > fabs(step)) {
1044       if (verbose) {
1045         origState.Print();
1046         state.Print();
1047
1048         std::cout << "extrapolatedLen-step = " << extrapolatedLen - step << "\n";
1049         std::cout << "started extrapolation from: ";
1050         posOrig.Print();
1051         std::cout << "extrapolated to ";
1052         posExt.Print();
1053         std::cout << "difference = " << (posOrig - posExt).Mag() << "; step = " << step << "; delta = "
1054                   << (posOrig - posExt).Mag() - fabs(step) << "\n";
1055       }
1056       delete rep;
1057       return kFailed;
1058     }
1059
1060     delete rep;
1061     return kPassed;
1062
1063 }
1064 //=====================================================================================================================
1065 //=====================================================================================================================
1066 //=====================================================================================================================
1067 //=====================================================================================================================
1068
1069 struct TestCase {
1070   TestCase(const std::string &name, e_testStatus(*function)(bool)) : name_(name), function_(function), nPassed_(0), nFailed_(0), nException_(0) {;}
1071   void Print() {std::cout << name_ << " \t" << nPassed_ << " \t" << nFailed_ << " \t" << nException_ << "\n";}
1072
1073   std::string name_;
1074   e_testStatus(*function_)(bool);
1075   unsigned int nPassed_;
1076   unsigned int nFailed_;
1077   unsigned int nException_;
1078 };
1079
1080
1081 int main() {
1082
1083   const double BField = 15.;       // kGauss
1084   //const bool debug = true;
1085
1086   gRandom->SetSeed(10);
1087   signal(SIGSEGV, handler);   // install our handler
1088
1089   // init geometry and mag. field
1090   new TGeoManager("Geometry", "Geane geometry");
1091   TGeoManager::Import("genfitGeom.root");
1092   genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0.,BField));
1093   genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
1094
1095   /*genfit::MaterialEffects::getInstance()->drawdEdx(2212);
1096   genfit::MaterialEffects::getInstance()->drawdEdx(11);
1097   genfit::MaterialEffects::getInstance()->drawdEdx(211);
1098   genfit::MaterialEffects::getInstance()->drawdEdx(13);
1099   return 0;*/
1100
1101   TDatabasePDG::Instance()->GetParticle(211);
1102
1103   const unsigned int nTests(1000);
1104
1105   std::vector<TestCase> testCases;
1106   testCases.push_back(TestCase(std::string("checkSetGetPosMom()            "), &checkSetGetPosMom));
1107   testCases.push_back(TestCase(std::string("compareForthBackExtrapolation()"), &compareForthBackExtrapolation));
1108   testCases.push_back(TestCase(std::string("checkStopAtBoundary()          "), &checkStopAtBoundary));
1109   testCases.push_back(TestCase(std::string("checkErrorPropagation()        "), &checkErrorPropagation));
1110   testCases.push_back(TestCase(std::string("checkExtrapolateToLine()       "), &checkExtrapolateToLine));
1111   testCases.push_back(TestCase(std::string("checkExtrapolateToPoint()      "), &checkExtrapolateToPoint));
1112   testCases.push_back(TestCase(std::string("checkExtrapolateToCylinder()   "), &checkExtrapolateToCylinder));
1113   testCases.push_back(TestCase(std::string("checkExtrapolateToSphere()     "), &checkExtrapolateToSphere));
1114   testCases.push_back(TestCase(std::string("checkExtrapolateBy()           "), &checkExtrapolateBy));
1115   testCases.push_back(TestCase(std::string("compareForthBackJacNoise()     "), &compareForthBackJacNoise));
1116
1117
1118   for (unsigned int i=0; i<nTests; ++i) {
1119
1120     for (unsigned int j=0; j<testCases.size(); ++j) {
1121       e_testStatus status = testCases[j].function_(false);
1122       switch (status) {
1123       case kPassed:
1124         testCases[j].nPassed_++;
1125         break;
1126       case kFailed:
1127         testCases[j].nFailed_++;
1128         std::cout << "failed " << testCases[j].name_ << " nr " << i << "\n";
1129         break;
1130       case kException:
1131         testCases[j].nException_++;
1132         std::cout << "exception at " << testCases[j].name_ << " nr " << i << "\n";
1133       }
1134     }
1135
1136   }
1137
1138   //CALLGRIND_STOP_INSTRUMENTATION;
1139   std::cout << "name                           " << " \t" << "pass" << " \t" << "fail" << " \t" << "exception" << "\n";
1140   for (unsigned int j=0; j<testCases.size(); ++j) {
1141     testCases[j].Print();
1142   }
1143
1144   for (unsigned int j=0; j<testCases.size(); ++j) {
1145     testCases[j].function_(true);
1146   }
1147
1148   return 0;
1149 }
1150
1151