18 #include <PlanarMeasurement.h>
23 #include <SpacepointMeasurement.h>
39 #include <TApplication.h>
41 #include <TDatabasePDG.h>
42 #include <TEveManager.h>
43 #include <TGeoManager.h>
58 #include "TDatabasePDG.h"
78 size = backtrace(array, 10);
81 fprintf(stderr,
"Error: signal %d:\n", sig);
82 backtrace_symbols_fd(array, size, 2);
89 switch(
int(gRandom->Uniform(8))) {
113 if (gRandom->Uniform(1) > 0.5)
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)
131 double maxAbsErr = maxRelErr*max;
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 ) {
140 std::cout <<
"compareMatrices: A("<<
i<<
","<<
j<<
") = " <<
A(
i,
j) <<
" B("<<
i<<
","<<
j<<
") = " << B(
i,
j) <<
" absErr = " << absErr <<
" relErr = " << relErr <<
"\n";
152 if (!(cov.IsSymmetric())) {
154 std::cout <<
"isCovMatrix: not symmetric\n";
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))) {
163 std::cout <<
"isCovMatrix: element isnan\n";
169 std::cout <<
"isCovMatrix: negative diagonal element\n";
186 double epsilonLen = 1.E-10;
187 double epsilonMom = 1.E-10;
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));
206 const TVector3&
u = plane->getU();
207 const TVector3&
v = plane->getV();
210 pos += gRandom->Gaus() *
u;
211 pos += gRandom->Gaus() *
v;
214 mom.SetXYZ(0,0.5,gRandom->Gaus(0,0.3));
223 std::cout <<
"plane has changed unexpectedly! \n";
232 if ((
pos - rep->
getPos(state)).Mag() > epsilonLen ||
233 (mom - rep->
getMom(state)).Mag() > epsilonMom) {
237 std::cout <<
"pos difference = " << (
pos - rep->
getPos(state)).Mag() <<
"\n";
238 std::cout <<
"mom difference = " << (mom - rep->
getMom(state)).Mag() <<
"\n";
240 std::cout << std::endl;
253 static std::map<int, std::vector<TH2D*> > histoMap;
255 static const int nPDGs(5);
256 int pdgs[nPDGs]={0, 211,13,11,2212};
257 static bool fill(
true);
259 for (
int i = 0;
i<nPDGs; ++
i) {
265 Result = convert.str();
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));
276 TFile
outfile(
"deviation.root",
"recreate");
279 for (
int i = 0;
i<nPDGs; ++
i) {
281 histoMap[
pdg][0]->Write();
282 histoMap[
pdg][1]->Write();
283 histoMap[
pdg][2]->Write();
292 double epsilonLen = 5.E-5;
293 double epsilonMom = 1.E-6;
305 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
306 double mass = part->Mass();
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 );
315 double betaGamma = log(mom.Mag()/
mass);
328 double momLoss1, momLoss2;
335 mom2 = state.getMom();
336 momLoss1 = mom.Mag()-mom2.Mag();
343 std::cerr <<
"Exception in forth Extrapolation. PDG = " << pdg <<
"; mom: \n";
346 std::cerr << e.
what();
355 double backExtrapLen(0);
359 momLoss2 = mom2.Mag()-state.getMom().Mag();
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);
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);
377 std::cerr <<
"Exception in back Extrapolation. PDG = " << pdg <<
"; mom: \n";
379 std::cerr <<
"mom2: \n";
382 std::cerr << e.
what();
389 if ((rep->
getPos(origState) - rep->
getPos(state)).Mag() > epsilonLen ||
390 (rep->
getMom(origState) - rep->
getMom(state)).Mag() > epsilonMom ||
391 fabs(extrapLen + backExtrapLen) > epsilonLen) {
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";
401 std::cout << std::endl;
423 double deltaJac = 0.005;
424 double deltaNoise = 0.00001;
425 double deltaState = 3.E-6;
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));
442 mom.SetMag(gRandom->Uniform(2)+0.3);
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;
454 static const double smear = 0.2;
455 TVector3 normal(mom);
457 normal.SetXYZ(gRandom->Gaus(normal.X(), smear),
458 gRandom->Gaus(normal.Y(), smear),
459 gRandom->Gaus(normal.Z(), smear));
462 double rotAngleOrig = gRandom->Uniform(2.*TMath::Pi());
463 origPlanePtr->
rotate(rotAngleOrig);
474 normal.SetXYZ(gRandom->Gaus(normal.X(), smear),
475 gRandom->Gaus(normal.Y(), smear),
476 gRandom->Gaus(normal.Z(), smear));
481 double rotAngle = gRandom->Uniform(2.*TMath::Pi());
482 planePtr->
rotate(rotAngle);
502 extrapolatedState =
state;
507 std::cerr << e.
what();
523 std::cerr << e.
what();
534 std::cout <<
"(origState.getState() - state.getState()) ";
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) ) {
548 std::cout <<
"(jac_f * origState.getState() + c_f - extrapolatedState.getState()) ";
550 std::cout <<
"(jac_bi * origState.getState() + c_bi - extrapolatedState.getState()) ";
552 std::cout <<
"(jac_b * extrapolatedState.getState() + c_b - origState.getState()) ";
554 std::cout <<
"(jac_fi * extrapolatedState.getState() + c_fi - origState.getState()) ";
568 std::cout <<
"jac_f = ";
570 std::cout <<
"jac_bi = ";
572 std::cout << std::endl;
580 std::cout <<
"jac_b = ";
582 std::cout <<
"jac_fi = ";
584 std::cout << std::endl;
592 std::cout <<
"noise_f = ";
594 std::cout <<
"noise_bi = ";
596 std::cout << std::endl;
604 std::cout <<
"noise_b = ";
606 std::cout <<
"noise_fi = ";
608 std::cout << std::endl;
618 std::cout <<
"jac_f = ";
620 std::cout <<
"jac_f_num = ";
622 std::cout << std::endl;
640 double epsilonLen = 1.E-4;
643 double matRadius(1.);
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));
668 std::cerr << e.
what();
676 if (fabs(rep->
getPos(state).Perp() - matRadius) > epsilonLen) {
681 std::cerr <<
"radius difference = " << rep->
getPos(state).Perp() - matRadius <<
"\n";
683 std::cerr << std::endl;
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));
726 std::cerr << e.
what();
754 double epsilonLen = 1.E-4;
755 double epsilonMom = 1.E-5;
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.);
773 TVector3 linePoint(gRandom->Gaus(),
randomSign()*10+gRandom->Gaus(),gRandom->Gaus());
774 TVector3 lineDirection(gRandom->Gaus(),gRandom->Gaus(),
randomSign()*10+gRandom->Gaus());
781 std::cerr << e.
what();
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) {
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";
817 double epsilonLen = 1.E-4;
818 double epsilonMom = 1.E-5;
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));
836 TVector3
point(gRandom->Gaus(),
randomSign()*10+gRandom->Gaus(),gRandom->Gaus());
843 std::cerr << e.
what();
852 fabs((rep->
getMom(state).Unit() * state.
getPlane()->getNormal())) - 1 > epsilonMom) {
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";
875 double epsilonLen = 1.E-4;
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.);
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);
903 std::cerr << e.
what();
910 TVector3 pocaOnLine(lineDirection);
911 double t = 1./(pocaOnLine.Mag2()) * ((rep->
getPos(state)*pocaOnLine) - (linePoint*pocaOnLine));
913 pocaOnLine += linePoint;
915 TVector3 radiusVec = rep->
getPos(state) - pocaOnLine;
918 if (fabs(state.
getPlane()->getNormal()*radiusVec.Unit())-1 > epsilonLen ||
919 fabs(lineDirection*radiusVec) > epsilonLen ||
920 fabs(radiusVec.Mag()-radius) > epsilonLen) {
925 std::cout <<
"lineDirection*radiusVec = " << lineDirection * radiusVec <<
"\n";
926 std::cout <<
"radiusVec.Mag()-radius = " << radiusVec.Mag() - radius <<
"\n";
943 double epsilonLen = 1.E-4;
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.);
962 const TVector3 centerPoint(gRandom->Gaus(0,10), gRandom->Gaus(0,10), gRandom->Gaus(0,10));
963 const double radius = gRandom->Uniform(10);
970 std::cerr << e.
what();
978 TVector3 radiusVec = rep->
getPos(state) - centerPoint;
981 if (fabs(state.
getPlane()->getNormal()*radiusVec.Unit())-1 > epsilonLen ||
982 fabs(radiusVec.Mag()-radius) > epsilonLen) {
987 std::cout <<
"state.getPlane()->getNormal()*radiusVec = " << state.
getPlane()->getNormal() * radiusVec <<
"\n";
988 std::cout <<
"radiusVec.Mag()-radius = " << radiusVec.Mag() - radius <<
"\n";
1005 double epsilonLen = 1.E-3;
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.);
1020 TVector3 posOrig(state.
getPos());
1025 double step = gRandom->Uniform(-15.,15.);
1026 double extrapolatedLen(0);
1036 TVector3 posExt(state.
getPos());
1042 if (fabs(extrapolatedLen-step) > epsilonLen ||
1043 (posOrig - posExt).Mag() > fabs(step)) {
1048 std::cout <<
"extrapolatedLen-step = " << extrapolatedLen - step <<
"\n";
1049 std::cout <<
"started extrapolation from: ";
1051 std::cout <<
"extrapolated to ";
1053 std::cout <<
"difference = " << (posOrig - posExt).Mag() <<
"; step = " << step <<
"; delta = "
1054 << (posOrig - posExt).Mag() - fabs(step) <<
"\n";
1083 const double BField = 15.;
1086 gRandom->SetSeed(10);
1090 new TGeoManager(
"Geometry",
"Geane geometry");
1091 TGeoManager::Import(
"genfitGeom.root");
1101 TDatabasePDG::Instance()->GetParticle(211);
1103 const unsigned int nTests(1000);
1105 std::vector<TestCase> testCases;
1118 for (
unsigned int i=0;
i<nTests; ++
i) {
1120 for (
unsigned int j=0;
j<testCases.size(); ++
j) {
1124 testCases[
j].nPassed_++;
1127 testCases[
j].nFailed_++;
1128 std::cout <<
"failed " << testCases[
j].name_ <<
" nr " <<
i <<
"\n";
1131 testCases[
j].nException_++;
1132 std::cout <<
"exception at " << testCases[
j].name_ <<
" nr " <<
i <<
"\n";
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();
1144 for (
unsigned int j=0;
j<testCases.size(); ++
j) {
1145 testCases[
j].function_(
true);