CbmRoot
Loading...
Searching...
No Matches
CbmRichCorrectionVector.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 "CbmRichCorrectionVector.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
54CbmRichCorrectionVector::CbmRichCorrectionVector()
55 : FairTask()
56 , fRichHits(NULL)
57 , fRichRings(NULL)
58 , fRichProjections(NULL)
59 , fRichMirrorPoints(NULL)
60 , fRichMCPoints(NULL)
61 , fMCTracks(NULL)
62 , fRichRingMatches(NULL)
63 , fRichRefPlanePoints(NULL)
64 , fRichPoints(NULL)
65 , fGlobalTracks(NULL)
66 , fHM(NULL)
67 , fHM2(NULL)
68 ,
69 //fGP(),
70 fEventNum(0)
71 , fOutputDir("")
72 , fRunTitle("")
73 , fAxisRotTitle("")
74 , fDrawAlignment(kTRUE)
75 , fDrawMapping(kFALSE)
76 , fDrawProjection(kTRUE)
77 , fIsMeanCenter(kFALSE)
78 , fIsReconstruction(kFALSE)
79 , fCopFit(NULL)
80 , fTauFit(NULL)
81 , fPathsMap()
82 , fPathsMapEllipse()
83 , fPhi()
84{
85 fMirrCounter = 0.;
86 for (int i = 0; i < 3; i++) {
87 fArray[i] = 0.;
88 }
89}
90
91CbmRichCorrectionVector::~CbmRichCorrectionVector() {}
92
93InitStatus CbmRichCorrectionVector::Init()
94{
95 FairRootManager* manager = FairRootManager::Instance();
96
97 fRichHits = (TClonesArray*) manager->GetObject("RichHit");
98 if (NULL == fRichHits) {
99 Fatal("CbmRichCorrectionVector::Init", "No RichHit array !");
100 }
101
102 fRichRings = (TClonesArray*) manager->GetObject("RichRing");
103 if (NULL == fRichRings) {
104 Fatal("CbmRichCorrectionVector::Init", "No RichRing array !");
105 }
106
107 fRichProjections = (TClonesArray*) manager->GetObject("RichProjection");
108 if (NULL == fRichProjections) {
109 Fatal("CbmRichCorrectionVector::Init", "No RichProjection array !");
110 }
111
112 fRichMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
113 if (NULL == fRichMirrorPoints) {
114 Fatal("CbmRichCorrectionVector::Init", "No RichMirrorPoints array !");
115 }
116
117 fRichMCPoints = (TClonesArray*) manager->GetObject("RichPoint");
118 if (NULL == fRichMCPoints) {
119 Fatal("CbmRichCorrectionVector::Init", "No RichMCPoints array !");
120 }
121
122 fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
123 if (NULL == fMCTracks) {
124 Fatal("CbmRichCorrectionVector::Init", "No MCTracks array !");
125 }
126
127 fRichRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
128 if (NULL == fRichRingMatches) {
129 Fatal("CbmRichCorrectionVector::Init", "No RichRingMatches array !");
130 }
131
132 fRichRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
133 if (NULL == fRichRefPlanePoints) {
134 Fatal("CbmRichCorrectionVector::Init", "No RichRefPlanePoint array !");
135 }
136
137 fRichPoints = (TClonesArray*) manager->GetObject("RichPoint");
138 if (NULL == fRichPoints) {
139 Fatal("CbmRichCorrectionVector::Init", "No RichPoint array !");
140 }
141
142 fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
143 if (NULL == fGlobalTracks) {
144 Fatal("CbmAnaDielectronTask::Init", "No GlobalTrack array!");
145 }
146
147 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
148 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_53"] = "2_53";
149 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
150 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_52"] = "2_52";
151 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
152 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_1_51"] = "1_51";
153 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
154 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_77"] = "2_77";
155 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
156 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_76"] = "2_76";
157 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
158 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_1_75"] = "1_75";
159 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
160 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_29"] = "2_29";
161 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
162 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_28"] = "2_28";
163 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
164 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_1_27"] = "1_27";
165 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
166 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_3_15"] = "3_15";
167 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
168 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_16"] = "2_16";
169 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
170 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_17"] = "2_17";
171
172 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
173 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_53"] = "2_53";
174 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
175 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_52"] = "2_52";
176 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
177 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_1_51"] = "1_51";
178 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
179 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_77"] = "2_77";
180 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
181 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_76"] = "2_76";
182 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
183 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_1_75"] = "1_75";
184 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
185 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_29"] = "2_29";
186 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
187 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_28"] = "2_28";
188 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
189 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_1_27"] = "1_27";
190 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
191 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_3_15"] = "3_15";
192 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
193 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_16"] = "2_16";
194 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
195 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_17"] = "2_17";
196
197 fCopFit = new CbmRichRingFitterCOP();
198 fTauFit = new CbmRichRingFitterEllipseTau();
199 CbmRichConverter::Init();
200
201 InitHistProjection();
202 InitHistAlignment();
203
204 return kSUCCESS;
205}
206
207void CbmRichCorrectionVector::InitHistProjection()
208{
209 fHM = new CbmHistManager();
210 for (std::map<string, string>::iterator it = fPathsMap.begin(); it != fPathsMap.end();
211 ++it) { // Initialize all the histograms, using map IDs as inputs.
212 string name = "fHMCPoints_" + it->second; // it->first gives the paths; it->second gives the ID.
213 fHM->Create2<TH2D>(name, name + ";X_Axis [];Y_Axis [];Entries", 2001, -100., 100., 2001, 60., 210.);
214 }
215
216 for (std::map<string, string>::iterator it = fPathsMapEllipse.begin(); it != fPathsMapEllipse.end(); ++it) {
217 string name = "fHPoints_Ellipse_" + it->second;
218 fHM->Create2<TH2D>(name, name + ";X_Axis [];Y_Axis [];Entries", 2001, -100., 100., 2001, 60., 210.);
219 }
220
221 fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrack",
222 "fhDistanceCenterToExtrapolatedTrack;Distance fitted "
223 "center to extrapolated track;Number of entries",
224 400, 0., 2.);
225 fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrackInPlane",
226 "fhDistanceCenterToExtrapolatedTrackInPlane;Distance fitted center to "
227 "extrapolated track in plane;Number of entries",
228 400, 0., 2.);
229 fHM->Create1<TH1D>("fhDifferenceX",
230 "fhDifferenceX;Difference in X (fitted center - "
231 "extrapolated track);Number of entries",
232 400, 0., 2.);
233 fHM->Create1<TH1D>("fhDifferenceY",
234 "fhDifferenceY;Difference in Y (fitted center - "
235 "extrapolated track);Number of entries",
236 400, 0., 2.);
237
238 fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrackInPlaneUncorrected",
239 "fhDistanceCenterToExtrapolatedTrackInPlaneUncorrected;Distance fitted "
240 "center to extrapolated track in plane;Number of entries",
241 300, 0., 3.);
242 fHM->Create1<TH1D>("fhDifferenceXUncorrected",
243 "fhDifferenceXUncorrected;Difference in X (fitted center "
244 "- extrapolated track);Number of entries",
245 300, 0., 3.);
246 fHM->Create1<TH1D>("fhDifferenceYUncorrected",
247 "fhDifferenceYUncorrected;Difference in Y (fitted center "
248 "- extrapolated track);Number of entries",
249 300, 0., 3.);
250}
251
252void CbmRichCorrectionVector::InitHistAlignment()
253{
254 fHM2 = new CbmHistManager();
255
256 fHM2->Create1<TH1D>("fHCenterDistance", "fHCenterDistance;Distance C-C';Nb of entries", 100, -0.1, 5.);
257 fHM2->Create1<TH1D>("fHPhi", "fHPhi;Phi_Ch [rad];Nb of entries", 200, -3.4, 3.4);
258 fHM2->Create1<TH1D>("fHThetaDiff", "fHThetaDiff;Th_Ch-Th_0 [cm];Nb of entries", 252, -5., 5.);
259 fHM2->Create2<TH2D>("fHCherenkovHitsDistribTheta0", "fHCherenkovHitsDistribTheta0;Phi_0 [rad];Theta_0 [cm];Entries",
260 200, -2., 2., 600, 2., 8.);
261 fHM2->Create2<TH2D>("fHCherenkovHitsDistribThetaCh",
262 "fHCherenkovHitsDistribThetaCh;Phi_Ch [rad];Theta_Ch [cm];Entries", 200, -3.4, 3.4, 600, 0., 20);
263 fHM2->Create2<TH2D>("fHCherenkovHitsDistribReduced",
264 "fHCherenkovHitsDistribReduced;Phi_Ch [rad];Th_Ch-Th_0 [cm];Entries", 200, -3.4, 3.4, 500, -5.,
265 5.);
266}
267
268void CbmRichCorrectionVector::Exec(Option_t* /*option*/)
269{
270 cout << endl
271 << "//"
272 "--------------------------------------------------------------------"
273 "------------------------------------------//"
274 << endl;
275 cout << "//---------------------------------------- EXEC Function - Defining "
276 "the entries ----------------------------------------//"
277 << endl;
278 cout << "//"
279 "--------------------------------------------------------------------"
280 "--------------------------------------------------//"
281 << endl;
282 fEventNum++;
283 //LOG(debug2) << "CbmRichCorrectionVector : Event #" << fEventNum;
284 cout << "CbmRichCorrectionVector : Event #" << fEventNum << endl;
285
286 Int_t nofRingsInEvent = fRichRings->GetEntriesFast();
287 Int_t nofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
288 Int_t nofHitsInEvent = fRichHits->GetEntriesFast();
289 Int_t NofMCPoints = fRichMCPoints->GetEntriesFast();
290 Int_t NofMCTracks = fMCTracks->GetEntriesFast();
291 cout << "Nb of rings in evt = " << nofRingsInEvent << ", nb of mirror points = " << nofMirrorPoints
292 << ", nb of hits in evt = " << nofHitsInEvent << ", nb of Monte-Carlo points = " << NofMCPoints
293 << " and nb of Monte-Carlo tracks = " << NofMCTracks << endl
294 << endl;
295
296 cout << "//"
297 "--------------------------------------------------------------------"
298 "--------------------------------------//"
299 << endl;
300 cout << "//-------------------------------- EXEC Function - Correction "
301 "Vector class ------------------------------------------//"
302 << endl;
303 cout << "//"
304 "--------------------------------------------------------------------"
305 "--------------------------------------//"
306 << endl
307 << endl;
308
309 TClonesArray* projectedPoint;
310
311 if (nofRingsInEvent == 0) {
312 cout << "Error no rings registered in event." << endl << endl;
313 }
314 else {
315 CalculateAnglesAndDrawDistrib();
316 //MatchFinder(); // PMT Mapping
317 //fGP = CbmRichHitProducer::InitGeometry();
318 //fGP.Print();
319 //ProjectionProducer(projectedPoint);
320 }
321}
322
323void CbmRichCorrectionVector::CalculateAnglesAndDrawDistrib()
324{
325 Double_t trackX = 0., trackY = 0.;
326 GetTrackPosition(trackX, trackY);
327
328 Int_t nofRingsInEvent = fRichRings->GetEntriesFast();
329 Float_t DistCenters, Theta_Ch, Theta_0, Angles_0;
330 Float_t Pi = 3.14159265;
331 Float_t TwoPi = 2. * 3.14159265;
332 fPhi.resize(kMAX_NOF_HITS);
333
334 // ------------------------- Loop to get ring center coordinates and photon hit coordinates per ring and per event ------------------------- //
335 if (nofRingsInEvent >= 1) {
336 cout << "Number of Rings in event: " << nofRingsInEvent << endl;
337 //sleep(2);
338 for (Int_t iR = 0; iR < nofRingsInEvent; iR++) {
339 // ----- Convert Ring to Ring Light ----- //
340 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
341 CbmRichRingLight ringL;
342 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
343 fCopFit->DoFit(&ringL);
344 // ----- Distance between mean center and fitted center calculation ----- //
345 DistCenters = sqrt(TMath::Power(ringL.GetCenterX() - trackX, 2) + TMath::Power(ringL.GetCenterY() - trackY, 2));
346 fHM2->H1("fHCenterDistance")->Fill(DistCenters);
347 // ----- Declaration of new variables ----- //
348 Int_t NofHits = ringL.GetNofHits();
349 Float_t xRing = ringL.GetCenterX();
350 Float_t yRing = ringL.GetCenterY();
351
352 for (Int_t iH = 0; iH < NofHits; iH++) {
353 // ----- Phi angle calculation ----- //
354 Int_t HitIndex = ring->GetHit(iH);
355 CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(HitIndex));
356 Float_t xHit = hit->GetX();
357 Float_t yHit = hit->GetY();
358 Angles_0 = TMath::ATan2((hit->GetX() - ringL.GetCenterX()),
359 (ringL.GetCenterY() - hit->GetY())); //* TMath::RadToDeg();
360 //cout << "Angles_0 = " << Angles_0[iH] << endl;
361
362 if (xRing - xHit == 0 || yRing - yHit == 0) continue;
363 fPhi[iH] = TMath::ATan2(yHit - yRing, xHit - xRing);
364 fHM2->H1("fHPhi")->Fill(fPhi[iH]);
365
366 // ----- Theta_Ch and Theta_0 calculations ----- //
367 Theta_Ch = sqrt(TMath::Power(trackX - hit->GetX(), 2) + TMath::Power(trackY - hit->GetY(), 2));
368 Theta_0 =
369 sqrt(TMath::Power(ringL.GetCenterX() - hit->GetX(), 2) + TMath::Power(ringL.GetCenterY() - hit->GetY(), 2));
370 //cout << "Theta_0 = " << Theta_0 << endl;
371 fHM2->H1("fHThetaDiff")->Fill(Theta_Ch - Theta_0);
372
373 // ----- Filling of final histograms ----- //
374 fHM2->H2("fHCherenkovHitsDistribTheta0")->Fill(Angles_0, Theta_0);
375 fHM2->H2("fHCherenkovHitsDistribThetaCh")->Fill(fPhi[iH], Theta_Ch);
376 fHM2->H2("fHCherenkovHitsDistribReduced")->Fill(fPhi[iH], (Theta_Ch - Theta_0));
377 }
378 //cout << endl;
379 }
380 }
381}
382
383void CbmRichCorrectionVector::GetTrackPosition(Double_t& x, Double_t& y)
384{
385 Int_t NofProjections = fRichProjections->GetEntriesFast();
386 //cout << "!!! NB PROJECTIONS !!! " << NofProjections << endl;
387 for (Int_t iP = 0; iP < NofProjections; iP++) {
388 FairTrackParam* pr = (FairTrackParam*) fRichProjections->At(iP);
389 if (NULL == pr) {
390 x = 0.;
391 y = 0.;
392 cout << "Error: CbmRichCorrectionVector::GetTrackPosition. No fair track "
393 "param found."
394 << endl;
395 }
396 x = pr->GetX();
397 y = pr->GetY();
398 //cout << "Center X: " << *x << " and Center y: " << *y << endl;
399 }
400}
401
402void CbmRichCorrectionVector::MatchFinder()
403{
404 cout << "//---------------------------------------- MATCH_FINDER Function "
405 "----------------------------------------//"
406 << endl;
407 Int_t NofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
408 Int_t NofRingsInEvent = fRichRings->GetEntriesFast();
409 Int_t NofMCPoints = fRichMCPoints->GetEntriesFast();
410 Int_t NofMCTracks = fMCTracks->GetEntriesFast();
411 //Int_t NofProjections = fRichProjections->GetEntriesFast();
412
413 Double_t x_Mirr = 0, y_Mirr = 0, z_Mirr = 0, x_PMT = 0, y_PMT = 0, z_PMT = 0;
414 Double_t CenterX = 0, CenterY = 0;
415 TGeoNode* mirr_node;
416 const Char_t* mirr_path;
417 UShort_t pHit;
418
419 for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
420 CbmRichPoint* MirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
421 Int_t trackID = MirrPoint->GetTrackID();
422 //cout << "Track ID = " << trackID << endl;
423
424 for (Int_t iMCPoint = 0; iMCPoint < NofMCPoints; iMCPoint++) {
425 CbmRichPoint* pPoint = (CbmRichPoint*) fRichMCPoints->At(iMCPoint);
426 if (NULL == pPoint) continue;
427 CbmMCTrack* pTrack = (CbmMCTrack*) fMCTracks->At(pPoint->GetTrackID());
428 if (NULL == pTrack) continue;
429 Int_t gcode = pTrack->GetPdgCode();
430 Int_t motherID = pTrack->GetMotherId();
431 if (motherID == -1) continue;
432
433 if (trackID == motherID) {
434 //cout << "MATCH BETWEEN TRACK ID AND MOTHER ID FOUND !" << endl;
435 //sleep(2);
436 //cout << "TrackID from mirror point = " << trackID << " and mother ID from MC point = " << motherID << endl;
437 x_Mirr = MirrPoint->GetX();
438 y_Mirr = MirrPoint->GetY();
439 z_Mirr = MirrPoint->GetZ();
440 //cout << "x_Mirr: " << x_Mirr << ", y_Mirr: " << y_Mirr << " and z_Mirr: " << z_Mirr << endl;
441 mirr_node = gGeoManager->FindNode(x_Mirr, y_Mirr, z_Mirr);
442 mirr_path = gGeoManager->GetPath();
443 FillPMTMap(mirr_path, pPoint);
444 //break;
445 }
446 }
447
448 // Loop filling the PMT map with ellipse points
449 for (Int_t iRing = 0; iRing < NofRingsInEvent; iRing++) {
450 CbmRichRing* ring = (CbmRichRing*) fRichRings->At(iRing);
451 if (NULL == ring) continue;
452 CbmRichRingLight ringL;
453 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
454 //fCopFit->DoFit(&ringL);
455 fTauFit->DoFit(&ringL);
456 CenterX = ringL.GetCenterX();
457 CenterY = ringL.GetCenterY();
458 CbmTrackMatchNew* richRingMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iRing));
459 Int_t richMCID = richRingMatch->GetMatchedLink().GetIndex();
460 //Int_t trackID2 = ring->GetTrackID();
461 //cout << "Track ID from ring = " << richMCID << endl;
462 if (richMCID == -1) continue;
463
464 if (trackID == richMCID) {
465
466 cout << "MATCH BETWEEN TRACK_ID AND RICH_MC_ID FOUND !" << endl;
467 cout << "Center X = " << CenterX << " and center Y = " << CenterY << endl;
468 x_Mirr = MirrPoint->GetX();
469 y_Mirr = MirrPoint->GetY();
470 z_Mirr = MirrPoint->GetZ();
471 cout << "x_Mirr: " << x_Mirr << ", y_Mirr: " << y_Mirr << " and z_Mirr: " << z_Mirr << endl;
472 mirr_node = gGeoManager->FindNode(x_Mirr, y_Mirr, z_Mirr);
473 mirr_path = gGeoManager->GetPath();
474 cout << "Mirror path: " << mirr_path << endl;
475
476 FillPMTMapEllipse(mirr_path, ringL.GetCenterX(), ringL.GetCenterY());
477 //break;
478 }
479 }
480 }
481}
482
483void CbmRichCorrectionVector::FillPMTMap(const Char_t* mirr_path, CbmRichPoint* pPoint)
484{
485 //cout << "//---------------------------------------- FILL_PMT_MAP Function ----------------------------------------//" << endl;
486 string name = string(mirr_path);
487 for (std::map<string, string>::iterator it = fPathsMap.begin(); it != fPathsMap.end(); ++it) {
488 if (name.compare(it->first) == 0) {
489 //cout << "IDENTICAL PATHS FOUND !" << endl;
490 //cout << "Mirror ID: " << it->second << endl;
491 Double_t x_PMT = pPoint->GetX();
492 Double_t y_PMT = pPoint->GetY();
493 Double_t z_PMT = pPoint->GetZ();
494 //cout << "x_PMT: " << x_PMT << ", y_PMT: " << y_PMT << " and z_PMT: " << z_PMT << endl << endl;
495 //sleep(1);
496 fHM->H2("fHMCPoints_" + it->second)->Fill(-x_PMT, y_PMT);
497 fMirrCounter++;
498 }
499 }
500}
501
502void CbmRichCorrectionVector::FillPMTMapEllipse(const Char_t* mirr_path, Float_t CenterX, Float_t CenterY)
503{
504 cout << "//---------------------------------------- FILL_PMT_MAP_ELLIPSE "
505 "Function ----------------------------------------//"
506 << endl;
507 string name = string(mirr_path);
508 for (std::map<string, string>::iterator it = fPathsMap.begin(); it != fPathsMap.end(); ++it) {
509 if (name.compare(it->first) == 0) {
510 cout << "IDENTICAL PATHS FOUND !" << endl;
511 //sleep(2);
512 cout << "Mirror ID: " << it->second << endl;
513 //cout << "Center X: " << CenterX << " and Center Y: " << CenterY << endl << endl;
514 fHM->H2("fHPoints_Ellipse_" + it->second)->Fill(-CenterX, CenterY);
515 }
516 }
517 cout << endl;
518}
519
520void CbmRichCorrectionVector::ProjectionProducer(TClonesArray* projectedPoint)
521{
522 cout << "//------------------------------ CbmRichCorrectionVector: "
523 "Projection Producer ------------------------------//"
524 << endl
525 << endl;
526
527 Int_t NofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
528 Int_t NofRingsInEvent = fRichRings->GetEntriesFast();
529 Int_t NofGTracks = fGlobalTracks->GetEntriesFast();
530 Int_t NofRefPlanePoints = fRichRefPlanePoints->GetEntriesFast();
531 Int_t NofPMTPoints = fRichPoints->GetEntriesFast();
532
533 if (fIsReconstruction) {
534 projectedPoint->Delete();
535 }
536 TMatrixFSym covMat(5);
537 for (Int_t iMatrix = 0; iMatrix < 5; iMatrix++) {
538 for (Int_t jMatrix = 0; jMatrix <= iMatrix; jMatrix++) {
539 covMat(iMatrix, jMatrix) = 0;
540 }
541 }
542 covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) = 1.e-4;
543
544 // Declaration of points coordinates.
545 Double_t sphereRadius = 300., constantePMT = 0.;
546 vector<Double_t> ptM(3), ptC(3), ptR1(3), momR1(3), normalPMT(3), ptR2Mirr(3), ptR2Center(3), ptPMirr(3), ptPR2(3);
547 vector<Double_t> ptMUnCorr(3), ptCUnCorr(3), ptR2CenterUnCorr(3), ptR2MirrUnCorr(3), ptPMirrUnCorr(3), ptPR2UnCorr(3);
548 Double_t reflectedPtCooVectSphereUnity[] = {0., 0., 0.};
549 TVector3 outPos, outPosUnCorr;
550 // Declaration of ring parameters.
551 Double_t ringCenter[] = {0., 0., 0.}, distToExtrapTrackHit = 0., distToExtrapTrackHitInPlane = 0.;
552 //Declarations related to geometry.
553 Int_t mirrTrackID = -1, pmtTrackID = -1, refPlaneTrackID = -1, motherID = -100, pmtMotherID = -100;
554 CbmMCTrack *track = NULL, *track2 = NULL;
555 TGeoNavigator* navi;
556 TGeoNode* mirrNode;
557 TGeoMatrix *mirrMatrix, *pmtMatrix, *richMatrix;
558
559 CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
560 double mirrorX = gp->fMirrorX;
561 Double_t pmtPlaneX = gp->fPmt.fPlaneX;
562 double pmtPlaneY = gp->fPmt.fPlaneY;
563 double pmtWidth = gp->fPmt.fWidth;
564 double pmtHeight = gp->fPmt.fHeight;
565
566 GetPmtNormal(NofPMTPoints, normalPMT, constantePMT);
567 cout << "Calculated normal vector to PMT plane = {" << normalPMT.at(0) << ", " << normalPMT.at(1) << ", "
568 << normalPMT.at(2) << "} and constante d = " << constantePMT << endl
569 << endl;
570
571 for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
572 //cout << "NofMirrorPoints = " << NofMirrorPoints << " and iMirr = " << iMirr << endl;
573 CbmRichPoint* mirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
574 mirrTrackID = mirrPoint->GetTrackID();
575 //cout << "Mirror track ID = " << mirrTrackID << endl;
576 if (mirrTrackID <= -1) {
577 cout << "Mirror track ID <= 1 !!!" << endl;
578 cout << "----------------------------------- End of loop N°" << iMirr + 1
579 << " on the mirror points. -----------------------------------" << endl
580 << endl;
581 continue;
582 }
583 track = (CbmMCTrack*) fMCTracks->At(mirrTrackID);
584 motherID = track->GetMotherId();
585 if (motherID == -1) {
586 //cout << "Mirror motherID == -1 !!!" << endl << endl;
587 ptM.at(0) = mirrPoint->GetX(), ptM.at(1) = mirrPoint->GetY(), ptM.at(2) = mirrPoint->GetZ();
588 ptMUnCorr.at(0) = -70.6622238, ptMUnCorr.at(1) = 55.1816025, ptMUnCorr.at(2) = 335.3621216;
589 //cout << "Mirror Point coordinates; x = " << ptM.at(0) << ", y = " << ptM.at(1) << " and z = " << ptM.at(2) << endl;
590 mirrNode = gGeoManager->FindNode(ptM.at(0), ptM.at(1), ptM.at(2));
591 if (mirrNode) {
592 cout << "Mirror node found! Mirror node name = " << mirrNode->GetName() << endl;
593 navi = gGeoManager->GetCurrentNavigator();
594 cout << "Navigator path: " << navi->GetPath() << endl;
595 cout << "Coordinates of sphere center: " << endl;
596 navi->GetCurrentMatrix()->Print();
597 if (fIsMeanCenter)
598 GetMeanSphereCenter(navi,
599 ptC); //IF NO INFORMATION ON MIRRORS ARE KNOWN (TO BE USED IN RECONSTRUCTION STEP) !!!
600 else {
601 ptC.at(0) = navi->GetCurrentMatrix()->GetTranslation()[0];
602 ptC.at(1) = navi->GetCurrentMatrix()->GetTranslation()[1];
603 ptC.at(2) = navi->GetCurrentMatrix()->GetTranslation()[2];
604 }
605 ptCUnCorr.at(0) = 0., ptCUnCorr.at(1) = 132.594000, ptCUnCorr.at(2) = 54.267226;
606 cout << "Coordinates of tile center: " << endl;
607 navi->GetMotherMatrix()->Print();
608 cout << endl
609 << "Sphere center coordinates of the rotated mirror tile = {" << ptC.at(0) << ", " << ptC.at(1) << ", "
610 << ptC.at(2) << "} and sphere inner radius = " << sphereRadius << endl;
611
612 for (Int_t iRefl = 0; iRefl < NofRefPlanePoints; iRefl++) {
613 //new((*projectedPoint)[iRefl]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
614 CbmRichPoint* refPlanePoint = (CbmRichPoint*) fRichRefPlanePoints->At(iRefl);
615 refPlaneTrackID = refPlanePoint->GetTrackID();
616 //cout << "Reflective plane track ID = " << refPlaneTrackID << endl;
617 if (mirrTrackID == refPlaneTrackID) {
618 //cout << "IDENTICAL TRACK ID FOUND !!!" << endl << endl;
619 ptR1.at(0) = refPlanePoint->GetX(), ptR1.at(1) = refPlanePoint->GetY(), ptR1.at(2) = refPlanePoint->GetZ();
620 momR1.at(0) = refPlanePoint->GetPx(), momR1.at(1) = refPlanePoint->GetPy(),
621 momR1.at(2) = refPlanePoint->GetPz();
622 cout << "Reflective Plane Point coordinates = {" << ptR1.at(0) << ", " << ptR1.at(1) << ", " << ptR1.at(2)
623 << "}" << endl;
624 cout << "And reflective Plane Point momenta = {" << momR1.at(0) << ", " << momR1.at(1) << ", "
625 << momR1.at(2) << "}" << endl;
626 cout << "Mirror Point coordinates = {" << ptM.at(0) << ", " << ptM.at(1) << ", " << ptM.at(2) << "}" << endl
627 << endl;
628 cout << "Mirror Point coordinates (Uncorrected) = {" << ptMUnCorr.at(0) << ", " << ptMUnCorr.at(1) << ", "
629 << ptMUnCorr.at(2) << "}" << endl
630 << endl;
631
632 if (fIsMeanCenter) {
633 GetMirrorIntersection(ptM, ptR1, momR1, ptC, sphereRadius);
634 // From ptM: how to retrieve tile ID ???
635 // => Compare distance of ptM to tile centers
636 }
637
638 ComputeR2(ptR2Center, ptR2Mirr, ptC, ptM, ptR1);
639 ComputeR2(ptR2CenterUnCorr, ptR2MirrUnCorr, ptCUnCorr, ptM, ptR1);
640
641 ComputeP(ptPMirr, ptPR2, normalPMT, ptM, ptR2Mirr, constantePMT);
642 ComputeP(ptPMirrUnCorr, ptPR2UnCorr, normalPMT, ptM, ptR2MirrUnCorr, constantePMT);
643
644 TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
645 CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
646 cout << "New mirror points coordinates = {" << outPos.x() << ", " << outPos.y() << ", " << outPos.z() << "}"
647 << endl;
648
649 TVector3 inPosUnCorr(ptPMirrUnCorr.at(0), ptPMirrUnCorr.at(1), ptPMirrUnCorr.at(2));
650 CbmRichGeoManager::GetInstance().RotatePoint(&inPosUnCorr, &outPosUnCorr);
651 cout << "New mirror points coordinates = {" << outPosUnCorr.x() << ", " << outPosUnCorr.y() << ", "
652 << outPosUnCorr.z() << "}" << endl
653 << endl;
654
655 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
656 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
657 pmtTrackID = pmtPoint->GetTrackID();
658 CbmMCTrack* track3 = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
659 pmtMotherID = track3->GetMotherId();
660 //cout << "pmt mother ID = " << pmtMotherID << endl;
661 if (pmtMotherID == mirrTrackID) {
662 ptP[0] = pmtPoint->GetX(), ptP[1] = pmtPoint->GetY(), ptP[2] = pmtPoint->GetZ();
663 //cout << "Identical mirror track ID and PMT mother ID !!!" << endl;
664 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
665 }
666 }
667 cout << "Looking for PMT hits: end." << endl << endl;*/
668 }
669 //else { cout << "No identical track ID between mirror point and reflective plane point found ..." << endl << endl; }
670 }
671
672 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
673 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
674 cout << "PMT points = {" << pmtPoint->GetX() << ", " << pmtPoint->GetY() << ", " << pmtPoint->GetZ() << "}" << endl;
675 TVector3 inputPoint(pmtPoint->GetX(), pmtPoint->GetY(), pmtPoint->GetZ());
676 TVector3 outputPoint;
677 CbmRichHitProducer::TiltPoint(&inputPoint, &outputPoint, fGP.fPmtPhi, fGP.fPmtTheta, fGP.fPmtZOrig);
678 cout << "New PMT points after rotation of the PMT plane: " << endl;
679 cout << outputPoint.X() << "\t" << outputPoint.Y() << "\t" << outputPoint.Z() << endl;
680 }*/
681
682 FillHistProjection(outPos, outPosUnCorr, NofGTracks, normalPMT, constantePMT);
683 }
684 else {
685 cout << "Not a mother particle ..." << endl;
686 }
687 cout << "----------------------------------- "
688 << "End of loop N°" << iMirr + 1 << " on the mirror points."
689 << " -----------------------------------" << endl
690 << endl;
691 }
692 }
693}
694
695void CbmRichCorrectionVector::GetPmtNormal(Int_t NofPMTPoints, vector<Double_t>& normalPMT, Double_t& normalCste)
696{
697 //cout << endl << "//------------------------------ CbmRichCorrectionVector: Calculate PMT Normal ------------------------------//" << endl << endl;
698
699 Int_t pmtTrackID, pmtMotherID;
700 Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0., scalarProd = 0.;
701 Double_t pmtPt[] = {0., 0., 0.};
702 Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
703 CbmMCTrack* track;
704
705 /*
706 * 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.
707 * Formula used is: vect(AB) x vect(AC) = normal.
708 * Normalize the normal vector obtained and check with three random points from the PMT plane, whether the scalar product is equal to zero.
709 */
710 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
711 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
712 pmtTrackID = pmtPoint->GetTrackID();
713 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
714 pmtMotherID = track->GetMotherId();
715 a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
716 //cout << "a[0] = " << a[0] << ", a[1] = " << a[1] << " et a[2] = " << a[2] << endl;
717 break;
718 }
719 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
720 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
721 pmtTrackID = pmtPoint->GetTrackID();
722 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
723 pmtMotherID = track->GetMotherId();
724 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
725 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
726 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
727 > 7) {
728 b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
729 //cout << "b[0] = " << b[0] << ", b[1] = " << b[1] << " et b[2] = " << b[2] << endl;
730 break;
731 }
732 }
733 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
734 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
735 pmtTrackID = pmtPoint->GetTrackID();
736 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
737 pmtMotherID = track->GetMotherId();
738 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
739 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
740 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
741 > 7
742 && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2) + TMath::Power(b[1] - pmtPoint->GetY(), 2)
743 + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
744 > 7) {
745 c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
746 //cout << "c[0] = " << c[0] << ", c[1] = " << c[1] << " et c[2] = " << c[2] << endl;
747 break;
748 }
749 }
750
751 k = (b[0] - a[0]) / (c[0] - a[0]);
752 if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0 || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
753 cout << "Error in normal calculation, vect_AB and vect_AC are collinear." << endl;
754 }
755 else {
756 buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
757 buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
758 buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
759 normalPMT.at(0) =
760 buffNormX / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
761 normalPMT.at(1) =
762 buffNormY / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
763 normalPMT.at(2) =
764 buffNormZ / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
765 }
766
767 CbmRichPoint* pmtPoint1 = (CbmRichPoint*) fRichPoints->At(20);
768 scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
769 + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
770 //cout << "1st scalar product between vectAM and normale = " << scalarProd << endl;
771 // 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.
772 normalCste =
773 -1
774 * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY() + normalPMT.at(2) * pmtPoint1->GetZ());
775 CbmRichPoint* pmtPoint2 = (CbmRichPoint*) fRichPoints->At(15);
776 scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
777 + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
778 //cout << "2nd scalar product between vectAM and normale = " << scalarProd << endl;
779 CbmRichPoint* pmtPoint3 = (CbmRichPoint*) fRichPoints->At(25);
780 scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0]) + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
781 + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
782 //cout << "3nd scalar product between vectAM and normale = " << scalarProd << endl;
783}
784
785void CbmRichCorrectionVector::GetMeanSphereCenter(TGeoNavigator* navi, vector<Double_t>& ptC)
786{
787 const Char_t* topNodePath;
788 topNodePath = gGeoManager->GetTopNode()->GetName();
789 cout << "Top node path: " << topNodePath << endl;
790 TGeoVolume* rootTop;
791 rootTop = gGeoManager->GetTopVolume();
792 rootTop->Print();
793
794 TGeoIterator nextNode(rootTop);
795 TGeoNode* curNode;
796 const TGeoMatrix* curMatrix;
797 const Double_t* curNodeTranslation; // 3 components - pointers to some memory which is provided by ROOT
798 const Double_t* curNodeRotationM; // 9 components - pointers to some memory which is provided by ROOT
799 TString filterName0("mirror_tile_type0");
800 TString filterName1("mirror_tile_type1");
801 TString filterName2("mirror_tile_type2");
802 TString filterName3("mirror_tile_type3");
803 TString filterName4("mirror_tile_type4");
804 TString filterName5("mirror_tile_type5");
805 Int_t counter = 0;
806 Double_t sphereXTot = 0., sphereYTot = 0., sphereZTot = 0.;
807
808 while ((curNode = nextNode())) {
809 TString nodeName(curNode->GetName());
810 TString nodePath;
811
812 // Filter using volume name, not node name
813 // But you can do 'if (nodeName.Contains("filter"))'
814 if (curNode->GetVolume()->GetName() == filterName0 || curNode->GetVolume()->GetName() == filterName1
815 || curNode->GetVolume()->GetName() == filterName2 || curNode->GetVolume()->GetName() == filterName3
816 || curNode->GetVolume()->GetName() == filterName4 || curNode->GetVolume()->GetName() == filterName5) {
817 if (curNode->GetNdaughters() == 0) {
818 // All deepest nodes of mirror tiles here (leaves)
819 // Thus we get spherical surface centers
820 nextNode.GetPath(nodePath);
821 curMatrix = nextNode.GetCurrentMatrix();
822 curNodeTranslation = curMatrix->GetTranslation();
823 curNodeRotationM = curMatrix->GetRotationMatrix();
824 printf("%s tr:\t", nodePath.Data());
825 printf("%08f\t%08f\t%08f\t\n", curNodeTranslation[0], curNodeTranslation[1], curNodeTranslation[2]);
826 if (curNodeTranslation[1] > 0) { // CONDITION FOR UPPER MIRROR WALL STUDY
827 sphereXTot += curNodeTranslation[0];
828 sphereYTot += curNodeTranslation[1];
829 sphereZTot += curNodeTranslation[2];
830 counter++;
831 }
832 }
833 }
834 }
835 ptC.at(0) = sphereXTot / counter;
836 ptC.at(1) = sphereYTot / counter;
837 ptC.at(2) = sphereZTot / counter;
838
839 counter = 0;
840 nextNode.Reset();
841}
842
843void CbmRichCorrectionVector::GetMirrorIntersection(vector<Double_t>& ptM, vector<Double_t> ptR1,
844 vector<Double_t> momR1, vector<Double_t> ptC, Double_t sphereRadius)
845{
846 Double_t a = 0., b = 0., c = 0., d = 0., k0 = 0., k1 = 0., k2 = 0.;
847
848 a = TMath::Power(momR1.at(0), 2) + TMath::Power(momR1.at(1), 2) + TMath::Power(momR1.at(2), 2);
849 b = 2
850 * (momR1.at(0) * (ptR1.at(0) - ptC.at(0)) + momR1.at(1) * (ptR1.at(1) - ptC.at(1))
851 + momR1.at(2) * (ptR1.at(2) - ptC.at(2)));
852 c = TMath::Power(ptR1.at(0) - ptC.at(0), 2) + TMath::Power(ptR1.at(1) - ptC.at(1), 2)
853 + TMath::Power(ptR1.at(2) - ptC.at(2), 2) - TMath::Power(sphereRadius, 2);
854 d = b * b - 4 * a * c;
855 cout << "d = " << d << endl;
856
857 if (d < 0) {
858 cout << "Error no solution to degree 2 equation found ; discriminant below 0." << endl;
859 ptM.at(0) = 0., ptM.at(1) = 0., ptM.at(2) = 0.;
860 }
861 else if (d == 0) {
862 cout << "One solution to degree 2 equation found." << endl;
863 k0 = -b / (2 * a);
864 ptM.at(0) = ptR1.at(0) + k0 * momR1.at(0);
865 ptM.at(1) = ptR1.at(1) + k0 * momR1.at(1);
866 ptM.at(2) = ptR1.at(2) + k0 * momR1.at(2);
867 }
868 else if (d > 0) {
869 cout << "Two solutions to degree 2 equation found." << endl;
870 k1 = ((-b - TMath::Sqrt(d)) / (2 * a));
871 k2 = ((-b + TMath::Sqrt(d)) / (2 * a));
872
873 if (ptR1.at(2) + k1 * momR1.at(2) > ptR1.at(2) + k2 * momR1.at(2)) {
874 ptM.at(0) = ptR1.at(0) + k1 * momR1.at(0);
875 ptM.at(1) = ptR1.at(1) + k1 * momR1.at(1);
876 ptM.at(2) = ptR1.at(2) + k1 * momR1.at(2);
877 }
878 else if (ptR1.at(2) + k1 * momR1.at(2) < ptR1.at(2) + k2 * momR1.at(2)) {
879 ptM.at(0) = ptR1.at(0) + k2 * momR1.at(0);
880 ptM.at(1) = ptR1.at(1) + k2 * momR1.at(1);
881 ptM.at(2) = ptR1.at(2) + k2 * momR1.at(2);
882 }
883 }
884}
885
886void CbmRichCorrectionVector::ComputeR2(vector<Double_t>& ptR2Center, vector<Double_t>& ptR2Mirr, vector<Double_t> ptM,
887 vector<Double_t> ptC, vector<Double_t> ptR1)
888{
889 vector<Double_t> normalMirr(3);
890 Double_t t1 = 0., t2 = 0., t3 = 0.;
891
892 normalMirr.at(0) = (ptC.at(0) - ptM.at(0))
893 / TMath::Sqrt(TMath::Power(ptC.at(0) - ptM.at(0), 2) + TMath::Power(ptC.at(1) - ptM.at(1), 2)
894 + TMath::Power(ptC.at(2) - ptM.at(2), 2));
895 normalMirr.at(1) = (ptC.at(1) - ptM.at(1))
896 / TMath::Sqrt(TMath::Power(ptC.at(0) - ptM.at(0), 2) + TMath::Power(ptC.at(1) - ptM.at(1), 2)
897 + TMath::Power(ptC.at(2) - ptM.at(2), 2));
898 normalMirr.at(2) = (ptC.at(2) - ptM.at(2))
899 / TMath::Sqrt(TMath::Power(ptC.at(0) - ptM.at(0), 2) + TMath::Power(ptC.at(1) - ptM.at(1), 2)
900 + TMath::Power(ptC.at(2) - ptM.at(2), 2));
901 //cout << "Calculated and normalized normal of mirror tile = {" << normalMirr.at(0) << ", " << normalMirr.at(1) << ", " << normalMirr.at(2) << "}" << endl;
902
903 t1 = ((ptR1.at(0) - ptM.at(0)) * (ptC.at(0) - ptM.at(0)) + (ptR1.at(1) - ptM.at(1)) * (ptC.at(1) - ptM.at(1))
904 + (ptR1.at(2) - ptM.at(2)) * (ptC.at(2) - ptM.at(2)))
905 / (TMath::Power(ptC.at(0) - ptM.at(0), 2) + TMath::Power(ptC.at(1) - ptM.at(1), 2)
906 + TMath::Power(ptC.at(2) - ptM.at(2), 2));
907 ptR2Center.at(0) = 2 * (ptM.at(0) + t1 * (ptC.at(0) - ptM.at(0))) - ptR1.at(0);
908 ptR2Center.at(1) = 2 * (ptM.at(1) + t1 * (ptC.at(1) - ptM.at(1))) - ptR1.at(1);
909 ptR2Center.at(2) = 2 * (ptM.at(2) + t1 * (ptC.at(2) - ptM.at(2))) - ptR1.at(2);
910 t2 = ((ptR1.at(0) - ptC.at(0)) * (ptC.at(0) - ptM.at(0)) + (ptR1.at(1) - ptC.at(1)) * (ptC.at(1) - ptM.at(1))
911 + (ptR1.at(2) - ptC.at(2)) * (ptC.at(2) - ptM.at(2)))
912 / (TMath::Power(ptC.at(0) - ptM.at(0), 2) + TMath::Power(ptC.at(1) - ptM.at(1), 2)
913 + TMath::Power(ptC.at(2) - ptM.at(2), 2));
914 ptR2Mirr.at(0) = 2 * (ptC.at(0) + t2 * (ptC.at(0) - ptM.at(0))) - ptR1.at(0);
915 ptR2Mirr.at(1) = 2 * (ptC.at(1) + t2 * (ptC.at(1) - ptM.at(1))) - ptR1.at(1);
916 ptR2Mirr.at(2) = 2 * (ptC.at(2) + t2 * (ptC.at(2) - ptM.at(2))) - ptR1.at(2);
917 /*//SAME AS calculation of t2 above
918 t3 = ((ptR1.at(0)-ptC.at(0))*(ptC.at(0)-ptM.at(0)) + (ptR1.at(1)-ptC.at(1))*(ptC.at(1)-ptM.at(1)) + (ptR1.at(2)-ptC.at(2))*(ptC.at(2)-ptM.at(2)))/TMath::Sqrt(TMath::Power(ptC.at(0) - ptM.at(0),2)+TMath::Power(ptC.at(1) - ptM.at(1),2)+TMath::Power(ptC.at(2) - ptM.at(2),2));
919 reflectedPtCooVectSphereUnity[0] = 2*(ptC.at(0)+t3*(normalMirr.at(0)))-ptR1.at(0);
920 reflectedPtCooVectSphereUnity[1] = 2*(ptC.at(1)+t3*(normalMirr.at(1)))-ptR1.at(1);
921 reflectedPtCooVectSphereUnity[2] = 2*(ptC.at(2)+t3*(normalMirr.at(2)))-ptR1.at(2);*/
922 cout << "Coordinates of point R2 on reflective plane after reflection on the "
923 "mirror tile:"
924 << endl;
925 //cout << "* using mirror point M to define \U00000394: {" << ptR2Center.at(0) << ", " << ptR2Center.at(1) << ", " << ptR2Center.at(2) << "}" << endl;
926 cout << "* using sphere center C to define \U00000394: {" << ptR2Mirr.at(0) << ", " << ptR2Mirr.at(1) << ", "
927 << ptR2Mirr.at(2) << "}" << endl
928 << endl;
929 //cout << "Ref Pt Coo using unity Mirror-Sphere vector & sphere pt = {" << reflectedPtCooVectSphereUnity[0] << ", " << reflectedPtCooVectSphereUnity[1] << ", " << reflectedPtCooVectSphereUnity[2] << "}" << endl << endl;
930 //cout << "NofPMTPoints = " << NofPMTPoints << endl;
931}
932
933void CbmRichCorrectionVector::ComputeP(vector<Double_t>& ptPMirr, vector<Double_t>& ptPR2, vector<Double_t> normalPMT,
934 vector<Double_t> ptM, vector<Double_t> ptR2Mirr, Double_t constantePMT)
935{
936 Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
937
938 k1 = -1
939 * ((normalPMT.at(0) * ptM.at(0) + normalPMT.at(1) * ptM.at(1) + normalPMT.at(2) * ptM.at(2) + constantePMT)
940 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0)) + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
941 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
942 ptPMirr.at(0) = ptM.at(0) + k1 * (ptR2Mirr.at(0) - ptM.at(0));
943 ptPMirr.at(1) = ptM.at(1) + k1 * (ptR2Mirr.at(1) - ptM.at(1));
944 ptPMirr.at(2) = ptM.at(2) + k1 * (ptR2Mirr.at(2) - ptM.at(2));
945 k2 = -1
946 * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1) + normalPMT.at(2) * ptR2Mirr.at(2)
947 + constantePMT)
948 / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0)) + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
949 + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
950 ptPR2.at(0) = ptR2Mirr.at(0) + k2 * (ptR2Mirr.at(0) - ptM.at(0));
951 ptPR2.at(1) = ptR2Mirr.at(1) + k2 * (ptR2Mirr.at(1) - ptM.at(1));
952 ptPR2.at(2) = ptR2Mirr.at(2) + k2 * (ptR2Mirr.at(2) - ptM.at(2));
953 cout << "Coordinates of point P on PMT plane, after reflection on the mirror "
954 "tile and extrapolation to the PMT plane:"
955 << endl;
956 cout << "* using mirror point M to define \U0001D49F ': {" << ptPMirr.at(0) << ", " << ptPMirr.at(1) << ", "
957 << ptPMirr.at(2) << "}" << endl;
958 //cout << "* using reflected point R2 to define \U0001D49F ': {" << ptPR2.at(0) << ", " << ptPR2.at(1) << ", " << ptPR2.at(2) << "}" << endl;
959 checkCalc1 =
960 ptPMirr.at(0) * normalPMT.at(0) + ptPMirr.at(1) * normalPMT.at(1) + ptPMirr.at(2) * normalPMT.at(2) + constantePMT;
961 cout << "Check whether extrapolated track point on PMT plane verifies its "
962 "equation (value should be 0.):"
963 << endl;
964 cout << "* using mirror point M, checkCalc = " << checkCalc1 << endl;
965 checkCalc2 =
966 ptPR2.at(0) * normalPMT.at(0) + ptPR2.at(1) * normalPMT.at(1) + ptPR2.at(2) * normalPMT.at(2) + constantePMT;
967 //cout << "* using reflected point R2, checkCalc = " << checkCalc2 << endl << endl;
968}
969
970void CbmRichCorrectionVector::FillHistProjection(TVector3 outPos, TVector3 outPosUnCorr, Int_t NofGTracks,
971 vector<Double_t> normalPMT, Double_t constantePMT)
972{
973 CbmMCTrack* track2 = NULL;
974 Double_t ringCenter[] = {0, 0, 0}, distToExtrapTrackHit = 0, distToExtrapTrackHitInPlane = 0;
975 Double_t distToExtrapTrackHitInPlaneUnCorr = 0;
976
977 for (Int_t iGlobalTrack = 0; iGlobalTrack < NofGTracks; iGlobalTrack++) {
978 //cout << "Nb of global tracks = " << NofGTracks << " and iGlobalTrack = " << iGlobalTrack << endl;
979 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack);
980 if (NULL == gTrack) continue;
981 Int_t richInd = gTrack->GetRichRingIndex();
982 //cout << "Rich index = " << richInd << endl;
983 if (richInd < 0) {
984 continue;
985 }
986 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richInd));
987 if (NULL == ring) {
988 continue;
989 }
990 Int_t ringTrackID = ring->GetTrackID();
991 track2 = (CbmMCTrack*) fMCTracks->At(ringTrackID);
992 Int_t ringMotherID = track2->GetMotherId();
993 //cout << "Ring mother ID = " << ringMotherID << ", ring track ID = " << ringTrackID << " and track2 pdg = " << track2->GetPdgCode() << endl;
994 CbmRichRingLight ringL;
995 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
996 //RotateAndCopyHitsToRingLight(ring, &ringL);
997 fCopFit->DoFit(&ringL);
998 //fTauFit->DoFit(&ringL);
999 ringCenter[0] = ringL.GetCenterX();
1000 ringCenter[1] = ringL.GetCenterY();
1001 ringCenter[2] =
1002 -1 * ((normalPMT.at(0) * ringCenter[0] + normalPMT.at(1) * ringCenter[1] + constantePMT) / normalPMT.at(2));
1003
1004 vector<Double_t> r(3),
1005 p(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1006 r.at(0) = TMath::Abs(ringCenter[0]), r.at(1) = TMath::Abs(ringCenter[1]), r.at(2) = TMath::Abs(ringCenter[2]);
1007 p.at(0) = TMath::Abs(outPos.x()), p.at(1) = TMath::Abs(outPos.y()), p.at(2) = TMath::Abs(outPos.z());
1008 cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}"
1009 << endl;
1010 cout << "Difference in X = " << TMath::Abs(r.at(0) - p.at(0)) << "\t ; \t"
1011 << "Difference in Y = " << TMath::Abs(r.at(1) - p.at(1)) << "\t ; \t"
1012 << "Difference in Z = " << TMath::Abs(r.at(2) - p.at(2)) << endl;
1013 fHM->H1("fhDifferenceX")->Fill(TMath::Abs(r.at(0) - p.at(0)));
1014 fHM->H1("fhDifferenceY")->Fill(TMath::Abs(r.at(1) - p.at(1)));
1015 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2)
1016 + TMath::Power(r.at(2) - p.at(2), 2));
1017 distToExtrapTrackHitInPlane = TMath::Sqrt(TMath::Power(r.at(0) - p.at(0), 2) + TMath::Power(r.at(1) - p.at(1), 2));
1018 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
1019 fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlane")->Fill(distToExtrapTrackHitInPlane);
1020 cout << "Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl;
1021 cout << "Distance between fitted ring center and extrapolated track hit in "
1022 "plane = "
1023 << distToExtrapTrackHitInPlane << endl;
1024
1025 vector<Double_t> pUnCorr(3); // Absolute coordinates of fitted ring Center r and PMT extrapolated point p
1026 pUnCorr.at(0) = TMath::Abs(outPosUnCorr.x()), pUnCorr.at(1) = TMath::Abs(outPosUnCorr.y()),
1027 pUnCorr.at(2) = TMath::Abs(outPosUnCorr.z());
1028 //cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}" << endl;
1029 cout << "Difference in X = " << TMath::Abs(r.at(0) - pUnCorr.at(0)) << "\t ; \t"
1030 << "Difference in Y = " << TMath::Abs(r.at(1) - pUnCorr.at(1)) << "\t ; \t"
1031 << "Difference in Z = " << TMath::Abs(r.at(2) - pUnCorr.at(2)) << endl;
1032 fHM->H1("fhDifferenceXUncorrected")->Fill(TMath::Abs(r.at(0) - pUnCorr.at(0)));
1033 fHM->H1("fhDifferenceYUncorrected")->Fill(TMath::Abs(r.at(1) - pUnCorr.at(1)));
1034 distToExtrapTrackHitInPlaneUnCorr =
1035 TMath::Sqrt(TMath::Power(r.at(0) - pUnCorr.at(0), 2) + TMath::Power(r.at(1) - pUnCorr.at(1), 2));
1036 fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlaneUncorrected")->Fill(distToExtrapTrackHitInPlaneUnCorr);
1037 cout << "Distance between fitted ring center and extrapolated track hit in "
1038 "plane = "
1039 << distToExtrapTrackHitInPlaneUnCorr << endl
1040 << endl;
1041 //}
1042 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
1043 }
1044 cout << "End of loop on global tracks;" << endl;
1045}
1046
1047void CbmRichCorrectionVector::RotateAndCopyHitsToRingLight(const CbmRichRing* ring1, CbmRichRingLight* ring2)
1048{
1049 Int_t nofHits = ring1->GetNofHits();
1050
1051 for (Int_t i = 0; i < nofHits; i++) {
1052 Int_t hitInd = ring1->GetHit(i);
1053 CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
1054 if (NULL == hit) continue;
1055 TVector3 inputHit(hit->GetX(), hit->GetY(), hit->GetZ());
1056 TVector3 outputHit(0, 0, 0);
1057 //CbmRichHitProducer::TiltPoint(&inputHit, &outputHit, fGP.fPmtPhi, fGP.fPmtTheta, fGP.fPmtZOrig);
1058 CbmRichHitLight hl(outputHit.X(), outputHit.Y());
1059 ring2->AddHit(hl);
1060 }
1061}
1062
1063void CbmRichCorrectionVector::DrawHistAlignment()
1064{
1065 TCanvas* c1 = new TCanvas(fRunTitle + "_Main_Analysis_" + fAxisRotTitle,
1066 fRunTitle + "_Main_Analysis_" + fAxisRotTitle, 1500, 600);
1067 c1->Divide(3, 2);
1068 c1->cd(1);
1069 /*c1->SetGridx(1);
1070 c1->SetGridy(1);
1071 fHM2->H1("fHCenterDistance")->GetXaxis()->SetLabelSize(0.03);
1072 fHM2->H1("fHCenterDistance")->GetXaxis()->SetTitleSize(0.03);
1073 fHM2->H1("fHCenterDistance")->GetXaxis()->CenterTitle();
1074 fHM2->H1("fHCenterDistance")->GetXaxis()->SetNdivisions(612,kTRUE);
1075 fHM2->H1("fHCenterDistance")->GetYaxis()->CenterTitle();
1076 fHM2->H1("fHCenterDistance")->GetYaxis()->SetTitleSize(0.03);
1077 fHM2->H1("fHCenterDistance")->GetYaxis()->SetTitleOffset(1.);
1078 fHM2->H1("fHCenterDistance")->Draw();*/
1079 DrawH1(fHM2->H1("fHCenterDistance"));
1080 c1->cd(2);
1081 DrawH1(fHM2->H1("fHThetaDiff"));
1082 c1->cd(3);
1083 DrawH2(fHM2->H2("fHCherenkovHitsDistribTheta0"));
1084 c1->cd(4);
1085 DrawH1(fHM2->H1("fHPhi"));
1086 c1->cd(5);
1087 DrawH2(fHM2->H2("fHCherenkovHitsDistribThetaCh"));
1088 c1->cd(6);
1089 DrawH2(fHM2->H2("fHCherenkovHitsDistribReduced"));
1090 Cbm::SaveCanvasAsImage(c1, string(fOutputDir.Data()), "png");
1091}
1092
1093void CbmRichCorrectionVector::DrawFit(vector<Double_t>& outputFit, Int_t thresh)
1094{
1095 //vector<Double_t> paramVect;
1096 //paramVect.reserve(5);
1097
1098 TCanvas* c3 =
1099 new TCanvas(fRunTitle + "_Fit_Slices_" + fAxisRotTitle, fRunTitle + "_Fit_Slices_" + fAxisRotTitle, 1100, 600);
1100 c3->SetFillColor(42);
1101 c3->Divide(4, 2);
1102 gPad->SetTopMargin(0.1);
1103 gPad->SetFillColor(33);
1104 c3->cd(1);
1105 // TH2D* CloneArr = (TH2D*)fHM2->H2("fHCherenkovHitsDistribThetaCh")->Clone();
1106 TH2D* CloneArr = (TH2D*) fHM2->H2("fHCherenkovHitsDistribReduced")->Clone();
1107 CloneArr->GetXaxis()->SetLabelSize(0.03);
1108 CloneArr->GetXaxis()->SetTitleSize(0.03);
1109 CloneArr->GetXaxis()->CenterTitle();
1110 CloneArr->GetXaxis()->SetNdivisions(612, kTRUE);
1111 CloneArr->GetYaxis()->SetLabelSize(0.03);
1112 CloneArr->GetYaxis()->SetTitleSize(0.03);
1113 CloneArr->GetYaxis()->SetNdivisions(612, kTRUE);
1114 CloneArr->GetYaxis()->CenterTitle();
1115 //CloneArr->GetYaxis()->SetRangeUser(-2.5,2.5);
1116 CloneArr->GetZaxis()->SetLabelSize(0.03);
1117 CloneArr->GetZaxis()->SetTitleSize(0.03);
1118 CloneArr->GetYaxis()->SetTitleOffset(1.0);
1119 CloneArr->Draw("colz");
1120 //Double_t ymax = CloneArr->GetYaxis()->GetXmax();
1121 //Double_t ymin = CloneArr->GetYaxis()->GetXmin();
1122 //TF1 *fgauss = TF1 fgauss("gauss", "[0]*exp(-0.5*((x-[1])/[2])**2)", 0, 100);
1123
1124 // ------------------------------ APPLY THRESHOLD TO 2D-HISTO ------------------------------ //
1125 // TH2D* CloneArr_2 = (TH2D*)fHM2->H2("fHCherenkovHitsDistribThetaCh")->Clone();
1126 TH2D* CloneArr_2 = (TH2D*) fHM2->H2("fHCherenkovHitsDistribReduced")->Clone();
1127 for (Int_t y_bin = 1; y_bin <= 500; y_bin++) {
1128 for (Int_t x_bin = 1; x_bin <= 200; x_bin++) {
1129 /*if (CloneArr_2->GetBinContent(x_bin, y_bin)!=0) {
1130 cout << "Bin Content: " << CloneArr_2->GetBinContent(x_bin, y_bin) << endl;
1131 sleep(1);
1132 }
1133 else;*/
1134 if (CloneArr_2->GetBinContent(x_bin, y_bin) < thresh) {
1135 CloneArr_2->SetBinContent(x_bin, y_bin, 0);
1136 }
1137 }
1138 }
1139 c3->cd(2);
1140 CloneArr_2->GetXaxis()->SetLabelSize(0.03);
1141 CloneArr_2->GetXaxis()->SetTitleSize(0.03);
1142 CloneArr_2->GetXaxis()->CenterTitle();
1143 CloneArr_2->GetXaxis()->SetNdivisions(612, kTRUE);
1144 CloneArr_2->GetYaxis()->SetLabelSize(0.03);
1145 CloneArr_2->GetYaxis()->SetTitleSize(0.03);
1146 CloneArr_2->GetYaxis()->SetNdivisions(612, kTRUE);
1147 CloneArr_2->GetYaxis()->CenterTitle();
1148 //CloneArr_2->GetYaxis()->SetRangeUser(-2.5,2.5);
1149 CloneArr_2->GetZaxis()->SetLabelSize(0.03);
1150 CloneArr_2->GetZaxis()->SetTitleSize(0.03);
1151 CloneArr_2->GetYaxis()->SetTitleOffset(1.0);
1152 CloneArr_2->Draw("colz");
1153 CloneArr_2->Write();
1154
1155 // -------------------- FIT SLICES AND FIT THE MEAN OF THE RESULT TO A SIN FUNCTION -------------------- //
1156 CloneArr_2->FitSlicesY(0, 0, -1, 1);
1157 c3->cd(3);
1158 TH1D* histo_0 = gDirectory->Get<TH1D>("fHCherenkovHitsDistribReduced_0");
1159 histo_0->Draw();
1160 c3->cd(4);
1161 TH1D* histo_1 = gDirectory->Get<TH1D>("fHCherenkovHitsDistribReduced_1");
1162 //histo_1->GetYaxis()->SetRangeUser(-2.5, 2.5);
1163 histo_1->Draw();
1164 c3->cd(5);
1165 TH1D* histo_2 = gDirectory->Get<TH1D>("fHCherenkovHitsDistribReduced_2");
1166 histo_2->Draw();
1167 c3->cd(6);
1168 TH1D* histo_chi2 = gDirectory->Get<TH1D>("fHCherenkovHitsDistribReduced_chi2");
1169 histo_chi2->Draw();
1170
1171 c3->cd(7);
1172 TF1* f1 = new TF1("f1", "[2]+[0]*cos(x)+[1]*sin(x)", -3.5, 3.5);
1173 f1->SetParameters(0, 0, 0);
1174 f1->SetParNames("Delta_phi", "Delta_lambda", "Offset");
1175 histo_1->Fit("f1", "", "");
1176 TF1* fit = histo_1->GetFunction("f1");
1177 Double_t p1 = fit->GetParameter("Delta_phi");
1178 Double_t p2 = fit->GetParameter("Delta_lambda");
1179 Double_t p3 = fit->GetParameter("Offset");
1180 Double_t chi2 = fit->GetChisquare();
1181 //cout << setprecision(6) << "Delta_phi = " << fit->GetParameter(0) << " and delta_lambda = " << fit->GetParameter(1) << endl;
1182 //cout << "Delta_phi error = " << fit->GetParError(0) << " and delta_lambda error = " << fit->GetParError(1) << endl;
1183 //cout << endl << "Chi2: " << chi2 << endl;
1184
1185 /* paramVect.push_back(fit->GetParameter("Delta_phi"));
1186 paramVect.push_back(fit->GetParameter("Delta_lambda"));
1187 paramVect.push_back(fit->GetParameter("Offset"));
1188 paramVect.push_back(fit->GetChisquare());
1189 //cout << "Vectors: Delta_phi = " << paramVect[0] << ", Delta_lambda = " << paramVect[1] << ", Offset = " << paramVect[2] << endl;
1190*/
1191 f1->SetParameters(fit->GetParameter(0), fit->GetParameter(1));
1192 char leg[128];
1193 f1->SetLineColor(2);
1194 f1->Draw();
1195 f1->Write();
1196
1197 // ------------------------------ CALCULATION OF MISALIGNMENT ANGLE ------------------------------ //
1198 Double_t Focal_length = 150., q = 0., A = 0., Alpha = 0., mis_x = 0., mis_y = 0.;
1199 q = TMath::ATan(fit->GetParameter(0) / fit->GetParameter(1));
1200 cout << endl << "fit_1 = " << fit->GetParameter(0) << " and fit_2 = " << fit->GetParameter(1) << endl;
1201 //cout << "q = " << q << endl;
1202 A = fit->GetParameter(1) / TMath::Cos(q);
1203 //cout << "Parameter a = " << A << endl;
1204 Alpha =
1205 TMath::ATan(A / 1.5) * 0.5
1206 * TMath::Power(
1207 10,
1208 3); // *0.5, because a mirror rotation of alpha implies a rotation in the particle trajectory of 2*alpha ; 1.5 meters = Focal length = Radius_of_curvature/2
1209 //cout << setprecision(6) << "Total angle of misalignment alpha = " << Alpha << endl; // setprecision(#) gives the number of digits in the cout.
1210 mis_x = TMath::ATan(fit->GetParameter(0) / Focal_length) * 0.5 * TMath::Power(10, 3);
1211 mis_y = TMath::ATan(fit->GetParameter(1) / Focal_length) * 0.5 * TMath::Power(10, 3);
1212 //cout << "Horizontal displacement = " << outputFit[0] << " [mrad] and vertical displacement = " << outputFit[1] << " [mrad]." << endl;
1213
1214 TLegend* LEG = new TLegend(0.30, 0.7, 0.72, 0.85); // Set legend position
1215 LEG->SetBorderSize(1);
1216 LEG->SetFillColor(0);
1217 LEG->SetMargin(0.2);
1218 LEG->SetTextSize(0.03);
1219 sprintf(leg, "Fitted sinusoid");
1220 LEG->AddEntry(f1, leg, "l");
1221 sprintf(leg, "Misalign in X = %f", mis_x);
1222 LEG->AddEntry("", leg, "l");
1223 sprintf(leg, "Misalign in Y = %f", mis_y);
1224 LEG->AddEntry("", leg, "l");
1225 sprintf(leg, "Offset = %f", fit->GetParameter(2));
1226 LEG->AddEntry("", leg, "l");
1227 LEG->Draw();
1228 Cbm::SaveCanvasAsImage(c3, string(fOutputDir.Data()), "png");
1229
1230 TCanvas* c4 = new TCanvas(fRunTitle + fAxisRotTitle, fRunTitle + fAxisRotTitle, 400, 400);
1231 CloneArr_2->Draw("colz");
1232 f1->Draw("same");
1233 TLegend* LEG1 = new TLegend(0.30, 0.7, 0.72, 0.85); // Set legend position
1234 LEG1->SetBorderSize(1);
1235 LEG1->SetFillColor(0);
1236 LEG1->SetMargin(0.2);
1237 LEG1->SetTextSize(0.03);
1238 sprintf(leg, "Fitted sinusoid");
1239 LEG1->AddEntry(f1, leg, "l");
1240 sprintf(leg, "Misalign in X = %f", mis_x);
1241 LEG1->AddEntry("", leg, "l");
1242 sprintf(leg, "Misalign in Y = %f", mis_y);
1243 LEG1->AddEntry("", leg, "l");
1244 sprintf(leg, "Offset = %f", fit->GetParameter(2));
1245 LEG1->AddEntry("", leg, "l");
1246 LEG1->Draw();
1247
1248 // ------------------------------ APPLY SECOND FIT USING LOG-LIKELIHOOD METHOD ------------------------------ //
1249 /* TCanvas* c4 = new TCanvas(fRunTitle + "_Second_Fit_" + fAxisRotTitle, fRunTitle + "_Second_Fit_" + fAxisRotTitle, 600, 600);
1250 c4->SetFillColor(42);
1251 gPad->SetTopMargin(0.1);
1252 gPad->SetFillColor(33);
1253 f1->SetParameters(fit->GetParameter(0), fit->GetParameter(1), fit->GetParameter(2));
1254 histo_1->Fit("f1","L","");
1255 TF1 *fit2 = histo_1->GetFunction("f1");
1256 f1->SetParameters(fit2->GetParameter(0), fit2->GetParameter(1), fit2->GetParameter(2));
1257 f1->Draw();
1258
1259 Double_t q_2 = TMath::ATan(fit2->GetParameter(0)/fit2->GetParameter(1));
1260 //cout << endl << "fit2_1 = " << fit2->GetParameter(0) << " and fit2_2 = " << fit2->GetParameter(1) << endl;
1261 //cout << "q_2 = " << q_2 << endl;
1262 Double_t A_2 = fit2->GetParameter(1)/TMath::Cos(q_2);
1263 //cout << "Parameter a_2 = " << A_2 << endl;
1264 Double_t Alpha_2 = TMath::ATan(A_2/1.5)*0.5*TMath::Power(10,3);
1265 //cout << setprecision(6) << "Total angle of misalignment alpha_2 = " << Alpha_2 << endl;
1266 Double_t mis_x_2 = TMath::ATan(fit2->GetParameter(0)/Focal_length)*0.5*TMath::Power(10,3);
1267 Double_t mis_y_2 = TMath::ATan(fit2->GetParameter(1)/Focal_length)*0.5*TMath::Power(10,3);
1268
1269 TLegend* LEG2= new TLegend(0.31,0.7,0.72,0.85); // Set legend position
1270 LEG2->SetBorderSize(1);
1271 LEG2->SetFillColor(0);
1272 LEG2->SetMargin(0.2);
1273 LEG2->SetTextSize(0.03);
1274 sprintf(leg, "Fitted sinusoid");
1275 LEG2->AddEntry(f1, leg, "l");
1276 sprintf(leg, "Misalign in X = %f", mis_x_2);
1277 LEG2->AddEntry("", leg, "l");
1278 sprintf(leg, "Misalign in Y = %f", mis_y_2);
1279 LEG2->AddEntry("", leg, "l");
1280 sprintf(leg, "Offset = %f", fit2->GetParameter(2));
1281 LEG2->AddEntry("", leg, "l");
1282 LEG2->Draw();
1283 Cbm::SaveCanvasAsImage(c4, string(fOutputDir.Data()), "png");*/
1284
1285 outputFit.at(0) = mis_x;
1286 outputFit.at(1) = mis_y;
1287}
1288
1289void CbmRichCorrectionVector::DrawHistMapping()
1290{
1291 TCanvas* can = new TCanvas(fRunTitle + "_Separated_Hits", fRunTitle + "_Separated_Hits", 1500, 900);
1292 can->SetGridx(1);
1293 can->SetGridy(1);
1294 can->Divide(4, 3);
1295 can->cd(9);
1296 DrawH2(fHM->H2("fHMCPoints_3_15"));
1297 can->cd(5);
1298 DrawH2(fHM->H2("fHMCPoints_2_16"));
1299 can->cd(1);
1300 DrawH2(fHM->H2("fHMCPoints_2_17"));
1301 can->cd(2);
1302 DrawH2(fHM->H2("fHMCPoints_2_29"));
1303 can->cd(3);
1304 DrawH2(fHM->H2("fHMCPoints_2_53"));
1305 can->cd(4);
1306 DrawH2(fHM->H2("fHMCPoints_2_77"));
1307 can->cd(6);
1308 DrawH2(fHM->H2("fHMCPoints_2_28"));
1309 can->cd(7);
1310 DrawH2(fHM->H2("fHMCPoints_2_52"));
1311 can->cd(8);
1312 DrawH2(fHM->H2("fHMCPoints_2_76"));
1313 can->cd(10);
1314 DrawH2(fHM->H2("fHMCPoints_1_27"));
1315 can->cd(11);
1316 DrawH2(fHM->H2("fHMCPoints_1_51"));
1317 can->cd(12);
1318 DrawH2(fHM->H2("fHMCPoints_1_75"));
1319 Cbm::SaveCanvasAsImage(can, string(fOutputDir.Data()), "png");
1320
1321 TCanvas* can2 = new TCanvas(fRunTitle + "_Separated_Ellipse", fRunTitle + "_Separated_Ellipse", 1500, 900);
1322 can2->SetGridx(1);
1323 can2->SetGridy(1);
1324 can2->Divide(4, 3);
1325 can2->cd(9);
1326 DrawH2(fHM->H2("fHPoints_Ellipse_3_15"));
1327 can2->cd(5);
1328 DrawH2(fHM->H2("fHPoints_Ellipse_2_16"));
1329 can2->cd(1);
1330 DrawH2(fHM->H2("fHPoints_Ellipse_2_17"));
1331 can2->cd(2);
1332 DrawH2(fHM->H2("fHPoints_Ellipse_2_29"));
1333 can2->cd(3);
1334 DrawH2(fHM->H2("fHPoints_Ellipse_2_53"));
1335 can2->cd(4);
1336 DrawH2(fHM->H2("fHPoints_Ellipse_2_77"));
1337 can2->cd(6);
1338 DrawH2(fHM->H2("fHPoints_Ellipse_2_28"));
1339 can2->cd(7);
1340 DrawH2(fHM->H2("fHPoints_Ellipse_2_52"));
1341 can2->cd(8);
1342 DrawH2(fHM->H2("fHPoints_Ellipse_2_76"));
1343 can2->cd(10);
1344 DrawH2(fHM->H2("fHPoints_Ellipse_1_27"));
1345 can2->cd(11);
1346 DrawH2(fHM->H2("fHPoints_Ellipse_1_51"));
1347 can2->cd(12);
1348 DrawH2(fHM->H2("fHPoints_Ellipse_1_75"));
1349 Cbm::SaveCanvasAsImage(can2, string(fOutputDir.Data()), "png");
1350}
1351
1352void CbmRichCorrectionVector::DrawHistProjection()
1353{
1354 TCanvas* can3 = new TCanvas(fRunTitle + "_Projected_Points_w/Corr_" + fAxisRotTitle,
1355 fRunTitle + "_Projected_Points_w/Corr_" + fAxisRotTitle, 1500, 900);
1356 can3->Divide(2, 2);
1357 can3->cd(1);
1358 DrawH1(fHM->H1("fhDifferenceX"));
1359 can3->cd(2);
1360 DrawH1(fHM->H1("fhDifferenceY"));
1361 can3->cd(3);
1362 /* DrawH1(fHM->H1("fhDistanceCenterToExtrapolatedTrack"));
1363 can3->cd(4);*/
1364 DrawH1(fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlane"));
1365 Cbm::SaveCanvasAsImage(can3, string(fOutputDir.Data()), "png");
1366
1367 TCanvas* can4 = new TCanvas(fRunTitle + "_Projected_Points_w/oCorr_" + fAxisRotTitle,
1368 fRunTitle + "_Projected_Points_w/oCorr_" + fAxisRotTitle, 1500, 900);
1369 can4->Divide(2, 2);
1370 can4->cd(1);
1371 DrawH1(fHM->H1("fhDifferenceXUncorrected"));
1372 can4->cd(2);
1373 DrawH1(fHM->H1("fhDifferenceYUncorrected"));
1374 can4->cd(3);
1375 DrawH1(fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlaneUncorrected"));
1376 can4->SaveAs("Uncorrected_Histograms_" + fAxisRotTitle + ".png");
1377}
1378
1379void CbmRichCorrectionVector::DrawHistFromFile(TString fileName)
1380{
1382 TFile* oldFile = gFile;
1383 TDirectory* oldDir = gDirectory;
1384
1385 fHM = new CbmHistManager();
1386 TFile* file = new TFile(fileName, "READ");
1387 LOG_IF(fatal, !file) << "Could not open file " << fileName;
1388
1389 fHM->ReadFromFile(file);
1390 DrawHistMapping();
1391
1393 gFile = oldFile;
1394 gDirectory = oldDir;
1395}
1396
1397void CbmRichCorrectionVector::Finish()
1398{
1399 // ---------------------------------------------------------------------------------------------------------------------------------------- //
1400 // -------------------------------------------------- Mapping for mirror - PMT relations -------------------------------------------------- //
1401 // ---------------------------------------------------------------------------------------------------------------------------------------- //
1402
1403 if (fDrawAlignment) {
1404 DrawHistAlignment();
1405 Int_t thresh = 5;
1406 vector<Double_t> outputFit(2);
1407 DrawFit(outputFit, thresh);
1408 cout << setprecision(6) << endl;
1409 cout << "Horizontal displacement = " << outputFit[0] << " [mrad] and vertical displacement = " << outputFit[1]
1410 << " [mrad]." << endl;
1411 /*Double_t Focal_length = 150;
1412 Double_t q = TMath::ATan(outputFit[0]/outputFit[1]);
1413 //cout << endl << "fit_1 = " << fit->GetParameter(0) << " and fit_2 = " << fit->GetParameter(1) << endl;
1414 //cout << "q = " << q << endl;
1415 Double_t A = outputFit[1]/TMath::Cos(q);
1416 //cout << "Parameter a = " << A << endl;
1417 Double_t Alpha = TMath::ATan(A/Focal_length)*0.5*TMath::Power(10,3); // *0.5, because a mirror rotation of alpha implies a rotation in the particle trajectory of 2*alpha ; 1.5 meters = Focal length = Radius_of_curvature/2
1418 cout << endl << setprecision(6) << "Angle of misalignment alpha [mrad] = " << Alpha << endl; // setprecision(#) gives the number of digits in the cout.
1419 Double_t mis_x = TMath::ATan(outputFit[0]/Focal_length)*0.5*TMath::Power(10,3);
1420 Double_t mis_y = TMath::ATan(outputFit[1]/Focal_length)*0.5*TMath::Power(10,3);
1421 cout << "Misalignment in X [mrad] = " << mis_x << " and misalignment in Y [mrad] = " << mis_y << endl;*/
1422
1423 fHM2->Create2<TH2D>("fHCherenkovHitsDistribReduced",
1424 "fHCherenkovHitsDistribReduced;Phi_Ch [rad];Th_Ch-Th_0 [cm];Entries", 200, -3.4, 3.4, 500, -5.,
1425 5.);
1426 }
1427
1428 if (fDrawMapping) {
1429 DrawHistMapping();
1430 }
1431
1432 if (fDrawProjection) {
1433 DrawHistProjection();
1434 fHM->H1("fhDifferenceX")->Write();
1435 fHM->H1("fhDifferenceY")->Write();
1436 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Write();
1437 fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlane")->Write();
1438 fHM->H1("fhDifferenceXUncorrected")->Write();
1439 fHM->H1("fhDifferenceYUncorrected")->Write();
1440 fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlaneUncorrected")->Write();
1441 }
1442
1443 //cout << endl << "Mirror counter = " << fMirrCounter << endl;
1444 //cout << setprecision(6) << endl;
1445}
1446ClassImp(CbmRichCorrectionVector)
ClassImp
ClassImp(CbmConverterManager)
DrawH1
void DrawH1(TH1 *hist, HistScale logx, HistScale logy, const string &drawOpt, Int_t color, Int_t lineWidth, Int_t lineStyle, Int_t markerSize, Int_t markerStyle)
Definition CbmDrawHist.cxx:54
DrawH2
void DrawH2(TH2 *hist, HistScale logx, HistScale logy, HistScale logz, const string &drawOpt)
Definition CbmDrawHist.cxx:83
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.
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmHistManager
Histogram manager.
Definition CbmHistManager.h:48
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition CbmHistManager.h:243
CbmHistManager::Create2
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...
Definition CbmHistManager.h:123
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
CbmHit::GetZ
double GetZ() const
Definition CbmHit.h:71
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
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
CbmRichCorrectionVector::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 CbmRichCorrectionVector.cxx:933
CbmRichCorrectionVector::fPathsMapEllipse
std::map< string, string > fPathsMapEllipse
Definition CbmRichCorrectionVector.h:219
CbmRichCorrectionVector::fPathsMap
std::map< string, string > fPathsMap
Definition CbmRichCorrectionVector.h:218
CbmRichCorrectionVector::DrawHistFromFile
void DrawHistFromFile(TString fileName)
Definition CbmRichCorrectionVector.cxx:1379
CbmRichCorrectionVector::GetTrackPosition
void GetTrackPosition(Double_t &x, Double_t &y)
Definition CbmRichCorrectionVector.cxx:383
CbmRichCorrectionVector::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition CbmRichCorrectionVector.cxx:268
CbmRichCorrectionVector::FillPMTMapEllipse
void FillPMTMapEllipse(const Char_t *mirr_path, Float_t CenterX, Float_t CenterY)
Definition CbmRichCorrectionVector.cxx:502
CbmRichCorrectionVector::FillHistProjection
void FillHistProjection(TVector3 outPos, TVector3 outPosUnCorr, Int_t NofGlobalTracks, vector< Double_t > normalPMT, Double_t constantePMT)
Definition CbmRichCorrectionVector.cxx:970
CbmRichCorrectionVector::fCopFit
CbmRichRingFitterCOP * fCopFit
Definition CbmRichCorrectionVector.h:225
CbmRichCorrectionVector::GetMeanSphereCenter
void GetMeanSphereCenter(TGeoNavigator *navi, vector< Double_t > &ptC)
Definition CbmRichCorrectionVector.cxx:785
CbmRichCorrectionVector::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition CbmRichCorrectionVector.cxx:93
CbmRichCorrectionVector::ComputeR2
void ComputeR2(vector< Double_t > &ptR2Center, vector< Double_t > &ptR2Mirr, vector< Double_t > ptM, vector< Double_t > ptC, vector< Double_t > ptR1)
Definition CbmRichCorrectionVector.cxx:886
CbmRichCorrectionVector::fTauFit
CbmRichRingFitterEllipseTau * fTauFit
Definition CbmRichCorrectionVector.h:226
CbmRichCorrectionVector::Finish
virtual void Finish()
Inherited from FairTask.
Definition CbmRichCorrectionVector.cxx:1397
CbmRichCorrectionVector::GetMirrorIntersection
void GetMirrorIntersection(vector< Double_t > &ptM, vector< Double_t > ptR1, vector< Double_t > momR1, vector< Double_t > ptC, Double_t sphereRadius)
Definition CbmRichCorrectionVector.cxx:843
CbmRichCorrectionVector::DrawFit
void DrawFit(vector< Double_t > &outputFit, Int_t thresh)
Definition CbmRichCorrectionVector.cxx:1093
CbmRichCorrectionVector::FillPMTMap
void FillPMTMap(const Char_t *mirr_path, CbmRichPoint *pPoint)
Definition CbmRichCorrectionVector.cxx:483
CbmRichCorrectionVector::ProjectionProducer
void ProjectionProducer(TClonesArray *projectedPoint)
Definition CbmRichCorrectionVector.cxx:520
CbmRichCorrectionVector::GetPmtNormal
void GetPmtNormal(Int_t NofPMTPoints, vector< Double_t > &normalPMT, Double_t &normalCste)
Definition CbmRichCorrectionVector.cxx:695
CbmRichCorrectionVector::RotateAndCopyHitsToRingLight
void RotateAndCopyHitsToRingLight(const CbmRichRing *ring1, CbmRichRingLight *ring2)
Definition CbmRichCorrectionVector.cxx:1047
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
CbmRichRingFitterEllipseTau::DoFit
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition CbmRichRingFitterEllipseTau.cxx:98
CbmRichRingLight::GetCenterX
float GetCenterX() const
Definition CbmRichRingLight.h:165
CbmRichRingLight::GetCenterY
float GetCenterY() const
Definition CbmRichRingLight.h:166
CbmRichRingLight::GetNofHits
int GetNofHits() const
Return number of hits in ring.
Definition CbmRichRingLight.h:112
CbmRichRingLight::AddHit
void AddHit(CbmRichHitLight hit)
Add new hit to the ring.
Definition CbmRichRingLight.h:90
CbmRichRing::GetHit
uint32_t GetHit(int32_t i) const
Definition CbmRichRing.h:39
CbmRichRing::GetNofHits
int32_t GetNofHits() const
Definition CbmRichRing.h:37
Cbm::SaveCanvasAsImage
void SaveCanvasAsImage(TCanvas *c, const string &dir, const string &option)
Definition CbmUtils.cxx:36