CbmRoot
alignment/CbmRichProjectionProducerAnalytical.cxx
Go to the documentation of this file.
1 
8 #include "CbmRichProjectionProducerAnalytical.h"
9 
10 #include "CbmMCTrack.h"
11 #include "FairGeoNode.h"
12 #include "FairGeoTransform.h"
13 #include "FairGeoVector.h"
14 #include "FairRootManager.h"
15 #include "FairRun.h"
16 #include "FairRunAna.h"
17 #include "FairRuntimeDb.h"
18 #include "FairTrackParam.h"
19 #include "TGeoManager.h"
20 
21 #include "TClonesArray.h"
22 #include "TMatrixFSym.h"
23 #include "TVector3.h"
24 
25 #include "CbmRichGeoManager.h"
26 #include "CbmRichNavigationUtil.h"
27 #include "CbmRichPoint.h"
28 
29 #include "FairLogger.h"
30 
31 #include <cmath>
32 #include <iostream>
33 
34 using std::cout;
35 using std::endl;
36 
37 #include "../alignment/CbmRichMirrorMisalignmentCorrectionUtils.h"
38 
39 
40 CbmRichProjectionProducerAnalytical::CbmRichProjectionProducerAnalytical()
41  : fTrackParams(NULL), fNHits(0), fEventNum(0) {}
42 
43 CbmRichProjectionProducerAnalytical::~CbmRichProjectionProducerAnalytical() {
44  FairRootManager* fManager = FairRootManager::Instance();
45  fManager->Write();
46 }
47 
48 
49 void CbmRichProjectionProducerAnalytical::Init() {
50  LOG(info) << "CbmRichProjectionProducerAnalytical::Init()";
51  FairRootManager* manager = FairRootManager::Instance();
52 
53  fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
54  if (NULL == fTrackParams) {
55  Fatal("CbmRichProjectionProducerAnalytical::Init:",
56  "No RichTrackParamZ array!");
57  }
58 
59  fMirrorCorrectionParameterFile =
60  new CbmRichMirrorMisalignmentCorrectionUtils();
61  fMirrorCorrectionParameterFile->Init(fPathToMirrorCorrectionParameterFile);
62 }
63 
64 void CbmRichProjectionProducerAnalytical::DoProjection(TClonesArray* richProj) {
65  fEventNum++;
66  cout << "CbmRichProjectionProducerAnalytical:: event " << fEventNum << endl;
67 
68  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
69  Double_t mirrorX = gp->fMirrorX;
70  Double_t mirrorY = gp->fMirrorY;
71  Double_t mirrorZ = gp->fMirrorZ;
72  Double_t mirrorR = gp->fMirrorR;
73 
74  richProj->Delete();
75  TMatrixFSym covMat(5);
76  for (Int_t i = 0; i < 5; i++) {
77  for (Int_t j = 0; j <= i; j++) {
78  covMat(i, j) = 0;
79  }
80  }
81  covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
82  1.e-4;
83 
84  for (Int_t j = 0; j < fTrackParams->GetEntriesFast(); j++) {
85  FairTrackParam* point = (FairTrackParam*) fTrackParams->At(j);
86  new ((*richProj)[j]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
87 
88  // check if Array was filled
89  if (point->GetX() == 0 && point->GetY() == 0 && point->GetZ() == 0
90  && point->GetTx() == 0 && point->GetTy() == 0)
91  continue;
92  if (point->GetQp() == 0) continue;
93 
94  Double_t rho1 = 0.;
95  TVector3 startP, momP, crossP, centerP;
96 
97 
98  Double_t p = 1. / TMath::Abs(point->GetQp());
99  Double_t pz;
100  if ((1 + point->GetTx() * point->GetTx() + point->GetTy() * point->GetTy())
101  > 0.)
102  pz = p
103  / TMath::Sqrt(1 + point->GetTx() * point->GetTx()
104  + point->GetTy() * point->GetTy());
105  else {
106  cout << " -E- RichProjectionProducer: strange value for calculating pz: "
107  << (1 + point->GetTx() * point->GetTx()
108  + point->GetTy() * point->GetTy())
109  << endl;
110  pz = 0.;
111  continue;
112  }
113  Double_t px = pz * point->GetTx();
114  Double_t py = pz * point->GetTy();
115  momP.SetXYZ(px, py, pz);
116  point->Position(startP);
117  if ((mirrorY * startP.y()) < 0)
118  mirrorY =
119  -mirrorY; // check that mirror center and startP are in same hemisphere
120 
121  // calculation of intersection of track with selected mirror
122  // corresponds to calculation of intersection between a straight line and a sphere:
123  // vector: r = startP - mirrorCenter
124  // RxP = r*momP
125  // normP2 = momP^2
126  // dist = r^2 - fR^2
127  // -> rho1 = (-RxP+sqrt(RxP^2-normP2*dist))/normP2 extrapolation factor for:
128  // intersection point crossP = startP + rho1 * momP
129  Double_t RxP =
130  (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
131  + momP.z() * (startP.z() - mirrorZ));
132  Double_t normP2 =
133  (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
134  Double_t dist =
135  (startP.x() * startP.x() + mirrorX * mirrorX + startP.y() * startP.y()
136  + mirrorY * mirrorY + startP.z() * startP.z() + mirrorZ * mirrorZ
137  - 2 * startP.x() * mirrorX - 2 * startP.y() * mirrorY
138  - 2 * startP.z() * mirrorZ - mirrorR * mirrorR);
139 
140  if ((RxP * RxP - normP2 * dist) > 0.) {
141  if (normP2 != 0.)
142  rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
143  if (normP2 == 0)
144  cout << " Error in track extrapolation: momentum = 0 " << endl;
145  } else {
146  //cout << " -E- RichProjectionProducer: RxP*RxP-normP2*dist = " << RxP*RxP-normP2*dist << endl;
147  }
148 
149  Double_t crossPx = startP.x() + rho1 * momP.x();
150  Double_t crossPy = startP.y() + rho1 * momP.y();
151  Double_t crossPz = startP.z() + rho1 * momP.z();
152  crossP.SetXYZ(crossPx, crossPy, crossPz);
153 
154  // check if crosspoint with mirror and chosen mirrorcenter (y) are in same hemisphere
155  // if not recalculate crossing point
156  if ((mirrorY * crossP.y()) < 0) {
157  mirrorY = -mirrorY;
158  RxP =
159  (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
160  + momP.z() * (startP.z() - mirrorZ));
161  normP2 =
162  (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
163  dist =
164  (startP.x() * startP.x() + mirrorX * mirrorX + startP.y() * startP.y()
165  + mirrorY * mirrorY + startP.z() * startP.z() + mirrorZ * mirrorZ
166  - 2 * startP.x() * mirrorX - 2 * startP.y() * mirrorY
167  - 2 * startP.z() * mirrorZ - mirrorR * mirrorR);
168 
169  if ((RxP * RxP - normP2 * dist) > 0.) {
170  if (normP2 != 0.)
171  rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
172  if (normP2 == 0)
173  cout << " Error in track extrapolation: momentum = 0 " << endl;
174  } else {
175  //cout << " -E- RichProjectionProducer: RxP*RxP-normP2*dist = " << RxP*RxP-normP2*dist << endl;
176  }
177 
178  crossPx = startP.x() + rho1 * momP.x();
179  crossPy = startP.y() + rho1 * momP.y();
180  crossPz = startP.z() + rho1 * momP.z();
181  crossP.SetXYZ(crossPx, crossPy, crossPz);
182  }
183 
184  centerP.SetXYZ(mirrorX, mirrorY, mirrorZ); // mirror center
185 
186  TVector3 crossPoint;
187  string volumeName = CbmRichNavigationUtil::FindIntersection(
188  point, crossPoint, "mirror_tile_");
189  //cout << "CbmRichProjectionAnalytical: volume Name = " << volumeName << endl;
190 
191  TVector3 newCenterP;
192  newCenterP = MirrorCenter(centerP, volumeName);
193 
194  // calculate normal on crosspoint with mirror
195  TVector3 normP(crossP.x() - newCenterP.x(),
196  crossP.y() - newCenterP.y(),
197  crossP.z() - newCenterP.z());
198  // TVector3 normP(crossP.x()-centerP.x(),crossP.y()-centerP.y(),crossP.z()-centerP.z());
199  normP = normP.Unit();
200  // check that normal has same z-direction as momentum
201  if ((normP.z() * momP.z()) < 0.)
202  normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
203 
204  // reflect track
205  Double_t np =
206  normP.x() * momP.x() + normP.y() * momP.y() + normP.z() * momP.z();
207 
208  TVector3 ref;
209  ref.SetXYZ(2 * np * normP.x() - momP.x(),
210  2 * np * normP.y() - momP.y(),
211  2 * np * normP.z() - momP.z());
212 
213  if (ref.z() != 0.) {
214 
215  TVector3 inPos(0., 0., 0.), outPos(0., 0., 0.);
216 
217  if (gp->fGeometryType == CbmRichGeometryTypeCylindrical) {
218  GetPmtIntersectionPointCyl(&newCenterP, &crossP, &ref, &inPos);
219  // cout << "inPos new center P: " << inPos.X() << ", " << inPos.Y() << ", " << inPos.Z() << endl;
220  // GetPmtIntersectionPointCyl(&centerP, &crossP, &ref, &inPos);
221  // Transform intersection point in same way as MCPoints were transformed in HitProducer before stored as Hit
222  // For the cylindrical geometry we also pass the path to the strip block
223  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
224 
225  //cout << "inPoint:" << inPos.X() << " " << inPos.Y() << " " << inPos.Z() << " outPoint:" << outPos.X() << " " << outPos.Y() << " " << outPos.Z() << endl;
226 
227 
228  } else if (gp->fGeometryType == CbmRichGeometryTypeTwoWings) {
229  GetPmtIntersectionPointTwoWings(&centerP, &crossP, &ref, &inPos);
230  // Transform intersection point in same way as MCPoints were transformed in HitProducer before stored as Hit
231  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
232  } else {
233  Fatal("CbmRichProjectionProducerAnalytical ", "unknown geometry type");
234  }
235 
236  Bool_t isInsidePmt =
237  CbmRichGeoManager::GetInstance().IsPointInsidePmt(&outPos);
238 
239  if (isInsidePmt) {
240  FairTrackParam richtrack(
241  outPos.x(), outPos.y(), outPos.z(), 0., 0., 0., covMat);
242  *(FairTrackParam*) (richProj->At(j)) = richtrack;
243  }
244  } // if (refZ!=0.)
245  } // j
246 }
247 
248 TVector3
249 CbmRichProjectionProducerAnalytical::MirrorCenter(const TVector3 centerP,
250  const string volumeName) {
251  if (volumeName == ""
252  || !(fMirrorCorrectionParameterFile->GetMirrorCorrectionParamBool())) {
253  //cout << "center P: " << centerP.x() << ", " << centerP.y() << ", " << centerP.z() << endl;
254  return (centerP);
255  }
256 
257  TVector3 newCenterP;
258 
259  string mirrorID = GetMirrorID(volumeName);
260 
261  std::map<string, std::pair<Double_t, Double_t>> mirrorCorrectionParamMap;
262  mirrorCorrectionParamMap =
263  fMirrorCorrectionParameterFile->GetMirrorCorrectionParamMap();
264  for (std::map<string, std::pair<Double_t, Double_t>>::iterator it =
265  mirrorCorrectionParamMap.begin();
266  it != mirrorCorrectionParamMap.end();
267  ++it) {
268  if (it->first == mirrorID) {
269  cout << "Map's first key: " << it->first
270  << " and pair: " << it->second.first << ", " << it->second.second
271  << endl;
272  newCenterP.SetX(centerP.X() + it->second.first);
273  newCenterP.SetY(centerP.Y() + it->second.second);
274  newCenterP.SetZ(centerP.Z());
275  }
276  }
277  //cout << "center P: " << centerP.x() << ", " << centerP.y() << ", " << centerP.z() << endl;
278  //cout << "new center P: " << newCenterP.x() << ", " << newCenterP.y() << ", " << newCenterP.z() << endl;
279 
280  return (newCenterP);
281 }
282 
283 /*
284 TVector3 CbmRichProjectionProducerAnalytical::MirrorCenter(
285  const TVector3 centerP,
286  const string volumeName)
287 {
288  if ( volumeName == "" ) {
289  return(centerP);
290  }
291 
292  TVector3 newCenterP;
293 
294  string mirrorID = GetMirrorID(volumeName);
295 
296  Int_t lineCounter=1, lineIndex=0;
297  //TString corrFileName = "/data/Sim_Outputs/Correction_test/new_code/corr_params/correction_param_array_Full_Correction.txt";
298  TString corrFileName = fPathToMirrorCorrectionParameterFile;
299  //TString corrFileName = fTxtFile + "correction_param_array_" + fStudyName + ".txt";
300  string fileLine = "", strMisX = "", strMisY = "";
301  Double_t misX=0., misY=0.;
302  ifstream corrFile;
303  corrFile.open(corrFileName);
304 
305  if (corrFile.is_open())
306  {
307  while (!corrFile.eof())
308  {
309  getline(corrFile, fileLine);
310  lineIndex = fileLine.find(mirrorID, 0);
311  if (lineIndex != string::npos)
312  {
313  cout << mirrorID << " has been found in the file at line: " << lineCounter << " and position: " << lineIndex << "." << endl;
314  break;
315  }
316  lineCounter++;
317  }
318  corrFile >> misY;
319  //cout << "number at line: " << lineCounter+2 << " = " << misX << "." << endl;
320  corrFile >> misX;
321  //cout << "number at line: " << lineCounter+3 << " = " << misY << "." << endl;
322 
323  corrFile.close();
324  }
325  else {
326  cout << "CbmRichProjectionProducerAnalytical:: unable to open file with correction parameters!" << endl;
327  cout << "Parameter file path located at: " << corrFileName << endl;
328  cout << "Using ideal tile center." << endl;
329  return(centerP);
330  }
331 
332  cout << "Correction parameters = " << misX << ", " << misY << endl;
333  newCenterP.SetX(centerP.X() + misX);
334  newCenterP.SetY(centerP.Y() + misY);
335  newCenterP.SetZ(centerP.Z());
336 
337  cout << "center P: " << centerP.x() << ", " << centerP.y() << ", " << centerP.z() << endl;
338  cout << "new center P: " << newCenterP.x() << ", " << newCenterP.y() << ", " << newCenterP.z() << endl;
339 
340  return(newCenterP);
341 }
342 */
343 
344 string
345 CbmRichProjectionProducerAnalytical::GetMirrorID(const string volumeName) {
346  string str1 = "mirror_tile_", str2 = "";
347  std::size_t found = volumeName.find(str1);
348  if (found != std::string::npos) {
349  //cout << "first 'mirror_tile_type' found at: " << found << '\n';
350  Int_t end = str1.length() + 3;
351  str2 = volumeName.substr(found, end);
352  }
353  cout << "Mirror ID: " << str2 << endl;
354 
355  return (str2);
356 }
357 
358 void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(
359  const TVector3* centerP,
360  const TVector3* crossP,
361  const TVector3* ref,
362  TVector3* outPoint) {
363  // crosspoint whith photodetector plane:
364  // calculate intersection between straight line and (tilted) plane:
365  // normal on plane tilted by theta around x-axis: (0,-sin(theta),cos(theta)) = n
366  // normal on plane tilted by phi around y-axis: (-sin(phi),0,cos(phi)) = n
367  // normal on plane tilted by theta around x-axis and phi around y-axis: (-sin(phi),-sin(theta)cos(phi),cos(theta)cos(phi)) = n
368  // point on plane is (fDetX,fDetY,fDetZ) = p as photodetector is tiled around its center
369  // equation of plane for r being point in plane: n(r-p) = 0
370  // calculate intersection point of reflected track with plane: r=intersection point
371  // intersection point = crossP + rho2 * refl_track
372  // take care for all 4 cases:
373  // -> first calculate for case x>0, then check
374 
375  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
376 
377  if (gp->fGeometryType != CbmRichGeometryTypeTwoWings) {
378  Fatal(
379  "CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings ",
380  "Wrong geometry, this method could be used only for "
381  "CbmRichGeometryTypeTwoWings");
382  }
383 
384  Double_t pmtPhi = gp->fPmt.fPhi;
385  Double_t pmtTheta = gp->fPmt.fTheta;
386  Double_t pmtPlaneX = gp->fPmt.fPlaneX;
387  Double_t pmtPlaneY = gp->fPmt.fPlaneY;
388  Double_t pmtPlaneZ = gp->fPmt.fPlaneZ;
389  Double_t rho2 = 0.;
390 
391  if (centerP->y() > 0) {
392  rho2 =
393  (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
394  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
395  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
396  / (-TMath::Sin(pmtPhi) * ref->x()
397  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
398  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
399  }
400  if (centerP->y() < 0) {
401  rho2 =
402  (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
403  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * (-pmtPlaneY - crossP->y())
404  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
405  / (-TMath::Sin(pmtPhi) * ref->x()
406  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
407  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
408  }
409 
410  //rho2 = -1*(crossP.z() - fDetZ)/refZ; // only for theta = 0, phi=0
411  Double_t xX = crossP->x() + ref->x() * rho2;
412  Double_t yY = crossP->y() + ref->y() * rho2;
413  Double_t zZ = crossP->z() + ref->z() * rho2;
414 
415  if (xX < 0) {
416  if (centerP->y() > 0) {
417  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
418  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi)
419  * (pmtPlaneY - crossP->y())
420  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi)
421  * (pmtPlaneZ - crossP->z()))
422  / (-TMath::Sin(-pmtPhi) * ref->x()
423  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
424  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
425  }
426  if (centerP->y() < 0) {
427  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
428  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi)
429  * (-pmtPlaneY - crossP->y())
430  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi)
431  * (pmtPlaneZ - crossP->z()))
432  / (-TMath::Sin(-pmtPhi) * ref->x()
433  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
434  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
435  }
436 
437  xX = crossP->x() + ref->x() * rho2;
438  yY = crossP->y() + ref->y() * rho2;
439  zZ = crossP->z() + ref->z() * rho2;
440  }
441 
442  outPoint->SetXYZ(xX, yY, zZ);
443 }
444 
445 
446 void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(
447  const TVector3* centerP,
448  const TVector3* crossP,
449  const TVector3* ref,
450  TVector3* outPoint) {
451  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
452  if (gp->fGeometryType != CbmRichGeometryTypeCylindrical) {
453  Fatal("CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl ",
454  "Wrong geometry, this method could be used only for "
455  "CbmRichGeometryTypeCylindrical");
456  }
457  Double_t xX;
458  Double_t yY;
459  Double_t zZ;
460  Double_t maxDist = 0.;
461 
462  for (map<string, CbmRichRecGeoParPmt>::iterator it = gp->fPmtMap.begin();
463  it != gp->fPmtMap.end();
464  it++) {
465  Double_t pmtPlaneX = it->second.fPlaneX;
466  Double_t pmtPlaneY = it->second.fPlaneY;
467  Double_t pmtPlaneZ = it->second.fPlaneZ;
468  Double_t pmtPhi = it->second.fPhi;
469  Double_t pmtTheta = it->second.fTheta;
470 
471  if (!(pmtPlaneX >= 0 && pmtPlaneY >= 0)) continue;
472 
473  double rho2 = 0.;
474 
475  if (centerP->y() > 0) {
476  rho2 =
477  (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
478  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
479  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi)
480  * (pmtPlaneZ - crossP->z()))
481  / (-TMath::Sin(pmtPhi) * ref->x()
482  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
483  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
484  }
485  if (centerP->y() < 0) {
486  rho2 = (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
487  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi)
488  * (-pmtPlaneY - crossP->y())
489  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi)
490  * (pmtPlaneZ - crossP->z()))
491  / (-TMath::Sin(pmtPhi) * ref->x()
492  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
493  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
494  }
495 
496  Double_t cxX = crossP->x() + ref->x() * rho2;
497  Double_t cyY = crossP->y() + ref->y() * rho2;
498  Double_t czZ = crossP->z() + ref->z() * rho2;
499 
500  if (cxX < 0) {
501  if (centerP->y() > 0) {
502  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
503  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi)
504  * (pmtPlaneY - crossP->y())
505  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi)
506  * (pmtPlaneZ - crossP->z()))
507  / (-TMath::Sin(-pmtPhi) * ref->x()
508  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
509  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
510  }
511  if (centerP->y() < 0) {
512  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
513  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi)
514  * (-pmtPlaneY - crossP->y())
515  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi)
516  * (pmtPlaneZ - crossP->z()))
517  / (-TMath::Sin(-pmtPhi) * ref->x()
518  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
519  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
520  }
521 
522  cxX = crossP->x() + ref->x() * rho2;
523  cyY = crossP->y() + ref->y() * rho2;
524  czZ = crossP->z() + ref->z() * rho2;
525  }
526 
527  Double_t dP = TMath::Sqrt((crossP->x() - cxX) * (crossP->X() - cxX)
528  + (crossP->y() - cyY) * (crossP->y() - cyY)
529  + (crossP->z() - czZ) * (crossP->Z() - czZ));
530 
531  if (dP >= maxDist) {
532  maxDist = dP;
533  xX = cxX;
534  yY = cyY;
535  zZ = czZ;
536  }
537  }
538 
539  outPoint->SetXYZ(xX, yY, zZ);
540 }
CbmRichPoint.h
CbmRichProjectionProducerBase::fPathToMirrorCorrectionParameterFile
string fPathToMirrorCorrectionParameterFile
Definition: alignment/CbmRichProjectionProducerBase.h:59
CbmRichRecGeoPar::fMirrorZ
Double_t fMirrorZ
Definition: CbmRichRecGeoPar.h:232
CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl
void GetPmtIntersectionPointCyl(const TVector3 *centerP, const TVector3 *crossP, const TVector3 *ref, TVector3 *outPoint)
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:446
CbmRichGeometryTypeTwoWings
@ CbmRichGeometryTypeTwoWings
Definition: CbmRichRecGeoPar.h:24
CbmRichRecGeoParPmt::fPhi
Double_t fPhi
Definition: CbmRichRecGeoPar.h:58
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
CbmRichProjectionProducerAnalytical::CbmRichProjectionProducerAnalytical
CbmRichProjectionProducerAnalytical()
Standard constructor.
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:40
CbmRichRecGeoParPmt::fPlaneY
Double_t fPlaneY
Definition: CbmRichRecGeoPar.h:67
CbmRichMirrorMisalignmentCorrectionUtils::Init
void Init(const string &s)
Initialization in case one needs to initialize some TCloneArrays.
Definition: CbmRichMirrorMisalignmentCorrectionUtils.h:37
CbmRichRecGeoPar::fPmtMap
std::map< std::string, CbmRichRecGeoParPmt > fPmtMap
Definition: CbmRichRecGeoPar.h:223
CbmRichRecGeoPar::fPmt
CbmRichRecGeoParPmt fPmt
Definition: CbmRichRecGeoPar.h:218
CbmRichRecGeoPar::fMirrorY
Double_t fMirrorY
Definition: CbmRichRecGeoPar.h:231
CbmRichProjectionProducerBase::fMirrorCorrectionParameterFile
CbmRichMirrorMisalignmentCorrectionUtils * fMirrorCorrectionParameterFile
Definition: alignment/CbmRichProjectionProducerBase.h:60
CbmRichProjectionProducerAnalytical::~CbmRichProjectionProducerAnalytical
virtual ~CbmRichProjectionProducerAnalytical()
Destructor.
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:43
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRecGeoPar::fMirrorX
Double_t fMirrorX
Definition: CbmRichRecGeoPar.h:230
CbmRichProjectionProducerAnalytical::GetMirrorID
string GetMirrorID(const string volumeName)
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:345
CbmRichRecGeoPar::fGeometryType
CbmRichGeometryType fGeometryType
Definition: CbmRichRecGeoPar.h:220
CbmRichNavigationUtil::FindIntersection
static string FindIntersection(const FairTrackParam *par, TVector3 &crossPoint, const string &volumeName)
Definition: CbmRichNavigationUtil.h:14
CbmRichRecGeoParPmt::fTheta
Double_t fTheta
Definition: CbmRichRecGeoPar.h:57
CbmRichGeoManager.h
CbmRichProjectionProducerAnalytical::Init
virtual void Init()
Initialization of the task.
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:49
CbmRichMirrorMisalignmentCorrectionUtils::GetMirrorCorrectionParamMap
std::map< string, std::pair< Double_t, Double_t > > GetMirrorCorrectionParamMap()
Definition: CbmRichMirrorMisalignmentCorrectionUtils.h:104
CbmRichRecGeoParPmt::fPlaneX
Double_t fPlaneX
Definition: CbmRichRecGeoPar.h:66
CbmRichRecGeoPar
PMT parameters for the RICH geometry.
Definition: CbmRichRecGeoPar.h:87
CbmRichMirrorMisalignmentCorrectionUtils
class checks correction parameter file containing mirror misalignment information.
Definition: CbmRichMirrorMisalignmentCorrectionUtils.h:20
CbmRichGeoManager::IsPointInsidePmt
Bool_t IsPointInsidePmt(const TVector3 *rotatedPoint)
Definition: CbmRichGeoManager.cxx:584
CbmRichGeometryTypeCylindrical
@ CbmRichGeometryTypeCylindrical
Definition: CbmRichRecGeoPar.h:25
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:407
CbmRichProjectionProducerAnalytical.h
Project track by straight line from imaginary plane to the mirror and reflect it to the photodetector...
CbmRichProjectionProducerAnalytical::fTrackParams
TClonesArray * fTrackParams
Definition: alignment/CbmRichProjectionProducerAnalytical.h:78
CbmRichProjectionProducerAnalytical::fEventNum
int fEventNum
Definition: alignment/CbmRichProjectionProducerAnalytical.h:81
richProj
TClonesArray * richProj
Definition: Compute_distance.h:18
CbmMCTrack.h
CbmRichMirrorMisalignmentCorrectionUtils::GetMirrorCorrectionParamBool
bool GetMirrorCorrectionParamBool()
Definition: CbmRichMirrorMisalignmentCorrectionUtils.h:108
CbmRichRecGeoPar::fMirrorR
Double_t fMirrorR
Definition: CbmRichRecGeoPar.h:233
CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings
void GetPmtIntersectionPointTwoWings(const TVector3 *centerP, const TVector3 *crossP, const TVector3 *ref, TVector3 *outPoint)
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:358
CbmRichProjectionProducerAnalytical::MirrorCenter
TVector3 MirrorCenter(const TVector3 centerP, const string volumeName)
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:249
CbmRichProjectionProducerAnalytical::DoProjection
virtual void DoProjection(TClonesArray *richProj)
Execute task.
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:64
CbmRichGeoManager::fGP
CbmRichRecGeoPar * fGP
Definition: CbmRichGeoManager.h:23
CbmRichNavigationUtil.h
CbmRichRecGeoParPmt::fPlaneZ
Double_t fPlaneZ
Definition: CbmRichRecGeoPar.h:68