test/fitterTests/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 <AbsKalmanFitter.h>
0016 #include <KalmanFitter.h>
0017 #include <KalmanFitterRefTrack.h>
0018 #include <KalmanFitterInfo.h>
0019 #include <KalmanFitStatus.h>
0020 #include <DAF.h>
0021 #include <GFGbl.h>
0022 #include <MeasuredStateOnPlane.h>
0023 #include <MeasurementOnPlane.h>
0024 #include <FullMeasurement.h>
0025 #include <PlanarMeasurement.h>
0026 #include <ProlateSpacepointMeasurement.h>
0027 #include <RectangularFinitePlane.h>
0028 #include <ReferenceStateOnPlane.h>
0029 #include <SharedPlanePtr.h>
0030 #include <SpacepointMeasurement.h>
0031 #include <StateOnPlane.h>
0032 #include <Tools.h>
0033 #include <TrackCand.h>
0034 #include <TrackCandHit.h>
0035 #include <Track.h>
0036 #include <TrackPoint.h>
0037 #include <WireMeasurement.h>
0038 #include <WirePointMeasurement.h>
0039
0040 #include <MaterialEffects.h>
0041 #include <RKTools.h>
0042 #include <RKTrackRep.h>
0043 #include <StepLimits.h>
0044 #include <TGeoMaterialInterface.h>
0045
0046 #include <EventDisplay.h>
0047
0048 #include <HelixTrackModel.h>
0049 #include <MeasurementCreator.h>
0050
0051 #include <TApplication.h>
0052 #include <TCanvas.h>
0053 #include <TDatabasePDG.h>
0054 #include <TEveManager.h>
0055 #include <TGeoManager.h>
0056 #include <TH1D.h>
0057 #include <TRandom.h>
0058 #include <TStyle.h>
0059 #include <TVector3.h>
0060 #include <vector>
0061
0062 #include <TROOT.h>
0063 #include <TFile.h>
0064 #include <TTree.h>
0065 #include "TDatabasePDG.h"
0066 #include <TMath.h>
0067 #include <TString.h>
0068
0069 #include <memory>
0070
0071
0072 void handler(int sig) {
0073   void *array[10];
0074   size_t size;
0075
0076   // get void*'s for all entries on the stack
0077   size = backtrace(array, 10);
0078
0079   // print out all the frames to stderr
0080   fprintf(stderr, "Error: signal %d:\n", sig);
0081   backtrace_symbols_fd(array, size, 2);
0082   exit(1);
0083 }
0084
0085 int randomPdg() {
0086   int pdg;
0087
0088   switch(int(gRandom->Uniform(8))) {
0089   case 1:
0090     pdg = -11; break;
0091   case 2:
0092     pdg = 11; break;
0093   case 3:
0094     pdg = 13; break;
0095   case 4:
0096     pdg = -13; break;
0097   case 5:
0098     pdg = 211; break;
0099   case 6:
0100     pdg = -211; break;
0101   case 7:
0102     pdg = 2212; break;
0103   default:
0104     pdg = 211;
0105   }
0106
0107   return pdg;
0108 }
0109
0110
0111 int randomSign() {
0112   if (gRandom->Uniform(1) > 0.5)
0113     return 1;
0114   return -1;
0115 }
0116
0117 //---------------------------------------------------------------------------------------------------------
0118 //---------------------------------------------------------------------------------------------------------
0119 //---------------------------------------------------------------------------------------------------------
0120 //---------------------------------------------------------------------------------------------------------
0121 //---------------------------------------------------------------------------------------------------------
0122
0123 //#define VALGRIND
0124
0125 #ifdef VALGRIND
0126   #include <valgrind/callgrind.h>
0127 #else
0128 #define CALLGRIND_START_INSTRUMENTATION
0129 #define CALLGRIND_STOP_INSTRUMENTATION
0130 #define CALLGRIND_DUMP_STATS
0131 #endif
0132
0133 int main() {
0134   std::cout<<"main"<<std::endl;
0135   gRandom->SetSeed(14);
0136
0137
0138   const unsigned int nEvents = 1000;
0139   const unsigned int nMeasurements = 10;
0140   const double BField = 20.;       // kGauss
0141   const double momentum = 0.1;     // GeV
0142   const double theta = 110;         // degree
0143   const double thetaDetPlane = 90;         // degree
0144   const double phiDetPlane = 0;         // degree
0145   const double pointDist = 3.;      // cm; approx. distance between measurements
0146   const double resolution = 0.05;   // cm; resolution of generated measurements
0147
0148   const double resolutionWire = 5*resolution;   // cm; resolution of generated measurements
0149   const TVector3 wireDir(0,0,1);
0150   const double skewAngle(5);
0151   const bool useSkew = true;
0152   const int nSuperLayer = 10;
0153   const double minDrift = 0.;
0154   const double maxDrift = 2;
0155   const bool idealLRResolution = false; // resolve the l/r ambiguities of the wire measurements
0156
0157   const double outlierProb = -0.1;
0158   const double outlierRange = 2;
0159
0160   const double hitSwitchProb = -0.1; // probability to give hits to fit in wrong order (flip two hits)
0161
0162   const int splitTrack = -5; //nMeasurements/2; // for track merging testing.
0163   const bool fullMeasurement = false; // put fit result of first tracklet as FullMeasurement into second tracklet, don't merge
0164
0165   //const genfit::eFitterType fitterId = genfit::SimpleKalman;
0166   const genfit::eFitterType fitterId = genfit::RefKalman;
0167   //const genfit::eFitterType fitterId = genfit::DafRef;
0168   //const genfit::eFitterType fitterId = genfit::DafSimple;
0169   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedAverage;
0170   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToReference;
0171   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToPrediction;
0172   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToReference;
0173   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToPrediction;
0174   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToReferenceWire;
0175   const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToPredictionWire;
0176   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToReferenceWire;
0177   //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToPredictionWire;
0178   const int nIter = 20; // max number of iterations
0179   const double dPVal = 1.E-3; // convergence criterion
0180
0181   const bool resort = false;
0182   const bool prefit = false; // make a simple Kalman iteration before the actual fit
0183   const bool refit  = false; // if fit did not converge, try to fit again
0184
0185   const bool twoReps = false; // test if everything works with more than one rep in the tracks
0186
0187   const bool checkPruning = true; // test pruning
0188
0189   const int pdg = 13;               // particle pdg code
0190
0191   const bool smearPosMom = true;     // init the Reps with smeared pos and mom
0192   const double chargeSwitchProb = -0.1; // probability to seed with wrong charge sign
0193   const double posSmear = 10*resolution;     // cm
0194   const double momSmear = 5. /180.*TMath::Pi();     // rad
0195   const double momMagSmear = 0.2;   // relative
0196   const double zSmearFac = 2;
0197
0198
0199   const bool matFX = false;         // include material effects; can only be disabled for RKTrackRep!
0200
0201   const bool debug = false;
0202   const bool onlyDisplayFailed = false; // only load non-converged tracks into the display
0203
0204   std::vector<genfit::eMeasurementType> measurementTypes;
0205   for (unsigned int i = 0; i<nMeasurements; ++i) {
0206     measurementTypes.push_back(genfit::eMeasurementType(i%8));
0207   }
0208
0209   signal(SIGSEGV, handler);   // install our handler
0210
0211   // init fitter
0212   genfit::AbsKalmanFitter* fitter = 0;
0213   switch (fitterId) {
0214     case genfit::SimpleKalman:
0215       fitter = new genfit::KalmanFitter(nIter, dPVal);
0216       fitter->setMultipleMeasurementHandling(mmHandling);
0217       break;
0218
0219     case genfit::RefKalman:
0220       fitter = new genfit::KalmanFitterRefTrack(nIter, dPVal);
0221       fitter->setMultipleMeasurementHandling(mmHandling);
0222       break;
0223
0224     case genfit::DafSimple:
0225       fitter = new genfit::DAF(false);
0226       break;
0227     case genfit::DafRef:
0228       fitter = new genfit::DAF();
0229       break;
0230   }
0231   fitter->setMaxIterations(nIter);
0232   //if (debug)
0233   //  fitter->setDebugLvl(10);
0234
0235   /*if (dynamic_cast<genfit::DAF*>(fitter) != nullptr) {
0236     //static_cast<genfit::DAF*>(fitter)->setBetas(100, 50, 25, 12, 6, 3, 1, 0.5, 0.1);
0237     //static_cast<genfit::DAF*>(fitter)->setBetas(81, 8, 4, 0.5, 0.1);
0238     static_cast<genfit::DAF*>(fitter)->setAnnealingScheme(100, 0.1, 5);
0239     //static_cast<genfit::DAF*>(fitter)->setConvergenceDeltaWeight(0.0001);
0240     //fitter->setMaxIterations(nIter);
0241   }*/
0242
0243
0244   // init MeasurementCreator
0245   genfit::MeasurementCreator measurementCreator;
0246   measurementCreator.setResolution(resolution);
0247   measurementCreator.setResolutionWire(resolutionWire);
0248   measurementCreator.setOutlierProb(outlierProb);
0249   measurementCreator.setOutlierRange(outlierRange);
0250   measurementCreator.setThetaDetPlane(thetaDetPlane);
0251   measurementCreator.setPhiDetPlane(phiDetPlane);
0252   measurementCreator.setWireDir(wireDir);
0253   measurementCreator.setMinDrift(minDrift);
0254   measurementCreator.setMaxDrift(maxDrift);
0255   measurementCreator.setIdealLRResolution(idealLRResolution);
0256   measurementCreator.setUseSkew(useSkew);
0257   measurementCreator.setSkewAngle(skewAngle);
0258   measurementCreator.setNSuperLayer(nSuperLayer);
0259   measurementCreator.setDebug(debug);
0260
0261
0262   // init geometry and mag. field
0263   new TGeoManager("Geometry", "Geane geometry");
0264   TGeoManager::Import("genfitGeom.root");
0265   genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0.,BField));
0266   genfit::FieldManager::getInstance()->useCache(true, 8);
0267   genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
0268
0269   const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
0270
0271   // init event display
0272 #ifndef VALGRIND
0273   genfit::EventDisplay* display = genfit::EventDisplay::getInstance();
0274   display->reset();
0275 #endif
0276
0277
0278 #ifndef VALGRIND
0279   // create histograms
0280   gROOT->SetStyle("Plain");
0281   gStyle->SetPalette(1);
0282   gStyle->SetOptFit(1111);
0283
0284   TH1D *hmomRes = new TH1D("hmomRes","mom res",500,-20*resolution*momentum/nMeasurements,20*resolution*momentum/nMeasurements);
0285   TH1D *hupRes = new TH1D("hupRes","u' res",500,-15*resolution/nMeasurements, 15*resolution/nMeasurements);
0286   TH1D *hvpRes = new TH1D("hvpRes","v' res",500,-15*resolution/nMeasurements, 15*resolution/nMeasurements);
0287   TH1D *huRes = new TH1D("huRes","u res",500,-15*resolution, 15*resolution);
0288   TH1D *hvRes = new TH1D("hvRes","v res",500,-15*resolution, 15*resolution);
0289
0290   TH1D *hqopPu = new TH1D("hqopPu","q/p pull",200,-6.,6.);
0291   TH1D *pVal = new TH1D("pVal","p-value",100,0.,1.00000001);
0292   pVal->SetMinimum(0);
0293   TH1D *hupPu = new TH1D("hupPu","u' pull",200,-6.,6.);
0294   TH1D *hvpPu = new TH1D("hvpPu","v' pull",200,-6.,6.);
0295   TH1D *huPu = new TH1D("huPu","u pull",200,-6.,6.);
0296   TH1D *hvPu = new TH1D("hvPu","v pull",200,-6.,6.);
0297
0298   TH1D *weights = new TH1D("weights","Daf vs true weights",500,-1.01,1.01);
0299
0300   TH1D *trackLenRes = new TH1D("trackLenRes","(trueLen - FittedLen) / trueLen",500,-0.01,0.01);
0301 #endif
0302
0303   double maxWeight(0);
0304   unsigned int nTotalIterConverged(0);
0305   unsigned int nTotalIterNotConverged(0);
0306   unsigned int nTotalIterSecondConverged(0);
0307   unsigned int nTotalIterSecondNotConverged(0);
0308   unsigned int nConvergedFits(0);
0309   unsigned int nUNConvergedFits(0);
0310   unsigned int nConvergedFitsSecond(0);
0311   unsigned int nUNConvergedFitsSecond(0);
0312
0313
0314   CALLGRIND_START_INSTRUMENTATION;
0315
0316   genfit::Track* fitTrack(nullptr);
0317   genfit::Track* secondTrack(nullptr);
0318
0319   // main loop
0320   for (unsigned int iEvent=0; iEvent<nEvents; ++iEvent){
0321
0322       /*measurementTypes.clear();
0323       for (unsigned int i = 0; i < nMeasurements; ++i)
0324         measurementTypes.push_back(genfit::eMeasurementType(gRandom->Uniform(8)));*/
0325
0326
0327       if (debug || (iEvent)%10==0)
0328         std::cout << "Event Nr. " << iEvent << " ";
0329       else std::cout << ". ";
0330       if (debug || (iEvent+1)%10==0)
0331         std::cout << "\n";
0332
0333
0334       // clean up
0335       delete fitTrack;
0336       fitTrack = nullptr;
0337       delete secondTrack;
0338       secondTrack = nullptr;
0339
0340
0341       // true start values
0342       TVector3 pos(0, 0, 0);
0343       TVector3 mom(1.,0,0);
0344       mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
0345       //mom.SetTheta(gRandom->Uniform(0.5*TMath::Pi(),0.9*TMath::Pi()));
0346       mom.SetTheta(theta*TMath::Pi()/180);
0347       mom.SetMag(momentum);
0348       TMatrixDSym covM(6);
0349       for (int i = 0; i < 3; ++i)
0350         covM(i,i) = resolution*resolution;
0351       for (int i = 3; i < 6; ++i)
0352         covM(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
0353
0354       if (debug) {
0355         std::cout << "start values \n";
0356         pos.Print();
0357         mom.Print();
0358       }
0359
0360       // calc helix parameters
0361       genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
0362       measurementCreator.setTrackModel(helix);
0363
0364       // smeared start values
0365       TVector3 posM(pos);
0366       TVector3 momM(mom);
0367       if (smearPosMom) {
0368         posM.SetX(gRandom->Gaus(posM.X(),posSmear));
0369         posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
0370         posM.SetZ(gRandom->Gaus(posM.Z(),zSmearFac*posSmear));
0371
0372         momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
0373         momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
0374         momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
0375       }
0376
0377
0378       // trackrep for creating measurements
0379       double sign(1.);
0380       if (chargeSwitchProb > gRandom->Uniform(1.))
0381         sign = -1.;
0382       genfit::AbsTrackRep* rep = new genfit::RKTrackRep(sign*pdg);
0383       sign = 1.;
0384       if (chargeSwitchProb > gRandom->Uniform(1.))
0385         sign = -1.;
0386       genfit::AbsTrackRep* secondRep = new genfit::RKTrackRep(sign*-211);
0387       genfit::MeasuredStateOnPlane stateRef(rep);
0388       rep->setPosMomCov(stateRef, pos, mom, covM);
0389
0390       // smeared start state
0391       genfit::MeasuredStateOnPlane stateSmeared(rep);
0392       rep->setPosMomCov(stateSmeared, posM, momM, covM);
0393
0394       //rep->setPropDir(1);
0395
0396       if (!matFX) genfit::MaterialEffects::getInstance()->setNoEffects();
0397
0398       // remember original initial state
0399       const genfit::StateOnPlane stateRefOrig(stateRef);
0400
0401       // create smeared measurements
0402       std::vector< std::vector<genfit::AbsMeasurement*> > measurements;
0403
0404       std::vector<bool> outlierTrue;
0405       bool outlier;
0406       // true values for left right. 0 for non wire measurements
0407       std::vector<int> leftRightTrue;
0408       int lr;
0409
0410       double trueLen(-1);
0411
0412       try{
0413         for (unsigned int i=0; i<measurementTypes.size(); ++i){
0414           trueLen = i*pointDist;
0415
0416           measurements.push_back(measurementCreator.create(measurementTypes[i], trueLen, outlier, lr));
0417           outlierTrue.push_back(outlier);
0418           leftRightTrue.push_back(lr);
0419         }
0420         assert(measurementTypes.size() == leftRightTrue.size());
0421         assert(measurementTypes.size() == outlierTrue.size());
0422       }
0423       catch(genfit::Exception& e){
0424         std::cerr<<"Exception, next track"<<std::endl;
0425         std::cerr << e.what();
0426         continue; // here is a memleak!
0427       }
0428
0429       if (debug) std::cout << "... done creating measurements \n";
0430
0431
0432
0433       // create track
0434       TVectorD seedState(6);
0435       TMatrixDSym seedCov(6);
0436       rep->get6DStateCov(stateSmeared, seedState, seedCov);
0437       fitTrack = new genfit::Track(rep, seedState, seedCov); //initialized with smeared rep
0438       secondTrack = new genfit::Track(rep->clone(), seedState, seedCov); //initialized with smeared rep
0439       if (twoReps) {
0440         fitTrack->addTrackRep(secondRep);
0441         secondTrack->addTrackRep(secondRep->clone());
0442       }
0443       else
0444         delete secondRep;
0445       //if (debug) fitTrack->Print("C");
0446
0447       fitTrack->checkConsistency();
0448       //fitTrack->addTrackRep(rep->clone()); // check if everything works fine with more than one rep
0449
0450       // add measurements
0451       for(unsigned int i=0; i<measurements.size(); ++i){
0452         if (splitTrack > 0 && (int)i >= splitTrack)
0453           break;
0454         if (i>0 && hitSwitchProb > gRandom->Uniform(1.))
0455           fitTrack->insertPoint(new genfit::TrackPoint(measurements[i], fitTrack), -2);
0456         else
0457           fitTrack->insertPoint(new genfit::TrackPoint(measurements[i], fitTrack));
0458
0459         fitTrack->checkConsistency();
0460         //if (debug) fitTrack->Print("C");
0461       }
0462
0463       if (splitTrack > 0) {
0464         for(unsigned int i=splitTrack; i<measurements.size(); ++i){
0465           if (i>0 && hitSwitchProb > gRandom->Uniform(1.))
0466             secondTrack->insertPoint(new genfit::TrackPoint(measurements[i], secondTrack), -2);
0467           else
0468             secondTrack->insertPoint(new genfit::TrackPoint(measurements[i], secondTrack));
0469
0470           //if (debug) fitTrack->Print("C");
0471         }
0472       }
0473
0474       fitTrack->checkConsistency();
0475       secondTrack->checkConsistency();
0476
0477       //if (debug) fitTrack->Print();
0478
0479       // do the fit
0480       try{
0481         if (debug) std::cout<<"Starting the fitter"<<std::endl;
0482
0483         if (prefit) {
0484           genfit::KalmanFitter prefitter(1, dPVal);
0485           prefitter.setMultipleMeasurementHandling(genfit::weightedClosestToPrediction);
0486           prefitter.processTrackWithRep(fitTrack, fitTrack->getCardinalRep());
0487         }
0488
0489         fitter->processTrack(fitTrack, resort);
0490         if (splitTrack > 0)
0491           fitter->processTrack(secondTrack, resort);
0492
0493         if (debug) std::cout<<"fitter is finished!"<<std::endl;
0494       }
0495       catch(genfit::Exception& e){
0496         std::cerr << e.what();
0497         std::cerr << "Exception, next track" << std::endl;
0498         continue;
0499       }
0500
0501       if (splitTrack > 0) {
0502         if (debug) fitTrack->Print("C");
0503         if (debug) secondTrack->Print("C");
0504
0505         if (fullMeasurement) {
0506           genfit::FullMeasurement* fullM = new genfit::FullMeasurement(secondTrack->getFittedState());
0507           fitTrack->insertPoint(new genfit::TrackPoint(fullM, fitTrack));
0508         }
0509         else
0510           fitTrack->mergeTrack(secondTrack);
0511
0512         if (debug) fitTrack->Print("C");
0513
0514         try{
0515           if (debug) std::cout<<"Starting the fitter"<<std::endl;
0516           fitter->processTrack(fitTrack, resort);
0517           if (debug) std::cout<<"fitter is finished!"<<std::endl;
0518         }
0519         catch(genfit::Exception& e){
0520           std::cerr << e.what();
0521           std::cerr << "Exception, next track" << std::endl;
0522           continue;
0523         }
0524       }
0525
0526
0527       if (refit && !fitTrack->getFitStatus(rep)->isFitConverged()) {
0528         std::cout<<"Trying to fit again "<<std::endl;
0529         fitter->processTrack(fitTrack, resort);
0530       }
0531
0532
0533
0534       if (debug) {
0535         fitTrack->Print("C");
0536         fitTrack->getFitStatus(rep)->Print();
0537       }
0538
0539       fitTrack->checkConsistency();
0540       secondTrack->checkConsistency();
0541
0542 #ifndef VALGRIND
0543       if (!onlyDisplayFailed && iEvent < 1000) {
0544         std::vector<genfit::Track*> event;
0545         event.push_back(fitTrack);
0546         if (splitTrack > 0)
0547           event.push_back(secondTrack);
0548         display->addEvent(event);
0549       }
0550       else if (onlyDisplayFailed &&
0551                (!fitTrack->getFitStatus(rep)->isFitConverged() ||
0552                 fitTrack->getFitStatus(rep)->getPVal() < 0.01)) {
0553         // add track to event display
0554         display->addEvent(fitTrack);
0555       }
0556 #endif
0557
0558
0559       if (fitTrack->getFitStatus(rep)->isFitConverged()) {
0560         nTotalIterConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(rep))->getNumIterations();
0561         nConvergedFits += 1;
0562       }
0563       else {
0564         nTotalIterNotConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(rep))->getNumIterations();
0565         nUNConvergedFits += 1;
0566       }
0567
0568       if (twoReps) {
0569         if (fitTrack->getFitStatus(secondRep)->isFitConverged()) {
0570           nTotalIterSecondConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(secondRep))->getNumIterations();
0571           nConvergedFitsSecond += 1;
0572         }
0573         else {
0574           nTotalIterSecondNotConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(secondRep))->getNumIterations();
0575           nUNConvergedFitsSecond += 1;
0576         }
0577       }
0578
0579
0580       // check if fit was successful
0581       if (! fitTrack->getFitStatus(rep)->isFitConverged()) {
0582         std::cout << "Track could not be fitted successfully! Fit is not converged! \n";
0583         continue;
0584       }
0585
0586
0587       genfit::TrackPoint* tp = fitTrack->getPointWithMeasurementAndFitterInfo(0, rep);
0588       if (tp == nullptr) {
0589         std::cout << "Track has no TrackPoint with fitterInfo! \n";
0590         continue;
0591       }
0592       genfit::KalmanFittedStateOnPlane kfsop(*(static_cast<genfit::KalmanFitterInfo*>(tp->getFitterInfo(rep))->getBackwardUpdate()));
0593       if (debug) {
0594         std::cout << "state before extrapolating back to reference plane \n";
0595         kfsop.Print();
0596       }
0597
0598       // extrapolate back to reference plane.
0599       try{
0600         rep->extrapolateToPlane(kfsop, stateRefOrig.getPlane());;
0601       }
0602       catch(genfit::Exception& e){
0603         std::cerr<<"Exception, next track"<<std::endl;
0604         std::cerr << e.what();
0605         continue;
0606       }
0607
0608 #ifndef VALGRIND
0609       // calculate pulls
0610       const TVectorD& referenceState = stateRefOrig.getState();
0611
0612       const TVectorD& state = kfsop.getState();
0613       const TMatrixDSym& cov = kfsop.getCov();
0614
0615       double pval = fitter->getPVal(fitTrack, rep);
0616       //assert( fabs(pval - static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(rep))->getBackwardPVal()) < 1E-10 );
0617
0618       hmomRes->Fill( (charge/state[0]-momentum));
0619       hupRes->Fill(  (state[1]-referenceState[1]));
0620       hvpRes->Fill(  (state[2]-referenceState[2]));
0621       huRes->Fill(   (state[3]-referenceState[3]));
0622       hvRes->Fill(   (state[4]-referenceState[4]));
0623
0624       hqopPu->Fill( (state[0]-referenceState[0]) / sqrt(cov[0][0]) );
0625       pVal->Fill(   pval);
0626       hupPu->Fill(  (state[1]-referenceState[1]) / sqrt(cov[1][1]) );
0627       hvpPu->Fill(  (state[2]-referenceState[2]) / sqrt(cov[2][2]) );
0628       huPu->Fill(   (state[3]-referenceState[3]) / sqrt(cov[3][3]) );
0629       hvPu->Fill(   (state[4]-referenceState[4]) / sqrt(cov[4][4]) );
0630
0631
0632       try {
0633         trackLenRes->Fill( (trueLen - fitTrack->getTrackLen(rep)) / trueLen );
0634
0635         if (debug) {
0636           std::cout << "true track length = " << trueLen << "; fitted length = " << fitTrack->getTrackLen(rep) << "\n";
0637           std::cout << "fitted tof = " << fitTrack->getTOF(rep) << " ns\n";
0638         }
0639       }
0640       catch (genfit::Exception& e) {
0641         std::cerr << e.what();
0642         std::cout << "could not get TrackLen or TOF! \n";
0643       }
0644
0645
0646
0647       // check l/r resolution and outlier rejection
0648       if (dynamic_cast<genfit::DAF*>(fitter) != nullptr) {
0649         for (unsigned int i=0; i<fitTrack->getNumPointsWithMeasurement(); ++i){
0650
0651           if (! fitTrack->getPointWithMeasurement(i)->hasFitterInfo(rep))
0652             continue;
0653
0654           if (debug) {
0655             std::vector<double> dafWeights = dynamic_cast<genfit::KalmanFitterInfo*>(fitTrack->getPointWithMeasurement(i)->getFitterInfo(rep))->getWeights();
0656             std::cout << "hit " << i;
0657             if (outlierTrue[i]) std::cout << " is an OUTLIER";
0658             std::cout << " weights: ";
0659             for (unsigned int j=0; j<dafWeights.size(); ++j){
0660               std::cout << dafWeights[j] << "  ";
0661             }
0662             std::cout << "   l/r: " << leftRightTrue[i];
0663             std::cout << "\n";
0664           }
0665           int trueSide = leftRightTrue[i];
0666           if (trueSide == 0) continue; // not a wire measurement
0667           if (outlierTrue[i]) continue; // an outlier
0668           std::vector<double> dafWeightLR = dynamic_cast<genfit::KalmanFitterInfo*>(fitTrack->getPointWithMeasurement(i)->getFitterInfo(rep))->getWeights();
0669           if(dafWeightLR.size() != 2)
0670             continue;
0671
0672           double weightCorrectSide, weightWrongSide;
0673
0674           if (trueSide < 0) {
0675             weightCorrectSide = dafWeightLR[0];
0676             weightWrongSide =  dafWeightLR[1];
0677           }
0678           else {
0679             weightCorrectSide = dafWeightLR[1];
0680             weightWrongSide =  dafWeightLR[0];
0681           }
0682           weightWrongSide -= 1.;
0683
0684           weights->Fill(weightCorrectSide);
0685           weights->Fill(weightWrongSide);
0686
0687           if (weightCorrectSide>maxWeight) maxWeight = weightCorrectSide;
0688         }
0689
0690         for (unsigned int i=0; i<fitTrack->getNumPointsWithMeasurement(); ++i){
0691           if (! fitTrack->getPointWithMeasurement(i)->hasFitterInfo(rep))
0692             continue;
0693
0694           std::vector<double> dafWeights = dynamic_cast<genfit::KalmanFitterInfo*>(fitTrack->getPointWithMeasurement(i)->getFitterInfo(rep))->getWeights();
0695
0696           if (outlierTrue[i]) { // an outlier
0697             for (unsigned int j=0; j<dafWeights.size(); ++j){
0698               weights->Fill(dafWeights[j]-1);
0699             }
0700           }
0701           else if (leftRightTrue[i] == 0) { // only for non wire hits
0702             for (unsigned int j=0; j<dafWeights.size(); ++j){
0703               weights->Fill(dafWeights[j]);
0704             }
0705           }
0706         }
0707
0708       }
0709
0710       if (checkPruning) { //check pruning
0711         //std::cout<<"\n";
0712         //std::cout<<"get stFirst ";
0713         genfit::MeasuredStateOnPlane stFirst = fitTrack->getFittedState();
0714         //std::cout<<"get stLast ";
0715         genfit::MeasuredStateOnPlane stLast = fitTrack->getFittedState(-1);
0716
0717         for (unsigned int i=0; i<1; ++i) {
0718           genfit::Track trClone(*fitTrack);
0719           trClone.checkConsistency();
0720
0721           bool first(false), last(false);
0722
0723           TString opt("");
0724           try {
0725             if (gRandom->Uniform() < 0.5) trClone.prune("C");
0726             if (gRandom->Uniform() < 0.5) {
0727               opt.Append("F");
0728               first = true;
0729             }
0730             if (gRandom->Uniform() < 0.5) {
0731               opt.Append("L");
0732               last = true;
0733             }
0734             if (gRandom->Uniform() < 0.5) opt.Append("W");
0735             if (gRandom->Uniform() < 0.5) opt.Append("R");
0736             if (gRandom->Uniform() < 0.5) opt.Append("M");
0737             if (gRandom->Uniform() < 0.5) opt.Append("I");
0738             if (gRandom->Uniform() < 0.5) opt.Append("U");
0739
0740             trClone.prune(opt);
0741
0742             try {
0743               trClone.checkConsistency();
0744             } catch (genfit::Exception& e) {
0745               trClone.getFitStatus()->getPruneFlags().Print();
0746             }
0747
0748             //std::cout<<"get stCloneFirst ";
0749             genfit::MeasuredStateOnPlane stCloneFirst = trClone.getFittedState();
0750             //std::cout<<"get stCloneLast ";
0751             genfit::MeasuredStateOnPlane stCloneLast = trClone.getFittedState(-1);
0752
0753             if (first and ! (stFirst.getState() == stCloneFirst.getState() and stFirst.getCov() == stCloneFirst.getCov() )) {
0754               //std::cout<<" failed first state ";
0755               //stFirst.Print();
0756               //stCloneFirst.Print();
0757
0758               if (debug)
0759                 trClone.getFitStatus()->getPruneFlags().Print();
0760             }
0761
0762             if (last  and ! (stLast.getState()  == stCloneLast.getState()  and stLast.getCov()  == stCloneLast.getCov() )) {
0763               //std::cout<<" failed last state ";
0764               //stLast.Print();
0765               //stCloneLast.Print();
0766
0767               if (debug)
0768                 trClone.getFitStatus()->getPruneFlags().Print();
0769             }
0770
0771             if (debug) {
0772               std::cout<<" pruned track: ";
0773               trClone.Print();
0774             }
0775           }
0776           catch (genfit::Exception &e) {
0777             std::cerr << e.what();
0778           }
0779         }
0780
0781       } // end check pruning
0782
0783 #endif
0784
0785
0786   }// end loop over events
0787
0788   delete fitTrack;
0789   delete secondTrack;
0790   delete fitter;
0791
0792   CALLGRIND_STOP_INSTRUMENTATION;
0793   CALLGRIND_DUMP_STATS;
0794
0795   std::cout<<"maxWeight = " << maxWeight << std::endl;
0796   std::cout<<"avg nr iterations =                     " << (double)(nTotalIterConverged + nTotalIterNotConverged)/(double)(nConvergedFits + nUNConvergedFits) << std::endl;
0797   std::cout<<"avg nr iterations of converged fits =   " << (double)(nTotalIterConverged)/(double)(nConvergedFits) << std::endl;
0798   std::cout<<"avg nr iterations of UNconverged fits = " << (double)(nTotalIterNotConverged)/(double)(nUNConvergedFits) << std::endl;
0799   std::cout<<"fit efficiency =                        " << (double)nConvergedFits/nEvents << std::endl;
0800
0801   if (twoReps) {
0802     std::cout<<"second rep: \navg nr iterations =                     " << (double)(nTotalIterSecondConverged + nTotalIterSecondNotConverged)/(double)(nConvergedFitsSecond + nUNConvergedFitsSecond) << std::endl;
0803     std::cout<<"avg nr iterations of converged fits =   " << (double)(nTotalIterSecondConverged)/(double)(nConvergedFitsSecond) << std::endl;
0804     std::cout<<"avg nr iterations of UNconverged fits = " << (double)(nTotalIterSecondNotConverged)/(double)(nUNConvergedFitsSecond) << std::endl;
0805     std::cout<<"fit efficiency =                        " << (double)nConvergedFitsSecond/nEvents << std::endl;
0806   }
0807
0808
0809   //std::cout<<"avg nr iterations (2nd rep) = " << (double)nTotalIterSecond/nSuccessfullFitsSecond << std::endl;
0810   //std::cout<<"fit efficiency (2nd rep) = " << (double)nConvergedFitsSecond/nEvents << std::endl;
0811
0812
0813 #ifndef VALGRIND
0814
0815   if (debug) std::cout<<"Draw histograms ...";
0816   // fit and draw histograms
0817   TCanvas* c1 = new TCanvas();
0818   c1->Divide(2,3);
0819
0820   c1->cd(1);
0821   hmomRes->Fit("gaus");
0822   hmomRes->Draw();
0823
0824   c1->cd(2);
0825   weights->Draw();
0826
0827   c1->cd(3);
0828   hupRes->Fit("gaus");
0829   hupRes->Draw();
0830
0831   c1->cd(4);
0832   hvpRes->Fit("gaus");
0833   hvpRes->Draw();
0834
0835   c1->cd(5);
0836   huRes->Fit("gaus");
0837   huRes->Draw();
0838
0839   c1->cd(6);
0840   hvRes->Fit("gaus");
0841   hvRes->Draw();
0842
0843   c1->Write();
0844
0845   TCanvas* c2 = new TCanvas();
0846   c2->Divide(2,3);
0847
0848   c2->cd(1);
0849   hqopPu->Fit("gaus");
0850   hqopPu->Draw();
0851
0852   c2->cd(2);
0853   pVal->Fit("pol1");
0854   pVal->Draw();
0855   c2->cd(3);
0856   hupPu->Fit("gaus");
0857   hupPu->Draw();
0858
0859   c2->cd(4);
0860   hvpPu->Fit("gaus");
0861   hvpPu->Draw();
0862
0863   c2->cd(5);
0864   huPu->Fit("gaus");
0865   huPu->Draw();
0866
0867   c2->cd(6);
0868   hvPu->Fit("gaus");
0869   hvPu->Draw();
0870
0871   c2->Write();
0872
0873
0874
0875   TCanvas* c3 = new TCanvas();
0876   //c3->Divide(2,3);
0877
0878   c3->cd(1);
0879   trackLenRes->Fit("gaus");
0880   trackLenRes->Draw();
0881
0882   c3->Write();
0883
0884   if (debug) std::cout<<"... done"<<std::endl;
0885
0886   // open event display
0887   display->setOptions("ABDEFHMPT"); // G show geometry
0888   if (matFX) display->setOptions("ABDEFGHMPT"); // G show geometry
0889   display->open();
0890
0891
0892 #endif
0893
0894
0895 }