2 #include "FairLogger.h"
3 #include "FairRootManager.h"
14 #include "TClonesArray.h"
23 #include "TGeoManager.h"
24 #include "TGeoNavigator.h"
29 #include "TMCProcess.h"
40 : FairTask(
"CbmRichMirrorSortingCorrection")
51 , fRefPlanePoints(NULL)
53 , fRichProjections(NULL)
55 , fRichRingMatches(NULL)
56 , fStsTrackMatches(NULL)
57 , fTrackCenterDistanceIdeal(0)
58 , fTrackCenterDistanceCorrected(0)
59 , fTrackCenterDistanceUncorrected(0)
60 , fCorrectionMatching(
"")
63 , fCorrectionTableDir(
"")
69 FairRootManager* manager = FairRootManager::Instance();
71 fGlobalTracks = (TClonesArray*) manager->GetObject(
"GlobalTrack");
73 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No GlobalTrack array!");
76 fRichRings = (TClonesArray*) manager->GetObject(
"RichRing");
78 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichRing array !");
81 fMCTracks = (TClonesArray*) manager->GetObject(
"MCTrack");
83 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No MCTracks array !");
86 fMirrorPoints = (TClonesArray*) manager->GetObject(
"RichMirrorPoint");
88 Fatal(
"CbmRichMirrorSortingCorrection::Init",
89 "No RichMirrorPoints array !");
94 Fatal(
"CbmRichMirrorSortingCorrection::Init",
95 "No RichRefPlanePoint array !");
98 fPmtPoints = (TClonesArray*) manager->GetObject(
"RichPoint");
100 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichPoint array !");
105 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichProjection array !");
108 fTrackParams = (TClonesArray*) manager->GetObject(
"RichTrackParamZ");
110 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichTrackParamZ array!");
115 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichRingMatch array!");
120 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No StsTrackMatch array!");
137 Double_t upperScaleLimit = 6., bin = 400.;
139 fHM->
Create1<TH1D>(
"fhDistanceCenterToExtrapolatedTrack",
140 "fhDistanceCenterToExtrapolatedTrack;Distance fitted "
141 "center to extrapolated track;Number of entries",
146 "fhDistanceCorrected;Distance a [cm];A.U.",
152 "fhDifferenceX;Difference in X (fitted center - extrapolated track);A.U.",
158 "fhDifferenceY;Difference in Y (fitted center - extrapolated track);A.U.",
164 "fhDistanceUncorrected;Distance a [cm];A.U.",
169 "fhDifferenceXUncorrected",
170 "fhDifferenceXUncorrected;Difference in X uncorrected [cm];A.U.",
175 "fhDifferenceYUncorrected",
176 "fhDifferenceYUncorrected;Difference in Y uncorrected [cm];A.U.",
182 "fhDistanceIdeal;Distance a [cm];A.U.",
187 "fhDifferenceXIdeal;Difference in X ideal [cm];A.U.",
192 "fhDifferenceYIdeal;Difference in Y ideal [cm];A.U.",
198 "fHistoDiffX;Histogram difference between corrected and "
199 "ideal X positions;A.U.",
204 "fHistoDiffY;Histogram difference between corrected and "
205 "ideal Y positions;A.U.",
211 "fHistoBoA;Histogram B axis over A axis;A.U.",
219 Double_t upperScaleLimit = 6., bin = 400.;
220 TString strCorrX =
"fhDifferenceCorrectedX_mirror_tile_",
221 strCorrY =
"fhDifferenceCorrectedY_mirror_tile_",
222 strDiffCorrX =
"DiffCorrX_mirror_tile_",
223 strDiffCorrY =
"DiffCorrY_mirror_tile_";
224 stringstream ssDiffCorrX, ssDiffCorrY, ssCorrX, ssCorrY;
225 TString strUncorrX =
"fhDifferenceUncorrectedX_mirror_tile_",
226 strUncorrY =
"fhDifferenceUncorrectedY_mirror_tile_",
227 strDiffUncorrX =
"DiffUncorrX_mirror_tile_",
228 strDiffUncorrY =
"DiffUncorrY_mirror_tile_";
229 stringstream ssDiffUncorrX, ssDiffUncorrY, ssUncorrX, ssUncorrY;
230 TString strIdealX =
"fhDifferenceIdealX_mirror_tile_",
231 strIdealY =
"fhDifferenceIdealY_mirror_tile_",
232 strDiffIdealX =
"DiffIdealX_mirror_tile_",
233 strDiffIdealY =
"DiffIdealY_mirror_tile_";
234 stringstream ssDiffIdealX, ssDiffIdealY, ssIdealX, ssIdealY;
236 cout << endl <<
"Init histo: " << endl << endl;
237 for (Int_t j = 0; j < 4; j++) {
238 for (Int_t
i = 0;
i < 10;
i++) {
239 ssDiffCorrX << strDiffCorrX << j <<
"_" <<
i;
240 ssDiffCorrY << strDiffCorrY << j <<
"_" <<
i;
241 ssCorrX << strCorrX << j <<
"_" <<
i;
242 ssCorrY << strCorrY << j <<
"_" <<
i;
244 new TH1D(ssCorrX.str().c_str(),
245 ";Difference in X (fitted center - extrapolated track);A.U.",
250 new TH1D(ssCorrY.str().c_str(),
251 ";Difference in Y (fitted center - extrapolated track);A.U.",
256 ssDiffUncorrX << strDiffUncorrX << j <<
"_" <<
i;
257 ssDiffUncorrY << strDiffUncorrY << j <<
"_" <<
i;
258 ssUncorrX << strUncorrX << j <<
"_" <<
i;
259 ssUncorrY << strUncorrY << j <<
"_" <<
i;
261 new TH1D(ssUncorrX.str().c_str(),
262 ";Difference in X (fitted center - extrapolated track);A.U.",
267 new TH1D(ssUncorrY.str().c_str(),
268 ";Difference in Y (fitted center - extrapolated track);A.U.",
273 ssDiffIdealX << strDiffIdealX << j <<
"_" <<
i;
274 ssDiffIdealY << strDiffIdealY << j <<
"_" <<
i;
275 ssIdealX << strIdealX << j <<
"_" <<
i;
276 ssIdealY << strIdealY << j <<
"_" <<
i;
278 new TH1D(ssIdealX.str().c_str(),
279 ";Difference in X (fitted center - extrapolated track);A.U.",
284 new TH1D(ssIdealY.str().c_str(),
285 ";Difference in Y (fitted center - extrapolated track);A.U.",
297 ssDiffUncorrX.str(
"");
299 ssDiffUncorrY.str(
"");
301 ssDiffIdealX.str(
"");
303 ssDiffIdealY.str(
"");
309 Double_t xMin = -120., xMax = 120., nBinsX1 = 60, yMax = 200., range = 3.;
312 for (Int_t k = 1; k < 3; k++) {
314 fHM->
Create3<TH3D>(
"fhRingTrackDistVsXYTruematch" + ss.str(),
315 "fhRingTrackDistVsXYTruematch" + ss.str()
316 +
";X [cm];Y [cm];Ring-track distance [cm]",
326 fHM->
Create2<TH2D>(
"fhRingTrackDistVsXTruematch" + ss.str(),
327 "fhRingTrackDistVsXTruematch" + ss.str()
328 +
";X [cm];Ring-track distance [cm]",
335 fHM->
Create2<TH2D>(
"fhRingTrackDistVsYTruematch" + ss.str(),
336 "fhRingTrackDistVsYTruematch" + ss.str()
337 +
";Abs(Y) [cm];Ring-track distance [cm]",
344 fHM->
Create3<TH3D>(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str(),
345 "fhRingTrackDistDiffXRingVsXYTruematch" + ss.str()
346 +
";X [cm];Y [cm];X Ring-track distance [cm]",
356 fHM->
Create3<TH3D>(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str(),
357 "fhRingTrackDistDiffYRingVsXYTruematch" + ss.str()
358 +
";X [cm];Y [cm];Y Ring-track distance [cm]",
372 "fDistUncorr;Uncorrected Distance;Number of entries",
377 "fDistCorr;Corrected Distance;Number of entries",
442 cout <<
"CbmRichMirrorSortingCorrection: Event #" <<
fEventNb << endl;
443 TVector3 momentum, outPos, outPosUnCorr, outPosIdeal;
444 Double_t constantePMT = 0., trackX = 0., trackY = 0.;
445 vector<Double_t> vect(2, 0), ptM(3, 0), ptC(3, 0), ptCIdeal(3, 0), ptR1(3, 0),
446 ptR2Center(3, 0), ptR2Mirr(3, 0), ptPR2(3, 0), ptPMirr(3, 0),
448 vector<Double_t> ptR2CenterUnCorr(3, 0), ptR2CenterIdeal(3, 0),
449 ptR2MirrUnCorr(3, 0), ptR2MirrIdeal(3, 0), ptPMirrUnCorr(3, 0),
450 ptPMirrIdeal(3, 0), ptPR2UnCorr(3, 0), ptPR2Ideal(3, 0);
452 ptC.at(0) = 0., ptC.at(1) = 132.594000, ptC.at(2) = 54.267226;
453 TVector3 mirrorPoint, dirCos,
pos;
454 Double_t nx = 0., ny = 0., nz = 0.;
459 cout <<
"Nb of rings in evt = " <<
fRichRings->GetEntries() << endl << endl;
463 for (Int_t iGlobalTrack = 0; iGlobalTrack <
fGlobalTracks->GetEntriesFast();
473 cout <<
"Error richInd < 0" << endl;
478 cout <<
"Error ring == NULL!" << endl;
487 if (NULL == cbmRichTrackMatch) {
continue; }
488 cout <<
"Nof true hits = " << cbmRichTrackMatch->
GetNofTrueHits() << endl;
494 if (mcRichTrackId < 0)
continue;
496 if (mcStsTrackId != mcRichTrackId) {
497 cout <<
"Error StsTrackIndex and TrackIndex from Ring do not match!"
502 if (!mcTrack)
continue;
509 cout <<
"ring Center Coo: " << ringL.
GetCenterX() <<
", "
514 if (pTrack == NULL) {
515 cout <<
"CbmRichMirrorSortingCorrection::Exec : pTrack = NULL." << endl;
518 trackX = pTrack->GetX(), trackY = pTrack->GetY();
519 cout <<
"Track: " << trackX <<
", " << trackY << endl;
524 Int_t pdg = TMath::Abs(mcTrack->
GetPdgCode());
525 if (trackMotherId == -1) {
528 for (Int_t iMirrPt = 0; iMirrPt <
fMirrorPoints->GetEntries();
531 if (mirrPoint == 0) {
continue; }
533 if (mirrPoint->GetTrackID() == mcRichTrackId) {
break; }
535 ptM.at(0) = mirrPoint->GetX(), ptM.at(1) = mirrPoint->GetY(),
536 ptM.at(2) = mirrPoint->GetZ();
538 mirrNode = gGeoManager->FindNode(ptM.at(0), ptM.at(1), ptM.at(2));
540 string str1 = gGeoManager->GetPath(), str2 =
"mirror_tile_",
542 std::size_t found = str1.find(str2);
543 if (found != std::string::npos) {
545 Int_t end = str2.length() + 3;
546 str3 = str1.substr(found, end);
548 cout <<
"Mirror ID: " << str3 << endl;
551 TGeoNavigator* navi = gGeoManager->GetCurrentNavigator();
553 ptCIdeal.at(0) = navi->GetCurrentMatrix()->GetTranslation()[0];
554 ptCIdeal.at(1) = navi->GetCurrentMatrix()->GetTranslation()[1];
555 ptCIdeal.at(2) = navi->GetCurrentMatrix()->GetTranslation()[2];
556 cout <<
"Sphere center coordinates of the aligned mirror tile, "
558 << ptCIdeal.at(0) <<
", " << ptCIdeal.at(1) <<
", "
559 << ptCIdeal.at(2) <<
"}" << endl;
564 if (refPlanePoint->GetTrackID() == mcRichTrackId) {
break; }
566 ptR1.at(0) = refPlanePoint->GetX(),
567 ptR1.at(1) = refPlanePoint->GetY(),
568 ptR1.at(2) = refPlanePoint->GetZ();
569 cout <<
"Refl plane point coo = {" << ptR1[0] <<
", " << ptR1[1]
570 <<
", " << ptR1[2] <<
"}" << endl;
572 ptR2Center, ptR2Mirr, ptM, ptC, ptR1, navi,
"Corrected", str3);
589 ComputeP(ptPMirr, ptPR2, normalPMT, ptM, ptR2Mirr, constantePMT);
602 cout <<
"PMT points mirr coordinates before rotation = {"
603 << ptPMirr[0] <<
", " << ptPMirr[1] <<
", " << ptPMirr[2]
605 cout <<
"PMT points mirr uncorr coordinates before rotation = {"
606 << ptPMirrUnCorr[0] <<
", " << ptPMirrUnCorr[1] <<
", "
607 << ptPMirrUnCorr[2] <<
"}" << endl;
608 cout <<
"PMT points mirr ideal coordinates before rotation = {"
609 << ptPMirrIdeal[0] <<
", " << ptPMirrIdeal[1] <<
", "
610 << ptPMirrIdeal[2] <<
"}" << endl;
612 TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
615 <<
"New PMT points coordinates = {" << outPos.x() <<
", "
616 << outPos.y() <<
", " << outPos.z() <<
"}" << endl;
617 TVector3 inPosUnCorr(
618 ptPMirrUnCorr.at(0), ptPMirrUnCorr.at(1), ptPMirrUnCorr.at(2));
621 cout <<
"New mirror points coordinates = {" << outPosUnCorr.x()
622 <<
", " << outPosUnCorr.y() <<
", " << outPosUnCorr.z() <<
"}"
625 ptPMirrIdeal.at(0), ptPMirrIdeal.at(1), ptPMirrIdeal.at(2));
628 cout <<
"New mirror points coordinates = {" << outPosIdeal.x()
629 <<
", " << outPosIdeal.y() <<
", " << outPosIdeal.z() <<
"}"
633 cout <<
"pTrack [X,Y]: " << pTrack->GetX() <<
", " << pTrack->GetY()
636 pTrack->SetX(outPosUnCorr.x());
637 pTrack->SetY(outPosUnCorr.y());
639 pTrack->SetX(outPos.x());
640 pTrack->SetY(outPos.y());
642 pTrack->SetX(outPosIdeal.x());
643 pTrack->SetY(outPosIdeal.y());
654 cout <<
"pTrack [X,Y]: " << pTrack->GetX() <<
", " << pTrack->GetY()
659 cout <<
"No mirror points registered." << endl;
662 cout <<
"Not a mother particle." << endl;
667 cout <<
"CbmRichMirrorSortingCorrection::Exec No rings in event were found."
675 vector<Double_t>& normalPMT,
676 Double_t& normalCste) {
679 Int_t pmtTrackID, pmtMotherID;
680 Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0.,
682 Double_t pmtPt[] = {0., 0., 0.};
683 Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
691 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
693 pmtTrackID = pmtPoint->GetTrackID();
696 a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
700 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
702 pmtTrackID = pmtPoint->GetTrackID();
706 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2)
707 + TMath::Power(a[1] - pmtPoint->GetY(), 2)
708 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
710 b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
715 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
717 pmtTrackID = pmtPoint->GetTrackID();
721 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2)
722 + TMath::Power(a[1] - pmtPoint->GetY(), 2)
723 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
725 && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2)
726 + TMath::Power(b[1] - pmtPoint->GetY(), 2)
727 + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
729 c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
735 k = (b[0] - a[0]) / (c[0] - a[0]);
736 if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0
737 || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
738 cout <<
"Error in normal calculation, vect_AB and vect_AC are collinear."
741 buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
742 buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
743 buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
746 / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
747 + TMath::Power(buffNormZ, 2));
750 / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
751 + TMath::Power(buffNormZ, 2));
754 / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
755 + TMath::Power(buffNormZ, 2));
759 scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0])
760 + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
761 + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
766 * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY()
767 + normalPMT.at(2) * pmtPoint1->GetZ());
769 scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0])
770 + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
771 + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
774 scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0])
775 + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
776 + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
781 vector<Double_t>& ptR2Mirr,
782 vector<Double_t> ptM,
783 vector<Double_t> ptC,
784 vector<Double_t> ptR1,
787 TString mirrorTileName) {
790 vector<Double_t> normalMirr(3), ptCNew(3), ptTileCenter(3);
791 Double_t t1 = 0., t2 = 0., t3 = 0.;
793 if (option ==
"Corrected") {
796 Int_t lineCounter = 1, lineIndex = 0;
800 string fileLine =
"", strMisX =
"", strMisY =
"";
801 Double_t misX = 0., misY = 0.;
817 if (corrFile.is_open()) {
818 while (!corrFile.eof()) {
819 getline(corrFile, fileLine);
820 lineIndex = fileLine.find(mirrorTileName, 0);
821 if (lineIndex != string::npos) {
843 cout <<
"Error in CbmRichCorrection: unable to open parameter file!"
845 cout <<
"Parameter file path: " << str << endl << endl;
852 cout <<
"Correction parameters = " << misX <<
", " << misY << endl;
853 ptCNew.at(0) = ptC.at(0) + misX;
854 ptCNew.at(1) = ptC.at(1) + misY;
855 ptCNew.at(2) = ptC.at(2);
856 ptTileCenter.at(0) = navi->GetMotherMatrix()->GetTranslation()[0];
857 ptTileCenter.at(1) = navi->GetMotherMatrix()->GetTranslation()[1];
858 ptTileCenter.at(2) = navi->GetMotherMatrix()->GetTranslation()[2];
860 Double_t
x = 0.,
y = 0., z = 0., dist = 0., dist2 = 0.,
z2 = 0.;
861 x = TMath::Power(ptCNew.at(0) - ptTileCenter.at(0), 2);
862 y = TMath::Power(ptCNew.at(1) - ptTileCenter.at(1), 2);
863 z = TMath::Power(ptCNew.at(2) - ptTileCenter.at(2), 2);
864 dist = TMath::Sqrt(
x +
y + z);
865 z2 = ptTileCenter.at(2) - TMath::Sqrt(TMath::Power(300, 2) -
x -
y)
868 dist2 = TMath::Sqrt(
x +
y + TMath::Power(
z2 - ptTileCenter.at(2), 2));
871 cout <<
"Sphere center coordinates of the rotated mirror tile, after "
873 << ptCNew.at(0) <<
", " << ptCNew.at(1) <<
", " << ptCNew.at(2) <<
"}"
875 }
else if (option ==
"Uncorrected") {
878 cout <<
"Sphere center coordinates of the rotated mirror tile, without "
880 << ptCNew.at(0) <<
", " << ptCNew.at(1) <<
", " << ptCNew.at(2) <<
"}"
888 normalMirr.at(0) = (ptCNew.at(0) - ptM.at(0))
889 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
890 + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
891 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
892 normalMirr.at(1) = (ptCNew.at(1) - ptM.at(1))
893 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
894 + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
895 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
896 normalMirr.at(2) = (ptCNew.at(2) - ptM.at(2))
897 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
898 + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
899 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
902 t1 = ((ptR1.at(0) - ptM.at(0)) * (ptCNew.at(0) - ptM.at(0))
903 + (ptR1.at(1) - ptM.at(1)) * (ptCNew.at(1) - ptM.at(1))
904 + (ptR1.at(2) - ptM.at(2)) * (ptCNew.at(2) - ptM.at(2)))
905 / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
906 + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
907 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
909 2 * (ptM.at(0) + t1 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
911 2 * (ptM.at(1) + t1 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
913 2 * (ptM.at(2) + t1 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
914 t2 = ((ptR1.at(0) - ptCNew.at(0)) * (ptCNew.at(0) - ptM.at(0))
915 + (ptR1.at(1) - ptCNew.at(1)) * (ptCNew.at(1) - ptM.at(1))
916 + (ptR1.at(2) - ptCNew.at(2)) * (ptCNew.at(2) - ptM.at(2)))
917 / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
918 + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
919 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
921 2 * (ptCNew.at(0) + t2 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
923 2 * (ptCNew.at(1) + t2 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
925 2 * (ptCNew.at(2) + t2 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
940 vector<Double_t>& ptPR2,
941 vector<Double_t> normalPMT,
942 vector<Double_t> ptM,
943 vector<Double_t> ptR2Mirr,
944 Double_t constantePMT) {
947 Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
950 * ((normalPMT.at(0) * ptM.at(0) + normalPMT.at(1) * ptM.at(1)
951 + normalPMT.at(2) * ptM.at(2) + constantePMT)
952 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0))
953 + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
954 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
955 ptPMirr.at(0) = ptM.at(0) + k1 * (ptR2Mirr.at(0) - ptM.at(0));
956 ptPMirr.at(1) = ptM.at(1) + k1 * (ptR2Mirr.at(1) - ptM.at(1));
957 ptPMirr.at(2) = ptM.at(2) + k1 * (ptR2Mirr.at(2) - ptM.at(2));
959 * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1)
960 + normalPMT.at(2) * ptR2Mirr.at(2) + constantePMT)
961 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0))
962 + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
963 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
964 ptPR2.at(0) = ptR2Mirr.at(0) + k2 * (ptR2Mirr.at(0) - ptM.at(0));
965 ptPR2.at(1) = ptR2Mirr.at(1) + k2 * (ptR2Mirr.at(1) - ptM.at(1));
966 ptPR2.at(2) = ptR2Mirr.at(2) + k2 * (ptR2Mirr.at(2) - ptM.at(2));
970 checkCalc1 = ptPMirr.at(0) * normalPMT.at(0) + ptPMirr.at(1) * normalPMT.at(1)
971 + ptPMirr.at(2) * normalPMT.at(2) + constantePMT;
972 cout <<
"Check whether extrapolated track point on PMT plane verifies its "
973 "equation (value should be 0.):"
976 checkCalc2 = ptPR2.at(0) * normalPMT.at(0) + ptPR2.at(1) * normalPMT.at(1)
977 + ptPR2.at(2) * normalPMT.at(2) + constantePMT;
978 cout <<
"* using reflected point R2, checkCalc = " << checkCalc2 << endl;
982 TVector3 outPosIdeal,
983 TVector3 outPosUnCorr,
986 vector<Double_t> normalPMT,
987 Double_t constantePMT,
989 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0,
990 distToExtrapTrackHitInPlane = 0,
991 distToExtrapTrackHitInPlaneUnCorr = 0,
992 distToExtrapTrackHitInPlaneIdeal = 0;
993 string histoNameX =
"", histoNameY =
"";
994 string nameX =
"", nameY =
"";
999 * ((normalPMT.at(0) * ringCenter[0]
1000 + normalPMT.at(1) * ringCenter[1] + constantePMT)
1004 vector<Double_t> r(3),
1006 r.at(0) = ringCenter[0], r.at(1) = ringCenter[1], r.at(2) = ringCenter[2];
1007 p.at(0) = outPos.x(), p.at(1) = outPos.y(), p.at(2) = outPos.z();
1008 cout <<
"Ring center coordinates = {" << ringCenter[0] <<
", "
1009 << ringCenter[1] <<
", " << ringCenter[2] <<
"}" << endl;
1010 cout <<
"Difference in X = " << TMath::Abs(r.at(0) - p.at(0)) <<
"; \t"
1011 <<
"Difference in Y = " << TMath::Abs(r.at(1) - p.at(1)) <<
"; \t"
1012 <<
"Difference in Z = " << TMath::Abs(r.at(2) - p.at(2)) << endl;
1014 nameX = string(
"DiffCorrX_") + str;
1015 nameY = string(
"DiffCorrY_") + str;
1016 for (std::map<string, TH1D*>::iterator it =
fDiffHistoMap.begin();
1019 if (nameX.compare(it->first) == 0) {
1020 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - p.at(0)));
1022 if (nameY.compare(it->first) == 0) {
1023 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - p.at(1)));
1027 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2)
1028 + TMath::Power(r.at(1) - p.at(1), 2)
1029 + TMath::Power(r.at(2) - p.at(2), 2));
1030 distToExtrapTrackHitInPlane = TMath::Sqrt(
1031 TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2));
1032 cout <<
"Distance between fitted ring center and extrapolated track hit = "
1033 << distToExtrapTrackHit << endl;
1034 cout <<
"Distance between fitted ring center and extrapolated track hit in "
1036 << distToExtrapTrackHitInPlane << endl;
1037 fHM->
H1(
"fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
1038 fHM->
H1(
"fhDistanceCorrected")->Fill(distToExtrapTrackHitInPlane);
1041 vector<Double_t> pUncorr(
1043 pUncorr.at(0) = outPosUnCorr.x(), pUncorr.at(1) = outPosUnCorr.y(),
1044 pUncorr.at(2) = outPosUnCorr.z();
1046 cout <<
"Difference in X w/o correction = "
1047 << TMath::Abs(r.at(0) - pUncorr.at(0)) <<
"; \t"
1048 <<
"Difference in Y = " << TMath::Abs(r.at(1) - pUncorr.at(1)) <<
"; \t"
1049 <<
"Difference in Z = " << TMath::Abs(r.at(2) - pUncorr.at(2)) << endl;
1051 nameX = string(
"DiffUncorrX_") + str;
1052 nameY = string(
"DiffUncorrY_") + str;
1053 for (std::map<string, TH1D*>::iterator it =
fDiffHistoMap.begin();
1056 if (nameX.compare(it->first) == 0) {
1057 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - pUncorr.at(0)));
1059 if (nameY.compare(it->first) == 0) {
1060 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - pUncorr.at(1)));
1064 distToExtrapTrackHitInPlaneUnCorr =
1065 TMath::Sqrt(TMath::Power(r.at(0) - pUncorr.at(0), 2)
1066 + TMath::Power(r.at(1) - pUncorr.at(1), 2));
1067 fHM->
H1(
"fhDistanceUncorrected")->Fill(distToExtrapTrackHitInPlaneUnCorr);
1068 cout <<
"Distance between fitted ring center and extrapolated track hit in "
1070 << distToExtrapTrackHitInPlaneUnCorr << endl;
1073 vector<Double_t> pIdeal(
1075 pIdeal.at(0) = outPosIdeal.x(), pIdeal.at(1) = outPosIdeal.y(),
1076 pIdeal.at(2) = outPosIdeal.z();
1078 cout <<
"Difference in X w/ ideal correction = "
1079 << TMath::Abs(r.at(0) - pIdeal.at(0)) <<
"; \t"
1080 <<
"Difference in Y = " << TMath::Abs(r.at(1) - pIdeal.at(1)) <<
"; \t"
1081 <<
"Difference in Z = " << TMath::Abs(r.at(2) - pIdeal.at(2)) << endl;
1083 nameX = string(
"DiffIdealX_") + str;
1084 nameY = string(
"DiffIdealY_") + str;
1085 for (std::map<string, TH1D*>::iterator it =
fDiffHistoMap.begin();
1088 if (nameX.compare(it->first) == 0) {
1089 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - pIdeal.at(0)));
1091 if (nameY.compare(it->first) == 0) {
1092 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - pIdeal.at(1)));
1096 distToExtrapTrackHitInPlaneIdeal =
1097 TMath::Sqrt(TMath::Power(r.at(0) - pIdeal.at(0), 2)
1098 + TMath::Power(r.at(1) - pIdeal.at(1), 2));
1099 fHM->
H1(
"fhDistanceIdeal")->Fill(distToExtrapTrackHitInPlaneIdeal);
1100 cout <<
"Distance between fitted ring center and extrapolated track hit in "
1102 << distToExtrapTrackHitInPlaneIdeal << endl
1107 if (distToExtrapTrackHitInPlane < 25.
1108 && distToExtrapTrackHitInPlaneUnCorr < 25.
1109 && distToExtrapTrackHitInPlaneIdeal < 25.) {
1115 cout <<
"Distance hit-ring too high!" << endl;
1122 Int_t counterX = 1, counterY = 1, canX = 1500, canY = 400;
1123 vector<TH1D*> histoVectX, histoVectY;
1124 vector<string> stringVectX, stringVectY;
1125 TString strX =
"X_mirror_tile_", strY =
"Y_mirror_tile_", str1 =
"",
1127 stringstream ssX, ssY, str3, str4;
1128 ssX << strX << axisY <<
"_" << axisX;
1129 ssY << strY << axisY <<
"_" << axisX;
1130 cout <<
"ssX: " << ssX.str().c_str() <<
" and ssY: " << ssY.str().c_str()
1134 for (std::map<string, TH1D*>::iterator it =
fDiffHistoMap.begin();
1137 if (it->first.find(ssX.str().c_str()) != std::string::npos
1140 cout <<
"Histo entries: " << it->second->GetEntries() << endl;
1141 stringVectX.push_back(it->first);
1142 histoVectX.push_back(it->second);
1143 }
else if (it->first.find(ssY.str().c_str()) != std::string::npos
1146 stringVectY.push_back(it->first);
1147 histoVectY.push_back(it->second);
1151 if (!histoVectX.empty()) {
1152 cout <<
"Vector size: " << histoVectX.size() << endl;
1153 TCanvas* c1 =
new TCanvas(ssX.str().c_str(), ssX.str().c_str(), canX, canY);
1155 for (Int_t
i = 0;
i < 3;
i++) {
1157 histoVectX[
i]->Draw();
1158 histoVectX[
i]->SetTitle(stringVectX[
i].c_str());
1159 histoVectX[
i]->SetLineColor(counterX + 1);
1160 histoVectX[
i]->SetLineWidth(2);
1161 histoVectX[
i]->Write();
1167 if (!histoVectY.empty()) {
1168 TCanvas* c2 =
new TCanvas(ssY.str().c_str(), ssY.str().c_str(), canX, canY);
1170 for (Int_t
i = 0;
i < 3;
i++) {
1172 histoVectY[
i]->Draw();
1173 histoVectY[
i]->SetTitle(stringVectY[
i].c_str());
1174 histoVectY[
i]->SetLineColor(counterY + 1);
1175 histoVectY[
i]->SetLineWidth(2);
1176 histoVectY[
i]->Write();
1187 int counter1 = 1, counter2 = 1, counter3 = 1, counter4 = 1, counter5 = 1,
1188 counter6 = 1, counter7 = 1, counter8 = 1, counter9 = 1, counter10 = 1,
1189 counter11 = 1, counter12 = 1;
1203 cout << endl <<
"CALLING FUNCTION 'DRAWMAP()' ... " << endl << endl;
1204 for (Int_t j = 0; j < 4; j++) {
1205 for (Int_t
i = 0;
i < 10;
i++) {
1380 for (Int_t iTrack = 0; iTrack <
fGlobalTracks->GetEntriesFast(); iTrack++) {
1385 if (stsId < 0 || richId < 0)
continue;
1388 if (stsTrackMatch == NULL)
continue;
1392 if (richRingMatch == NULL)
continue;
1396 if (NULL == ring)
continue;
1404 if (mctrack == NULL)
continue;
1407 if (isEl && stsMcTrackId == richMcTrackId) {
1411 fHM->
H3(
"fhRingTrackDistVsXYTruematch" + ss.str())
1412 ->Fill(xc, yc, rtDistance);
1413 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str())->Fill(xc, rtDistance);
1414 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str())
1415 ->Fill(abs(yc), rtDistance);
1416 fHM->
H3(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str())
1417 ->Fill(xc, yc, rtDistanceX);
1418 fHM->
H3(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str())
1419 ->Fill(xc, yc, rtDistanceY);
1420 if (rtDistance >= -10 && rtDistance <= 10) {
1421 fHM->
H1(
"fDistUncorr")->Fill(rtDistance);
1430 const FairTrackParam* pTrack,
1434 Double_t dx = richRing->
GetCenterX() - pTrack->GetX();
1435 Double_t dy = richRing->
GetCenterY() - pTrack->GetY();
1436 Double_t dist = TMath::Sqrt(dx * dx + dy * dy);
1443 fHM->
H3(
"fhRingTrackDistVsXYTruematch" + ss.str())
1444 ->Fill(xRing, yRing, dist);
1445 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str())->Fill(xRing, dist);
1446 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str())->Fill(abs(yRing), dist);
1447 fHM->
H3(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str())
1448 ->Fill(xRing, yRing, dist);
1449 fHM->
H3(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str())
1450 ->Fill(xRing, yRing, dist);
1451 if (dist >= -10 && dist <= 10) {
fHM->
H1(
"fDistCorr")->Fill(dist); }
1458 if (mctrack == NULL)
return false;
1468 TCanvas* c =
fHM->
CreateCanvas(
"fh_ring_track_distance_vs_xy_truematch" + k,
1469 "fh_ring_track_distance_vs_xy_truematch" + k,
1475 fHM->
H3(
"fhRingTrackDistVsXYTruematch" + ss.str()),
true,
false, 0., 2.);
1478 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str()),
false,
true);
1479 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str())
1481 ->SetRangeUser(0., 1.5);
1484 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str()),
false,
true);
1485 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str())
1487 ->SetRangeUser(0., 1.5);
1490 Double_t range = 2.5;
1493 "fh_ring_track_distance_diff_x_y_ring_vs_xy_truematch" + k,
1494 "fh_ring_track_distance_diff_x_y_ring_vs_xy_truematch" + k,
1515 "fh_distance_distribution_corr",
"fh_distance_distribution_corr", 600, 900);
1525 TDirectory* oldir = gDirectory;
1526 TFile* outFile = FairRootManager::Instance()->GetOutFile();
1528 gDirectory->cd(oldir->GetPath());
1532 string str = str1 +
"/" + str2;
1533 cout << endl << endl <<
"output string: " << str << endl << endl;
1537 TString s =
fOutputDir +
"/correction_table/track_ring_distances.txt";
1539 corrFile.open(s, std::ofstream::trunc);
1540 if (corrFile.is_open()) {
1541 corrFile <<
"number of events: " <<
fEventNb << endl;
1544 <<
"Mean distance between track hit and ring center ; corrected case = "
1548 <<
"Mean distance between track hit and ring center ; uncorrected case = "
1552 <<
"Mean distance between track hit and ring center ; ideal case = "
1556 cout <<
"Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
1561 cout <<
"number of events: " <<
fEventNb << endl;
1562 cout << setprecision(9)
1563 <<
"Mean distance between track hit and ring center ; corrected case = "
1567 <<
"Mean distance between track hit and ring center ; uncorrected case = "
1570 cout <<
"Mean distance between track hit and ring center ; ideal case = "