25 #include <PlanarMeasurement.h>
30 #include <SpacepointMeasurement.h>
51 #include <TApplication.h>
53 #include <TDatabasePDG.h>
54 #include <TEveManager.h>
55 #include <TGeoManager.h>
65 #include "TDatabasePDG.h"
77 size = backtrace(array, 10);
80 fprintf(stderr,
"Error: signal %d:\n", sig);
81 backtrace_symbols_fd(array, size, 2);
88 switch(
int(gRandom->Uniform(8))) {
112 if (gRandom->Uniform(1) > 0.5)
126 #include <valgrind/callgrind.h>
128 #define CALLGRIND_START_INSTRUMENTATION
129 #define CALLGRIND_STOP_INSTRUMENTATION
130 #define CALLGRIND_DUMP_STATS
134 std::cout<<
"main"<<std::endl;
135 gRandom->SetSeed(14);
138 const unsigned int nEvents = 1000;
139 const unsigned int nMeasurements = 10;
140 const double BField = 20.;
142 const double theta = 110;
143 const double thetaDetPlane = 90;
144 const double phiDetPlane = 0;
145 const double pointDist = 3.;
146 const double resolution = 0.05;
148 const double resolutionWire = 5*resolution;
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;
157 const double outlierProb = -0.1;
158 const double outlierRange = 2;
160 const double hitSwitchProb = -0.1;
162 const int splitTrack = -5;
163 const bool fullMeasurement =
false;
178 const int nIter = 20;
179 const double dPVal = 1.E-3;
181 const bool resort =
false;
182 const bool prefit =
false;
183 const bool refit =
false;
185 const bool twoReps =
false;
187 const bool checkPruning =
true;
191 const bool smearPosMom =
true;
192 const double chargeSwitchProb = -0.1;
193 const double posSmear = 10*resolution;
194 const double momSmear = 5. /180.*TMath::Pi();
195 const double momMagSmear = 0.2;
196 const double zSmearFac = 2;
199 const bool matFX =
false;
201 const bool debug =
false;
202 const bool onlyDisplayFailed =
false;
204 std::vector<genfit::eMeasurementType> measurementTypes;
205 for (
unsigned int i = 0;
i<nMeasurements; ++
i) {
263 new TGeoManager(
"Geometry",
"Geane geometry");
264 TGeoManager::Import(
"genfitGeom.root");
269 const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
280 gROOT->SetStyle(
"Plain");
281 gStyle->SetPalette(1);
282 gStyle->SetOptFit(1111);
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);
290 TH1D *hqopPu =
new TH1D(
"hqopPu",
"q/p pull",200,-6.,6.);
291 TH1D *pVal =
new TH1D(
"pVal",
"p-value",100,0.,1.00000001);
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.);
298 TH1D *weights =
new TH1D(
"weights",
"Daf vs true weights",500,-1.01,1.01);
300 TH1D *trackLenRes =
new TH1D(
"trackLenRes",
"(trueLen - FittedLen) / trueLen",500,-0.01,0.01);
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);
320 for (
unsigned int iEvent=0; iEvent<
nEvents; ++iEvent){
327 if (debug || (iEvent)%10==0)
328 std::cout <<
"Event Nr. " << iEvent <<
" ";
329 else std::cout <<
". ";
330 if (debug || (iEvent+1)%10==0)
338 secondTrack =
nullptr;
342 TVector3
pos(0, 0, 0);
343 TVector3 mom(1.,0,0);
344 mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
346 mom.SetTheta(theta*TMath::Pi()/180);
347 mom.SetMag(momentum);
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);
355 std::cout <<
"start values \n";
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));
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()));
380 if (chargeSwitchProb > gRandom->Uniform(1.))
384 if (chargeSwitchProb > gRandom->Uniform(1.))
402 std::vector< std::vector<genfit::AbsMeasurement*> > measurements;
404 std::vector<bool> outlierTrue;
407 std::vector<int> leftRightTrue;
413 for (
unsigned int i=0;
i<measurementTypes.size(); ++
i){
414 trueLen =
i*pointDist;
416 measurements.push_back(measurementCreator.
create(measurementTypes[
i], trueLen, outlier, lr));
417 outlierTrue.push_back(outlier);
418 leftRightTrue.push_back(lr);
420 assert(measurementTypes.size() == leftRightTrue.size());
421 assert(measurementTypes.size() == outlierTrue.size());
424 std::cerr<<
"Exception, next track"<<std::endl;
425 std::cerr << e.
what();
429 if (debug) std::cout <<
"... done creating measurements \n";
434 TVectorD seedState(6);
435 TMatrixDSym seedCov(6);
451 for(
unsigned int i=0;
i<measurements.size(); ++
i){
452 if (splitTrack > 0 && (
int)
i >= splitTrack)
454 if (
i>0 && hitSwitchProb > gRandom->Uniform(1.))
463 if (splitTrack > 0) {
464 for(
unsigned int i=splitTrack;
i<measurements.size(); ++
i){
465 if (
i>0 && hitSwitchProb > gRandom->Uniform(1.))
481 if (debug) std::cout<<
"Starting the fitter"<<std::endl;
493 if (debug) std::cout<<
"fitter is finished!"<<std::endl;
496 std::cerr << e.
what();
497 std::cerr <<
"Exception, next track" << std::endl;
501 if (splitTrack > 0) {
502 if (debug) fitTrack->
Print(
"C");
503 if (debug) secondTrack->
Print(
"C");
505 if (fullMeasurement) {
512 if (debug) fitTrack->
Print(
"C");
515 if (debug) std::cout<<
"Starting the fitter"<<std::endl;
517 if (debug) std::cout<<
"fitter is finished!"<<std::endl;
520 std::cerr << e.
what();
521 std::cerr <<
"Exception, next track" << std::endl;
528 std::cout<<
"Trying to fit again "<<std::endl;
535 fitTrack->
Print(
"C");
543 if (!onlyDisplayFailed && iEvent < 1000) {
544 std::vector<genfit::Track*>
event;
545 event.push_back(fitTrack);
547 event.push_back(secondTrack);
550 else if (onlyDisplayFailed &&
565 nUNConvergedFits += 1;
571 nConvergedFitsSecond += 1;
575 nUNConvergedFitsSecond += 1;
582 std::cout <<
"Track could not be fitted successfully! Fit is not converged! \n";
589 std::cout <<
"Track has no TrackPoint with fitterInfo! \n";
594 std::cout <<
"state before extrapolating back to reference plane \n";
603 std::cerr<<
"Exception, next track"<<std::endl;
604 std::cerr << e.
what();
610 const TVectorD& referenceState = stateRefOrig.
getState();
615 double pval = fitter->
getPVal(fitTrack, rep);
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]));
624 hqopPu->Fill( (state[0]-referenceState[0]) / sqrt(cov[0][0]) );
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]) );
633 trackLenRes->Fill( (trueLen - fitTrack->
getTrackLen(rep)) / trueLen );
636 std::cout <<
"true track length = " << trueLen <<
"; fitted length = " << fitTrack->
getTrackLen(rep) <<
"\n";
637 std::cout <<
"fitted tof = " << fitTrack->
getTOF(rep) <<
" ns\n";
641 std::cerr << e.
what();
642 std::cout <<
"could not get TrackLen or TOF! \n";
648 if (dynamic_cast<genfit::DAF*>(fitter) !=
nullptr) {
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] <<
" ";
662 std::cout <<
" l/r: " << leftRightTrue[
i];
665 int trueSide = leftRightTrue[
i];
666 if (trueSide == 0)
continue;
667 if (outlierTrue[
i])
continue;
669 if(dafWeightLR.size() != 2)
672 double weightCorrectSide, weightWrongSide;
675 weightCorrectSide = dafWeightLR[0];
676 weightWrongSide = dafWeightLR[1];
679 weightCorrectSide = dafWeightLR[1];
680 weightWrongSide = dafWeightLR[0];
682 weightWrongSide -= 1.;
684 weights->Fill(weightCorrectSide);
685 weights->Fill(weightWrongSide);
687 if (weightCorrectSide>maxWeight) maxWeight = weightCorrectSide;
696 if (outlierTrue[
i]) {
697 for (
unsigned int j=0;
j<dafWeights.size(); ++
j){
698 weights->Fill(dafWeights[
j]-1);
701 else if (leftRightTrue[i] == 0) {
702 for (
unsigned int j=0;
j<dafWeights.size(); ++
j){
703 weights->Fill(dafWeights[
j]);
717 for (
unsigned int i=0;
i<1; ++
i) {
721 bool first(
false), last(
false);
725 if (gRandom->Uniform() < 0.5) trClone.
prune(
"C");
726 if (gRandom->Uniform() < 0.5) {
730 if (gRandom->Uniform() < 0.5) {
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");
772 std::cout<<
" pruned track: ";
777 std::cerr << e.
what();
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;
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;
815 if (debug) std::cout<<
"Draw histograms ...";
817 TCanvas* c1 =
new TCanvas();
821 hmomRes->Fit(
"gaus");
845 TCanvas*
c2 =
new TCanvas();
875 TCanvas* c3 =
new TCanvas();
879 trackLenRes->Fit(
"gaus");
884 if (debug) std::cout<<
"... done"<<std::endl;