7#include "FairRootManager.h"
20#include "TClonesArray.h"
30#include "TGeoManager.h"
31#include "TGeoNavigator.h"
36#include "TMCProcess.h"
48 : FairTask(
"CbmRichMirrorSortingCorrection")
59 , fRefPlanePoints(NULL)
61 , fRichProjections(NULL)
63 , fRichRingMatches(NULL)
64 , fStsTrackMatches(NULL)
65 , fTrackCenterDistanceIdeal(0)
66 , fTrackCenterDistanceCorrected(0)
67 , fTrackCenterDistanceUncorrected(0)
68 , fCorrectionMatching(
"")
71 , fCorrectionTableDir(
"")
80 FairRootManager* manager = FairRootManager::Instance();
82 fGlobalTracks = (TClonesArray*) manager->GetObject(
"GlobalTrack");
84 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No GlobalTrack array!");
87 fRichRings = (TClonesArray*) manager->GetObject(
"RichRing");
89 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichRing array !");
92 fMCTracks = (TClonesArray*) manager->GetObject(
"MCTrack");
94 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No MCTracks array !");
97 fMirrorPoints = (TClonesArray*) manager->GetObject(
"RichMirrorPoint");
99 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichMirrorPoints array !");
104 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichRefPlanePoint array !");
107 fPmtPoints = (TClonesArray*) manager->GetObject(
"RichPoint");
109 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichPoint array !");
114 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichProjection array !");
117 fTrackParams = (TClonesArray*) manager->GetObject(
"RichTrackParamZ");
119 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichTrackParamZ array!");
124 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No RichRingMatch array!");
129 Fatal(
"CbmRichMirrorSortingCorrection::Init",
"No StsTrackMatch array!");
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;
411 if (trackMotherId == -1) {
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;
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();
965 int counter1 = 1, counter2 = 1, counter3 = 1, counter4 = 1, counter5 = 1, counter6 = 1, counter7 = 1, counter8 = 1,
966 counter9 = 1, counter10 = 1, counter11 = 1, counter12 = 1;
980 cout << endl <<
"CALLING FUNCTION 'DRAWMAP()' ... " << endl << endl;
981 for (
Int_t j = 0; j < 4; j++) {
982 for (
Int_t i = 0; i < 10; i++) {
1162 if (stsId < 0 || richId < 0)
continue;
1164 if (stsTrackMatch == NULL)
continue;
1167 if (richRingMatch == NULL)
continue;
1170 if (NULL == ring)
continue;
1178 if (mctrack == NULL)
continue;
1181 if (isEl && stsMcTrackId == richMcTrackId) {
1185 fHM->
H3(
"fhRingTrackDistVsXYTruematch" + ss.str())->Fill(xc, yc, rtDistance);
1186 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str())->Fill(xc, rtDistance);
1187 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str())->Fill(abs(yc), rtDistance);
1188 fHM->
H3(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str())->Fill(xc, yc, rtDistanceX);
1189 fHM->
H3(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str())->Fill(xc, yc, rtDistanceY);
1190 if (rtDistance >= -10 && rtDistance <= 10) {
1191 fHM->
H1(
"fDistUncorr")->Fill(rtDistance);
1199 const FairTrackParam* pTrack,
const CbmMCTrack* mcTrack)
1203 Double_t dx = richRing->
GetCenterX() - pTrack->GetX();
1204 Double_t dy = richRing->
GetCenterY() - pTrack->GetY();
1205 Double_t dist = TMath::Sqrt(dx * dx + dy * dy);
1212 fHM->
H3(
"fhRingTrackDistVsXYTruematch" + ss.str())->Fill(xRing, yRing, dist);
1213 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str())->Fill(xRing, dist);
1214 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str())->Fill(abs(yRing), dist);
1215 fHM->
H3(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str())->Fill(xRing, yRing, dist);
1216 fHM->
H3(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str())->Fill(xRing, yRing, dist);
1217 if (dist >= -10 && dist <= 10) {
1218 fHM->
H1(
"fDistCorr")->Fill(dist);
1226 if (mctrack == NULL)
return false;
1237 TCanvas* c =
fHM->
CreateCanvas(
"fh_ring_track_distance_vs_xy_truematch" + k,
1238 "fh_ring_track_distance_vs_xy_truematch" + k, 1800, 600);
1241 DrawH3Profile(
fHM->
H3(
"fhRingTrackDistVsXYTruematch" + ss.str()),
true,
false, 0., 2.);
1244 fHM->
H2(
"fhRingTrackDistVsXTruematch" + ss.str())->GetYaxis()->SetRangeUser(0., 1.5);
1247 fHM->
H2(
"fhRingTrackDistVsYTruematch" + ss.str())->GetYaxis()->SetRangeUser(0., 1.5);
1250 Double_t range = 2.5;
1252 TCanvas* c =
fHM->
CreateCanvas(
"fh_ring_track_distance_diff_x_y_ring_vs_xy_truematch" + k,
1253 "fh_ring_track_distance_diff_x_y_ring_vs_xy_truematch" + k, 1800, 600);
1256 DrawH3Profile(
fHM->
H3(
"fhRingTrackDistDiffXRingVsXYTruematch" + ss.str()),
true,
false, -range, range);
1258 DrawH3Profile(
fHM->
H3(
"fhRingTrackDistDiffYRingVsXYTruematch" + ss.str()),
true,
false, -range, range);
1264 TCanvas* c1 =
fHM->
CreateCanvas(
"fh_distance_distribution_corr",
"fh_distance_distribution_corr", 600, 900);
1275 TDirectory* oldir = gDirectory;
1276 TFile* outFile = FairRootManager::Instance()->GetOutFile();
1277 if (outFile != NULL) {
1280 gDirectory->cd(oldir->GetPath());
1284 string str = str1 +
"/" + str2;
1285 cout << endl << endl <<
"output string: " << str << endl << endl;
1289 TString s =
fOutputDir +
"/correction_table/track_ring_distances.txt";
1291 corrFile.open(s, std::ofstream::trunc);
1292 if (corrFile.is_open()) {
1293 corrFile <<
"number of events: " <<
fEventNb << endl;
1294 corrFile << setprecision(9) <<
"Mean distance between track hit and ring center ; corrected case = "
1297 corrFile <<
"Mean distance between track hit and ring center ; uncorrected case = "
1300 corrFile <<
"Mean distance between track hit and ring center ; ideal case = "
1304 cout <<
"Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
1309 cout <<
"number of events: " <<
fEventNb << endl;
1310 cout << setprecision(9) <<
"Mean distance between track hit and ring center ; corrected case = "
1312 cout <<
"Mean distance between track hit and ring center ; uncorrected case = "
ClassImp(CbmConverterManager)
void DrawH1andFitGauss(TH1 *hist, Bool_t drawResults, Bool_t doScale, Double_t userRangeMin, Double_t userRangeMax)
void DrawH2WithProfile(TH2 *hist, Bool_t doGaussFit, Bool_t drawOnlyMean, const string &drawOpt2D, Int_t profileColor, Int_t profileLineWidth)
TH2D * DrawH3Profile(TH3 *h, Bool_t drawMean, Bool_t doGaussFit, Double_t zMin, Double_t zMax)
Helper functions for drawing 1D and 2D histograms and graphs.
Convert internal data classes to cbmroot common data classes.
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
int32_t GetStsTrackIndex() const
int32_t GetRichRingIndex() const
TCanvas * CreateCanvas(const std::string &name, const std::string &title, Int_t width, Int_t height)
Create and draw TCanvas and store pointer to it.
void SaveCanvasToImage(const std::string &outputDir, const std::string &options="png,eps")
Save all stored canvases to images.
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
void Create2(const std::string &name, const std::string &title, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY)
Helper function for creation of 2-dimensional histograms and profiles. Template argument is a real ob...
void WriteToFile()
Write all objects to current opened file.
TH3 * H3(const std::string &name) const
Return pointer to TH3 histogram.
void Create1(const std::string &name, const std::string &title, Int_t nofBins, Double_t minBin, Double_t maxBin)
Helper function for creation of 1-dimensional histograms and profiles. Template argument is a real ob...
void Create3(const std::string &name, const std::string &title, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY, Int_t nofBinsZ, Double_t minBinZ, Double_t maxBinZ)
Helper function for creation of 3-dimensional histograms and profiles. Template argument is a real ob...
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
uint32_t GetGeantProcessId() const
void GetMomentum(TVector3 &momentum) const
int32_t GetMotherId() const
int32_t GetPdgCode() const
const CbmLink & GetMatchedLink() const
static void Init()
Initialize array of RICH hits.
static void CopyHitsToRingLight(const CbmRichRing *ring1, CbmRichRingLight *ring2)
Copy hits from CbmRichRing to CbmRichRingLight.
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
static CbmRichGeoManager & GetInstance()
void GetPmtNormal(Int_t NofPMTPoints, vector< Double_t > &normalPMT, Double_t &normalCste)
void FillHistProjection(TVector3 outPosIdeal, TVector3 outPosUnCorr, TVector3 outPos, CbmRichRingLight ringLight, vector< Double_t > normalPMT, Double_t constantePMT, string str)
TClonesArray * fRefPlanePoints
void DrawMap(Int_t strX, Int_t strY)
virtual InitStatus Init()
Inherited from FairTask.
TClonesArray * fTrackParams
TClonesArray * fMirrorPoints
TString fCorrectionMatching
TString fCorrectionTableDir
CbmRichRingFitterCOP * fCopFit
TClonesArray * fGlobalTracks
CbmRichRingFitterEllipseTau * fTauFit
void DrawHistProjection()
TClonesArray * fStsTrackMatches
void FillRingTrackDistanceCorr(const CbmRichRing *richRing, const FairTrackParam *pTrack, const CbmMCTrack *mcTrack)
virtual void Exec(Option_t *option)
Inherited from FairTask.
TClonesArray * fRichRingMatches
void DrawRingTrackDistance(Int_t k)
void ComputeR2(vector< Double_t > &ptR2Center, vector< Double_t > &ptR2Mirr, vector< Double_t > ptM, vector< Double_t > ptC, vector< Double_t > ptR1, TGeoNavigator *navi, TString option, TString mirrorTileName)
TClonesArray * fRichProjections
virtual ~CbmRichMirrorSortingCorrection()
void ComputeP(vector< Double_t > &ptPMirr, vector< Double_t > &ptPR2, vector< Double_t > normalPMT, vector< Double_t > ptM, vector< Double_t > ptR2Mirr, Double_t constantePMT)
CbmRichMirrorSortingCorrection()
virtual void Finish()
Inherited from FairTask.
TClonesArray * fPmtPoints
Double_t fTrackCenterDistanceCorrected
bool IsMcPrimaryElectron(const CbmMCTrack *mctrack)
TClonesArray * fRichRings
void InitHistProjection()
std::map< string, TH1D * > fDiffHistoMap
void FillRingTrackDistance()
Double_t fTrackCenterDistanceIdeal
Double_t fTrackCenterDistanceUncorrected
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
static double GetRingTrackDistanceY(int globalTrackId)
static double GetRingTrackDistance(int globalTrackId)
static double GetRingTrackDistanceX(int globalTrackId)
int32_t GetNofWrongHits() const
int32_t GetNofTrueHits() const
void SaveCanvasAsImage(TCanvas *c, const string &dir, const string &option)