CbmRoot
Loading...
Searching...
No Matches
CbmRichMirrorSortingCorrection.cxx
Go to the documentation of this file.
1/* Copyright (C) 2016-2021 Justus-Liebig-Universitaet Giessen, Giessen
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Jordan Bendarouach [committer] */
4
6
7#include "FairRootManager.h"
9
10#include <Logger.h>
11
12// ----- PART 1 ----- //
13#include "CbmGlobalTrack.h"
14#include "CbmRichConverter.h"
17#include "CbmRichUtil.h"
18#include "CbmUtils.h"
19#include "TCanvas.h"
20#include "TClonesArray.h"
21#include "TH2D.h"
22#include "TLegend.h"
23#include "TVector3.h"
24
25#include <vector>
26
27// ----- PART 2 ----- //
28#include "CbmRichGeoManager.h"
29#include "CbmRichPoint.h"
30#include "TGeoManager.h"
31#include "TGeoNavigator.h"
32class TGeoNode;
33class TGeoMatrix;
34#include "CbmDrawHist.h"
35#include "CbmTrackMatchNew.h"
36#include "TMCProcess.h"
37#include "TStyle.h"
38#include "string.h"
39
40#include <fstream>
41#include <iostream>
42#include <sstream>
43
44using namespace std;
45
46
48 : FairTask("CbmRichMirrorSortingCorrection")
49 , fEventNb(0)
50 , fHM(NULL)
51 , fHM2(NULL)
52 , fDiffHistoMap()
53 , fCopFit(NULL)
54 , fTauFit(NULL)
55 , fGlobalTracks(NULL)
56 , fRichRings(NULL)
57 , fMCTracks(NULL)
58 , fMirrorPoints(NULL)
59 , fRefPlanePoints(NULL)
60 , fPmtPoints(NULL)
61 , fRichProjections(NULL)
62 , fTrackParams(NULL)
63 , fRichRingMatches(NULL)
64 , fStsTrackMatches(NULL)
65 , fTrackCenterDistanceIdeal(0)
66 , fTrackCenterDistanceCorrected(0)
67 , fTrackCenterDistanceUncorrected(0)
68 , fCorrectionMatching("")
69 , fThreshold(0)
70 , fOutputDir("")
71 , fCorrectionTableDir("")
72 , fStudyName("")
73{
74}
75
77
79{
80 FairRootManager* manager = FairRootManager::Instance();
81
82 fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
83 if (NULL == fGlobalTracks) {
84 Fatal("CbmRichMirrorSortingCorrection::Init", "No GlobalTrack array!");
85 }
86
87 fRichRings = (TClonesArray*) manager->GetObject("RichRing");
88 if (NULL == fRichRings) {
89 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichRing array !");
90 }
91
92 fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
93 if (NULL == fMCTracks) {
94 Fatal("CbmRichMirrorSortingCorrection::Init", "No MCTracks array !");
95 }
96
97 fMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
98 if (NULL == fMirrorPoints) {
99 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichMirrorPoints array !");
100 }
101
102 fRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
103 if (NULL == fRefPlanePoints) {
104 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichRefPlanePoint array !");
105 }
106
107 fPmtPoints = (TClonesArray*) manager->GetObject("RichPoint");
108 if (NULL == fPmtPoints) {
109 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichPoint array !");
110 }
111
112 fRichProjections = (TClonesArray*) manager->GetObject("RichProjection");
113 if (NULL == fRichProjections) {
114 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichProjection array !");
115 }
116
117 fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
118 if (NULL == fTrackParams) {
119 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichTrackParamZ array!");
120 }
121
122 fRichRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
123 if (NULL == fRichRingMatches) {
124 Fatal("CbmRichMirrorSortingCorrection::Init", "No RichRingMatch array!");
125 }
126
127 fStsTrackMatches = (TClonesArray*) manager->GetObject("StsTrackMatch");
128 if (NULL == fStsTrackMatches) {
129 Fatal("CbmRichMirrorSortingCorrection::Init", "No StsTrackMatch array!");
130 }
131
135
137
138 InitHistoMap();
139
140 return kSUCCESS;
141}
142
144{
145 fHM = new CbmHistManager();
146
147 Double_t upperScaleLimit = 6., bin = 400.;
148 // fhDistance => fhDistanceCenterToExtrapolatedTrack.
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);
158
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);
164
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,
167 upperScaleLimit);
168 fHM->Create1<TH1D>("fhDifferenceYIdeal", "fhDifferenceYIdeal;Difference in Y ideal [cm];A.U.", bin, -upperScaleLimit,
169 upperScaleLimit);
170
171 fHM->Create1<TH1D>("fHistoDiffX",
172 "fHistoDiffX;Histogram difference between corrected and "
173 "ideal X positions;A.U.",
174 bin, 0., upperScaleLimit);
175 fHM->Create1<TH1D>("fHistoDiffY",
176 "fHistoDiffY;Histogram difference between corrected and "
177 "ideal Y positions;A.U.",
178 bin, 0., upperScaleLimit);
179
180 fHM->Create1<TH1D>("fHistoBoA", "fHistoBoA;Histogram B axis over A axis;A.U.", bin, 0., upperScaleLimit);
181}
182
184{
185 //fHM2 = new CbmHistManager();
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;
196
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;
204 fDiffHistoMap[ssDiffCorrX.str().c_str()] = new TH1D(
205 ssCorrX.str().c_str(), ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
206 fDiffHistoMap[ssDiffCorrY.str().c_str()] = new TH1D(
207 ssCorrY.str().c_str(), ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
208
209 ssDiffUncorrX << strDiffUncorrX << j << "_" << i;
210 ssDiffUncorrY << strDiffUncorrY << j << "_" << i;
211 ssUncorrX << strUncorrX << j << "_" << i;
212 ssUncorrY << strUncorrY << j << "_" << i;
213 fDiffHistoMap[ssDiffUncorrX.str().c_str()] =
214 new TH1D(ssUncorrX.str().c_str(), ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0.,
215 upperScaleLimit);
216 fDiffHistoMap[ssDiffUncorrY.str().c_str()] =
217 new TH1D(ssUncorrY.str().c_str(), ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0.,
218 upperScaleLimit);
219
220 ssDiffIdealX << strDiffIdealX << j << "_" << i;
221 ssDiffIdealY << strDiffIdealY << j << "_" << i;
222 ssIdealX << strIdealX << j << "_" << i;
223 ssIdealY << strIdealY << j << "_" << i;
224 fDiffHistoMap[ssDiffIdealX.str().c_str()] = new TH1D(
225 ssIdealX.str().c_str(), ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
226 fDiffHistoMap[ssDiffIdealY.str().c_str()] = new TH1D(
227 ssIdealY.str().c_str(), ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
228
229 // cout << "CorrX: " << ssDiffCorrX.str().c_str() << ", " << ssCorrX.str().c_str() << endl << "UncorrX: " << ssDiffUncorrX.str().c_str() << ", " <<
230 // ssUncorrX.str().c_str() << endl << "IdealX: " << ssDiffIdealX.str().c_str() << ", " << ssIdealX.str().c_str() << endl << endl;
231
232 ssDiffCorrX.str("");
233 ssCorrX.str("");
234 ssDiffCorrY.str("");
235 ssCorrY.str("");
236 ssDiffUncorrX.str("");
237 ssUncorrX.str("");
238 ssDiffUncorrY.str("");
239 ssUncorrY.str("");
240 ssDiffIdealX.str("");
241 ssIdealX.str("");
242 ssDiffIdealY.str("");
243 ssIdealY.str("");
244 }
245 }
246
247 //fHM = new CbmHistManager();
248 Double_t xMin = -120., xMax = 120., nBinsX1 = 60, yMax = 200., range = 3.;
249 stringstream ss;
250
251 for (Int_t k = 1; k < 3; k++) {
252 ss << 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,
258 xMax, 100, 0., 5.);
259 fHM->Create2<TH2D>("fhRingTrackDistVsYTruematch" + ss.str(),
260 "fhRingTrackDistVsYTruematch" + ss.str() + ";Abs(Y) [cm];Ring-track distance [cm]", 34, 110.,
261 yMax, 100, 0., 5.);
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);
268 ss.str("");
269 }
270
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);
273
274 /*
275// Double_t upperScaleLimit = 6., bin = 400.;
276 fDiffHistoMap["DiffCorrX_mirror_tile_2_8"] = new TH1D("fhDifferenceCorrectedX_mirror_tile_2_8", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
277 fDiffHistoMap["DiffCorrY_mirror_tile_2_8"] = new TH1D("fhDifferenceCorrectedY_mirror_tile_2_8", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
278 fDiffHistoMap["DiffUncorrX_mirror_tile_2_8"] = new TH1D("fhDifferenceUncorrectedX_mirror_tile_2_8", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
279 fDiffHistoMap["DiffUncorrY_mirror_tile_2_8"] = new TH1D("fhDifferenceUncorrectedY_mirror_tile_2_8", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
280 fDiffHistoMap["DiffIdealX_mirror_tile_2_8"] = new TH1D("fhDifferenceIdealX_mirror_tile_2_8", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
281 fDiffHistoMap["DiffIdealY_mirror_tile_2_8"] = new TH1D("fhDifferenceIdealY_mirror_tile_2_8", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
282
283 fDiffHistoMap["DiffCorrX_mirror_tile_1_3"] = new TH1D("fhDifferenceCorrectedX_mirror_tile_1_3", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
284 fDiffHistoMap["DiffCorrY_mirror_tile_1_3"] = new TH1D("fhDifferenceCorrectedY_mirror_tile_1_3", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
285 fDiffHistoMap["DiffUncorrX_mirror_tile_1_3"] = new TH1D("fhDifferenceUncorrectedX_mirror_tile_1_3", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
286 fDiffHistoMap["DiffUncorrY_mirror_tile_1_3"] = new TH1D("fhDifferenceUncorrectedY_mirror_tile_1_3", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
287 fDiffHistoMap["DiffIdealX_mirror_tile_1_3"] = new TH1D("fhDifferenceIdealX_mirror_tile_1_3", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
288 fDiffHistoMap["DiffIdealY_mirror_tile_1_3"] = new TH1D("fhDifferenceIdealY_mirror_tile_1_3", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
289
290 fDiffHistoMap["DiffCorrX_mirror_tile_1_4"] = new TH1D("fhDifferenceCorrectedX_mirror_tile_1_4", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
291 fDiffHistoMap["DiffCorrY_mirror_tile_1_4"] = new TH1D("fhDifferenceCorrectedY_mirror_tile_1_4", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
292 fDiffHistoMap["DiffUncorrX_mirror_tile_1_4"] = new TH1D("fhDifferenceUncorrectedX_mirror_tile_1_4", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
293 fDiffHistoMap["DiffUncorrY_mirror_tile_1_4"] = new TH1D("fhDifferenceUncorrectedY_mirror_tile_1_4", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
294 fDiffHistoMap["DiffIdealX_mirror_tile_1_4"] = new TH1D("fhDifferenceIdealX_mirror_tile_1_4", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
295 fDiffHistoMap["DiffIdealY_mirror_tile_1_4"] = new TH1D("fhDifferenceIdealY_mirror_tile_1_4", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
296
297 fDiffHistoMap["DiffCorrX_mirror_tile_0_8"] = new TH1D("fhDifferenceCorrectedX_mirror_tile_0_8", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
298 fDiffHistoMap["DiffCorrY_mirror_tile_0_8"] = new TH1D("fhDifferenceCorrectedY_mirror_tile_0_8", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
299 fDiffHistoMap["DiffUncorrX_mirror_tile_0_8"] = new TH1D("fhDifferenceUncorrectedX_mirror_tile_0_8", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
300 fDiffHistoMap["DiffUncorrY_mirror_tile_0_8"] = new TH1D("fhDifferenceUncorrectedY_mirror_tile_0_8", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
301 fDiffHistoMap["DiffIdealX_mirror_tile_0_8"] = new TH1D("fhDifferenceIdealX_mirror_tile_0_8", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
302 fDiffHistoMap["DiffIdealY_mirror_tile_0_8"] = new TH1D("fhDifferenceIdealY_mirror_tile_0_8", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
303
304 fDiffHistoMap["DiffCorrX_mirror_tile_1_5"] = new TH1D("fhDifferenceCorrectedX_mirror_tile_1_5", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
305 fDiffHistoMap["DiffCorrY_mirror_tile_1_5"] = new TH1D("fhDifferenceCorrectedY_mirror_tile_1_5", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
306 fDiffHistoMap["DiffUncorrX_mirror_tile_1_5"] = new TH1D("fhDifferenceUncorrectedX_mirror_tile_1_5", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
307 fDiffHistoMap["DiffUncorrY_mirror_tile_1_5"] = new TH1D("fhDifferenceUncorrectedY_mirror_tile_1_5", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
308 fDiffHistoMap["DiffIdealX_mirror_tile_1_5"] = new TH1D("fhDifferenceIdealX_mirror_tile_1_5", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
309 fDiffHistoMap["DiffIdealY_mirror_tile_1_5"] = new TH1D("fhDifferenceIdealY_mirror_tile_1_5", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
310
311 fDiffHistoMap["DiffCorrX_mirror_tile_0_1"] = new TH1D("fhDifferenceCorrectedX_mirror_tile_0_1", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
312 fDiffHistoMap["DiffCorrY_mirror_tile_0_1"] = new TH1D("fhDifferenceCorrectedY_mirror_tile_0_1", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
313 fDiffHistoMap["DiffUncorrX_mirror_tile_0_1"] = new TH1D("fhDifferenceUncorrectedX_mirror_tile_0_1", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
314 fDiffHistoMap["DiffUncorrY_mirror_tile_0_1"] = new TH1D("fhDifferenceUncorrectedY_mirror_tile_0_1", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
315 fDiffHistoMap["DiffIdealX_mirror_tile_0_1"] = new TH1D("fhDifferenceIdealX_mirror_tile_0_1", ";Difference in X (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
316 fDiffHistoMap["DiffIdealY_mirror_tile_0_1"] = new TH1D("fhDifferenceIdealY_mirror_tile_0_1", ";Difference in Y (fitted center - extrapolated track);A.U.", bin, 0., upperScaleLimit);
317*/
318
319 /*fDiffHistoMap["fhDifferenceX_mirror_tile_2_8"] = "X_mirror_tile_2_8";
320 fDiffHistoMap["fhDifferenceY_mirror_tile_2_8"] = "Y_mirror_tile_2_8";
321 fDiffHistoMap["fhDifferenceX_mirror_tile_1_3"] = "X_mirror_tile_1_3";
322 fDiffHistoMap["fhDifferenceY_mirror_tile_1_3"] = "Y_mirror_tile_1_3";
323 fDiffHistoMap["fhDifferenceX_mirror_tile_1_4"] = "X_mirror_tile_1_4";
324 fDiffHistoMap["fhDifferenceY_mirror_tile_1_4"] = "Y_mirror_tile_1_4";
325
326 for (std::map<string,string>::iterator it=fDiffHistoMap.begin(); it!=fDiffHistoMap.end(); ++it) { // Initialize all the histograms, using map IDs as inputs.
327 cout << "first: " << it->first << " and second: " << it->second << endl;
328 fHM2->Create1<TH1D>(it->first, it->first + ";Difference in X (fitted center - extrapolated track);A.U.", bin, -upperScaleLimit, upperScaleLimit);
329 }*/
330}
331
333{
334 fEventNb++;
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);
342 ptC[0] = 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.;
346 TGeoNode* mirrNode;
347 CbmRichPoint *mirrPoint, *refPlanePoint;
348
349 if (fRichRings->GetEntriesFast() != 0) {
350 cout << "Nb of rings in evt = " << fRichRings->GetEntriesFast() << endl << endl;
351 GetPmtNormal(fPmtPoints->GetEntriesFast(), normalPMT, constantePMT);
352 //cout << "Calculated normal vector to PMT plane = {" << normalPMT.at(0) << ", " << normalPMT.at(1) << ", " << normalPMT.at(2) << "} and constante d = " << constantePMT << endl;
353
354 for (Int_t iGlobalTrack = 0; iGlobalTrack < fGlobalTracks->GetEntriesFast(); iGlobalTrack++) {
355 // ----- PART 1 ----- //
356 // Ring-Track matching + Ring fit + Track momentum:
357 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack);
358 Int_t richInd = gTrack->GetRichRingIndex();
359 Int_t stsInd = gTrack->GetStsTrackIndex();
360 //cout << "richInd: " << richInd << endl;
361 if (richInd < 0) {
362 cout << "Error richInd < 0" << endl;
363 continue;
364 }
365 CbmRichRing* ring = (CbmRichRing*) fRichRings->At(richInd);
366 if (ring == NULL) {
367 cout << "Error ring == NULL!" << endl;
368 continue;
369 }
370 //Int_t ringTrackID = ring->GetTrackID(); //Old code not working with changes for new matching method.
371 //cout << "ringTrackID: " << ringTrackID << endl;
372 CbmTrackMatchNew* cbmRichTrackMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
373 CbmTrackMatchNew* cbmStsTrackMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
374 if (NULL == cbmRichTrackMatch) {
375 continue;
376 }
377 cout << "Nof true hits = " << cbmRichTrackMatch->GetNofTrueHits() << endl;
378 cout << "Nof wrong hits = " << cbmRichTrackMatch->GetNofWrongHits() << endl;
379 Int_t mcRichTrackId = cbmRichTrackMatch->GetMatchedLink().GetIndex();
380 Int_t mcStsTrackId = cbmStsTrackMatch->GetMatchedLink().GetIndex();
381 //cout << "mcTrackId: " << mcRichTrackId << endl;
382 if (mcRichTrackId < 0) continue;
383 //if (mcStsTrackId != ringTrackID) { //Old code not working with changes for new matching method.
384 if (mcStsTrackId != mcRichTrackId) {
385 cout << "Error StsTrackIndex and TrackIndex from Ring do not match!" << endl;
386 continue;
387 }
388 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcRichTrackId);
389 if (!mcTrack) continue;
390 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMCTracks->At(mcStsTrackId);
391
392 CbmRichRingLight ringL;
394 fCopFit->DoFit(&ringL);
395 //fTauFit->DoFit(&ringL);
396 cout << "ring Center Coo: " << ringL.GetCenterX() << ", " << ringL.GetCenterY() << endl;
397 mcTrack->GetMomentum(momentum);
398 //FairTrackParam* pTrack = (FairTrackParam*) fRichProjections->At(ringTrackID); //Old code not working with changes for new matching method.
399 FairTrackParam* pTrack = (FairTrackParam*) fRichProjections->At(stsInd);
400 if (pTrack == NULL) {
401 cout << "CbmRichMirrorSortingCorrection::Exec : pTrack = NULL." << endl;
402 continue;
403 }
404 trackX = pTrack->GetX(), trackY = pTrack->GetY();
405 cout << "Track: " << trackX << ", " << trackY << endl;
406
407 // ----- PART 2 ----- //
408 // Mirror ID via TGeoNavigator + Extrap hit:
409 Int_t trackMotherId = mcTrack->GetMotherId();
410 Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
411 if (trackMotherId == -1) {
412 if (fMirrorPoints->GetEntriesFast() > 0) {
413 //loop on mirrorPoint and compare w/ TrackID->GetTrackId to get correct one
414 for (Int_t iMirrPt = 0; iMirrPt < fMirrorPoints->GetEntriesFast(); iMirrPt++) {
415 mirrPoint = (CbmRichPoint*) fMirrorPoints->At(iMirrPt);
416 if (mirrPoint == 0) {
417 continue;
418 }
419 //cout << "Mirror point track ID: " << mirrPoint->GetTrackID() << endl;
420 if (mirrPoint->GetTrackID() == mcRichTrackId) {
421 break;
422 }
423 }
424 ptM.at(0) = mirrPoint->GetX(), ptM.at(1) = mirrPoint->GetY(), ptM.at(2) = mirrPoint->GetZ();
425 //cout << "mirrPoint: {" << mirrPoint->GetX() << ", " << mirrPoint->GetY() << ", " << mirrPoint->GetZ() << "}" << endl;
426 mirrNode = gGeoManager->FindNode(ptM.at(0), ptM.at(1), ptM.at(2));
427 //cout << "Mirror node name: " << mirrNode->GetName() << " and full path " << gGeoManager->GetPath() << endl;
428 string str1 = gGeoManager->GetPath(), str2 = "mirror_tile_", str3 = "";
429 std::size_t found = str1.find(str2);
430 if (found != std::string::npos) {
431 //cout << "first 'mirror_tile_type' found at: " << found << '\n';
432 Int_t end = str2.length() + 3;
433 str3 = str1.substr(found, end);
434 }
435 cout << "Mirror ID: " << str3 << endl;
436
437 if (mirrNode) {
438 TGeoNavigator* navi = gGeoManager->GetCurrentNavigator();
439 //cout << "Navigator path: " << navi->GetPath() << endl;
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, "
444 "ideal = {"
445 << ptCIdeal.at(0) << ", " << ptCIdeal.at(1) << ", " << ptCIdeal.at(2) << "}" << endl;
446 for (Int_t iRefPt = 0; iRefPt < fRefPlanePoints->GetEntriesFast(); iRefPt++) {
447 refPlanePoint = (CbmRichPoint*) fRefPlanePoints->At(iRefPt);
448 //cout << "Refl plane point track ID: " << refPlanePoint->GetTrackID() << endl;
449 if (refPlanePoint->GetTrackID() == mcRichTrackId) {
450 break;
451 }
452 }
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;
467
468 TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
470 cout << endl
471 << "New PMT points coordinates = {" << outPos.x() << ", " << outPos.y() << ", " << outPos.z() << "}"
472 << endl;
473 TVector3 inPosUnCorr(ptPMirrUnCorr.at(0), ptPMirrUnCorr.at(1), ptPMirrUnCorr.at(2));
474 CbmRichGeoManager::GetInstance().RotatePoint(&inPosUnCorr, &outPosUnCorr);
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));
478 CbmRichGeoManager::GetInstance().RotatePoint(&inPosIdeal, &outPosIdeal);
479 cout << "New mirror points coordinates = {" << outPosIdeal.x() << ", " << outPosIdeal.y() << ", "
480 << outPosIdeal.z() << "}" << endl
481 << endl;
482
483 cout << "pTrack [X,Y]: " << pTrack->GetX() << ", " << pTrack->GetY() << endl;
484 if (fCorrectionMatching == "standard") {
485 pTrack->SetX(outPosUnCorr.x());
486 pTrack->SetY(outPosUnCorr.y());
487 }
488 else if (fCorrectionMatching == "correction") {
489 pTrack->SetX(outPos.x());
490 pTrack->SetY(outPos.y());
491 }
492 else if (fCorrectionMatching == "ideal") {
493 pTrack->SetX(outPosIdeal.x());
494 pTrack->SetY(outPosIdeal.y());
495 }
496
497 FillHistProjection(outPosIdeal, outPosUnCorr, outPos, ringL, normalPMT, constantePMT, str3);
498
499 cout << "pTrack [X,Y]: " << pTrack->GetX() << ", " << pTrack->GetY() << endl;
500 FillRingTrackDistanceCorr(ring, pTrack, mcTrack2);
501 }
502 }
503 else {
504 cout << "No mirror points registered." << endl;
505 }
506 }
507 else {
508 cout << "Not a mother particle." << endl;
509 }
510 //ComputeAngles();
511 }
512 }
513 else {
514 cout << "CbmRichMirrorSortingCorrection::Exec No rings in event were found." << endl;
515 }
516
518}
519
520void CbmRichMirrorSortingCorrection::GetPmtNormal(Int_t NofPMTPoints, vector<Double_t>& normalPMT, Double_t& normalCste)
521{
522 //cout << endl << "//------------------------------ CbmRichMirrorSortingAlignment: Calculate PMT Normal ------------------------------//" << endl << endl;
523
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.};
528 CbmMCTrack* track;
529
530 /*
531 * Selection of three points (A, B, C), which form a plan and from which the calculation of the normal of the plan can be computed.
532 * Formula used is: vect(AB) x vect(AC) = normal.
533 * Normalize the normal vector obtained and check with three random points from the PMT plane, whether the scalar product is equal to zero.
534 */
535 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
536 CbmRichPoint* pmtPoint = (CbmRichPoint*) fPmtPoints->At(iPmt);
537 pmtTrackID = pmtPoint->GetTrackID();
538 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
539 pmtMotherID = track->GetMotherId();
540 a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
541 //cout << "a[0] = " << a[0] << ", a[1] = " << a[1] << " et a[2] = " << a[2] << endl;
542 break;
543 }
544 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
545 CbmRichPoint* pmtPoint = (CbmRichPoint*) fPmtPoints->At(iPmt);
546 pmtTrackID = pmtPoint->GetTrackID();
547 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
548 pmtMotherID = track->GetMotherId();
549 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
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))
552 > 7) {
553 b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
554 //cout << "b[0] = " << b[0] << ", b[1] = " << b[1] << " et b[2] = " << b[2] << endl;
555 break;
556 }
557 }
558 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
559 CbmRichPoint* pmtPoint = (CbmRichPoint*) fPmtPoints->At(iPmt);
560 pmtTrackID = pmtPoint->GetTrackID();
561 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
562 pmtMotherID = track->GetMotherId();
563 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
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))
566 > 7
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))
569 > 7) {
570 c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
571 //cout << "c[0] = " << c[0] << ", c[1] = " << c[1] << " et c[2] = " << c[2] << endl;
572 break;
573 }
574 }
575
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;
579 }
580 else {
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]);
584 normalPMT.at(0) =
585 buffNormX / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
586 normalPMT.at(1) =
587 buffNormY / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
588 normalPMT.at(2) =
589 buffNormZ / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
590 }
591
592 CbmRichPoint* pmtPoint1 = (CbmRichPoint*) fPmtPoints->At(20);
593 scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
594 + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
595 //cout << "1st scalar product between vectAM and normale = " << scalarProd << endl;
596 // To determine the constant term of the plane equation, inject the coordinates of a pmt point, which should solve it: a*x+b*y+c*z+d=0.
597 normalCste =
598 -1
599 * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY() + normalPMT.at(2) * pmtPoint1->GetZ());
600 CbmRichPoint* pmtPoint2 = (CbmRichPoint*) fPmtPoints->At(15);
601 scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
602 + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
603 //cout << "2nd scalar product between vectAM and normale = " << scalarProd << endl;
604 CbmRichPoint* pmtPoint3 = (CbmRichPoint*) fPmtPoints->At(25);
605 scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
606 + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
607 //cout << "3nd scalar product between vectAM and normale = " << scalarProd << endl;
608}
609
610void CbmRichMirrorSortingCorrection::ComputeR2(vector<Double_t>& ptR2Center, vector<Double_t>& ptR2Mirr,
611 vector<Double_t> ptM, vector<Double_t> ptC, vector<Double_t> ptR1,
612 TGeoNavigator* navi, TString option, TString mirrorTileName)
613{
614 //cout << endl << "//------------------------------ CbmRichCorrection: ComputeR2 ------------------------------//" << endl << endl;
615
616 vector<Double_t> normalMirr(3), ptCNew(3), ptTileCenter(3);
617 Double_t t1 = 0., t2 = 0., t3 = 0.;
618
619 if (option == "Corrected") {
620 // Use the correction information from text file, to the tile sphere center:
621 // Reading misalignment information from correction_param.txt text file.
622 Int_t lineCounter = 1, lineIndex = 0;
623 //TString str = fOutputDir + "correction_param_array_" + fStudyName + ".txt";
624 TString str = fCorrectionTableDir + "/correction_table/correction_param_array.txt";
625 string fileLine = "", strMisX = "", strMisY = "";
626 Double_t misX = 0., misY = 0.;
627 ifstream corrFile;
628 corrFile.open(str);
629
630 /*std::ifstream inFile(corrFile);
631 if(!inFile) {
632 cout << endl << "Failed to open file " << corrFile;
633 return;
634 }
635 double d1 = 0.;
636 double d2 = 0.;
637 while(!inFile.eof()) {
638 inFile >> d1 >> d2;
639 cout << d1 << " " << d2 << endl;;
640 }*/
641
642 if (corrFile.is_open()) {
643 while (!corrFile.eof()) {
644 getline(corrFile, fileLine);
645 lineIndex = fileLine.find(mirrorTileName, 0);
646 if (lineIndex != string::npos) {
647 //cout << mirrorTileName << " has been found in the file at line: " << lineCounter << " and position: " << lineIndex << "." << endl;
648 break;
649 }
650 lineCounter++;
651 }
652 // getline(corrFile, strMisX);
653 corrFile >> misY;
654 //cout << "number at line: " << lineCounter+2 << " = " << misX << "." << endl;
655 // getline(corrFile, strMisY);
656 corrFile >> misX;
657 //cout << "number at line: " << lineCounter+3 << " = " << misY << "." << endl;
658
659 /*std::istringstream i1(strMisX);
660 i1 >> misX;
661 std::istringstream i2(strMisY);
662 i2 >> misY;
663 double sum = misX + misY;
664 cout << "x1 = " << misX << ", x2 = " << misY << ", sum = " << sum << endl;*/
665
666 corrFile.close();
667 }
668 else {
669 cout << "Error in CbmRichCorrection: unable to open parameter file!" << endl;
670 cout << "Parameter file path: " << str << endl << endl;
671 //sleep(5);
672 }
673 //cout << "Misalignment parameters read from file = [" << outputFit.at(0) << " ; " << outputFit.at(1) << " ; " << outputFit.at(2) << " ; " << outputFit.at(3) << "]" << endl;
674
675 //ptCNew.at(0) = TMath::Abs(ptC.at(0) - TMath::Abs(outputFit.at(3)));
676 //ptCNew.at(1) = TMath::Abs(ptC.at(1) - TMath::Abs(outputFit.at(2)));
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];
684 //cout << "Mirror tile center coordinates = {" << ptTileCenter.at(0) << ", " << ptTileCenter.at(1) << ", " << ptTileCenter.at(2) << "}" << endl;
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);
691 //cout << "{x, y, z} = {" << x << ", " << y << ", " << z << "}, dist = " << dist << " and z2 = " << z2 << endl;
692 dist2 = TMath::Sqrt(x + y + TMath::Power(z2 - ptTileCenter.at(2), 2));
693 //cout << "dist2 = " << dist2 << endl;
694 ptCNew.at(2) += z2;
695 cout << "Sphere center coordinates of the rotated mirror tile, after "
696 "correction, = {"
697 << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
698 }
699 else if (option == "Uncorrected") {
700 // Keep the same tile sphere center, with no correction information.
701 ptCNew = ptC;
702 cout << "Sphere center coordinates of the rotated mirror tile, without "
703 "correction = {"
704 << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
705 }
706 else {
707 //cout << "No input given in function ComputeR2! Uncorrected parameters for the sphere center of the tile will be used!" << endl;
708 ptCNew = ptC;
709 //cout << "Sphere center coordinates of the rotated mirror tile, without correction = {" << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
710 }
711
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));
721 //cout << "Calculated and normalized normal of mirror tile = {" << normalMirr.at(0) << ", " << normalMirr.at(1) << ", " << normalMirr.at(2) << "}" << endl;
722
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);
730 t2 =
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);
738 /*//SAME AS calculation of t2 above
739 t3 = ((ptR1.at(0)-ptCNew.at(0))*(ptCNew.at(0)-ptM.at(0)) + (ptR1.at(1)-ptCNew.at(1))*(ptCNew.at(1)-ptM.at(1)) + (ptR1.at(2)-ptCNew.at(2))*(ptCNew.at(2)-ptM.at(2)))/TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0),2)+TMath::Power(ptCNew.at(1) - ptM.at(1),2)+TMath::Power(ptCNew.at(2) - ptM.at(2),2));
740 reflectedPtCooVectSphereUnity[0] = 2*(ptCNew.at(0)+t3*(normalMirr.at(0)))-ptR1.at(0);
741 reflectedPtCooVectSphereUnity[1] = 2*(ptCNew.at(1)+t3*(normalMirr.at(1)))-ptR1.at(1);
742 reflectedPtCooVectSphereUnity[2] = 2*(ptCNew.at(2)+t3*(normalMirr.at(2)))-ptR1.at(2);*/
743 //cout << "* using mirror point M to define \U00000394: {" << ptR2Center.at(0) << ", " << ptR2Center.at(1) << ", " << ptR2Center.at(2) << "}" << endl;
744 //cout << "Ref Pt Coo using unity Mirror-Sphere vector & sphere pt = {" << reflectedPtCooVectSphereUnity[0] << ", " << reflectedPtCooVectSphereUnity[1] << ", " << reflectedPtCooVectSphereUnity[2] << "}" << endl << endl;
745 //cout << "NofPMTPoints = " << NofPMTPoints << endl;
746
747 //cout << "Coordinates of point R2 on reflective plane after reflection on the mirror tile:" << endl;
748 //cout << "* using sphere center C to define \U00000394: {" << ptR2Mirr.at(0) << ", " << ptR2Mirr.at(1) << ", " << ptR2Mirr.at(2) << "}" << endl;
749}
750
751void CbmRichMirrorSortingCorrection::ComputeP(vector<Double_t>& ptPMirr, vector<Double_t>& ptPR2,
752 vector<Double_t> normalPMT, vector<Double_t> ptM,
753 vector<Double_t> ptR2Mirr, Double_t constantePMT)
754{
755 //cout << endl << "//------------------------------ CbmRichCorrection: ComputeP ------------------------------//" << endl << endl;
756
757 Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
758
759 k1 = -1
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));
766 k2 = -1
767 * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1) + normalPMT.at(2) * ptR2Mirr.at(2)
768 + constantePMT)
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));
774 //cout << "Coordinates of point P on PMT plane, after reflection on the mirror tile and extrapolation to the PMT plane:" << endl;
775 //cout << "* using mirror point M to define \U0001D49F ': {" << ptPMirr.at(0) << ", " << ptPMirr.at(1) << ", " << ptPMirr.at(2) << "}" << endl;
776 //cout << "* using reflected point R2 to define \U0001D49F ': {" << ptPR2.at(0) << ", " << ptPR2.at(1) << ", " << ptPR2.at(2) << "}" << endl;
777 checkCalc1 =
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.):"
781 << endl;
782 //cout << "* using mirror point M, checkCalc = " << checkCalc1 << endl;
783 checkCalc2 =
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;
786}
787
788void CbmRichMirrorSortingCorrection::FillHistProjection(TVector3 outPosIdeal, TVector3 outPosUnCorr, TVector3 outPos,
789 CbmRichRingLight ringL, vector<Double_t> normalPMT,
790 Double_t constantePMT, string str)
791{
792 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0, distToExtrapTrackHitInPlane = 0,
793 distToExtrapTrackHitInPlaneUnCorr = 0, distToExtrapTrackHitInPlaneIdeal = 0;
794 string histoNameX = "", histoNameY = "";
795 string nameX = "", nameY = "";
796
797 ringCenter[0] = ringL.GetCenterX();
798 ringCenter[1] = ringL.GetCenterY();
799 ringCenter[2] =
800 -1 * ((normalPMT.at(0) * ringCenter[0] + normalPMT.at(1) * ringCenter[1] + constantePMT) / normalPMT.at(2));
801
802 // Calculation using the corrected mirror hit/position, using the correction method
803 vector<Double_t> r(3),
804 p(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
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] << "}"
808 << endl;
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;
812
813 nameX = string("DiffCorrX_") + str;
814 nameY = string("DiffCorrY_") + str;
815 for (std::map<string, TH1D*>::iterator it = fDiffHistoMap.begin(); it != fDiffHistoMap.end(); ++it) {
816 if (nameX.compare(it->first) == 0) {
817 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - p.at(0)));
818 }
819 if (nameY.compare(it->first) == 0) {
820 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - p.at(1)));
821 }
822 }
823
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 "
829 "plane = "
830 << distToExtrapTrackHitInPlane << endl;
831 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
832 fHM->H1("fhDistanceCorrected")->Fill(distToExtrapTrackHitInPlane);
833
834 // Calculation using the uncorrected mirror hit/position
835 vector<Double_t> pUncorr(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
836 pUncorr.at(0) = outPosUnCorr.x(), pUncorr.at(1) = outPosUnCorr.y(), pUncorr.at(2) = outPosUnCorr.z();
837 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
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;
841
842 nameX = string("DiffUncorrX_") + str;
843 nameY = string("DiffUncorrY_") + str;
844 for (std::map<string, TH1D*>::iterator it = fDiffHistoMap.begin(); it != fDiffHistoMap.end(); ++it) {
845 if (nameX.compare(it->first) == 0) {
846 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - pUncorr.at(0)));
847 }
848 if (nameY.compare(it->first) == 0) {
849 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - pUncorr.at(1)));
850 }
851 }
852
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 "
857 "plane = "
858 << distToExtrapTrackHitInPlaneUnCorr << endl;
859
860 // Calculation using the ideally corrected mirror hit/position
861 vector<Double_t> pIdeal(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
862 pIdeal.at(0) = outPosIdeal.x(), pIdeal.at(1) = outPosIdeal.y(), pIdeal.at(2) = outPosIdeal.z();
863 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
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;
867
868 nameX = string("DiffIdealX_") + str;
869 nameY = string("DiffIdealY_") + str;
870 for (std::map<string, TH1D*>::iterator it = fDiffHistoMap.begin(); it != fDiffHistoMap.end(); ++it) {
871 if (nameX.compare(it->first) == 0) {
872 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(0) - pIdeal.at(0)));
873 }
874 if (nameY.compare(it->first) == 0) {
875 fDiffHistoMap[it->first]->Fill(TMath::Abs(r.at(1) - pIdeal.at(1)));
876 }
877 }
878
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 "
883 "plane = "
884 << distToExtrapTrackHitInPlaneIdeal << endl
885 << endl;
886 //}
887 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
888
889 if (distToExtrapTrackHitInPlane < 25. && distToExtrapTrackHitInPlaneUnCorr < 25.
890 && distToExtrapTrackHitInPlaneIdeal < 25.) {
891 //cout << "SEVERAL: " << distToExtrapTrackHitInPlane << endl << distToExtrapTrackHitInPlaneUnCorr << endl << distToExtrapTrackHitInPlaneIdeal << endl << endl;
892 fTrackCenterDistanceCorrected += distToExtrapTrackHitInPlane;
893 fTrackCenterDistanceUncorrected += distToExtrapTrackHitInPlaneUnCorr;
894 fTrackCenterDistanceIdeal += distToExtrapTrackHitInPlaneIdeal;
895 }
896 else {
897 cout << "Distance hit-ring too high!" << endl;
898 //sleep(5);
899 }
900}
901
902void CbmRichMirrorSortingCorrection::DrawMap(Int_t axisX, Int_t axisY)
903{
904 //vector<TCanvas> Can(40);
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;
913
914 for (std::map<string, TH1D*>::iterator it = fDiffHistoMap.begin(); it != fDiffHistoMap.end(); ++it) {
915 if (it->first.find(ssX.str().c_str()) != std::string::npos && it->second->GetEntries() > fThreshold) {
916 // cout << "ssX: " << ssX.str().c_str() << " and histo key: " << it->first << endl;
917 cout << "Histo entries: " << it->second->GetEntries() << endl;
918 stringVectX.push_back(it->first);
919 histoVectX.push_back(it->second);
920 }
921 else if (it->first.find(ssY.str().c_str()) != std::string::npos && it->second->GetEntries() > fThreshold) {
922 // cout << "ssY: " << ssY.str().c_str() << " and histo key: " << it->first << endl;
923 stringVectY.push_back(it->first);
924 histoVectY.push_back(it->second);
925 }
926 }
927
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);
931 c1->Divide(3, 1);
932 for (Int_t i = 0; i < 3; i++) {
933 c1->cd(counterX);
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();
939 counterX++;
940 }
941 Cbm::SaveCanvasAsImage(c1, string(fOutputDir.Data() + fStudyName), "png");
942 }
943
944 if (!histoVectY.empty()) {
945 TCanvas* c2 = new TCanvas(ssY.str().c_str(), ssY.str().c_str(), canX, canY);
946 c2->Divide(3, 1);
947 for (Int_t i = 0; i < 3; i++) {
948 c2->cd(counterY);
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();
954 counterY++;
955 }
956 Cbm::SaveCanvasAsImage(c2, string(fOutputDir.Data() + fStudyName), "png");
957 }
958
959 histoVectX.clear();
960 histoVectY.clear();
961}
962
964{
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;
967
968 // vector<TCanvas> Can;
969 // std::stringstream ss;
970 // ss << "X_mirror_" << X << "_" << Y;
971 // ss.str();
972 // ss.str().c_str();
973 // for ( Int_t i=0; i<10; i++ ) {
974 // ss << can << i;
975 // TCan* c = new TCan(ss.str().c_str(), )
976 // draw
977 // }
978
979 // vector<TCanvas> CanVect(80);
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++) {
983 DrawMap(i, j);
984 }
985 }
986
987 /* TCanvas* can1 = new TCanvas("X_mirror_tile_0_1","X_mirror_tile_0_1",1500,400);
988 can1->Divide(3,1);
989 TCanvas* can2 = new TCanvas("Y_mirror_tile_0_1","Y_mirror_tile_0_1",1500,400);
990 can2->Divide(3,1);
991 TCanvas* can3 = new TCanvas("X_mirror_tile_1_5","X_mirror_tile_1_5",1500,400);
992 can3->Divide(3,1);
993 TCanvas* can4 = new TCanvas("Y_mirror_tile_1_5","Y_mirror_tile_1_5",1500,400);
994 can4->Divide(3,1);
995 TCanvas* can5 = new TCanvas("X_mirror_tile_2_8","X_mirror_tile_2_8",1500,400);
996 can5->Divide(3,1);
997 TCanvas* can6 = new TCanvas("Y_mirror_tile_2_8","Y_mirror_tile_2_8",1500,400);
998 can6->Divide(3,1);
999 TCanvas* can7 = new TCanvas("X_mirror_tile_1_3","X_mirror_tile_1_3",1500,400);
1000 can7->Divide(3,1);
1001 TCanvas* can8 = new TCanvas("Y_mirror_tile_1_3","Y_mirror_tile_1_3",1500,400);
1002 can8->Divide(3,1);
1003 TCanvas* can9 = new TCanvas("X_mirror_tile_1_4","X_mirror_tile_1_4",1500,400);
1004 can9->Divide(3,1);
1005 TCanvas* can10 = new TCanvas("Y_mirror_tile_1_4","Y_mirror_tile_1_4",1500,400);
1006 can10->Divide(3,1);
1007 TCanvas* can11 = new TCanvas("X_mirror_tile_0_8","X_mirror_tile_0_8",1500,400);
1008 can11->Divide(3,1);
1009 TCanvas* can12 = new TCanvas("Y_mirror_tile_0_8","Y_mirror_tile_0_8",1500,400);
1010 can12->Divide(3,1);
1011
1012
1013 for (std::map<string,TH1D*>::iterator it=fDiffHistoMap.begin(); it!=fDiffHistoMap.end(); ++it) {
1014 if ( it->first.find("X_mirror_tile_0_1")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1015 can1->cd(counter1);
1016 fDiffHistoMap[it->first]->Draw();
1017 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1018 fDiffHistoMap[it->first]->SetLineColor(counter1+1);
1019 fDiffHistoMap[it->first]->SetLineWidth(2);
1020 fDiffHistoMap[it->first]->Write();
1021 counter1++;
1022 }
1023 else if ( it->first.find("Y_mirror_tile_0_1")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1024 can2->cd(counter2);
1025 fDiffHistoMap[it->first]->Draw();
1026 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1027 fDiffHistoMap[it->first]->SetLineColor(counter2+1);
1028 fDiffHistoMap[it->first]->SetLineWidth(2);
1029 fDiffHistoMap[it->first]->Write();
1030 counter2++;
1031 }
1032 else if ( it->first.find("X_mirror_tile_1_5")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1033 can3->cd(counter3);
1034 fDiffHistoMap[it->first]->Draw();
1035 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1036 fDiffHistoMap[it->first]->SetLineColor(counter3+1);
1037 fDiffHistoMap[it->first]->SetLineWidth(2);
1038 fDiffHistoMap[it->first]->Write();
1039 counter3++;
1040 }
1041 else if ( it->first.find("Y_mirror_tile_1_5")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1042 can4->cd(counter4);
1043 fDiffHistoMap[it->first]->Draw();
1044 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1045 fDiffHistoMap[it->first]->SetLineColor(counter4+1);
1046 fDiffHistoMap[it->first]->SetLineWidth(2);
1047 fDiffHistoMap[it->first]->Write();
1048 counter4++;
1049 }
1050 else if ( it->first.find("X_mirror_tile_2_8")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1051 can5->cd(counter5);
1052 fDiffHistoMap[it->first]->Draw();
1053 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1054 fDiffHistoMap[it->first]->SetLineColor(counter5+1);
1055 fDiffHistoMap[it->first]->SetLineWidth(2);
1056 fDiffHistoMap[it->first]->Write();
1057 counter5++;
1058 }
1059 else if ( it->first.find("Y_mirror_tile_2_8")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1060 can6->cd(counter6);
1061 fDiffHistoMap[it->first]->Draw();
1062 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1063 fDiffHistoMap[it->first]->SetLineColor(counter6+1);
1064 fDiffHistoMap[it->first]->SetLineWidth(2);
1065 fDiffHistoMap[it->first]->Write();
1066 counter6++;
1067 }
1068 else if ( it->first.find("X_mirror_tile_1_3")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1069 can7->cd(counter7);
1070 fDiffHistoMap[it->first]->Draw();
1071 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1072 fDiffHistoMap[it->first]->SetLineColor(counter7+1);
1073 fDiffHistoMap[it->first]->SetLineWidth(2);
1074 counter7++;
1075 }
1076 else if ( it->first.find("Y_mirror_tile_1_3")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1077 can8->cd(counter8);
1078 fDiffHistoMap[it->first]->Draw();
1079 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1080 fDiffHistoMap[it->first]->SetLineColor(counter8+1);
1081 fDiffHistoMap[it->first]->SetLineWidth(2);
1082 counter8++;
1083 }
1084 else if ( it->first.find("X_mirror_tile_1_4")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1085 can9->cd(counter9);
1086 fDiffHistoMap[it->first]->Draw();
1087 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1088 fDiffHistoMap[it->first]->SetLineColor(counter9+1);
1089 fDiffHistoMap[it->first]->SetLineWidth(2);
1090 fDiffHistoMap[it->first]->Write();
1091 counter9++;
1092 }
1093 else if ( it->first.find("Y_mirror_tile_1_4")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1094 can10->cd(counter10);
1095 //fDiffHistoMap[it->first]->Draw();
1096 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1097 fDiffHistoMap[it->first]->SetLineColor(counter10+1);
1098 fDiffHistoMap[it->first]->SetLineWidth(2);
1099 fDiffHistoMap[it->first]->Write();
1100 counter10++;
1101 }
1102 else if ( it->first.find("X_mirror_tile_0_8")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1103 can11->cd(counter11);
1104 fDiffHistoMap[it->first]->Draw();
1105 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1106 fDiffHistoMap[it->first]->SetLineColor(counter11+1);
1107 fDiffHistoMap[it->first]->SetLineWidth(2);
1108 fDiffHistoMap[it->first]->Write();
1109 counter11++;
1110 }
1111 else if ( it->first.find("Y_mirror_tile_0_8")!=std::string::npos && it->second->GetEntries() > fThreshold ) {
1112 can12->cd(counter12);
1113 fDiffHistoMap[it->first]->Draw();
1114 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1115 fDiffHistoMap[it->first]->SetLineColor(counter12+1);
1116 fDiffHistoMap[it->first]->SetLineWidth(2);
1117 fDiffHistoMap[it->first]->Write();
1118 counter12++;
1119 }
1120 else if ( it->second->GetEntries() > fThreshold ) {
1121 TCanvas* can = new TCanvas();
1122 fDiffHistoMap[it->first]->Draw();
1123 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1124 fDiffHistoMap[it->first]->SetLineColor(4);
1125 fDiffHistoMap[it->first]->SetLineWidth(2);
1126 fDiffHistoMap[it->first]->Write();
1127 }
1128 }
1129
1130 Cbm::SaveCanvasAsImage(can1, string(fOutputDir.Data()+fStudyName), "png");
1131 Cbm::SaveCanvasAsImage(can2, string(fOutputDir.Data()+fStudyName), "png");
1132 Cbm::SaveCanvasAsImage(can3, string(fOutputDir.Data()+fStudyName), "png");
1133 Cbm::SaveCanvasAsImage(can4, string(fOutputDir.Data()+fStudyName), "png");*/
1134 //Cbm::SaveCanvasAsImage(can5, string(fOutputDir.Data()+fStudyName), "png");
1135 //Cbm::SaveCanvasAsImage(can6, string(fOutputDir.Data()+fStudyName), "png");
1136 //Cbm::SaveCanvasAsImage(can7, string(fOutputDir.Data()+fStudyName), "png");
1137 //Cbm::SaveCanvasAsImage(can8, string(fOutputDir.Data()+fStudyName), "png");
1138 //Cbm::SaveCanvasAsImage(can9, string(fOutputDir.Data()+fStudyName), "png");
1139 //Cbm::SaveCanvasAsImage(can10, string(fOutputDir.Data()+fStudyName), "png");
1140 //Cbm::SaveCanvasAsImage(can11, string(fOutputDir.Data()+fStudyName), "png");
1141 //Cbm::SaveCanvasAsImage(can12, string(fOutputDir.Data()+fStudyName), "png");
1142
1143 /* char title[128];
1144 for (std::map<string,TH1D*>::iterator it=fDiffHistoMap.begin(); it!=fDiffHistoMap.end(); ++it) {
1145 TCanvas* can = new TCanvas();
1146 fDiffHistoMap[it->first]->Draw();
1147 //title = (it->first).c_str();
1148 fDiffHistoMap[it->first]->SetTitle((it->first).c_str());
1149 fDiffHistoMap[it->first]->SetLineColor(4);
1150 fDiffHistoMap[it->first]->SetLineWidth(2);
1151 //Cbm::SaveCanvasAsImage(can, string(fOutputDir.Data()), "png");
1152 }
1153*/
1154}
1155
1157{
1158 for (Int_t iTrack = 0; iTrack < fGlobalTracks->GetEntriesFast(); iTrack++) {
1159 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
1160 Int_t stsId = globalTrack->GetStsTrackIndex();
1161 Int_t richId = globalTrack->GetRichRingIndex();
1162 if (stsId < 0 || richId < 0) continue;
1163 const CbmTrackMatchNew* stsTrackMatch = static_cast<const CbmTrackMatchNew*>(fStsTrackMatches->At(stsId));
1164 if (stsTrackMatch == NULL) continue;
1165 int stsMcTrackId = stsTrackMatch->GetMatchedLink().GetIndex();
1166 const CbmTrackMatchNew* richRingMatch = static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
1167 if (richRingMatch == NULL) continue;
1168 int richMcTrackId = richRingMatch->GetMatchedLink().GetIndex();
1169 const CbmRichRing* ring = static_cast<const CbmRichRing*>(fRichRings->At(richId));
1170 if (NULL == ring) continue;
1171 Double_t rtDistance = CbmRichUtil::GetRingTrackDistance(iTrack);
1172 Double_t rtDistanceX = CbmRichUtil::GetRingTrackDistanceX(iTrack);
1173 Double_t rtDistanceY = CbmRichUtil::GetRingTrackDistanceY(iTrack);
1174 Double_t xc = ring->GetCenterX();
1175 Double_t yc = ring->GetCenterY();
1176
1177 CbmMCTrack* mctrack = static_cast<CbmMCTrack*>(fMCTracks->At(stsMcTrackId));
1178 if (mctrack == NULL) continue;
1179 bool isEl = IsMcPrimaryElectron(mctrack);
1180
1181 if (isEl && stsMcTrackId == richMcTrackId) {
1182 Int_t k = 1;
1183 stringstream ss;
1184 ss << k;
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);
1192 }
1193 ss.str("");
1194 }
1195 }
1196}
1197
1199 const FairTrackParam* pTrack, const CbmMCTrack* mcTrack)
1200{
1201 Double_t xRing = richRing->GetCenterX();
1202 Double_t yRing = richRing->GetCenterY();
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);
1206
1207 bool isEl = IsMcPrimaryElectron(mcTrack);
1208 if (isEl) {
1209 Int_t k = 2;
1210 stringstream ss;
1211 ss << k;
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);
1219 }
1220 ss.str("");
1221 }
1222}
1223
1225{
1226 if (mctrack == NULL) return false;
1227 int pdg = TMath::Abs(mctrack->GetPdgCode());
1228 if (mctrack->GetGeantProcessId() == kPPrimary && pdg == 11) return true;
1229 return false;
1230}
1231
1233{
1234 stringstream ss;
1235 ss << k;
1236 {
1237 TCanvas* c = fHM->CreateCanvas("fh_ring_track_distance_vs_xy_truematch" + k,
1238 "fh_ring_track_distance_vs_xy_truematch" + k, 1800, 600);
1239 c->Divide(3, 1);
1240 c->cd(1);
1241 DrawH3Profile(fHM->H3("fhRingTrackDistVsXYTruematch" + ss.str()), true, false, 0., 2.);
1242 c->cd(2);
1243 DrawH2WithProfile(fHM->H2("fhRingTrackDistVsXTruematch" + ss.str()), false, true);
1244 fHM->H2("fhRingTrackDistVsXTruematch" + ss.str())->GetYaxis()->SetRangeUser(0., 1.5);
1245 c->cd(3);
1246 DrawH2WithProfile(fHM->H2("fhRingTrackDistVsYTruematch" + ss.str()), false, true);
1247 fHM->H2("fhRingTrackDistVsYTruematch" + ss.str())->GetYaxis()->SetRangeUser(0., 1.5);
1248 }
1249
1250 Double_t range = 2.5;
1251 {
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);
1254 c->Divide(2, 1);
1255 c->cd(1);
1256 DrawH3Profile(fHM->H3("fhRingTrackDistDiffXRingVsXYTruematch" + ss.str()), true, false, -range, range);
1257 c->cd(2);
1258 DrawH3Profile(fHM->H3("fhRingTrackDistDiffYRingVsXYTruematch" + ss.str()), true, false, -range, range);
1259 }
1260}
1261
1263{
1264 TCanvas* c1 = fHM->CreateCanvas("fh_distance_distribution_corr", "fh_distance_distribution_corr", 600, 900);
1265 DrawH1andFitGauss(fHM->H1("fDistCorr"), true, false, -2., 4.);
1266}
1267
1269{
1270 //DrawHistProjection();
1274
1275 TDirectory* oldir = gDirectory;
1276 TFile* outFile = FairRootManager::Instance()->GetOutFile();
1277 if (outFile != NULL) {
1278 fHM->WriteToFile();
1279 }
1280 gDirectory->cd(oldir->GetPath());
1281
1282 string str1 = fOutputDir.Data();
1283 string str2 = fStudyName.Data();
1284 string str = str1 + "/" + str2;
1285 cout << endl << endl << "output string: " << str << endl << endl;
1286 fHM->SaveCanvasToImage(str);
1287
1288 // TString s = fOutputDir + "track_ring_distances_" + fStudyName + ".txt";
1289 TString s = fOutputDir + "/correction_table/track_ring_distances.txt";
1290 ofstream corrFile;
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 = "
1296 << endl;
1297 corrFile << "Mean distance between track hit and ring center ; uncorrected case = "
1299 << endl;
1300 corrFile << "Mean distance between track hit and ring center ; ideal case = "
1301 << fTrackCenterDistanceIdeal / fEventNb << " and total sum = " << fTrackCenterDistanceIdeal << endl;
1302 }
1303 else {
1304 cout << "Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
1305 "parameter file!"
1306 << endl;
1307 }
1308
1309 cout << "number of events: " << fEventNb << endl;
1310 cout << setprecision(9) << "Mean distance between track hit and ring center ; corrected case = "
1311 << fTrackCenterDistanceCorrected / fEventNb << " and total sum = " << fTrackCenterDistanceCorrected << endl;
1312 cout << "Mean distance between track hit and ring center ; uncorrected case = "
1313 << fTrackCenterDistanceUncorrected / fEventNb << " and total sum = " << fTrackCenterDistanceUncorrected << endl;
1314 cout << "Mean distance between track hit and ring center ; ideal case = " << fTrackCenterDistanceIdeal / fEventNb
1315 << " and total sum = " << fTrackCenterDistanceIdeal << endl;
1316}
1318
1319
1320 /*void CbmRichMirrorSortingCorrection::FillHistProjection(TVector3 outPosIdeal, TVector3 outPosUnCorr, TVector3 outPos, CbmRichRingLight ringL, vector<Double_t> normalPMT, Double_t constantePMT, string str)
1321{
1322 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0, distToExtrapTrackHitInPlane = 0, distToExtrapTrackHitInPlaneUnCorr = 0, distToExtrapTrackHitInPlaneIdeal = 0;
1323 string histoNameX = "", histoNameY = "";
1324 string nameX = "", nameY = "";
1325
1326 ringCenter[0] = ringL.GetCenterX();
1327 ringCenter[1] = ringL.GetCenterY();
1328 ringCenter[2] = -1*((normalPMT.at(0)*ringCenter[0] + normalPMT.at(1)*ringCenter[1] + constantePMT)/normalPMT.at(2));
1329
1330 vector<Double_t> r(3), p(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1331 r.at(0) = ringCenter[0], r.at(1) = ringCenter[1], r.at(2) = ringCenter[2];
1332 p.at(0) = outPos.x(), p.at(1) = outPos.y(), p.at(2) = outPos.z();
1333 cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1334 cout << "Difference in X = " << TMath::Abs(r.at(0) - p.at(0)) << "; \t" << "Difference in Y = " << TMath::Abs(r.at(1) - p.at(1)) << "; \t" << "Difference in Z = " << TMath::Abs(r.at(2) - p.at(2)) << endl;
1335
1336 nameX = string("X_") + str;
1337 nameY = string("Y_") + str;
1338 for (std::map<string,string>::iterator it=fDiffHistoMap.begin(); it!=fDiffHistoMap.end(); ++it) {
1339 if ( nameX.compare(it->second) == 0 ) {
1340 fHM2->H1(it->first)->Fill(TMath::Abs(r.at(0) - p.at(0)));
1341 }
1342 if ( nameY.compare(it->second) == 0 ) {
1343 fHM2->H1(it->first)->Fill(TMath::Abs(r.at(1) - p.at(1)));
1344 }
1345 }
1346
1347 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0),2) + TMath::Power(r.at(1) - p.at(1),2) + TMath::Power(r.at(2) - p.at(2),2));
1348 distToExtrapTrackHitInPlane = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0),2) + TMath::Power(r.at(1) - p.at(1),2));
1349 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
1350 fHM->H1("fhDistanceCorrected")->Fill(distToExtrapTrackHitInPlane);
1351 //cout << "Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl;
1352 cout << "Distance between fitted ring center and extrapolated track hit in plane = " << distToExtrapTrackHitInPlane << endl;
1353
1354 vector<Double_t> pUnCorr(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1355 pUnCorr.at(0) = outPosUnCorr.x(), pUnCorr.at(1) = outPosUnCorr.y(), pUnCorr.at(2) = outPosUnCorr.z();
1356 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1357 cout << "Difference in X w/o correction = " << TMath::Abs(r.at(0) - pUnCorr.at(0)) << "; \t" << "Difference in Y = " << TMath::Abs(r.at(1) - pUnCorr.at(1)) << "; \t" << "Difference in Z = " << TMath::Abs(r.at(2) - pUnCorr.at(2)) << endl;
1358 fHM->H1("fhDifferenceXUncorrected")->Fill(TMath::Abs(r.at(0) - pUnCorr.at(0)));
1359 fHM->H1("fhDifferenceYUncorrected")->Fill(TMath::Abs(r.at(1) - pUnCorr.at(1)));
1360 distToExtrapTrackHitInPlaneUnCorr = TMath::Sqrt(TMath::Power(r.at(0) - pUnCorr.at(0),2) + TMath::Power(r.at(1) - pUnCorr.at(1),2));
1361 fHM->H1("fhDistanceUncorrected")->Fill(distToExtrapTrackHitInPlaneUnCorr);
1362 cout << "Distance between fitted ring center and extrapolated track hit in plane = " << distToExtrapTrackHitInPlaneUnCorr << endl;
1363
1364 vector<Double_t> pIdeal(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1365 pIdeal.at(0) = outPosIdeal.x(), pIdeal.at(1) = outPosIdeal.y(), pIdeal.at(2) = outPosIdeal.z();
1366 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1367 cout << "Difference in X w/ ideal correction = " << TMath::Abs(r.at(0) - pIdeal.at(0)) << "; \t" << "Difference in Y = " << TMath::Abs(r.at(1) - pIdeal.at(1)) << "; \t" << "Difference in Z = " << TMath::Abs(r.at(2) - pIdeal.at(2)) << endl;
1368 fHM->H1("fhDifferenceXIdeal")->Fill(TMath::Abs(r.at(0) - pIdeal.at(0)));
1369 fHM->H1("fhDifferenceYIdeal")->Fill(TMath::Abs(r.at(1) - pIdeal.at(1)));
1370 distToExtrapTrackHitInPlaneIdeal = TMath::Sqrt(TMath::Power(r.at(0) - pIdeal.at(0),2) + TMath::Power(r.at(1) - pIdeal.at(1),2));
1371 fHM->H1("fhDistanceIdeal")->Fill(distToExtrapTrackHitInPlaneIdeal);
1372 cout << "Distance between fitted ring center and extrapolated track hit in plane = " << distToExtrapTrackHitInPlaneIdeal << endl << endl;
1373 //}
1374 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
1375}
1376
1377/*void CbmRichMirrorSortingCorrection::FillHistProjection(TVector3 outPosIdeal, TVector3 outPosUnCorr, TVector3 outPos, CbmRichRingLight ringL, vector<Double_t> normalPMT, Double_t constantePMT)
1378{
1379 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0, distToExtrapTrackHitInPlane = 0, distToExtrapTrackHitInPlaneUnCorr = 0, distToExtrapTrackHitInPlaneIdeal = 0;
1380
1381 ringCenter[0] = ringL.GetCenterX();
1382 ringCenter[1] = ringL.GetCenterY();
1383 ringCenter[2] = -1*((normalPMT.at(0)*ringCenter[0] + normalPMT.at(1)*ringCenter[1] + constantePMT)/normalPMT.at(2));
1384
1385 vector<Double_t> r(3), p(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1386 r.at(0) = TMath::Abs(ringCenter[0]), r.at(1) = TMath::Abs(ringCenter[1]), r.at(2) = TMath::Abs(ringCenter[2]);
1387 p.at(0) = TMath::Abs(outPos.x()), p.at(1) = TMath::Abs(outPos.y()), p.at(2) = TMath::Abs(outPos.z());
1388 cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1389 cout << "Difference in X = " << TMath::Abs(r.at(0) - p.at(0)) << "; \t" << "Difference in Y = " << TMath::Abs(r.at(1) - p.at(1)) << "; \t" << "Difference in Z = " << TMath::Abs(r.at(2) - p.at(2)) << endl;
1390 fHM->H1("fhDifferenceX")->Fill(TMath::Abs(r.at(0) - p.at(0)));
1391 fHM->H1("fhDifferenceY")->Fill(TMath::Abs(r.at(1) - p.at(1)));
1392 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0),2) + TMath::Power(r.at(1) - p.at(1),2) + TMath::Power(r.at(2) - p.at(2),2));
1393 distToExtrapTrackHitInPlane = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0),2) + TMath::Power(r.at(1) - p.at(1),2));
1394 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
1395 fHM->H1("fhDistanceCorrected")->Fill(distToExtrapTrackHitInPlane);
1396 //cout << "Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl;
1397 cout << "Distance between fitted ring center and extrapolated track hit in plane = " << distToExtrapTrackHitInPlane << endl;
1398
1399 vector<Double_t> pUnCorr(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1400 pUnCorr.at(0) = TMath::Abs(outPosUnCorr.x()), pUnCorr.at(1) = TMath::Abs(outPosUnCorr.y()), pUnCorr.at(2) = TMath::Abs(outPosUnCorr.z());
1401 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1402 cout << "Difference in X w/o correction = " << TMath::Abs(r.at(0) - pUnCorr.at(0)) << "; \t" << "Difference in Y = " << TMath::Abs(r.at(1) - pUnCorr.at(1)) << "; \t" << "Difference in Z = " << TMath::Abs(r.at(2) - pUnCorr.at(2)) << endl;
1403 fHM->H1("fhDifferenceXUncorrected")->Fill(TMath::Abs(r.at(0) - pUnCorr.at(0)));
1404 fHM->H1("fhDifferenceYUncorrected")->Fill(TMath::Abs(r.at(1) - pUnCorr.at(1)));
1405 distToExtrapTrackHitInPlaneUnCorr = TMath::Sqrt(TMath::Power(r.at(0) - pUnCorr.at(0),2) + TMath::Power(r.at(1) - pUnCorr.at(1),2));
1406 fHM->H1("fhDistanceUncorrected")->Fill(distToExtrapTrackHitInPlaneUnCorr);
1407 cout << "Distance between fitted ring center and extrapolated track hit in plane = " << distToExtrapTrackHitInPlaneUnCorr << endl;
1408
1409 vector<Double_t> pIdeal(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1410 pIdeal.at(0) = TMath::Abs(outPosIdeal.x()), pIdeal.at(1) = TMath::Abs(outPosIdeal.y()), pIdeal.at(2) = TMath::Abs(outPosIdeal.z());
1411 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1412 cout << "Difference in X w/ ideal correction = " << TMath::Abs(r.at(0) - pIdeal.at(0)) << "; \t" << "Difference in Y = " << TMath::Abs(r.at(1) - pIdeal.at(1)) << "; \t" << "Difference in Z = " << TMath::Abs(r.at(2) - pIdeal.at(2)) << endl;
1413 fHM->H1("fhDifferenceXIdeal")->Fill(TMath::Abs(r.at(0) - pIdeal.at(0)));
1414 fHM->H1("fhDifferenceYIdeal")->Fill(TMath::Abs(r.at(1) - pIdeal.at(1)));
1415 distToExtrapTrackHitInPlaneIdeal = TMath::Sqrt(TMath::Power(r.at(0) - pIdeal.at(0),2) + TMath::Power(r.at(1) - pIdeal.at(1),2));
1416 fHM->H1("fhDistanceIdeal")->Fill(distToExtrapTrackHitInPlaneIdeal);
1417 cout << "Distance between fitted ring center and extrapolated track hit in plane = " << distToExtrapTrackHitInPlaneIdeal << endl << endl;
1418 //}
1419 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
1420}
1421
1422
1423void CbmRichMirrorSortingCorrection::DrawHistProjection()
1424{
1425 char leg[128];
1426 int colorInd = 1;
1427 string diffCorrX = "DiffX_mirror_tile_2_8", diffCorrY = "DiffY_mirror_tile_2_8";
1428
1429 TCanvas* can1 = new TCanvas("Distance_Histos_Difference_X", "Distance_Histos_Difference_X", 1500, 400);
1430 can1->SetGrid(1,1);
1431 can1->Divide(3,1);
1432 can1->cd(1);
1433 fHM->H1("fhDifferenceXIdeal")->Draw();
1434 fHM->H1("fhDifferenceXIdeal")->SetTitle("Difference in X ideal");
1435 fHM->H1("fhDifferenceXIdeal")->SetLineColor(kBlue);
1436 fHM->H1("fhDifferenceXIdeal")->SetLineWidth(2);
1437 can1->cd(2);
1438 fDiffHistoMap[diffCorrX]->Draw();
1439 fDiffHistoMap[diffCorrX]->SetTitle("Difference in X corrected");
1440 fDiffHistoMap[diffCorrX]->SetLineColor(kRed);
1441 fDiffHistoMap[diffCorrX]->SetLineWidth(2);
1442 can1->cd(3);
1443 fHM->H1("fhDifferenceXUncorrected")->Draw();
1444 fHM->H1("fhDifferenceXUncorrected")->SetTitle("Difference in X uncorrected");
1445 fHM->H1("fhDifferenceXUncorrected")->SetLineColor(kGreen);
1446 fHM->H1("fhDifferenceXUncorrected")->SetLineWidth(2);
1447
1448 TCanvas* can2 = new TCanvas("Distance_Histos_Difference_Y", "Distance_Histos_Difference_Y", 1500, 400);
1449 can2->SetGrid(1,1);
1450 can2->Divide(3,1);
1451 can2->cd(1);
1452 fHM->H1("fhDifferenceYIdeal")->Draw();
1453 fHM->H1("fhDifferenceYIdeal")->SetTitle("Difference in Y ideal");
1454 fHM->H1("fhDifferenceYIdeal")->SetLineColor(kBlue);
1455 fHM->H1("fhDifferenceYIdeal")->SetLineWidth(2);
1456 can2->cd(2);
1457 fDiffHistoMap[diffCorrY]->Draw();
1458 fDiffHistoMap[diffCorrY]->SetTitle("Difference in Y corrected");
1459 fDiffHistoMap[diffCorrY]->SetLineColor(kRed);
1460 fDiffHistoMap[diffCorrY]->SetLineWidth(2);
1461 can2->cd(3);
1462 fHM->H1("fhDifferenceYUncorrected")->Draw();
1463 fHM->H1("fhDifferenceYUncorrected")->SetTitle("Difference in Y uncorrected");
1464 fHM->H1("fhDifferenceYUncorrected")->SetLineColor(kGreen);
1465 fHM->H1("fhDifferenceYUncorrected")->SetLineWidth(2);
1466
1467 /*can3->SetGrid(1,1);
1468 can3->Divide(2,1);
1469 can3->cd(1)->SetGrid(1,1);
1470 can3->cd(2)->SetGrid(1,1);
1471 can3->cd(1);
1472 TH1D* Clone1 = (TH1D*)fHM->H1("fhDifferenceXIdeal")->Clone();
1473 Clone1->GetXaxis()->SetTitleSize(0.04);
1474 Clone1->GetYaxis()->SetTitleSize(0.04);
1475 Clone1->GetXaxis()->SetLabelSize(0.03);
1476 Clone1->GetYaxis()->SetLabelSize(0.03);
1477 Clone1->GetXaxis()->CenterTitle();
1478 Clone1->GetYaxis()->CenterTitle();
1479 Clone1->SetTitle("Difference in X");
1480 Clone1->SetLineColor(kBlue);
1481 Clone1->SetLineWidth(2);
1482 Clone1->Rebin(2);
1483 Clone1->Draw();
1484 TH1D* Clone2 = (TH1D*)fHM->H1("fhDifferenceX")->Clone();
1485 Clone2->SetTitleSize(0.04);
1486 Clone2->SetLabelSize(0.03);
1487 Clone2->SetLineColor(kGreen);
1488 Clone2->SetLineWidth(2);
1489 Clone2->Rebin(2);
1490 Clone2->Draw("same");
1491 TH1D* Clone3 = (TH1D*)fHM->H1("fhDifferenceXUncorrected")->Clone();
1492 Clone3->SetTitleSize(0.04);
1493 Clone3->SetLabelSize(0.03);
1494 Clone3->SetLineColor(kRed);
1495 Clone3->SetLineWidth(2);
1496 Clone3->Rebin(2);
1497 Clone3->Draw("same");
1498 gStyle->Clear();
1499
1500 TLegend* LEG = new TLegend(0.3,0.78,0.5,0.88); // Set legend position
1501 LEG->SetBorderSize(1);
1502 LEG->SetFillColor(0);
1503 LEG->SetMargin(0.2);
1504 LEG->SetTextSize(0.02);
1505 sprintf(leg, "X diff uncorr");
1506 LEG->AddEntry(Clone3, leg, "l");
1507 sprintf(leg, "X diff corr");
1508 LEG->AddEntry(Clone2, leg, "l");
1509 sprintf(leg, "X diff ideal");
1510 LEG->AddEntry(Clone1, leg, "l");
1511 LEG->Draw();
1512
1513 can3->cd(2);
1514 can3->SetGrid(1,1);
1515 TH1D* Clone4 = (TH1D*)fHM->H1("fhDifferenceYIdeal")->Clone();
1516 Clone4->GetXaxis()->SetTitleSize(0.04);
1517 Clone4->GetYaxis()->SetTitleSize(0.04);
1518 Clone4->GetXaxis()->SetLabelSize(0.03);
1519 Clone4->GetYaxis()->SetLabelSize(0.03);
1520 Clone4->GetXaxis()->CenterTitle();
1521 Clone4->GetYaxis()->CenterTitle();
1522 Clone4->SetTitle("Difference in Y");
1523 Clone4->SetLineColor(kBlue);
1524 Clone4->SetLineWidth(2);
1525 Clone4->Rebin(2);
1526 Clone4->Draw();
1527 TH1D* Clone5 = (TH1D*)fHM->H1("fhDifferenceY")->Clone();
1528 Clone5->SetTitleSize(0.04);
1529 Clone5->SetLabelSize(0.03);
1530 Clone5->SetLineColor(kGreen);
1531 Clone5->SetLineWidth(2);
1532 Clone5->Rebin(2);
1533 Clone5->Draw("same");
1534 TH1D* Clone6 = (TH1D*)fHM->H1("fhDifferenceYUncorrected")->Clone();
1535 Clone6->SetTitleSize(0.04);
1536 Clone6->SetLabelSize(0.03);
1537 Clone6->SetLineColor(kRed);
1538 Clone6->SetLineWidth(2);
1539 Clone6->Rebin(2);
1540 Clone6->Draw("same");
1541
1542 TLegend* LEG1 = new TLegend(0.3,0.78,0.5,0.88); // Set legend position
1543 LEG1->SetBorderSize(1);
1544 LEG1->SetFillColor(0);
1545 LEG1->SetMargin(0.2);
1546 LEG1->SetTextSize(0.02);
1547 sprintf(leg, "Y diff uncorr");
1548 LEG1->AddEntry(Clone6, leg, "l");
1549 sprintf(leg, "Y diff corr");
1550 LEG1->AddEntry(Clone5, leg, "l");
1551 sprintf(leg, "Y diff ideal");
1552 LEG1->AddEntry(Clone4, leg, "l");
1553 LEG1->Draw();
1554
1555 /*can3->cd(3);
1556 //DrawH1(list_of(fHM->H1("fhDistanceCorrected"))(fHM->H1("fhDistanceUncorrected"))(fHM->H1("fhDistanceIdeal")), list_of("a corrected")("a uncorrected")("a ideal"), kLinear, kLog, true, 0.7, 0.7, 0.99, 0.99);
1557 TH1D* CloneDist1 = (TH1D*)fHM->H1("fhDistanceIdeal")->Clone();
1558 CloneDist1->GetXaxis()->SetTitleSize(0.04);
1559 CloneDist1->GetYaxis()->SetTitleSize(0.04);
1560 CloneDist1->GetXaxis()->SetLabelSize(0.03);
1561 CloneDist1->GetYaxis()->SetLabelSize(0.03);
1562 CloneDist1->GetXaxis()->CenterTitle();
1563 CloneDist1->GetYaxis()->CenterTitle();
1564 CloneDist1->SetTitle("Distance between extrapolated track center and fitted ring center");
1565 CloneDist1->SetLineColor(kBlue);
1566 CloneDist1->SetLineWidth(2);
1567 CloneDist1->Rebin(2);
1568 CloneDist1->Draw();
1569 TH1D* CloneDist2 = (TH1D*)fHM->H1("fhDistanceUncorrected")->Clone();
1570 CloneDist2->SetTitleSize(0.04);
1571 CloneDist2->SetTitleSize(0.04);
1572 CloneDist2->SetLineColor(kRed);
1573 CloneDist2->SetLineWidth(2);
1574 CloneDist2->Rebin(2);
1575 CloneDist2->Draw("same");
1576 TH1D* CloneDist3 = (TH1D*)fHM->H1("fhDistanceCorrected")->Clone();
1577 CloneDist3->SetTitleSize(0.04);
1578 CloneDist3->SetTitleSize(0.04);
1579 CloneDist3->SetLineColor(kGreen);
1580 CloneDist3->SetLineWidth(2);
1581 CloneDist3->Rebin(2);
1582 CloneDist3->Draw("same");
1583
1584 TLegend* LEG2 = new TLegend(0.3,0.78,0.5,0.88); // Set legend position
1585 LEG2->SetBorderSize(1);
1586 LEG2->SetFillColor(0);
1587 LEG2->SetMargin(0.2);
1588 LEG2->SetTextSize(0.02);
1589 sprintf(leg, "Distance uncorr");
1590 LEG2->AddEntry(CloneDist2, leg, "l");
1591 sprintf(leg, "Distance corr");
1592 LEG2->AddEntry(CloneDist3, leg, "l");
1593 sprintf(leg, "Distance ideal");
1594 LEG2->AddEntry(CloneDist1, leg, "l");
1595 LEG2->Draw();*/
1596 /* gStyle->SetOptStat(000000);
1597
1598 Cbm::SaveCanvasAsImage(can1, string(fOutputDir.Data()), "png");
1599 Cbm::SaveCanvasAsImage(can2, string(fOutputDir.Data()), "png");
1600}
1601*/
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
Histogram manager.
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
Definition CbmMCTrack.h:67
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:170
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
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)
virtual InitStatus Init()
Inherited from FairTask.
void FillRingTrackDistanceCorr(const CbmRichRing *richRing, const FairTrackParam *pTrack, const CbmMCTrack *mcTrack)
virtual void Exec(Option_t *option)
Inherited from FairTask.
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)
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)
virtual void Finish()
Inherited from FairTask.
bool IsMcPrimaryElectron(const CbmMCTrack *mctrack)
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.
float GetCenterX() const
float GetCenterY() const
float GetCenterX() const
Definition CbmRichRing.h:77
float GetCenterY() const
Definition CbmRichRing.h:78
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)
Definition CbmUtils.cxx:36
Hash for CbmL1LinkKey.