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