CbmRoot
Loading...
Searching...
No Matches
CbmRichCorrection.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
5// ---------- Original Headers ---------- //
6#include "CbmRichCorrection.h"
7
8#include "CbmDrawHist.h"
9#include "CbmRichHit.h"
10#include "FairRootManager.h"
11#include "TCanvas.h"
12#include "TClonesArray.h"
13#include "TF1.h"
14#include "TH1D.h"
15#include "TH2D.h"
16
17#include <Logger.h>
18
19#include <iostream>
20#include <vector>
21
22// ---------- Included Headers ---------- //
23#include "CbmMCTrack.h"
24#include "CbmRichConverter.h"
25#include "CbmRichPoint.h"
26#include "CbmRichRingFitterCOP.h"
27#include "CbmRichRingFitterEllipseTau.h"
28#include "CbmRichRingLight.h"
29#include "CbmTrackMatchNew.h"
30#include "CbmUtils.h"
31#include "FairTrackParam.h"
32#include "FairVolume.h"
33#include "TEllipse.h"
34#include "TGeoManager.h"
35
36#include <algorithm>
37#include <fstream>
38#include <iomanip>
39
40#include <stdlib.h>
41//#include <stdio.h>
42#include "CbmGlobalTrack.h"
43#include "CbmRichGeoManager.h"
44#include "CbmRichHitProducer.h"
45
46//#include "TLorentzVector.h"
47#include "TGeoSphere.h"
48#include "TVirtualMC.h"
49class TGeoNode;
50class TGeoVolume;
51class TGeoShape;
52class TGeoMatrix;
53
54#include <boost/assign/list_of.hpp>
55using boost::assign::list_of;
56#include "TStyle.h"
57
58#include <sstream>
59
60CbmRichCorrection::CbmRichCorrection()
61 : FairTask()
62 , fRichHits(NULL)
63 , fRichRings(NULL)
64 , fRichProjections(NULL)
65 , fRichMirrorPoints(NULL)
66 , fRichMCPoints(NULL)
67 , fMCTracks(NULL)
68 , fRichRingMatches(NULL)
69 , fRichRefPlanePoints(NULL)
70 , fRichPoints(NULL)
71 , fGlobalTracks(NULL)
72 , fHM(NULL)
73 ,
74 //fGP(),
75 fNumbAxis(0)
76 , fTile(0)
77 , fEventNum(0)
78 , fOutputDir("")
79 , fRunTitle("")
80 , fAxisRotTitle("")
81 , fDrawProjection(kTRUE)
82 , fIsMeanCenter(kFALSE)
83 , fIsReconstruction(kFALSE)
84 , fCopFit(NULL)
85 , fTauFit(NULL)
86 , fPhi()
87{
88}
89
90CbmRichCorrection::~CbmRichCorrection() {}
91
92InitStatus CbmRichCorrection::Init()
93{
94 FairRootManager* manager = FairRootManager::Instance();
95
96 fRichHits = (TClonesArray*) manager->GetObject("RichHit");
97 if (NULL == fRichHits) {
98 Fatal("CbmRichCorrection::Init", "No RichHit array !");
99 }
100
101 fRichRings = (TClonesArray*) manager->GetObject("RichRing");
102 if (NULL == fRichRings) {
103 Fatal("CbmRichCorrection::Init", "No RichRing array !");
104 }
105
106 fRichProjections = (TClonesArray*) manager->GetObject("RichProjection");
107 if (NULL == fRichProjections) {
108 Fatal("CbmRichCorrection::Init", "No RichProjection array !");
109 }
110
111 fRichMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
112 if (NULL == fRichMirrorPoints) {
113 Fatal("CbmRichCorrection::Init", "No RichMirrorPoints array !");
114 }
115
116 fRichMCPoints = (TClonesArray*) manager->GetObject("RichPoint");
117 if (NULL == fRichMCPoints) {
118 Fatal("CbmRichCorrection::Init", "No RichMCPoints array !");
119 }
120
121 fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
122 if (NULL == fMCTracks) {
123 Fatal("CbmRichCorrection::Init", "No MCTracks array !");
124 }
125
126 fRichRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
127 if (NULL == fRichRingMatches) {
128 Fatal("CbmRichCorrection::Init", "No RichRingMatches array !");
129 }
130
131 fRichRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
132 if (NULL == fRichRefPlanePoints) {
133 Fatal("CbmRichCorrection::Init", "No RichRefPlanePoint array !");
134 }
135
136 fRichPoints = (TClonesArray*) manager->GetObject("RichPoint");
137 if (NULL == fRichPoints) {
138 Fatal("CbmRichCorrection::Init", "No RichPoint array !");
139 }
140
141 fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
142 if (NULL == fGlobalTracks) {
143 Fatal("CbmRichCorrection::Init", "No GlobalTrack array!");
144 }
145
146 fCopFit = new CbmRichRingFitterCOP();
147 fTauFit = new CbmRichRingFitterEllipseTau();
148 CbmRichConverter::Init();
149
150 InitHistProjection();
151
152 return kSUCCESS;
153}
154
155void CbmRichCorrection::InitHistProjection()
156{
157 fHM = new CbmHistManager();
158 /* for (std::map<string,string>::iterator it=fPathsMap.begin(); it!=fPathsMap.end(); ++it) { // Initialize all the histograms, using map IDs as inputs.
159 string name = "fHMCPoints_" + it->second; // it->first gives the paths; it->second gives the ID.
160 fHM->Create2<TH2D>(name, name + ";X_Axis [];Y_Axis [];Entries", 2001, -100., 100.,2001, 60., 210.);
161 }
162*/
163 Double_t upperScaleLimit = 3.5, bin = 400.;
164 // fhDistance => fhDistanceCenterToExtrapolatedTrack.
165 fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrack",
166 "fhDistanceCenterToExtrapolatedTrack;Distance fitted "
167 "center to extrapolated track;Number of entries",
168 bin, 0., 2.);
169 fHM->Create1<TH1D>("fhDistanceCorrected", "fhDistanceCorrected;Distance a [cm];A.U.", bin, 0., upperScaleLimit);
170 fHM->Create1<TH1D>("fhDifferenceX", "fhDifferenceX;Difference in X (fitted center - extrapolated track);A.U.", bin,
171 0., upperScaleLimit);
172 fHM->Create1<TH1D>("fhDifferenceY", "fhDifferenceY;Difference in Y (fitted center - extrapolated track);A.U.", bin,
173 0., upperScaleLimit);
174
175 fHM->Create1<TH1D>("fhDistanceUncorrected", "fhDistanceUncorrected;Distance a [cm];A.U.", bin, 0., upperScaleLimit);
176 fHM->Create1<TH1D>("fhDifferenceXUncorrected", "fhDifferenceXUncorrected;Difference in X uncorrected [cm];A.U.", bin,
177 0., upperScaleLimit);
178 fHM->Create1<TH1D>("fhDifferenceYUncorrected", "fhDifferenceYUncorrected;Difference in Y uncorrected [cm];A.U.", bin,
179 0., upperScaleLimit);
180
181 fHM->Create1<TH1D>("fhDistanceIdeal", "fhDistanceIdeal;Distance a [cm];A.U.", bin, 0., upperScaleLimit);
182 fHM->Create1<TH1D>("fhDifferenceXIdeal", "fhDifferenceXIdeal;Difference in X ideal [cm];A.U.", bin, 0.,
183 upperScaleLimit);
184 fHM->Create1<TH1D>("fhDifferenceYIdeal", "fhDifferenceYIdeal;Difference in Y ideal [cm];A.U.", bin, 0.,
185 upperScaleLimit);
186
187 fHM->Create1<TH1D>("fHistoDiffX",
188 "fHistoDiffX;Histogram difference between corrected and "
189 "ideal X positions;A.U.",
190 bin, 0., upperScaleLimit);
191 fHM->Create1<TH1D>("fHistoDiffY",
192 "fHistoDiffY;Histogram difference between corrected and "
193 "ideal Y positions;A.U.",
194 bin, 0., upperScaleLimit);
195
196 fHM->Create1<TH1D>("fHistoBoA", "fHistoBoA;Histogram B axis over A axis;A.U.", bin, 0., upperScaleLimit);
197}
198
199void CbmRichCorrection::Exec(Option_t* /*option*/)
200{
201 cout << endl
202 << "//"
203 "--------------------------------------------------------------------"
204 "------------------------------------------//"
205 << endl;
206 cout << "//---------------------------------------- CbmRichCorrection - EXEC "
207 "Function ----------------------------------------//"
208 << endl;
209 cout << "//"
210 "--------------------------------------------------------------------"
211 "--------------------------------------------------//"
212 << endl;
213 fEventNum++;
214 //LOG(debug2) << "CbmRichCorrection : Event #" << fEventNum;
215 cout << "CbmRichCorrection : Event #" << fEventNum << endl;
216
217 Int_t nofRingsInEvent = fRichRings->GetEntriesFast();
218 Int_t nofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
219 Int_t nofHitsInEvent = fRichHits->GetEntriesFast();
220 Int_t NofMCPoints = fRichMCPoints->GetEntriesFast();
221 Int_t NofMCTracks = fMCTracks->GetEntriesFast();
222 cout << "Nb of rings in evt = " << nofRingsInEvent << ", nb of mirror points = " << nofMirrorPoints
223 << ", nb of hits in evt = " << nofHitsInEvent << ", nb of Monte-Carlo points = " << NofMCPoints
224 << " and nb of Monte-Carlo tracks = " << NofMCTracks << endl
225 << endl;
226
227 TClonesArray* projectedPoint;
228
229 if (nofRingsInEvent == 0) {
230 cout << "Error no rings registered in event." << endl << endl;
231 }
232 else {
233 ProjectionProducer();
234 }
235}
236
237void CbmRichCorrection::ProjectionProducer()
238{
239 cout << "//------------------------------ CbmRichCorrection: Projection "
240 "Producer ------------------------------//"
241 << endl
242 << endl;
243
244 Int_t NofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
245 Int_t NofRingsInEvent = fRichRings->GetEntriesFast();
246 Int_t NofGTracks = fGlobalTracks->GetEntriesFast();
247 Int_t NofRefPlanePoints = fRichRefPlanePoints->GetEntriesFast();
248 Int_t NofPMTPoints = fRichPoints->GetEntriesFast();
249
250 // Declaration of points coordinates.
251 Double_t sphereRadius = 300., constantePMT = 0.;
252 vector<Double_t> ptM(3), ptMNew(3), ptC(3), ptCNew(3), ptR1(3), momR1(3), normalPMT(3), ptR2Mirr(3), ptR2Center(3),
253 ptPMirr(3), ptPR2(3), ptTileCenter(3);
254 vector<Double_t> ptCIdeal(3), ptR2CenterUnCorr(3), ptR2CenterIdeal(3), ptR2MirrUnCorr(3), ptR2MirrIdeal(3),
255 ptPMirrUnCorr(3), ptPMirrIdeal(3), ptPR2UnCorr(3), ptPR2Ideal(3);
256 Double_t reflectedPtCooVectSphereUnity[] = {0., 0., 0.};
257 TVector3 outPos, outPosUnCorr, outPosIdeal;
258 // Declaration of ring parameters.
259 Double_t ringCenter[] = {0., 0., 0.}, distToExtrapTrackHit = 0., distToExtrapTrackHitInPlane = 0.;
260 //Declarations related to geometry.
261 Int_t mirrTrackID = -1, pmtTrackID = -1, refPlaneTrackID = -1, motherID = -100, pmtMotherID = -100;
262 CbmMCTrack* track = NULL;
263 TGeoNavigator* navi;
264 TGeoNode* mirrNode;
265 TGeoMatrix *mirrMatrix, *pmtMatrix, *richMatrix;
266
267 CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
268 Double_t pmtPlaneX = gp->fPmt.fPlaneX;
269 Double_t pmtPlaneY = gp->fPmt.fPlaneY;
270 Double_t pmtWidth = gp->fPmt.fWidth;
271 Double_t pmtHeight = gp->fPmt.fHeight;
272
273 GetPmtNormal(NofPMTPoints, normalPMT, constantePMT);
274 cout << "Calculated normal vector to PMT plane = {" << normalPMT.at(0) << ", " << normalPMT.at(1) << ", "
275 << normalPMT.at(2) << "} and constante d = " << constantePMT << endl
276 << endl;
277
278 for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
279 //cout << "NofMirrorPoints = " << NofMirrorPoints << " and iMirr = " << iMirr << endl;
280 CbmRichPoint* mirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
281 mirrTrackID = mirrPoint->GetTrackID();
282 //cout << "Mirror track ID = " << mirrTrackID << endl;
283 if (mirrTrackID <= -1) {
284 cout << "Mirror track ID <= 1 !!!" << endl;
285 cout << "----------------------------------- End of loop N°" << iMirr + 1
286 << " on the mirror points. -----------------------------------" << endl
287 << endl;
288 continue;
289 }
290 track = (CbmMCTrack*) fMCTracks->At(mirrTrackID);
291 motherID = track->GetMotherId();
292 if (motherID == -1) {
293 //cout << "Mirror motherID == -1 !!!" << endl << endl;
294 ptM.at(0) = mirrPoint->GetX(), ptM.at(1) = mirrPoint->GetY(), ptM.at(2) = mirrPoint->GetZ();
295 //cout << "Mirror Point coordinates; x = " << ptM.at(0) << ", y = " << ptM.at(1) << " and z = " << ptM.at(2) << endl;
296 mirrNode = gGeoManager->FindNode(ptM.at(0), ptM.at(1), ptM.at(2));
297 if (mirrNode) {
298 //cout << "Mirror node found! Mirror node name = " << mirrNode->GetName() << endl;
299 navi = gGeoManager->GetCurrentNavigator();
300 cout << "Navigator path: " << navi->GetPath() << endl;
301 cout << "Coordinates of sphere center: " << endl;
302 navi->GetCurrentMatrix()->Print();
303 if (fIsMeanCenter)
304 GetMeanSphereCenter(navi,
305 ptC); //IF NO INFORMATION ON MIRRORS ARE KNOWN (TO BE USED IN RECONSTRUCTION STEP) !!!
306 else {
307 ptCIdeal.at(0) = navi->GetCurrentMatrix()->GetTranslation()[0];
308 ptCIdeal.at(1) = navi->GetCurrentMatrix()->GetTranslation()[1];
309 ptCIdeal.at(2) = navi->GetCurrentMatrix()->GetTranslation()[2];
310 }
311 cout << "Coordinates of tile center: " << endl;
312 navi->GetMotherMatrix()->Print();
313 ptC.at(0) = 0., ptC.at(1) = 132.594000, ptC.at(2) = 54.267226;
314 cout << "Sphere center coordinates of the aligned mirror tile, ideal = {" << ptCIdeal.at(0) << ", "
315 << ptCIdeal.at(1) << ", " << ptCIdeal.at(2) << "}" << endl;
316 cout << "Sphere center coordinates of the rotated mirror tile, w/ "
317 "GeoManager, = {"
318 << ptC.at(0) << ", " << ptC.at(1) << ", " << ptC.at(2) << "} and sphere inner radius = " << sphereRadius
319 << endl
320 << endl;
321 //ptCNew = RotateSphereCenter(ptTileCenter, ptC, navi);
322
323 for (Int_t iRefl = 0; iRefl < NofRefPlanePoints; iRefl++) {
324 //new((*projectedPoint)[iRefl]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
325 CbmRichPoint* refPlanePoint = (CbmRichPoint*) fRichRefPlanePoints->At(iRefl);
326 refPlaneTrackID = refPlanePoint->GetTrackID();
327 //cout << "Reflective plane track ID = " << refPlaneTrackID << endl;
328 if (mirrTrackID == refPlaneTrackID) {
329 //cout << "IDENTICAL TRACK ID FOUND !!!" << endl << endl;
330 ptR1.at(0) = refPlanePoint->GetX(), ptR1.at(1) = refPlanePoint->GetY(), ptR1.at(2) = refPlanePoint->GetZ();
331 momR1.at(0) = refPlanePoint->GetPx(), momR1.at(1) = refPlanePoint->GetPy(),
332 momR1.at(2) = refPlanePoint->GetPz();
333 cout << "Reflective Plane Point coordinates = {" << ptR1.at(0) << ", " << ptR1.at(1) << ", " << ptR1.at(2)
334 << "}" << endl;
335 cout << "And reflective Plane Point momenta = {" << momR1.at(0) << ", " << momR1.at(1) << ", "
336 << momR1.at(2) << "}" << endl;
337 cout << "Mirror Point coordinates = {" << ptM.at(0) << ", " << ptM.at(1) << ", " << ptM.at(2) << "}"
338 << endl;
339 CalculateMirrorIntersection(ptM, ptCIdeal, ptMNew);
340
341 if (fIsMeanCenter) {
342 GetMirrorIntersection(ptM, ptR1, momR1, ptC, sphereRadius);
343 // From ptM: how to retrieve tile ID ???
344 // => Compare distance of ptM to tile centers
345 }
346
347 ComputeR2(ptR2CenterUnCorr, ptR2MirrUnCorr, ptM, ptC, ptR1, navi, "Uncorrected");
348 ComputeR2(ptR2Center, ptR2Mirr, ptM, ptC, ptR1, navi, "Corrected");
349 ComputeR2(ptR2CenterIdeal, ptR2MirrIdeal, ptM, ptCIdeal, ptR1, navi, "Uncorrected");
350
351 ComputeP(ptPMirrUnCorr, ptPR2UnCorr, normalPMT, ptM, ptR2MirrUnCorr, constantePMT);
352 ComputeP(ptPMirr, ptPR2, normalPMT, ptM, ptR2Mirr, constantePMT);
353 ComputeP(ptPMirrIdeal, ptPR2Ideal, normalPMT, ptM, ptR2MirrIdeal, constantePMT);
354
355 TVector3 inPosUnCorr(ptPMirrUnCorr.at(0), ptPMirrUnCorr.at(1), ptPMirrUnCorr.at(2));
356 CbmRichGeoManager::GetInstance().RotatePoint(&inPosUnCorr, &outPosUnCorr);
357 cout << endl
358 << "New mirror points coordinates = {" << outPosUnCorr.x() << ", " << outPosUnCorr.y() << ", "
359 << outPosUnCorr.z() << "}" << endl;
360 TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
361 CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
362 cout << "New mirror points coordinates = {" << outPos.x() << ", " << outPos.y() << ", " << outPos.z() << "}"
363 << endl;
364 TVector3 inPosIdeal(ptPMirrIdeal.at(0), ptPMirrIdeal.at(1), ptPMirrIdeal.at(2));
365 CbmRichGeoManager::GetInstance().RotatePoint(&inPosIdeal, &outPosIdeal);
366 cout << "New mirror points coordinates = {" << outPosIdeal.x() << ", " << outPosIdeal.y() << ", "
367 << outPosIdeal.z() << "}" << endl
368 << endl;
369
370 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
371 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
372 pmtTrackID = pmtPoint->GetTrackID();
373 CbmMCTrack* track3 = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
374 pmtMotherID = track3->GetMotherId();
375 //cout << "pmt mother ID = " << pmtMotherID << endl;
376 if (pmtMotherID == mirrTrackID) {
377 ptP[0] = pmtPoint->GetX(), ptP[1] = pmtPoint->GetY(), ptP[2] = pmtPoint->GetZ();
378 //cout << "Identical mirror track ID and PMT mother ID !!!" << endl;
379 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
380 }
381 }
382 cout << "Looking for PMT hits: end." << endl << endl;*/
383 }
384 //else { cout << "No identical track ID between mirror point and reflective plane point found ..." << endl << endl; }
385 }
386
387 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
388 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
389 cout << "PMT points = {" << pmtPoint->GetX() << ", " << pmtPoint->GetY() << ", " << pmtPoint->GetZ() << "}" << endl;
390 TVector3 inputPoint(pmtPoint->GetX(), pmtPoint->GetY(), pmtPoint->GetZ());
391 TVector3 outputPoint;
392 CbmRichHitProducer::TiltPoint(&inputPoint, &outputPoint, fGP.fPmtPhi, fGP.fPmtTheta, fGP.fPmtZOrig);
393 cout << "New PMT points after rotation of the PMT plane: " << endl;
394 cout << outputPoint.X() << "\t" << outputPoint.Y() << "\t" << outputPoint.Z() << endl;
395 }*/
396
397 FillHistProjection(outPosIdeal, outPosUnCorr, outPos, NofGTracks, normalPMT, constantePMT);
398 }
399 else {
400 cout << "Not a mother particle ..." << endl;
401 }
402 cout << "----------------------------------- "
403 << "End of loop N°" << iMirr + 1 << " on the mirror points."
404 << " -----------------------------------" << endl
405 << endl;
406 }
407 }
408}
409
410void CbmRichCorrection::GetPmtNormal(Int_t NofPMTPoints, vector<Double_t>& normalPMT, Double_t& normalCste)
411{
412 //cout << endl << "//------------------------------ CbmRichCorrection: Calculate PMT Normal ------------------------------//" << endl << endl;
413
414 Int_t pmtTrackID, pmtMotherID;
415 Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0., scalarProd = 0.;
416 Double_t pmtPt[] = {0., 0., 0.};
417 Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
418 CbmMCTrack* track;
419
420 /*
421 * 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.
422 * Formula used is: vect(AB) x vect(AC) = normal.
423 * Normalize the normal vector obtained and check with three random points from the PMT plane, whether the scalar product is equal to zero.
424 */
425 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
426 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
427 pmtTrackID = pmtPoint->GetTrackID();
428 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
429 pmtMotherID = track->GetMotherId();
430 a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
431 //cout << "a[0] = " << a[0] << ", a[1] = " << a[1] << " et a[2] = " << a[2] << endl;
432 break;
433 }
434 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
435 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
436 pmtTrackID = pmtPoint->GetTrackID();
437 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
438 pmtMotherID = track->GetMotherId();
439 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
440 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
441 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
442 > 7) {
443 b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
444 //cout << "b[0] = " << b[0] << ", b[1] = " << b[1] << " et b[2] = " << b[2] << endl;
445 break;
446 }
447 }
448 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
449 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
450 pmtTrackID = pmtPoint->GetTrackID();
451 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
452 pmtMotherID = track->GetMotherId();
453 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
454 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
455 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
456 > 7
457 && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2) + TMath::Power(b[1] - pmtPoint->GetY(), 2)
458 + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
459 > 7) {
460 c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
461 //cout << "c[0] = " << c[0] << ", c[1] = " << c[1] << " et c[2] = " << c[2] << endl;
462 break;
463 }
464 }
465
466 k = (b[0] - a[0]) / (c[0] - a[0]);
467 if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0 || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
468 cout << "Error in normal calculation, vect_AB and vect_AC are collinear." << endl;
469 }
470 else {
471 buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
472 buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
473 buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
474 normalPMT.at(0) =
475 buffNormX / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
476 normalPMT.at(1) =
477 buffNormY / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
478 normalPMT.at(2) =
479 buffNormZ / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
480 }
481
482 CbmRichPoint* pmtPoint1 = (CbmRichPoint*) fRichPoints->At(20);
483 scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
484 + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
485 //cout << "1st scalar product between vectAM and normale = " << scalarProd << endl;
486 // 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.
487 normalCste =
488 -1
489 * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY() + normalPMT.at(2) * pmtPoint1->GetZ());
490 CbmRichPoint* pmtPoint2 = (CbmRichPoint*) fRichPoints->At(15);
491 scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
492 + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
493 //cout << "2nd scalar product between vectAM and normale = " << scalarProd << endl;
494 CbmRichPoint* pmtPoint3 = (CbmRichPoint*) fRichPoints->At(25);
495 scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
496 + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
497 //cout << "3nd scalar product between vectAM and normale = " << scalarProd << endl;
498}
499
500void CbmRichCorrection::GetMeanSphereCenter(TGeoNavigator* navi, vector<Double_t>& ptC)
501{
502 const Char_t* topNodePath;
503 topNodePath = gGeoManager->GetTopNode()->GetName();
504 cout << "Top node path: " << topNodePath << endl;
505 TGeoVolume* rootTop;
506 rootTop = gGeoManager->GetTopVolume();
507 rootTop->Print();
508
509 TGeoIterator nextNode(rootTop);
510 TGeoNode* curNode;
511 const TGeoMatrix* curMatrix;
512 const Double_t* curNodeTranslation; // 3 components - pointers to some memory which is provided by ROOT
513 const Double_t* curNodeRotationM; // 9 components - pointers to some memory which is provided by ROOT
514 TString filterName0("mirror_tile_type0");
515 TString filterName1("mirror_tile_type1");
516 TString filterName2("mirror_tile_type2");
517 TString filterName3("mirror_tile_type3");
518 TString filterName4("mirror_tile_type4");
519 TString filterName5("mirror_tile_type5");
520 Int_t counter = 0;
521 Double_t sphereXTot = 0., sphereYTot = 0., sphereZTot = 0.;
522
523 while ((curNode = nextNode())) {
524 TString nodeName(curNode->GetName());
525 TString nodePath;
526
527 // Filter using volume name, not node name
528 // But you can do 'if (nodeName.Contains("filter"))'
529 if (curNode->GetVolume()->GetName() == filterName0 || curNode->GetVolume()->GetName() == filterName1
530 || curNode->GetVolume()->GetName() == filterName2 || curNode->GetVolume()->GetName() == filterName3
531 || curNode->GetVolume()->GetName() == filterName4 || curNode->GetVolume()->GetName() == filterName5) {
532 if (curNode->GetNdaughters() == 0) {
533 // All deepest nodes of mirror tiles here (leaves)
534 // Thus we get spherical surface centers
535 nextNode.GetPath(nodePath);
536 curMatrix = nextNode.GetCurrentMatrix();
537 curNodeTranslation = curMatrix->GetTranslation();
538 curNodeRotationM = curMatrix->GetRotationMatrix();
539 printf("%s tr:\t", nodePath.Data());
540 printf("%08f\t%08f\t%08f\t\n", curNodeTranslation[0], curNodeTranslation[1], curNodeTranslation[2]);
541 if (curNodeTranslation[1] > 0) { // CONDITION FOR UPPER MIRROR WALL STUDY
542 sphereXTot += curNodeTranslation[0];
543 sphereYTot += curNodeTranslation[1];
544 sphereZTot += curNodeTranslation[2];
545 counter++;
546 }
547 }
548 }
549 }
550 ptC.at(0) = sphereXTot / counter;
551 ptC.at(1) = sphereYTot / counter;
552 ptC.at(2) = sphereZTot / counter;
553
554 counter = 0;
555 nextNode.Reset();
556}
557
558void CbmRichCorrection::GetMirrorIntersection(vector<Double_t>& ptM, vector<Double_t> ptR1, vector<Double_t> momR1,
559 vector<Double_t> ptC, Double_t sphereRadius)
560{
561 Double_t a = 0., b = 0., c = 0., d = 0., k0 = 0., k1 = 0., k2 = 0.;
562
563 a = TMath::Power(momR1.at(0), 2) + TMath::Power(momR1.at(1), 2) + TMath::Power(momR1.at(2), 2);
564 b = 2
565 * (momR1.at(0) * (ptR1.at(0) - ptC.at(0)) + momR1.at(1) * (ptR1.at(1) - ptC.at(1))
566 + momR1.at(2) * (ptR1.at(2) - ptC.at(2)));
567 c = TMath::Power(ptR1.at(0) - ptC.at(0), 2) + TMath::Power(ptR1.at(1) - ptC.at(1), 2)
568 + TMath::Power(ptR1.at(2) - ptC.at(2), 2) - TMath::Power(sphereRadius, 2);
569 d = b * b - 4 * a * c;
570 cout << "d = " << d << endl;
571
572 if (d < 0) {
573 cout << "Error no solution to degree 2 equation found ; discriminant below 0." << endl;
574 ptM.at(0) = 0., ptM.at(1) = 0., ptM.at(2) = 0.;
575 }
576 else if (d == 0) {
577 cout << "One solution to degree 2 equation found." << endl;
578 k0 = -b / (2 * a);
579 ptM.at(0) = ptR1.at(0) + k0 * momR1.at(0);
580 ptM.at(1) = ptR1.at(1) + k0 * momR1.at(1);
581 ptM.at(2) = ptR1.at(2) + k0 * momR1.at(2);
582 }
583 else if (d > 0) {
584 cout << "Two solutions to degree 2 equation found." << endl;
585 k1 = ((-b - TMath::Sqrt(d)) / (2 * a));
586 k2 = ((-b + TMath::Sqrt(d)) / (2 * a));
587
588 if (ptR1.at(2) + k1 * momR1.at(2) > ptR1.at(2) + k2 * momR1.at(2)) {
589 ptM.at(0) = ptR1.at(0) + k1 * momR1.at(0);
590 ptM.at(1) = ptR1.at(1) + k1 * momR1.at(1);
591 ptM.at(2) = ptR1.at(2) + k1 * momR1.at(2);
592 }
593 else if (ptR1.at(2) + k1 * momR1.at(2) < ptR1.at(2) + k2 * momR1.at(2)) {
594 ptM.at(0) = ptR1.at(0) + k2 * momR1.at(0);
595 ptM.at(1) = ptR1.at(1) + k2 * momR1.at(1);
596 ptM.at(2) = ptR1.at(2) + k2 * momR1.at(2);
597 }
598 }
599}
600
601vector<Double_t> CbmRichCorrection::RotateSphereCenter(vector<Double_t> ptM, vector<Double_t> ptC, TGeoNavigator* navi)
602{
603 vector<Double_t> ptCNew(3), ptCNew2(3), ptCNew3(3);
604 Double_t cosPhi = 0., sinPhi = 0., cosTheta = 0., sinTheta = 0., phi2 = 0., theta2 = 0.;
605 Double_t diff[3], transfoMat[3][3], invMat[3][3], corrMat[3][3], buff1[3][3], buff2[3][3], buff3[3][3], buff4[3][3],
606 buff5[3][3], RotX[3][3], RotY[3][3];
607 Double_t corrMat2[3][3], RotX2[3][3], RotY2[3][3];
608 InvertMatrix(transfoMat, invMat, navi);
609 /*for (Int_t i=0; i<3; i++) {
610 for (Int_t j=0; j<3; j++) {
611 cout << invMat[i][j] << "\t";
612 }
613 cout << endl;
614 }*/
615
616 // Reading misalignment information from correction_param.txt text file.
617 vector<Double_t> outputFit(4);
618 ifstream corr_file;
619 corr_file.open("correction_param.txt");
620 if (corr_file.is_open()) {
621 for (Int_t i = 0; i < 4; i++) {
622 corr_file >> outputFit.at(i);
623 }
624 corr_file.close();
625 }
626 else {
627 cout << "Error in CbmRichCorrection: unable to open parameter file!" << endl << endl;
628 sleep(5);
629 }
630 cout << "Misalignment parameters read from file = [" << outputFit.at(0) << " ; " << outputFit.at(1) << " ; "
631 << outputFit.at(2) << " ; " << outputFit.at(3) << "]" << endl;
632
633 // Initializing the matrices used for further calculations.
634 for (Int_t i = 0; i < 3; i++) {
635 for (Int_t j = 0; j < 3; j++) {
636 corrMat[i][j] = 0.;
637 corrMat2[i][j] = 0.;
638 buff1[i][j] = 0.;
639 buff2[i][j] = 0.;
640 buff3[i][j] = 0.;
641 buff4[i][j] = 0.;
642 buff5[i][j] = 0.;
643 RotX[i][j] = 0.;
644 RotX2[i][j] = 0.;
645 RotY[i][j] = 0.;
646 RotY2[i][j] = 0.;
647 }
648 ptCNew.at(i) = 0.;
649 ptCNew2.at(i) = 0.;
650 ptCNew3.at(i) = 0.;
651 }
652
653 //Initializing the cosine and sine functions, using inputs from text file. Phi (resp. Theta) angle corresponds to calculated rotation around X (resp. Y) axis.
654 cosPhi = TMath::Cos(outputFit.at(0));
655 sinPhi = TMath::Sin(outputFit.at(0));
656 cosTheta = TMath::Cos(outputFit.at(1));
657 sinTheta = TMath::Sin(outputFit.at(1));
658
659 //Initializing the rotation matrices for rotations around X and Y axes.
660 // Y towards Z is the rotation direction defined in the tile frame, which is the same definition as for the global
661 // frame. But as the X axis direction is opposite in the two frames, to correct for misalignment, the rotation
662 // direction remains unchanged in the global frame.
663 // So define rotation around X axis, with Y towards Z, in the global frame:
664 RotX[0][0] = 1;
665 RotX[0][1] = 0;
666 RotX[0][2] = 0;
667 RotX[1][0] = 0;
668 RotX[1][1] = cosPhi;
669 RotX[1][2] = -1 * sinPhi;
670 RotX[2][0] = 0;
671 RotX[2][1] = sinPhi;
672 RotX[2][2] = cosPhi;
673 // X towards Z is the rotation direction defined in the tile frame, which is the same definition as for the global
674 // frame. And as the Y axes direction in the tile and global frames are identical, the rotation direction for the
675 // correction is the opposite in the global frame - in order to correct for the misalignment.
676 // So define rotation around Y axis, with Z towards X, in the global frame:
677 RotY[0][0] = cosTheta;
678 RotY[0][1] = 0;
679 RotY[0][2] = sinTheta;
680 RotY[1][0] = 0;
681 RotY[1][1] = 1;
682 RotY[1][2] = 0;
683 RotY[2][0] = -1 * sinTheta;
684 RotY[2][1] = 0;
685 RotY[2][2] = cosTheta;
686
687 // Translation of the sphere center of the tile: S.
688 for (Int_t i = 0; i < 3; i++) {
689 diff[i] = ptC.at(i) - ptM.at(i);
690 }
691
692 // Calculation of the correction matrix, from rotation matrices.
693 for (Int_t i = 0; i < 3; i++) {
694 for (Int_t j = 0; j < 3; j++) {
695 //corrMat[i][j] = RotZ[i][j];
696 for (Int_t k = 0; k < 3; k++) {
697 corrMat[i][j] = RotY[i][k] * RotX[k][j] + corrMat[i][j];
698 }
699 }
700 }
701 // Calculation of the new coordinates of the translated point S: newS = transfoMat*corrMat*invMat*diff.
702 for (Int_t i = 0; i < 3; i++) {
703 for (Int_t j = 0; j < 3; j++) {
704 for (Int_t k = 0; k < 3; k++) {
705 buff1[i][j] = corrMat[i][k] * invMat[k][j] + buff1[i][j];
706 }
707 }
708 }
709 for (Int_t i = 0; i < 3; i++) {
710 for (Int_t j = 0; j < 3; j++) {
711 for (Int_t k = 0; k < 3; k++) {
712 buff2[i][j] = transfoMat[i][k] * buff1[k][j] + buff2[i][j];
713 }
714 }
715 }
716 for (Int_t i = 0; i < 3; i++) {
717 for (Int_t j = 0; j < 3; j++) {
718 ptCNew.at(i) = buff2[i][j] * diff[j] + ptCNew.at(i);
719 //ptCNew.at(i) = invMat[i][j]*diff[j] + ptCNew.at(i);
720 }
721 }
722
723 // Calculating the theoretical rotation angles to be applied to the translated point S, to get the point along the Z axis (should obtain a {0; 0; -300} coo).
724 phi2 = TMath::ATan2(diff[1], -1 * diff[2]);
725 theta2 = TMath::ATan2(diff[0], -1 * diff[2]);
726 cout << "Calculated Phi (= arctan(y/-z)), in degrees: " << TMath::RadToDeg() * phi2
727 << " and calculated Theta (= arctan(x/-z)), in degrees: " << TMath::RadToDeg() * theta2 << endl;
728 // Defining the rotation matrices accordingly:
729 // Rotation around X axis, with Z towards Y.
730 RotX2[0][0] = 1;
731 RotX2[0][1] = 0;
732 RotX2[0][2] = 0;
733 RotX2[1][0] = 0;
734 RotX2[1][1] = TMath::Cos(phi2);
735 RotX2[1][2] = TMath::Sin(phi2);
736 RotX2[2][0] = 0;
737 RotX2[2][1] = -1 * TMath::Sin(phi2);
738 RotX2[2][2] = TMath::Cos(phi2);
739 // Rotation around Y axis, with Z towards X.
740 RotY2[0][0] = TMath::Cos(theta2);
741 RotY2[0][1] = 0;
742 RotY2[0][2] = TMath::Sin(theta2);
743 RotY2[1][0] = 0;
744 RotY2[1][1] = 1;
745 RotY2[1][2] = 0;
746 RotY2[2][0] = -1 * TMath::Sin(theta2);
747 RotY2[2][1] = 0;
748 RotY2[2][2] = TMath::Cos(theta2);
749
750 // Calculation of the correction matrix, from new rotation matrices.
751 for (Int_t i = 0; i < 3; i++) {
752 for (Int_t j = 0; j < 3; j++) {
753 for (Int_t k = 0; k < 3; k++) {
754 corrMat2[i][j] = RotY2[i][k] * RotX2[k][j] + corrMat2[i][j];
755 }
756 }
757 }
758 // Calculation of the new coordinates of the translated point S: newS = transfoMat*corrMat2*invMat*diff.
759 for (Int_t i = 0; i < 3; i++) {
760 for (Int_t j = 0; j < 3; j++) {
761 for (Int_t k = 0; k < 3; k++) {
762 buff3[i][j] = corrMat2[i][k] * invMat[k][j] + buff3[i][j];
763 }
764 }
765 }
766 for (Int_t i = 0; i < 3; i++) {
767 for (Int_t j = 0; j < 3; j++) {
768 for (Int_t k = 0; k < 3; k++) {
769 buff4[i][j] = transfoMat[i][k] * buff3[k][j] + buff4[i][j];
770 }
771 }
772 }
773 for (Int_t i = 0; i < 3; i++) {
774 for (Int_t j = 0; j < 3; j++) {
775 for (Int_t k = 0; k < 3; k++) {
776 buff5[i][j] = RotX[i][k] * invMat[k][j] + buff5[i][j];
777 //buff5[i][j] = RotY[i][k]*invMat[k][j] + buff5[i][j];
778 }
779 }
780 }
781 for (Int_t i = 0; i < 3; i++) {
782 for (Int_t j = 0; j < 3; j++) {
783 ptCNew2.at(i) = buff4[i][j] * diff[j] + ptCNew2.at(i);
784 //ptCNew3.at(i) = buff3[i][j]*diff[j] + ptCNew3.at(i);
785 ptCNew3.at(i) = buff5[i][j] * diff[j] + ptCNew3.at(i);
786 }
787 }
788
789 for (Int_t i = 0; i < 3; i++) {
790 ptCNew.at(i) += ptM.at(i);
791 ptCNew2.at(i) += ptM.at(i);
792 ptCNew3.at(i) += ptM.at(i);
793 }
794 cout << "diff = {" << diff[0] << ", " << diff[1] << ", " << diff[2] << "}" << endl;
795 cout << "New coordinates of the rotated tile sphere center (using angles "
796 "from text file) = {"
797 << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
798 cout << "New coordinates of the rotated tile sphere center (using calculated "
799 "angles) = {"
800 << ptCNew2.at(0) << ", " << ptCNew2.at(1) << ", " << ptCNew2.at(2) << "}" << endl;
801 cout << "Tile coordinates after translation, invMat and rotations around X "
802 "and Y axes = {"
803 << ptCNew3.at(0) << ", " << ptCNew3.at(1) << ", " << ptCNew3.at(2) << "}" << endl
804 << endl;
805
806 return ptCNew;
807}
808
809void CbmRichCorrection::InvertMatrix(Double_t mat[3][3], Double_t invMat[3][3], TGeoNavigator* navi)
810{
811 Double_t deter = 0., det11 = 0., det12 = 0., det13 = 0., det21 = 0., det22 = 0., det23 = 0., det31 = 0., det32 = 0.,
812 det33 = 0., buff[3][3], prodMat[3][3];
813
814 //Filling the transformation matrix of the tile.
815 // STANDARD FILL:
816 mat[0][0] = navi->GetMotherMatrix()->GetRotationMatrix()[0];
817 mat[0][1] = navi->GetMotherMatrix()->GetRotationMatrix()[1];
818 mat[0][2] = navi->GetMotherMatrix()->GetRotationMatrix()[2];
819 mat[1][0] = navi->GetMotherMatrix()->GetRotationMatrix()[3];
820 mat[1][1] = navi->GetMotherMatrix()->GetRotationMatrix()[4];
821 mat[1][2] = navi->GetMotherMatrix()->GetRotationMatrix()[5];
822 mat[2][0] = navi->GetMotherMatrix()->GetRotationMatrix()[6];
823 mat[2][1] = navi->GetMotherMatrix()->GetRotationMatrix()[7];
824 mat[2][2] = navi->GetMotherMatrix()->GetRotationMatrix()[8];
825
826 /* // TEST FILL:
827 mat[0][0] = navi->GetMotherMatrix()->GetRotationMatrix()[2];
828 mat[0][1] = navi->GetMotherMatrix()->GetRotationMatrix()[1];
829 mat[0][2] = navi->GetMotherMatrix()->GetRotationMatrix()[0];
830 mat[1][0] = navi->GetMotherMatrix()->GetRotationMatrix()[5];
831 mat[1][1] = navi->GetMotherMatrix()->GetRotationMatrix()[4];
832 mat[1][2] = navi->GetMotherMatrix()->GetRotationMatrix()[3];
833 mat[2][0] = navi->GetMotherMatrix()->GetRotationMatrix()[8];
834 mat[2][1] = navi->GetMotherMatrix()->GetRotationMatrix()[7];
835 mat[2][2] = navi->GetMotherMatrix()->GetRotationMatrix()[6];
836*/
837 /*for (Int_t i=0; i<3; i++) {
838 for (Int_t j=0; j<3; j++) {
839 mat[i][j] = navi->GetMotherMatrix()->GetRotationMatrix()[i*3+j];
840 }
841 }
842 cout << "Matrix index [2][0] = " << mat[2][0] << endl;*/
843
844 //Computing the inverse of the matrix: inv = (1/det)*adjugate(mat) = (1/det)*transpose(cofactor_matrix(mat)).
845 deter = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1])
846 - mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0])
847 + mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
848 if (deter == 0) {
849 cout << "Error in CbmRichCorrection::InvertMatrix; determinant of input "
850 "matrix equals to zero !!!"
851 << endl;
852 sleep(5);
853 for (Int_t i = 0; i < 3; i++) {
854 for (Int_t j = 0; j < 3; j++) {
855 invMat[i][j] = 0;
856 }
857 }
858 }
859 else {
860 buff[0][0] = TMath::Power(-1, 2) * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]);
861 buff[0][1] = TMath::Power(-1, 3) * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]);
862 buff[0][2] = TMath::Power(-1, 4) * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
863 buff[1][0] = TMath::Power(-1, 3) * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]);
864 buff[1][1] = TMath::Power(-1, 4) * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]);
865 buff[1][2] = TMath::Power(-1, 5) * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
866 buff[2][0] = TMath::Power(-1, 4) * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
867 buff[2][1] = TMath::Power(-1, 5) * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]);
868 buff[2][2] = TMath::Power(-1, 6) * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
869
870 for (Int_t i = 0; i < 3; i++) {
871 for (Int_t j = 0; j < 3; j++) {
872 invMat[i][j] = buff[i][j] / deter;
873 prodMat[i][j] = 0;
874 }
875 }
876
877 /* for (Int_t i=0; i<3; i++) {
878 for (Int_t j=0; j<3; j++) {
879 for (Int_t k=0; k<3; k++) {
880 prodMat[i][j] = mat[i][k]*invMat[k][j] + prodMat[i][j];
881 }
882 }
883 }
884
885 cout << "Matrix calculation to check whether inverse matrix is correct:" << endl;
886 for (Int_t i=0; i<3; i++) {
887 for (Int_t j=0; j<3; j++) {
888 cout << "Resulting matrix = [" << prodMat[i][j] << "] \t";
889 }
890 cout << endl;
891 }*/
892 }
893}
894
895void CbmRichCorrection::CalculateMirrorIntersection(vector<Double_t> ptM, vector<Double_t> ptCIdeal,
896 vector<Double_t>& ptMNew)
897{
898 Double_t t = 0., diffX = 0., diffY = 0., diffZ = 0.;
899 diffX = ptM.at(0) - ptCIdeal.at(0);
900 diffY = ptM.at(1) - ptCIdeal.at(1);
901 diffZ = ptM.at(2) - ptCIdeal.at(2);
902 t = TMath::Sqrt(300 * 300 / (diffX * diffX + diffY * diffY + diffZ * diffZ));
903
904 ptMNew.at(0) = t * diffX + ptCIdeal.at(0);
905 ptMNew.at(1) = t * diffY + ptCIdeal.at(1);
906 ptMNew.at(2) = t * diffZ + ptCIdeal.at(2);
907 cout << "New coordinates of point M = {" << ptMNew.at(0) << ", " << ptMNew.at(1) << ", " << ptMNew.at(2) << "}"
908 << endl;
909}
910
911void CbmRichCorrection::ComputeR2(vector<Double_t>& ptR2Center, vector<Double_t>& ptR2Mirr, vector<Double_t> ptM,
912 vector<Double_t> ptC, vector<Double_t> ptR1, TGeoNavigator* navi, TString s)
913{
914 cout << endl
915 << "//------------------------------ CbmRichCorrection: ComputeR2 "
916 "------------------------------//"
917 << endl
918 << endl;
919
920 vector<Double_t> normalMirr(3), ptCNew(3), ptTileCenter(3);
921 Double_t t1 = 0., t2 = 0., t3 = 0.;
922
923 if (s == "Corrected") {
924 // Use the correction information from text file, to the tile sphere center:
925 // Reading misalignment information from correction_param.txt text file.
926 vector<Double_t> outputFit(4);
927 ifstream corr_file;
928 TString str = fOutputDir + "correction_param_array___" + fNumbAxis + fTile + ".txt";
929 corr_file.open(str);
930 if (corr_file.is_open()) {
931 for (Int_t i = 0; i < 4; i++) {
932 corr_file >> outputFit.at(i);
933 }
934 //for (Int_t i=0; i<2; i++) {corr_file >> outputFit.at(i);}
935 corr_file.close();
936 }
937 else {
938 cout << "Error in CbmRichCorrection: unable to open parameter file!" << endl;
939 cout << "Parameter file path = " << str << endl << endl;
940 sleep(5);
941 }
942 cout << "Misalignment parameters read from file = [" << outputFit.at(0) << " ; " << outputFit.at(1) << " ; "
943 << outputFit.at(2) << " ; " << outputFit.at(3) << "]" << endl;
944
945 //ptCNew.at(0) = TMath::Abs(ptC.at(0) - TMath::Abs(outputFit.at(3)));
946 //ptCNew.at(1) = TMath::Abs(ptC.at(1) - TMath::Abs(outputFit.at(2)));
947 ptCNew.at(0) = ptC.at(0) + outputFit.at(3);
948 ptCNew.at(1) = ptC.at(1) + outputFit.at(2);
949 ptCNew.at(2) = ptC.at(2);
950 ptTileCenter.at(0) = navi->GetMotherMatrix()->GetTranslation()[0];
951 ptTileCenter.at(1) = navi->GetMotherMatrix()->GetTranslation()[1];
952 ptTileCenter.at(2) = navi->GetMotherMatrix()->GetTranslation()[2];
953 cout << "Mirror tile center coordinates = {" << ptTileCenter.at(0) << ", " << ptTileCenter.at(1) << ", "
954 << ptTileCenter.at(2) << "}" << endl;
955 Double_t x = 0., y = 0., z = 0., dist = 0., dist2 = 0., z2 = 0.;
956 x = TMath::Power(ptCNew.at(0) - ptTileCenter.at(0), 2);
957 y = TMath::Power(ptCNew.at(1) - ptTileCenter.at(1), 2);
958 z = TMath::Power(ptCNew.at(2) - ptTileCenter.at(2), 2);
959 dist = TMath::Sqrt(x + y + z);
960 z2 = ptTileCenter.at(2) - TMath::Sqrt(TMath::Power(300, 2) - x - y) - ptCNew.at(2);
961 cout << "{x, y, z} = {" << x << ", " << y << ", " << z << "}, dist = " << dist << " and z2 = " << z2 << endl;
962 dist2 = TMath::Sqrt(x + y + TMath::Power(z2 - ptTileCenter.at(2), 2));
963 cout << "dist2 = " << dist2 << endl;
964 ptCNew.at(2) += z2;
965 cout << "Sphere center coordinates of the rotated mirror tile, after "
966 "correction, = {"
967 << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
968 }
969 else if (s == "Uncorrected") {
970 // Keep the same tile sphere center, with no correction information.
971 ptCNew = ptC;
972 cout << "Sphere center coordinates of the rotated mirror tile, without "
973 "correction = {"
974 << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
975 }
976 else {
977 cout << "No input given in function ComputeR2! Uncorrected parameters for "
978 "the sphere center of the tile will be used!"
979 << endl;
980 ptCNew = ptC;
981 cout << "Sphere center coordinates of the rotated mirror tile, without "
982 "correction = {"
983 << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
984 }
985
986 normalMirr.at(0) = (ptCNew.at(0) - ptM.at(0))
987 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
988 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
989 normalMirr.at(1) = (ptCNew.at(1) - ptM.at(1))
990 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
991 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
992 normalMirr.at(2) = (ptCNew.at(2) - ptM.at(2))
993 / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
994 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
995 //cout << "Calculated and normalized normal of mirror tile = {" << normalMirr.at(0) << ", " << normalMirr.at(1) << ", " << normalMirr.at(2) << "}" << endl;
996
997 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))
998 + (ptR1.at(2) - ptM.at(2)) * (ptCNew.at(2) - ptM.at(2)))
999 / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
1000 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
1001 ptR2Center.at(0) = 2 * (ptM.at(0) + t1 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
1002 ptR2Center.at(1) = 2 * (ptM.at(1) + t1 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
1003 ptR2Center.at(2) = 2 * (ptM.at(2) + t1 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
1004 t2 =
1005 ((ptR1.at(0) - ptCNew.at(0)) * (ptCNew.at(0) - ptM.at(0)) + (ptR1.at(1) - ptCNew.at(1)) * (ptCNew.at(1) - ptM.at(1))
1006 + (ptR1.at(2) - ptCNew.at(2)) * (ptCNew.at(2) - ptM.at(2)))
1007 / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2) + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
1008 + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
1009 ptR2Mirr.at(0) = 2 * (ptCNew.at(0) + t2 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
1010 ptR2Mirr.at(1) = 2 * (ptCNew.at(1) + t2 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
1011 ptR2Mirr.at(2) = 2 * (ptCNew.at(2) + t2 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
1012 /*//SAME AS calculation of t2 above
1013 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));
1014 reflectedPtCooVectSphereUnity[0] = 2*(ptCNew.at(0)+t3*(normalMirr.at(0)))-ptR1.at(0);
1015 reflectedPtCooVectSphereUnity[1] = 2*(ptCNew.at(1)+t3*(normalMirr.at(1)))-ptR1.at(1);
1016 reflectedPtCooVectSphereUnity[2] = 2*(ptCNew.at(2)+t3*(normalMirr.at(2)))-ptR1.at(2);*/
1017 cout << "Coordinates of point R2 on reflective plane after reflection on the "
1018 "mirror tile:"
1019 << endl;
1020 //cout << "* using mirror point M to define \U00000394: {" << ptR2Center.at(0) << ", " << ptR2Center.at(1) << ", " << ptR2Center.at(2) << "}" << endl;
1021 cout << "* using sphere center C to define \U00000394: {" << ptR2Mirr.at(0) << ", " << ptR2Mirr.at(1) << ", "
1022 << ptR2Mirr.at(2) << "}" << endl;
1023 //cout << "Ref Pt Coo using unity Mirror-Sphere vector & sphere pt = {" << reflectedPtCooVectSphereUnity[0] << ", " << reflectedPtCooVectSphereUnity[1] << ", " << reflectedPtCooVectSphereUnity[2] << "}" << endl << endl;
1024 //cout << "NofPMTPoints = " << NofPMTPoints << endl;
1025}
1026
1027void CbmRichCorrection::ComputeP(vector<Double_t>& ptPMirr, vector<Double_t>& ptPR2, vector<Double_t> normalPMT,
1028 vector<Double_t> ptM, vector<Double_t> ptR2Mirr, Double_t constantePMT)
1029{
1030 cout << endl
1031 << "//------------------------------ CbmRichCorrection: ComputeP "
1032 "------------------------------//"
1033 << endl
1034 << endl;
1035
1036 Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
1037
1038 k1 = -1
1039 * ((normalPMT.at(0) * ptM.at(0) + normalPMT.at(1) * ptM.at(1) + normalPMT.at(2) * ptM.at(2) + constantePMT)
1040 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0)) + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
1041 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
1042 ptPMirr.at(0) = ptM.at(0) + k1 * (ptR2Mirr.at(0) - ptM.at(0));
1043 ptPMirr.at(1) = ptM.at(1) + k1 * (ptR2Mirr.at(1) - ptM.at(1));
1044 ptPMirr.at(2) = ptM.at(2) + k1 * (ptR2Mirr.at(2) - ptM.at(2));
1045 k2 = -1
1046 * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1) + normalPMT.at(2) * ptR2Mirr.at(2)
1047 + constantePMT)
1048 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0)) + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
1049 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
1050 ptPR2.at(0) = ptR2Mirr.at(0) + k2 * (ptR2Mirr.at(0) - ptM.at(0));
1051 ptPR2.at(1) = ptR2Mirr.at(1) + k2 * (ptR2Mirr.at(1) - ptM.at(1));
1052 ptPR2.at(2) = ptR2Mirr.at(2) + k2 * (ptR2Mirr.at(2) - ptM.at(2));
1053 cout << "Coordinates of point P on PMT plane, after reflection on the mirror "
1054 "tile and extrapolation to the PMT plane:"
1055 << endl;
1056 cout << "* using mirror point M to define \U0001D49F ': {" << ptPMirr.at(0) << ", " << ptPMirr.at(1) << ", "
1057 << ptPMirr.at(2) << "}" << endl;
1058 //cout << "* using reflected point R2 to define \U0001D49F ': {" << ptPR2.at(0) << ", " << ptPR2.at(1) << ", " << ptPR2.at(2) << "}" << endl;
1059 checkCalc1 =
1060 ptPMirr.at(0) * normalPMT.at(0) + ptPMirr.at(1) * normalPMT.at(1) + ptPMirr.at(2) * normalPMT.at(2) + constantePMT;
1061 cout << "Check whether extrapolated track point on PMT plane verifies its "
1062 "equation (value should be 0.):"
1063 << endl;
1064 cout << "* using mirror point M, checkCalc = " << checkCalc1 << endl;
1065 checkCalc2 =
1066 ptPR2.at(0) * normalPMT.at(0) + ptPR2.at(1) * normalPMT.at(1) + ptPR2.at(2) * normalPMT.at(2) + constantePMT;
1067 //cout << "* using reflected point R2, checkCalc = " << checkCalc2 << endl;
1068}
1069
1070void CbmRichCorrection::FillHistProjection(TVector3 outPosIdeal, TVector3 outPosUnCorr, TVector3 outPos,
1071 Int_t NofGTracks, vector<Double_t> normalPMT, Double_t constantePMT)
1072{
1073 CbmMCTrack* track2 = NULL;
1074 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0, distToExtrapTrackHitInPlane = 0,
1075 distToExtrapTrackHitInPlaneUnCorr = 0, distToExtrapTrackHitInPlaneIdeal = 0;
1076
1077 for (Int_t iGlobalTrack = 0; iGlobalTrack < NofGTracks; iGlobalTrack++) {
1078 //cout << "Nb of global tracks = " << NofGTracks << " and iGlobalTrack = " << iGlobalTrack << endl;
1079 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack);
1080 if (NULL == gTrack) continue;
1081 Int_t richInd = gTrack->GetRichRingIndex();
1082 //cout << "Rich index = " << richInd << endl;
1083 if (richInd < 0) {
1084 continue;
1085 }
1086 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richInd));
1087 if (NULL == ring) {
1088 continue;
1089 }
1090 Int_t ringTrackID = ring->GetTrackID();
1091 track2 = (CbmMCTrack*) fMCTracks->At(ringTrackID);
1092 Int_t ringMotherID = track2->GetMotherId();
1093 //cout << "Ring mother ID = " << ringMotherID << ", ring track ID = " << ringTrackID << " and track2 pdg = " << track2->GetPdgCode() << endl;
1094 CbmRichRingLight ringL;
1095 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
1096 fCopFit->DoFit(&ringL);
1097 //fTauFit->DoFit(&ringL);
1098 ringCenter[0] = ringL.GetCenterX();
1099 ringCenter[1] = ringL.GetCenterY();
1100 ringCenter[2] =
1101 -1 * ((normalPMT.at(0) * ringCenter[0] + normalPMT.at(1) * ringCenter[1] + constantePMT) / normalPMT.at(2));
1102
1103 vector<Double_t> r(3),
1104 p(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1105 r.at(0) = TMath::Abs(ringCenter[0]), r.at(1) = TMath::Abs(ringCenter[1]), r.at(2) = TMath::Abs(ringCenter[2]);
1106 p.at(0) = TMath::Abs(outPos.x()), p.at(1) = TMath::Abs(outPos.y()), p.at(2) = TMath::Abs(outPos.z());
1107 cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}"
1108 << endl;
1109 cout << "Difference in X = " << TMath::Abs(r.at(0) - p.at(0)) << "; \t"
1110 << "Difference in Y = " << TMath::Abs(r.at(1) - p.at(1)) << "; \t"
1111 << "Difference in Z = " << TMath::Abs(r.at(2) - p.at(2)) << endl;
1112 fHM->H1("fhDifferenceX")->Fill(TMath::Abs(r.at(0) - p.at(0)));
1113 fHM->H1("fhDifferenceY")->Fill(TMath::Abs(r.at(1) - p.at(1)));
1114 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2)
1115 + TMath::Power(r.at(2) - p.at(2), 2));
1116 distToExtrapTrackHitInPlane = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2));
1117 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
1118 fHM->H1("fhDistanceCorrected")->Fill(distToExtrapTrackHitInPlane);
1119 //cout << "Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl;
1120 cout << "Distance between fitted ring center and extrapolated track hit in "
1121 "plane = "
1122 << distToExtrapTrackHitInPlane << endl;
1123
1124 vector<Double_t> pUnCorr(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1125 pUnCorr.at(0) = TMath::Abs(outPosUnCorr.x()), pUnCorr.at(1) = TMath::Abs(outPosUnCorr.y()),
1126 pUnCorr.at(2) = TMath::Abs(outPosUnCorr.z());
1127 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1128 cout << "Difference in X = " << TMath::Abs(r.at(0) - pUnCorr.at(0)) << "; \t"
1129 << "Difference in Y = " << TMath::Abs(r.at(1) - pUnCorr.at(1)) << "; \t"
1130 << "Difference in Z = " << TMath::Abs(r.at(2) - pUnCorr.at(2)) << endl;
1131 fHM->H1("fhDifferenceXUncorrected")->Fill(TMath::Abs(r.at(0) - pUnCorr.at(0)));
1132 fHM->H1("fhDifferenceYUncorrected")->Fill(TMath::Abs(r.at(1) - pUnCorr.at(1)));
1133 distToExtrapTrackHitInPlaneUnCorr =
1134 TMath::Sqrt(TMath::Power(r.at(0) - pUnCorr.at(0), 2) + TMath::Power(r.at(1) - pUnCorr.at(1), 2));
1135 fHM->H1("fhDistanceUncorrected")->Fill(distToExtrapTrackHitInPlaneUnCorr);
1136 cout << "Distance between fitted ring center and extrapolated track hit in "
1137 "plane = "
1138 << distToExtrapTrackHitInPlaneUnCorr << endl;
1139
1140 vector<Double_t> pIdeal(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1141 pIdeal.at(0) = TMath::Abs(outPosIdeal.x()), pIdeal.at(1) = TMath::Abs(outPosIdeal.y()),
1142 pIdeal.at(2) = TMath::Abs(outPosIdeal.z());
1143 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1144 cout << "Difference in X = " << TMath::Abs(r.at(0) - pIdeal.at(0)) << "; \t"
1145 << "Difference in Y = " << TMath::Abs(r.at(1) - pIdeal.at(1)) << "; \t"
1146 << "Difference in Z = " << TMath::Abs(r.at(2) - pIdeal.at(2)) << endl;
1147 fHM->H1("fhDifferenceXIdeal")->Fill(TMath::Abs(r.at(0) - pIdeal.at(0)));
1148 fHM->H1("fhDifferenceYIdeal")->Fill(TMath::Abs(r.at(1) - pIdeal.at(1)));
1149 distToExtrapTrackHitInPlaneIdeal =
1150 TMath::Sqrt(TMath::Power(r.at(0) - pIdeal.at(0), 2) + TMath::Power(r.at(1) - pIdeal.at(1), 2));
1151 fHM->H1("fhDistanceIdeal")->Fill(distToExtrapTrackHitInPlaneIdeal);
1152 cout << "Distance between fitted ring center and extrapolated track hit in "
1153 "plane = "
1154 << distToExtrapTrackHitInPlaneIdeal << endl
1155 << endl;
1156 //}
1157 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
1158 }
1159 cout << "End of loop on global tracks;" << endl;
1160}
1161
1162void CbmRichCorrection::DrawHistProjection()
1163{
1164 char leg[128];
1165 int colorInd = 1;
1166
1167 TCanvas* can3 = new TCanvas(fRunTitle + "_Distance_Histos_" + fAxisRotTitle,
1168 fRunTitle + "_Distance_Histos_" + fAxisRotTitle, 1500, 400);
1169 can3->SetGrid(1, 1);
1170 can3->Divide(2, 1);
1171 can3->cd(1)->SetGrid(1, 1);
1172 can3->cd(2)->SetGrid(1, 1);
1173 can3->cd(1);
1174 TH1D* Clone1 = (TH1D*) fHM->H1("fhDifferenceXIdeal")->Clone();
1175 Clone1->GetXaxis()->SetTitleSize(0.04);
1176 Clone1->GetYaxis()->SetTitleSize(0.04);
1177 Clone1->GetXaxis()->SetLabelSize(0.03);
1178 Clone1->GetYaxis()->SetLabelSize(0.03);
1179 Clone1->GetXaxis()->CenterTitle();
1180 Clone1->GetYaxis()->CenterTitle();
1181 Clone1->SetTitle("Difference in X");
1182 Clone1->SetLineColor(kBlue);
1183 Clone1->SetLineWidth(2);
1184 Clone1->Rebin(2);
1185 Clone1->Draw();
1186 TH1D* Clone2 = (TH1D*) fHM->H1("fhDifferenceX")->Clone();
1187 Clone2->SetTitleSize(0.04);
1188 Clone2->SetLabelSize(0.03);
1189 Clone2->SetLineColor(kGreen);
1190 Clone2->SetLineWidth(2);
1191 Clone2->Rebin(2);
1192 Clone2->Draw("same");
1193 TH1D* Clone3 = (TH1D*) fHM->H1("fhDifferenceXUncorrected")->Clone();
1194 Clone3->SetTitleSize(0.04);
1195 Clone3->SetLabelSize(0.03);
1196 Clone3->SetLineColor(kRed);
1197 Clone3->SetLineWidth(2);
1198 Clone3->Rebin(2);
1199 Clone3->Draw("same");
1200 gStyle->Clear();
1201
1202 TLegend* LEG = new TLegend(0.3, 0.78, 0.5, 0.88); // Set legend position
1203 LEG->SetBorderSize(1);
1204 LEG->SetFillColor(0);
1205 LEG->SetMargin(0.2);
1206 LEG->SetTextSize(0.02);
1207 sprintf(leg, "X diff uncorr");
1208 LEG->AddEntry(Clone3, leg, "l");
1209 sprintf(leg, "X diff corr");
1210 LEG->AddEntry(Clone2, leg, "l");
1211 sprintf(leg, "X diff ideal");
1212 LEG->AddEntry(Clone1, leg, "l");
1213 LEG->Draw();
1214
1215 can3->cd(2);
1216 can3->SetGrid(1, 1);
1217 TH1D* Clone4 = (TH1D*) fHM->H1("fhDifferenceYIdeal")->Clone();
1218 Clone4->GetXaxis()->SetTitleSize(0.04);
1219 Clone4->GetYaxis()->SetTitleSize(0.04);
1220 Clone4->GetXaxis()->SetLabelSize(0.03);
1221 Clone4->GetYaxis()->SetLabelSize(0.03);
1222 Clone4->GetXaxis()->CenterTitle();
1223 Clone4->GetYaxis()->CenterTitle();
1224 Clone4->SetTitle("Difference in Y");
1225 Clone4->SetLineColor(kBlue);
1226 Clone4->SetLineWidth(2);
1227 Clone4->Rebin(2);
1228 Clone4->Draw();
1229 TH1D* Clone5 = (TH1D*) fHM->H1("fhDifferenceY")->Clone();
1230 Clone5->SetTitleSize(0.04);
1231 Clone5->SetLabelSize(0.03);
1232 Clone5->SetLineColor(kGreen);
1233 Clone5->SetLineWidth(2);
1234 Clone5->Rebin(2);
1235 Clone5->Draw("same");
1236 TH1D* Clone6 = (TH1D*) fHM->H1("fhDifferenceYUncorrected")->Clone();
1237 Clone6->SetTitleSize(0.04);
1238 Clone6->SetLabelSize(0.03);
1239 Clone6->SetLineColor(kRed);
1240 Clone6->SetLineWidth(2);
1241 Clone6->Rebin(2);
1242 Clone6->Draw("same");
1243
1244 TLegend* LEG1 = new TLegend(0.3, 0.78, 0.5, 0.88); // Set legend position
1245 LEG1->SetBorderSize(1);
1246 LEG1->SetFillColor(0);
1247 LEG1->SetMargin(0.2);
1248 LEG1->SetTextSize(0.02);
1249 sprintf(leg, "Y diff uncorr");
1250 LEG1->AddEntry(Clone6, leg, "l");
1251 sprintf(leg, "Y diff corr");
1252 LEG1->AddEntry(Clone5, leg, "l");
1253 sprintf(leg, "Y diff ideal");
1254 LEG1->AddEntry(Clone4, leg, "l");
1255 LEG1->Draw();
1256
1257 /*can3->cd(3);
1258 //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);
1259 TH1D* CloneDist1 = (TH1D*)fHM->H1("fhDistanceIdeal")->Clone();
1260 CloneDist1->GetXaxis()->SetTitleSize(0.04);
1261 CloneDist1->GetYaxis()->SetTitleSize(0.04);
1262 CloneDist1->GetXaxis()->SetLabelSize(0.03);
1263 CloneDist1->GetYaxis()->SetLabelSize(0.03);
1264 CloneDist1->GetXaxis()->CenterTitle();
1265 CloneDist1->GetYaxis()->CenterTitle();
1266 CloneDist1->SetTitle("Distance between extrapolated track center and fitted ring center");
1267 CloneDist1->SetLineColor(kBlue);
1268 CloneDist1->SetLineWidth(2);
1269 CloneDist1->Rebin(2);
1270 CloneDist1->Draw();
1271 TH1D* CloneDist2 = (TH1D*)fHM->H1("fhDistanceUncorrected")->Clone();
1272 CloneDist2->SetTitleSize(0.04);
1273 CloneDist2->SetTitleSize(0.04);
1274 CloneDist2->SetLineColor(kRed);
1275 CloneDist2->SetLineWidth(2);
1276 CloneDist2->Rebin(2);
1277 CloneDist2->Draw("same");
1278 TH1D* CloneDist3 = (TH1D*)fHM->H1("fhDistanceCorrected")->Clone();
1279 CloneDist3->SetTitleSize(0.04);
1280 CloneDist3->SetTitleSize(0.04);
1281 CloneDist3->SetLineColor(kGreen);
1282 CloneDist3->SetLineWidth(2);
1283 CloneDist3->Rebin(2);
1284 CloneDist3->Draw("same");
1285
1286 TLegend* LEG2 = new TLegend(0.3,0.78,0.5,0.88); // Set legend position
1287 LEG2->SetBorderSize(1);
1288 LEG2->SetFillColor(0);
1289 LEG2->SetMargin(0.2);
1290 LEG2->SetTextSize(0.02);
1291 sprintf(leg, "Distance uncorr");
1292 LEG2->AddEntry(CloneDist2, leg, "l");
1293 sprintf(leg, "Distance corr");
1294 LEG2->AddEntry(CloneDist3, leg, "l");
1295 sprintf(leg, "Distance ideal");
1296 LEG2->AddEntry(CloneDist1, leg, "l");
1297 LEG2->Draw();*/
1298 gStyle->SetOptStat(000000);
1299
1300 Cbm::SaveCanvasAsImage(can3, string(fOutputDir.Data()), "png");
1301
1302 /*TCanvas* can4 = new TCanvas(fRunTitle + "_Difference_Fits_" + fAxisRotTitle, fRunTitle + "_Difference_Fits_" + fAxisRotTitle, 800, 800);
1303 int colorInd3 = 1;
1304 can4->Divide(2,2);
1305 can4->cd(1);
1306 //fHM->H1("fHistoDiffX")->Add(fHM->H1("fhDifferenceX"), fHM->H1("fhDifferenceXIdeal"), -1, 1);
1307 Clone2->GetXaxis()->SetTitleSize(0.04);
1308 Clone2->GetYaxis()->SetTitleSize(0.04);
1309 Clone2->GetXaxis()->SetLabelSize(0.03);
1310 Clone2->GetYaxis()->SetLabelSize(0.03);
1311 Clone2->GetXaxis()->CenterTitle();
1312 Clone2->GetYaxis()->CenterTitle();
1313 Clone2->SetTitle("Corrected difference in X");
1314 //Clone2->SetAxisRange(0., 0.6);
1315 Clone2->SetLineColor(kGreen);
1316 //Clone2->SetLineColor(colorInd3);
1317 //colorInd3++;
1318 Clone2->Draw();
1319 Clone2->Fit("gaus", "0", "", 0., 1.5);
1320 gStyle->SetOptFit(1111);
1321 TF1 *fit1 = Clone2->GetFunction("gaus");
1322 can4->cd(2);
1323 //Clone3->SetLineColor(colorInd3);
1324 //colorInd3++;
1325 //Clone1->SetAxisRange(0., 0.6);
1326 Clone1->SetTitle("Difference in X ideal");
1327 Clone1->Draw();
1328 Clone1->Fit("gaus", "0", "", 0., 1.5);
1329 gStyle->SetOptFit(1111);
1330 TF1 *fit2 = Clone1->GetFunction("gaus");
1331
1332 can4->cd(3);
1333 Clone5->GetXaxis()->SetTitleSize(0.04);
1334 Clone5->GetYaxis()->SetTitleSize(0.04);
1335 Clone5->GetXaxis()->SetLabelSize(0.03);
1336 Clone5->GetYaxis()->SetLabelSize(0.03);
1337 Clone5->GetXaxis()->CenterTitle();
1338 Clone5->GetYaxis()->CenterTitle();
1339 Clone5->SetTitle("Corrected difference in Y");
1340 //Clone5->SetAxisRange(0., 0.6);
1341 //Clone5->SetLineColor(colorInd3);
1342 //colorInd3++;
1343 Clone5->Draw();
1344 Clone5->Fit("gaus", "0", "", 0., 1.5);
1345 gStyle->SetOptFit(1111);
1346 TF1 *fit3 = Clone5->GetFunction("gaus");
1347 can4->cd(4);
1348 //Clone4->SetAxisRange(0., 0.6);
1349 //Clone6->SetLineColor(colorInd3);
1350 //colorInd3++;
1351 Clone4->SetTitle("Difference in Y ideal");
1352 Clone4->Draw();
1353 Clone4->Fit("gaus", "0", "", 0., 1.5);
1354 gStyle->SetOptFit(1111);
1355 TF1 *fit4 = Clone4->GetFunction("gaus");
1356
1357 //Cbm::SaveCanvasAsImage(can4, string(fOutputDir.Data()), "png");
1358
1359 Double_t meanX_corr = fit1->GetParameter(1);
1360 Double_t meanX_ideal = fit2->GetParameter(1);
1361 cout << endl;
1362 cout << "Mean X corrected = " << meanX_corr << " and mean X ideal = " << meanX_ideal << endl;
1363 cout << "Fitted parameters of differenceX histo = " << fit1->GetParameter(0) << ", " << fit1->GetParameter(1) << " and " << fit1->GetParameter(2) << endl;
1364 cout << "Fitted parameters of differenceXIdeal histo = " << fit2->GetParameter(0) << ", " << fit2->GetParameter(1) << " and " << fit2->GetParameter(2) << endl;
1365 Double_t meanY_corr = fit3->GetParameter(1);
1366 Double_t meanY_ideal = fit4->GetParameter(1);
1367 cout << "Fitted parameters of differenceY histo = " << fit3->GetParameter(0) << ", " << fit3->GetParameter(1) << " and " << fit3->GetParameter(2) << endl;
1368 cout << "Fitted parameters of differenceYIdeal histo = " << fit4->GetParameter(0) << ", " << fit4->GetParameter(1) << " and " << fit4->GetParameter(2) << endl;
1369 cout << "Mean Y corrected = " << meanY_corr << " and mean Y ideal = " << meanY_ideal << endl;*/
1370}
1371
1372void CbmRichCorrection::DrawHistFromFile(TString fileName)
1373{
1375 TFile* oldFile = gFile;
1376 TDirectory* oldDir = gDirectory;
1377
1378 fHM = new CbmHistManager();
1379 TFile* file = new TFile(fileName, "READ");
1380 fHM->ReadFromFile(file);
1381 DrawHistProjection();
1382
1384 gFile = oldFile;
1385 gDirectory = oldDir;
1386}
1387
1388void CbmRichCorrection::Finish()
1389{
1390 // ---------------------------------------------------------------------------------------------------------------------------------------- //
1391 // -------------------------------------------------- Mapping for mirror - PMT relations -------------------------------------------------- //
1392 // ---------------------------------------------------------------------------------------------------------------------------------------- //
1393
1394 if (fDrawProjection) {
1395 DrawHistProjection();
1396 fHM->H1("fhDifferenceXUncorrected")->Write();
1397 fHM->H1("fhDifferenceYUncorrected")->Write();
1398 fHM->H1("fhDistanceUncorrected")->Write();
1399 fHM->H1("fhDifferenceX")->Write();
1400 fHM->H1("fhDifferenceY")->Write();
1401 fHM->H1("fhDistanceCorrected")->Write();
1402 fHM->H1("fhDifferenceXIdeal")->Write();
1403 fHM->H1("fhDifferenceYIdeal")->Write();
1404 fHM->H1("fhDistanceIdeal")->Write();
1405 }
1406
1407 //cout << setprecision(6) << endl;
1408}
1409ClassImp(CbmRichCorrection)
1410
1411 /*
1412vector<Double_t> CbmRichCorrection::RotateSphereCenter(vector<Double_t> ptM, vector<Double_t> ptC, TGeoNavigator* navi)
1413{
1414 vector<Double_t> ptCNew(3);
1415 Double_t cosPhi=0., sinPhi=0., cosTheta=0., sinTheta=0., diff[3], mat[3][3], invMat[3][3], corrMat[3][3], buff1[3][3], buff2[3][3];
1416
1417 InvertMatrix(mat, invMat, navi);
1418 /*for (Int_t i=0; i<3; i++) {
1419 for (Int_t j=0; j<3; j++) {
1420 cout << invMat[i][j] << "\t";
1421 }
1422 cout << endl;
1423 }*/
1424 /*
1425 // Reading misalignment information from correction_param.txt text file.
1426 vector<Float_t> outputFit(4);
1427 ifstream corr_file;
1428 corr_file.open("correction_param.txt");
1429 if (corr_file.is_open())
1430 {
1431 for (Int_t i=0; i<4; i++) {corr_file >> outputFit.at(i);}
1432 corr_file.close();
1433 }
1434 else {
1435 cout << "Error in CbmRichCorrection: unable to open parameter file!" << endl << endl;
1436 sleep(5);
1437 }
1438 cout << "Misalignment parameters read from file = [" << outputFit.at(0) << " ; " << outputFit.at(1) << " ; " << outputFit.at(2) << " ; " << outputFit.at(3) << "]" << endl;
1439
1440 // Calculating the new coordinates of sphere center C, according to the correction parameters.
1441 for (Int_t i=0; i<3; i++) {
1442 for (Int_t j=0; j<3; j++) {
1443 corrMat[i][j] = 0.;
1444 buff1[i][j] = 0.;
1445 buff2[i][j] = 0.;
1446 }
1447 }
1448
1449 // Correction Matrix = Id(3).
1450 //corrMat[0][0] = corrMat[1][1] = corrMat[2][2] = 1;
1451
1452 cosPhi = TMath::Cos(outputFit.at(0));
1453 sinPhi = TMath::Sin(outputFit.at(0)); // BEWARE !!! the horizontal rotation axis of the tile is opposite to the global horizontal rotation axis !!!
1454 cosTheta = TMath::Cos(outputFit.at(1));
1455 sinTheta = -1*TMath::Sin(outputFit.at(1));
1456
1457 Double_t RotX[3][3], RotY[3][3], RotZ[3][3];
1458 // Define rotation around X axis, with y towards z.
1459 // y towards z is the rotation direction defined in the tile frame, which is the same definition as for the global
1460 // frame. But as the X axis direction is opposite in the two frames, to correct for misalignment, the rotation
1461 // direction remains unchanged in the global frame.
1462 RotX[0][0] = 1;
1463 RotX[0][1] = 0;
1464 RotX[0][2] = 0;
1465 RotX[1][0] = 0;
1466 RotX[1][1] = cosPhi;
1467 RotX[1][2] = -1*sinPhi;
1468 RotX[2][0] = 0;
1469 RotX[2][1] = sinPhi;
1470 RotX[2][2] = cosPhi;
1471 // Define rotation around Y axis, with z towards x.
1472 // x towards z is the rotation direction defined in the tile frame, which is the same definition as for the global
1473 // frame. And as the tile and global frames are identical, the rotation direction for the correction is the opposite
1474 // in the global frame - in order to correct for the mibuff3[i][j] = 0.;salignment.
1475 RotY[0][0] = cosTheta;
1476 RotY[0][1] = 0;
1477 RotY[0][2] = sinTheta;
1478 RotY[1][0] = 0;
1479 RotY[1][1] = 1;
1480 RotY[1][2] = 0;
1481 RotY[2][0] = -1*sinTheta;
1482 RotY[2][1] = 0;
1483 RotY[2][2] = cosTheta;
1484 // Define rotation around Z axis.
1485 RotZ[0][0] = cosTheta;
1486 RotZ[0][1] = -1*sinTheta;
1487 RotZ[0][2] = 0;
1488 RotZ[1][0] = sinTheta;
1489 RotZ[1][1] = cosTheta;
1490 RotZ[1][2] = 0;
1491 RotZ[2][0] = 0;
1492 RotZ[2][1] = 0;
1493 RotZ[2][2] = 1;
1494 cout << "RotY[0][0] = " << RotY[0][0] << " and RotY[0][2] = " << RotY[0][2] << endl;
1495
1496 for (Int_t i=0; i<3; i++) {
1497 for (Int_t j=0; j<3; j++) {
1498 //corrMat[i][j] = RotZ[i][j];
1499 for (Int_t k=0; k<3; k++) {
1500 corrMat[i][j] = RotY[i][k]*RotX[k][j] + corrMat[i][j];
1501 }
1502 }
1503 }
1504
1505/* corrMat[0][0] = cosTheta;
1506 corrMat[0][1] = -1*sinTheta*sinPhi;
1507 corrMat[0][2] = -1*sinTheta*cosPhi;
1508 corrMat[1][0] = 0;
1509 corrMat[1][1] = cosPhi;
1510 corrMat[1][2] = -1*sinPhi;
1511 corrMat[2][0] = sinTheta;
1512 corrMat[2][1] = cosPhi;
1513 corrMat[2][2] = cosTheta*cosPhi;
1514 corrMat[0][0] = 1;
1515 corrMat[0][1] = 0;
1516 corrMat[0][2] = 0;
1517 corrMat[1][0] = 0;
1518 corrMat[1][1] = cosPhi;
1519 corrMat[1][2] = -1*sinPhi;
1520 corrMat[2][0] = 0;
1521 corrMat[2][1] = sinPhi;
1522 corrMat[2][2] = cosPhi;
1523 corrMat[0][0] = cosTheta;
1524 corrMat[0][1] = 0;
1525 corrMat[0][2] = -1*sinTheta;
1526 corrMat[1][0] = -1*sinTheta*sinPhi;
1527 corrMat[1][1] = cosPhi;
1528 corrMat[1][2] = -1*sinPhi*cosTheta;
1529 corrMat[2][0] = cosPhi*sinTheta;
1530 corrMat[2][1] = sinPhi;
1531 corrMat[2][2] = cosTheta*cosPhi;
1532
1533 for (Int_t i=0; i<3; i++) {
1534 for (Int_t j=0; j<3; j++) {
1535 for (Int_t k=0; k<3; k++) {
1536 buff1[i][j] = corrMat[i][k]*invMat[k][j] + buff1[i][j];
1537 }
1538 }
1539 }
1540 for (Int_t i=0; i<3; i++) {
1541 for (Int_t j=0; j<3; j++) {
1542 for (Int_t k=0; k<3; k++) {
1543 buff2[i][j] = mat[i][k]*buff1[k][j] + buff2[i][j];
1544 }
1545 }
1546 }
1547
1548 diff[0] = ptC.at(0) - ptM.at(0);
1549 diff[1] = ptC.at(1) - ptM.at(1);
1550 diff[2] = ptC.at(2) - ptM.at(2);
1551 cout << "diff = {" << diff[0] << ", " << diff[1] << ", " << diff[2] << "}" << endl;
1552 Double_t phi2 =TMath::ATan2(diff[1], -1*diff[2]), theta2 =TMath::ATan2(diff[0], -1*diff[2]);
1553 cout << "Calculated phi (= arctan(y/-z)): " << TMath::RadToDeg()*phi2 << " and calculated theta (= arctan(x/-z)): " << TMath::RadToDeg()*theta2 << endl;
1554 cout << "cosPhi = " << TMath::Cos(phi2) << ", sinPhi = " << TMath::Sin(phi2) << " and cosTheta = " << TMath::Cos(theta2) << ", sinTheta = " << TMath::Sin(theta2) << endl;
1555 Double_t RotX2[3][3], RotY2[3][3], corrMat2[3][3];
1556 RotX2[0][0] = 1;
1557 RotX2[0][1] = 0;
1558 RotX2[0][2] = 0;
1559 RotX2[1][0] = 0;
1560 RotX2[1][1] = TMath::Cos(phi2);
1561 RotX2[1][2] = TMath::Sin(phi2);
1562 RotX2[2][0] = 0;
1563 RotX2[2][1] = -1*TMath::Sin(phi2);
1564 RotX2[2][2] = TMath::Cos(phi2);
1565 RotY2[0][0] = TMath::Cos(theta2);
1566 RotY2[0][1] = 0;
1567 RotY2[0][2] = TMath::Sin(theta2);
1568 RotY2[1][0] = 0;
1569 RotY2[1][1] = 1;
1570 RotY2[1][2] = 0;
1571 RotY2[2][0] = -1*TMath::Sin(theta2);
1572 RotY2[2][1] = 0;
1573 RotY2[2][2] = TMath::Cos(theta2);
1574
1575 for (Int_t i=0; i<3; i++) {
1576 for (Int_t j=0; j<3; j++) {
1577 for (Int_t k=0; k<3; k++) {
1578 corrMat2[i][j] = RotX2[i][k]*RotY2[k][j] + corrMat2[i][j];
1579 }
1580 }
1581 }
1582
1583 vector<Double_t> ptCNew2(3), ptCNew3(3);
1584 for (Int_t i=0; i<3; i++) {
1585 for (Int_t j=0; j<3; j++) {
1586 //ptCNew.at(i) = buff2[i][j]*diff[j] + ptCNew.at(i);
1587 //ptCNew.at(i) = invMat[i][j]*diff[j] + ptCNew.at(i);
1588 ptCNew.at(i) = corrMat2[i][j]*diff[j] + ptCNew.at(i);
1589 ptCNew2.at(i) = RotY[i][j]*diff[j] + ptCNew2.at(i);
1590 }
1591 }
1592
1593 phi2 = TMath::ATan2(ptCNew2.at(1), -ptCNew2.at(2));
1594 cout << "New phi2 = " << TMath::RadToDeg()*phi2 << endl;
1595 RotX2[0][0] = 1;
1596 RotX2[0][1] = 0;
1597 RotX2[0][2] = 0;
1598 RotX2[1][0] = 0;
1599 RotX2[1][1] = TMath::Cos(phi2);
1600 RotX2[1][2] = TMath::Sin(phi2);
1601 RotX2[2][0] = 0;
1602 RotX2[2][1] = -1*TMath::Sin(phi2);
1603 RotX2[2][2] = TMath::Cos(phi2);
1604 for (Int_t i=0; i<3; i++) {
1605 for (Int_t j=0; j<3; j++) {
1606 ptCNew3.at(i) = RotX[i][j]*diff[j] + ptCNew3.at(i);
1607 }
1608 }
1609
1610 //ptCNew.at(0) += ptM.at(0);
1611 //ptCNew.at(1) += ptM.at(1);
1612 //ptCNew.at(2) += ptM.at(2);
1613
1614 /*ptCNew.at(0) = ptM.at(0) + diffX*cosTheta - diffY*sinTheta*sinPhi - diffZ*sinTheta*cosPhi;
1615 ptCNew.at(1) = ptM.at(1) + diffY*cosPhi - diffZ*sinPhi;
1616 ptCNew.at(2) = ptM.at(2) + diffX*sinTheta + diffY*cosPhi + diffZ*cosTheta*cosPhi;*/
1617 /* cout << "Sphere center coordinates of the rotated mirror tile, w/ calculation, = {" << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
1618 cout << "Sphere center coordinates of the rotated mirror tile (after rotation around Y), w/ calculation, = {" << ptCNew2.at(0) << ", " << ptCNew2.at(1) << ", " << ptCNew2.at(2) << "}" << endl;
1619 cout << "Sphere center coordinates of the rotated mirror tile (around unrotated axes), w/ calculation, = {" << ptCNew3.at(0) << ", " << ptCNew3.at(1) << ", " << ptCNew3.at(2) << "}" << endl << endl;
1620
1621 return ptCNew;
1622}
1623*/
ClassImp
ClassImp(CbmConverterManager)
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichConverter.h
Convert internal data classes to cbmroot common data classes.
CbmRichHitProducer.h
Class for producing RICH hits directly from MCPoints.
CbmRichRingFitterCOP.h
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
CbmRichRingFitterEllipseTau.h
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmHistManager
Histogram manager.
Definition CbmHistManager.h:48
CbmHistManager::ReadFromFile
void ReadFromFile(TFile *file)
Read histograms from file.
Definition core/base/CbmHistManager.cxx:131
CbmHistManager::Create1
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...
Definition CbmHistManager.h:88
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition CbmHistManager.h:207
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmRichConverter::Init
static void Init()
Initialize array of RICH hits.
Definition CbmRichConverter.h:98
CbmRichConverter::CopyHitsToRingLight
static void CopyHitsToRingLight(const CbmRichRing *ring1, CbmRichRingLight *ring2)
Copy hits from CbmRichRing to CbmRichRingLight.
Definition CbmRichConverter.h:45
CbmRichCorrection::fGlobalTracks
TClonesArray * fGlobalTracks
Definition CbmRichCorrection.h:176
CbmRichCorrection::fMCTracks
TClonesArray * fMCTracks
Definition CbmRichCorrection.h:172
CbmRichCorrection::InvertMatrix
void InvertMatrix(Double_t mat[3][3], Double_t invMat[3][3], TGeoNavigator *navi)
Definition CbmRichCorrection.cxx:809
CbmRichCorrection::RotateSphereCenter
vector< Double_t > RotateSphereCenter(vector< Double_t > ptM, vector< Double_t > ptC, TGeoNavigator *navi)
Definition CbmRichCorrection.cxx:601
CbmRichCorrection::Finish
virtual void Finish()
Inherited from FairTask.
Definition CbmRichCorrection.cxx:1388
CbmRichCorrection::FillHistProjection
void FillHistProjection(TVector3 outPosIdeal, TVector3 outPosUnCorr, TVector3 outPos, Int_t NofGlobalTracks, vector< Double_t > normalPMT, Double_t constantePMT)
Definition CbmRichCorrection.cxx:1070
CbmRichCorrection::fRichRingMatches
TClonesArray * fRichRingMatches
Definition CbmRichCorrection.h:173
CbmRichCorrection::fRichPoints
TClonesArray * fRichPoints
Definition CbmRichCorrection.h:175
CbmRichCorrection::fRichMirrorPoints
TClonesArray * fRichMirrorPoints
Definition CbmRichCorrection.h:169
CbmRichCorrection::fRichHits
TClonesArray * fRichHits
Definition CbmRichCorrection.h:167
CbmRichCorrection::fRichProjections
TClonesArray * fRichProjections
Definition CbmRichCorrection.h:170
CbmRichCorrection::fRichRefPlanePoints
TClonesArray * fRichRefPlanePoints
Definition CbmRichCorrection.h:174
CbmRichCorrection::ComputeP
void ComputeP(vector< Double_t > &ptPMirr, vector< Double_t > &ptPR2, vector< Double_t > normalPMT, vector< Double_t > ptM, vector< Double_t > ptR2Mirr, Double_t normalCste)
Definition CbmRichCorrection.cxx:1027
CbmRichCorrection::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition CbmRichCorrection.cxx:199
CbmRichCorrection::fRichRings
TClonesArray * fRichRings
Definition CbmRichCorrection.h:168
CbmRichCorrection::CalculateMirrorIntersection
void CalculateMirrorIntersection(vector< Double_t > ptM, vector< Double_t > ptCUnCorr, vector< Double_t > &ptMNew)
Definition CbmRichCorrection.cxx:895
CbmRichCorrection::GetPmtNormal
void GetPmtNormal(Int_t NofPMTPoints, vector< Double_t > &normalPMT, Double_t &normalCste)
Definition CbmRichCorrection.cxx:410
CbmRichCorrection::GetMeanSphereCenter
void GetMeanSphereCenter(TGeoNavigator *navi, vector< Double_t > &ptC)
Definition CbmRichCorrection.cxx:500
CbmRichCorrection::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition CbmRichCorrection.cxx:92
CbmRichCorrection::GetMirrorIntersection
void GetMirrorIntersection(vector< Double_t > &ptM, vector< Double_t > ptR1, vector< Double_t > momR1, vector< Double_t > ptC, Double_t sphereRadius)
Definition CbmRichCorrection.cxx:558
CbmRichCorrection::ComputeR2
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 s)
Definition CbmRichCorrection.cxx:911
CbmRichCorrection::fHM
CbmHistManager * fHM
Definition CbmRichCorrection.h:177
CbmRichCorrection::fCopFit
CbmRichRingFitterCOP * fCopFit
Definition CbmRichCorrection.h:192
CbmRichCorrection::DrawHistFromFile
void DrawHistFromFile(TString fileName)
Definition CbmRichCorrection.cxx:1372
CbmRichCorrection::fRichMCPoints
TClonesArray * fRichMCPoints
Definition CbmRichCorrection.h:171
CbmRichCorrection::fTauFit
CbmRichRingFitterEllipseTau * fTauFit
Definition CbmRichCorrection.h:193
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition CbmRichGeoManager.cxx:399
CbmRichGeoManager::fGP
CbmRichRecGeoPar * fGP
Definition CbmRichGeoManager.h:28
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition CbmRichGeoManager.h:34
CbmRichRecGeoPar
PMT parameters for the RICH geometry.
Definition CbmRichRecGeoPar.h:94
CbmRichRecGeoPar::fPmt
CbmRichRecGeoParPmt fPmt
Definition CbmRichRecGeoPar.h:218
CbmRichRingFitterCOP
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
Definition CbmRichRingFitterCOP.h:32
CbmRichRingFitterCOP::DoFit
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition CbmRichRingFitterCOP.cxx:24
CbmRichRingFitterEllipseTau
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
Definition CbmRichRingFitterEllipseTau.h:39
CbmRichRingLight::GetCenterX
float GetCenterX() const
Definition CbmRichRingLight.h:165
CbmRichRingLight::GetCenterY
float GetCenterY() const
Definition CbmRichRingLight.h:166
Cbm::SaveCanvasAsImage
void SaveCanvasAsImage(TCanvas *c, const string &dir, const string &option)
Definition CbmUtils.cxx:36