11 #include "FairGeoNode.h"
12 #include "FairGeoTransform.h"
13 #include "FairGeoVector.h"
14 #include "FairRootManager.h"
16 #include "FairRunAna.h"
17 #include "FairRuntimeDb.h"
18 #include "FairTrackParam.h"
19 #include "TGeoManager.h"
21 #include "TClonesArray.h"
22 #include "TMatrixFSym.h"
28 #include "FairLogger.h"
38 : fTrackParams(NULL), fNHits(0), fEventNum(0) {}
41 FairRootManager* fManager = FairRootManager::Instance();
47 LOG(info) <<
"CbmRichProjectionProducerAnalytical::Init()";
48 FairRootManager* manager = FairRootManager::Instance();
50 fTrackParams = (TClonesArray*) manager->GetObject(
"RichTrackParamZ");
51 if (NULL ==
fTrackParams) { LOG(fatal) <<
"No RichTrackParamZ array!"; }
56 cout <<
"CbmRichProjectionProducerAnalytical:: event " <<
fEventNum << endl;
65 TMatrixFSym covMat(5);
66 for (Int_t
i = 0;
i < 5;
i++) {
67 for (Int_t j = 0; j <=
i; j++) {
71 covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
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);
79 if (point->GetX() == 0 && point->GetY() == 0 && point->GetZ() == 0
80 && point->GetTx() == 0 && point->GetTy() == 0)
82 if (point->GetQp() == 0)
continue;
85 TVector3 startP, momP, crossP, centerP;
88 Double_t p = 1. / TMath::Abs(point->GetQp());
90 if ((1 + point->GetTx() * point->GetTx() + point->GetTy() * point->GetTy())
93 / TMath::Sqrt(1 + point->GetTx() * point->GetTx()
94 + point->GetTy() * point->GetTy());
96 cout <<
" -E- RichProjectionProducer: strange value for calculating pz: "
97 << (1 + point->GetTx() * point->GetTx()
98 + point->GetTy() * point->GetTy())
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)
119 (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
120 + momP.z() * (startP.z() - mirrorZ));
122 (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
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);
129 if ((RxP * RxP - normP2 * dist) > 0.) {
131 rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
133 cout <<
" Error in track extrapolation: momentum = 0 " << endl;
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);
145 if ((mirrorY * crossP.y()) < 0) {
148 (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
149 + momP.z() * (startP.z() - mirrorZ));
151 (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
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);
158 if ((RxP * RxP - normP2 * dist) > 0.) {
160 rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
162 cout <<
" Error in track extrapolation: momentum = 0 " << endl;
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);
173 centerP.SetXYZ(mirrorX, mirrorY, mirrorZ);
177 TVector3 normP(crossP.x() - centerP.x(),
178 crossP.y() - centerP.y(),
179 crossP.z() - centerP.z());
180 normP = normP.Unit();
182 if ((normP.z() * momP.z()) < 0.)
183 normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
187 normP.x() * momP.x() + normP.y() * momP.y() + normP.z() * momP.z();
190 ref.SetXYZ(2 * np * normP.x() - momP.x(),
191 2 * np * normP.y() - momP.y(),
192 2 * np * normP.z() - momP.z());
196 TVector3 inPos(0., 0., 0.), outPos(0., 0., 0.);
213 Fatal(
"CbmRichProjectionProducerAnalytical ",
"unnown geometry type");
220 FairTrackParam richtrack(
221 outPos.x(), outPos.y(), outPos.z(), 0., 0., 0., covMat);
222 *(FairTrackParam*) (
richProj->At(j)) = richtrack;
229 const TVector3* centerP,
230 const TVector3* crossP,
232 TVector3* outPoint) {
249 "CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings ",
250 "Wrong geometry, this method could be used only for "
251 "CbmRichGeometryTypeTwoWings");
261 if (centerP->y() > 0) {
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());
270 if (centerP->y() < 0) {
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());
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;
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());
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());
307 xX = crossP->x() + ref->x() * rho2;
308 yY = crossP->y() + ref->y() * rho2;
309 zZ = crossP->z() + ref->z() * rho2;
312 outPoint->SetXYZ(xX, yY, zZ);
317 const TVector3* centerP,
318 const TVector3* crossP,
320 TVector3* outPoint) {
323 Fatal(
"CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl ",
324 "Wrong geometry, this method could be used only for "
325 "CbmRichGeometryTypeCylindrical");
330 Double_t maxDist = 0.;
332 for (map<string, CbmRichRecGeoParPmt>::iterator it = gp->
fPmtMap.begin();
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;
341 if (!(pmtPlaneX >= 0 && pmtPlaneY >= 0))
continue;
345 if (centerP->y() > 0) {
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());
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());
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;
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());
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());
392 cxX = crossP->x() + ref->x() * rho2;
393 cyY = crossP->y() + ref->y() * rho2;
394 czZ = crossP->z() + ref->z() * rho2;
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));
409 outPoint->SetXYZ(xX, yY, zZ);