CbmRoot
CbmRichGeoManager.cxx
Go to the documentation of this file.
1 /*
2  * CbmRichGeoManager.cxx
3  *
4  * Created on: Dec 16, 2015
5  * Author: slebedev
6  */
7 #include "CbmRichGeoManager.h"
8 
9 #include "CbmRichRecGeoPar.h" // for CbmRichRecGeoParPmt, CbmRichRecGeoPar
10 
11 #include <FairLogger.h> // for LOG
12 
13 #include <TGeoBBox.h> // for TGeoBBox
14 #include <TGeoManager.h> // for TGeoManager, gGeoManager
15 #include <TGeoMatrix.h> // for TGeoMatrix
16 #include <TGeoNode.h> // for TGeoIterator, TGeoNode
17 #include <TGeoShape.h> // for TGeoShape
18 #include <TGeoSphere.h> // for TGeoSphere
19 #include <TGeoVolume.h> // for TGeoVolume
20 #include <TMath.h> // for ASin, Cos, ACos, Sin, DegToRad, Sqrt
21 #include <TMathBase.h> // for Abs, Max, Min
22 #include <TString.h> // for TString, operator==
23 #include <TVector3.h> // for TVector3
24 
25 
26 #include <algorithm> // for sort
27 #include <map> // for map, __map_iterator, map<>::iterator
28 #include <stddef.h> // for size_t
29 #include <string> // for operator<, operator+
30 #include <utility> // for pair
31 #include <vector> // for vector
32 
33 using namespace std;
34 
35 CbmRichGeoManager::CbmRichGeoManager() : fGP(nullptr) { InitGeometry(); }
36 
37 void CbmRichGeoManager::InitGeometry() {
38 
39  LOG(info) << "CbmRichGeoManager::InitGeometry";
40 
41  fGP = new CbmRichRecGeoPar();
42  //TODO: get refractive index from material
43  fGP->fNRefrac = 1.000446242;
44 
45  DetectGeometryType();
46 
47  if (fGP->fGeometryType == CbmRichGeometryTypeTwoWings) {
48  LOG(info) << "CbmRichGeoManager: Init CbmRichGeometryTypeTwoWings";
49  InitPmt();
50  } else if (fGP->fGeometryType == CbmRichGeometryTypeCylindrical) {
51  LOG(info) << "CbmRichGeoManager: Init CbmRichGeometryTypeCylindrical";
52  InitPmtCyl();
53  } else if (fGP->fGeometryType == CbmRichGeometryTypeNotDefined) {
54  // Fatal("CbmRichGeoManager::InitGeometry()", " Geometry type is CbmRichGeometryTypeNotDefined. Geometry could not be defined automatically.");
55  LOG(info);
56  LOG(error) << "CbmRichGeoManager::InitGeometry(): Geometry type is "
57  "CbmRichGeometryTypeNotDefined. Geometry could not be "
58  "defined automatically.";
59  }
60 
61  InitMirror();
62 
63  //fGP->Print();
64 }
65 
66 void CbmRichGeoManager::DetectGeometryType() {
67  TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
68  geoIterator.SetTopName("/cave_1");
69  TGeoNode* curNode;
70 
71  geoIterator.Reset();
72  while ((curNode = geoIterator())) {
73  TString nodeName(curNode->GetName());
74  TString nodePath;
75  // if (curNode->GetVolume()->GetName() == TString("camera_two_quarters")) {
76  if (curNode->GetVolume()->GetName() == TString("camera_strip")) {
77  fGP->fGeometryType = CbmRichGeometryTypeCylindrical;
78  return;
79  }
80  }
81 
82  geoIterator.Reset();
83  while ((curNode = geoIterator())) {
84  TString nodeName(curNode->GetName());
85  TString nodePath;
86  if (curNode->GetVolume()->GetName() == TString("camera_quarter")) {
87  fGP->fGeometryType = CbmRichGeometryTypeTwoWings;
88  return;
89  }
90  }
91 
92  fGP->fGeometryType = CbmRichGeometryTypeNotDefined;
93 }
94 
95 void CbmRichGeoManager::InitPmtCyl() {
96  TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
97  geoIterator.SetTopName("/cave_1");
98  TGeoNode* curNode;
99 
100  vector<Double_t> xCoord;
101  geoIterator.Reset();
102  while ((curNode = geoIterator())) {
103  TString nodeName(curNode->GetName());
104  TString nodePath;
105  // if (curNode->GetVolume()->GetName() == TString("pmt_block_strip")) {
106  if (curNode->GetVolume()->GetName() == TString("camera_strip")) {
107  geoIterator.GetPath(nodePath);
108  //cout << "nodePath:" << nodePath << endl;
109  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
110  const Double_t* curNodeTr = curMatrix->GetTranslation();
111  const Double_t* curNodeRot = curMatrix->GetRotationMatrix();
112 
113  //curMatrix->Print();
114 
115  double pmtX = curNodeTr[0];
116  double pmtY = curNodeTr[1];
117  double pmtZ = curNodeTr[2];
118 
119  double rotY = TMath::ASin(-curNodeRot[2]); // around Y
120  // double rotZ = TMath::ASin(curNodeRot[1]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around Z
121  //double rotX = TMath::ASin(curNodeRot[5]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X
122  double rotX = TMath::ACos(
123  curNodeRot[8] / TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X
124 
125 
126  //double rotX = TMath::ATan2(curNodeRot[5], curNodeRot[8]); // around Y
127  //double c2 = TMath::Sqrt(curNodeRot[0]*curNodeRot[0] + curNodeRot[1]*curNodeRot[1]);
128  //double rotY = TMath::ATan2(-curNodeRot[2], c2); // around X
129 
130  string nodePathStr = string(nodePath.Data()) + "/";
131 
132  fGP->fPmtMap[nodePathStr].fTheta = rotX;
133  fGP->fPmtMap[nodePathStr].fPhi = rotY;
134  const TGeoBBox* shape =
135  (const TGeoBBox*) (curNode->GetVolume()->GetShape());
136 
137  fGP->fPmtMap[nodePathStr].fWidth = 2. * shape->GetDX();
138  fGP->fPmtMap[nodePathStr].fHeight = 2. * shape->GetDY();
139  fGP->fPmtMap[nodePathStr].fZ = pmtZ;
140  fGP->fPmtMap[nodePathStr].fX = pmtX;
141  fGP->fPmtMap[nodePathStr].fY = pmtY;
142 
143  if (pmtX >= 0 && pmtY >= 0) { xCoord.push_back(pmtX); }
144  }
145  }
146  std::sort(xCoord.begin(), xCoord.end());
147 
148  for (map<string, CbmRichRecGeoParPmt>::iterator it = fGP->fPmtMap.begin();
149  it != fGP->fPmtMap.end();
150  it++) {
151  Double_t curX = TMath::Abs(it->second.fX);
152  int pos = -1;
153  //int pos = find(xCoord.begin(), xCoord.end(), curX) - xCoord.begin();
154  for (unsigned int i = 0; i < xCoord.size(); i++) {
155  if (TMath::Abs(curX - xCoord[i]) < 0.1) {
156  pos = i;
157  break;
158  }
159  }
160  it->second.fPmtPositionIndexX = pos;
161  }
162 
163  // We also need to find pmt plane center
164  map<string, CbmRichPmtPlaneMinMax> mapPmtPlaneMinMax;
165  TString filterNamePixel("pmt_pixel");
166  geoIterator.Reset();
167  CbmRichPmtPlaneMinMax pmtPlaneMinMax;
168 
169  while ((curNode = geoIterator())) {
170  TString nodeName(curNode->GetName());
171  TString nodePath;
172  if (curNode->GetVolume()->GetName() == filterNamePixel) {
173 
174  geoIterator.GetPath(nodePath);
175 
176 
177  string nodePathStr = string(nodePath.Data());
178  //size_t foundIndex1 = nodePathStr.find("pmt_block_strip");
179  size_t foundIndex1 = nodePathStr.find("camera_strip");
180 
181  if (foundIndex1 == string::npos) continue;
182  size_t foundIndex2 = nodePathStr.find("/", foundIndex1 + 1);
183  if (foundIndex2 == string::npos) continue;
184 
185  string mapKey = nodePathStr.substr(0, foundIndex2) + "/";
186 
187  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
188  const Double_t* curNodeTr = curMatrix->GetTranslation();
189 
190  double pmtX = curNodeTr[0];
191  double pmtY = curNodeTr[1];
192  double pmtZ = curNodeTr[2];
193 
194  mapPmtPlaneMinMax[mapKey].AddPoint(pmtX, pmtY, pmtZ);
195  }
196  }
197 
198  for (map<string, CbmRichRecGeoParPmt>::iterator it = fGP->fPmtMap.begin();
199  it != fGP->fPmtMap.end();
200  it++) {
201  it->second.fPlaneX = mapPmtPlaneMinMax[it->first].GetMeanX();
202  it->second.fPlaneY = mapPmtPlaneMinMax[it->first].GetMeanY();
203  it->second.fPlaneZ = mapPmtPlaneMinMax[it->first].GetMeanZ();
204 
205  // cout << "name:" << it->first << " strip(x,y,z):" <<it->second.fX << "," << it->second.fY << "," << it->second.fZ <<
206  // " pmtPlane(x,y,z):" <<it->second.fPlaneX << "," << it->second.fPlaneY << "," << it->second.fPlaneZ << ", " <<
207  // "theta:" << it->second.fTheta << ", phi:" << it->second.fPhi << endl;
208  }
209 
210  // Calculate gap between camera_strip
211  geoIterator.Reset();
212  double master1[3], master2[3];
213  while ((curNode = geoIterator())) {
214  TString nodeName(curNode->GetName());
215  TString nodePath;
216  if (curNode->GetVolume()->GetName() == TString("camera_strip")) {
217 
218  geoIterator.GetPath(nodePath);
219  string nodePathStr = string(nodePath.Data()) + "/";
220  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
221  const Double_t* curNodeTr = curMatrix->GetTranslation();
222  // const Double_t* curNodeRot = curMatrix->GetRotationMatrix();
223 
224  double pmtX = curNodeTr[0];
225  double pmtY = curNodeTr[1];
226  // double pmtZ = curNodeTr[2];
227 
228  if (pmtX < 0 || pmtY < 0) continue;
229  const TGeoBBox* shape =
230  (const TGeoBBox*) (curNode->GetVolume()->GetShape());
231 
232 
233  double loc[3];
234  if (fGP->fPmtMap[nodePathStr].fPmtPositionIndexX == 1) {
235  loc[0] = shape->GetDX() + shape->GetOrigin()[0];
236  loc[1] = shape->GetDY() + shape->GetOrigin()[1];
237  loc[2] = shape->GetDZ() + shape->GetOrigin()[2];
238  curMatrix->LocalToMaster(loc, master1);
239  } else if (fGP->fPmtMap[nodePathStr].fPmtPositionIndexX == 2) {
240  loc[0] = -shape->GetDX() + shape->GetOrigin()[0];
241  loc[1] = shape->GetDY() + shape->GetOrigin()[1];
242  loc[2] = shape->GetDZ() + shape->GetOrigin()[2];
243  curMatrix->LocalToMaster(loc, master2);
244  }
245  }
246  }
247  //cout << master1[0] << " "<< master1[1] << " "<< master1[2]<< endl;
248  //cout << master2[0] << " "<< master2[1] << " "<< master2[2]<< endl;
249  double dist =
250  TMath::Sqrt((master1[0] - master2[0]) * (master1[0] - master2[0])
251  + (master1[1] - master2[1]) * (master1[1] - master2[1])
252  + (master1[2] - master2[2]) * (master1[2] - master2[2]));
253  //cout << "Gap:" << dist << endl;
254  fGP->fPmtStripGap = dist; // v17a = 0.3083394
255 }
256 
257 void CbmRichGeoManager::InitPmt() {
258  TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
259  TGeoNode* curNode;
260 
261  // PMT plane position\rotation
262  TString filterName_pixel("pmt_pixel");
263  geoIterator.Reset();
264  CbmRichPmtPlaneMinMax pmtPlaneMinMax;
265  while ((curNode = geoIterator())) {
266  TString nodeName(curNode->GetName());
267  TString nodePath;
268  if (curNode->GetVolume()->GetName() == filterName_pixel) {
269 
270  geoIterator.GetPath(nodePath);
271  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
272  const Double_t* curNodeTr = curMatrix->GetTranslation();
273  const Double_t* curNodeRot = curMatrix->GetRotationMatrix();
274 
275  double pmtX = curNodeTr[0];
276  double pmtY = curNodeTr[1];
277  double pmtZ = curNodeTr[2];
278 
279  if (pmtX > 0. && pmtY > 0) {
280  //printf ("%08f\t%08f\t%08f\t\n", curNodeTranslation[0], curNodeTranslation[1], curNodeTranslation[2]);
281  double rotY = TMath::ASin(-curNodeRot[2]); // around Y
282  // double rotZ = TMath::ASin(curNodeRot[1]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around Z
283  //double rotX = TMath::ASin(curNodeRot[5]/TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X
284  double rotX = TMath::ACos(
285  curNodeRot[8] / TMath::Cos(TMath::ASin(-curNodeRot[2]))); // around X
286 
287  fGP->fPmt.fTheta = rotX;
288  fGP->fPmt.fPhi = rotY;
289 
290  pmtPlaneMinMax.AddPoint(pmtX, pmtY, pmtZ);
291  }
292  }
293  }
294 
295  fGP->fPmt.fPlaneX = pmtPlaneMinMax.GetMeanX();
296  fGP->fPmt.fPlaneY = pmtPlaneMinMax.GetMeanY();
297  fGP->fPmt.fPlaneZ = pmtPlaneMinMax.GetMeanZ();
298 
299  geoIterator.Reset();
300  while ((curNode = geoIterator())) {
301  TString nodeName(curNode->GetName());
302  TString nodePath;
303  if (curNode->GetVolume()->GetName() == TString("camera_quarter")) {
304 
305  geoIterator.GetPath(nodePath);
306  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
307  const Double_t* curNodeTr = curMatrix->GetTranslation();
308  //const Double_t* curNodeRot = curMatrix->GetRotationMatrix();
309 
310 
311  double pmtX = curNodeTr[0];
312  double pmtY = curNodeTr[1];
313  double pmtZ = curNodeTr[2];
314 
315  if (pmtX > 0. && pmtY > 0) {
316  const TGeoBBox* shape =
317  (const TGeoBBox*) (curNode->GetVolume()->GetShape());
318  fGP->fPmt.fWidth = shape->GetDX();
319  fGP->fPmt.fHeight = shape->GetDY();
320  fGP->fPmt.fZ = pmtZ;
321  fGP->fPmt.fX = pmtX;
322  fGP->fPmt.fY = pmtY;
323  }
324  }
325  }
326 }
327 
328 void CbmRichGeoManager::InitMirror() {
329 
330  TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
331  TGeoNode* curNode;
332  geoIterator.Reset();
333 
334  //mirror position\rotation
335  TString mirrorName0("mirror_tile_type0");
336  TString mirrorName1("mirror_tile_type1");
337  TString mirrorName2("mirror_tile_type2");
338  TString mirrorName3("mirror_tile_type3");
339  //TString mirrorName4("mirror_tile_type4");
340  //TString mirrorName5("mirror_tile_type5");
341 
342  // these names are needed for misaligned geometry
343  TString mirrorMisAlignName0("mirror_tile_0");
344  TString mirrorMisAlignName1("mirror_tile_1");
345  TString mirrorMisAlignName2("mirror_tile_2");
346  TString mirrorMisAlignName3("mirror_tile_3");
347  TString mirrorMisAlignName4("mirror_tile_4");
348  TString mirrorMisAlignName5("mirror_tile_5");
349  TString mirrorMisAlignName6("mirror_tile_6");
350  TString mirrorMisAlignName7("mirror_tile_7");
351 
352  geoIterator.Reset();
353  double maxTheta = 0.;
354  double minTheta = 999999999.;
355  double mirrorX = 0.;
356  double mirrorY = 0.;
357  double mirrorZ = 0.;
358  double mirrorRadius = 0.;
359  while ((curNode = geoIterator())) {
360  TString nodeName(curNode->GetName());
361  TString nodePath;
362 
363  if (nodeName.Contains(mirrorName0) || nodeName.Contains(mirrorName1)
364  || nodeName.Contains(mirrorName2) || nodeName.Contains(mirrorName3)
365  || nodeName.Contains(mirrorMisAlignName0)
366  || nodeName.Contains(mirrorMisAlignName1)
367  || nodeName.Contains(mirrorMisAlignName2)
368  || nodeName.Contains(mirrorMisAlignName3)
369  || nodeName.Contains(mirrorMisAlignName4)
370  || nodeName.Contains(mirrorMisAlignName5)
371  || nodeName.Contains(mirrorMisAlignName6)
372  || nodeName.Contains(mirrorMisAlignName7)) {
373  geoIterator.GetPath(nodePath);
374  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
375  const Double_t* curNodeTr = curMatrix->GetTranslation();
376  mirrorX = TMath::Abs(curNodeTr[0]);
377  mirrorY = TMath::Abs(curNodeTr[1]);
378  mirrorZ = TMath::Abs(curNodeTr[2]);
379 
380  const TGeoSphere* shape =
381  dynamic_cast<const TGeoSphere*>(curNode->GetVolume()->GetShape());
382  //const TGeoSphere* shape = (const TGeoSphere*)(curNode->GetVolume()->GetShape());
383  if (shape != nullptr) {
384  mirrorRadius = shape->GetRmin();
385 
386  double theta1 = shape->GetTheta1();
387  double theta2 = shape->GetTheta2();
388  if (maxTheta < theta1 || maxTheta < theta2) {
389  maxTheta = TMath::Max(theta1, theta2);
390  }
391  if (minTheta > theta1 || minTheta > theta2) {
392  minTheta = TMath::Min(theta1, theta2);
393  }
394  }
395  }
396  }
397 
398  fGP->fMirrorTheta = -((maxTheta + minTheta) / 2. - 90.)
399  * TMath::DegToRad(); // rotation angle around x-axis
400  fGP->fMirrorX = mirrorX;
401  fGP->fMirrorY = mirrorY;
402  fGP->fMirrorZ = mirrorZ;
403  fGP->fMirrorR = mirrorRadius;
404 }
405 
406 
407 void CbmRichGeoManager::RotatePoint(TVector3* inPos,
408  TVector3* outPos,
409  Bool_t noTilting) {
410  if (fGP == nullptr) {
411  LOG(error) << "CbmRichGeoManager::RotatePoint RICH geometry is not "
412  "initialized. fGP == nullptr";
413  }
414 
415  if (fGP->fGeometryType == CbmRichGeometryTypeTwoWings) {
416  RotatePointTwoWings(inPos, outPos, noTilting);
417  } else if (fGP->fGeometryType == CbmRichGeometryTypeCylindrical) {
418  RotatePointCyl(inPos, outPos, noTilting);
419  } else {
420  outPos->SetXYZ(inPos->X(), inPos->Y(), inPos->Z());
421  }
422 }
423 
424 void CbmRichGeoManager::RotatePointTwoWings(TVector3* inPos,
425  TVector3* outPos,
426  Bool_t noTilting) {
427  if (noTilting == false) {
428  RotatePointImpl(inPos,
429  outPos,
430  fGP->fPmt.fPhi,
431  fGP->fPmt.fTheta,
432  fGP->fPmt.fX,
433  fGP->fPmt.fY,
434  fGP->fPmt.fZ);
435  } else {
436  outPos->SetXYZ(inPos->X(), inPos->Y(), inPos->Z());
437  }
438 }
439 
440 
441 void CbmRichGeoManager::RotatePointCyl(TVector3* inPos,
442  TVector3* outPos,
443  Bool_t noTilting,
444  Bool_t noShift) {
445  if (noTilting == false) {
446  // TGeoNode* node = gGeoManager->FindNode(inPos->X(), inPos->Y(), inPos->Z());
447  gGeoManager->FindNode(inPos->X(), inPos->Y(), inPos->Z());
448  string path(gGeoManager->GetPath());
449 
450 
451  CbmRichRecGeoParPmt pmtPar =
452  fGP->GetGeoRecPmtByBlockPathOrClosest(path, inPos);
453 
454  RotatePointImpl(inPos,
455  outPos,
456  -TMath::Abs(pmtPar.fPhi),
457  TMath::Abs(pmtPar.fTheta),
458  TMath::Abs(pmtPar.fX),
459  TMath::Abs(pmtPar.fY),
460  TMath::Abs(pmtPar.fZ));
461 
462  // After rotation wee need to shift point
463  if (!noShift) {
464  // All pmt blocks centers will be move to this Y position
465  // TODO: We can also take smallest Y from all pmt blocks
466  Double_t baseLineY = (outPos->Y() >= 0) ? 160. : -160.; //cm
467  Double_t dY = pmtPar.fY - baseLineY;
468  outPos->SetY(outPos->Y() - dY);
469 
470  // Calculate pmt block center after rotation
471  TVector3 inPosPmt(pmtPar.fX, pmtPar.fY, pmtPar.fZ);
472  TVector3 outPosPmt;
473  RotatePointImpl(&inPosPmt,
474  &outPosPmt,
475  -TMath::Abs(pmtPar.fPhi),
476  TMath::Abs(pmtPar.fTheta),
477  TMath::Abs(pmtPar.fX),
478  TMath::Abs(pmtPar.fY),
479  TMath::Abs(pmtPar.fZ));
480  //RotatePointCyl(&inPosPmt, &outPosPmt, false, true);
481 
482  // calculate ideal position assuming the same width for all pmt blocks
483  //TODO:Actually we need to implement general solution if pmt-block widths are not the same
484  Double_t gap = fGP->fPmtStripGap;
485  Double_t padding = gap / 2.;
486  Double_t idealX = padding + 0.5 * pmtPar.fWidth
487  + pmtPar.fPmtPositionIndexX * (pmtPar.fWidth + gap);
488  if (outPos->X() < 0) idealX = -idealX;
489  Double_t dX = idealX - outPosPmt.X();
490 
491  outPos->SetX(outPos->X() + dX);
492  }
493  } else {
494  outPos->SetXYZ(inPos->X(), inPos->Y(), inPos->Z());
495  }
496 }
497 
498 
499 void CbmRichGeoManager::RotatePointImpl(TVector3* inPos,
500  TVector3* outPos,
501  Double_t phi,
502  Double_t theta,
503  Double_t pmtX,
504  Double_t pmtY,
505  Double_t pmtZ) {
506  Double_t xDet = 0., yDet = 0., zDet = 0.;
507  Double_t x = inPos->X();
508  Double_t y = inPos->Y();
509  Double_t z = inPos->Z();
510 
511  Double_t sinTheta = TMath::Sin(theta);
512  Double_t cosTheta = TMath::Cos(theta);
513  Double_t sinPhi = TMath::Sin(phi);
514  Double_t cosPhi = TMath::Cos(phi);
515 
516  if (x > 0 && y > 0) {
517  y -= pmtY;
518  x -= pmtX;
519  z -= pmtZ;
520  //xDet = x*cosPhi + z*sinPhi;// - detZOrig*sinPhi;
521  //yDet = -x*sinTheta*sinPhi + y*cosTheta + z*sinTheta*cosPhi;
522  //zDet = -x*cosTheta*sinPhi - y*sinTheta + z*cosTheta*cosPhi;
523 
524  xDet = x * cosPhi - y * sinPhi * sinTheta + z * cosTheta * sinPhi;
525  yDet = y * cosTheta + z * sinTheta;
526  zDet = -x * sinPhi - y * sinTheta * cosPhi + z * cosTheta * cosPhi;
527 
528  yDet += pmtY;
529  xDet += pmtX;
530  zDet += pmtZ;
531 
532  } else if (x > 0 && y < 0) {
533  y += pmtY;
534  x -= pmtX;
535  z -= pmtZ;
536  // xDet = x*cosPhi + z*sinPhi;// - detZOrig*sinPhi;
537  //yDet = x*sinTheta*sinPhi + y*cosTheta - z*sinTheta*cosPhi;
538  //zDet = -x*cosTheta*sinPhi + y*sinTheta + z*cosTheta*cosPhi;
539 
540  xDet = x * cosPhi + y * sinPhi * sinTheta + z * cosTheta * sinPhi;
541  yDet = y * cosTheta - z * sinTheta;
542  zDet = -x * sinPhi + y * sinTheta * cosPhi + z * cosTheta * cosPhi;
543 
544  yDet -= pmtY;
545  xDet += pmtX;
546  zDet += pmtZ;
547  } else if (x < 0 && y < 0) {
548  y += pmtY;
549  x += pmtX;
550  z -= pmtZ;
551  // xDet = x*cosPhi - z*sinPhi;// + detZOrig*sinPhi;
552  //yDet = -x*sinTheta*sinPhi + y*cosTheta - z*sinTheta*cosPhi;
553  //zDet = x*cosTheta*sinPhi + y*sinTheta + z*cosTheta*cosPhi;
554 
555  xDet = x * cosPhi - y * sinPhi * sinTheta - z * cosTheta * sinPhi;
556  yDet = y * cosTheta - z * sinTheta;
557  zDet = x * sinPhi + y * sinTheta * cosPhi + z * cosTheta * cosPhi;
558 
559  yDet -= pmtY;
560  xDet -= pmtX;
561  zDet += pmtZ;
562  } else if (x < 0 && y > 0) {
563  y -= pmtY;
564  x += pmtX;
565  z -= pmtZ;
566  //xDet = x*cosPhi - z*sinPhi;// + detZOrig*sinPhi;
567  //yDet = x*sinTheta*sinPhi + y*cosTheta + z*sinTheta*cosPhi;
568  //zDet = x*cosTheta*sinPhi - y*sinTheta + z*cosTheta*cosPhi;
569 
570  xDet = x * cosPhi + y * sinPhi * sinTheta - z * cosTheta * sinPhi;
571  yDet = y * cosTheta + z * sinTheta;
572  zDet = x * sinPhi - y * sinTheta * cosPhi + z * cosTheta * cosPhi;
573 
574 
575  yDet += pmtY;
576  xDet -= pmtX;
577  zDet += pmtZ;
578  } else {
579  outPos->SetXYZ(x, y, z);
580  }
581  outPos->SetXYZ(xDet, yDet, zDet);
582 }
583 
584 Bool_t CbmRichGeoManager::IsPointInsidePmt(const TVector3* rotatedPoint) {
585  if (fGP->fGeometryType == CbmRichGeometryTypeTwoWings) {
586  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
587  Double_t pmtPlaneX = gp->fPmt.fPlaneX;
588  Double_t pmtPlaneY = gp->fPmt.fPlaneY;
589  Double_t pmtWidth = gp->fPmt.fWidth; // half width
590  Double_t pmtHeight = gp->fPmt.fHeight; // half height
591 
592  Double_t marginX = 2.; // [cm]
593  Double_t marginY = 2.; // [cm]
594  // upper pmt planes
595  Double_t pmtYTop = TMath::Abs(pmtPlaneY) + pmtHeight + marginY;
596  Double_t pmtYBottom = TMath::Abs(pmtPlaneY) - pmtHeight - marginY;
597  Double_t absYDet = TMath::Abs(rotatedPoint->y());
598  Bool_t isYOk = (absYDet <= pmtYTop && absYDet >= pmtYBottom);
599 
600  Double_t pmtXMin = -TMath::Abs(pmtPlaneX) - pmtWidth - marginX;
601  Double_t pmtXMax = TMath::Abs(pmtPlaneX) + pmtWidth + marginX;
602  Bool_t isXOk =
603  (rotatedPoint->x() >= pmtXMin && rotatedPoint->x() <= pmtXMax);
604 
605  return (isXOk && isYOk);
606  } else if (fGP->fGeometryType == CbmRichGeometryTypeCylindrical) {
607  CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
608  map<string, CbmRichRecGeoParPmt> pmtMap = gp->fPmtMap;
609 
610  int maxIndex = -1;
611  string maxKey = "";
612  for (auto const& entPmt : pmtMap) {
613  if (maxIndex < entPmt.second.fPmtPositionIndexX && entPmt.second.fX > 0
614  && entPmt.second.fY > 0) {
615  maxKey = entPmt.first;
616  maxIndex = entPmt.second.fPmtPositionIndexX;
617  }
618  }
619  CbmRichRecGeoParPmt pmtPar = pmtMap[maxKey];
620  Double_t pmtPlaneX = pmtPar.fPlaneX;
621  Double_t pmtPlaneY = pmtPar.fPlaneY;
622  Double_t pmtPlaneZ = pmtPar.fPlaneZ;
623  Double_t pmtWidth = pmtPar.fWidth; // full width of the strip
624  Double_t pmtHeight = pmtPar.fHeight; // full height
625 
626  //cout << "IsPointInsidePmt" << endl;
627  //cout << "maxIndex:" << maxIndex << " maxKey:" << maxKey << endl;
628  //cout << "pmtPlaneX:"<< pmtPlaneX << " pmtPlaneY:" << pmtPlaneY << " pmtPlaneZ:" << pmtPlaneZ<< " pmtWidth:" << pmtWidth << " pmtHeight:" <<pmtHeight << endl;
629  TVector3 inPmt(pmtPlaneX, pmtPlaneY, pmtPlaneZ), outPmt(0., 0., 0.);
630  CbmRichGeoManager::GetInstance().RotatePoint(&inPmt, &outPmt);
631  //cout << "Rotated pmtPlaneX:"<< outPmt.X() << " pmtPlaneY:" << outPmt.Y() << " pmtPlaneZ:" << outPmt.Z()<< endl;
632 
633  Double_t marginX = 2.; // [cm]
634  Double_t marginY = 2.; // [cm]
635  // upper pmt planes
636  Double_t pmtYTop = TMath::Abs(outPmt.Y()) + 0.5 * pmtHeight + marginY;
637  Double_t pmtYBottom = TMath::Abs(outPmt.Y()) - 0.5 * pmtHeight - marginY;
638  //cout << "pmtYTop:"<< pmtYTop << " pmtYBottom:" << pmtYBottom << endl;
639  Double_t absYDet = TMath::Abs(rotatedPoint->y());
640  Bool_t isYOk = (absYDet <= pmtYTop && absYDet >= pmtYBottom);
641 
642  Double_t pmtXMin = -TMath::Abs(outPmt.X()) - 0.5 * pmtWidth - marginX;
643  Double_t pmtXMax = TMath::Abs(outPmt.X()) + 0.5 * pmtWidth + marginX;
644  //cout << "pmtXMin:"<< pmtXMin << " pmtXMax:" << pmtXMax << endl;
645  Bool_t isXOk =
646  (rotatedPoint->x() >= pmtXMin && rotatedPoint->x() <= pmtXMax);
647 
648  return (isXOk && isYOk);
649  // return true;
650  } else {
651  return false;
652  }
653 }
CbmRichGeoManager::CbmRichGeoManager
CbmRichGeoManager()
Definition: CbmRichGeoManager.cxx:35
CbmRichRecGeoPar::fMirrorZ
Double_t fMirrorZ
Definition: CbmRichRecGeoPar.h:232
CbmRichRecGeoPar::fNRefrac
Double_t fNRefrac
Definition: CbmRichRecGeoPar.h:228
CbmRichRecGeoParPmt::fHeight
Double_t fHeight
Definition: CbmRichRecGeoPar.h:71
CbmRichGeometryTypeTwoWings
@ CbmRichGeometryTypeTwoWings
Definition: CbmRichRecGeoPar.h:24
CbmRichRecGeoParPmt::fPhi
Double_t fPhi
Definition: CbmRichRecGeoPar.h:58
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
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
CbmRichPmtPlaneMinMax
This class is used to store pmt_pixel min and max positions.
Definition: CbmRichRecGeoPar.h:246
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRecGeoPar::fMirrorX
Double_t fMirrorX
Definition: CbmRichRecGeoPar.h:230
CbmRichRecGeoPar.h
RICH geometry parameters for the reconstruction. This class is used for convinient storing of the bas...
CbmRichGeometryTypeNotDefined
@ CbmRichGeometryTypeNotDefined
Definition: CbmRichRecGeoPar.h:23
CbmRichGeoManager::RotatePointCyl
void RotatePointCyl(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false, Bool_t noShift=false)
Definition: CbmRichGeoManager.cxx:441
CbmRichGeoManager::InitGeometry
void InitGeometry()
Definition: CbmRichGeoManager.cxx:37
CbmRichGeoManager::DetectGeometryType
void DetectGeometryType()
Definition: CbmRichGeoManager.cxx:66
CbmRichRecGeoPar::fGeometryType
CbmRichGeometryType fGeometryType
Definition: CbmRichRecGeoPar.h:220
CbmRichRecGeoParPmt::fTheta
Double_t fTheta
Definition: CbmRichRecGeoPar.h:57
CbmRichRecGeoParPmt
Definition: CbmRichRecGeoPar.h:36
CbmRichGeoManager::InitMirror
void InitMirror()
Definition: CbmRichGeoManager.cxx:328
CbmRichGeoManager.h
CbmRichGeoManager::InitPmtCyl
void InitPmtCyl()
Definition: CbmRichGeoManager.cxx:95
CbmRichRecGeoPar::fMirrorTheta
Double_t fMirrorTheta
Definition: CbmRichRecGeoPar.h:235
CbmRichRecGeoParPmt::fY
Double_t fY
Definition: CbmRichRecGeoPar.h:62
CbmRichRecGeoParPmt::fPlaneX
Double_t fPlaneX
Definition: CbmRichRecGeoPar.h:66
CbmRichPmtPlaneMinMax::AddPoint
void AddPoint(Double_t x, Double_t y, Double_t z)
Definition: CbmRichRecGeoPar.h:256
CbmRichRecGeoPar
PMT parameters for the RICH geometry.
Definition: CbmRichRecGeoPar.h:87
CbmRichPmtPlaneMinMax::GetMeanX
Double_t GetMeanX()
Definition: CbmRichRecGeoPar.h:265
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
CbmRichGeoManager::InitPmt
void InitPmt()
Definition: CbmRichGeoManager.cxx:257
CbmRichPmtPlaneMinMax::GetMeanY
Double_t GetMeanY()
Definition: CbmRichRecGeoPar.h:267
CbmRichRecGeoParPmt::fX
Double_t fX
Definition: CbmRichRecGeoPar.h:61
CbmRichRecGeoParPmt::fZ
Double_t fZ
Definition: CbmRichRecGeoPar.h:63
shape
UInt_t shape
Definition: CbmMvdSensorDigiToHitTask.cxx:73
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichRecGeoPar::fMirrorR
Double_t fMirrorR
Definition: CbmRichRecGeoPar.h:233
CbmRichGeoManager::RotatePointTwoWings
void RotatePointTwoWings(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:424
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichRecGeoParPmt::fWidth
Double_t fWidth
Definition: CbmRichRecGeoPar.h:70
CbmRichRecGeoParPmt::fPmtPositionIndexX
Int_t fPmtPositionIndexX
Definition: CbmRichRecGeoPar.h:74
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmRichRecGeoPar::GetGeoRecPmtByBlockPathOrClosest
CbmRichRecGeoParPmt GetGeoRecPmtByBlockPathOrClosest(const std::string &path, TVector3 *pos)
Definition: CbmRichRecGeoPar.h:181
CbmRichGeoManager::fGP
CbmRichRecGeoPar * fGP
Definition: CbmRichGeoManager.h:23
CbmRichRecGeoPar::fPmtStripGap
Double_t fPmtStripGap
Definition: CbmRichRecGeoPar.h:225
CbmRichPmtPlaneMinMax::GetMeanZ
Double_t GetMeanZ()
Definition: CbmRichRecGeoPar.h:269
CbmRichRecGeoParPmt::fPlaneZ
Double_t fPlaneZ
Definition: CbmRichRecGeoPar.h:68
CbmRichGeoManager::RotatePointImpl
void RotatePointImpl(TVector3 *inPos, TVector3 *outPos, Double_t phi, Double_t theta, Double_t pmtX, Double_t pmtY, Double_t pmtZ)
Definition: CbmRichGeoManager.cxx:499