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