CbmRoot
CbmRichMirrorSortingAlignment.cxx
Go to the documentation of this file.
1 #include "CbmRichMirrorSortingAlignment.h"
2 #include "FairLogger.h"
3 #include "FairRootManager.h"
4 
5 // ----- PART 1 ----- //
6 #include "CbmGlobalTrack.h"
7 #include "CbmMCTrack.h"
8 #include "CbmRichConverter.h"
9 #include "CbmRichRing.h"
10 #include "CbmRichRingFitterCOP.h"
11 #include "CbmRichRingFitterEllipseTau.h"
12 #include "CbmUtils.h"
13 #include "TCanvas.h"
14 #include "TClonesArray.h"
15 #include "TF1.h"
16 #include "TH2D.h"
17 #include "TLegend.h"
18 #include "TVector3.h"
19 #include <vector>
20 
21 // ----- PART 2 ----- //
22 #include "CbmRichGeoManager.h"
23 #include "CbmRichPoint.h"
24 #include "TGeoManager.h"
25 #include "TGeoNavigator.h"
26 class TGeoNode;
27 class TGeoMatrix;
28 #include "CbmRichNavigationUtil.h"
29 #include "CbmTrackMatchNew.h"
30 #include "TStyle.h"
31 #include "string.h"
32 
33 #include "TSystem.h"
34 #include <iostream>
35 
36 CbmRichMirrorSortingAlignment::CbmRichMirrorSortingAlignment()
37  : FairTask("CbmRichMirrorSortingAlignment")
38  , fEventNb(0)
39  , fGlobalTracks(NULL)
40  , fRichRings(NULL)
41  , fMCTracks(NULL)
42  , fCopFit(NULL)
43  , fTauFit(NULL)
44  , fOutputDir("")
45  , fStudyName("")
46  , fMirrorPoints(NULL)
47  , fRefPlanePoints(NULL)
48  , fPmtPoints(NULL)
49  , fRichProjections(NULL)
50  , fTrackParams(NULL)
51  , fRichRingMatches(NULL)
52  , fStsTrackMatches(NULL)
53  , fMirrorMap()
54  , fThreshold(0) {}
55 
56 CbmRichMirrorSortingAlignment::~CbmRichMirrorSortingAlignment() {}
57 
58 InitStatus CbmRichMirrorSortingAlignment::Init() {
59  FairRootManager* manager = FairRootManager::Instance();
60 
61  fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
62  if (NULL == fGlobalTracks) {
63  Fatal("CbmRichMirrorSortingAlignment::Init", "No GlobalTrack array!");
64  }
65 
66  fRichRings = (TClonesArray*) manager->GetObject("RichRing");
67  if (NULL == fRichRings) {
68  Fatal("CbmRichMirrorSortingAlignment::Init", "No RichRing array !");
69  }
70 
71  fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
72  if (NULL == fMCTracks) {
73  Fatal("CbmRichMirrorSortingAlignment::Init", "No MCTracks array !");
74  }
75 
76  fMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
77  if (NULL == fMirrorPoints) {
78  Fatal("CbmRichMirrorSortingAlignment::Init", "No RichMirrorPoints array !");
79  }
80 
81  fRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
82  if (NULL == fRefPlanePoints) {
83  Fatal("CbmRichMirrorSortingAlignment::Init",
84  "No RichRefPlanePoint array !");
85  }
86 
87  fPmtPoints = (TClonesArray*) manager->GetObject("RichPoint");
88  if (NULL == fPmtPoints) {
89  Fatal("CbmRichMirrorSortingAlignment::Init", "No RichPoint array !");
90  }
91 
92  fRichProjections = (TClonesArray*) manager->GetObject("RichProjection");
93  if (NULL == fRichProjections) {
94  Fatal("CbmRichMirrorSortingAlignment::Init", "No RichProjection array !");
95  }
96 
97  fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
98  if (NULL == fTrackParams) {
99  Fatal("CbmRichMirrorSortingAlignment::Init", "No RichTrackParamZ array!");
100  }
101 
102  fRichRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
103  if (NULL == fRichRingMatches) {
104  Fatal("CbmRichMirrorSortingAlignment::Init", "No RichRingMatch array!");
105  }
106 
107  fStsTrackMatches = (TClonesArray*) manager->GetObject("StsTrackMatch");
108  if (NULL == fStsTrackMatches) {
109  Fatal("CbmRichMirrorSortingAlignment::Init", "No StsTrackMatch array!");
110  }
111 
112  fCopFit = new CbmRichRingFitterCOP();
113  fTauFit = new CbmRichRingFitterEllipseTau();
114  CbmRichConverter::Init();
115 
116  return kSUCCESS;
117 }
118 
119 void CbmRichMirrorSortingAlignment::Exec(Option_t* Option) {
120  fEventNb++;
121  cout << "CbmRichMirrorSortingAlignment: Event #" << fEventNb << endl;
122  TVector3 momentum, outPos;
123  Double_t constantePMT = 0., trackX = 0., trackY = 0.;
124  vector<Double_t> vect(2, 0), ptM(3, 0), ptC(3, 0), ptCIdeal(3, 0), ptR1(3, 0),
125  ptR2Center(3, 0), ptR2Mirr(3, 0), ptPR2(3, 0), ptPMirr(3, 0),
126  normalPMT(3, 0);
127  //ptC[0]=0.;
128  ptC.at(0) = 0., ptC.at(1) = 132.594000, ptC.at(2) = 54.267226;
129  TVector3 mirrorPoint, dirCos, pos;
130  Double_t nx = 0., ny = 0., nz = 0.;
131  TGeoNode* mirrNode;
132  CbmRichPoint *mirrPoint, *refPlanePoint;
133 
134  if (fRichRings->GetEntries() != 0) {
135  cout << "Nb of rings in evt = " << fRichRings->GetEntries() << endl << endl;
136  GetPmtNormal(fPmtPoints->GetEntries(), normalPMT, constantePMT);
137  //cout << "Calculated normal vector to PMT plane = {" << normalPMT.at(0) << ", " << normalPMT.at(1) << ", " << normalPMT.at(2) << "} and constante d = " << constantePMT << endl;
138 
139  for (Int_t iGlobalTrack = 0; iGlobalTrack < fGlobalTracks->GetEntriesFast();
140  iGlobalTrack++) {
141  CbmRichMirror* mirrorObject =
142  new CbmRichMirror(); // Create CbmRichMirror object, to be later filled.
143  // ----- PART 1 ----- //
144  // Ring-Track matching + Ring fit + Track momentum:
145  CbmGlobalTrack* gTrack =
146  (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack);
147  Int_t richInd = gTrack->GetRichRingIndex();
148  Int_t stsInd = gTrack->GetStsTrackIndex();
149  //cout << "richInd: " << richInd << endl;
150  if (richInd < 0) {
151  cout << "Error richInd < 0" << endl;
152  continue;
153  }
154  CbmRichRing* ring = (CbmRichRing*) fRichRings->At(richInd);
155  if (ring == NULL) {
156  cout << "Error ring == NULL!" << endl;
157  continue;
158  }
159  //Int_t ringTrackID = ring->GetTrackID(); //Old code not working with changes for new matching method.
160  //cout << "ringTrackID: " << ringTrackID << endl;
161  CbmTrackMatchNew* cbmRichTrackMatch =
162  (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
163  CbmTrackMatchNew* cbmStsTrackMatch =
164  (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
165  if (NULL == cbmRichTrackMatch) { continue; }
166  cout << "Nof true hits = " << cbmRichTrackMatch->GetNofTrueHits() << endl;
167  cout << "Nof wrong hits = " << cbmRichTrackMatch->GetNofWrongHits()
168  << endl;
169  Int_t mcRichTrackId = cbmRichTrackMatch->GetMatchedLink().GetIndex();
170  Int_t mcStsTrackId = cbmStsTrackMatch->GetMatchedLink().GetIndex();
171  //cout << "mcTrackId: " << mcRichTrackId << endl;
172  if (mcRichTrackId < 0) continue;
173  //if (mcStsTrackId != ringTrackID) { //Old code not working with changes for new matching method.
174  if (mcStsTrackId != mcRichTrackId) {
175  cout << "Error StsTrackIndex and TrackIndex from Ring do not match!"
176  << endl;
177  continue;
178  }
179  CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcRichTrackId);
180  if (!mcTrack) continue;
181 
182  CbmRichRingLight ringL;
183  CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
184  fCopFit->DoFit(&ringL);
185  //fTauFit->DoFit(&ringL);
186  mirrorObject->setRingLight(
187  ringL); // Fill Cbm Rich Ring Light inside mirrorObject.
188  cout << "ring Center Coo: " << ringL.GetCenterX() << ", "
189  << ringL.GetCenterY() << endl;
190  mcTrack->GetMomentum(momentum);
191  mirrorObject->setMomentum(
192  momentum); // Fill track momentum inside mirrorObject.
193  //FairTrackParam* pr = (FairTrackParam*) fRichProjections->At(ringTrackID); //Old code not working with changes for new matching method.
194  FairTrackParam* pr = (FairTrackParam*) fRichProjections->At(stsInd);
195  FairTrackParam* point = (FairTrackParam*) fTrackParams->At(stsInd);
196  if (pr == NULL) {
197  cout << "CbmRichMirrorSortingAlignment::Exec : pr = NULL." << endl;
198  continue;
199  }
200  trackX = pr->GetX(), trackY = pr->GetY();
201  cout << "Track: " << trackX << ", " << trackY << endl;
202  mirrorObject->setProjHit(trackX, trackY);
203 
204  // ----- PART 2 ----- //
205  // Mirror ID via TGeoNavigator + Extrap hit:
206  Int_t trackMotherId = mcTrack->GetMotherId();
207  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
208  if (trackMotherId == -1) {
209  if (fMirrorPoints->GetEntries() > 0) {
210  //loop on mirrorPoint and compare w/ TrackID->GetTrackId to get correct one
211  for (Int_t iMirrPt = 0; iMirrPt < fMirrorPoints->GetEntries();
212  iMirrPt++) {
213  mirrPoint = (CbmRichPoint*) fMirrorPoints->At(iMirrPt);
214  if (mirrPoint == 0) { continue; }
215  //cout << "Mirror point track ID: " << mirrPoint->GetTrackID() << endl;
216  if (mirrPoint->GetTrackID() == mcRichTrackId) { break; }
217  }
218  ptM.at(0) = mirrPoint->GetX(), ptM.at(1) = mirrPoint->GetY(),
219  ptM.at(2) = mirrPoint->GetZ();
220  cout << "mirrPoint: {" << mirrPoint->GetX() << ", "
221  << mirrPoint->GetY() << ", " << mirrPoint->GetZ() << "}" << endl;
222 
223  mirrNode = gGeoManager->FindNode(ptM.at(0), ptM.at(1), ptM.at(2));
224  cout << "Mirror node name: " << mirrNode->GetName()
225  << " and full path " << gGeoManager->GetPath() << endl;
226  string str1 = gGeoManager->GetPath(), str2 = "mirror_tile_",
227  str3 = "";
228  cout << "str1 before CbmRichNavigationUtil::FindIntersection: "
229  << str1 << endl;
230  TVector3 crossP;
231  str1 = CbmRichNavigationUtil::FindIntersection(
232  point, crossP, "mirror_tile_");
233  cout << "str1 after CbmRichNavigationUtil::FindIntersection: " << str1
234  << endl;
235 
236  std::size_t found = str1.find(str2);
237  if (found != std::string::npos) {
238  //cout << "first 'mirror_tile_type' found at: " << found << '\n';
239  Int_t end = str2.length() + 3;
240  str3 = str1.substr(found, end);
241  }
242  cout << "Mirror ID: " << str3 << endl;
243  mirrorObject->setMirrorId(str3);
244 
245  if (mirrNode) {
246  TGeoNavigator* navi = gGeoManager->GetCurrentNavigator();
247  //cout << "Navigator path: " << navi->GetPath() << endl;
248  ptCIdeal.at(0) = navi->GetCurrentMatrix()->GetTranslation()[0];
249  ptCIdeal.at(1) = navi->GetCurrentMatrix()->GetTranslation()[1];
250  ptCIdeal.at(2) = navi->GetCurrentMatrix()->GetTranslation()[2];
251  cout << "Sphere center coordinates of the aligned mirror tile, "
252  "ideal = {"
253  << ptCIdeal.at(0) << ", " << ptCIdeal.at(1) << ", "
254  << ptCIdeal.at(2) << "}" << endl;
255  for (Int_t iRefPt = 0; iRefPt < fRefPlanePoints->GetEntries();
256  iRefPt++) {
257  refPlanePoint = (CbmRichPoint*) fRefPlanePoints->At(iRefPt);
258  //cout << "Refl plane point track ID: " << refPlanePoint->GetTrackID() << endl;
259  if (refPlanePoint->GetTrackID() == mcRichTrackId) { break; }
260  }
261  ptR1.at(0) = refPlanePoint->GetX(),
262  ptR1.at(1) = refPlanePoint->GetY(),
263  ptR1.at(2) = refPlanePoint->GetZ();
264  cout << "Refl plane point coo = {" << ptR1[0] << ", " << ptR1[1]
265  << ", " << ptR1[2] << "}" << endl;
266  ComputeR2(
267  ptR2Center, ptR2Mirr, ptM, ptC, ptR1, navi, "Uncorrected");
268  ComputeP(ptPMirr, ptPR2, normalPMT, ptM, ptR2Mirr, constantePMT);
269  TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
270  cout << "inPos vector: " << inPos.x() << ", " << inPos.y() << ", "
271  << inPos.z() << endl;
272  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
273  cout << "New PMT points coordinates = {" << outPos.x() << ", "
274  << outPos.y() << ", " << outPos.z() << "}" << endl;
275  mirrorObject->setExtrapHit(outPos.x(), outPos.y());
276  }
277  } else {
278  cout << "No mirror points registered." << endl;
279  }
280  } else {
281  cout << "Not a mother particle." << endl;
282  }
283  //ComputeAngles();
284  fMirrorMap[mirrorObject->getMirrorId()].push_back(mirrorObject);
285  cout << "Key str: " << mirrorObject->getMirrorId() << endl
286  << "Mirror map: " << fMirrorMap[mirrorObject->getMirrorId()].size()
287  << endl
288  << endl;
289  //mirrNode->Clear();
290  //gGeoManager->Clear();
291  }
292  } else {
293  cout << "CbmRichMirrorSortingAlignment::Exec No rings in event were found."
294  << endl;
295  }
296 }
297 
298 void CbmRichMirrorSortingAlignment::GetPmtNormal(Int_t NofPMTPoints,
299  vector<Double_t>& normalPMT,
300  Double_t& normalCste) {
301  //cout << endl << "//------------------------------ CbmRichMirrorSortingAlignment: Calculate PMT Normal ------------------------------//" << endl << endl;
302 
303  Int_t pmtTrackID, pmtMotherID;
304  Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0.,
305  scalarProd = 0.;
306  Double_t pmtPt[] = {0., 0., 0.};
307  Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
308  CbmMCTrack* track;
309 
310  /*
311  * 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.
312  * Formula used is: vect(AB) x vect(AC) = normal.
313  * Normalize the normal vector obtained and check with three random points from the PMT plane, whether the scalar product is equal to zero.
314  */
315  for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
316  CbmRichPoint* pmtPoint = (CbmRichPoint*) fPmtPoints->At(iPmt);
317  pmtTrackID = pmtPoint->GetTrackID();
318  track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
319  pmtMotherID = track->GetMotherId();
320  a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
321  //cout << "a[0] = " << a[0] << ", a[1] = " << a[1] << " et a[2] = " << a[2] << endl;
322  break;
323  }
324  for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
325  CbmRichPoint* pmtPoint = (CbmRichPoint*) fPmtPoints->At(iPmt);
326  pmtTrackID = pmtPoint->GetTrackID();
327  track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
328  pmtMotherID = track->GetMotherId();
329  //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
330  if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2)
331  + TMath::Power(a[1] - pmtPoint->GetY(), 2)
332  + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
333  > 7) {
334  b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
335  //cout << "b[0] = " << b[0] << ", b[1] = " << b[1] << " et b[2] = " << b[2] << endl;
336  break;
337  }
338  }
339  for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
340  CbmRichPoint* pmtPoint = (CbmRichPoint*) fPmtPoints->At(iPmt);
341  pmtTrackID = pmtPoint->GetTrackID();
342  track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
343  pmtMotherID = track->GetMotherId();
344  //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
345  if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2)
346  + TMath::Power(a[1] - pmtPoint->GetY(), 2)
347  + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
348  > 7
349  && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2)
350  + TMath::Power(b[1] - pmtPoint->GetY(), 2)
351  + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
352  > 7) {
353  c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
354  //cout << "c[0] = " << c[0] << ", c[1] = " << c[1] << " et c[2] = " << c[2] << endl;
355  break;
356  }
357  }
358 
359  k = (b[0] - a[0]) / (c[0] - a[0]);
360  if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0
361  || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
362  cout << "Error in normal calculation, vect_AB and vect_AC are collinear."
363  << endl;
364  } else {
365  buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
366  buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
367  buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
368  normalPMT.at(0) =
369  buffNormX
370  / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
371  + TMath::Power(buffNormZ, 2));
372  normalPMT.at(1) =
373  buffNormY
374  / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
375  + TMath::Power(buffNormZ, 2));
376  normalPMT.at(2) =
377  buffNormZ
378  / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
379  + TMath::Power(buffNormZ, 2));
380  }
381 
382  CbmRichPoint* pmtPoint1 = (CbmRichPoint*) fPmtPoints->At(20);
383  scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0])
384  + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
385  + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
386  //cout << "1st scalar product between vectAM and normale = " << scalarProd << endl;
387  // 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.
388  normalCste =
389  -1
390  * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY()
391  + normalPMT.at(2) * pmtPoint1->GetZ());
392  CbmRichPoint* pmtPoint2 = (CbmRichPoint*) fPmtPoints->At(15);
393  scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0])
394  + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
395  + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
396  //cout << "2nd scalar product between vectAM and normale = " << scalarProd << endl;
397  CbmRichPoint* pmtPoint3 = (CbmRichPoint*) fPmtPoints->At(25);
398  scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0])
399  + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
400  + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
401  //cout << "3nd scalar product between vectAM and normale = " << scalarProd << endl;
402 }
403 
404 void CbmRichMirrorSortingAlignment::ComputeR2(vector<Double_t>& ptR2Center,
405  vector<Double_t>& ptR2Mirr,
406  vector<Double_t> ptM,
407  vector<Double_t> ptC,
408  vector<Double_t> ptR1,
409  TGeoNavigator* navi,
410  TString s) {
411  //cout << endl << "//------------------------------ CbmRichCorrection: ComputeR2 ------------------------------//" << endl << endl;
412 
413  vector<Double_t> normalMirr(3), ptCNew(3), ptTileCenter(3);
414  Double_t t1 = 0., t2 = 0., t3 = 0.;
415 
416  if (s == "Corrected") {
417  // Use the correction information from text file, to the tile sphere center:
418  // Reading misalignment information from correction_param.txt text file.
419  vector<Double_t> outputFit(4);
420  ifstream corrFile;
421  //TString str = fOutputDir + "correction_param_" + fNumbAxis + fTile + ".txt";
422  TString str = fOutputDir + "/correction_table/correction_param.txt";
423  corrFile.open(str);
424  if (corrFile.is_open()) {
425  for (Int_t i = 0; i < 4; i++) {
426  corrFile >> outputFit.at(i);
427  }
428  corrFile.close();
429  } else {
430  cout << "Error in CbmRichCorrection: unable to open parameter file!"
431  << endl;
432  cout << "Parameter file path = " << str << endl << endl;
433  sleep(5);
434  }
435  //cout << "Misalignment parameters read from file = [" << outputFit.at(0) << " ; " << outputFit.at(1) << " ; " << outputFit.at(2) << " ; " << outputFit.at(3) << "]" << endl;
436 
437  //ptCNew.at(0) = TMath::Abs(ptC.at(0) - TMath::Abs(outputFit.at(3)));
438  //ptCNew.at(1) = TMath::Abs(ptC.at(1) - TMath::Abs(outputFit.at(2)));
439  ptCNew.at(0) = ptC.at(0) + outputFit.at(3);
440  ptCNew.at(1) = ptC.at(1) + outputFit.at(2);
441  ptCNew.at(2) = ptC.at(2);
442  ptTileCenter.at(0) = navi->GetMotherMatrix()->GetTranslation()[0];
443  ptTileCenter.at(1) = navi->GetMotherMatrix()->GetTranslation()[1];
444  ptTileCenter.at(2) = navi->GetMotherMatrix()->GetTranslation()[2];
445  //cout << "Mirror tile center coordinates = {" << ptTileCenter.at(0) << ", " << ptTileCenter.at(1) << ", " << ptTileCenter.at(2) << "}" << endl;
446  Double_t x = 0., y = 0., z = 0., dist = 0., dist2 = 0., z2 = 0.;
447  x = TMath::Power(ptCNew.at(0) - ptTileCenter.at(0), 2);
448  y = TMath::Power(ptCNew.at(1) - ptTileCenter.at(1), 2);
449  z = TMath::Power(ptCNew.at(2) - ptTileCenter.at(2), 2);
450  dist = TMath::Sqrt(x + y + z);
451  z2 = ptTileCenter.at(2) - TMath::Sqrt(TMath::Power(300, 2) - x - y)
452  - ptCNew.at(2);
453  //cout << "{x, y, z} = {" << x << ", " << y << ", " << z << "}, dist = " << dist << " and z2 = " << z2 << endl;
454  dist2 = TMath::Sqrt(x + y + TMath::Power(z2 - ptTileCenter.at(2), 2));
455  //cout << "dist2 = " << dist2 << endl;
456  ptCNew.at(2) += z2;
457  //cout << "Sphere center coordinates of the rotated mirror tile, after correction, = {" << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
458  } else if (s == "Uncorrected") {
459  // Keep the same tile sphere center, with no correction information.
460  ptCNew = ptC;
461  //cout << "Sphere center coordinates of the rotated mirror tile, without correction = {" << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
462  } else {
463  //cout << "No input given in function ComputeR2! Uncorrected parameters for the sphere center of the tile will be used!" << endl;
464  ptCNew = ptC;
465  //cout << "Sphere center coordinates of the rotated mirror tile, without correction = {" << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}" << endl;
466  }
467 
468  normalMirr.at(0) = (ptCNew.at(0) - ptM.at(0))
469  / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
470  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
471  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
472  normalMirr.at(1) = (ptCNew.at(1) - ptM.at(1))
473  / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
474  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
475  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
476  normalMirr.at(2) = (ptCNew.at(2) - ptM.at(2))
477  / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
478  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
479  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
480  //cout << "Calculated and normalized normal of mirror tile = {" << normalMirr.at(0) << ", " << normalMirr.at(1) << ", " << normalMirr.at(2) << "}" << endl;
481 
482  t1 = ((ptR1.at(0) - ptM.at(0)) * (ptCNew.at(0) - ptM.at(0))
483  + (ptR1.at(1) - ptM.at(1)) * (ptCNew.at(1) - ptM.at(1))
484  + (ptR1.at(2) - ptM.at(2)) * (ptCNew.at(2) - ptM.at(2)))
485  / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
486  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
487  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
488  ptR2Center.at(0) =
489  2 * (ptM.at(0) + t1 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
490  ptR2Center.at(1) =
491  2 * (ptM.at(1) + t1 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
492  ptR2Center.at(2) =
493  2 * (ptM.at(2) + t1 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
494  t2 = ((ptR1.at(0) - ptCNew.at(0)) * (ptCNew.at(0) - ptM.at(0))
495  + (ptR1.at(1) - ptCNew.at(1)) * (ptCNew.at(1) - ptM.at(1))
496  + (ptR1.at(2) - ptCNew.at(2)) * (ptCNew.at(2) - ptM.at(2)))
497  / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
498  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
499  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
500  ptR2Mirr.at(0) =
501  2 * (ptCNew.at(0) + t2 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
502  ptR2Mirr.at(1) =
503  2 * (ptCNew.at(1) + t2 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
504  ptR2Mirr.at(2) =
505  2 * (ptCNew.at(2) + t2 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
506  /*//SAME AS calculation of t2 above
507  t3 = ((ptR1.at(0)-ptCNew.at(0))*(ptCNew.at(0)-ptM.at(0)) + (ptR1.at(1)-ptCNew.at(1))*(ptCNew.at(1)-ptM.at(1)) + (ptR1.at(2)-ptCNew.at(2))*(ptCNew.at(2)-ptM.at(2)))/TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0),2)+TMath::Power(ptCNew.at(1) - ptM.at(1),2)+TMath::Power(ptCNew.at(2) - ptM.at(2),2));
508  reflectedPtCooVectSphereUnity[0] = 2*(ptCNew.at(0)+t3*(normalMirr.at(0)))-ptR1.at(0);
509  reflectedPtCooVectSphereUnity[1] = 2*(ptCNew.at(1)+t3*(normalMirr.at(1)))-ptR1.at(1);
510  reflectedPtCooVectSphereUnity[2] = 2*(ptCNew.at(2)+t3*(normalMirr.at(2)))-ptR1.at(2);*/
511  //cout << "* using mirror point M to define \U00000394: {" << ptR2Center.at(0) << ", " << ptR2Center.at(1) << ", " << ptR2Center.at(2) << "}" << endl;
512  //cout << "Ref Pt Coo using unity Mirror-Sphere vector & sphere pt = {" << reflectedPtCooVectSphereUnity[0] << ", " << reflectedPtCooVectSphereUnity[1] << ", " << reflectedPtCooVectSphereUnity[2] << "}" << endl << endl;
513  //cout << "NofPMTPoints = " << NofPMTPoints << endl;
514 
515  //cout << "Coordinates of point R2 on reflective plane after reflection on the mirror tile:" << endl;
516  //cout << "* using sphere center C to define \U00000394: {" << ptR2Mirr.at(0) << ", " << ptR2Mirr.at(1) << ", " << ptR2Mirr.at(2) << "}" << endl;
517 }
518 
519 void CbmRichMirrorSortingAlignment::ComputeP(vector<Double_t>& ptPMirr,
520  vector<Double_t>& ptPR2,
521  vector<Double_t> normalPMT,
522  vector<Double_t> ptM,
523  vector<Double_t> ptR2Mirr,
524  Double_t constantePMT) {
525  //cout << endl << "//------------------------------ CbmRichCorrection: ComputeP ------------------------------//" << endl << endl;
526 
527  Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
528 
529  k1 = -1
530  * ((normalPMT.at(0) * ptM.at(0) + normalPMT.at(1) * ptM.at(1)
531  + normalPMT.at(2) * ptM.at(2) + constantePMT)
532  / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0))
533  + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
534  + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
535  ptPMirr.at(0) = ptM.at(0) + k1 * (ptR2Mirr.at(0) - ptM.at(0));
536  ptPMirr.at(1) = ptM.at(1) + k1 * (ptR2Mirr.at(1) - ptM.at(1));
537  ptPMirr.at(2) = ptM.at(2) + k1 * (ptR2Mirr.at(2) - ptM.at(2));
538  k2 = -1
539  * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1)
540  + normalPMT.at(2) * ptR2Mirr.at(2) + constantePMT)
541  / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0))
542  + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
543  + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
544  ptPR2.at(0) = ptR2Mirr.at(0) + k2 * (ptR2Mirr.at(0) - ptM.at(0));
545  ptPR2.at(1) = ptR2Mirr.at(1) + k2 * (ptR2Mirr.at(1) - ptM.at(1));
546  ptPR2.at(2) = ptR2Mirr.at(2) + k2 * (ptR2Mirr.at(2) - ptM.at(2));
547  //cout << "Coordinates of point P on PMT plane, after reflection on the mirror tile and extrapolation to the PMT plane:" << endl;
548  //cout << "* using mirror point M to define \U0001D49F ': {" << ptPMirr.at(0) << ", " << ptPMirr.at(1) << ", " << ptPMirr.at(2) << "}" << endl;
549  //cout << "* using reflected point R2 to define \U0001D49F ': {" << ptPR2.at(0) << ", " << ptPR2.at(1) << ", " << ptPR2.at(2) << "}" << endl;
550  checkCalc1 = ptPMirr.at(0) * normalPMT.at(0) + ptPMirr.at(1) * normalPMT.at(1)
551  + ptPMirr.at(2) * normalPMT.at(2) + constantePMT;
552  //cout << "Check whether extrapolated track point on PMT plane verifies its equation (value should be 0.):" << endl;
553  //cout << "* using mirror point M, checkCalc = " << checkCalc1 << endl;
554  checkCalc2 = ptPR2.at(0) * normalPMT.at(0) + ptPR2.at(1) * normalPMT.at(1)
555  + ptPR2.at(2) * normalPMT.at(2) + constantePMT;
556  //cout << "* using reflected point R2, checkCalc = " << checkCalc2 << endl;
557 }
558 
559 void CbmRichMirrorSortingAlignment::CreateHistoMap(
560  std::map<string, vector<CbmRichMirror*>> mirrorMap,
561  std::map<string, TH2D*>& histoMap) {
562  Int_t NofHits = 0;
563  Double_t phi = 0., theta0 = 0., thetaCh = 0.;
564 
565  for (std::map<string, vector<CbmRichMirror*>>::iterator it =
566  mirrorMap.begin();
567  it != mirrorMap.end();
568  ++it) {
569  // to get mirror Id:
570  string curMirrorId = it->first;
571  cout << "curMirrorId: '" << curMirrorId
572  << "' and vector size: " << it->second.size() << endl;
573  vector<CbmRichMirror*> mirror = it->second;
574  if (curMirrorId != "" && it->second.size() > fThreshold) {
575  histoMap[it->first] =
576  new TH2D(string("CherenkovHitsDistribReduced_" + it->first).c_str(),
577  "CherenkovHitsDistribReduced;#Phi_{Ch} "
578  "[rad];#theta_{Ch}-#theta_{0} [cm];Entries",
579  200,
580  -3.4,
581  3.4,
582  500,
583  -6.,
584  6.);
585  histoMap[it->first]->GetXaxis()->SetTitleSize(0.05);
586  histoMap[it->first]->GetXaxis()->SetTitleOffset(0.75);
587  histoMap[it->first]->GetYaxis()->SetTitleSize(0.04);
588  histoMap[it->first]->GetYaxis()->SetTitleOffset(1.2);
589  histoMap[it->first]->GetZaxis()->SetTitleSize(0.03);
590  histoMap[it->first]->GetZaxis()->SetTitleOffset(0.6);
591  for (int i = 0; i < it->second.size(); i++) {
592  CbmRichMirror* mirr = mirror.at(i);
593  string str = mirr->getMirrorId();
594  TVector3 mom = mirr->getMomentum();
595  vector<Double_t> projHit = mirr->getProjHit();
596  vector<Double_t> extrapHit = mirr->getExtrapHit();
597  CbmRichRingLight ringL = mirr->getRingLight();
598  // cout << "mirror: " << str << endl;
599  // cout << "momentum: {" << mom.X() << ", " << mom.Y() << ", " << mom.Z() << "}" << endl;
600  // cout << "Proj hit coo: {" << projHit[0] << ", " << projHit[1] << "}" << endl;
601  // cout << "Extrap hit coo: {" << extrapHit[0] << ", " << extrapHit[1] << "}" << endl;
602 
603  for (Int_t iH = 0; iH < ringL.GetNofHits(); iH++) {
604  // ----- Phi angle calculation ----- //
605  CbmRichHitLight hit = ringL.GetHit(iH);
606  //CbmRichHit* hit = static_cast<CbmRichHit*>(fPmtPoints->At(HitIndex));
607  phi = TMath::ATan2(hit.fY - ringL.GetCenterY(),
608  hit.fX - ringL.GetCenterX());
609 
610  // ----- Theta_Ch and Theta_0 calculations ----- //
611  thetaCh = sqrt(TMath::Power(projHit[0] - hit.fX, 2)
612  + TMath::Power(projHit[1] - hit.fY, 2));
613  theta0 = sqrt(TMath::Power(ringL.GetCenterX() - hit.fX, 2)
614  + TMath::Power(ringL.GetCenterY() - hit.fY, 2));
615  //cout << "Theta_0 = " << Theta_0 << endl;
616 
617  histoMap[it->first]->Fill(phi, thetaCh - theta0);
618  }
619  }
620  }
621  }
622  cout << endl;
623  // For loop to draw all histograms collected for all misaligned mirror tiles
624  // for (std::map<string, TH2D*>::iterator it=histoMap.begin(); it!=histoMap.end(); ++it) {
625  // cout << "Key str: " << it->first << " and nb of entries: " << it->second->GetEntries() << endl << endl;
626  // TCanvas* can = new TCanvas();
627  // it->second->Draw("colz");
628  // }
629 }
630 
631 void CbmRichMirrorSortingAlignment::DrawFitAndExtractAngles(
632  std::map<string, vector<Double_t>>& anglesMap,
633  std::map<string, TH2D*> histoMap) {
634  Int_t thresh = 3;
635  Double_t p1 = 0, p2 = 0, p3 = 0, chi2 = 0, focalLength = 150., q = 0., A = 0.,
636  alpha = 0., mis_x = 0., mis_y = 0.;
637  // mis_x && mis_y corresponds respect. to rotation angles around the Y and X axes.
638  // !!! BEWARE: AXES INDEXES ARE SWITCHED !!!
639 
640  for (std::map<string, TH2D*>::iterator it = histoMap.begin();
641  it != histoMap.end();
642  ++it) {
643  if (it->first != "") {
644  TCanvas* can = new TCanvas("can");
645  // can->Divide(3,1);
646  can->Divide(2, 1);
647  can->SetGrid();
648  gStyle->SetOptStat(0);
649  can->cd(1);
650  TH2D* histo = it->second;
651  for (Int_t y_bin = 1; y_bin <= 500; y_bin++) {
652  for (Int_t x_bin = 1; x_bin <= 200; x_bin++) {
653  if (histo->GetBinContent(x_bin, y_bin) < thresh) {
654  histo->SetBinContent(x_bin, y_bin, 0);
655  }
656  }
657  }
658  histo->Draw("colz");
659  histo->FitSlicesY(0, 0, -1, 1);
660  histo->Write();
661 
662  can->cd(2);
663  string histoName = "CherenkovHitsDistribReduced_" + it->first + "_1";
664  //cout << "HistoName: " << histoName << endl;
665  TH1D* histo_1 = (TH1D*) gDirectory->Get((histoName).c_str());
666  histo_1->GetXaxis()->SetTitle("#Phi_{Ch} [rad]");
667  histo_1->GetYaxis()->SetTitle("#theta_{Ch}-#theta_{0} [cm]");
668  histo_1->GetXaxis()->SetTitleSize(0.05);
669  histo_1->GetXaxis()->SetTitleOffset(0.75);
670  histo_1->GetYaxis()->SetTitleSize(0.04);
671  histo_1->GetYaxis()->SetTitleOffset(1.2);
672  histo_1->Draw("HIST");
673  histo_1->Write();
674  //
675  // if ( histo_1->GetSumOfWeights() == 0 || histo_1->Integral() == 0 ) {
676  // cout << "For mirror tile: " << it->first << ":" << endl;
677  // cout << "Sum of weights: " << histo_1->GetSumOfWeights() << endl;
678  // cout << "Integral: " << histo_1->Integral() << endl;
679  // continue;
680  // }
681  //
682  // can->cd(3);
683  TF1* f1 = new TF1("f1", "[2]+[0]*cos(x)+[1]*sin(x)", -3.5, 3.5);
684  f1->SetParameters(0, 0, 0);
685  f1->SetParNames("Delta_phi", "Delta_lambda", "Offset");
686  histo_1->Fit("f1", "", "");
687  TF1* fit = histo_1->GetFunction("f1");
688  p1 = fit->GetParameter("Delta_phi"),
689  p2 = fit->GetParameter("Delta_lambda"), p3 = fit->GetParameter("Offset"),
690  chi2 = fit->GetChisquare();
691  f1->SetParameters(fit->GetParameter(0), fit->GetParameter(1));
692  char leg[128];
693  f1->SetLineColor(2);
694  f1->Draw("same");
695  f1->Write();
696  // ------------------------------ CALCULATION OF MISALIGNMENT ANGLE ------------------------------ //
697  cout << setprecision(6) << endl;
698  q = TMath::ATan(fit->GetParameter(0) / fit->GetParameter(1));
699  //cout << endl << "fit_1 = " << fit->GetParameter(0) << " and fit_2 = " << fit->GetParameter(1) << endl;
700  //cout << "q = " << q << endl;
701  A = fit->GetParameter(1) / TMath::Cos(q);
702  //cout << "Parameter a = " << A << endl;
703  alpha =
704  TMath::ATan(A / 1.5) * 0.5
705  * TMath::Power(
706  10,
707  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
708  //cout << setprecision(6) << "Total angle of misalignment alpha = " << alpha << endl;
709  mis_x = TMath::ATan(fit->GetParameter(1) / focalLength) * 0.5
710  * TMath::Power(10, 3);
711  mis_y = TMath::ATan(fit->GetParameter(0) / focalLength) * 0.5
712  * TMath::Power(10, 3);
713  cout << "Horizontal displacement = " << mis_x
714  << " [mrad] and vertical displacement = " << mis_y << " [mrad]."
715  << endl
716  << endl;
717  TLegend* LEG = new TLegend(0.15, 0.25, 0.84, 0.4);
718  LEG->SetBorderSize(1);
719  LEG->SetFillColor(0);
720  LEG->SetMargin(0.2);
721  LEG->SetTextSize(0.04);
722  sprintf(leg, "Fitted sinusoid");
723  LEG->AddEntry(f1, leg, "l");
724  sprintf(leg, "Rotation angle around X = %.3f", -1 * mis_x);
725  LEG->AddEntry("", leg, "l");
726  sprintf(leg, "Rotation angle around Y = %.3f", mis_y);
727  LEG->AddEntry("", leg, "l");
728  sprintf(leg, "Offset = %.3f", fit->GetParameter(2));
729  LEG->AddEntry("", leg, "l");
730  LEG->Draw();
731  Cbm::SaveCanvasAsImage(
732  can, string(fOutputDir.Data() + fStudyName), "png");
733 
734  anglesMap[it->first].push_back(fit->GetParameter(1));
735  anglesMap[it->first].push_back(fit->GetParameter(0));
736  anglesMap[it->first].push_back(mis_x);
737  anglesMap[it->first].push_back(mis_y);
738  }
739  }
740 }
741 
742 void CbmRichMirrorSortingAlignment::Finish() {
743  std::map<string, TH2D*> histoMap;
744  histoMap.clear();
745  std::map<string, vector<Double_t>> anglesMap;
746  anglesMap.clear();
747 
748  // Filling the reduced thetaCh VS phi histogram and writing the resulting histogram in histoMap:
749  CreateHistoMap(fMirrorMap, histoMap);
750 
751  // Drawing the obtained thetaCh VS phi histogram ; fitting with sinusoid ; write calculated misalignment angles in anglesMap:
752  DrawFitAndExtractAngles(anglesMap, histoMap);
753 
754  //TString str_correction = fOutputDir + "/corr_params/correction_param_array_" + fStudyName + ".txt";
755  TString str_correction =
756  fOutputDir + "/correction_table/correction_param_array.txt";
757 
758  // Check if the directory exists, otherwise creates it.
759  TString pathToCorrectionTable = fOutputDir + "/correction_table/";
760  gSystem->mkdir(pathToCorrectionTable, true);
761 
762  ofstream corrFile;
763  corrFile.open(str_correction, std::ofstream::trunc);
764  if (corrFile.is_open()) {
765  corrFile.close();
766  } else {
767  cout << "Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
768  "parameter file!"
769  << endl;
770  }
771 
772  // Write correction values in output file:
773  for (std::map<string, vector<Double_t>>::iterator it = anglesMap.begin();
774  it != anglesMap.end();
775  ++it) {
776  string mirrorId = it->first;
777  cout << "curMirrorId: " << mirrorId << endl;
778  vector<Double_t> misAngles = it->second;
779  cout << "mirror correction parameters infos: {" << misAngles[0] << ", "
780  << misAngles[1] << "}" << endl;
781  corrFile.open(str_correction, std::ofstream::app);
782  if (corrFile.is_open()) {
783  corrFile << mirrorId << "\n";
784  corrFile << setprecision(7) << misAngles[0] << "\n";
785  corrFile << setprecision(7) << misAngles[1] << "\n";
786  corrFile.close();
787  cout << "Wrote correction parameters to: " << str_correction << endl;
788  } else {
789  cout << "Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
790  "parameter file!"
791  << endl;
792  }
793  }
794 
795  //TString str_angles = fOutputDir + "/corr_params/reconstructed_angles_array_" + fStudyName + ".txt";
796  TString str_angles =
797  fOutputDir + "/correction_table/reconstructed_angles_array.txt";
798  ofstream anglesFile;
799  anglesFile.open(str_angles, std::ofstream::trunc);
800  if (anglesFile.is_open()) {
801  anglesFile.close();
802  } else {
803  cout << "Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
804  "parameter file!"
805  << endl;
806  }
807 
808  // Write misalignment angles in output file:
809  for (std::map<string, vector<Double_t>>::iterator it = anglesMap.begin();
810  it != anglesMap.end();
811  ++it) {
812  string mirrorId = it->first;
813  cout << "curMirrorId: " << mirrorId << endl;
814  vector<Double_t> misAngles = it->second;
815  cout << "mirror reconstructed angles infos: {" << misAngles[2] << ", "
816  << misAngles[3] << "}" << endl;
817  anglesFile.open(str_angles, std::ofstream::app);
818  if (anglesFile.is_open()) {
819  anglesFile << mirrorId << "\n";
820  anglesFile << setprecision(7) << misAngles[2] << "\n";
821  anglesFile << setprecision(7) << misAngles[3] << "\n";
822  anglesFile.close();
823  cout << "Wrote reconstructed angles to: " << str_angles << endl;
824  } else {
825  cout << "Error in CbmRichMirrorSortingAlignment::Finish ; unable to open "
826  "parameter file!"
827  << endl;
828  }
829  }
830 }
831 ClassImp(CbmRichMirrorSortingAlignment)
CbmRichPoint.h
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition: CbmMCTrack.h:172
CbmRichMirror
Definition: CbmRichMirror.h:20
CbmRichMirrorSortingAlignment::fMirrorPoints
TClonesArray * fMirrorPoints
Definition: CbmRichMirrorSortingAlignment.h:94
CbmTrackMatchNew::GetNofWrongHits
Int_t GetNofWrongHits() const
Definition: CbmTrackMatchNew.h:33
CbmRichMirror::getMirrorId
string getMirrorId()
Definition: CbmRichMirror.h:53
CbmRichMirrorSortingAlignment.h
CbmRichMirror::getRingLight
CbmRichRingLight getRingLight()
Definition: CbmRichMirror.h:58
CbmRichMirrorSortingAlignment::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmRichMirrorSortingAlignment.h:99
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
CbmRichMirrorSortingAlignment::CbmRichMirrorSortingAlignment
CbmRichMirrorSortingAlignment()
Definition: CbmRichMirrorSortingAlignment.cxx:36
CbmRichMirrorSortingAlignment::fRichRings
TClonesArray * fRichRings
Definition: CbmRichMirrorSortingAlignment.h:92
CbmRichMirrorSortingAlignment::fThreshold
Int_t fThreshold
Definition: CbmRichMirrorSortingAlignment.h:88
CbmRichRingFitterEllipseTau.h
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmRichMirrorSortingAlignment::fTauFit
CbmRichRingFitterEllipseTau * fTauFit
Definition: CbmRichMirrorSortingAlignment.h:85
CbmRichMirror::getProjHit
vector< Double_t > getProjHit()
Definition: CbmRichMirror.h:55
CbmGlobalTrack::GetRichRingIndex
Int_t GetRichRingIndex() const
Definition: CbmGlobalTrack.h:41
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
CbmRichMirror::setRingLight
void setRingLight(CbmRichRingLight ringL)
Definition: CbmRichMirror.h:50
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichMirrorSortingAlignment::CreateHistoMap
void CreateHistoMap(std::map< string, vector< CbmRichMirror * >> mirrorMap, std::map< string, TH2D * > &histoMap)
Definition: CbmRichMirrorSortingAlignment.cxx:559
CbmGlobalTrack.h
CbmRichMirrorSortingAlignment::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition: CbmRichMirrorSortingAlignment.h:100
CbmRichRing
Definition: CbmRichRing.h:17
CbmRichMirrorSortingAlignment::fPmtPoints
TClonesArray * fPmtPoints
Definition: CbmRichMirrorSortingAlignment.h:96
CbmRichHitLight::fX
float fX
Definition: CbmRichRingLight.h:34
CbmRichRing.h
CbmRichMirrorSortingAlignment::fTrackParams
TClonesArray * fTrackParams
Definition: CbmRichMirrorSortingAlignment.h:98
CbmRichMirror::setExtrapHit
void setExtrapHit(Double_t xX, Double_t yY)
Definition: CbmRichMirror.h:44
CbmRichMirrorSortingAlignment::ComputeR2
void ComputeR2(vector< Double_t > &ptR2Center, vector< Double_t > &ptR2Mirr, vector< Double_t > ptM, vector< Double_t > ptC, vector< Double_t > ptR1, TGeoNavigator *navi, TString s)
Definition: CbmRichMirrorSortingAlignment.cxx:404
CbmRichNavigationUtil::FindIntersection
static string FindIntersection(const FairTrackParam *par, TVector3 &crossPoint, const string &volumeName)
Definition: CbmRichNavigationUtil.h:14
CbmRichRingFitterCOP
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
Definition: CbmRichRingFitterCOP.h:28
CbmRichGeoManager.h
CbmRichMirrorSortingAlignment::fOutputDir
TString fOutputDir
Definition: CbmRichMirrorSortingAlignment.h:86
CbmRichMirrorSortingAlignment::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichMirrorSortingAlignment.cxx:119
CbmRichRingLight::GetNofHits
int GetNofHits() const
Return number of hits in ring.
Definition: CbmRichRingLight.h:108
CbmRichMirrorSortingAlignment::fMCTracks
TClonesArray * fMCTracks
Definition: CbmRichMirrorSortingAlignment.h:93
CbmRichMirrorSortingAlignment::fMirrorMap
std::map< string, vector< CbmRichMirror * > > fMirrorMap
Definition: CbmRichMirrorSortingAlignment.h:89
CbmGlobalTrack::GetStsTrackIndex
Int_t GetStsTrackIndex() const
Definition: CbmGlobalTrack.h:38
CbmRichMirrorSortingAlignment::fRichProjections
TClonesArray * fRichProjections
Definition: CbmRichMirrorSortingAlignment.h:97
z2
Double_t z2[nSects2]
Definition: pipe_v16a_mvdsts100.h:11
CbmRichConverter.h
Convert internal data classes to cbmroot common data classes.
CbmRichRingLight::GetCenterY
float GetCenterY() const
Definition: CbmRichRingLight.h:160
CbmRichMirror::setMirrorId
void setMirrorId(string s)
Definition: CbmRichMirror.h:39
CbmTrackMatchNew.h
CbmRichMirrorSortingAlignment::~CbmRichMirrorSortingAlignment
virtual ~CbmRichMirrorSortingAlignment()
Definition: CbmRichMirrorSortingAlignment.cxx:56
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:407
CbmRichMirrorSortingAlignment::fGlobalTracks
TClonesArray * fGlobalTracks
Definition: CbmRichMirrorSortingAlignment.h:91
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmRichMirror::setProjHit
void setProjHit(Double_t xX, Double_t yY)
Definition: CbmRichMirror.h:41
CbmRichMirrorSortingAlignment::GetPmtNormal
void GetPmtNormal(Int_t NofPMTPoints, vector< Double_t > &normalPMT, Double_t &normalCste)
Definition: CbmRichMirrorSortingAlignment.cxx:298
CbmRichHitLight
Definition: CbmRichRingLight.h:14
CbmUtils.h
CbmRichRingFitterCOP.h
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
CbmRichMirrorSortingAlignment::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichMirrorSortingAlignment.cxx:58
CbmRichConverter::Init
static void Init()
Initialize array of RICH hits.
Definition: CbmRichConverter.h:93
CbmRichMirrorSortingAlignment::fEventNb
UInt_t fEventNb
Definition: CbmRichMirrorSortingAlignment.h:83
CbmGlobalTrack
Definition: CbmGlobalTrack.h:26
CbmRichMirrorSortingAlignment
Definition: CbmRichMirrorSortingAlignment.h:21
CbmRichRingFitterCOP::DoFit
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition: CbmRichRingFitterCOP.cxx:19
CbmRichMirrorSortingAlignment::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichMirrorSortingAlignment.cxx:742
CbmMCTrack.h
CbmTrackMatchNew::GetNofTrueHits
Int_t GetNofTrueHits() const
Definition: CbmTrackMatchNew.h:32
CbmRichHitLight::fY
float fY
Definition: CbmRichRingLight.h:35
CbmMCTrack
Definition: CbmMCTrack.h:34
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichMirror::getExtrapHit
vector< Double_t > getExtrapHit()
Definition: CbmRichMirror.h:56
CbmRichConverter::CopyHitsToRingLight
static void CopyHitsToRingLight(const CbmRichRing *ring1, CbmRichRingLight *ring2)
Copy hits from CbmRichRing to CbmRichRingLight.
Definition: CbmRichConverter.h:41
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichRingLight::GetHit
CbmRichHitLight GetHit(int ind)
Return hit by the index.
Definition: CbmRichRingLight.h:114
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmRichMirror::getMomentum
TVector3 getMomentum()
Definition: CbmRichMirror.h:54
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmRichMirrorSortingAlignment::fCopFit
CbmRichRingFitterCOP * fCopFit
Definition: CbmRichMirrorSortingAlignment.h:84
Cbm::SaveCanvasAsImage
void SaveCanvasAsImage(TCanvas *c, const std::string &dir, const std::string &option)
Definition: CbmUtils.cxx:20
CbmRichRingLight::GetCenterX
float GetCenterX() const
Definition: CbmRichRingLight.h:159
CbmRichMirrorSortingAlignment::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 constantePMT)
Definition: CbmRichMirrorSortingAlignment.cxx:519
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmRichNavigationUtil.h
CbmRichRingLight
Definition: CbmRichRingLight.h:39
CbmRichMirror::setMomentum
void setMomentum(TVector3 v)
Definition: CbmRichMirror.h:40
CbmRichMirrorSortingAlignment::fRefPlanePoints
TClonesArray * fRefPlanePoints
Definition: CbmRichMirrorSortingAlignment.h:95
CbmRichMirrorSortingAlignment::fStudyName
TString fStudyName
Definition: CbmRichMirrorSortingAlignment.h:87
CbmRichMirrorSortingAlignment::DrawFitAndExtractAngles
void DrawFitAndExtractAngles(std::map< string, vector< Double_t >> &anglesMap, std::map< string, TH2D * > histoMap)
Definition: CbmRichMirrorSortingAlignment.cxx:631