51 #include "PlanarMeasurement.h"
63 #include <Math/SMatrix.h>
65 #include <TVectorDfwd.h>
90 TH1F* downWeightsHistosU[12];
91 TH1F* downWeightsHistosV[12];
102 TH1F* localCov12[12];
103 TH1F* localCov13[12];
104 TH1F* localCov14[12];
105 TH1F* localCov15[12];
106 TH1F* localCov23[12];
107 TH1F* localCov24[12];
108 TH1F* localCov25[12];
109 TH1F* localCov34[12];
110 TH1F* localCov35[12];
111 TH1F* localCov45[12];
120 bool writeHistoDataForLabel(
double label, TVectorD res, TVectorD measErr, TVectorD resErr, TVectorD downWeights, TVectorD localPar, TMatrixDSym localCov)
122 if (label < 1.)
return false;
124 unsigned int id = floor(label);
127 unsigned int sensor =
id & 7;
129 unsigned int ladder =
id & 31;
131 unsigned int layer =
id & 7;
132 if (layer == 7 && ladder == 2) {
134 }
else if (layer == 7 && ladder == 3) {
135 label = sensor + 9 - 3;
140 if (label > 12.)
return false;
145 resHistosU[i - 1]->Fill(res[0]);
146 resHistosV[i - 1]->Fill(res[1]);
147 mhistosU[i - 1]->Fill(res[0] / measErr[0]);
148 mhistosV[i - 1]->Fill(res[1] / measErr[1]);
149 ghistosU[i - 1]->Fill(res[0] / resErr[0]);
150 ghistosV[i - 1]->Fill(res[1] / resErr[1]);
151 downWeightsHistosU[i - 1]->Fill(downWeights[0]);
152 downWeightsHistosV[i - 1]->Fill(downWeights[1]);
153 localPar1[i - 1]->Fill(localPar(0));
154 localPar2[i - 1]->Fill(localPar(1));
155 localPar3[i - 1]->Fill(localPar(2));
156 localPar4[i - 1]->Fill(localPar(3));
157 localPar5[i - 1]->Fill(localPar(4));
173 using namespace genfit;
176 AbsFitter(), m_milleFileName(
"millefile.dat"), m_gblInternalIterations(
"THC"), m_pValueCut(0.), m_minNdf(1),
178 m_enableScatterers(
true),
179 m_enableIntermediateScatterer(
true)
189 diag =
new TFile(rootFileName.c_str(),
"RECREATE");
192 for (
int i = 0; i < 12; i++) {
193 sprintf(name,
"res_u_%i", i + 1);
194 resHistosU[
i] =
new TH1F(name,
"Residual (U)", 1000, -0.1, 0.1);
195 sprintf(name,
"res_v_%i", i + 1);
196 resHistosV[
i] =
new TH1F(name,
"Residual (V)", 1000, -0.1, 0.1);
197 sprintf(name,
"meas_pull_u_%i", i + 1);
198 mhistosU[
i] =
new TH1F(name,
"Res/Meas.Err. (U)", 1000, -20., 20.);
199 sprintf(name,
"meas_pull_v_%i", i + 1);
200 mhistosV[
i] =
new TH1F(name,
"Res/Meas.Err. (V)", 1000, -20., 20.);
201 sprintf(name,
"pull_u_%i", i + 1);
202 ghistosU[
i] =
new TH1F(name,
"Res/Res.Err. (U)", 1000, -20., 20.);
203 sprintf(name,
"pull_v_%i", i + 1);
204 ghistosV[
i] =
new TH1F(name,
"Res/Res.Err. (V)", 1000, -20., 20.);
205 sprintf(name,
"downWeights_u_%i", i + 1);
206 downWeightsHistosU[
i] =
new TH1F(name,
"Down-weights (U)", 1000, 0., 1.);
207 sprintf(name,
"downWeights_v_%i", i + 1);
208 downWeightsHistosV[
i] =
new TH1F(name,
"Down-weights (V)", 1000, 0., 1.);
209 sprintf(name,
"localPar1_%i", i + 1);
211 localPar1[
i] =
new TH1F(name,
"Residual (U)", 1000, -0.1, 0.1);
212 sprintf(name,
"localPar2_%i", i + 1);
213 localPar2[
i] =
new TH1F(name,
"Residual (U)", 1000, -0.1, 0.1);
214 sprintf(name,
"localPar3_%i", i + 1);
215 localPar3[
i] =
new TH1F(name,
"Residual (U)", 1000, -0.1, 0.1);
216 sprintf(name,
"localPar4_%i", i + 1);
217 localPar4[
i] =
new TH1F(name,
"Residual (U)", 1000, -0.1, 0.1);
218 sprintf(name,
"localPar5_%i", i + 1);
219 localPar5[
i] =
new TH1F(name,
"Residual (U)", 1000, -0.1, 0.1);
221 fitResHisto =
new TH1I(
"fit_result",
"GBL Fit Result", 21, -1, 20);
222 ndfHisto =
new TH1I(
"ndf",
"GBL Track NDF", 41, -1, 40);
223 chi2OndfHisto =
new TH1F(
"chi2_ndf",
"Track Chi2/NDF", 100, 0., 10.);
224 pValueHisto =
new TH1F(
"p_value",
"Track P-value", 100, 0., 1.);
226 stats =
new TH1I(
"stats",
"0: NDF>0 | 1: fTel&VXD | 2: all | 3: VXD | 4: SVD | 5: all - PXD | 6: fTel&SVD | 7: bTel", 10, 0, 10);
265 theta = 0.; s = 0.; ds = 0.;
266 if (steps.empty())
return;
281 for (
unsigned int i = 0; i < steps.size(); i++) {
291 sumxx += rho * (xmax -
xmin);
293 sumx2x2 += rho * (xmax * xmax - xmin *
xmin) / 2.;
295 sumx3x3 += rho * (xmax * xmax * xmax - xmin * xmin *
xmin) / 3.;
298 if (sumxx < 1.0
e-10)
return;
302 double beta = p / sqrt(p * p + mass * mass);
303 theta = (0.0136 / p / beta) * fabs(charge) * sqrt(sumxx) * (1. + 0.038 * log(sumxx));
309 double N = 1. / sumxx;
314 double ds_2 = N * (sumx3x3 - 2. * sumx2x2 * s + sumxx * s *
s);
331 bool skipMeasurement =
false;
336 bool fitQoverP =
true;
350 TVectorD scatResidual(2);
360 cout <<
"WARNING: Hits resorting in GBL interface not supported." << endl;
362 cout <<
"-------------------------------------------------------" << endl;
363 cout <<
" GBL processing genfit::Track " << endl;
364 cout <<
"-------------------------------------------------------" << endl;
365 cout <<
" # Ref. Track Points : " << npoints_all << endl;
366 cout <<
" # Meas. Points : " << npoints_meas << endl;
370 std::vector<GblPoint> listOfPoints;
372 std::vector<double> listOfLayers;
375 unsigned int seedLabel = 0;
379 TMatrixDSym clSeed(dim);
382 TMatrixD jacPointToPoint(dim, dim);
383 jacPointToPoint.UnitMatrix();
385 int n_gbl_points = 0;
386 int n_gbl_meas_points = 0;
389 double lostWeight = 0.;
392 double trackMomMag = 0.;
394 double particleCharge = 1.;
396 for (
int ipoint_meas = 0; ipoint_meas < npoints_meas; ipoint_meas++) {
400 cout <<
" ERROR: Measurement point does not have a fitter info. Track will be skipped." << endl;
406 cout <<
" ERROR: KalmanFitterInfo (with reference state) for measurement does not exist. Track will be skipped." << endl;
412 cout <<
" ERROR: Fitter info does not contain reference state. Track will be skipped." << endl;
420 TVector3 trackDir = rep->
getDir(*reference);
422 trackMomMag = rep->
getMomMag(*reference);
424 particleCharge = rep->
getCharge(*reference);
426 double particleMass = rep->
getMass(*reference);
429 double trackLen = 0.;
430 double scatTheta = 0.;
431 double scatSMean = 0.;
432 double scatDeltaS = 0.;
442 TMatrixD jacMeas2Scat(dim, dim);
443 jacMeas2Scat.UnitMatrix();
452 if (measPlanar) sensorId = measPlanar->
getPlaneId();
468 if (measPlanar2) sensorId2 = measPlanar->
getPlaneId();
472 if (sensorId != sensorId2) {
473 skipMeasurement =
true;
474 cout <<
" ERROR: Incompatible sensorIDs at measurement point " << ipoint_meas <<
". Measurement will be skipped." << endl;
484 raw_measU = raw_meas1;
485 raw_measV = raw_meas2;
489 raw_measU = raw_meas2;
490 raw_measV = raw_meas1;
493 cout <<
" ERROR: Incompatible 1D measurements at meas. point " << ipoint_meas <<
". Track will be skipped." << endl;
497 TVectorD _raw_coor(2);
501 TMatrixDSym _raw_cov(2);
513 skipMeasurement =
true;
515 cout <<
" WARNING: Measurement " << (ipoint_meas + 1) <<
" is not 2D. Measurement Will be skipped. " << endl;
521 if (!skipMeasurement) {
527 std::unique_ptr<const AbsHMatrix> HitHMatrix(raw_meas->
constructHMatrix(rep));
529 TVectorD residual = -1. * (raw_coor - HitHMatrix->Hv(state));
531 trkChi2 += residual(0) * residual(0) / raw_cov(0, 0) + residual(1) * residual(1) / raw_cov(1, 1);
534 GblPoint measPoint(jacPointToPoint);
538 TMatrixD proL2m(2, 2);
540 proL2m(0, 0) = HitHMatrix->getMatrix()(0, 3);
541 proL2m(0, 1) = HitHMatrix->getMatrix()(0, 4);
542 proL2m(1, 1) = HitHMatrix->getMatrix()(1, 4);
543 proL2m(1, 0) = HitHMatrix->getMatrix()(1, 3);
552 int label = sensorId * 10;
555 TMatrixD derGlobal(2, 12);
558 std::vector<int> labGlobal;
561 TVector3 tDir = trackDir;
563 TVector3 uDir = plane->getU();
565 TVector3 vDir = plane->getV();
567 TVector3 nDir = plane->getNormal();
573 TVector3 tLoc = TVector3(uDir.Dot(tDir), vDir.Dot(tDir), nDir.Dot(tDir));
576 double uSlope = tLoc[0] / tLoc[2];
578 double vSlope = tLoc[1] / tLoc[2];
581 double uPos = raw_coor[0];
583 double vPos = raw_coor[1];
586 derGlobal(0, 0) = 1.0;
587 derGlobal(0, 1) = 0.0;
588 derGlobal(0, 2) = - uSlope;
589 derGlobal(0, 3) = vPos * uSlope;
590 derGlobal(0, 4) = -uPos * uSlope;
591 derGlobal(0, 5) = vPos;
593 derGlobal(1, 0) = 0.0;
594 derGlobal(1, 1) = 1.0;
595 derGlobal(1, 2) = - vSlope;
596 derGlobal(1, 3) = vPos * vSlope;
597 derGlobal(1, 4) = -uPos * vSlope;
598 derGlobal(1, 5) = -uPos;
600 labGlobal.push_back(label + 1);
601 labGlobal.push_back(label + 2);
602 labGlobal.push_back(label + 3);
603 labGlobal.push_back(label + 4);
604 labGlobal.push_back(label + 5);
605 labGlobal.push_back(label + 6);
611 TVector3 detPos = plane->getO();
616 TVector3 pred = detPos + uPos * uDir + vPos * vDir;
620 double xPred = pred[0];
621 double yPred = pred[1];
622 double zPred = pred[2];
625 double tn = tDir.Dot(nDir);
632 for (
int row = 0; row < 3; row++)
634 drdm(row,
col) -= tDir[row] * nDir[
col] / tn;
642 dmdg(0, 0) = 1.; dmdg(0, 4) = -zPred; dmdg(0, 5) = yPred;
643 dmdg(1, 1) = 1.; dmdg(1, 3) = zPred; dmdg(1, 5) = -xPred;
644 dmdg(2, 2) = 1.; dmdg(2, 3) = -yPred; dmdg(2, 4) = xPred;
650 TMatrixD drldrg(3, 3);
652 drldrg(0, 0) = uDir[0]; drldrg(0, 1) = uDir[1]; drldrg(0, 2) = uDir[2];
653 drldrg(1, 0) = vDir[0]; drldrg(1, 1) = vDir[1]; drldrg(1, 2) = vDir[2];
662 TMatrixD drldg(3, 6);
663 drldg = drldrg * (drdm * dmdg);
674 derGlobal(0, 6) = drldg(0, 0); labGlobal.push_back(offset + 1);
675 derGlobal(0, 7) = drldg(0, 1); labGlobal.push_back(offset + 2);
676 derGlobal(0, 8) = drldg(0, 2); labGlobal.push_back(offset + 3);
677 derGlobal(0, 9) = drldg(0, 3); labGlobal.push_back(offset + 4);
678 derGlobal(0, 10) = drldg(0, 4); labGlobal.push_back(offset + 5);
679 derGlobal(0, 11) = drldg(0, 5); labGlobal.push_back(offset + 6);
681 derGlobal(1, 6) = drldg(1, 0);
682 derGlobal(1, 7) = drldg(1, 1);
683 derGlobal(1, 8) = drldg(1, 2);
684 derGlobal(1, 9) = drldg(1, 3);
685 derGlobal(1, 10) = drldg(1, 4);
686 derGlobal(1, 11) = drldg(1, 5);
691 listOfPoints.push_back(measPoint);
692 listOfLayers.push_back((
unsigned int) sensorId);
697 GblPoint dummyPoint(jacPointToPoint);
698 listOfPoints.push_back(dummyPoint);
699 listOfLayers.push_back((
unsigned int) sensorId);
701 skipMeasurement =
false;
703 cout <<
" Dummy point inserted. " << endl;
712 if (ipoint_meas < npoints_meas - 1) {
715 cout <<
" ERROR: Measurement point does not have a fitter info. Track will be skipped." << endl;
721 cout <<
" ERROR: KalmanFitterInfo (with reference state) for measurement does not exist. Track will be skipped." << endl;
738 s2 = scatSMean + scatDeltaS * scatDeltaS / (scatSMean - s1);
739 theta1 = sqrt(scatTheta * scatTheta * scatDeltaS * scatDeltaS / (scatDeltaS * scatDeltaS + (scatSMean - s1) * (scatSMean - s1)));
740 theta2 = sqrt(scatTheta * scatTheta * (scatSMean - s1) * (scatSMean - s1) / (scatDeltaS * scatDeltaS + (scatSMean - s1) * (scatSMean - s1)));
750 if (s2 < 1.e-4 || s2 >= trackLen - 1.
e-4 || s2 <= 1.
e-4) {
751 cout <<
" WARNING: GBL points will be too close. GBLTrajectory construction might fail. Let's try it..." << endl;
760 if (ipoint_meas < npoints_meas - 1) {
761 if (theta2 > scatEpsilon) {
771 cout <<
" ERROR: The extrapolation to measurement point " << (ipoint_meas + 2) <<
" stepped back by " << nextStep <<
"cm !!! Track will be skipped." << endl;
786 cout <<
" ERROR: Extrapolation failed. Track will be skipped." << endl;
795 if (ipoint_meas < npoints_meas - 1) {
797 if (theta1 > scatEpsilon) {
802 double c1 = trackDir.Dot(plane->getU());
803 double c2 = trackDir.Dot(plane->getV());
804 TMatrixDSym scatCov(2);
805 scatCov(0, 0) = 1. - c2 *
c2;
806 scatCov(1, 1) = 1. - c1 * c1;
807 scatCov(0, 1) = c1 *
c2;
808 scatCov(1, 0) = c1 *
c2;
809 scatCov *= theta1 * theta1 / (1. - c1 * c1 - c2 *
c2) / (1. - c1 * c1 - c2 * c2) ;
812 GblPoint& lastPoint = listOfPoints.at(n_gbl_points - 1);
817 if (theta2 > scatEpsilon) {
821 TMatrixDSym scatCov(2);
823 scatCov(0, 0) = theta2 * theta2;
824 scatCov(1, 1) = theta2 * theta2;
828 listOfPoints.push_back(scatPoint);
829 listOfLayers.push_back(((
unsigned int) sensorId) + 0.5);
839 if (n_gbl_meas_points > 1) {
845 traj =
new GblTrajectory(listOfPoints, seedLabel, clSeed, fitQoverP);
860 pvalue = TMath::Prob(Chi2, Ndf);
868 fitResHisto->Fill(fitRes);
873 cout <<
" Ref. Track Chi2 : " << trkChi2 << endl;
874 cout <<
" Ref. end momentum : " << trackMomMag <<
" GeV/c ";
876 if (particleCharge < 0.)
877 cout <<
"(electron)";
879 cout <<
"(positron)";
882 cout <<
"------------------ GBL Fit Results --------------------" << endl;
883 cout <<
" Fit q/p parameter : " << ((fitQoverP) ? (
"True") : (
"False")) << endl;
884 cout <<
" Valid trajectory : " << ((traj->
isValid()) ? (
"True") : (
"False")) << endl;
885 cout <<
" Fit result : " << fitRes <<
" (0 for success)" << endl;
886 cout <<
" # GBL meas. points : " << n_gbl_meas_points << endl;
887 cout <<
" # GBL all points : " << n_gbl_points << endl;
888 cout <<
" GBL track NDF : " << Ndf <<
" (-1 for failure)" << endl;
889 cout <<
" GBL track Chi2 : " << Chi2 << endl;
890 cout <<
" GBL track P-value : " << pvalue;
894 cout <<
"-------------------------------------------------------" << endl;
898 bool hittedLayers[12];
899 for (
int hl = 0; hl < 12; hl++) {
900 hittedLayers[hl] =
false;
914 for (
unsigned int p = 0;
p < listOfPoints.size();
p++) {
915 unsigned int label =
p + 1;
917 TVectorD residuals(2);
918 TVectorD measErrors(2);
919 TVectorD resErrors(2);
920 TVectorD downWeights(2);
922 if (!listOfPoints.at(
p).hasMeasurement())
928 unsigned int id = listOfLayers.at(
p);
931 unsigned int sensor =
id & 7;
933 unsigned int ladder =
id & 31;
935 unsigned int layer =
id & 7;
937 if (layer == 7 && ladder == 2) {
939 }
else if (layer == 7 && ladder == 3) {
945 hittedLayers[l - 1] =
true;
947 TVectorD localPar(5);
948 TMatrixDSym localCov(5);
951 traj->
getMeasResults(label, numRes, residuals, measErrors, resErrors, downWeights);
956 if (!writeHistoDataForLabel(listOfLayers.at(
p), residuals, measErrors, resErrors, downWeights, localPar, localCov))
966 cout <<
" GBL Track written to Millepede II binary file." << endl;
967 cout <<
"-------------------------------------------------------" << endl;
973 chi2OndfHisto->Fill(Chi2 / Ndf);
974 pValueHisto->Fill(TMath::Prob(Chi2, Ndf));