Analysis Software
Documentation for sPHENIX simulation 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>
15 #include <KalmanFitterInfo.h>
16 #include <MeasuredStateOnPlane.h>
17 #include <MeasurementOnPlane.h>
18 #include <PlanarMeasurement.h>
20 #include <RectangularFinitePlane.h>
21 #include <ReferenceStateOnPlane.h>
22 #include <SharedPlanePtr.h>
23 #include <SpacepointMeasurement.h>
24 #include <StateOnPlane.h>
25 #include <Tools.h>
26 #include <TrackCand.h>
27 #include <TrackCandHit.h>
28 #include <Track.h>
29 #include <TrackPoint.h>
30 #include <WireMeasurement.h>
31 #include <WirePointMeasurement.h>
32 
33 #include <MaterialEffects.h>
34 #include <RKTools.h>
35 #include <RKTrackRep.h>
36 #include <StepLimits.h>
37 #include <TGeoMaterialInterface.h>
38 
39 #include <TApplication.h>
40 #include <TCanvas.h>
41 #include <TDatabasePDG.h>
42 #include <TEveManager.h>
43 #include <TGeoManager.h>
44 #include <TH1D.h>
45 #include <TH2D.h>
46 #include <TRandom.h>
47 #include <TStyle.h>
48 #include <TVector3.h>
49 
50 #include <vector>
51 #include <map>
52 #include <stdio.h>
53 #include <stdlib.h>
54 
55 #include <TROOT.h>
56 #include <TFile.h>
57 #include <TTree.h>
58 #include "TDatabasePDG.h"
59 #include <TMath.h>
60 #include <TString.h>
61 
62 //#include <callgrind/callgrind.h>
63 
64 
69 };
70 
71 constexpr bool verbose = false;
72 
73 void handler(int sig) {
74  void *array[10];
75  size_t size;
76 
77  // get void*'s for all entries on the stack
78  size = backtrace(array, 10);
79 
80  // print out all the frames to stderr
81  fprintf(stderr, "Error: signal %d:\n", sig);
82  backtrace_symbols_fd(array, size, 2);
83  exit(1);
84 }
85 
86 int randomPdg() {
87  int pdg;
88 
89  switch(int(gRandom->Uniform(8))) {
90  case 1:
91  pdg = -11; break;
92  case 2:
93  pdg = 11; break;
94  case 3:
95  pdg = 13; break;
96  case 4:
97  pdg = -13; break;
98  case 5:
99  pdg = 211; break;
100  case 6:
101  pdg = -211; break;
102  case 7:
103  pdg = 2212; break;
104  default:
105  pdg = 211;
106  }
107 
108  return pdg;
109 }
110 
111 
112 int randomSign() {
113  if (gRandom->Uniform(1) > 0.5)
114  return 1;
115  return -1;
116 }
117 
118 
119 e_testStatus compareMatrices(const TMatrixTBase<double>& A, const TMatrixTBase<double>& B, double maxRelErr) {
120  e_testStatus retVal = kPassed;
121 
122  // search max abs value
123  double max(0);
124  for (int i=0; i<A.GetNrows(); ++i) {
125  for (int j=0; j<A.GetNcols(); ++j) {
126  if (fabs(A(i,j)) > max)
127  max = fabs(A(i,j));
128  }
129  }
130 
131  double maxAbsErr = maxRelErr*max;
132 
133  for (int i=0; i<A.GetNrows(); ++i) {
134  for (int j=0; j<A.GetNcols(); ++j) {
135  double absErr = A(i,j) - B(i,j);
136  if ( fabs(absErr) > maxAbsErr ) {
137  double relErr = A(i,j)/B(i,j) - 1;
138  if ( fabs(relErr) > maxRelErr ) {
139  if (verbose) {
140  std::cout << "compareMatrices: A("<<i<<","<<j<<") = " << A(i,j) << " B("<<i<<","<<j<<") = " << B(i,j) << " absErr = " << absErr << " relErr = " << relErr << "\n";
141  }
142  retVal = kFailed;
143  }
144  }
145  }
146  }
147  return retVal;
148 }
149 
150 e_testStatus isCovMatrix(TMatrixTBase<double>& cov) {
151 
152  if (!(cov.IsSymmetric())) {
153  if (verbose) {
154  std::cout << "isCovMatrix: not symmetric\n";
155  }
156  return kFailed;
157  }
158 
159  for (int i=0; i<cov.GetNrows(); ++i) {
160  for (int j=0; j<cov.GetNcols(); ++j) {
161  if (std::isnan(cov(i,j))) {
162  if (verbose) {
163  std::cout << "isCovMatrix: element isnan\n";
164  }
165  return kFailed;
166  }
167  if (i==j && cov(i,j) < 0) {
168  if (verbose) {
169  std::cout << "isCovMatrix: negative diagonal element\n";
170  }
171  return kFailed;
172  }
173  }
174  }
175 
176  return kPassed;
177 }
178 
179 
180 
181 e_testStatus checkSetGetPosMom(bool writeHisto = false) {
182 
183  if (writeHisto)
184  return kPassed;
185 
186  double epsilonLen = 1.E-10;
187  double epsilonMom = 1.E-10;
188 
189  int pdg = randomPdg();
190  genfit::AbsTrackRep* rep;
191  rep = new genfit::RKTrackRep(pdg);
192 
193  //TVector3 pos(0,0,0);
194  TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
195  TVector3 mom(0,0.5,gRandom->Gaus(0,0.3));
196  mom.SetMag(0.5);
197  mom *= randomSign();
198 
199 
201  rep->setPosMom(state, pos, mom);
202 
203  // check if we can set another position in the same plane
204  if (randomSign() == 1) {
206  const TVector3& u = plane->getU();
207  const TVector3& v = plane->getV();
208 
209  // random position on plane
210  pos += gRandom->Gaus() * u;
211  pos += gRandom->Gaus() * v;
212 
213  // new random momentum
214  mom.SetXYZ(0,0.5,gRandom->Gaus(0,0.3));
215  mom.SetMag(0.5);
216  mom *= randomSign();
217 
218  rep->setPosMom(state, pos, mom);
219 
220  // check if plane has changed
221  if (state.getPlane() != plane) {
222  if (verbose) {
223  std::cout << "plane has changed unexpectedly! \n";
224  }
225  delete rep;
226  return kFailed;
227  }
228  }
229 
230 
231  // compare
232  if ((pos - rep->getPos(state)).Mag() > epsilonLen ||
233  (mom - rep->getMom(state)).Mag() > epsilonMom) {
234 
235  if (verbose) {
236  state.Print();
237  std::cout << "pos difference = " << (pos - rep->getPos(state)).Mag() << "\n";
238  std::cout << "mom difference = " << (mom - rep->getMom(state)).Mag() << "\n";
239 
240  std::cout << std::endl;
241  }
242  delete rep;
243  return kFailed;
244  }
245 
246  delete rep;
247  return kPassed;
248 
249 }
250 
251 
252 e_testStatus compareForthBackExtrapolation(bool writeHisto = false) {
253  static std::map<int, std::vector<TH2D*> > histoMap;
254 
255  static const int nPDGs(5);
256  int pdgs[nPDGs]={0, 211,13,11,2212};
257  static bool fill(true);
258  if (fill) {
259  for (int i = 0; i<nPDGs; ++i) {
260  int pdg = pdgs[i];
261 
262  std::string Result;//string which will contain the result
263  std::stringstream convert; // stringstream used for the conversion
264  convert << pdg;//add the value of Number to the characters in the stream
265  Result = convert.str();//set Result to the content of the stream
266 
267  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));
268  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));
269  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));
270  }
271  fill = false;
272  }
273 
274 
275  if (writeHisto) {
276  TFile outfile("deviation.root", "recreate");
277  outfile.cd();
278 
279  for (int i = 0; i<nPDGs; ++i) {
280  int pdg = pdgs[i];
281  histoMap[pdg][0]->Write();
282  histoMap[pdg][1]->Write();
283  histoMap[pdg][2]->Write();
284  }
285 
286  outfile.Close();
287 
288  return kPassed;
289  }
290 
291 
292  double epsilonLen = 5.E-5; // 0.5 mu
293  double epsilonMom = 1.E-6; // 1 keV
294 
295  int pdg = randomPdg();
296  //pdg = 211;
297  genfit::AbsTrackRep* rep;
298  rep = new genfit::RKTrackRep(pdg);
299 
300  //rep->setDebugLvl(1);
301  //genfit::MaterialEffects::getInstance()->setDebugLvl(1);
302 
303 
304 
305  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
306  double mass = part->Mass(); // GeV
307 
308  //TVector3 pos(0,0,0);
309  TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
310  TVector3 mom(0,0.5,gRandom->Gaus(0,0.3));
311  mom.SetMag( exp(gRandom->Uniform(-4, 8)) * mass );
312  mom *= randomSign();
313 
314 
315  double betaGamma = log(mom.Mag()/mass);
316 
317  //mom.Print();
318 
320  rep->setPosMom(state, pos, mom);
321 
322  genfit::SharedPlanePtr origPlane = state.getPlane();
323  genfit::SharedPlanePtr plane(new genfit::DetPlane(TVector3(0,randomSign()*10,0), TVector3(0,randomSign()*1,0)));
324 
325  genfit::StateOnPlane origState(state);
326 
327  TVector3 mom2;
328  double momLoss1, momLoss2;
329 
330  // forth
331  double extrapLen(0);
332  try {
333  extrapLen = rep->extrapolateToPlane(state, plane);
334 
335  mom2 = state.getMom();
336  momLoss1 = mom.Mag()-mom2.Mag();
337 
338  //mom2.Print();
339  //std::cout << "MomLoss = " << momLoss1 << "\n";
340  }
341  catch (genfit::Exception& e) {
342  if (verbose) {
343  std::cerr << "Exception in forth Extrapolation. PDG = " << pdg << "; mom: \n";
344  mom.Print();
345 
346  std::cerr << e.what();
347  }
348  delete rep;
349  return kException;
350  }
351 
352 
353 
354  // back
355  double backExtrapLen(0);
356  try {
357  backExtrapLen = rep->extrapolateToPlane(state, origPlane);
358 
359  momLoss2 = mom2.Mag()-state.getMom().Mag();
360 
361  //state.getMom().Print();
362  //std::cout << "MomLoss = " << momLoss2 << "\n";
363 
364  double deviation = 1. + momLoss1/momLoss2;
365  histoMap[abs(pdg)][0]->Fill(deviation, betaGamma);
366  histoMap[abs(pdg)][1]->Fill((mom.Mag() - state.getMom().Mag())*1e6, betaGamma);
367  histoMap[abs(pdg)][2]->Fill(deviation, extrapLen+backExtrapLen);
368 
369  histoMap[0][0]->Fill(deviation, betaGamma);
370  histoMap[0][1]->Fill((mom.Mag() - state.getMom().Mag())*1e6, betaGamma);
371  histoMap[0][2]->Fill(deviation, extrapLen+backExtrapLen);
372 
373  //std::cout << "deviation = " << deviation << "\n";
374  }
375  catch (genfit::Exception& e) {
376  if (verbose) {
377  std::cerr << "Exception in back Extrapolation. PDG = " << pdg << "; mom: \n";
378  mom.Print();
379  std::cerr << "mom2: \n";
380  mom2.Print();
381  }
382  std::cerr << e.what();
383 
384  delete rep;
385  return kException;
386  }
387 
388  // compare
389  if ((rep->getPos(origState) - rep->getPos(state)).Mag() > epsilonLen ||
390  (rep->getMom(origState) - rep->getMom(state)).Mag() > epsilonMom ||
391  fabs(extrapLen + backExtrapLen) > epsilonLen) {
392 
393  if (verbose) {
394  origState.Print();
395  state.Print();
396 
397  std::cout << "pos difference = " << (rep->getPos(origState) - rep->getPos(state)).Mag() << "\n";
398  std::cout << "mom difference = " << (rep->getMom(origState) - rep->getMom(state)).Mag() << "\n";
399  std::cout << "len difference = " << extrapLen + backExtrapLen << "\n";
400 
401  std::cout << std::endl;
402  }
403  delete rep;
404  return kFailed;
405  }
406 
407  delete rep;
408  return kPassed;
409 
410 }
411 
412 
413 e_testStatus compareForthBackJacNoise(bool writeHisto = false) {
414 
415  if (writeHisto)
416  return kPassed;
417 
418  e_testStatus retVal(kPassed);
419 
420  bool fx( randomSign() > 0);
422 
423  double deltaJac = 0.005; // relative
424  double deltaNoise = 0.00001;
425  double deltaState = 3.E-6; // absolute
426 
427  if (fx) {
428  deltaJac = 0.1; // relative
429  deltaNoise = 0.1;
430  deltaState = 5.E-4; // absolute
431  }
432 
433  int pdg = randomPdg();
434  genfit::AbsTrackRep* rep;
435  rep = new genfit::RKTrackRep(pdg);
436 
437  //TVector3 pos(0,0,0);
438  //TVector3 mom(0,1,2);
439  TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
440  TVector3 mom(0, 0.5, gRandom->Gaus(0, 1));
441  mom *= randomSign();
442  mom.SetMag(gRandom->Uniform(2)+0.3);
443  //mom.SetMag(3);
444 
445  TMatrixD jac_f, jac_fi, jac_b, jac_bi;
446  TMatrixDSym noise_f, noise_fi, noise_b, noise_bi;
447  TVectorD c_f, c_fi, c_b, c_bi;
448  TVectorD state_man, stateOrig_man;
449 
450  // original state and plane
452  rep->setPosMom(state, pos, mom);
453 
454  static const double smear = 0.2;
455  TVector3 normal(mom);
456  normal.SetMag(1);
457  normal.SetXYZ(gRandom->Gaus(normal.X(), smear),
458  gRandom->Gaus(normal.Y(), smear),
459  gRandom->Gaus(normal.Z(), smear));
460  genfit::DetPlane* origPlanePtr = new genfit::DetPlane (pos, normal);
461  //genfit::DetPlane* origPlanePtr = new genfit::DetPlane (pos, TVector3(1,0,0), TVector3(0,0,1));
462  double rotAngleOrig = gRandom->Uniform(2.*TMath::Pi());
463  origPlanePtr->rotate(rotAngleOrig);
464  genfit::SharedPlanePtr origPlane(origPlanePtr);
465  //genfit::SharedPlanePtr origPlane = state.getPlane();
466  rep->extrapolateToPlane(state, origPlane);
467 
468  const genfit::StateOnPlane origState(state);
469 
470 
471  // dest plane
472  normal = mom;
473  normal.SetMag(1);
474  normal.SetXYZ(gRandom->Gaus(normal.X(), smear),
475  gRandom->Gaus(normal.Y(), smear),
476  gRandom->Gaus(normal.Z(), smear));
477  TVector3 dest(mom);
478  dest.SetMag(10);
479  genfit::DetPlane* planePtr = new genfit::DetPlane (dest, normal);
480  //genfit::DetPlane* planePtr = new genfit::DetPlane (dest, TVector3(1,0,0), TVector3(0,0,1));
481  double rotAngle = gRandom->Uniform(2.*TMath::Pi());
482  planePtr->rotate(rotAngle);
483  genfit::SharedPlanePtr plane(planePtr);
484 
485  /* genfit::DetPlane* planePtr = new genfit::DetPlane (*origPlane);
486  planePtr->setO(TVector3(0,randomSign()*10,0));
487  //planePtr->rotate(rotAngle);
488  genfit::SharedPlanePtr plane(planePtr);
489 */
490 
491  // numerical calculation
492  TMatrixD jac_f_num;
493  rep->calcJacobianNumerically(origState, plane, jac_f_num);
494 
495 
496  // forth
497  genfit::StateOnPlane extrapolatedState;
498  try {
499  //std::cout << "DO FORTH EXTRAPOLATION \n";
500  rep->extrapolateToPlane(state, plane);
501  //std::cout << "GET INFO FOR FORTH EXTRAPOLATION \n";
502  extrapolatedState = state;
503  rep->getForwardJacobianAndNoise(jac_f, noise_f, c_f);
504  rep->getBackwardJacobianAndNoise(jac_fi, noise_fi, c_fi);
505  }
506  catch (genfit::Exception& e) {
507  std::cerr << e.what();
508 
509  delete rep;
511  return kException;
512  }
513 
514  // back
515  try {
516  //std::cout << "DO BACK EXTRAPOLATION \n";
517  rep->extrapolateToPlane(state, origPlane);
518  //std::cout << "GET INFO FOR BACK EXTRAPOLATION \n";
519  rep->getForwardJacobianAndNoise(jac_b, noise_b, c_b);
520  rep->getBackwardJacobianAndNoise(jac_bi, noise_bi, c_bi);
521  }
522  catch (genfit::Exception& e) {
523  std::cerr << e.what();
524 
525  delete rep;
527  return kException;
528  }
529 
530 
531  // compare
532  if (!((origState.getState() - state.getState()).Abs() < deltaState) ) {
533  if (verbose) {
534  std::cout << "(origState.getState() - state.getState()) ";
535  (origState.getState() - state.getState()).Print();
536  }
537 
538  retVal = kFailed;
539  }
540 
541  // check c
542  if (!((jac_f * origState.getState() + c_f - extrapolatedState.getState()).Abs() < deltaState) ||
543  !((jac_bi * origState.getState() + c_bi - extrapolatedState.getState()).Abs() < deltaState) ||
544  !((jac_b * extrapolatedState.getState() + c_b - origState.getState()).Abs() < deltaState) ||
545  !((jac_fi * extrapolatedState.getState() + c_fi - origState.getState()).Abs() < deltaState) ) {
546 
547  if (verbose) {
548  std::cout << "(jac_f * origState.getState() + c_f - extrapolatedState.getState()) ";
549  (jac_f * origState.getState() + c_f - extrapolatedState.getState()).Print();
550  std::cout << "(jac_bi * origState.getState() + c_bi - extrapolatedState.getState()) ";
551  (jac_bi * origState.getState() + c_bi - extrapolatedState.getState()).Print();
552  std::cout << "(jac_b * extrapolatedState.getState() + c_b - origState.getState()) ";
553  (jac_b * extrapolatedState.getState() + c_b - origState.getState()).Print();
554  std::cout << "(jac_fi * extrapolatedState.getState() + c_fi - origState.getState()) ";
555  (jac_fi * extrapolatedState.getState() + c_fi - origState.getState()).Print();
556  }
557  retVal = kFailed;
558  }
559 
560  if (isCovMatrix(state.getCov()) == kFailed) {
561  retVal = kFailed;
562  }
563 
564 
565  // compare
566  if (compareMatrices(jac_f, jac_bi, deltaJac) == kFailed) {
567  if (verbose) {
568  std::cout << "jac_f = ";
569  jac_f.Print();
570  std::cout << "jac_bi = ";
571  jac_bi.Print();
572  std::cout << std::endl;
573  }
574  retVal = kFailed;
575  }
576 
577  // compare
578  if (compareMatrices(jac_b, jac_fi, deltaJac) == kFailed) {
579  if (verbose) {
580  std::cout << "jac_b = ";
581  jac_b.Print();
582  std::cout << "jac_fi = ";
583  jac_fi.Print();
584  std::cout << std::endl;
585  }
586  retVal = kFailed;
587  }
588 
589  // compare
590  if (compareMatrices(noise_f, noise_bi, deltaNoise) == kFailed) {
591  if (verbose) {
592  std::cout << "noise_f = ";
593  noise_f.Print();
594  std::cout << "noise_bi = ";
595  noise_bi.Print();
596  std::cout << std::endl;
597  }
598  retVal = kFailed;
599  }
600 
601  // compare
602  if (compareMatrices(noise_b, noise_fi, deltaNoise) == kFailed) {
603  if (verbose) {
604  std::cout << "noise_b = ";
605  noise_b.Print();
606  std::cout << "noise_fi = ";
607  noise_fi.Print();
608  std::cout << std::endl;
609  }
610  retVal = kFailed;
611  }
612 
613 
614  if (!fx) {
615  // compare
616  if (compareMatrices(jac_f, jac_f_num, deltaJac) == kFailed) {
617  if (verbose) {
618  std::cout << "jac_f = ";
619  jac_f.Print();
620  std::cout << "jac_f_num = ";
621  jac_f_num.Print();
622  std::cout << std::endl;
623  }
624  retVal = kFailed;
625  }
626  }
627 
628  delete rep;
630 
631  return retVal;
632 }
633 
634 
635 e_testStatus checkStopAtBoundary(bool writeHisto = false) {
636 
637  if (writeHisto)
638  return kPassed;
639 
640  double epsilonLen = 1.E-4; // 1 mu
641  //double epsilonMom = 1.E-5; // 10 keV
642 
643  double matRadius(1.);
644 
645  int pdg = randomPdg();
646  genfit::AbsTrackRep* rep;
647  rep = new genfit::RKTrackRep(pdg);
648 
649  //TVector3 pos(0,0,0);
650  TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
651  TVector3 mom(0,0.5,gRandom->Gaus(0,0.1));
652  mom *= randomSign();
653 
654 
656  rep->setPosMom(state, pos, mom);
657 
658  genfit::SharedPlanePtr origPlane = state.getPlane();
659  genfit::SharedPlanePtr plane(new genfit::DetPlane(TVector3(0,randomSign()*10,0), TVector3(0,randomSign()*1,0)));
660 
661  genfit::StateOnPlane origState(state);
662 
663  // forth
664  try {
665  rep->extrapolateToPlane(state, plane, true);
666  }
667  catch (genfit::Exception& e) {
668  std::cerr << e.what();
669 
670  delete rep;
671  return kException;
672  }
673 
674 
675  // compare
676  if (fabs(rep->getPos(state).Perp() - matRadius) > epsilonLen) {
677  if (verbose) {
678  origState.Print();
679  state.Print();
680 
681  std::cerr << "radius difference = " << rep->getPos(state).Perp() - matRadius << "\n";
682 
683  std::cerr << std::endl;
684  }
685  delete rep;
686  return kFailed;
687  }
688 
689  delete rep;
690  return kPassed;
691 
692 }
693 
694 
695 e_testStatus checkErrorPropagation(bool writeHisto = false) {
696 
697  if (writeHisto)
698  return kPassed;
699 
700  //double epsilonLen = 1.E-4; // 1 mu
701  //double epsilonMom = 1.E-5; // 10 keV
702 
703  int pdg = randomPdg();
704  genfit::AbsTrackRep* rep;
705  rep = new genfit::RKTrackRep(pdg);
706 
707  //TVector3 pos(0,0,0);
708  TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
709  TVector3 mom(0,0.5,gRandom->Gaus(0,0.1));
710  mom *= randomSign();
711 
712 
714  rep->setPosMom(state, pos, mom);
715 
716  genfit::SharedPlanePtr origPlane = state.getPlane();
717  genfit::SharedPlanePtr plane(new genfit::DetPlane(TVector3(0,randomSign()*50,0), TVector3(0,randomSign()*1,0)));
718 
719  genfit::MeasuredStateOnPlane origState(state);
720 
721  // forth
722  try {
723  rep->extrapolateToPlane(state, plane);
724  }
725  catch (genfit::Exception& e) {
726  std::cerr << e.what();
727 
728  delete rep;
729  return kException;
730  }
731 
732 
733  // check
734  if (isCovMatrix(state.getCov()) == kFailed) {
735  if (verbose) {
736  origState.Print();
737  state.Print();
738  }
739  delete rep;
740  return kFailed;
741  }
742 
743  delete rep;
744  return kPassed;
745 
746 }
747 
748 
749 e_testStatus checkExtrapolateToLine(bool writeHisto = false) {
750 
751  if (writeHisto)
752  return kPassed;
753 
754  double epsilonLen = 1.E-4; // 1 mu
755  double epsilonMom = 1.E-5; // 10 keV
756 
757  int pdg = randomPdg();
758  genfit::AbsTrackRep* rep;
759  rep = new genfit::RKTrackRep(pdg);
760 
761  //TVector3 pos(0,0,0);
762  TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
763  TVector3 mom(0,0.5,0.);
764  mom *= randomSign();
765 
766 
768  rep->setPosMom(state, pos, mom);
769 
770  genfit::SharedPlanePtr origPlane = state.getPlane();
771  genfit::StateOnPlane origState(state);
772 
773  TVector3 linePoint(gRandom->Gaus(),randomSign()*10+gRandom->Gaus(),gRandom->Gaus());
774  TVector3 lineDirection(gRandom->Gaus(),gRandom->Gaus(),randomSign()*10+gRandom->Gaus());
775 
776  // forth
777  try {
778  rep->extrapolateToLine(state, linePoint, lineDirection, false);
779  }
780  catch (genfit::Exception& e) {
781  std::cerr << e.what();
782 
783  delete rep;
784  return kException;
785  }
786 
787 
788  // compare
789  if (fabs(state.getPlane()->distance(linePoint)) > epsilonLen ||
790  fabs(state.getPlane()->distance(linePoint+lineDirection)) > epsilonLen ||
791  (rep->getMom(state).Unit() * state.getPlane()->getNormal()) > epsilonMom) {
792 
793  if (verbose) {
794  origState.Print();
795  state.Print();
796 
797  std::cout << "distance of linePoint to plane = " << state.getPlane()->distance(linePoint) << "\n";
798  std::cout << "distance of linePoint+lineDirection to plane = "
799  << state.getPlane()->distance(linePoint + lineDirection) << "\n";
800  std::cout << "direction * plane normal = " << rep->getMom(state).Unit() * state.getPlane()->getNormal() << "\n";
801  }
802  delete rep;
803  return kFailed;
804  }
805 
806  delete rep;
807  return kPassed;
808 
809 }
810 
811 
812 e_testStatus checkExtrapolateToPoint(bool writeHisto = false) {
813 
814  if (writeHisto)
815  return kPassed;
816 
817  double epsilonLen = 1.E-4; // 1 mu
818  double epsilonMom = 1.E-5; // 10 keV
819 
820  int pdg = randomPdg();
821  genfit::AbsTrackRep* rep;
822  rep = new genfit::RKTrackRep(pdg);
823 
824  //TVector3 pos(0,0,0);
825  TVector3 pos(gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1),gRandom->Gaus(0,0.1));
826  TVector3 mom(0,0.5,gRandom->Gaus(0,0.1));
827  mom *= randomSign();
828 
829 
831  rep->setPosMom(state, pos, mom);
832 
833  genfit::SharedPlanePtr origPlane = state.getPlane();
834  genfit::StateOnPlane origState(state);
835 
836  TVector3 point(gRandom->Gaus(),randomSign()*10+gRandom->Gaus(),gRandom->Gaus());
837 
838  // forth
839  try {
840  rep->extrapolateToPoint(state, point, false);
841  }
842  catch (genfit::Exception& e) {
843  std::cerr << e.what();
844 
845  delete rep;
846  return kException;
847  }
848 
849 
850  // compare
851  if (fabs(state.getPlane()->distance(point)) > epsilonLen ||
852  fabs((rep->getMom(state).Unit() * state.getPlane()->getNormal())) - 1 > epsilonMom) {
853  if (verbose) {
854  origState.Print();
855  state.Print();
856 
857  std::cout << "distance of point to plane = " << state.getPlane()->distance(point) << "\n";
858  std::cout << "direction * plane normal = " << rep->getMom(state).Unit() * state.getPlane()->getNormal() << "\n";
859  }
860  delete rep;
861  return kFailed;
862  }
863 
864  delete rep;
865  return kPassed;
866 
867 }
868 
869 
870 e_testStatus checkExtrapolateToCylinder(bool writeHisto = false) {
871 
872  if (writeHisto)
873  return kPassed;
874 
875  double epsilonLen = 1.E-4; // 1 mu
876  //double epsilonMom = 1.E-5; // 10 keV
877 
878  int pdg = randomPdg();
879  genfit::AbsTrackRep* rep;
880  rep = new genfit::RKTrackRep(pdg);
881 
882  //TVector3 pos(0,0,0);
883  TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
884  TVector3 mom(0,0.5,0.);
885  mom *= randomSign();
886 
887 
889  rep->setPosMom(state, pos, mom);
890 
891  genfit::SharedPlanePtr origPlane = state.getPlane();
892  genfit::StateOnPlane origState(state);
893 
894  const TVector3 linePoint(gRandom->Gaus(0,5), gRandom->Gaus(0,5), gRandom->Gaus(0,5));
895  const TVector3 lineDirection(gRandom->Gaus(),gRandom->Gaus(),2+gRandom->Gaus());
896  const double radius = gRandom->Uniform(10);
897 
898  // forth
899  try {
900  rep->extrapolateToCylinder(state, radius, linePoint, lineDirection, false);
901  }
902  catch (genfit::Exception& e) {
903  std::cerr << e.what();
904 
905  delete rep;
906 
907  return kException;
908  }
909 
910  TVector3 pocaOnLine(lineDirection);
911  double t = 1./(pocaOnLine.Mag2()) * ((rep->getPos(state)*pocaOnLine) - (linePoint*pocaOnLine));
912  pocaOnLine *= t;
913  pocaOnLine += linePoint;
914 
915  TVector3 radiusVec = rep->getPos(state) - pocaOnLine;
916 
917  // compare
918  if (fabs(state.getPlane()->getNormal()*radiusVec.Unit())-1 > epsilonLen ||
919  fabs(lineDirection*radiusVec) > epsilonLen ||
920  fabs(radiusVec.Mag()-radius) > epsilonLen) {
921  if (verbose) {
922  origState.Print();
923  state.Print();
924 
925  std::cout << "lineDirection*radiusVec = " << lineDirection * radiusVec << "\n";
926  std::cout << "radiusVec.Mag()-radius = " << radiusVec.Mag() - radius << "\n";
927  }
928  delete rep;
929  return kFailed;
930  }
931 
932  delete rep;
933  return kPassed;
934 
935 }
936 
937 
938 e_testStatus checkExtrapolateToSphere(bool writeHisto = false) {
939 
940  if (writeHisto)
941  return kPassed;
942 
943  double epsilonLen = 1.E-4; // 1 mu
944  //double epsilonMom = 1.E-5; // 10 keV
945 
946  int pdg = randomPdg();
947  genfit::AbsTrackRep* rep;
948  rep = new genfit::RKTrackRep(pdg);
949 
950  //TVector3 pos(0,0,0);
951  TVector3 pos(0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1),0+gRandom->Gaus(0,0.1));
952  TVector3 mom(0,0.5,0.);
953  mom *= randomSign();
954 
955 
957  rep->setPosMom(state, pos, mom);
958 
959  genfit::SharedPlanePtr origPlane = state.getPlane();
960  genfit::StateOnPlane origState(state);
961 
962  const TVector3 centerPoint(gRandom->Gaus(0,10), gRandom->Gaus(0,10), gRandom->Gaus(0,10));
963  const double radius = gRandom->Uniform(10);
964 
965  // forth
966  try {
967  rep->extrapolateToSphere(state, radius, centerPoint, false);
968  }
969  catch (genfit::Exception& e) {
970  std::cerr << e.what();
971 
972  delete rep;
973 
974  return kException;
975  }
976 
977 
978  TVector3 radiusVec = rep->getPos(state) - centerPoint;
979 
980  // compare
981  if (fabs(state.getPlane()->getNormal()*radiusVec.Unit())-1 > epsilonLen ||
982  fabs(radiusVec.Mag()-radius) > epsilonLen) {
983  if (verbose) {
984  origState.Print();
985  state.Print();
986 
987  std::cout << "state.getPlane()->getNormal()*radiusVec = " << state.getPlane()->getNormal() * radiusVec << "\n";
988  std::cout << "radiusVec.Mag()-radius = " << radiusVec.Mag() - radius << "\n";
989  }
990  delete rep;
991  return kFailed;
992  }
993 
994  delete rep;
995  return kPassed;
996 
997 }
998 
999 
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 
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 
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");
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