147 Double_t upperScaleLimit = 6., bin = 400.;
149 fHM->
Create1<TH1D>(
"fhDistanceCenterToExtrapolatedTrack",
150 "fhDistanceCenterToExtrapolatedTrack;Distance fitted "
151 "center to extrapolated track;Number of entries",
152 bin, 0., upperScaleLimit);
153 fHM->
Create1<TH1D>(
"fhDistanceCorrected",
"fhDistanceCorrected;Distance a [cm];A.U.", bin, 0., upperScaleLimit);
154 fHM->
Create1<TH1D>(
"fhDifferenceX",
"fhDifferenceX;Difference in X (fitted center - extrapolated track);A.U.", bin,
155 -upperScaleLimit, upperScaleLimit);
156 fHM->
Create1<TH1D>(
"fhDifferenceY",
"fhDifferenceY;Difference in Y (fitted center - extrapolated track);A.U.", bin,
157 -upperScaleLimit, upperScaleLimit);
159 fHM->
Create1<TH1D>(
"fhDistanceUncorrected",
"fhDistanceUncorrected;Distance a [cm];A.U.", bin, 0., upperScaleLimit);
160 fHM->
Create1<TH1D>(
"fhDifferenceXUncorrected",
"fhDifferenceXUncorrected;Difference in X uncorrected [cm];A.U.", bin,
161 -upperScaleLimit, upperScaleLimit);
162 fHM->
Create1<TH1D>(
"fhDifferenceYUncorrected",
"fhDifferenceYUncorrected;Difference in Y uncorrected [cm];A.U.", bin,
163 -upperScaleLimit, upperScaleLimit);
165 fHM->
Create1<TH1D>(
"fhDistanceIdeal",
"fhDistanceIdeal;Distance a [cm];A.U.", bin, 0., upperScaleLimit);
166 fHM->
Create1<TH1D>(
"fhDifferenceXIdeal",
"fhDifferenceXIdeal;Difference in X ideal [cm];A.U.", bin, -upperScaleLimit,
168 fHM->
Create1<TH1D>(
"fhDifferenceYIdeal",
"fhDifferenceYIdeal;Difference in Y ideal [cm];A.U.", bin, -upperScaleLimit,
172 "fHistoDiffX;Histogram difference between corrected and "
173 "ideal X positions;A.U.",
174 bin, 0., upperScaleLimit);
176 "fHistoDiffY;Histogram difference between corrected and "
177 "ideal Y positions;A.U.",
178 bin, 0., upperScaleLimit);
180 fHM->
Create1<TH1D>(
"fHistoBoA",
"fHistoBoA;Histogram B axis over A axis;A.U.", bin, 0., upperScaleLimit);
186 Double_t upperScaleLimit = 6., bin = 400.;
187 TString strCorrX =
"fhDifferenceCorrectedX_mirror_tile_", strCorrY =
"fhDifferenceCorrectedY_mirror_tile_",
188 strDiffCorrX =
"DiffCorrX_mirror_tile_", strDiffCorrY =
"DiffCorrY_mirror_tile_";
189 stringstream ssDiffCorrX, ssDiffCorrY, ssCorrX, ssCorrY;
190 TString strUncorrX =
"fhDifferenceUncorrectedX_mirror_tile_", strUncorrY =
"fhDifferenceUncorrectedY_mirror_tile_",
191 strDiffUncorrX =
"DiffUncorrX_mirror_tile_", strDiffUncorrY =
"DiffUncorrY_mirror_tile_";
192 stringstream ssDiffUncorrX, ssDiffUncorrY, ssUncorrX, ssUncorrY;
193 TString strIdealX =
"fhDifferenceIdealX_mirror_tile_", strIdealY =
"fhDifferenceIdealY_mirror_tile_",
194 strDiffIdealX =
"DiffIdealX_mirror_tile_", strDiffIdealY =
"DiffIdealY_mirror_tile_";
195 stringstream ssDiffIdealX, ssDiffIdealY, ssIdealX, ssIdealY;
197 cout << endl <<
"Init histo: " << endl << endl;
198 for (Int_t j = 0; j < 4; j++) {
199 for (Int_t i = 0; i < 10; i++) {
200 ssDiffCorrX << strDiffCorrX << j <<
"_" << i;
201 ssDiffCorrY << strDiffCorrY << j <<
"_" << i;
202 ssCorrX << strCorrX << j <<
"_" << i;
203 ssCorrY << strCorrY << j <<
"_" << i;
205 ssCorrX.str().c_str(),
";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
207 ssCorrY.str().c_str(),
";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
209 ssDiffUncorrX << strDiffUncorrX << j <<
"_" << i;
210 ssDiffUncorrY << strDiffUncorrY << j <<
"_" << i;
211 ssUncorrX << strUncorrX << j <<
"_" << i;
212 ssUncorrY << strUncorrY << j <<
"_" << i;
214 new TH1D(ssUncorrX.str().c_str(),
";Difference in X (fitted center - extrapolated track);A.U.", bin, 0.,
217 new TH1D(ssUncorrY.str().c_str(),
";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0.,
220 ssDiffIdealX << strDiffIdealX << j <<
"_" << i;
221 ssDiffIdealY << strDiffIdealY << j <<
"_" << i;
222 ssIdealX << strIdealX << j <<
"_" << i;
223 ssIdealY << strIdealY << j <<
"_" << i;
225 ssIdealX.str().c_str(),
";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
227 ssIdealY.str().c_str(),
";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
236 ssDiffUncorrX.str(
"");
238 ssDiffUncorrY.str(
"");
240 ssDiffIdealX.str(
"");
242 ssDiffIdealY.str(
"");
248 Double_t xMin = -120., xMax = 120., nBinsX1 = 60, yMax = 200., range = 3.;
251 for (Int_t k = 1; k < 3; k++) {
253 fHM->
Create3<TH3D>(
"fhRingTrackDistVsXYTruematch" + ss.str(),
254 "fhRingTrackDistVsXYTruematch" + ss.str() +
";X [cm];Y [cm];Ring-track distance [cm]", 60, xMin,
255 xMax, 104, 110., 200., 100, 0., 5.);
256 fHM->
Create2<TH2D>(
"fhRingTrackDistVsXTruematch" + ss.str(),
257 "fhRingTrackDistVsXTruematch" + ss.str() +
";X [cm];Ring-track distance [cm]", nBinsX1, xMin,
259 fHM->
Create2<TH2D>(
"fhRingTrackDistVsYTruematch" + ss.str(),
260 "fhRingTrackDistVsYTruematch" + ss.str() +
";Abs(Y) [cm];Ring-track distance [cm]", 34, 110.,
262 fHM->
Create3<TH3D>(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str(),
263 "fhRingTrackDistDiffXRingVsXYTruematch" + ss.str() +
";X [cm];Y [cm];X Ring-track distance [cm]",
264 60, xMin, xMax, 104, 110, 200, 200, -range, range);
265 fHM->
Create3<TH3D>(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str(),
266 "fhRingTrackDistDiffYRingVsXYTruematch" + ss.str() +
";X [cm];Y [cm];Y Ring-track distance [cm]",
267 60, xMin, xMax, 104, 110, 200, 200, -range, range);
271 fHM->
Create1<TH1D>(
"fDistUncorr",
"fDistUncorr;Uncorrected Distance;Number of entries", 600, -10., 10);
272 fHM->
Create1<TH1D>(
"fDistCorr",
"fDistCorr;Corrected Distance;Number of entries", 600, -10., 10);
335 cout <<
"CbmRichMirrorSortingCorrection: Event #" <<
fEventNb << endl;
336 TVector3 momentum, outPos, outPosUnCorr, outPosIdeal;
337 Double_t constantePMT = 0., trackX = 0., trackY = 0.;
338 vector<Double_t> vect(2, 0), ptM(3, 0), ptC(3, 0), ptCIdeal(3, 0), ptR1(3, 0), ptR2Center(3, 0), ptR2Mirr(3, 0),
339 ptPR2(3, 0), ptPMirr(3, 0), normalPMT(3, 0);
340 vector<Double_t> ptR2CenterUnCorr(3, 0), ptR2CenterIdeal(3, 0), ptR2MirrUnCorr(3, 0), ptR2MirrIdeal(3, 0),
341 ptPMirrUnCorr(3, 0), ptPMirrIdeal(3, 0), ptPR2UnCorr(3, 0), ptPR2Ideal(3, 0);
343 ptC.at(0) = 0., ptC.at(1) = 132.594000, ptC.at(2) = 54.267226;
344 TVector3 mirrorPoint, dirCos,
pos;
345 Double_t nx = 0., ny = 0., nz = 0.;
350 cout <<
"Nb of rings in evt = " <<
fRichRings->GetEntriesFast() << endl << endl;
354 for (Int_t iGlobalTrack = 0; iGlobalTrack <
fGlobalTracks->GetEntriesFast(); iGlobalTrack++) {
362 cout <<
"Error richInd < 0" << endl;
367 cout <<
"Error ring == NULL!" << endl;
374 if (NULL == cbmRichTrackMatch) {
377 cout <<
"Nof true hits = " << cbmRichTrackMatch->
GetNofTrueHits() << endl;
378 cout <<
"Nof wrong hits = " << cbmRichTrackMatch->
GetNofWrongHits() << endl;
382 if (mcRichTrackId < 0)
continue;
384 if (mcStsTrackId != mcRichTrackId) {
385 cout <<
"Error StsTrackIndex and TrackIndex from Ring do not match!" << endl;
389 if (!mcTrack)
continue;
400 if (pTrack == NULL) {
401 cout <<
"CbmRichMirrorSortingCorrection::Exec : pTrack = NULL." << endl;
404 trackX = pTrack->GetX(), trackY = pTrack->GetY();
405 cout <<
"Track: " << trackX <<
", " << trackY << endl;
410 Int_t pdg = TMath::Abs(mcTrack->
GetPdgCode());
411 if (trackMotherId == -1) {
414 for (Int_t iMirrPt = 0; iMirrPt <
fMirrorPoints->GetEntriesFast(); iMirrPt++) {
416 if (mirrPoint == 0) {
420 if (mirrPoint->GetTrackID() == mcRichTrackId) {
424 ptM.at(0) = mirrPoint->GetX(), ptM.at(1) = mirrPoint->GetY(), ptM.at(2) = mirrPoint->GetZ();
426 mirrNode = gGeoManager->FindNode(ptM.at(0), ptM.at(1), ptM.at(2));
428 string str1 = gGeoManager->GetPath(), str2 =
"mirror_tile_", str3 =
"";
429 std::size_t found = str1.find(str2);
430 if (found != std::string::npos) {
432 Int_t end = str2.length() + 3;
433 str3 = str1.substr(found, end);
435 cout <<
"Mirror ID: " << str3 << endl;
438 TGeoNavigator* navi = gGeoManager->GetCurrentNavigator();
440 ptCIdeal.at(0) = navi->GetCurrentMatrix()->GetTranslation()[0];
441 ptCIdeal.at(1) = navi->GetCurrentMatrix()->GetTranslation()[1];
442 ptCIdeal.at(2) = navi->GetCurrentMatrix()->GetTranslation()[2];
443 cout <<
"Sphere center coordinates of the aligned mirror tile, "
445 << ptCIdeal.at(0) <<
", " << ptCIdeal.at(1) <<
", " << ptCIdeal.at(2) <<
"}" << endl;
446 for (Int_t iRefPt = 0; iRefPt <
fRefPlanePoints->GetEntriesFast(); iRefPt++) {
449 if (refPlanePoint->GetTrackID() == mcRichTrackId) {
453 ptR1.at(0) = refPlanePoint->GetX(), ptR1.at(1) = refPlanePoint->GetY(), ptR1.at(2) = refPlanePoint->GetZ();
454 cout <<
"Refl plane point coo = {" << ptR1[0] <<
", " << ptR1[1] <<
", " << ptR1[2] <<
"}" << endl;
455 ComputeR2(ptR2Center, ptR2Mirr, ptM, ptC, ptR1, navi,
"Corrected", str3);
456 ComputeR2(ptR2CenterUnCorr, ptR2MirrUnCorr, ptM, ptC, ptR1, navi,
"Uncorrected", str3);
457 ComputeR2(ptR2CenterIdeal, ptR2MirrIdeal, ptM, ptCIdeal, ptR1, navi,
"Uncorrected", str3);
458 ComputeP(ptPMirr, ptPR2, normalPMT, ptM, ptR2Mirr, constantePMT);
459 ComputeP(ptPMirrUnCorr, ptPR2UnCorr, normalPMT, ptM, ptR2MirrUnCorr, constantePMT);
460 ComputeP(ptPMirrIdeal, ptPR2Ideal, normalPMT, ptM, ptR2MirrIdeal, constantePMT);
461 cout <<
"PMT points mirr coordinates before rotation = {" << ptPMirr[0] <<
", " << ptPMirr[1] <<
", "
462 << ptPMirr[2] <<
"}" << endl;
463 cout <<
"PMT points mirr uncorr coordinates before rotation = {" << ptPMirrUnCorr[0] <<
", "
464 << ptPMirrUnCorr[1] <<
", " << ptPMirrUnCorr[2] <<
"}" << endl;
465 cout <<
"PMT points mirr ideal coordinates before rotation = {" << ptPMirrIdeal[0] <<
", "
466 << ptPMirrIdeal[1] <<
", " << ptPMirrIdeal[2] <<
"}" << endl;
468 TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
471 <<
"New PMT points coordinates = {" << outPos.x() <<
", " << outPos.y() <<
", " << outPos.z() <<
"}"
473 TVector3 inPosUnCorr(ptPMirrUnCorr.at(0), ptPMirrUnCorr.at(1), ptPMirrUnCorr.at(2));
475 cout <<
"New mirror points coordinates = {" << outPosUnCorr.x() <<
", " << outPosUnCorr.y() <<
", "
476 << outPosUnCorr.z() <<
"}" << endl;
477 TVector3 inPosIdeal(ptPMirrIdeal.at(0), ptPMirrIdeal.at(1), ptPMirrIdeal.at(2));
479 cout <<
"New mirror points coordinates = {" << outPosIdeal.x() <<
", " << outPosIdeal.y() <<
", "
480 << outPosIdeal.z() <<
"}" << endl
483 cout <<
"pTrack [X,Y]: " << pTrack->GetX() <<
", " << pTrack->GetY() << endl;
485 pTrack->SetX(outPosUnCorr.x());
486 pTrack->SetY(outPosUnCorr.y());
489 pTrack->SetX(outPos.x());
490 pTrack->SetY(outPos.y());
493 pTrack->SetX(outPosIdeal.x());
494 pTrack->SetY(outPosIdeal.y());
497 FillHistProjection(outPosIdeal, outPosUnCorr, outPos, ringL, normalPMT, constantePMT, str3);
499 cout <<
"pTrack [X,Y]: " << pTrack->GetX() <<
", " << pTrack->GetY() << endl;
504 cout <<
"No mirror points registered." << endl;
508 cout <<
"Not a mother particle." << endl;
514 cout <<
"CbmRichMirrorSortingCorrection::Exec No rings in event were found." << endl;
524 Int_t pmtTrackID, pmtMotherID;
525 Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0., scalarProd = 0.;
526 Double_t pmtPt[] = {0., 0., 0.};
527 Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
535 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
537 pmtTrackID = pmtPoint->GetTrackID();
540 a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
544 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
546 pmtTrackID = pmtPoint->GetTrackID();
550 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
551 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
553 b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
558 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
560 pmtTrackID = pmtPoint->GetTrackID();
564 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
565 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
567 && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2) + TMath::Power(b[1] - pmtPoint->GetY(), 2)
568 + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
570 c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
576 k = (b[0] - a[0]) / (c[0] - a[0]);
577 if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0 || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
578 cout <<
"Error in normal calculation, vect_AB and vect_AC are collinear." << endl;
581 buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
582 buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
583 buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
585 buffNormX / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
587 buffNormY / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
589 buffNormZ / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
593 scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
594 + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
599 * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY() + normalPMT.at(2) * pmtPoint1->GetZ());
601 scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
602 + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
605 scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
606 + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
611 vector<Double_t> ptM, vector<Double_t> ptC, vector<Double_t> ptR1,
612 TGeoNavigator* navi, TString option, TString mirrorTileName)
616 vector<Double_t> normalMirr(3), ptCNew(3), ptTileCenter(3);
617 Double_t t1 = 0., t2 = 0., t3 = 0.;
619 if (option ==
"Corrected") {
622 Int_t lineCounter = 1, lineIndex = 0;
625 string fileLine =
"", strMisX =
"", strMisY =
"";
626 Double_t misX = 0., misY = 0.;
642 if (corrFile.is_open()) {
643 while (!corrFile.eof()) {
644 getline(corrFile, fileLine);
645 lineIndex = fileLine.find(mirrorTileName, 0);
646 if (lineIndex != string::npos) {
669 cout <<
"Error in CbmRichCorrection: unable to open parameter file!" << endl;
670 cout <<
"Parameter file path: " << str << endl << endl;
677 cout <<
"Correction parameters = " << misX <<
", " << misY << endl;
678 ptCNew.at(0) = ptC.at(0) + misX;
679 ptCNew.at(1) = ptC.at(1) + misY;
680 ptCNew.at(2) = ptC.at(2);
681 ptTileCenter.at(0) = navi->GetMotherMatrix()->GetTranslation()[0];
682 ptTileCenter.at(1) = navi->GetMotherMatrix()->GetTranslation()[1];
683 ptTileCenter.at(2) = navi->GetMotherMatrix()->GetTranslation()[2];
685 Double_t
x = 0.,
y = 0., z = 0., dist = 0., dist2 = 0., z2 = 0.;
686 x = TMath::Power(ptCNew.at(0) - ptTileCenter.at(0), 2);
687 y = TMath::Power(ptCNew.at(1) - ptTileCenter.at(1), 2);
688 z = TMath::Power(ptCNew.at(2) - ptTileCenter.at(2), 2);
689 dist = TMath::Sqrt(
x +
y + z);
690 z2 = ptTileCenter.at(2) - TMath::Sqrt(TMath::Power(300, 2) -
x -
y) - ptCNew.at(2);
692 dist2 = TMath::Sqrt(
x +
y + TMath::Power(z2 - ptTileCenter.at(2), 2));
695 cout <<
"Sphere center coordinates of the rotated mirror tile, after "
697 << ptCNew.at(0) <<
", " << ptCNew.at(1) <<
", " << ptCNew.at(2) <<
"}" << endl;
699 else if (option ==
"Uncorrected") {
702 cout <<
"Sphere center coordinates of the rotated mirror tile, without "
704 << ptCNew.at(0) <<
", " << ptCNew.at(1) <<
", " << ptCNew.at(2) <<
"}" << endl;
712 normalMirr.at(0) = (ptCNew.at(0) - ptM.at(0))
713 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
714 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
715 normalMirr.at(1) = (ptCNew.at(1) - ptM.at(1))
716 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
717 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
718 normalMirr.at(2) = (ptCNew.at(2) - ptM.at(2))
719 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
720 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
723 t1 = ((ptR1.at(0) - ptM.at(0)) * (ptCNew.at(0) - ptM.at(0)) + (ptR1.at(1) - ptM.at(1)) * (ptCNew.at(1) - ptM.at(1))
724 + (ptR1.at(2) - ptM.at(2)) * (ptCNew.at(2) - ptM.at(2)))
725 / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
726 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
727 ptR2Center.at(0) = 2 * (ptM.at(0) + t1 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
728 ptR2Center.at(1) = 2 * (ptM.at(1) + t1 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
729 ptR2Center.at(2) = 2 * (ptM.at(2) + t1 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
731 ((ptR1.at(0) - ptCNew.at(0)) * (ptCNew.at(0) - ptM.at(0)) + (ptR1.at(1) - ptCNew.at(1)) * (ptCNew.at(1) - ptM.at(1))
732 + (ptR1.at(2) - ptCNew.at(2)) * (ptCNew.at(2) - ptM.at(2)))
733 / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
734 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
735 ptR2Mirr.at(0) = 2 * (ptCNew.at(0) + t2 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
736 ptR2Mirr.at(1) = 2 * (ptCNew.at(1) + t2 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
737 ptR2Mirr.at(2) = 2 * (ptCNew.at(2) + t2 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
752 vector<Double_t> normalPMT, vector<Double_t> ptM,
753 vector<Double_t> ptR2Mirr, Double_t constantePMT)
757 Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
760 * ((normalPMT.at(0) * ptM.at(0) + normalPMT.at(1) * ptM.at(1) + normalPMT.at(2) * ptM.at(2) + constantePMT)
761 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0)) + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
762 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
763 ptPMirr.at(0) = ptM.at(0) + k1 * (ptR2Mirr.at(0) - ptM.at(0));
764 ptPMirr.at(1) = ptM.at(1) + k1 * (ptR2Mirr.at(1) - ptM.at(1));
765 ptPMirr.at(2) = ptM.at(2) + k1 * (ptR2Mirr.at(2) - ptM.at(2));
767 * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1) + normalPMT.at(2) * ptR2Mirr.at(2)
769 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0)) + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
770 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
771 ptPR2.at(0) = ptR2Mirr.at(0) + k2 * (ptR2Mirr.at(0) - ptM.at(0));
772 ptPR2.at(1) = ptR2Mirr.at(1) + k2 * (ptR2Mirr.at(1) - ptM.at(1));
773 ptPR2.at(2) = ptR2Mirr.at(2) + k2 * (ptR2Mirr.at(2) - ptM.at(2));
778 ptPMirr.at(0) * normalPMT.at(0) + ptPMirr.at(1) * normalPMT.at(1) + ptPMirr.at(2) * normalPMT.at(2) + constantePMT;
779 cout <<
"Check whether extrapolated track point on PMT plane verifies its "
780 "equation (value should be 0.):"
784 ptPR2.at(0) * normalPMT.at(0) + ptPR2.at(1) * normalPMT.at(1) + ptPR2.at(2) * normalPMT.at(2) + constantePMT;
785 cout <<
"* using reflected point R2, checkCalc = " << checkCalc2 << endl;
790 Double_t constantePMT,
string str)
792 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0, distToExtrapTrackHitInPlane = 0,
793 distToExtrapTrackHitInPlaneUnCorr = 0, distToExtrapTrackHitInPlaneIdeal = 0;
794 string histoNameX =
"", histoNameY =
"";
795 string nameX =
"", nameY =
"";
800 -1 * ((normalPMT.at(0) * ringCenter[0] + normalPMT.at(1) * ringCenter[1] + constantePMT) / normalPMT.at(2));
803 vector<Double_t> r(3),
805 r.at(0) = ringCenter[0], r.at(1) = ringCenter[1], r.at(2) = ringCenter[2];
806 p.at(0) = outPos.x(), p.at(1) = outPos.y(), p.at(2) = outPos.z();
807 cout <<
"Ring center coordinates = {" << ringCenter[0] <<
", " << ringCenter[1] <<
", " << ringCenter[2] <<
"}"
809 cout <<
"Difference in X = " << TMath::Abs(r.at(0) - p.at(0)) <<
"; \t"
810 <<
"Difference in Y = " << TMath::Abs(r.at(1) - p.at(1)) <<
"; \t"
811 <<
"Difference in Z = " << TMath::Abs(r.at(2) - p.at(2)) << endl;
813 nameX = string(
"DiffCorrX_") + str;
814 nameY = string(
"DiffCorrY_") + str;
816 if (nameX.compare(it->first) == 0) {
817 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - p.at(0)));
819 if (nameY.compare(it->first) == 0) {
820 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - p.at(1)));
824 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2)
825 + TMath::Power(r.at(2) - p.at(2), 2));
826 distToExtrapTrackHitInPlane = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2));
827 cout <<
"Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl;
828 cout <<
"Distance between fitted ring center and extrapolated track hit in "
830 << distToExtrapTrackHitInPlane << endl;
831 fHM->
H1(
"fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
832 fHM->
H1(
"fhDistanceCorrected")->Fill(distToExtrapTrackHitInPlane);
835 vector<Double_t> pUncorr(3);
836 pUncorr.at(0) = outPosUnCorr.x(), pUncorr.at(1) = outPosUnCorr.y(), pUncorr.at(2) = outPosUnCorr.z();
838 cout <<
"Difference in X w/o correction = " << TMath::Abs(r.at(0) - pUncorr.at(0)) <<
"; \t"
839 <<
"Difference in Y = " << TMath::Abs(r.at(1) - pUncorr.at(1)) <<
"; \t"
840 <<
"Difference in Z = " << TMath::Abs(r.at(2) - pUncorr.at(2)) << endl;
842 nameX = string(
"DiffUncorrX_") + str;
843 nameY = string(
"DiffUncorrY_") + str;
845 if (nameX.compare(it->first) == 0) {
846 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - pUncorr.at(0)));
848 if (nameY.compare(it->first) == 0) {
849 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - pUncorr.at(1)));
853 distToExtrapTrackHitInPlaneUnCorr =
854 TMath::Sqrt(TMath::Power(r.at(0) - pUncorr.at(0), 2) + TMath::Power(r.at(1) - pUncorr.at(1), 2));
855 fHM->
H1(
"fhDistanceUncorrected")->Fill(distToExtrapTrackHitInPlaneUnCorr);
856 cout <<
"Distance between fitted ring center and extrapolated track hit in "
858 << distToExtrapTrackHitInPlaneUnCorr << endl;
861 vector<Double_t> pIdeal(3);
862 pIdeal.at(0) = outPosIdeal.x(), pIdeal.at(1) = outPosIdeal.y(), pIdeal.at(2) = outPosIdeal.z();
864 cout <<
"Difference in X w/ ideal correction = " << TMath::Abs(r.at(0) - pIdeal.at(0)) <<
"; \t"
865 <<
"Difference in Y = " << TMath::Abs(r.at(1) - pIdeal.at(1)) <<
"; \t"
866 <<
"Difference in Z = " << TMath::Abs(r.at(2) - pIdeal.at(2)) << endl;
868 nameX = string(
"DiffIdealX_") + str;
869 nameY = string(
"DiffIdealY_") + str;
871 if (nameX.compare(it->first) == 0) {
872 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - pIdeal.at(0)));
874 if (nameY.compare(it->first) == 0) {
875 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - pIdeal.at(1)));
879 distToExtrapTrackHitInPlaneIdeal =
880 TMath::Sqrt(TMath::Power(r.at(0) - pIdeal.at(0), 2) + TMath::Power(r.at(1) - pIdeal.at(1), 2));
881 fHM->
H1(
"fhDistanceIdeal")->Fill(distToExtrapTrackHitInPlaneIdeal);
882 cout <<
"Distance between fitted ring center and extrapolated track hit in "
884 << distToExtrapTrackHitInPlaneIdeal << endl
889 if (distToExtrapTrackHitInPlane < 25. && distToExtrapTrackHitInPlaneUnCorr < 25.
890 && distToExtrapTrackHitInPlaneIdeal < 25.) {
897 cout <<
"Distance hit-ring too high!" << endl;
905 Int_t counterX = 1, counterY = 1, canX = 1500, canY = 400;
906 vector<TH1D*> histoVectX, histoVectY;
907 vector<string> stringVectX, stringVectY;
908 TString strX =
"X_mirror_tile_", strY =
"Y_mirror_tile_", str1 =
"", str2 =
"";
909 stringstream ssX, ssY, str3, str4;
910 ssX << strX << axisY <<
"_" << axisX;
911 ssY << strY << axisY <<
"_" << axisX;
912 cout <<
"ssX: " << ssX.str().c_str() <<
" and ssY: " << ssY.str().c_str() << endl << endl;
915 if (it->first.find(ssX.str().c_str()) != std::string::npos && it->second->GetEntries() >
fThreshold) {
917 cout <<
"Histo entries: " << it->second->GetEntries() << endl;
918 stringVectX.push_back(it->first);
919 histoVectX.push_back(it->second);
921 else if (it->first.find(ssY.str().c_str()) != std::string::npos && it->second->GetEntries() >
fThreshold) {
923 stringVectY.push_back(it->first);
924 histoVectY.push_back(it->second);
928 if (!histoVectX.empty()) {
929 cout <<
"Vector size: " << histoVectX.size() << endl;
930 TCanvas* c1 =
new TCanvas(ssX.str().c_str(), ssX.str().c_str(), canX, canY);
932 for (Int_t i = 0; i < 3; i++) {
934 histoVectX[i]->Draw();
935 histoVectX[i]->SetTitle(stringVectX[i].c_str());
936 histoVectX[i]->SetLineColor(counterX + 1);
937 histoVectX[i]->SetLineWidth(2);
938 histoVectX[i]->Write();
944 if (!histoVectY.empty()) {
945 TCanvas* c2 =
new TCanvas(ssY.str().c_str(), ssY.str().c_str(), canX, canY);
947 for (Int_t i = 0; i < 3; i++) {
949 histoVectY[i]->Draw();
950 histoVectY[i]->SetTitle(stringVectY[i].c_str());
951 histoVectY[i]->SetLineColor(counterY + 1);
952 histoVectY[i]->SetLineWidth(2);
953 histoVectY[i]->Write();