CbmRoot
tracks/CbmRichProjectionProducerAnalytical.cxx
Go to the documentation of this file.
1 
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 
27 
28 #include "FairLogger.h"
29 
30 #include <cmath>
31 #include <iostream>
32 
33 using std::cout;
34 using std::endl;
35 
36 
38  : fTrackParams(NULL), fNHits(0), fEventNum(0) {}
39 
41  FairRootManager* fManager = FairRootManager::Instance();
42  fManager->Write();
43 }
44 
45 
47  LOG(info) << "CbmRichProjectionProducerAnalytical::Init()";
48  FairRootManager* manager = FairRootManager::Instance();
49 
50  fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
51  if (NULL == fTrackParams) { LOG(fatal) << "No RichTrackParamZ array!"; }
52 }
53 
55  fEventNum++;
56  cout << "CbmRichProjectionProducerAnalytical:: event " << fEventNum << endl;
57 
59  Double_t mirrorX = gp->fMirrorX;
60  Double_t mirrorY = gp->fMirrorY;
61  Double_t mirrorZ = gp->fMirrorZ;
62  Double_t mirrorR = gp->fMirrorR;
63 
64  richProj->Delete();
65  TMatrixFSym covMat(5);
66  for (Int_t i = 0; i < 5; i++) {
67  for (Int_t j = 0; j <= i; j++) {
68  covMat(i, j) = 0;
69  }
70  }
71  covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
72  1.e-4;
73 
74  for (Int_t j = 0; j < fTrackParams->GetEntriesFast(); j++) {
75  FairTrackParam* point = (FairTrackParam*) fTrackParams->At(j);
76  new ((*richProj)[j]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
77 
78  // check if Array was filled
79  if (point->GetX() == 0 && point->GetY() == 0 && point->GetZ() == 0
80  && point->GetTx() == 0 && point->GetTy() == 0)
81  continue;
82  if (point->GetQp() == 0) continue;
83 
84  Double_t rho1 = 0.;
85  TVector3 startP, momP, crossP, centerP;
86 
87 
88  Double_t p = 1. / TMath::Abs(point->GetQp());
89  Double_t pz;
90  if ((1 + point->GetTx() * point->GetTx() + point->GetTy() * point->GetTy())
91  > 0.)
92  pz = p
93  / TMath::Sqrt(1 + point->GetTx() * point->GetTx()
94  + point->GetTy() * point->GetTy());
95  else {
96  cout << " -E- RichProjectionProducer: strange value for calculating pz: "
97  << (1 + point->GetTx() * point->GetTx()
98  + point->GetTy() * point->GetTy())
99  << endl;
100  pz = 0.;
101  }
102  Double_t px = pz * point->GetTx();
103  Double_t py = pz * point->GetTy();
104  momP.SetXYZ(px, py, pz);
105  point->Position(startP);
106  if ((mirrorY * startP.y()) < 0)
107  mirrorY =
108  -mirrorY; // check that mirror center and startP are in same hemisphere
109 
110  // calculation of intersection of track with selected mirror
111  // corresponds to calculation of intersection between a straight line and a sphere:
112  // vector: r = startP - mirrorCenter
113  // RxP = r*momP
114  // normP2 = momP^2
115  // dist = r^2 - fR^2
116  // -> rho1 = (-RxP+sqrt(RxP^2-normP2*dist))/normP2 extrapolation factor for:
117  // intersection point crossP = startP + rho1 * momP
118  Double_t RxP =
119  (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
120  + momP.z() * (startP.z() - mirrorZ));
121  Double_t normP2 =
122  (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
123  Double_t dist =
124  (startP.x() * startP.x() + mirrorX * mirrorX + startP.y() * startP.y()
125  + mirrorY * mirrorY + startP.z() * startP.z() + mirrorZ * mirrorZ
126  - 2 * startP.x() * mirrorX - 2 * startP.y() * mirrorY
127  - 2 * startP.z() * mirrorZ - mirrorR * mirrorR);
128 
129  if ((RxP * RxP - normP2 * dist) > 0.) {
130  if (normP2 != 0.)
131  rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
132  if (normP2 == 0)
133  cout << " Error in track extrapolation: momentum = 0 " << endl;
134  } else {
135  //cout << " -E- RichProjectionProducer: RxP*RxP-normP2*dist = " << RxP*RxP-normP2*dist << endl;
136  }
137 
138  Double_t crossPx = startP.x() + rho1 * momP.x();
139  Double_t crossPy = startP.y() + rho1 * momP.y();
140  Double_t crossPz = startP.z() + rho1 * momP.z();
141  crossP.SetXYZ(crossPx, crossPy, crossPz);
142 
143  // check if crosspoint with mirror and chosen mirrorcenter (y) are in same hemisphere
144  // if not recalculate crossing point
145  if ((mirrorY * crossP.y()) < 0) {
146  mirrorY = -mirrorY;
147  RxP =
148  (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
149  + momP.z() * (startP.z() - mirrorZ));
150  normP2 =
151  (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
152  dist =
153  (startP.x() * startP.x() + mirrorX * mirrorX + startP.y() * startP.y()
154  + mirrorY * mirrorY + startP.z() * startP.z() + mirrorZ * mirrorZ
155  - 2 * startP.x() * mirrorX - 2 * startP.y() * mirrorY
156  - 2 * startP.z() * mirrorZ - mirrorR * mirrorR);
157 
158  if ((RxP * RxP - normP2 * dist) > 0.) {
159  if (normP2 != 0.)
160  rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
161  if (normP2 == 0)
162  cout << " Error in track extrapolation: momentum = 0 " << endl;
163  } else {
164  //cout << " -E- RichProjectionProducer: RxP*RxP-normP2*dist = " << RxP*RxP-normP2*dist << endl;
165  }
166 
167  crossPx = startP.x() + rho1 * momP.x();
168  crossPy = startP.y() + rho1 * momP.y();
169  crossPz = startP.z() + rho1 * momP.z();
170  crossP.SetXYZ(crossPx, crossPy, crossPz);
171  }
172 
173  centerP.SetXYZ(mirrorX, mirrorY, mirrorZ); // mirror center
174 
175 
176  // calculate normal on crosspoint with mirror
177  TVector3 normP(crossP.x() - centerP.x(),
178  crossP.y() - centerP.y(),
179  crossP.z() - centerP.z());
180  normP = normP.Unit();
181  // check that normal has same z-direction as momentum
182  if ((normP.z() * momP.z()) < 0.)
183  normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
184 
185  // reflect track
186  Double_t np =
187  normP.x() * momP.x() + normP.y() * momP.y() + normP.z() * momP.z();
188 
189  TVector3 ref;
190  ref.SetXYZ(2 * np * normP.x() - momP.x(),
191  2 * np * normP.y() - momP.y(),
192  2 * np * normP.z() - momP.z());
193 
194  if (ref.z() != 0.) {
195 
196  TVector3 inPos(0., 0., 0.), outPos(0., 0., 0.);
197 
198 
200  GetPmtIntersectionPointCyl(&centerP, &crossP, &ref, &inPos);
201  // Transform intersection point in same way as MCPoints were transformed in HitProducer before stored as Hit
202  // For the cylindrical geometry we also pass the path to the strip block
203  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
204 
205  //cout << "inPoint:" << inPos.X() << " " << inPos.Y() << " " << inPos.Z() << " outPoint:" << outPos.X() << " " << outPos.Y() << " " << outPos.Z() << endl;
206 
207 
208  } else if (gp->fGeometryType == CbmRichGeometryTypeTwoWings) {
209  GetPmtIntersectionPointTwoWings(&centerP, &crossP, &ref, &inPos);
210  // Transform intersection point in same way as MCPoints were transformed in HitProducer before stored as Hit
211  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
212  } else {
213  Fatal("CbmRichProjectionProducerAnalytical ", "unnown geometry type");
214  }
215 
216  Bool_t isInsidePmt =
218 
219  if (isInsidePmt) {
220  FairTrackParam richtrack(
221  outPos.x(), outPos.y(), outPos.z(), 0., 0., 0., covMat);
222  *(FairTrackParam*) (richProj->At(j)) = richtrack;
223  }
224  } // if (refZ!=0.)
225  } // j
226 }
227 
229  const TVector3* centerP,
230  const TVector3* crossP,
231  const TVector3* ref,
232  TVector3* outPoint) {
233  // crosspoint whith photodetector plane:
234  // calculate intersection between straight line and (tilted) plane:
235  // normal on plane tilted by theta around x-axis: (0,-sin(theta),cos(theta)) = n
236  // normal on plane tilted by phi around y-axis: (-sin(phi),0,cos(phi)) = n
237  // 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
238  // point on plane is (fDetX,fDetY,fDetZ) = p as photodetector is tiled around its center
239  // equation of plane for r being point in plane: n(r-p) = 0
240  // calculate intersection point of reflected track with plane: r=intersection point
241  // intersection point = crossP + rho2 * refl_track
242  // take care for all 4 cases:
243  // -> first calculate for case x>0, then check
244 
246 
248  Fatal(
249  "CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings ",
250  "Wrong geometry, this method could be used only for "
251  "CbmRichGeometryTypeTwoWings");
252  }
253 
254  Double_t pmtPhi = gp->fPmt.fPhi;
255  Double_t pmtTheta = gp->fPmt.fTheta;
256  Double_t pmtPlaneX = gp->fPmt.fPlaneX;
257  Double_t pmtPlaneY = gp->fPmt.fPlaneY;
258  Double_t pmtPlaneZ = gp->fPmt.fPlaneZ;
259  Double_t rho2 = 0.;
260 
261  if (centerP->y() > 0) {
262  rho2 =
263  (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
264  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
265  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
266  / (-TMath::Sin(pmtPhi) * ref->x()
267  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
268  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
269  }
270  if (centerP->y() < 0) {
271  rho2 =
272  (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
273  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * (-pmtPlaneY - crossP->y())
274  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
275  / (-TMath::Sin(pmtPhi) * ref->x()
276  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
277  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
278  }
279 
280  //rho2 = -1*(crossP.z() - fDetZ)/refZ; // only for theta = 0, phi=0
281  Double_t xX = crossP->x() + ref->x() * rho2;
282  Double_t yY = crossP->y() + ref->y() * rho2;
283  Double_t zZ = crossP->z() + ref->z() * rho2;
284 
285  if (xX < 0) {
286  if (centerP->y() > 0) {
287  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
288  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi)
289  * (pmtPlaneY - crossP->y())
290  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi)
291  * (pmtPlaneZ - crossP->z()))
292  / (-TMath::Sin(-pmtPhi) * ref->x()
293  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
294  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
295  }
296  if (centerP->y() < 0) {
297  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
298  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi)
299  * (-pmtPlaneY - crossP->y())
300  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi)
301  * (pmtPlaneZ - crossP->z()))
302  / (-TMath::Sin(-pmtPhi) * ref->x()
303  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
304  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
305  }
306 
307  xX = crossP->x() + ref->x() * rho2;
308  yY = crossP->y() + ref->y() * rho2;
309  zZ = crossP->z() + ref->z() * rho2;
310  }
311 
312  outPoint->SetXYZ(xX, yY, zZ);
313 }
314 
315 
317  const TVector3* centerP,
318  const TVector3* crossP,
319  const TVector3* ref,
320  TVector3* outPoint) {
323  Fatal("CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl ",
324  "Wrong geometry, this method could be used only for "
325  "CbmRichGeometryTypeCylindrical");
326  }
327  Double_t xX = 0.;
328  Double_t yY = 0.;
329  Double_t zZ = 0.;
330  Double_t maxDist = 0.;
331 
332  for (map<string, CbmRichRecGeoParPmt>::iterator it = gp->fPmtMap.begin();
333  it != gp->fPmtMap.end();
334  it++) {
335  Double_t pmtPlaneX = it->second.fPlaneX;
336  Double_t pmtPlaneY = it->second.fPlaneY;
337  Double_t pmtPlaneZ = it->second.fPlaneZ;
338  Double_t pmtPhi = it->second.fPhi;
339  Double_t pmtTheta = it->second.fTheta;
340 
341  if (!(pmtPlaneX >= 0 && pmtPlaneY >= 0)) continue;
342 
343  double rho2 = 0.;
344 
345  if (centerP->y() > 0) {
346  rho2 =
347  (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
348  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
349  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi)
350  * (pmtPlaneZ - crossP->z()))
351  / (-TMath::Sin(pmtPhi) * ref->x()
352  - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
353  + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
354  }
355  if (centerP->y() < 0) {
356  rho2 = (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
357  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi)
358  * (-pmtPlaneY - crossP->y())
359  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi)
360  * (pmtPlaneZ - crossP->z()))
361  / (-TMath::Sin(pmtPhi) * ref->x()
362  - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
363  + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
364  }
365 
366  Double_t cxX = crossP->x() + ref->x() * rho2;
367  Double_t cyY = crossP->y() + ref->y() * rho2;
368  Double_t czZ = crossP->z() + ref->z() * rho2;
369 
370  if (cxX < 0) {
371  if (centerP->y() > 0) {
372  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
373  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi)
374  * (pmtPlaneY - crossP->y())
375  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi)
376  * (pmtPlaneZ - crossP->z()))
377  / (-TMath::Sin(-pmtPhi) * ref->x()
378  - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
379  + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
380  }
381  if (centerP->y() < 0) {
382  rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
383  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi)
384  * (-pmtPlaneY - crossP->y())
385  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi)
386  * (pmtPlaneZ - crossP->z()))
387  / (-TMath::Sin(-pmtPhi) * ref->x()
388  - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
389  + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
390  }
391 
392  cxX = crossP->x() + ref->x() * rho2;
393  cyY = crossP->y() + ref->y() * rho2;
394  czZ = crossP->z() + ref->z() * rho2;
395  }
396 
397  Double_t dP = TMath::Sqrt((crossP->x() - cxX) * (crossP->X() - cxX)
398  + (crossP->y() - cyY) * (crossP->y() - cyY)
399  + (crossP->z() - czZ) * (crossP->Z() - czZ));
400 
401  if (dP >= maxDist) {
402  maxDist = dP;
403  xX = cxX;
404  yY = cyY;
405  zZ = czZ;
406  }
407  }
408 
409  outPoint->SetXYZ(xX, yY, zZ);
410 }
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
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
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
CbmRichRecGeoPar::fGeometryType
CbmRichGeometryType fGeometryType
Definition: CbmRichRecGeoPar.h:220
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
CbmRichRecGeoParPmt::fPlaneX
Double_t fPlaneX
Definition: CbmRichRecGeoPar.h:66
CbmRichRecGeoPar
PMT parameters for the RICH geometry.
Definition: CbmRichRecGeoPar.h:87
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
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::DoProjection
virtual void DoProjection(TClonesArray *richProj)
Execute task.
Definition: alignment/CbmRichProjectionProducerAnalytical.cxx:64
CbmRichGeoManager::fGP
CbmRichRecGeoPar * fGP
Definition: CbmRichGeoManager.h:23
CbmRichRecGeoParPmt::fPlaneZ
Double_t fPlaneZ
Definition: CbmRichRecGeoPar.h:68