CbmRoot
Loading...
Searching...
No Matches
CbmRichPMTMapping.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 "CbmRichPMTMapping.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 "CbmRichHitProducer.h"
44#include "TGeoSphere.h"
45#include "TLorentzVector.h"
46#include "TVirtualMC.h"
47class TGeoNode;
48class TGeoVolume;
49class TGeoShape;
50class TGeoMatrix;
51
52CbmRichPMTMapping::CbmRichPMTMapping()
53 : FairTask()
54 , fRichHits(NULL)
55 , fRichRings(NULL)
56 , fRichProjections(NULL)
57 , fRichMirrorPoints(NULL)
58 , fRichMCPoints(NULL)
59 , fMCTracks(NULL)
60 , fRichRingMatches(NULL)
61 , fRichRefPlanePoints(NULL)
62 , fRichPoints(NULL)
63 , fGlobalTracks(NULL)
64 , fHM(NULL)
65 , fGP()
66 , fEventNum(0)
67 , fOutputDir("")
68 , fRunTitle("")
69 , fDrawHist(kFALSE)
70 , fIsMirrorUpperHalf(NULL)
71 , fCopFit(NULL)
72 , fTauFit(NULL)
73 , fPathsMap()
74 , fPathsMapEllipse()
75{
76 fCounterMapping = 0.;
77 fMirrCounter = 0.;
78 for (int i = 0; i < 3; i++) {
79 fArray[i] = 0.;
80 }
81}
82
83CbmRichPMTMapping::~CbmRichPMTMapping() {}
84
85InitStatus CbmRichPMTMapping::Init()
86{
87 FairRootManager* manager = FairRootManager::Instance();
88
89 fRichHits = (TClonesArray*) manager->GetObject("RichHit");
90 if (NULL == fRichHits) {
91 Fatal("CbmRichPMTMapping::Init", "No RichHit array !");
92 }
93
94 fRichRings = (TClonesArray*) manager->GetObject("RichRing");
95 if (NULL == fRichRings) {
96 Fatal("CbmRichPMTMapping::Init", "No RichRing array !");
97 }
98
99 fRichProjections = (TClonesArray*) manager->GetObject("RichProjection");
100 if (NULL == fRichProjections) {
101 Fatal("CbmRichPMTMapping::Init", "No RichProjection array !");
102 }
103
104 fRichMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
105 if (NULL == fRichMirrorPoints) {
106 Fatal("CbmRichPMTMapping::Init", "No RichMirrorPoints array !");
107 }
108
109 fRichMCPoints = (TClonesArray*) manager->GetObject("RichPoint");
110 if (NULL == fRichMCPoints) {
111 Fatal("CbmRichPMTMapping::Init", "No RichMCPoints array !");
112 }
113
114 fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
115 if (NULL == fMCTracks) {
116 Fatal("CbmRichPMTMapping::Init", "No MCTracks array !");
117 }
118
119 fRichRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
120 if (NULL == fRichRingMatches) {
121 Fatal("CbmRichPMTMapping::Init", "No RichRingMatches array !");
122 }
123
124 fRichRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
125 if (NULL == fRichRefPlanePoints) {
126 Fatal("CbmRichPMTMapping::Init", "No RichRefPlanePoint array !");
127 }
128
129 fRichPoints = (TClonesArray*) manager->GetObject("RichPoint");
130 if (NULL == fRichPoints) {
131 Fatal("CbmRichPMTMapping::Init", "No RichPoint array !");
132 }
133
134 fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
135 if (NULL == fGlobalTracks) {
136 Fatal("CbmAnaDielectronTask::Init", "No GlobalTrack array!");
137 }
138
139 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
140 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_53"] = "2_53";
141 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
142 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_52"] = "2_52";
143 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
144 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_1_51"] = "1_51";
145 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
146 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_77"] = "2_77";
147 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
148 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_76"] = "2_76";
149 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
150 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_1_75"] = "1_75";
151 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
152 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_29"] = "2_29";
153 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
154 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_28"] = "2_28";
155 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
156 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_1_27"] = "1_27";
157 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
158 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_3_15"] = "3_15";
159 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
160 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_16"] = "2_16";
161 fPathsMap["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
162 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_17"] = "2_17";
163
164 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
165 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_53"] = "2_53";
166 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
167 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_2_52"] = "2_52";
168 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
169 "RICH_mirror_and_support_belt_strip3_126/RICH_mirror_1_51"] = "1_51";
170 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
171 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_77"] = "2_77";
172 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
173 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_2_76"] = "2_76";
174 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
175 "RICH_mirror_and_support_belt_strip5_128/RICH_mirror_1_75"] = "1_75";
176 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
177 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_29"] = "2_29";
178 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
179 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_2_28"] = "2_28";
180 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
181 "RICH_mirror_and_support_belt_strip1_124/RICH_mirror_1_27"] = "1_27";
182 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
183 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_3_15"] = "3_15";
184 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
185 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_16"] = "2_16";
186 fPathsMapEllipse["/cave_1/rich1_0/RICH_gas_221/RICH_mirror_half_total_208/"
187 "RICH_mirror_and_support_belt_strip_cut_123/RICH_mirror_2_17"] = "2_17";
188
189 fCopFit = new CbmRichRingFitterCOP();
190 fTauFit = new CbmRichRingFitterEllipseTau();
191 CbmRichConverter::Init();
192
193 InitHist();
194
195 return kSUCCESS;
196}
197
198void CbmRichPMTMapping::InitHist()
199{
200 fHM = new CbmHistManager();
201 for (std::map<string, string>::iterator it = fPathsMap.begin(); it != fPathsMap.end();
202 ++it) { // Initialize all the histograms, using map IDs as inputs.
203 string name = "fHMCPoints_" + it->second; // it->first gives the paths; it->second gives the ID.
204 fHM->Create2<TH2D>(name, name + ";X_Axis [];Y_Axis [];Entries", 2001, -100., 100., 2001, 60., 210.);
205 }
206
207 for (std::map<string, string>::iterator it = fPathsMapEllipse.begin(); it != fPathsMapEllipse.end(); ++it) {
208 string name = "fHPoints_Ellipse_" + it->second;
209 fHM->Create2<TH2D>(name, name + ";X_Axis [];Y_Axis [];Entries", 2001, -100., 100., 2001, 60., 210.);
210 }
211
212 fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrack",
213 "fhDistanceCenterToExtrapolatedTrack;Distance fitted "
214 "center to extrapolated track;Number of entries",
215 750, 0., 20.);
216 // fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrackInPlane", "fhDistanceCenterToExtrapolatedTrack;Distance fitted center to extrapolated track;Number of entries", 400, 0., 50.);
217 fHM->Create1<TH1D>("fhDistanceCenterToExtrapolatedTrackInPlane",
218 "fhDistanceCenterToExtrapolatedTrackInPlane;Distance fitted center to "
219 "extrapolated track plane;Number of entries",
220 750, 0., 10.);
221 fHM->Create1<TH1D>("fhDifferenceX",
222 "fhDifferenceX;Difference in X (fitted center - "
223 "extrapolated track);Number of entries",
224 750, 0., 10.);
225 fHM->Create1<TH1D>("fhDifferenceY",
226 "fhDifferenceY;Difference in Y (fitted center - "
227 "extrapolated track);Number of entries",
228 750, 0., 10.);
229}
230
231void CbmRichPMTMapping::Exec(Option_t* /*option*/)
232{
233 cout << endl
234 << "//"
235 "--------------------------------------------------------------------"
236 "------------------------------------------//"
237 << endl;
238 cout << "//---------------------------------------- EXEC Function - Defining "
239 "the entries ----------------------------------------//"
240 << endl;
241 cout << "//"
242 "--------------------------------------------------------------------"
243 "--------------------------------------------------//"
244 << endl;
245 fEventNum++;
246 //LOG(debug2) << "CbmRichPMTMapping : Event #" << fEventNum;
247 cout << "CbmRichPMTMapping : Event #" << fEventNum << endl;
248
249 Int_t nofRingsInEvent = fRichRings->GetEntriesFast();
250 Int_t nofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
251 Int_t nofHitsInEvent = fRichHits->GetEntriesFast();
252 Int_t NofMCPoints = fRichMCPoints->GetEntriesFast();
253 Int_t NofMCTracks = fMCTracks->GetEntriesFast();
254 cout << "Nb of rings in evt = " << nofRingsInEvent << ", nb of mirror points = " << nofMirrorPoints
255 << ", nb of hits in evt = " << nofHitsInEvent << ", nb of Monte-Carlo points = " << NofMCPoints
256 << " and nb of Monte-Carlo tracks = " << NofMCTracks << endl
257 << endl;
258
259 cout << "//"
260 "--------------------------------------------------------------------"
261 "--------------------------------------//"
262 << endl;
263 cout << "//---------------------------------------- EXEC Function "
264 "---------------------------------------------------//"
265 << endl;
266 cout << "//------------------------------ Mapping for mirror - PMT relations "
267 "----------------------------------------//"
268 << endl;
269 cout << "//"
270 "--------------------------------------------------------------------"
271 "--------------------------------------//"
272 << endl
273 << endl;
274
275 if (nofRingsInEvent == 0) {
276 cout << "Error no rings registered in event." << endl << endl;
277 }
278 else {
279 //MatchFinder();
280 fGP = CbmRichHitProducer::Init();
281 fGP.Print();
282 //ProjectionProducer();
283 ProjectionProducer2();
284 }
285}
286
287void CbmRichPMTMapping::MatchFinder()
288{
289 cout << "//---------------------------------------- MATCH_FINDER Function "
290 "----------------------------------------//"
291 << endl;
292 Int_t NofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
293 Int_t NofRingsInEvent = fRichRings->GetEntriesFast();
294 Int_t NofMCPoints = fRichMCPoints->GetEntriesFast();
295 Int_t NofMCTracks = fMCTracks->GetEntriesFast();
296 //Int_t NofProjections = fRichProjections->GetEntriesFast();
297
298 Double_t x_Mirr = 0, y_Mirr = 0, z_Mirr = 0, x_PMT = 0, y_PMT = 0, z_PMT = 0;
299 Double_t CenterX = 0, CenterY = 0;
300 TGeoNode* mirr_node;
301 const Char_t* mirr_path;
302 UShort_t pHit;
303
304 for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
305 CbmRichPoint* MirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
306 Int_t trackID = MirrPoint->GetTrackID();
307 //cout << "Track ID = " << trackID << endl;
308
309 for (Int_t iMCPoint = 0; iMCPoint < NofMCPoints; iMCPoint++) {
310 CbmRichPoint* pPoint = (CbmRichPoint*) fRichMCPoints->At(iMCPoint);
311 if (NULL == pPoint) continue;
312 CbmMCTrack* pTrack = (CbmMCTrack*) fMCTracks->At(pPoint->GetTrackID());
313 if (NULL == pTrack) continue;
314 Int_t gcode = pTrack->GetPdgCode();
315 Int_t motherID = pTrack->GetMotherId();
316 if (motherID == -1) continue;
317
318 if (trackID == motherID) {
319 //cout << "MATCH BETWEEN TRACK ID AND MOTHER ID FOUND !" << endl;
320 //sleep(2);
321 //cout << "TrackID from mirror point = " << trackID << " and mother ID from MC point = " << motherID << endl;
322 x_Mirr = MirrPoint->GetX();
323 y_Mirr = MirrPoint->GetY();
324 z_Mirr = MirrPoint->GetZ();
325 //cout << "x_Mirr: " << x_Mirr << ", y_Mirr: " << y_Mirr << " and z_Mirr: " << z_Mirr << endl;
326 mirr_node = gGeoManager->FindNode(x_Mirr, y_Mirr, z_Mirr);
327 mirr_path = gGeoManager->GetPath();
328 FillPMTMap(mirr_path, pPoint);
329 //break;
330 }
331 }
332
333 // Loop filling the PMT map with ellipse points
334 for (Int_t iRing = 0; iRing < NofRingsInEvent; iRing++) {
335 CbmRichRing* ring = (CbmRichRing*) fRichRings->At(iRing);
336 if (NULL == ring) continue;
337 CbmRichRingLight ringL;
338 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
339 //fCopFit->DoFit(&ringL);
340 fTauFit->DoFit(&ringL);
341 CenterX = ringL.GetCenterX();
342 CenterY = ringL.GetCenterY();
343 CbmTrackMatchNew* richRingMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iRing));
344 Int_t richMCID = richRingMatch->GetMatchedLink().GetIndex();
345 //Int_t trackID2 = ring->GetTrackID();
346 //cout << "Track ID from ring = " << richMCID << endl;
347 if (richMCID == -1) continue;
348
349 if (trackID == richMCID) {
350
351 cout << "MATCH BETWEEN TRACK_ID AND RICH_MC_ID FOUND !" << endl;
352 cout << "Center X = " << CenterX << " and center Y = " << CenterY << endl;
353 x_Mirr = MirrPoint->GetX();
354 y_Mirr = MirrPoint->GetY();
355 z_Mirr = MirrPoint->GetZ();
356 cout << "x_Mirr: " << x_Mirr << ", y_Mirr: " << y_Mirr << " and z_Mirr: " << z_Mirr << endl;
357 mirr_node = gGeoManager->FindNode(x_Mirr, y_Mirr, z_Mirr);
358 mirr_path = gGeoManager->GetPath();
359 cout << "Mirror path: " << mirr_path << endl;
360
361 FillPMTMapEllipse(mirr_path, ringL.GetCenterX(), ringL.GetCenterY());
362 //break;
363 }
364 }
365 }
366}
367
368void CbmRichPMTMapping::FillPMTMap(const Char_t* mirr_path, CbmRichPoint* pPoint)
369{
370 //cout << "//---------------------------------------- FILL_PMT_MAP Function ----------------------------------------//" << endl;
371 string name = string(mirr_path);
372 for (std::map<string, string>::iterator it = fPathsMap.begin(); it != fPathsMap.end(); ++it) {
373 if (name.compare(it->first) == 0) {
374 //cout << "IDENTICAL PATHS FOUND !" << endl;
375 //cout << "Mirror ID: " << it->second << endl;
376 Double_t x_PMT = pPoint->GetX();
377 Double_t y_PMT = pPoint->GetY();
378 Double_t z_PMT = pPoint->GetZ();
379 //cout << "x_PMT: " << x_PMT << ", y_PMT: " << y_PMT << " and z_PMT: " << z_PMT << endl << endl;
380 //sleep(1);
381 fHM->H2("fHMCPoints_" + it->second)->Fill(-x_PMT, y_PMT);
382 fMirrCounter++;
383 }
384 }
385}
386
387void CbmRichPMTMapping::FillPMTMapEllipse(const Char_t* mirr_path, Float_t CenterX, Float_t CenterY)
388{
389 cout << "//---------------------------------------- FILL_PMT_MAP_ELLIPSE "
390 "Function ----------------------------------------//"
391 << endl;
392 string name = string(mirr_path);
393 for (std::map<string, string>::iterator it = fPathsMap.begin(); it != fPathsMap.end(); ++it) {
394 if (name.compare(it->first) == 0) {
395 cout << "IDENTICAL PATHS FOUND !" << endl;
396 //sleep(2);
397 cout << "Mirror ID: " << it->second << endl;
398 //cout << "Center X: " << CenterX << " and Center Y: " << CenterY << endl << endl;
399 fHM->H2("fHPoints_Ellipse_" + it->second)->Fill(-CenterX, CenterY);
400 }
401 }
402 cout << endl;
403}
404
405void CbmRichPMTMapping::ProjectionProducer2()
406{
407 cout << "//------------------------------ CbmRichPMTMapping: Projection "
408 "Producer ------------------------------//"
409 << endl
410 << endl;
411
412 Int_t NofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
413 Int_t NofRingsInEvent = fRichRings->GetEntriesFast();
414 Int_t NofGTracks = fGlobalTracks->GetEntriesFast();
415 Int_t NofRefPlanePoints = fRichRefPlanePoints->GetEntriesFast();
416 Int_t NofPMTPoints = fRichPoints->GetEntriesFast();
417
418 // Declarations of intermediate calculated variables.
419 Double_t t1 = 0., t2 = 0., t3 = 0., k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
420 // Declaration of points coordinates.
421 Double_t sphereRadius = 0., constantePMT = 0.;
422 Double_t ptMirr[] = {0., 0., 0.}, ptC[] = {0., 0., 0.}, ptR1[] = {0., 0., 0.}, normalPMT[] = {0., 0., 0.},
423 normalMirr[] = {0., 0., 0.};
424 Double_t ptR2Mirr[] = {0., 0., 0.}, ptR2Center[] = {0., 0., 0.}, ptPMirr[] = {0., 0., 0.}, ptPR2[] = {0., 0., 0.};
425 Double_t reflectedPtCooVectSphereUnity[] = {0., 0., 0.};
426 // Declaration of ring parameters.
427 Double_t ringCenter[] = {0., 0., 0.}, distToExtrapTrackHit = 0., distToExtrapTrackHitInPlane = 0.;
428 //Declarations related to geometry.
429 Int_t mirrTrackID = -1, pmtTrackID = -1, refPlaneTrackID = -1, motherID = -100, pmtMotherID = -100;
430 const Char_t *mirrPath, *topNodePath;
431 CbmMCTrack *track = NULL, *track2 = NULL;
432 TGeoNode* mirrNode;
433 TGeoMatrix *mirrMatrix, *pmtMatrix, *richMatrix;
434 TGeoShape* ptrShape;
435
436 GetPmtNormal(NofPMTPoints, normalPMT[0], normalPMT[1], normalPMT[2], constantePMT);
437 cout << "Calculated normal vector to PMT plane = {" << normalPMT[0] << ", " << normalPMT[1] << ", " << normalPMT[2]
438 << "} and constante d = " << constantePMT << endl
439 << endl;
440
441 for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
442 //cout << "NofMirrorPoints = " << NofMirrorPoints << " and iMirr = " << iMirr << endl;
443 CbmRichPoint* mirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
444 mirrTrackID = mirrPoint->GetTrackID();
445 //cout << "Mirror track ID = " << mirrTrackID << endl;
446 if (mirrTrackID <= -1) {
447 cout << "Mirror track ID <= 1 !!!" << endl;
448 cout << "----------------------------------- End of loop N°" << iMirr + 1
449 << " on the mirror points. -----------------------------------" << endl
450 << endl;
451 continue;
452 }
453 track = (CbmMCTrack*) fMCTracks->At(mirrTrackID);
454 motherID = track->GetMotherId();
455 if (motherID == -1) {
456 //cout << "Mirror motherID == -1 !!!" << endl << endl;
457 ptMirr[0] = mirrPoint->GetX(), ptMirr[1] = mirrPoint->GetY(), ptMirr[2] = mirrPoint->GetZ();
458 //cout << "Mirror Point coordinates; x = " << ptMirr[0] << ", y = " << ptMirr[1] << " and z = " << ptMirr[2] << endl;
459 mirrNode = gGeoManager->FindNode(ptMirr[0], ptMirr[1], ptMirr[2]);
460 //mirrPath = gGeoManager->GetPath();
461 mirrPath = mirrNode->GetName();
462 topNodePath = gGeoManager->GetTopNode()->GetName();
463 cout << "Top node path: " << topNodePath << " and mirror path: " << mirrPath << endl;
464 mirrMatrix = mirrNode->GetMatrix();
465 cout << "Mirror matrix parameters: " << endl;
466 mirrMatrix->Print();
467 ptrShape = mirrNode->GetVolume()->GetShape();
468 cout << "Shape of the mirror tile:" << endl;
469 ptrShape->Dump();
470
471 if (ptMirr[1] > 0) {
472 fIsMirrorUpperHalf = true;
473 }
474 else {
475 fIsMirrorUpperHalf = false;
476 }
477 CalculateSphereParameters2(mirrPath, ptC[0], ptC[1], ptC[2], sphereRadius);
478 cout << endl
479 << "Sphere center coordinates of the rotated mirror tile = {" << ptC[0] << ", " << ptC[1] << ", " << ptC[2]
480 << "} and sphere inner radius = " << sphereRadius << endl;
481
482 for (Int_t iRefl = 0; iRefl < NofRefPlanePoints; iRefl++) {
483 CbmRichPoint* refPlanePoint = (CbmRichPoint*) fRichRefPlanePoints->At(iRefl);
484 refPlaneTrackID = refPlanePoint->GetTrackID();
485 //cout << "Reflective plane track ID = " << refPlaneTrackID << endl;
486 if (mirrTrackID == refPlaneTrackID) {
487 //cout << "IDENTICAL TRACK ID FOUND !!!" << endl << endl;
488 ptR1[0] = refPlanePoint->GetX(), ptR1[1] = refPlanePoint->GetY(), ptR1[2] = refPlanePoint->GetZ();
489 cout << "Reflective Plane Point coordinates = {" << ptR1[0] << ", " << ptR1[1] << ", " << ptR1[2] << "}"
490 << endl;
491 cout << "Mirror Point coordinates = {" << ptMirr[0] << ", " << ptMirr[1] << ", " << ptMirr[2] << "}" << endl
492 << endl;
493 normalMirr[0] = (ptC[0] - ptMirr[0])
494 / TMath::Sqrt(TMath::Power(ptC[0] - ptMirr[0], 2) + TMath::Power(ptC[1] - ptMirr[1], 2)
495 + TMath::Power(ptC[2] - ptMirr[2], 2));
496 normalMirr[1] = (ptC[1] - ptMirr[1])
497 / TMath::Sqrt(TMath::Power(ptC[0] - ptMirr[0], 2) + TMath::Power(ptC[1] - ptMirr[1], 2)
498 + TMath::Power(ptC[2] - ptMirr[2], 2));
499 normalMirr[2] = (ptC[2] - ptMirr[2])
500 / TMath::Sqrt(TMath::Power(ptC[0] - ptMirr[0], 2) + TMath::Power(ptC[1] - ptMirr[1], 2)
501 + TMath::Power(ptC[2] - ptMirr[2], 2));
502 cout << "Calculated and normalized normal of mirror tile = {" << normalMirr[0] << ", " << normalMirr[1]
503 << ", " << normalMirr[2] << "}" << endl;
504
505 t1 = ((ptR1[0] - ptMirr[0]) * (ptC[0] - ptMirr[0]) + (ptR1[1] - ptMirr[1]) * (ptC[1] - ptMirr[1])
506 + (ptR1[2] - ptMirr[2]) * (ptC[2] - ptMirr[2]))
507 / (TMath::Power(ptC[0] - ptMirr[0], 2) + TMath::Power(ptC[1] - ptMirr[1], 2)
508 + TMath::Power(ptC[2] - ptMirr[2], 2));
509 ptR2Center[0] = 2 * (ptMirr[0] + t1 * (ptC[0] - ptMirr[0])) - ptR1[0];
510 ptR2Center[1] = 2 * (ptMirr[1] + t1 * (ptC[1] - ptMirr[1])) - ptR1[1];
511 ptR2Center[2] = 2 * (ptMirr[2] + t1 * (ptC[2] - ptMirr[2])) - ptR1[2];
512 t2 = ((ptR1[0] - ptC[0]) * (ptC[0] - ptMirr[0]) + (ptR1[1] - ptC[1]) * (ptC[1] - ptMirr[1])
513 + (ptR1[2] - ptC[2]) * (ptC[2] - ptMirr[2]))
514 / (TMath::Power(ptC[0] - ptMirr[0], 2) + TMath::Power(ptC[1] - ptMirr[1], 2)
515 + TMath::Power(ptC[2] - ptMirr[2], 2));
516 ptR2Mirr[0] = 2 * (ptC[0] + t2 * (ptC[0] - ptMirr[0])) - ptR1[0];
517 ptR2Mirr[1] = 2 * (ptC[1] + t2 * (ptC[1] - ptMirr[1])) - ptR1[1];
518 ptR2Mirr[2] = 2 * (ptC[2] + t2 * (ptC[2] - ptMirr[2])) - ptR1[2];
519 /*//SAME AS calculation of t2 above
520 t3 = ((ptR1[0]-ptC[0])*(ptC[0]-ptMirr[0]) + (ptR1[1]-ptC[1])*(ptC[1]-ptMirr[1]) + (ptR1[2]-ptC[2])*(ptC[2]-ptMirr[2]))/TMath::Sqrt(TMath::Power(ptC[0] - ptMirr[0],2)+TMath::Power(ptC[1] - ptMirr[1],2)+TMath::Power(ptC[2] - ptMirr[2],2));
521 reflectedPtCooVectSphereUnity[0] = 2*(ptC[0]+t3*(normalMirr[0]))-ptR1[0];
522 reflectedPtCooVectSphereUnity[1] = 2*(ptC[1]+t3*(normalMirr[1]))-ptR1[1];
523 reflectedPtCooVectSphereUnity[2] = 2*(ptC[2]+t3*(normalMirr[2]))-ptR1[2];*/
524 cout << "Coordinates of point R2 on reflective plane after "
525 "reflection on the mirror tile:"
526 << endl;
527 cout << "* using mirror point M to define \U00000394: {" << ptR2Center[0] << ", " << ptR2Center[1] << ", "
528 << ptR2Center[2] << "}" << endl;
529 cout << "* using sphere center C to define \U00000394: {" << ptR2Mirr[0] << ", " << ptR2Mirr[1] << ", "
530 << ptR2Mirr[2] << "}" << endl
531 << endl;
532 //cout << "Ref Pt Coo using unity Mirror-Sphere vector & sphere pt = {" << reflectedPtCooVectSphereUnity[0] << ", " << reflectedPtCooVectSphereUnity[1] << ", " << reflectedPtCooVectSphereUnity[2] << "}" << endl << endl;
533 //cout << "NofPMTPoints = " << NofPMTPoints << endl;
534
535 k1 = -1
536 * ((normalPMT[0] * ptMirr[0] + normalPMT[1] * ptMirr[1] + normalPMT[2] * ptMirr[2] + constantePMT)
537 / (normalPMT[0] * (ptR2Mirr[0] - ptMirr[0]) + normalPMT[1] * (ptR2Mirr[1] - ptMirr[1])
538 + normalPMT[2] * (ptR2Mirr[2] - ptMirr[2])));
539 ptPMirr[0] = ptMirr[0] + k1 * (ptR2Mirr[0] - ptMirr[0]);
540 ptPMirr[1] = ptMirr[1] + k1 * (ptR2Mirr[1] - ptMirr[1]);
541 ptPMirr[2] = ptMirr[2] + k1 * (ptR2Mirr[2] - ptMirr[2]);
542 k2 = -1
543 * ((normalPMT[0] * ptR2Mirr[0] + normalPMT[1] * ptR2Mirr[1] + normalPMT[2] * ptR2Mirr[2] + constantePMT)
544 / (normalPMT[0] * (ptR2Mirr[0] - ptMirr[0]) + normalPMT[1] * (ptR2Mirr[1] - ptMirr[1])
545 + normalPMT[2] * (ptR2Mirr[2] - ptMirr[2])));
546 ptPR2[0] = ptR2Mirr[0] + k2 * (ptR2Mirr[0] - ptMirr[0]);
547 ptPR2[1] = ptR2Mirr[1] + k2 * (ptR2Mirr[1] - ptMirr[1]);
548 ptPR2[2] = ptR2Mirr[2] + k2 * (ptR2Mirr[2] - ptMirr[2]);
549 cout << "Coordinates of point P on PMT plane, after reflection on "
550 "the mirror tile and extrapolation to the PMT plane:"
551 << endl;
552 cout << "* using mirror point M to define \U0001D49F ': {" << ptPMirr[0] << ", " << ptPMirr[1] << ", "
553 << ptPMirr[2] << "}" << endl;
554 cout << "* using reflected point R2 to define \U0001D49F ': {" << ptPR2[0] << ", " << ptPR2[1] << ", "
555 << ptPR2[2] << "}" << endl;
556 checkCalc1 = ptPMirr[0] * normalPMT[0] + ptPMirr[1] * normalPMT[1] + ptPMirr[2] * normalPMT[2] + constantePMT;
557 cout << "Check whether extrapolated track point on PMT plane "
558 "verifies its equation (value should be 0.):"
559 << endl;
560 cout << "* using mirror point M, checkCalc = " << checkCalc1 << endl;
561 checkCalc2 = ptPR2[0] * normalPMT[0] + ptPR2[1] * normalPMT[1] + ptPR2[2] * normalPMT[2] + constantePMT;
562 cout << "* using reflected point R2, checkCalc = " << checkCalc2 << endl;
563
564 TVector3 pmtVector(ptPMirr[0], ptPMirr[1], ptPMirr[2]);
565 TVector3 pmtVectorNew;
566 CbmRichHitProducer::TiltPoint(&pmtVector, &pmtVectorNew, fGP.fPmt.fPhi, fGP.fPmt.fTheta, fGP.fPmtZOrig);
567 cout << "New coordinates of point P on PMT plane, after PMT plane "
568 "rotation = {"
569 << pmtVectorNew.X() << ", " << pmtVectorNew.Y() << ", " << pmtVectorNew.Z() << "}" << endl
570 << endl;
571 ptPMirr[0] = pmtVectorNew.X(), ptPMirr[1] = pmtVectorNew.Y(), ptPMirr[2] = pmtVectorNew.Z();
572
573 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
574 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
575 pmtTrackID = pmtPoint->GetTrackID();
576 CbmMCTrack* track3 = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
577 pmtMotherID = track3->GetMotherId();
578 //cout << "pmt mother ID = " << pmtMotherID << endl;
579 if (pmtMotherID == mirrTrackID) {
580 ptP[0] = pmtPoint->GetX(), ptP[1] = pmtPoint->GetY(), ptP[2] = pmtPoint->GetZ();
581 //cout << "Identical mirror track ID and PMT mother ID !!!" << endl;
582 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
583 }
584 }
585 cout << "Looking for PMT hits: end." << endl << endl;*/
586 }
587 //else { cout << "No identical track ID between mirror point and reflective plane point found ..." << endl << endl; }
588 }
589
590 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
591 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
592 cout << "PMT points = {" << pmtPoint->GetX() << ", " << pmtPoint->GetY() << ", " << pmtPoint->GetZ() << "}" << endl;
593 TVector3 inputPoint(pmtPoint->GetX(), pmtPoint->GetY(), pmtPoint->GetZ());
594 TVector3 outputPoint;
595 CbmRichHitProducer::TiltPoint(&inputPoint, &outputPoint, fGP.fPmtPhi, fGP.fPmtTheta, fGP.fPmtZOrig);
596 cout << "New PMT points after rotation of the PMT plane: " << endl;
597 cout << outputPoint.X() << "\t" << outputPoint.Y() << "\t" << outputPoint.Z() << endl;
598 }*/
599
600 for (Int_t iGlobalTrack = 0; iGlobalTrack < NofGTracks; iGlobalTrack++) {
601 //cout << "Nb of global tracks = " << NofGTracks << " and iGlobalTrack = " << iGlobalTrack << endl;
602 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack);
603 if (NULL == gTrack) continue;
604 Int_t richInd = gTrack->GetRichRingIndex();
605 //cout << "Rich index = " << richInd << endl;
606 if (richInd < 0) {
607 continue;
608 }
609 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richInd));
610 if (NULL == ring) {
611 continue;
612 }
613 Int_t ringTrackID = ring->GetTrackID();
614 track2 = (CbmMCTrack*) fMCTracks->At(ringTrackID);
615 Int_t ringMotherID = track2->GetMotherId();
616 //cout << "Ring mother ID = " << ringMotherID << ", ring track ID = " << ringTrackID << " and track2 pdg = " << track2->GetPdgCode() << endl;
617 CbmRichRingLight ringL;
618 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
619 //RotateAndCopyHitsToRingLight(ring, &ringL);
620 fCopFit->DoFit(&ringL);
621 fTauFit->DoFit(&ringL);
622 ringCenter[0] = ringL.GetCenterX();
623 ringCenter[1] = ringL.GetCenterY();
624 ringCenter[2] =
625 -1 * ((normalPMT[0] * ringCenter[0] + normalPMT[1] * ringCenter[1] + constantePMT) / normalPMT[2]);
626 cout << "Ring center coordinates = {" << ringCenter[0] << ", " << ringCenter[1] << ", " << ringCenter[2] << "}"
627 << endl;
628 cout << "Difference in X = " << TMath::Abs(ringCenter[0] - ptPMirr[0]) << "\t"
629 << "Difference in Y = " << TMath::Abs(ringCenter[1] - ptPMirr[1]) << "\t"
630 << "Difference in Z = " << TMath::Abs(ringCenter[2] - ptPMirr[2]) << endl;
631 fHM->H1("fhDifferenceX")->Fill(TMath::Abs(ringCenter[0] - ptPMirr[0]));
632 fHM->H1("fhDifferenceY")->Fill(TMath::Abs(ringCenter[1] - ptPMirr[1]));
633 distToExtrapTrackHit =
634 TMath::Sqrt(TMath::Power(ringCenter[0] - ptPMirr[0], 2) + TMath::Power(ringCenter[1] - ptPMirr[1], 2)
635 + TMath::Power(ringCenter[2] - ptPMirr[2], 2));
636 distToExtrapTrackHitInPlane =
637 TMath::Sqrt(TMath::Power(ringCenter[0] - ptPMirr[0], 2) + TMath::Power(ringCenter[1] - ptPMirr[1], 2));
638 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
639 fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlane")->Fill(distToExtrapTrackHitInPlane);
640 cout << "Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl;
641 cout << "Distance between fitted ring center and extrapolated track "
642 "hit in plane = "
643 << distToExtrapTrackHitInPlane << endl
644 << endl;
645 //}
646 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
647 }
648 cout << "End of loop on global tracks;" << endl;
649 }
650 else {
651 cout << "Not a mother particle ..." << endl;
652 }
653 cout << "----------------------------------- "
654 << "End of loop N°" << iMirr + 1 << " on the mirror points."
655 << " -----------------------------------" << endl
656 << endl;
657 }
658}
659
660void CbmRichPMTMapping::ProjectionProducer()
661{
662 cout << "//------------------------------ CbmRichPMTMapping: Projection "
663 "Producer ------------------------------//"
664 << endl
665 << endl;
666
667 Int_t NofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
668 Int_t NofRingsInEvent = fRichRings->GetEntriesFast();
669 Int_t NofGTracks = fGlobalTracks->GetEntriesFast();
670 Int_t NofRefPlanePoints = fRichRefPlanePoints->GetEntriesFast();
671 Int_t NofPMTPoints = fRichPoints->GetEntriesFast();
672
673 //Declarations of intermediate and calculated variables.
674 Double_t sphereX = 0., sphereY = 0., sphereZ = 0., sphereR = 0., normalX = 0., normalY = 0., normalZ = 0.,
675 normalCste = 0.;
676 Double_t CenterX = 0., CenterY = 0., CenterZ = 0., distToExtrapTrackHit = 0.;
677 Double_t a1 = 0., a2 = 0., a3 = 0., a4 = 0., a5 = 0., t1 = 0., t2 = 0.;
678 // Declaration of points coordinates.
679 Double_t refPointCoo[] = {0., 0., 0.}, refPointMom[] = {0., 0., 0.};
680 Double_t reflectedPtCooNormMirr[] = {0., 0., 0.}, reflectedPtCooNormSphere[] = {0., 0., 0.};
681 Double_t reflectedPtCooVectMirr[] = {0., 0., 0.}, reflectedPtCooVectSphere[] = {0., 0., 0.};
682 Double_t reflectedPtCooVectSphereUnity[] = {0., 0., 0.}, vectMSUnity[] = {0., 0., 0.};
683 Double_t mirrPt[] = {0., 0., 0.}, mirrMom[] = {0., 0., 0.}, pmtPt[] = {0., 0., 0.};
684 Double_t computedNormal[] = {0., 0., 0.}, computedNormal2[] = {0., 0., 0.};
685 Double_t extrapolatedTrackCoo[] = {0., 0., 0.}, extrapolatedTrackCooComputedNormal[] = {0., 0., 0.};
686 Double_t checkCalc = 0.;
687 //Declarations related to geometries.
688 Int_t mirrTrackID = -1, pmtTrackID = -1, refPlaneTrackID = -1, motherID = -100, pmtMotherID = -100;
689 const Char_t *mirrPath, *topNodePath;
690 CbmMCTrack *track = NULL, *track2 = NULL;
691 TGeoNode* mirrNode;
692 TGeoMatrix *mirrMatrix, *pmtMatrix, *richMatrix;
693 TGeoShape* ptrShape;
694
695 /*From the point you can get the information using the following methods:
696 Double_t GetPx()
697 Double_t GetPy()
698 Double_t GetPz()
699 or if you want TVector - > void Momentum(TVector3& mom)
700 Double_t GetX()
701 Double_t GetY()
702 Double_t GetZ()
703 or if you want TVector void Position(TVector3& pos)*/
704
705 GetPmtNormal(NofPMTPoints, normalX, normalY, normalZ, normalCste);
706 cout << "Calculated normal vector to PMT plane = {" << normalX << ", " << normalY << ", " << normalZ
707 << "} and constante d = " << normalCste << endl
708 << endl;
709
710 for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
711 //cout << "NofMirrorPoints = " << NofMirrorPoints << " and iMirr = " << iMirr << endl;
712 CbmRichPoint* mirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
713 mirrTrackID = mirrPoint->GetTrackID();
714 //cout << "Mirror track ID = " << mirrTrackID << endl;
715 if (mirrTrackID <= -1) {
716 cout << "Mirror track ID <= 1 !!!" << endl;
717 cout << "----------------------------------- End of loop N°" << iMirr + 1
718 << " on the mirror points. -----------------------------------" << endl
719 << endl;
720 continue;
721 }
722 track = (CbmMCTrack*) fMCTracks->At(mirrTrackID);
723 motherID = track->GetMotherId();
724 if (motherID == -1) {
725 //cout << "Mirror motherID == -1 !!!" << endl << endl;
726 mirrPt[0] = mirrPoint->GetX(), mirrPt[1] = mirrPoint->GetY(), mirrPt[2] = mirrPoint->GetZ();
727 mirrMom[0] = mirrPoint->GetPx(), mirrMom[1] = mirrPoint->GetPy(), mirrMom[2] = mirrPoint->GetPz();
728 //cout << "Mirror Point coordinates; x = " << mirrPt[0] << ", y = " << mirrPt[1] << " and z = " << mirrPt[2] << endl;
729 mirrNode = gGeoManager->FindNode(mirrPt[0], mirrPt[1], mirrPt[2]);
730 //mirrPath = gGeoManager->GetPath();
731 mirrPath = mirrNode->GetName();
732 topNodePath = gGeoManager->GetTopNode()->GetName();
733 cout << "Mirror path: " << mirrPath << " and top node path: " << topNodePath << endl;
734 mirrMatrix = mirrNode->GetMatrix();
735 cout << "Mirror matrix parameters: " << endl;
736 mirrMatrix->Print();
737 ptrShape = mirrNode->GetVolume()->GetShape();
738 ptrShape->Dump();
739
740 if (mirrPt[1] > 0) {
741 fIsMirrorUpperHalf = true;
742 }
743 else {
744 fIsMirrorUpperHalf = false;
745 }
746 Double_t sphere2X = 0., sphere2Y = 0., sphere2Z = 0., sphere2R = 0.;
747 CalculateSphereParameters(mirrPath, sphere2X, sphere2Y, sphere2Z, sphere2R);
748 cout << endl
749 << "Old sphere coordinates = {" << sphere2X << ", " << sphere2Y << ", " << sphere2Z
750 << "} and sphere inner radius = " << sphere2R << endl;
751 CalculateSphereParameters2(mirrPath, sphereX, sphereY, sphereZ, sphereR);
752 cout << "New sphere coordinates = {" << sphereX << ", " << sphereY << ", " << sphereZ
753 << "} and sphere inner radius = " << sphereR << endl;
754
755 for (Int_t iRefl = 0; iRefl < NofRefPlanePoints; iRefl++) {
756 CbmRichPoint* refPlanePoint = (CbmRichPoint*) fRichRefPlanePoints->At(iRefl);
757 refPlaneTrackID = refPlanePoint->GetTrackID();
758 //cout << "Reflective plane track ID = " << refPlaneTrackID << endl;
759 if (mirrTrackID == refPlaneTrackID) {
760 //cout << "IDENTICAL TRACK ID FOUND !!!" << endl << endl;
761 refPointCoo[0] = refPlanePoint->GetX(), refPointCoo[1] = refPlanePoint->GetY(),
762 refPointCoo[2] = refPlanePoint->GetZ();
763 refPointMom[0] = refPlanePoint->GetPx(), refPointMom[1] = refPlanePoint->GetPy(),
764 refPointMom[2] = refPlanePoint->GetPz();
765 cout << "Reflective Plane Point coordinates = {" << refPointCoo[0] << ", " << refPointCoo[1] << ", "
766 << refPointCoo[2] << "} and momentum = {" << refPointMom[0] << ", " << refPointMom[1] << ", "
767 << refPointMom[2] << "}" << endl;
768 cout << "Mirror Point coordinates = {" << mirrPt[0] << ", " << mirrPt[1] << ", " << mirrPt[2]
769 << "} and momentum = {" << mirrMom[0] << ", " << mirrMom[1] << ", " << mirrMom[2] << "}" << endl
770 << endl;
771 ptrShape->ComputeNormal(refPointCoo, refPointMom, computedNormal);
772 cout << "Computed normal to mirror tile coordinates = {" << computedNormal[0] << ", " << computedNormal[1]
773 << ", " << computedNormal[2] << "}" << endl;
774 /*ptrShape->ComputeNormal(mirrPt, mirrMom, computedNormal2);
775 cout << "Computed normal 2 to mirror tile coordinates = {" << computedNormal2[0] << ", " << computedNormal2[1] << ", " << computedNormal2[2] << "}" << endl;*/
776 vectMSUnity[0] = (sphereX - mirrPt[0])
777 / TMath::Sqrt(TMath::Power(sphereX - mirrPt[0], 2) + TMath::Power(sphereY - mirrPt[1], 2)
778 + TMath::Power(sphereZ - mirrPt[2], 2));
779 vectMSUnity[1] = (sphereY - mirrPt[1])
780 / TMath::Sqrt(TMath::Power(sphereX - mirrPt[0], 2) + TMath::Power(sphereY - mirrPt[1], 2)
781 + TMath::Power(sphereZ - mirrPt[2], 2));
782 vectMSUnity[2] = (sphereZ - mirrPt[2])
783 / TMath::Sqrt(TMath::Power(sphereX - mirrPt[0], 2) + TMath::Power(sphereY - mirrPt[1], 2)
784 + TMath::Power(sphereZ - mirrPt[2], 2));
785 cout << "Calculated unity Mirror-Sphere vector = {" << vectMSUnity[0] << ", " << vectMSUnity[1] << ", "
786 << vectMSUnity[2] << "}" << endl;
787
788 a1 = (computedNormal[0] * (refPointCoo[0] - mirrPt[0]) + computedNormal[1] * (refPointCoo[1] - mirrPt[1])
789 + computedNormal[2] * (refPointCoo[2] - mirrPt[2]))
790 / (TMath::Power(computedNormal[0], 2) + TMath::Power(computedNormal[1], 2)
791 + TMath::Power(computedNormal[2], 2));
792 reflectedPtCooNormMirr[0] = 2 * (mirrPt[0] + a1 * computedNormal[0]) - refPointCoo[0];
793 reflectedPtCooNormMirr[1] = 2 * (mirrPt[1] + a1 * computedNormal[1]) - refPointCoo[1];
794 reflectedPtCooNormMirr[2] = 2 * (mirrPt[2] + a1 * computedNormal[2]) - refPointCoo[2];
795 /*a2 = (computedNormal[0]*(refPointCoo[0]-sphereX) + computedNormal[1]*(refPointCoo[1]-sphereY) + computedNormal[2]*(refPointCoo[2]-sphereZ))/(TMath::Power(computedNormal[0],2) + TMath::Power(computedNormal[1],2) + TMath::Power(computedNormal[2],2));
796 reflectedPtCooNormSphere[0] = 2*(mirrPt[0]+a2*computedNormal[0])-refPointCoo[0];
797 reflectedPtCooNormSphere[1] = 2*(mirrPt[1]+a2*computedNormal[1])-refPointCoo[1];
798 reflectedPtCooNormSphere[2] = 2*(mirrPt[2]+a2*computedNormal[2])-refPointCoo[2];
799 a3 = ((refPointCoo[0]-mirrPt[0])*(sphereX-mirrPt[0]) + (refPointCoo[1]-mirrPt[1])*(sphereY-mirrPt[1]) + (refPointCoo[2]-mirrPt[2])*(sphereZ-mirrPt[2]))/(TMath::Power(sphereX-mirrPt[0],2) + TMath::Power(sphereY-mirrPt[1],2) + TMath::Power(sphereZ-mirrPt[2],2));
800 reflectedPtCooVectMirr[0] = 2*(sphereX+a3*(sphereX-mirrPt[0]))-refPointCoo[0];
801 reflectedPtCooVectMirr[1] = 2*(sphereY+a3*(sphereY-mirrPt[1]))-refPointCoo[1];
802 reflectedPtCooVectMirr[2] = 2*(sphereZ+a3*(sphereZ-mirrPt[2]))-refPointCoo[2];*/
803 a4 = ((refPointCoo[0] - sphereX) * (sphereX - mirrPt[0]) + (refPointCoo[1] - sphereY) * (sphereY - mirrPt[1])
804 + (refPointCoo[2] - sphereZ) * (sphereZ - mirrPt[2]))
805 / (TMath::Power(sphereX - mirrPt[0], 2) + TMath::Power(sphereY - mirrPt[1], 2)
806 + TMath::Power(sphereZ - mirrPt[2], 2));
807 reflectedPtCooVectSphere[0] = 2 * (sphereX + a4 * (sphereX - mirrPt[0])) - refPointCoo[0];
808 reflectedPtCooVectSphere[1] = 2 * (sphereY + a4 * (sphereY - mirrPt[1])) - refPointCoo[1];
809 reflectedPtCooVectSphere[2] = 2 * (sphereZ + a4 * (sphereZ - mirrPt[2])) - refPointCoo[2];
810 /*//SAME AS calculation of a4 above
811 a5 = ((refPointCoo[0]-sphereX)*(sphereX-mirrPt[0]) + (refPointCoo[1]-sphereY)*(sphereY-mirrPt[1]) + (refPointCoo[2]-sphereZ)*(sphereZ-mirrPt[2]))/TMath::Sqrt(TMath::Power(sphereX - mirrPt[0],2)+TMath::Power(sphereY - mirrPt[1],2)+TMath::Power(sphereZ - mirrPt[2],2));
812 reflectedPtCooVectSphereUnity[0] = 2*(sphereX+a5*(vectMSUnity[0]))-refPointCoo[0];
813 reflectedPtCooVectSphereUnity[1] = 2*(sphereY+a5*(vectMSUnity[1]))-refPointCoo[1];
814 reflectedPtCooVectSphereUnity[2] = 2*(sphereZ+a5*(vectMSUnity[2]))-refPointCoo[2];*/
815
816 cout << "Ref Pt Coo using computed normal & mirror pt = {" << reflectedPtCooNormMirr[0] << ", "
817 << reflectedPtCooNormMirr[1] << ", " << reflectedPtCooNormMirr[2] << "}" << endl;
818 //cout << "Ref Pt Coo using normal & sphere pt = {" << reflectedPtCooNormSphere[0] << ", " << reflectedPtCooNormSphere[1] << ", " << reflectedPtCooNormSphere[2] << "}" << endl;
819 //cout << "Ref Pt Coo using MS vector & mirror pt = {" << reflectedPtCooVectMirr[0] << ", " << reflectedPtCooVectMirr[1] << ", " << reflectedPtCooVectMirr[2] << "}" << endl;
820 cout << "Ref Pt Coo using MS vector & sphere pt = {" << reflectedPtCooVectSphere[0] << ", "
821 << reflectedPtCooVectSphere[1] << ", " << reflectedPtCooVectSphere[2] << "}" << endl
822 << endl;
823 //cout << "Ref Pt Coo using unity Mirror-Sphere vector & sphere pt = {" << reflectedPtCooVectSphereUnity[0] << ", " << reflectedPtCooVectSphereUnity[1] << ", " << reflectedPtCooVectSphereUnity[2] << "}" << endl << endl;
824 //cout << "NofPMTPoints = " << NofPMTPoints << endl;
825
826 t1 = -1
827 * ((normalX * reflectedPtCooVectSphere[0] + normalY * reflectedPtCooVectSphere[1]
828 + normalZ * reflectedPtCooVectSphere[2] + normalCste)
829 / (normalX * (reflectedPtCooVectSphere[0] - mirrPt[0])
830 + normalY * (reflectedPtCooVectSphere[1] - mirrPt[1])
831 + normalZ * (reflectedPtCooVectSphere[2] - mirrPt[2])));
832 extrapolatedTrackCoo[0] = reflectedPtCooVectSphere[0] + t1 * (reflectedPtCooVectSphere[0] - mirrPt[0]);
833 extrapolatedTrackCoo[1] = reflectedPtCooVectSphere[1] + t1 * (reflectedPtCooVectSphere[1] - mirrPt[1]);
834 extrapolatedTrackCoo[2] = reflectedPtCooVectSphere[2] + t1 * (reflectedPtCooVectSphere[2] - mirrPt[2]);
835 cout << "Extrapolated track point on PMT plane using MS vector = {" << extrapolatedTrackCoo[0] << ", "
836 << extrapolatedTrackCoo[1] << ", " << extrapolatedTrackCoo[2] << "}" << endl;
837 checkCalc = extrapolatedTrackCoo[0] * normalX + extrapolatedTrackCoo[1] * normalY
838 + extrapolatedTrackCoo[2] * normalZ + normalCste;
839 cout << "Check whether extrapolated track point on PMT plane "
840 "verifies its equation (extrapolation with MS vector method):"
841 << endl;
842 cout << "Check calculation = " << checkCalc << endl;
843
844 t2 =
845 -1
846 * ((normalX * reflectedPtCooNormMirr[0] + normalY * reflectedPtCooNormMirr[1]
847 + normalZ * reflectedPtCooNormMirr[2] + normalCste)
848 / (normalX * (reflectedPtCooNormMirr[0] - mirrPt[0]) + normalY * (reflectedPtCooNormMirr[1] - mirrPt[1])
849 + normalZ * (reflectedPtCooNormMirr[2] - mirrPt[2])));
850 extrapolatedTrackCooComputedNormal[0] =
851 reflectedPtCooNormMirr[0] + t1 * (reflectedPtCooNormMirr[0] - mirrPt[0]);
852 extrapolatedTrackCooComputedNormal[1] =
853 reflectedPtCooNormMirr[1] + t1 * (reflectedPtCooNormMirr[1] - mirrPt[1]);
854 extrapolatedTrackCooComputedNormal[2] =
855 reflectedPtCooNormMirr[2] + t1 * (reflectedPtCooNormMirr[2] - mirrPt[2]);
856 cout << "Extrapolated track point on PMT plane using computed normal = {"
857 << extrapolatedTrackCooComputedNormal[0] << ", " << extrapolatedTrackCooComputedNormal[1] << ", "
858 << extrapolatedTrackCooComputedNormal[2] << "}" << endl;
859 checkCalc = extrapolatedTrackCooComputedNormal[0] * normalX + extrapolatedTrackCooComputedNormal[1] * normalY
860 + extrapolatedTrackCooComputedNormal[2] * normalZ + normalCste;
861 cout << "Check whether extrapolated track point on PMT plane verifies "
862 "its equation (extrapolation with computed normal method):"
863 << endl;
864 cout << "Check calculation = " << checkCalc << endl << endl;
865
866 /*for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
867 CbmRichPoint *pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
868 pmtTrackID = pmtPoint->GetTrackID();
869 CbmMCTrack* track3 = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
870 pmtMotherID = track3->GetMotherId();
871 //cout << "pmt mother ID = " << pmtMotherID << endl;
872 if (pmtMotherID == mirrTrackID) {
873 pmtPt[0] = pmtPoint->GetX(), pmtPt[1] = pmtPoint->GetY(), pmtPt[2] = pmtPoint->GetZ();
874 //cout << "Identical mirror track ID and PMT mother ID !!!" << endl;
875 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
876 }
877 }
878 cout << "Looking for PMT hits: end." << endl << endl;*/
879 }
880 //else { cout << "No identical track ID between mirror point and reflective plane point found ..." << endl << endl; }
881 }
882 for (Int_t iGlobalTrack = 0; iGlobalTrack < NofGTracks; iGlobalTrack++) {
883 //cout << "Nb of global tracks = " << NofGTracks << " and iGlobalTrack = " << iGlobalTrack << endl;
884 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack);
885 if (NULL == gTrack) continue;
886 Int_t richInd = gTrack->GetRichRingIndex();
887 //cout << "Rich index = " << richInd << endl;
888 if (richInd < 0) {
889 continue;
890 }
891 CbmRichRing* ring = (CbmRichRing*) fRichRings->At(richInd);
892 if (NULL == ring) {
893 continue;
894 }
895 Int_t ringTrackID = ring->GetTrackID();
896 track2 = (CbmMCTrack*) fMCTracks->At(ringTrackID);
897 Int_t ringMotherID = track2->GetMotherId();
898 //cout << "Ring mother ID = " << ringMotherID << ", ring track ID = " << ringTrackID << " and track2 pdg = " << track2->GetPdgCode() << endl;
899
900 CbmRichRingLight ringL;
901 CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
902 //fCopFit->DoFit(&ringL);
903 fTauFit->DoFit(&ringL);
904 CenterX = ringL.GetCenterX();
905 CenterY = ringL.GetCenterY();
906 CenterZ = -1 * ((normalX * CenterX + normalY * CenterY + normalCste) / normalZ);
907 cout << "Ring center coordinates = {" << CenterX << ", " << CenterY << ", " << CenterZ << "}" << endl;
908 cout << "Difference in X = " << TMath::Abs(CenterX - extrapolatedTrackCoo[0]) << endl;
909 cout << "Difference in Y = " << TMath::Abs(CenterY - extrapolatedTrackCoo[1]) << endl;
910 cout << "Difference in Z = " << TMath::Abs(CenterZ - extrapolatedTrackCoo[2]) << endl;
911 fHM->H1("fhDifferenceX")->Fill(TMath::Abs(CenterX - extrapolatedTrackCoo[0]));
912 fHM->H1("fhDifferenceY")->Fill(TMath::Abs(CenterY - extrapolatedTrackCoo[1]));
913 distToExtrapTrackHit = TMath::Sqrt(TMath::Power(CenterX - extrapolatedTrackCoo[0], 2)
914 + TMath::Power(CenterY - extrapolatedTrackCoo[1], 2)
915 + TMath::Power(CenterZ - extrapolatedTrackCoo[2], 2));
916 fHM->H1("fhDistanceCenterToExtrapolatedTrack")->Fill(distToExtrapTrackHit);
917 fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlane")
918 ->Fill(TMath::Sqrt(TMath::Power(CenterX - extrapolatedTrackCoo[0], 2)
919 + TMath::Power(CenterY - extrapolatedTrackCoo[1], 2)));
920 cout << "Distance between fitted ring center and extrapolated track hit = " << distToExtrapTrackHit << endl
921 << endl;
922 //}
923 //else { cout << "No identical ring mother ID and mirror track ID ..." << endl;}
924 }
925 cout << "End of loop on global tracks;" << endl;
926 }
927 else {
928 cout << "Not a mother particle ..." << endl;
929 }
930 cout << "----------------------------------- "
931 << "End of loop N°" << iMirr + 1 << " on the mirror points."
932 << " -----------------------------------" << endl
933 << endl;
934 }
935}
936
937void CbmRichPMTMapping::GetPmtNormal(Int_t NofPMTPoints, Double_t& normalX, Double_t& normalY, Double_t& normalZ,
938 Double_t& normalCste)
939{
940 //cout << endl << "//------------------------------ CbmRichPMTMapping: Calculate PMT Normal ------------------------------//" << endl << endl;
941
942 Int_t pmtTrackID, pmtMotherID;
943 Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0., scalarProd = 0.;
944 Double_t pmtPt[] = {0., 0., 0.};
945 Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
946 CbmMCTrack* track;
947
948 /*
949 * 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.
950 * Formula used is: vect(AB) x vect(AC) = normal.
951 * Normalize the normal vector obtained and check with three random points from the PMT plane, whether the scalar product is equal to zero.
952 */
953 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
954 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
955 pmtTrackID = pmtPoint->GetTrackID();
956 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
957 pmtMotherID = track->GetMotherId();
958 a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
959 //cout << "a[0] = " << a[0] << ", a[1] = " << a[1] << " et a[2] = " << a[2] << endl;
960 break;
961 }
962 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
963 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
964 pmtTrackID = pmtPoint->GetTrackID();
965 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
966 pmtMotherID = track->GetMotherId();
967 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
968 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
969 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
970 > 7) {
971 b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
972 //cout << "b[0] = " << b[0] << ", b[1] = " << b[1] << " et b[2] = " << b[2] << endl;
973 break;
974 }
975 }
976 for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
977 CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
978 pmtTrackID = pmtPoint->GetTrackID();
979 track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
980 pmtMotherID = track->GetMotherId();
981 //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
982 if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2) + TMath::Power(a[1] - pmtPoint->GetY(), 2)
983 + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
984 > 7
985 && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2) + TMath::Power(b[1] - pmtPoint->GetY(), 2)
986 + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
987 > 7) {
988 c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
989 //cout << "c[0] = " << c[0] << ", c[1] = " << c[1] << " et c[2] = " << c[2] << endl;
990 break;
991 }
992 }
993
994 k = (b[0] - a[0]) / (c[0] - a[0]);
995 if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0 || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
996 cout << "Error in normal calculation, vect_AB and vect_AC are collinear." << endl;
997 }
998 else {
999 buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
1000 buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
1001 buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
1002 normalX =
1003 buffNormX / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
1004 normalY =
1005 buffNormY / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
1006 normalZ =
1007 buffNormZ / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2) + TMath::Power(buffNormZ, 2));
1008 }
1009
1010 CbmRichPoint* pmtPoint1 = (CbmRichPoint*) fRichPoints->At(20);
1011 scalarProd =
1012 normalX * (pmtPoint1->GetX() - a[0]) + normalY * (pmtPoint1->GetY() - a[1]) + normalZ * (pmtPoint1->GetZ() - a[2]);
1013 //cout << "1st scalar product between vectAM and normale = " << scalarProd << endl;
1014 // 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.
1015 normalCste = -1 * (normalX * pmtPoint1->GetX() + normalY * pmtPoint1->GetY() + normalZ * pmtPoint1->GetZ());
1016 CbmRichPoint* pmtPoint2 = (CbmRichPoint*) fRichPoints->At(15);
1017 scalarProd =
1018 normalX * (pmtPoint2->GetX() - a[0]) + normalY * (pmtPoint2->GetY() - a[1]) + normalZ * (pmtPoint2->GetZ() - a[2]);
1019 //cout << "2nd scalar product between vectAM and normale = " << scalarProd << endl;
1020 CbmRichPoint* pmtPoint3 = (CbmRichPoint*) fRichPoints->At(25);
1021 scalarProd =
1022 normalX * (pmtPoint3->GetX() - a[0]) + normalY * (pmtPoint3->GetY() - a[1]) + normalZ * (pmtPoint3->GetZ() - a[2]);
1023 //cout << "3nd scalar product between vectAM and normale = " << scalarProd << endl;
1024}
1025
1026void CbmRichPMTMapping::CalculateSphereParameters(const Char_t* mirrID, Double_t& sphereX, Double_t& sphereY,
1027 Double_t& sphereZ, Double_t& sphereR)
1028{
1029 //cout << endl << "//------------------------------ CbmRichPMTMapping: Calculate Sphere Parameters ------------------------------//" << endl << endl;
1030
1031 const Char_t* mirrorHalf;
1032 if (fIsMirrorUpperHalf) {
1033 mirrorHalf = "RICH_mirror_half_total_208";
1034 }
1035 else {
1036 mirrorHalf = "RICH_mirror_half_total_207";
1037 }
1038 //cout << "Mirror half: " << mirrorHalf << " and mirrID = " << mirrID << endl;
1039
1040 TObjArray* nodesTop = gGeoManager->GetTopNode()->GetNodes();
1041 for (Int_t i1 = 0; i1 < nodesTop->GetEntriesFast(); i1++) {
1042 TGeoNode* richNode = (TGeoNode*) nodesTop->At(i1);
1043 if (TString(richNode->GetName()).Contains("rich")) {
1044 const Double_t* trRich = richNode->GetMatrix()->GetTranslation();
1045 TObjArray* nodes2 = richNode->GetNodes();
1046 for (Int_t i2 = 0; i2 < nodes2->GetEntriesFast(); i2++) {
1047 TGeoNode* gasNode = (TGeoNode*) nodes2->At(i2);
1048 if (TString(gasNode->GetName()).Contains("RICH_gas")) {
1049 const Double_t* trGas = gasNode->GetMatrix()->GetTranslation();
1050 TObjArray* nodes3 = gasNode->GetNodes();
1051 for (Int_t i3 = 0; i3 < nodes3->GetEntriesFast(); i3++) {
1052 TGeoNode* mirrorHalfNode = (TGeoNode*) nodes3->At(i3);
1053 if (TString(mirrorHalfNode->GetName()).Contains(mirrorHalf)) {
1054 const Double_t* rotMirror = mirrorHalfNode->GetMatrix()->GetRotationMatrix();
1055 //gp.fMirrorTheta = TMath::ASin(rm[3]); // tilting angle around x-axis
1056 //gp.fPmtPhi = -1.*TMath::ASin(rm[2]); // tilting angle around y-axis
1057 const Double_t* trHalfMirror = mirrorHalfNode->GetMatrix()->GetTranslation();
1058 const TGeoBBox* mirrorShape = (const TGeoBBox*) (mirrorHalfNode->GetVolume()->GetShape());
1059 TObjArray* nodes4 = mirrorHalfNode->GetNodes();
1060 for (Int_t i4 = 0; i4 < nodes4->GetEntriesFast(); i4++) {
1061 TGeoNode* suppBeltStripNode = (TGeoNode*) nodes4->At(i4);
1062 if (TString(suppBeltStripNode->GetName()).Contains("RICH_mirror_and_support_belt_strip")) {
1063 const Double_t* trSuppBeltStrip = suppBeltStripNode->GetMatrix()->GetTranslation();
1064 TObjArray* nodes5 = suppBeltStripNode->GetNodes();
1065 for (Int_t i5 = 0; i5 < nodes5->GetEntriesFast(); i5++) {
1066 TGeoNode* mirrorTileNode = (TGeoNode*) nodes5->At(i5);
1067 //if( TString(mirrorTileNode->GetName()).Contains("RICH_mirror_1") || TString(mirrorTileNode->GetName()).Contains("RICH_mirror_2") || TString(mirrorTileNode->GetName()).Contains("RICH_mirror_3") ) {
1068 if (TString(mirrorTileNode->GetName()).Contains(mirrID)) {
1069 //cout << "mirrorTileNode->GetName() => " << mirrorTileNode->GetName() << endl;
1070 const Double_t* trMirrorTile = mirrorTileNode->GetMatrix()->GetTranslation();
1071 sphereX = trRich[0] + trGas[0] + trHalfMirror[0] + trSuppBeltStrip[0] + trMirrorTile[0];
1072 sphereY = trRich[1] + trGas[1] + trHalfMirror[1] + trSuppBeltStrip[1] + trMirrorTile[1];
1073 sphereZ = trRich[2] + trGas[2] + trHalfMirror[2] + trSuppBeltStrip[2]
1074 + trMirrorTile[2]; // + mirrorShape->GetDZ();
1075 /*
1076 * The actual translation, using corrected transformation matrices.
1077 * sphereX = trRich[0] + trGas[0] + trHalfMirror[0] + trSuppBeltStrip[0] + trMirrorTile[1];
1078 * sphereY = trRich[1] + trGas[1] + trHalfMirror[1] + trSuppBeltStrip[1] + trMirrorTile[2];
1079 * sphereZ = trRich[2] + trGas[2] + trHalfMirror[2] + trSuppBeltStrip[2] + trMirrorTile[0]; // + mirrorShape->GetDZ();
1080 */
1081 TGeoShape* ptrShape = mirrorTileNode->GetVolume()->GetShape();
1082 TGeoSphere* ptrSphere = static_cast<TGeoSphere*>(ptrShape);
1083 sphereR = ptrSphere->GetRmin();
1084 }
1085 }
1086 }
1087 }
1088 }
1089 }
1090 }
1091 }
1092 }
1093 }
1094}
1095
1096void CbmRichPMTMapping::CalculateSphereParameters2(const Char_t* mirrID, Double_t& sphereX, Double_t& sphereY,
1097 Double_t& sphereZ, Double_t& sphereR)
1098{
1099 //cout << endl << "//------------------------------ CbmRichPMTMapping: Calculate Sphere Parameters ------------------------------//" << endl << endl;
1100
1101 const Char_t* mirrorHalf;
1102 if (fIsMirrorUpperHalf) {
1103 mirrorHalf = "RICH_mirror_half_total_208";
1104 }
1105 else {
1106 mirrorHalf = "RICH_mirror_half_total_207";
1107 }
1108 //cout << "Mirror half: " << mirrorHalf << " and mirrID = " << mirrID << endl;
1109
1110 TObjArray* nodesTop = gGeoManager->GetTopNode()->GetNodes();
1111 for (Int_t i1 = 0; i1 < nodesTop->GetEntriesFast(); i1++) {
1112 TGeoNode* richNode = (TGeoNode*) nodesTop->At(i1);
1113 if (TString(richNode->GetName()).Contains("rich")) {
1114 TGeoMatrix* matRich = richNode->GetMatrix();
1115 /*cout << "Matrix richNode:" << endl;
1116 matRich->Print();*/
1117 const Double_t* trRich = richNode->GetMatrix()->GetTranslation();
1118 TObjArray* nodes2 = richNode->GetNodes();
1119 for (Int_t i2 = 0; i2 < nodes2->GetEntriesFast(); i2++) {
1120 TGeoNode* gasNode = (TGeoNode*) nodes2->At(i2);
1121 if (TString(gasNode->GetName()).Contains("RICH_gas")) {
1122 TGeoMatrix* matRichGas = gasNode->GetMatrix();
1123 /*cout << "Matrix gasNode:" << endl;
1124 matRichGas->Print();*/
1125 const Double_t* trGas = gasNode->GetMatrix()->GetTranslation();
1126 TObjArray* nodes3 = gasNode->GetNodes();
1127 for (Int_t i3 = 0; i3 < nodes3->GetEntriesFast(); i3++) {
1128 TGeoNode* mirrorHalfNode = (TGeoNode*) nodes3->At(i3);
1129 if (TString(mirrorHalfNode->GetName()).Contains(mirrorHalf)) {
1130 TGeoMatrix* matMirrorHalf = mirrorHalfNode->GetMatrix();
1131 /*cout << "Matrix mirrorHalfNode:" << endl;
1132 matMirrorHalf->Print();*/
1133 const Double_t* rotMirrorHalf = mirrorHalfNode->GetMatrix()->GetRotationMatrix();
1134 //gp.fMirrorTheta = TMath::ASin(rm[3]); // tilting angle around x-axis
1135 //gp.fPmtPhi = -1.*TMath::ASin(rm[2]); // tilting angle around y-axis
1136 const Double_t* trHalfMirror = mirrorHalfNode->GetMatrix()->GetTranslation();
1137 const TGeoBBox* mirrorShape = (const TGeoBBox*) (mirrorHalfNode->GetVolume()->GetShape());
1138 TObjArray* nodes4 = mirrorHalfNode->GetNodes();
1139 for (Int_t i4 = 0; i4 < nodes4->GetEntriesFast(); i4++) {
1140 TGeoNode* suppBeltStripNode = (TGeoNode*) nodes4->At(i4);
1141 if (TString(suppBeltStripNode->GetName()).Contains("RICH_mirror_and_support_belt_strip")) {
1142 TGeoMatrix* matSuppBeltStrip = suppBeltStripNode->GetMatrix();
1143 /*cout << "Matrix suppBeltStripNode:" << suppBeltStripNode->GetName() << endl;
1144 matSuppBeltStrip->Print();*/
1145 const Double_t* trSuppBeltStrip = suppBeltStripNode->GetMatrix()->GetTranslation();
1146 TObjArray* nodes5 = suppBeltStripNode->GetNodes();
1147 for (Int_t i5 = 0; i5 < nodes5->GetEntriesFast(); i5++) {
1148 TGeoNode* mirrorTileNode = (TGeoNode*) nodes5->At(i5);
1149 //if( TString(mirrorTileNode->GetName()).Contains("RICH_mirror_1") || TString(mirrorTileNode->GetName()).Contains("RICH_mirror_2") || TString(mirrorTileNode->GetName()).Contains("RICH_mirror_3") ) {
1150 if (TString(mirrorTileNode->GetName()).Contains(mirrID)) {
1151 //cout << "mirrorTileNode->GetName() => " << mirrorTileNode->GetName() << endl;
1152 /*TGeoMatrix* matMirrorTile = mirrorTileNode->GetMatrix();
1153 cout << "Matrix mirrorTileNode:" << endl;
1154 matMirrorTile->Print();*/
1155 const Double_t* trMirrorTile = mirrorTileNode->GetMatrix()->GetTranslation();
1156 const Double_t* rotMirrorTile = mirrorTileNode->GetMatrix()->GetRotationMatrix();
1157 sphereX = trRich[0] + trGas[0] + trHalfMirror[0] + trSuppBeltStrip[0] + trMirrorTile[1];
1158 sphereY = trRich[1] + trGas[1] + trHalfMirror[1] + trSuppBeltStrip[1] + trMirrorTile[2];
1159 sphereZ = trRich[2] + trGas[2] + trHalfMirror[2] + trSuppBeltStrip[2]
1160 + trMirrorTile[0]; // + mirrorShape->GetDZ();
1161 TGeoShape* ptrShape = mirrorTileNode->GetVolume()->GetShape();
1162 TGeoSphere* ptrSphere = static_cast<TGeoSphere*>(ptrShape);
1163 sphereR = ptrSphere->GetRmin();
1164 }
1165 }
1166 }
1167 }
1168 }
1169 }
1170 }
1171 }
1172 }
1173 }
1174}
1175
1176void CbmRichPMTMapping::RotateAndCopyHitsToRingLight(const CbmRichRing* ring1, CbmRichRingLight* ring2)
1177{
1178 Int_t nofHits = ring1->GetNofHits();
1179
1180 for (Int_t i = 0; i < nofHits; i++) {
1181 Int_t hitInd = ring1->GetHit(i);
1182 CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
1183 if (NULL == hit) continue;
1184 TVector3 inputHit(hit->GetX(), hit->GetY(), hit->GetZ());
1185 TVector3 outputHit;
1186 CbmRichHitProducer::TiltPoint(&inputHit, &outputHit, fGP.fPmt.fPhi, fGP.fPmt.fTheta, fGP.fPmtZOrig);
1187 CbmRichHitLight hl(outputHit.X(), outputHit.Y());
1188 ring2->AddHit(hl);
1189 }
1190}
1191
1192void CbmRichPMTMapping::DrawHist()
1193{
1194 TCanvas* can = new TCanvas(fRunTitle + "_Separated_Hits", fRunTitle + "_Separated_Hits", 1500, 900);
1195 can->SetGridx(1);
1196 can->SetGridy(1);
1197 can->Divide(4, 3);
1198 can->cd(9);
1199 DrawH2(fHM->H2("fHMCPoints_3_15"));
1200 can->cd(5);
1201 DrawH2(fHM->H2("fHMCPoints_2_16"));
1202 can->cd(1);
1203 DrawH2(fHM->H2("fHMCPoints_2_17"));
1204 can->cd(2);
1205 DrawH2(fHM->H2("fHMCPoints_2_29"));
1206 can->cd(3);
1207 DrawH2(fHM->H2("fHMCPoints_2_53"));
1208 can->cd(4);
1209 DrawH2(fHM->H2("fHMCPoints_2_77"));
1210 can->cd(6);
1211 DrawH2(fHM->H2("fHMCPoints_2_28"));
1212 can->cd(7);
1213 DrawH2(fHM->H2("fHMCPoints_2_52"));
1214 can->cd(8);
1215 DrawH2(fHM->H2("fHMCPoints_2_76"));
1216 can->cd(10);
1217 DrawH2(fHM->H2("fHMCPoints_1_27"));
1218 can->cd(11);
1219 DrawH2(fHM->H2("fHMCPoints_1_51"));
1220 can->cd(12);
1221 DrawH2(fHM->H2("fHMCPoints_1_75"));
1222 Cbm::SaveCanvasAsImage(can, string(fOutputDir.Data()), "png");
1223
1224 TCanvas* can2 = new TCanvas(fRunTitle + "_Separated_Ellipse", fRunTitle + "_Separated_Ellipse", 1500, 900);
1225 can2->SetGridx(1);
1226 can2->SetGridy(1);
1227 can2->Divide(4, 3);
1228 can2->cd(9);
1229 DrawH2(fHM->H2("fHPoints_Ellipse_3_15"));
1230 can2->cd(5);
1231 DrawH2(fHM->H2("fHPoints_Ellipse_2_16"));
1232 can2->cd(1);
1233 DrawH2(fHM->H2("fHPoints_Ellipse_2_17"));
1234 can2->cd(2);
1235 DrawH2(fHM->H2("fHPoints_Ellipse_2_29"));
1236 can2->cd(3);
1237 DrawH2(fHM->H2("fHPoints_Ellipse_2_53"));
1238 can2->cd(4);
1239 DrawH2(fHM->H2("fHPoints_Ellipse_2_77"));
1240 can2->cd(6);
1241 DrawH2(fHM->H2("fHPoints_Ellipse_2_28"));
1242 can2->cd(7);
1243 DrawH2(fHM->H2("fHPoints_Ellipse_2_52"));
1244 can2->cd(8);
1245 DrawH2(fHM->H2("fHPoints_Ellipse_2_76"));
1246 can2->cd(10);
1247 DrawH2(fHM->H2("fHPoints_Ellipse_1_27"));
1248 can2->cd(11);
1249 DrawH2(fHM->H2("fHPoints_Ellipse_1_51"));
1250 can2->cd(12);
1251 DrawH2(fHM->H2("fHPoints_Ellipse_1_75"));
1252 Cbm::SaveCanvasAsImage(can2, string(fOutputDir.Data()), "png");
1253
1254 TCanvas* can3 = new TCanvas(fRunTitle + "_Separated_Ellipse", fRunTitle + "_Separated_Ellipse", 1500, 900);
1255 can3->Divide(2, 2);
1256 can3->cd(1);
1257 DrawH1(fHM->H1("fhDifferenceX"));
1258 can3->cd(2);
1259 DrawH1(fHM->H1("fhDifferenceY"));
1260 can3->cd(3);
1261 DrawH1(fHM->H1("fhDistanceCenterToExtrapolatedTrack"));
1262 can3->cd(4);
1263 DrawH1(fHM->H1("fhDistanceCenterToExtrapolatedTrackInPlane"));
1264}
1265
1266void CbmRichPMTMapping::DrawHistFromFile(TString fileName)
1267{
1269 TFile* oldFile = gFile;
1270 TDirectory* oldDir = gDirectory;
1271
1272 fHM = new CbmHistManager();
1273 TFile* file = new TFile(fileName, "READ");
1274 fHM->ReadFromFile(file);
1275 DrawHist();
1276
1278 gFile = oldFile;
1279 gDirectory = oldDir;
1280}
1281
1282void CbmRichPMTMapping::Finish()
1283{
1284 // ---------------------------------------------------------------------------------------------------------------------------------------- //
1285 // -------------------------------------------------- Mapping for mirror - PMT relations -------------------------------------------------- //
1286 // ---------------------------------------------------------------------------------------------------------------------------------------- //
1287
1288 if (fDrawHist) {
1289 DrawHist();
1290 }
1291 cout << endl << "Mirror counter = " << fMirrCounter << endl;
1292 //cout << setprecision(6) << endl;
1293}
1294ClassImp(CbmRichPMTMapping)
1295
1296 /* Old Code:
1297Double_t x_PMT, y_PMT, z_PMT;
1298Int_t nofMirrorPoints = fRichMirrorPoints->GetEntriesFast();
1299Int_t NofProjections = fRichProjections->GetEntriesFast();
1300cout << "Nb of Mirr_Pts: " << nofMirrorPoints << " and nb of Projections: " << NofProjections << endl;
1301//FairVolume* vol;
1302//vol->GetName();
1303if (nofMirrorPoints >= 1) {
1304 for (Int_t iMP = 0; iMP < nofMirrorPoints; iMP++) {
1305 CbmRichPoint *mirrorPoint = (CbmRichPoint*)fRichMirrorPoints->At(iMP);
1306 Double_t xMirr = mirrorPoint->GetX();
1307 Double_t yMirr = mirrorPoint->GetY();
1308 Double_t zMirr = mirrorPoint->GetZ();
1309 cout << "Particle hit coordinates on mirror: X = " << xMirr << ", Y = " << yMirr << " and Z = " << zMirr << endl;
1310
1311 TGeoNode *current_node = gGeoManager->GetCurrentNode();
1312 TGeoVolume *current_vol = current_node->GetVolume();
1313 // or: TGeoVolume *cvol = gGeoManager->GetCurrentVolume();
1314 TGeoMaterial *current_mat = current_vol->GetMedium()->GetMaterial();
1315
1316 TGeoNode *mirr_node = gGeoManager->FindNode(xMirr, yMirr, zMirr);
1317 TGeoVolume *mirr_vol = mirr_node->GetVolume();
1318 TGeoMaterial *mirr_mat = mirr_vol->GetMedium()->GetMaterial();
1319
1320 const Char_t *c1, *c2, *v1, *m1, *c3, *v2, *m2;
1321 c1 = mirr_node->GetName();
1322 c2 = current_node->GetName();
1323 cout << "NAMES:" << endl << "Node name mirr: " << c1 << " and current_node name: " << c2 << endl;
1324 v1 = mirr_vol->GetName();
1325 m1 = mirr_mat->GetName();
1326 v2 = current_vol->GetName();
1327 m2 = current_mat->GetName();
1328 cout << "Volume mirr: " << v1 << " and material mirr: " << m1 << endl;
1329 cout << "Current volume: " << v2 << " and current material: " << m2 << endl;
1330 Int_t Index_1 = mirr_node->GetIndex();
1331 Int_t Index_2 = current_node->GetIndex();
1332 cout << "Index mirr: " << Index_1 << " and current index: " << Index_2 << endl;
1333 //current->Draw("");
1334 //node->Draw("");
1335
1336 const Char_t *path = gGeoManager->GetPath();
1337 cout << "Current path is: " << path << endl;
1338
1339 for (Int_t iP = 0; iP < NofProjections; iP++) {
1340 FairTrackParam* pr = (FairTrackParam*) fRichProjections->At(iP);
1341 x_PMT = pr->GetX();
1342 y_PMT = pr->GetY();
1343 z_PMT = pr->GetZ();
1344 cout << "Center x_PMT: " << x_PMT << ", center y_PMT: " << y_PMT << " and z_PMT: " << z_PMT << endl;
1345 }
1346
1347 TGeoNode *node_pmt = gGeoManager->FindNode(x_PMT, y_PMT, z_PMT);
1348 c3 = node_pmt->GetName();
1349 cout << "Node name pmt: " << c3 << endl << endl;
1350 //TGeoManager* Node = FindNode(xMirr, yMirr, zMirr);
1351 //sleep(2);
1352 }
1353}
1354for (Int_t iMirr = 0; iMirr < NofMirrorPoints; iMirr++) {
1355 CbmRichPoint* MirrPoint = (CbmRichPoint*) fRichMirrorPoints->At(iMirr);
1356 Int_t trackID = MirrPoint->GetTrackID();
1357 for (Int_t iMCPoint = 0; iMCPoint < NofMCPoints; iMCPoint++) {
1358 CbmRichPoint* pPoint = (CbmRichPoint*) fRichMCPoints->At(iMCPoint);
1359 if ( NULL == pPoint)
1360 continue;
1361 CbmMCTrack* pTrack = (CbmMCTrack*) fMCTracks->At(pPoint->GetTrackID());
1362 if ( NULL == pTrack)
1363 continue;
1364 Int_t gcode = pTrack->GetPdgCode();
1365 Int_t motherID = pTrack->GetMotherId();
1366 if (motherID == -1)
1367 continue;
1368 if (trackID == motherID) {
1369 //cout << "MATCH BETWEEN TRACK ID AND MOTHER ID FOUND !" << endl << "TrackID from mirror point = " << trackID << " and mother ID from MC point = " << motherID << endl;
1370 //sleep(2);
1371 // Get transformation matrix of mirror from mirrNode
1372 xMirr = MirrPoint->GetX();
1373 yMirr = MirrPoint->GetY();
1374 zMirr = MirrPoint->GetZ();
1375 //cout << "x Mirr: " << xMirr << ", y Mirr: " << yMirr << " and z Mirr: " << zMirr << endl;
1376 mirrNode = gGeoManager->FindNode(xMirr, yMirr, zMirr);
1377 mirrPath = gGeoManager->GetPath();
1378 mirrMatrix = mirrNode->GetMatrix();
1379 const Double_t trans = *mirrMatrix->GetTranslation();
1380 cout << endl << " !!! HERE 1 !!! " << endl << "Translation matrix = " << trans << endl;
1381 cout << "Rotation matrix = " << *mirrMatrix->GetRotationMatrix() << endl;
1382 //if (mirrMatrix->IsCombi()){mirrMatrix->Print();}
1383 // Get shape for local mirror rotations
1384 vol = mirrNode->GetVolume();
1385 ptrShape = vol->GetShape();
1386 ptrSphere = static_cast<TGeoSphere*> (ptrShape);
1387 phi1 = ptrSphere->GetPhi1();
1388 phi2 = ptrSphere->GetPhi2();
1389 theta1 = ptrSphere->GetTheta1();
1390 theta2 = ptrSphere->GetTheta2();
1391 // Get transformation matrix of PMT plane from pmtNode
1392 xPMT = pPoint->GetX();
1393 yPMT = pPoint->GetY();
1394 zPMT = pPoint->GetZ();
1395 pmtNode = gGeoManager->FindNode(xPMT, yPMT, zPMT);
1396 pmtMatrix = pmtNode->GetMatrix();
1397 //fGP = CbmRichHitProducer::InitGeometry();
1398 //fGP.Print();
1399 }
1400 }
1401}*/
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.
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
CbmRichHitProducer::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition CbmRichHitProducer.cxx:44
CbmRichPMTMapping::RotateAndCopyHitsToRingLight
void RotateAndCopyHitsToRingLight(const CbmRichRing *ring1, CbmRichRingLight *ring2)
Definition CbmRichPMTMapping.cxx:1176
CbmRichPMTMapping::FillPMTMapEllipse
void FillPMTMapEllipse(const Char_t *mirr_path, Float_t CenterX, Float_t CenterY)
Definition CbmRichPMTMapping.cxx:387
CbmRichPMTMapping::fRichPoints
TClonesArray * fRichPoints
Definition CbmRichPMTMapping.h:138
CbmRichPMTMapping::fGlobalTracks
TClonesArray * fGlobalTracks
Definition CbmRichPMTMapping.h:139
CbmRichPMTMapping::fHM
CbmHistManager * fHM
Definition CbmRichPMTMapping.h:140
CbmRichPMTMapping::fRichRingMatches
TClonesArray * fRichRingMatches
Definition CbmRichPMTMapping.h:136
CbmRichPMTMapping::CalculateSphereParameters
void CalculateSphereParameters(const Char_t *mirrID, Double_t &sphereX, Double_t &sphereY, Double_t &sphereZ, Double_t &sphereR)
Definition CbmRichPMTMapping.cxx:1026
CbmRichPMTMapping::fRichMirrorPoints
TClonesArray * fRichMirrorPoints
Definition CbmRichPMTMapping.h:132
CbmRichPMTMapping::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition CbmRichPMTMapping.cxx:231
CbmRichPMTMapping::fMCTracks
TClonesArray * fMCTracks
Definition CbmRichPMTMapping.h:135
CbmRichPMTMapping::fRichMCPoints
TClonesArray * fRichMCPoints
Definition CbmRichPMTMapping.h:134
CbmRichPMTMapping::fPathsMap
std::map< string, string > fPathsMap
Definition CbmRichPMTMapping.h:150
CbmRichPMTMapping::Finish
virtual void Finish()
Inherited from FairTask.
Definition CbmRichPMTMapping.cxx:1282
CbmRichPMTMapping::fTauFit
CbmRichRingFitterEllipseTau * fTauFit
Definition CbmRichPMTMapping.h:157
CbmRichPMTMapping::CalculateSphereParameters2
void CalculateSphereParameters2(const Char_t *mirrID, Double_t &sphereX, Double_t &sphereY, Double_t &sphereZ, Double_t &sphereR)
Definition CbmRichPMTMapping.cxx:1096
CbmRichPMTMapping::fRichRings
TClonesArray * fRichRings
Definition CbmRichPMTMapping.h:131
CbmRichPMTMapping::GetPmtNormal
void GetPmtNormal(Int_t NofPMTPoints, Double_t &normalX, Double_t &normalY, Double_t &normalZ, Double_t &normalCste)
Definition CbmRichPMTMapping.cxx:937
CbmRichPMTMapping::fRichRefPlanePoints
TClonesArray * fRichRefPlanePoints
Definition CbmRichPMTMapping.h:137
CbmRichPMTMapping::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition CbmRichPMTMapping.cxx:85
CbmRichPMTMapping::fPathsMapEllipse
std::map< string, string > fPathsMapEllipse
Definition CbmRichPMTMapping.h:151
CbmRichPMTMapping::fGP
CbmRichRecGeoPar fGP
Definition CbmRichPMTMapping.h:141
CbmRichPMTMapping::DrawHistFromFile
void DrawHistFromFile(TString fileName)
Definition CbmRichPMTMapping.cxx:1266
CbmRichPMTMapping::FillPMTMap
void FillPMTMap(const Char_t *mirr_path, CbmRichPoint *pPoint)
Definition CbmRichPMTMapping.cxx:368
CbmRichPMTMapping::fRichHits
TClonesArray * fRichHits
Definition CbmRichPMTMapping.h:130
CbmRichPMTMapping::fRichProjections
TClonesArray * fRichProjections
Definition CbmRichPMTMapping.h:133
CbmRichPMTMapping::fCopFit
CbmRichRingFitterCOP * fCopFit
Definition CbmRichPMTMapping.h:156
CbmRichPoint::Print
virtual void Print(const Option_t *opt) const
Definition CbmRichPoint.cxx:38
CbmRichRecGeoPar::fPmt
CbmRichRecGeoParPmt fPmt
Definition CbmRichRecGeoPar.h:218
CbmRichRecGeoPar::Print
void Print()
Print geometry parameters.
Definition CbmRichRecGeoPar.h:123
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::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