CbmRoot
CbmRichProjectionProducer2.cxx
Go to the documentation of this file.
1 // ---------- Original Headers ---------- //
2 #include "CbmRichProjectionProducer2.h"
3 #include "FairLogger.h"
4 #include "FairRootManager.h"
5 
6 #include "CbmRichHit.h"
7 #include "TClonesArray.h"
8 
9 // ---------- Included Headers ---------- //
10 #include "CbmRichPoint.h"
11 #include "TGeoManager.h"
12 
13 #include "CbmMCTrack.h"
14 #include "CbmRichGeoManager.h"
15 #include "CbmRichHitProducer.h"
16 #include "CbmTrackMatchNew.h"
17 #include "CbmUtils.h"
18 #include "FairTrackParam.h"
19 #include "TMatrixFSym.h"
20 
21 #include "CbmRichNavigationUtil2.h"
22 
23 class TGeoNode;
24 class TGeoVolume;
25 
26 CbmRichProjectionProducer2::CbmRichProjectionProducer2()
27  : fTrackParams(NULL)
28  , fMCTracks(NULL)
29  , fRichPoints(NULL)
30  , fOutputDir("")
31  , fNumbAxis("")
32  , fTile("")
33  , fEventNum(0) {}
34 
35 CbmRichProjectionProducer2::~CbmRichProjectionProducer2() {
36  FairRootManager* fmanager = FairRootManager::Instance();
37  fmanager->Write();
38 }
39 
40 void CbmRichProjectionProducer2::Init() {
41  LOG(info) << "CbmRichProjectionProducerAnalytical::Init()";
42  FairRootManager* manager = FairRootManager::Instance();
43 
44  fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
45  if (NULL == fTrackParams) {
46  Fatal("CbmRichProjectionProducerAnalytical::Init",
47  "No RichTrackParamZ array !");
48  }
49 
50  fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
51  if (NULL == fMCTracks) {
52  Fatal("CbmRichCorrectionVector::Init", "No MCTracks array !");
53  }
54 
55  fRichPoints = (TClonesArray*) manager->GetObject("RichPoint");
56  if (NULL == fRichPoints) {
57  Fatal("CbmRichCorrectionVector::Init", "No RichPoint array !");
58  }
59 }
60 
61 void CbmRichProjectionProducer2::DoProjection(TClonesArray* projectedPoint) {
62  cout << endl
63  << "//"
64  "--------------------------------------------------------------------"
65  "------------------------------------------//"
66  << endl;
67  cout << "//-------------------------------- CbmRichProjectionProducer2: Do "
68  "Projection ---------------------------------//"
69  << endl
70  << endl;
71  cout << "//"
72  "--------------------------------------------------------------------"
73  "--------------------------------------------------//"
74  << endl;
75 
76  fEventNum++;
77  //LOG(debug2) << "CbmRichProjectionProducer2 : Event #" << fEventNum;
78  cout << "CbmRichProjectionProducer2 : Event #" << fEventNum << endl;
79 
80  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
81  Double_t pmtPlaneX = gp->fPmt.fPlaneX;
82  Double_t pmtPlaneY = gp->fPmt.fPlaneY;
83  Double_t pmtWidth = gp->fPmt.fWidth;
84  Double_t pmtHeight = gp->fPmt.fHeight;
85 
86  projectedPoint->Delete();
87  TMatrixFSym covMat(5);
88  for (Int_t iMatrix = 0; iMatrix < 5; iMatrix++) {
89  for (Int_t jMatrix = 0; jMatrix <= iMatrix; jMatrix++) {
90  covMat(iMatrix, jMatrix) = 0;
91  }
92  }
93  covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
94  1.e-4;
95 
96  for (Int_t j = 0; j < fTrackParams->GetEntriesFast(); j++) {
97  FairTrackParam* point = (FairTrackParam*) fTrackParams->At(j);
98  new ((*projectedPoint)[j]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
99 
100  // check if Array was filled
101  if (point->GetX() == 0 && point->GetY() == 0 && point->GetZ() == 0
102  && point->GetTx() == 0 && point->GetTy() == 0)
103  continue;
104  if (point->GetQp() == 0) continue;
105 
106  Double_t xDet = 0., yDet = 0., zDet = 0.;
107  Double_t* pmtPt;
108  pmtPt = ProjectionProducer(point);
109  xDet = pmtPt[0];
110  yDet = pmtPt[1];
111  zDet = pmtPt[2];
112 
113  //check that crosspoint inside the plane
114  Double_t marginX = 2.; // [cm]
115  Double_t marginY = 2.; // [cm]
116  // upper pmt planes
117  Double_t pmtYTop = TMath::Abs(pmtPlaneY) + pmtHeight + marginY;
118  Double_t pmtYBottom = TMath::Abs(pmtPlaneY) - pmtHeight - marginY;
119  Double_t absYDet = TMath::Abs(yDet);
120  Bool_t isYOk = (absYDet <= pmtYTop && absYDet >= pmtYBottom);
121 
122  Double_t pmtXMin = -TMath::Abs(pmtPlaneX) - pmtWidth - marginX;
123  Double_t pmtXMax = TMath::Abs(pmtPlaneX) + pmtWidth + marginX;
125  Bool_t isXOk = (xDet >= pmtXMin && xDet <= pmtXMax);
126 
127  if (isYOk && isXOk) {
128  FairTrackParam richtrack(xDet, yDet, zDet, 0., 0., 0., covMat);
129  *(FairTrackParam*) (projectedPoint->At(j)) = richtrack;
130  }
131  }
132 }
133 
134 Double_t*
135 CbmRichProjectionProducer2::ProjectionProducer(FairTrackParam* trackParam) {
136  cout << "//------------------------------ CbmRichProjectionProducer2: "
137  "Projection Producer ------------------------------//"
138  << endl
139  << endl;
140 
141  static Double_t pmtPt[3];
142 
143  Int_t NofPMTPoints = fRichPoints->GetEntriesFast();
144 
145  // Declaration of points coordinates.
146  Double_t sphereRadius = 300., constantePMT = 0.;
147  vector<Double_t> ptM(3), ptC(3), ptCNew(3), momR1(3), ptR1(3), normalPMT(3),
148  ptR2Mirr(3), ptR2Center(3), ptPMirr(3), ptPR2(3), ptTileCenter(3);
149  vector<Double_t> ptCIdeal(3), ptR2CenterUnCorr(3), ptR2CenterIdeal(3),
150  ptR2MirrUnCorr(3), ptR2MirrIdeal(3), ptPMirrUnCorr(3), ptPMirrIdeal(3),
151  ptPR2UnCorr(3), ptPR2Ideal(3);
152  TVector3 outPos, outPosUnCorr, outPosIdeal;
153  // Declaration of ring parameters.
154  Double_t ringCenter[] = {0., 0., 0.}, distToExtrapTrackHit = 0.,
155  distToExtrapTrackHitInPlane = 0.;
156  //Declarations related to geometry.
157  Int_t pmtTrackID = -1, pmtMotherID = -100;
158  TGeoNavigator* navi;
159 
160  ptR1.at(0) = trackParam->GetX();
161  ptR1.at(1) = trackParam->GetY();
162  ptR1.at(2) = trackParam->GetZ();
163  Double_t nx = 0., ny = 0., nz = 0.;
164  CbmRichNavigationUtil2::GetDirCos(trackParam, nx, ny, nz);
165  TVector3 dirCos;
166  dirCos.SetXYZ(nx, ny, nz);
167  GetPmtNormal(NofPMTPoints, normalPMT, constantePMT);
168  cout << "Calculated normal vector to PMT plane = {" << normalPMT.at(0) << ", "
169  << normalPMT.at(1) << ", " << normalPMT.at(2)
170  << "} and constante d = " << constantePMT << endl
171  << endl;
172 
173  TVector3 mirrorPoint;
174  TString mirrorIntersection1 = CbmRichNavigationUtil2::FindIntersection(
175  trackParam, mirrorPoint, "mirror_tile_type", navi);
176  ptM.at(0) = mirrorPoint.x();
177  ptM.at(1) = mirrorPoint.y();
178  ptM.at(2) = mirrorPoint.z();
179  cout << "mirrorIntersection1: " << mirrorIntersection1 << endl;
180  cout << "Mirror point coordinates = {" << mirrorPoint.x() << ", "
181  << mirrorPoint.y() << ", " << mirrorPoint.z() << "}" << endl;
182  TString mirrorIntersection =
183  "/cave_1/rich1_0/richContainer_333/rich_gas_329/mirror_full_half_321/"
184  "mirror_tile_2_0_148/mirror_tile_type2"
185  "_inter_108/mirror_tile_type2_94/"
186  + mirrorIntersection1;
187  cout << "mirrorIntersection: " << mirrorIntersection << endl;
188 
189  if (mirrorIntersection1) {
190  //navi->cd(mirrorIntersection1);
191  //navi = gGeoManager->GetCurrentNavigator();
192  navi->Dump();
193  cout << "Navigator path: " << navi->GetPath() << endl;
194  cout << "Coordinates of sphere center: " << endl;
195  navi->GetCurrentMatrix()->Print();
196  ptCIdeal.at(0) = navi->GetCurrentMatrix()->GetTranslation()[0];
197  ptCIdeal.at(1) = navi->GetCurrentMatrix()->GetTranslation()[1];
198  ptCIdeal.at(2) = navi->GetCurrentMatrix()->GetTranslation()[2];
199  cout << "Coordinates of tile center: " << endl;
200  navi->GetMotherMatrix()->Print();
201  ptC.at(0) = 0., ptC.at(1) = 132.594000,
202  ptC.at(2) = 54.267226; // Theoretical coordinates of point C.
203  cout << "Sphere center coordinates of the aligned mirror tile, ideal = {"
204  << ptCIdeal.at(0) << ", " << ptCIdeal.at(1) << ", " << ptCIdeal.at(2)
205  << "}" << endl;
206  cout << "Sphere center coordinates of the rotated mirror tile, w/ "
207  "GeoManager, = {"
208  << ptC.at(0) << ", " << ptC.at(1) << ", " << ptC.at(2)
209  << "} and sphere inner radius = " << sphereRadius << endl
210  << endl;
211  //ptCNew = RotateSphereCenter(ptTileCenter, ptC, navi);
212 
213  cout << "FairTrackParam Point coordinates = {" << ptR1.at(0) << ", "
214  << ptR1.at(1) << ", " << ptR1.at(2) << "}" << endl;
215  cout << "And FairTrackParam Point direction cosines = {" << ptR1.at(3)
216  << ", " << ptR1.at(4) << ", " << ptR1.at(5) << "}" << endl;
217  cout << "Mirror Point coordinates = {" << ptM.at(0) << ", " << ptM.at(1)
218  << ", " << ptM.at(2) << "}" << endl;
219 
220  ComputeR2(
221  ptR2CenterUnCorr, ptR2MirrUnCorr, ptM, ptC, ptR1, navi, "Uncorrected");
222  ComputeR2(ptR2Center, ptR2Mirr, ptM, ptC, ptR1, navi, "Corrected");
223  ComputeR2(
224  ptR2CenterIdeal, ptR2MirrIdeal, ptM, ptCIdeal, ptR1, navi, "Uncorrected");
225 
226  ComputeP(
227  ptPMirrUnCorr, ptPR2UnCorr, normalPMT, ptM, ptR2MirrUnCorr, constantePMT);
228  ComputeP(ptPMirr, ptPR2, normalPMT, ptM, ptR2Mirr, constantePMT);
229  ComputeP(
230  ptPMirrIdeal, ptPR2Ideal, normalPMT, ptM, ptR2MirrIdeal, constantePMT);
231 
232  TVector3 inPosUnCorr(
233  ptPMirrUnCorr.at(0), ptPMirrUnCorr.at(1), ptPMirrUnCorr.at(2));
234  CbmRichGeoManager::GetInstance().RotatePoint(&inPosUnCorr, &outPosUnCorr);
235  cout << "New mirror points coordinates = {" << outPosUnCorr.x() << ", "
236  << outPosUnCorr.y() << ", " << outPosUnCorr.z() << "}" << endl;
237  TVector3 inPos(ptPMirr.at(0), ptPMirr.at(1), ptPMirr.at(2));
238  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
239  cout << "New mirror points coordinates = {" << outPos.x() << ", "
240  << outPos.y() << ", " << outPos.z() << "}" << endl;
241  TVector3 inPosIdeal(
242  ptPMirrIdeal.at(0), ptPMirrIdeal.at(1), ptPMirrIdeal.at(2));
243  CbmRichGeoManager::GetInstance().RotatePoint(&inPosIdeal, &outPosIdeal);
244  cout << endl
245  << "New mirror points coordinates = {" << outPosIdeal.x() << ", "
246  << outPosIdeal.y() << ", " << outPosIdeal.z() << "}" << endl
247  << endl;
248  } else {
249  cout << "No mirror intersection found ..." << endl;
250  outPos.SetXYZ(0, 0, 0);
251  }
252 
253  pmtPt[0] = outPos.x();
254  pmtPt[1] = outPos.y();
255  pmtPt[2] = outPos.z();
256  return pmtPt;
257 }
258 
259 void CbmRichProjectionProducer2::GetPmtNormal(Int_t NofPMTPoints,
260  vector<Double_t>& normalPMT,
261  Double_t& normalCste) {
262  //cout << endl << "//------------------------------ CbmRichProjectionProducer2: Calculate PMT Normal ------------------------------//" << endl << endl;
263 
264  Int_t pmtTrackID, pmtMotherID;
265  Double_t buffNormX = 0., buffNormY = 0., buffNormZ = 0., k = 0.,
266  scalarProd = 0.;
267  Double_t pmtPt[] = {0., 0., 0.};
268  Double_t a[] = {0., 0., 0.}, b[] = {0., 0., 0.}, c[] = {0., 0., 0.};
269  CbmMCTrack* track;
270 
271  /*
272  * 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.
273  * Formula used is: vect(AB) x vect(AC) = normal.
274  * Normalize the normal vector obtained and check with three random points from the PMT plane, whether the scalar product is equal to zero.
275  */
276  for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
277  CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
278  pmtTrackID = pmtPoint->GetTrackID();
279  track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
280  pmtMotherID = track->GetMotherId();
281  //a[0] = pmtPoint->GetX(), a[1] = pmtPoint->GetY(), a[2] = pmtPoint->GetZ();
282  cout << "a[0] = " << a[0] << ", a[1] = " << a[1] << " et a[2] = " << a[2]
283  << endl;
284  break;
285  }
286  for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
287  CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
288  pmtTrackID = pmtPoint->GetTrackID();
289  track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
290  pmtMotherID = track->GetMotherId();
291  //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
292  if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2)
293  + TMath::Power(a[1] - pmtPoint->GetY(), 2)
294  + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
295  > 7) {
296  //b[0] = pmtPoint->GetX(), b[1] = pmtPoint->GetY(), b[2] = pmtPoint->GetZ();
297  cout << "b[0] = " << b[0] << ", b[1] = " << b[1] << " et b[2] = " << b[2]
298  << endl;
299  break;
300  }
301  }
302  for (Int_t iPmt = 0; iPmt < NofPMTPoints; iPmt++) {
303  CbmRichPoint* pmtPoint = (CbmRichPoint*) fRichPoints->At(iPmt);
304  pmtTrackID = pmtPoint->GetTrackID();
305  track = (CbmMCTrack*) fMCTracks->At(pmtTrackID);
306  pmtMotherID = track->GetMotherId();
307  //cout << "PMT Point coordinates; x = " << pmtPoint->GetX() << ", y = " << pmtPoint->GetY() << " and z = " << pmtPoint->GetZ() << endl;
308  if (TMath::Sqrt(TMath::Power(a[0] - pmtPoint->GetX(), 2)
309  + TMath::Power(a[1] - pmtPoint->GetY(), 2)
310  + TMath::Power(a[2] - pmtPoint->GetZ(), 2))
311  > 7
312  && TMath::Sqrt(TMath::Power(b[0] - pmtPoint->GetX(), 2)
313  + TMath::Power(b[1] - pmtPoint->GetY(), 2)
314  + TMath::Power(b[2] - pmtPoint->GetZ(), 2))
315  > 7) {
316  c[0] = pmtPoint->GetX(), c[1] = pmtPoint->GetY(), c[2] = pmtPoint->GetZ();
317  //cout << "c[0] = " << c[0] << ", c[1] = " << c[1] << " et c[2] = " << c[2] << endl;
318  break;
319  }
320  }
321 
322  k = (b[0] - a[0]) / (c[0] - a[0]);
323  if ((b[1] - a[1]) - (k * (c[1] - a[1])) == 0
324  || (b[2] - a[2]) - (k * (c[2] - a[2])) == 0) {
325  cout << "Error in normal calculation, vect_AB and vect_AC are collinear."
326  << endl;
327  } else {
328  buffNormX = (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]);
329  buffNormY = (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]);
330  buffNormZ = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]);
331  normalPMT.at(0) =
332  buffNormX
333  / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
334  + TMath::Power(buffNormZ, 2));
335  normalPMT.at(1) =
336  buffNormY
337  / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
338  + TMath::Power(buffNormZ, 2));
339  normalPMT.at(2) =
340  buffNormZ
341  / TMath::Sqrt(TMath::Power(buffNormX, 2) + TMath::Power(buffNormY, 2)
342  + TMath::Power(buffNormZ, 2));
343  }
344 
345  CbmRichPoint* pmtPoint1 = (CbmRichPoint*) fRichPoints->At(20);
346  scalarProd = normalPMT.at(0) * (pmtPoint1->GetX() - a[0])
347  + normalPMT.at(1) * (pmtPoint1->GetY() - a[1])
348  + normalPMT.at(2) * (pmtPoint1->GetZ() - a[2]);
349  //cout << "1st scalar product between vectAM and normale = " << scalarProd << endl;
350  // 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.
351  normalCste =
352  -1
353  * (normalPMT.at(0) * pmtPoint1->GetX() + normalPMT.at(1) * pmtPoint1->GetY()
354  + normalPMT.at(2) * pmtPoint1->GetZ());
355  CbmRichPoint* pmtPoint2 = (CbmRichPoint*) fRichPoints->At(15);
356  scalarProd = normalPMT.at(0) * (pmtPoint2->GetX() - a[0])
357  + normalPMT.at(1) * (pmtPoint2->GetY() - a[1])
358  + normalPMT.at(2) * (pmtPoint2->GetZ() - a[2]);
359  //cout << "2nd scalar product between vectAM and normale = " << scalarProd << endl;
360  CbmRichPoint* pmtPoint3 = (CbmRichPoint*) fRichPoints->At(25);
361  scalarProd = normalPMT.at(0) * (pmtPoint3->GetX() - a[0])
362  + normalPMT.at(1) * (pmtPoint3->GetY() - a[1])
363  + normalPMT.at(2) * (pmtPoint3->GetZ() - a[2]);
364  //cout << "3nd scalar product between vectAM and normale = " << scalarProd << endl;
365 }
366 
367 void CbmRichProjectionProducer2::ComputeR2(vector<Double_t>& ptR2Center,
368  vector<Double_t>& ptR2Mirr,
369  vector<Double_t> ptM,
370  vector<Double_t> ptC,
371  vector<Double_t> ptR1,
372  TGeoNavigator* navi,
373  TString s) {
374  cout << endl
375  << "//------------------------------ CbmRichCorrection: ComputeR2 "
376  "------------------------------//"
377  << endl
378  << endl;
379 
380  vector<Double_t> normalMirr(3), ptCNew(3), ptTileCenter(3);
381  Double_t t1 = 0., t2 = 0., t3 = 0.;
382 
383  if (s == "Corrected") {
384  // Use the correction information from text file, to the tile sphere center:
385  // Reading misalignment information from correction_param.txt text file.
386  vector<Double_t> outputFit(4);
387  ifstream corr_file;
388  TString str = fOutputDir + "correction_param_" + fNumbAxis + fTile + ".txt";
389  corr_file.open(str);
390  if (corr_file.is_open()) {
391  for (Int_t i = 0; i < 4; i++) {
392  corr_file >> outputFit.at(i);
393  }
394  corr_file.close();
395  } else {
396  cout << "Error in CbmRichCorrection: unable to open parameter file!"
397  << endl
398  << endl;
399  sleep(5);
400  }
401  cout << "Misalignment parameters read from file = [" << outputFit.at(0)
402  << " ; " << outputFit.at(1) << " ; " << outputFit.at(2) << " ; "
403  << outputFit.at(3) << "]" << endl;
404 
405  ptCNew.at(0) = TMath::Abs(ptC.at(0) - TMath::Abs(outputFit.at(3)));
406  ptCNew.at(1) = TMath::Abs(ptC.at(1) - TMath::Abs(outputFit.at(2)));
407  //ptCNew.at(0) = ptC.at(0) - outputFit.at(3);
408  //ptCNew.at(1) = ptC.at(1) - outputFit.at(2);
409  ptCNew.at(2) = ptC.at(2);
410  ptTileCenter.at(0) = navi->GetMotherMatrix()->GetTranslation()[0];
411  ptTileCenter.at(1) = navi->GetMotherMatrix()->GetTranslation()[1];
412  ptTileCenter.at(2) = navi->GetMotherMatrix()->GetTranslation()[2];
413  cout << "Mirror tile center coordinates = {" << ptTileCenter.at(0) << ", "
414  << ptTileCenter.at(1) << ", " << ptTileCenter.at(2) << "}" << endl;
415  Double_t x = 0., y = 0., z = 0., dist = 0., dist2 = 0., z2 = 0.;
416  x = TMath::Power(ptCNew.at(0) - ptTileCenter.at(0), 2);
417  y = TMath::Power(ptCNew.at(1) - ptTileCenter.at(1), 2);
418  z = TMath::Power(ptCNew.at(2) - ptTileCenter.at(2), 2);
419  dist = TMath::Sqrt(x + y + z);
420  z2 = ptTileCenter.at(2) - TMath::Sqrt(TMath::Power(300, 2) - x - y)
421  - ptCNew.at(2);
422  cout << "{x, y, z} = {" << x << ", " << y << ", " << z
423  << "}, dist = " << dist << " and z2 = " << z2 << endl;
424  dist2 = TMath::Sqrt(x + y + TMath::Power(z2 - ptTileCenter.at(2), 2));
425  cout << "dist2 = " << dist2 << endl;
426  ptCNew.at(2) += z2;
427  cout << "Sphere center coordinates of the rotated mirror tile, after "
428  "correction, = {"
429  << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}"
430  << endl;
431  } else if (s == "Uncorrected") {
432  // Keep the same tile sphere center, with no correction information.
433  ptCNew = ptC;
434  cout << "Sphere center coordinates of the rotated mirror tile, without "
435  "correction = {"
436  << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}"
437  << endl;
438  } else {
439  cout << "No input given in function ComputeR2! Uncorrected parameters for "
440  "the sphere center of the tile will be used!"
441  << endl;
442  ptCNew = ptC;
443  cout << "Sphere center coordinates of the rotated mirror tile, without "
444  "correction = {"
445  << ptCNew.at(0) << ", " << ptCNew.at(1) << ", " << ptCNew.at(2) << "}"
446  << endl;
447  }
448 
449  normalMirr.at(0) = (ptCNew.at(0) - ptM.at(0))
450  / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
451  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
452  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
453  normalMirr.at(1) = (ptCNew.at(1) - ptM.at(1))
454  / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
455  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
456  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
457  normalMirr.at(2) = (ptCNew.at(2) - ptM.at(2))
458  / TMath::Sqrt(TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
459  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
460  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
461  //cout << "Calculated and normalized normal of mirror tile = {" << normalMirr.at(0) << ", " << normalMirr.at(1) << ", " << normalMirr.at(2) << "}" << endl;
462 
463  t1 = ((ptR1.at(0) - ptM.at(0)) * (ptCNew.at(0) - ptM.at(0))
464  + (ptR1.at(1) - ptM.at(1)) * (ptCNew.at(1) - ptM.at(1))
465  + (ptR1.at(2) - ptM.at(2)) * (ptCNew.at(2) - ptM.at(2)))
466  / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
467  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
468  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
469  ptR2Center.at(0) =
470  2 * (ptM.at(0) + t1 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
471  ptR2Center.at(1) =
472  2 * (ptM.at(1) + t1 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
473  ptR2Center.at(2) =
474  2 * (ptM.at(2) + t1 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
475  t2 = ((ptR1.at(0) - ptCNew.at(0)) * (ptCNew.at(0) - ptM.at(0))
476  + (ptR1.at(1) - ptCNew.at(1)) * (ptCNew.at(1) - ptM.at(1))
477  + (ptR1.at(2) - ptCNew.at(2)) * (ptCNew.at(2) - ptM.at(2)))
478  / (TMath::Power(ptCNew.at(0) - ptM.at(0), 2)
479  + TMath::Power(ptCNew.at(1) - ptM.at(1), 2)
480  + TMath::Power(ptCNew.at(2) - ptM.at(2), 2));
481  ptR2Mirr.at(0) =
482  2 * (ptCNew.at(0) + t2 * (ptCNew.at(0) - ptM.at(0))) - ptR1.at(0);
483  ptR2Mirr.at(1) =
484  2 * (ptCNew.at(1) + t2 * (ptCNew.at(1) - ptM.at(1))) - ptR1.at(1);
485  ptR2Mirr.at(2) =
486  2 * (ptCNew.at(2) + t2 * (ptCNew.at(2) - ptM.at(2))) - ptR1.at(2);
487 
488  cout << "Coordinates of point R2 on reflective plane after reflection on the "
489  "mirror tile:"
490  << endl;
491  //cout << "* using mirror point M to define \U00000394: {" << ptR2Center.at(0) << ", " << ptR2Center.at(1) << ", " << ptR2Center.at(2) << "}" << endl;
492  cout << "* using sphere center C to define \U00000394: {" << ptR2Mirr.at(0)
493  << ", " << ptR2Mirr.at(1) << ", " << ptR2Mirr.at(2) << "}" << endl;
494 }
495 
496 void CbmRichProjectionProducer2::ComputeP(vector<Double_t>& ptPMirr,
497  vector<Double_t>& ptPR2,
498  vector<Double_t> normalPMT,
499  vector<Double_t> ptM,
500  vector<Double_t> ptR2Mirr,
501  Double_t constantePMT) {
502  cout << endl
503  << "//------------------------------ CbmRichCorrection: ComputeP "
504  "------------------------------//"
505  << endl
506  << endl;
507 
508  Double_t k1 = 0., k2 = 0., checkCalc1 = 0., checkCalc2 = 0.;
509 
510  k1 = -1
511  * ((normalPMT.at(0) * ptM.at(0) + normalPMT.at(1) * ptM.at(1)
512  + normalPMT.at(2) * ptM.at(2) + constantePMT)
513  / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0))
514  + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
515  + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
516  ptPMirr.at(0) = ptM.at(0) + k1 * (ptR2Mirr.at(0) - ptM.at(0));
517  ptPMirr.at(1) = ptM.at(1) + k1 * (ptR2Mirr.at(1) - ptM.at(1));
518  ptPMirr.at(2) = ptM.at(2) + k1 * (ptR2Mirr.at(2) - ptM.at(2));
519  k2 = -1
520  * ((normalPMT.at(0) * ptR2Mirr.at(0) + normalPMT.at(1) * ptR2Mirr.at(1)
521  + normalPMT.at(2) * ptR2Mirr.at(2) + constantePMT)
522  / (normalPMT.at(0) * (ptR2Mirr.at(0) - ptM.at(0))
523  + normalPMT.at(1) * (ptR2Mirr.at(1) - ptM.at(1))
524  + normalPMT.at(2) * (ptR2Mirr.at(2) - ptM.at(2))));
525  ptPR2.at(0) = ptR2Mirr.at(0) + k2 * (ptR2Mirr.at(0) - ptM.at(0));
526  ptPR2.at(1) = ptR2Mirr.at(1) + k2 * (ptR2Mirr.at(1) - ptM.at(1));
527  ptPR2.at(2) = ptR2Mirr.at(2) + k2 * (ptR2Mirr.at(2) - ptM.at(2));
528  cout << "Coordinates of point P on PMT plane, after reflection on the mirror "
529  "tile and extrapolation to the PMT plane:"
530  << endl;
531  cout << "* using mirror point M to define \U0001D49F ': {" << ptPMirr.at(0)
532  << ", " << ptPMirr.at(1) << ", " << ptPMirr.at(2) << "}" << endl;
533  //cout << "* using reflected point R2 to define \U0001D49F ': {" << ptPR2.at(0) << ", " << ptPR2.at(1) << ", " << ptPR2.at(2) << "}" << endl;
534  checkCalc1 = ptPMirr.at(0) * normalPMT.at(0) + ptPMirr.at(1) * normalPMT.at(1)
535  + ptPMirr.at(2) * normalPMT.at(2) + constantePMT;
536  cout << "Check whether extrapolated track point on PMT plane verifies its "
537  "equation (value should be 0.):"
538  << endl;
539  cout << "* using mirror point M, checkCalc = " << checkCalc1 << endl;
540  checkCalc2 = ptPR2.at(0) * normalPMT.at(0) + ptPR2.at(1) * normalPMT.at(1)
541  + ptPR2.at(2) * normalPMT.at(2) + constantePMT;
542  //cout << "* using reflected point R2, checkCalc = " << checkCalc2 << endl;
543 }
544 
545 ClassImp(CbmRichProjectionProducer2)
CbmRichPoint.h
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmRichRecGeoParPmt::fHeight
Double_t fHeight
Definition: CbmRichRecGeoPar.h:71
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
CbmRichRecGeoParPmt::fPlaneY
Double_t fPlaneY
Definition: CbmRichRecGeoPar.h:67
CbmRichProjectionProducer2
Definition: CbmRichProjectionProducer2.h:18
CbmRichRecGeoPar::fPmt
CbmRichRecGeoParPmt fPmt
Definition: CbmRichRecGeoPar.h:218
CbmRichProjectionProducer2::~CbmRichProjectionProducer2
virtual ~CbmRichProjectionProducer2()
Definition: CbmRichProjectionProducer2.cxx:35
CbmRichProjectionProducer2::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: CbmRichProjectionProducer2.cxx:367
CbmRichProjectionProducer2::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: CbmRichProjectionProducer2.cxx:496
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichProjectionProducer2::fMCTracks
TClonesArray * fMCTracks
Definition: CbmRichProjectionProducer2.h:92
CbmRichProjectionProducer2::fOutputDir
TString fOutputDir
Definition: CbmRichProjectionProducer2.h:97
CbmRichProjectionProducer2::DoProjection
virtual void DoProjection(TClonesArray *projectedPoint)
Definition: CbmRichProjectionProducer2.cxx:61
CbmRichProjectionProducer2.h
CbmRichProjectionProducer2::GetPmtNormal
void GetPmtNormal(Int_t NofPMTPoints, vector< Double_t > &normalPMT, Double_t &normalCste)
Definition: CbmRichProjectionProducer2.cxx:259
CbmRichGeoManager.h
CbmRichProjectionProducer2::fTile
TString fTile
Definition: CbmRichProjectionProducer2.h:96
CbmRichNavigationUtil2.h
CbmRichRecGeoParPmt::fPlaneX
Double_t fPlaneX
Definition: CbmRichRecGeoPar.h:66
z2
Double_t z2[nSects2]
Definition: pipe_v16a_mvdsts100.h:11
CbmRichRecGeoPar
PMT parameters for the RICH geometry.
Definition: CbmRichRecGeoPar.h:87
CbmRichProjectionProducer2::fTrackParams
TClonesArray * fTrackParams
Definition: CbmRichProjectionProducer2.h:91
CbmTrackMatchNew.h
CbmRichProjectionProducer2::Init
virtual void Init()
Initialization in case one needs to initialize some TCloneArrays.
Definition: CbmRichProjectionProducer2.cxx:40
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:407
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmRichNavigationUtil2::GetDirCos
static void GetDirCos(const FairTrackParam *par, Double_t &nx, Double_t &ny, Double_t &nz)
Definition: CbmRichNavigationUtil2.h:79
CbmUtils.h
CbmRichProjectionProducer2::CbmRichProjectionProducer2
CbmRichProjectionProducer2()
Definition: CbmRichProjectionProducer2.cxx:26
CbmRichNavigationUtil2::FindIntersection
static string FindIntersection(const FairTrackParam *par, TVector3 &crossPoint, const string &volumeName, TGeoNavigator *navi)
Definition: CbmRichNavigationUtil2.h:14
CbmRichProjectionProducer2::fRichPoints
TClonesArray * fRichPoints
Definition: CbmRichProjectionProducer2.h:93
CbmMCTrack.h
CbmMCTrack
Definition: CbmMCTrack.h:34
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichHitProducer.h
Class for producing RICH hits directly from MCPoints.
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichRecGeoParPmt::fWidth
Double_t fWidth
Definition: CbmRichRecGeoPar.h:70
CbmRichHit.h
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmRichGeoManager::fGP
CbmRichRecGeoPar * fGP
Definition: CbmRichGeoManager.h:23
CbmRichProjectionProducer2::fNumbAxis
TString fNumbAxis
Definition: CbmRichProjectionProducer2.h:95
CbmRichProjectionProducer2::fEventNum
Int_t fEventNum
Definition: CbmRichProjectionProducer2.h:98
CbmRichProjectionProducer2::ProjectionProducer
Double_t * ProjectionProducer(FairTrackParam *point)
Definition: CbmRichProjectionProducer2.cxx:135