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