EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
main.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file main.cc
1 #include <iostream>
2 #include <execinfo.h>
3 #include <signal.h>
4 #include <stdlib.h>
5 
6 #include <AbsFinitePlane.h>
7 #include <AbsFitterInfo.h>
8 #include <AbsMeasurement.h>
9 #include <AbsTrackRep.h>
10 #include <ConstField.h>
11 #include <DetPlane.h>
12 #include <Exception.h>
13 #include <FieldManager.h>
14 #include <KalmanFittedStateOnPlane.h>
15 #include <AbsKalmanFitter.h>
16 #include <KalmanFitter.h>
17 #include <KalmanFitterRefTrack.h>
18 #include <KalmanFitterInfo.h>
19 #include <KalmanFitStatus.h>
20 #include <DAF.h>
21 #include <GFGbl.h>
22 #include <MeasuredStateOnPlane.h>
23 #include <MeasurementOnPlane.h>
24 #include <FullMeasurement.h>
25 #include <PlanarMeasurement.h>
26 #include <ProlateSpacepointMeasurement.h>
27 #include <RectangularFinitePlane.h>
28 #include <ReferenceStateOnPlane.h>
29 #include <SharedPlanePtr.h>
30 #include <SpacepointMeasurement.h>
31 #include <StateOnPlane.h>
32 #include <Tools.h>
33 #include <TrackCand.h>
34 #include <TrackCandHit.h>
35 #include <Track.h>
36 #include <TrackPoint.h>
37 #include <WireMeasurement.h>
38 #include <WirePointMeasurement.h>
39 
40 #include <MaterialEffects.h>
41 #include <RKTools.h>
42 #include <RKTrackRep.h>
43 #include <StepLimits.h>
44 #include <TGeoMaterialInterface.h>
45 
46 #include <EventDisplay.h>
47 
48 #include <HelixTrackModel.h>
49 #include <MeasurementCreator.h>
50 
51 #include <TApplication.h>
52 #include <TCanvas.h>
53 #include <TDatabasePDG.h>
54 #include <TEveManager.h>
55 #include <TGeoManager.h>
56 #include <TH1D.h>
57 #include <TRandom.h>
58 #include <TStyle.h>
59 #include <TVector3.h>
60 #include <vector>
61 
62 #include <TROOT.h>
63 #include <TFile.h>
64 #include <TTree.h>
65 #include "TDatabasePDG.h"
66 #include <TMath.h>
67 #include <TString.h>
68 
69 #include <memory>
70 
71 
72 void handler(int sig) {
73  void *array[10];
74  size_t size;
75 
76  // get void*'s for all entries on the stack
77  size = backtrace(array, 10);
78 
79  // print out all the frames to stderr
80  fprintf(stderr, "Error: signal %d:\n", sig);
81  backtrace_symbols_fd(array, size, 2);
82  exit(1);
83 }
84 
85 int randomPdg() {
86  int pdg;
87 
88  switch(int(gRandom->Uniform(8))) {
89  case 1:
90  pdg = -11; break;
91  case 2:
92  pdg = 11; break;
93  case 3:
94  pdg = 13; break;
95  case 4:
96  pdg = -13; break;
97  case 5:
98  pdg = 211; break;
99  case 6:
100  pdg = -211; break;
101  case 7:
102  pdg = 2212; break;
103  default:
104  pdg = 211;
105  }
106 
107  return pdg;
108 }
109 
110 
111 int randomSign() {
112  if (gRandom->Uniform(1) > 0.5)
113  return 1;
114  return -1;
115 }
116 
117 //---------------------------------------------------------------------------------------------------------
118 //---------------------------------------------------------------------------------------------------------
119 //---------------------------------------------------------------------------------------------------------
120 //---------------------------------------------------------------------------------------------------------
121 //---------------------------------------------------------------------------------------------------------
122 
123 //#define VALGRIND
124 
125 #ifdef VALGRIND
126  #include <valgrind/callgrind.h>
127 #else
128 #define CALLGRIND_START_INSTRUMENTATION
129 #define CALLGRIND_STOP_INSTRUMENTATION
130 #define CALLGRIND_DUMP_STATS
131 #endif
132 
133 int main() {
134  std::cout<<"main"<<std::endl;
135  gRandom->SetSeed(14);
136 
137 
138  const unsigned int nEvents = 1000;
139  const unsigned int nMeasurements = 10;
140  const double BField = 20.; // kGauss
141  const double momentum = 0.1; // GeV
142  const double theta = 110; // degree
143  const double thetaDetPlane = 90; // degree
144  const double phiDetPlane = 0; // degree
145  const double pointDist = 3.; // cm; approx. distance between measurements
146  const double resolution = 0.05; // cm; resolution of generated measurements
147 
148  const double resolutionWire = 5*resolution; // cm; resolution of generated measurements
149  const TVector3 wireDir(0,0,1);
150  const double skewAngle(5);
151  const bool useSkew = true;
152  const int nSuperLayer = 10;
153  const double minDrift = 0.;
154  const double maxDrift = 2;
155  const bool idealLRResolution = false; // resolve the l/r ambiguities of the wire measurements
156 
157  const double outlierProb = -0.1;
158  const double outlierRange = 2;
159 
160  const double hitSwitchProb = -0.1; // probability to give hits to fit in wrong order (flip two hits)
161 
162  const int splitTrack = -5; //nMeasurements/2; // for track merging testing.
163  const bool fullMeasurement = false; // put fit result of first tracklet as FullMeasurement into second tracklet, don't merge
164 
165  //const genfit::eFitterType fitterId = genfit::SimpleKalman;
166  const genfit::eFitterType fitterId = genfit::RefKalman;
167  //const genfit::eFitterType fitterId = genfit::DafRef;
168  //const genfit::eFitterType fitterId = genfit::DafSimple;
169  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedAverage;
170  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToReference;
171  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToPrediction;
172  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToReference;
173  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToPrediction;
174  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToReferenceWire;
175  const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToPredictionWire;
176  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToReferenceWire;
177  //const genfit::eMultipleMeasurementHandling mmHandling = genfit::weightedClosestToPredictionWire;
178  const int nIter = 20; // max number of iterations
179  const double dPVal = 1.E-3; // convergence criterion
180 
181  const bool resort = false;
182  const bool prefit = false; // make a simple Kalman iteration before the actual fit
183  const bool refit = false; // if fit did not converge, try to fit again
184 
185  const bool twoReps = false; // test if everything works with more than one rep in the tracks
186 
187  const bool checkPruning = true; // test pruning
188 
189  const int pdg = 13; // particle pdg code
190 
191  const bool smearPosMom = true; // init the Reps with smeared pos and mom
192  const double chargeSwitchProb = -0.1; // probability to seed with wrong charge sign
193  const double posSmear = 10*resolution; // cm
194  const double momSmear = 5. /180.*TMath::Pi(); // rad
195  const double momMagSmear = 0.2; // relative
196  const double zSmearFac = 2;
197 
198 
199  const bool matFX = false; // include material effects; can only be disabled for RKTrackRep!
200 
201  const bool debug = false;
202  const bool onlyDisplayFailed = false; // only load non-converged tracks into the display
203 
204  std::vector<genfit::eMeasurementType> measurementTypes;
205  for (unsigned int i = 0; i<nMeasurements; ++i) {
206  measurementTypes.push_back(genfit::eMeasurementType(i%8));
207  }
208 
209  signal(SIGSEGV, handler); // install our handler
210 
211  // init fitter
212  genfit::AbsKalmanFitter* fitter = 0;
213  switch (fitterId) {
214  case genfit::SimpleKalman:
215  fitter = new genfit::KalmanFitter(nIter, dPVal);
216  fitter->setMultipleMeasurementHandling(mmHandling);
217  break;
218 
219  case genfit::RefKalman:
220  fitter = new genfit::KalmanFitterRefTrack(nIter, dPVal);
221  fitter->setMultipleMeasurementHandling(mmHandling);
222  break;
223 
224  case genfit::DafSimple:
225  fitter = new genfit::DAF(false);
226  break;
227  case genfit::DafRef:
228  fitter = new genfit::DAF();
229  break;
230  }
231  fitter->setMaxIterations(nIter);
232  //if (debug)
233  // fitter->setDebugLvl(10);
234 
235  /*if (dynamic_cast<genfit::DAF*>(fitter) != nullptr) {
236  //static_cast<genfit::DAF*>(fitter)->setBetas(100, 50, 25, 12, 6, 3, 1, 0.5, 0.1);
237  //static_cast<genfit::DAF*>(fitter)->setBetas(81, 8, 4, 0.5, 0.1);
238  static_cast<genfit::DAF*>(fitter)->setAnnealingScheme(100, 0.1, 5);
239  //static_cast<genfit::DAF*>(fitter)->setConvergenceDeltaWeight(0.0001);
240  //fitter->setMaxIterations(nIter);
241  }*/
242 
243 
244  // init MeasurementCreator
245  genfit::MeasurementCreator measurementCreator;
246  measurementCreator.setResolution(resolution);
247  measurementCreator.setResolutionWire(resolutionWire);
248  measurementCreator.setOutlierProb(outlierProb);
249  measurementCreator.setOutlierRange(outlierRange);
250  measurementCreator.setThetaDetPlane(thetaDetPlane);
251  measurementCreator.setPhiDetPlane(phiDetPlane);
252  measurementCreator.setWireDir(wireDir);
253  measurementCreator.setMinDrift(minDrift);
254  measurementCreator.setMaxDrift(maxDrift);
255  measurementCreator.setIdealLRResolution(idealLRResolution);
256  measurementCreator.setUseSkew(useSkew);
257  measurementCreator.setSkewAngle(skewAngle);
258  measurementCreator.setNSuperLayer(nSuperLayer);
259  measurementCreator.setDebug(debug);
260 
261 
262  // init geometry and mag. field
263  new TGeoManager("Geometry", "Geane geometry");
264  TGeoManager::Import("genfitGeom.root");
265  genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0.,BField));
266  genfit::FieldManager::getInstance()->useCache(true, 8);
267  genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
268 
269  const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
270 
271  // init event display
272 #ifndef VALGRIND
273  genfit::EventDisplay* display = genfit::EventDisplay::getInstance();
274  display->reset();
275 #endif
276 
277 
278 #ifndef VALGRIND
279  // create histograms
280  gROOT->SetStyle("Plain");
281  gStyle->SetPalette(1);
282  gStyle->SetOptFit(1111);
283 
284  TH1D *hmomRes = new TH1D("hmomRes","mom res",500,-20*resolution*momentum/nMeasurements,20*resolution*momentum/nMeasurements);
285  TH1D *hupRes = new TH1D("hupRes","u' res",500,-15*resolution/nMeasurements, 15*resolution/nMeasurements);
286  TH1D *hvpRes = new TH1D("hvpRes","v' res",500,-15*resolution/nMeasurements, 15*resolution/nMeasurements);
287  TH1D *huRes = new TH1D("huRes","u res",500,-15*resolution, 15*resolution);
288  TH1D *hvRes = new TH1D("hvRes","v res",500,-15*resolution, 15*resolution);
289 
290  TH1D *hqopPu = new TH1D("hqopPu","q/p pull",200,-6.,6.);
291  TH1D *pVal = new TH1D("pVal","p-value",100,0.,1.00000001);
292  pVal->SetMinimum(0);
293  TH1D *hupPu = new TH1D("hupPu","u' pull",200,-6.,6.);
294  TH1D *hvpPu = new TH1D("hvpPu","v' pull",200,-6.,6.);
295  TH1D *huPu = new TH1D("huPu","u pull",200,-6.,6.);
296  TH1D *hvPu = new TH1D("hvPu","v pull",200,-6.,6.);
297 
298  TH1D *weights = new TH1D("weights","Daf vs true weights",500,-1.01,1.01);
299 
300  TH1D *trackLenRes = new TH1D("trackLenRes","(trueLen - FittedLen) / trueLen",500,-0.01,0.01);
301 #endif
302 
303  double maxWeight(0);
304  unsigned int nTotalIterConverged(0);
305  unsigned int nTotalIterNotConverged(0);
306  unsigned int nTotalIterSecondConverged(0);
307  unsigned int nTotalIterSecondNotConverged(0);
308  unsigned int nConvergedFits(0);
309  unsigned int nUNConvergedFits(0);
310  unsigned int nConvergedFitsSecond(0);
311  unsigned int nUNConvergedFitsSecond(0);
312 
313 
314  CALLGRIND_START_INSTRUMENTATION;
315 
316  genfit::Track* fitTrack(nullptr);
317  genfit::Track* secondTrack(nullptr);
318 
319  // main loop
320  for (unsigned int iEvent=0; iEvent<nEvents; ++iEvent){
321 
322  /*measurementTypes.clear();
323  for (unsigned int i = 0; i < nMeasurements; ++i)
324  measurementTypes.push_back(genfit::eMeasurementType(gRandom->Uniform(8)));*/
325 
326 
327  if (debug || (iEvent)%10==0)
328  std::cout << "Event Nr. " << iEvent << " ";
329  else std::cout << ". ";
330  if (debug || (iEvent+1)%10==0)
331  std::cout << "\n";
332 
333 
334  // clean up
335  delete fitTrack;
336  fitTrack = nullptr;
337  delete secondTrack;
338  secondTrack = nullptr;
339 
340 
341  // true start values
342  TVector3 pos(0, 0, 0);
343  TVector3 mom(1.,0,0);
344  mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
345  //mom.SetTheta(gRandom->Uniform(0.5*TMath::Pi(),0.9*TMath::Pi()));
346  mom.SetTheta(theta*TMath::Pi()/180);
347  mom.SetMag(momentum);
348  TMatrixDSym covM(6);
349  for (int i = 0; i < 3; ++i)
350  covM(i,i) = resolution*resolution;
351  for (int i = 3; i < 6; ++i)
352  covM(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
353 
354  if (debug) {
355  std::cout << "start values \n";
356  pos.Print();
357  mom.Print();
358  }
359 
360  // calc helix parameters
361  genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
362  measurementCreator.setTrackModel(helix);
363 
364  // smeared start values
365  TVector3 posM(pos);
366  TVector3 momM(mom);
367  if (smearPosMom) {
368  posM.SetX(gRandom->Gaus(posM.X(),posSmear));
369  posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
370  posM.SetZ(gRandom->Gaus(posM.Z(),zSmearFac*posSmear));
371 
372  momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
373  momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
374  momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
375  }
376 
377 
378  // trackrep for creating measurements
379  double sign(1.);
380  if (chargeSwitchProb > gRandom->Uniform(1.))
381  sign = -1.;
382  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(sign*pdg);
383  sign = 1.;
384  if (chargeSwitchProb > gRandom->Uniform(1.))
385  sign = -1.;
386  genfit::AbsTrackRep* secondRep = new genfit::RKTrackRep(sign*-211);
387  genfit::MeasuredStateOnPlane stateRef(rep);
388  rep->setPosMomCov(stateRef, pos, mom, covM);
389 
390  // smeared start state
391  genfit::MeasuredStateOnPlane stateSmeared(rep);
392  rep->setPosMomCov(stateSmeared, posM, momM, covM);
393 
394  //rep->setPropDir(1);
395 
396  if (!matFX) genfit::MaterialEffects::getInstance()->setNoEffects();
397 
398  // remember original initial state
399  const genfit::StateOnPlane stateRefOrig(stateRef);
400 
401  // create smeared measurements
402  std::vector< std::vector<genfit::AbsMeasurement*> > measurements;
403 
404  std::vector<bool> outlierTrue;
405  bool outlier;
406  // true values for left right. 0 for non wire measurements
407  std::vector<int> leftRightTrue;
408  int lr;
409 
410  double trueLen(-1);
411 
412  try{
413  for (unsigned int i=0; i<measurementTypes.size(); ++i){
414  trueLen = i*pointDist;
415 
416  measurements.push_back(measurementCreator.create(measurementTypes[i], trueLen, outlier, lr));
417  outlierTrue.push_back(outlier);
418  leftRightTrue.push_back(lr);
419  }
420  assert(measurementTypes.size() == leftRightTrue.size());
421  assert(measurementTypes.size() == outlierTrue.size());
422  }
423  catch(genfit::Exception& e){
424  std::cerr<<"Exception, next track"<<std::endl;
425  std::cerr << e.what();
426  continue; // here is a memleak!
427  }
428 
429  if (debug) std::cout << "... done creating measurements \n";
430 
431 
432 
433  // create track
434  TVectorD seedState(6);
435  TMatrixDSym seedCov(6);
436  rep->get6DStateCov(stateSmeared, seedState, seedCov);
437  fitTrack = new genfit::Track(rep, seedState, seedCov); //initialized with smeared rep
438  secondTrack = new genfit::Track(rep->clone(), seedState, seedCov); //initialized with smeared rep
439  if (twoReps) {
440  fitTrack->addTrackRep(secondRep);
441  secondTrack->addTrackRep(secondRep->clone());
442  }
443  else
444  delete secondRep;
445  //if (debug) fitTrack->Print("C");
446 
447  fitTrack->checkConsistency();
448  //fitTrack->addTrackRep(rep->clone()); // check if everything works fine with more than one rep
449 
450  // add measurements
451  for(unsigned int i=0; i<measurements.size(); ++i){
452  if (splitTrack > 0 && (int)i >= splitTrack)
453  break;
454  if (i>0 && hitSwitchProb > gRandom->Uniform(1.))
455  fitTrack->insertPoint(new genfit::TrackPoint(measurements[i], fitTrack), -2);
456  else
457  fitTrack->insertPoint(new genfit::TrackPoint(measurements[i], fitTrack));
458 
459  fitTrack->checkConsistency();
460  //if (debug) fitTrack->Print("C");
461  }
462 
463  if (splitTrack > 0) {
464  for(unsigned int i=splitTrack; i<measurements.size(); ++i){
465  if (i>0 && hitSwitchProb > gRandom->Uniform(1.))
466  secondTrack->insertPoint(new genfit::TrackPoint(measurements[i], secondTrack), -2);
467  else
468  secondTrack->insertPoint(new genfit::TrackPoint(measurements[i], secondTrack));
469 
470  //if (debug) fitTrack->Print("C");
471  }
472  }
473 
474  fitTrack->checkConsistency();
475  secondTrack->checkConsistency();
476 
477  //if (debug) fitTrack->Print();
478 
479  // do the fit
480  try{
481  if (debug) std::cout<<"Starting the fitter"<<std::endl;
482 
483  if (prefit) {
484  genfit::KalmanFitter prefitter(1, dPVal);
485  prefitter.setMultipleMeasurementHandling(genfit::weightedClosestToPrediction);
486  prefitter.processTrackWithRep(fitTrack, fitTrack->getCardinalRep());
487  }
488 
489  fitter->processTrack(fitTrack, resort);
490  if (splitTrack > 0)
491  fitter->processTrack(secondTrack, resort);
492 
493  if (debug) std::cout<<"fitter is finished!"<<std::endl;
494  }
495  catch(genfit::Exception& e){
496  std::cerr << e.what();
497  std::cerr << "Exception, next track" << std::endl;
498  continue;
499  }
500 
501  if (splitTrack > 0) {
502  if (debug) fitTrack->Print("C");
503  if (debug) secondTrack->Print("C");
504 
505  if (fullMeasurement) {
506  genfit::FullMeasurement* fullM = new genfit::FullMeasurement(secondTrack->getFittedState());
507  fitTrack->insertPoint(new genfit::TrackPoint(fullM, fitTrack));
508  }
509  else
510  fitTrack->mergeTrack(secondTrack);
511 
512  if (debug) fitTrack->Print("C");
513 
514  try{
515  if (debug) std::cout<<"Starting the fitter"<<std::endl;
516  fitter->processTrack(fitTrack, resort);
517  if (debug) std::cout<<"fitter is finished!"<<std::endl;
518  }
519  catch(genfit::Exception& e){
520  std::cerr << e.what();
521  std::cerr << "Exception, next track" << std::endl;
522  continue;
523  }
524  }
525 
526 
527  if (refit && !fitTrack->getFitStatus(rep)->isFitConverged()) {
528  std::cout<<"Trying to fit again "<<std::endl;
529  fitter->processTrack(fitTrack, resort);
530  }
531 
532 
533 
534  if (debug) {
535  fitTrack->Print("C");
536  fitTrack->getFitStatus(rep)->Print();
537  }
538 
539  fitTrack->checkConsistency();
540  secondTrack->checkConsistency();
541 
542 #ifndef VALGRIND
543  if (!onlyDisplayFailed && iEvent < 1000) {
544  std::vector<genfit::Track*> event;
545  event.push_back(fitTrack);
546  if (splitTrack > 0)
547  event.push_back(secondTrack);
548  display->addEvent(event);
549  }
550  else if (onlyDisplayFailed &&
551  (!fitTrack->getFitStatus(rep)->isFitConverged() ||
552  fitTrack->getFitStatus(rep)->getPVal() < 0.01)) {
553  // add track to event display
554  display->addEvent(fitTrack);
555  }
556 #endif
557 
558 
559  if (fitTrack->getFitStatus(rep)->isFitConverged()) {
560  nTotalIterConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(rep))->getNumIterations();
561  nConvergedFits += 1;
562  }
563  else {
564  nTotalIterNotConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(rep))->getNumIterations();
565  nUNConvergedFits += 1;
566  }
567 
568  if (twoReps) {
569  if (fitTrack->getFitStatus(secondRep)->isFitConverged()) {
570  nTotalIterSecondConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(secondRep))->getNumIterations();
571  nConvergedFitsSecond += 1;
572  }
573  else {
574  nTotalIterSecondNotConverged += static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(secondRep))->getNumIterations();
575  nUNConvergedFitsSecond += 1;
576  }
577  }
578 
579 
580  // check if fit was successful
581  if (! fitTrack->getFitStatus(rep)->isFitConverged()) {
582  std::cout << "Track could not be fitted successfully! Fit is not converged! \n";
583  continue;
584  }
585 
586 
587  genfit::TrackPoint* tp = fitTrack->getPointWithMeasurementAndFitterInfo(0, rep);
588  if (tp == nullptr) {
589  std::cout << "Track has no TrackPoint with fitterInfo! \n";
590  continue;
591  }
592  genfit::KalmanFittedStateOnPlane kfsop(*(static_cast<genfit::KalmanFitterInfo*>(tp->getFitterInfo(rep))->getBackwardUpdate()));
593  if (debug) {
594  std::cout << "state before extrapolating back to reference plane \n";
595  kfsop.Print();
596  }
597 
598  // extrapolate back to reference plane.
599  try{
600  rep->extrapolateToPlane(kfsop, stateRefOrig.getPlane());;
601  }
602  catch(genfit::Exception& e){
603  std::cerr<<"Exception, next track"<<std::endl;
604  std::cerr << e.what();
605  continue;
606  }
607 
608 #ifndef VALGRIND
609  // calculate pulls
610  const TVectorD& referenceState = stateRefOrig.getState();
611 
612  const TVectorD& state = kfsop.getState();
613  const TMatrixDSym& cov = kfsop.getCov();
614 
615  double pval = fitter->getPVal(fitTrack, rep);
616  //assert( fabs(pval - static_cast<genfit::KalmanFitStatus*>(fitTrack->getFitStatus(rep))->getBackwardPVal()) < 1E-10 );
617 
618  hmomRes->Fill( (charge/state[0]-momentum));
619  hupRes->Fill( (state[1]-referenceState[1]));
620  hvpRes->Fill( (state[2]-referenceState[2]));
621  huRes->Fill( (state[3]-referenceState[3]));
622  hvRes->Fill( (state[4]-referenceState[4]));
623 
624  hqopPu->Fill( (state[0]-referenceState[0]) / sqrt(cov[0][0]) );
625  pVal->Fill( pval);
626  hupPu->Fill( (state[1]-referenceState[1]) / sqrt(cov[1][1]) );
627  hvpPu->Fill( (state[2]-referenceState[2]) / sqrt(cov[2][2]) );
628  huPu->Fill( (state[3]-referenceState[3]) / sqrt(cov[3][3]) );
629  hvPu->Fill( (state[4]-referenceState[4]) / sqrt(cov[4][4]) );
630 
631 
632  try {
633  trackLenRes->Fill( (trueLen - fitTrack->getTrackLen(rep)) / trueLen );
634 
635  if (debug) {
636  std::cout << "true track length = " << trueLen << "; fitted length = " << fitTrack->getTrackLen(rep) << "\n";
637  std::cout << "fitted tof = " << fitTrack->getTOF(rep) << " ns\n";
638  }
639  }
640  catch (genfit::Exception& e) {
641  std::cerr << e.what();
642  std::cout << "could not get TrackLen or TOF! \n";
643  }
644 
645 
646 
647  // check l/r resolution and outlier rejection
648  if (dynamic_cast<genfit::DAF*>(fitter) != nullptr) {
649  for (unsigned int i=0; i<fitTrack->getNumPointsWithMeasurement(); ++i){
650 
651  if (! fitTrack->getPointWithMeasurement(i)->hasFitterInfo(rep))
652  continue;
653 
654  if (debug) {
655  std::vector<double> dafWeights = dynamic_cast<genfit::KalmanFitterInfo*>(fitTrack->getPointWithMeasurement(i)->getFitterInfo(rep))->getWeights();
656  std::cout << "hit " << i;
657  if (outlierTrue[i]) std::cout << " is an OUTLIER";
658  std::cout << " weights: ";
659  for (unsigned int j=0; j<dafWeights.size(); ++j){
660  std::cout << dafWeights[j] << " ";
661  }
662  std::cout << " l/r: " << leftRightTrue[i];
663  std::cout << "\n";
664  }
665  int trueSide = leftRightTrue[i];
666  if (trueSide == 0) continue; // not a wire measurement
667  if (outlierTrue[i]) continue; // an outlier
668  std::vector<double> dafWeightLR = dynamic_cast<genfit::KalmanFitterInfo*>(fitTrack->getPointWithMeasurement(i)->getFitterInfo(rep))->getWeights();
669  if(dafWeightLR.size() != 2)
670  continue;
671 
672  double weightCorrectSide, weightWrongSide;
673 
674  if (trueSide < 0) {
675  weightCorrectSide = dafWeightLR[0];
676  weightWrongSide = dafWeightLR[1];
677  }
678  else {
679  weightCorrectSide = dafWeightLR[1];
680  weightWrongSide = dafWeightLR[0];
681  }
682  weightWrongSide -= 1.;
683 
684  weights->Fill(weightCorrectSide);
685  weights->Fill(weightWrongSide);
686 
687  if (weightCorrectSide>maxWeight) maxWeight = weightCorrectSide;
688  }
689 
690  for (unsigned int i=0; i<fitTrack->getNumPointsWithMeasurement(); ++i){
691  if (! fitTrack->getPointWithMeasurement(i)->hasFitterInfo(rep))
692  continue;
693 
694  std::vector<double> dafWeights = dynamic_cast<genfit::KalmanFitterInfo*>(fitTrack->getPointWithMeasurement(i)->getFitterInfo(rep))->getWeights();
695 
696  if (outlierTrue[i]) { // an outlier
697  for (unsigned int j=0; j<dafWeights.size(); ++j){
698  weights->Fill(dafWeights[j]-1);
699  }
700  }
701  else if (leftRightTrue[i] == 0) { // only for non wire hits
702  for (unsigned int j=0; j<dafWeights.size(); ++j){
703  weights->Fill(dafWeights[j]);
704  }
705  }
706  }
707 
708  }
709 
710  if (checkPruning) { //check pruning
711  //std::cout<<"\n";
712  //std::cout<<"get stFirst ";
713  genfit::MeasuredStateOnPlane stFirst = fitTrack->getFittedState();
714  //std::cout<<"get stLast ";
715  genfit::MeasuredStateOnPlane stLast = fitTrack->getFittedState(-1);
716 
717  for (unsigned int i=0; i<1; ++i) {
718  genfit::Track trClone(*fitTrack);
719  trClone.checkConsistency();
720 
721  bool first(false), last(false);
722 
723  TString opt("");
724  try {
725  if (gRandom->Uniform() < 0.5) trClone.prune("C");
726  if (gRandom->Uniform() < 0.5) {
727  opt.Append("F");
728  first = true;
729  }
730  if (gRandom->Uniform() < 0.5) {
731  opt.Append("L");
732  last = true;
733  }
734  if (gRandom->Uniform() < 0.5) opt.Append("W");
735  if (gRandom->Uniform() < 0.5) opt.Append("R");
736  if (gRandom->Uniform() < 0.5) opt.Append("M");
737  if (gRandom->Uniform() < 0.5) opt.Append("I");
738  if (gRandom->Uniform() < 0.5) opt.Append("U");
739 
740  trClone.prune(opt);
741 
742  try {
743  trClone.checkConsistency();
744  } catch (genfit::Exception& e) {
745  trClone.getFitStatus()->getPruneFlags().Print();
746  }
747 
748  //std::cout<<"get stCloneFirst ";
749  genfit::MeasuredStateOnPlane stCloneFirst = trClone.getFittedState();
750  //std::cout<<"get stCloneLast ";
751  genfit::MeasuredStateOnPlane stCloneLast = trClone.getFittedState(-1);
752 
753  if (first and ! (stFirst.getState() == stCloneFirst.getState() and stFirst.getCov() == stCloneFirst.getCov() )) {
754  //std::cout<<" failed first state ";
755  //stFirst.Print();
756  //stCloneFirst.Print();
757 
758  if (debug)
759  trClone.getFitStatus()->getPruneFlags().Print();
760  }
761 
762  if (last and ! (stLast.getState() == stCloneLast.getState() and stLast.getCov() == stCloneLast.getCov() )) {
763  //std::cout<<" failed last state ";
764  //stLast.Print();
765  //stCloneLast.Print();
766 
767  if (debug)
768  trClone.getFitStatus()->getPruneFlags().Print();
769  }
770 
771  if (debug) {
772  std::cout<<" pruned track: ";
773  trClone.Print();
774  }
775  }
776  catch (genfit::Exception &e) {
777  std::cerr << e.what();
778  }
779  }
780 
781  } // end check pruning
782 
783 #endif
784 
785 
786  }// end loop over events
787 
788  delete fitTrack;
789  delete secondTrack;
790  delete fitter;
791 
792  CALLGRIND_STOP_INSTRUMENTATION;
793  CALLGRIND_DUMP_STATS;
794 
795  std::cout<<"maxWeight = " << maxWeight << std::endl;
796  std::cout<<"avg nr iterations = " << (double)(nTotalIterConverged + nTotalIterNotConverged)/(double)(nConvergedFits + nUNConvergedFits) << std::endl;
797  std::cout<<"avg nr iterations of converged fits = " << (double)(nTotalIterConverged)/(double)(nConvergedFits) << std::endl;
798  std::cout<<"avg nr iterations of UNconverged fits = " << (double)(nTotalIterNotConverged)/(double)(nUNConvergedFits) << std::endl;
799  std::cout<<"fit efficiency = " << (double)nConvergedFits/nEvents << std::endl;
800 
801  if (twoReps) {
802  std::cout<<"second rep: \navg nr iterations = " << (double)(nTotalIterSecondConverged + nTotalIterSecondNotConverged)/(double)(nConvergedFitsSecond + nUNConvergedFitsSecond) << std::endl;
803  std::cout<<"avg nr iterations of converged fits = " << (double)(nTotalIterSecondConverged)/(double)(nConvergedFitsSecond) << std::endl;
804  std::cout<<"avg nr iterations of UNconverged fits = " << (double)(nTotalIterSecondNotConverged)/(double)(nUNConvergedFitsSecond) << std::endl;
805  std::cout<<"fit efficiency = " << (double)nConvergedFitsSecond/nEvents << std::endl;
806  }
807 
808 
809  //std::cout<<"avg nr iterations (2nd rep) = " << (double)nTotalIterSecond/nSuccessfullFitsSecond << std::endl;
810  //std::cout<<"fit efficiency (2nd rep) = " << (double)nConvergedFitsSecond/nEvents << std::endl;
811 
812 
813 #ifndef VALGRIND
814 
815  if (debug) std::cout<<"Draw histograms ...";
816  // fit and draw histograms
817  TCanvas* c1 = new TCanvas();
818  c1->Divide(2,3);
819 
820  c1->cd(1);
821  hmomRes->Fit("gaus");
822  hmomRes->Draw();
823 
824  c1->cd(2);
825  weights->Draw();
826 
827  c1->cd(3);
828  hupRes->Fit("gaus");
829  hupRes->Draw();
830 
831  c1->cd(4);
832  hvpRes->Fit("gaus");
833  hvpRes->Draw();
834 
835  c1->cd(5);
836  huRes->Fit("gaus");
837  huRes->Draw();
838 
839  c1->cd(6);
840  hvRes->Fit("gaus");
841  hvRes->Draw();
842 
843  c1->Write();
844 
845  TCanvas* c2 = new TCanvas();
846  c2->Divide(2,3);
847 
848  c2->cd(1);
849  hqopPu->Fit("gaus");
850  hqopPu->Draw();
851 
852  c2->cd(2);
853  pVal->Fit("pol1");
854  pVal->Draw();
855  c2->cd(3);
856  hupPu->Fit("gaus");
857  hupPu->Draw();
858 
859  c2->cd(4);
860  hvpPu->Fit("gaus");
861  hvpPu->Draw();
862 
863  c2->cd(5);
864  huPu->Fit("gaus");
865  huPu->Draw();
866 
867  c2->cd(6);
868  hvPu->Fit("gaus");
869  hvPu->Draw();
870 
871  c2->Write();
872 
873 
874 
875  TCanvas* c3 = new TCanvas();
876  //c3->Divide(2,3);
877 
878  c3->cd(1);
879  trackLenRes->Fit("gaus");
880  trackLenRes->Draw();
881 
882  c3->Write();
883 
884  if (debug) std::cout<<"... done"<<std::endl;
885 
886  // open event display
887  display->setOptions("ABDEFHMPT"); // G show geometry
888  if (matFX) display->setOptions("ABDEFGHMPT"); // G show geometry
889  display->open();
890 
891 
892 #endif
893 
894 
895 }