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