er  dev
ERTelescopeSetup.cxx
1 /********************************************************************************
2  * Copyright (C) Joint Institute for Nuclear Research *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence version 3 (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 
9 #include "ERTelescopeSetup.h"
10 
11 #include <ostream>
12 #include <iostream>
13 
14 #include "TGeoSphere.h"
15 #include "TError.h"
16 #include "TMath.h"
17 #include "TGeoManager.h"
18 #include "TGeoMatrix.h"
19 #include "TROOT.h"
20 #include <Riostream.h>
21 #include <TDOMParser.h>
22 #include <TXMLAttr.h>
23 #include <TXMLNode.h>
24 #include <TList.h>
25 
26 #include "FairRootManager.h"
27 #include "FairRunAna.h"
28 #include "FairRuntimeDb.h"
29 #include "FairGeoLoader.h"
30 #include "FairGeoMedium.h"
31 #include "FairGeoInterface.h"
32 #include "FairGeoBuilder.h"
33 #include "FairGeoMedia.h"
34 #include "FairLogger.h"
35 
36 ERTelescopeSetup* ERTelescopeSetup::fInstance = nullptr;
37 //--------------------------------------------------------------------------------------------------
38 ERTelescopeStrip::ERTelescopeStrip(Double_t globalX, Double_t globalY, Double_t globalZ,
39  Double_t localX, Double_t localY, Double_t localZ,
40  Double_t width)
41  : fGlobalX(globalX), fGlobalY(globalY), fGlobalZ(globalZ),
42  fLocalX(localX), fLocalY(localY), fLocalZ(localZ), fWidth(width)
43 {
44 }
45 //--------------------------------------------------------------------------------------------------
46 ERTelescopeStrip::ERTelescopeStrip(Double_t* globTrans, Double_t* localTrans, Double_t width)
47  : fGlobalX(globTrans[0]), fGlobalY(globTrans[1]), fGlobalZ(globTrans[2]),
48  fLocalX(localTrans[0]), fLocalY(localTrans[1]), fLocalZ(localTrans[2]), fWidth(width)
49 {
50 }
51 //--------------------------------------------------------------------------------------------------
52 ERTelescopeSetup* ERTelescopeSetup::Instance(){
53  if (fInstance == NULL) {
54  fInstance = new ERTelescopeSetup();
55  return fInstance;
56  }
57  else
58  return fInstance;
59 }
60 //--------------------------------------------------------------------------------------------------
61 Double_t ERTelescopeSetup::GetStripGlobalX(const TString& componentBranchName, Int_t stripNb) const {
62  return fStrips.at(componentBranchName).at(stripNb).fGlobalX;
63 }
64 //--------------------------------------------------------------------------------------------------
65 Double_t ERTelescopeSetup::GetStripGlobalY(const TString& componentBranchName, Int_t stripNb) const {
66  return fStrips.at(componentBranchName).at(stripNb).fGlobalY;
67 }
68 //--------------------------------------------------------------------------------------------------
69 Double_t ERTelescopeSetup::GetStripGlobalZ(const TString& componentBranchName, Int_t stripNb) const {
70  return fStrips.at(componentBranchName).at(stripNb).fGlobalZ;
71 }
72 //--------------------------------------------------------------------------------------------------
73 Double_t ERTelescopeSetup::GetStripLocalX(const TString& componentBranchName, Int_t stripNb) const {
74  return fStrips.at(componentBranchName).at(stripNb).fLocalX;
75 }
76 //--------------------------------------------------------------------------------------------------
77 Double_t ERTelescopeSetup::GetStripLocalY(const TString& componentBranchName, Int_t stripNb) const {
78  return fStrips.at(componentBranchName).at(stripNb).fLocalY;
79 }
80 //--------------------------------------------------------------------------------------------------
81 Double_t ERTelescopeSetup::GetStripLocalZ(const TString& componentBranchName, Int_t stripNb) const {
82  return fStrips.at(componentBranchName).at(stripNb).fLocalZ;
83 }
84 //--------------------------------------------------------------------------------------------------
85 Double_t ERTelescopeSetup::GetStripWidth(TString componentBranchName, Int_t stripNb) const {
86  return fStrips.at(componentBranchName).at(stripNb).fWidth;
87 }
88 //--------------------------------------------------------------------------------------------------
89 Double_t ERTelescopeSetup::GetStripPhi(const TString& componentBranchName, const Int_t stripNb) const {
90  if (fRStrips.find(componentBranchName) == fRStrips.end()) {
91  LOG(FATAL) << "R strips has not been filled for " << componentBranchName << FairLogger::endl;
92  }
93  if (stripNb >= fRStrips.at(componentBranchName).size()) {
94  LOG(FATAL) << "Wrong strip number " << stripNb << " for station " << componentBranchName << FairLogger::endl;
95  }
96  if (fRStrips.at(componentBranchName).at(stripNb).fPhi == -1) {
97  LOG(FATAL) << "Strip of station " << componentBranchName << " has not attribute phi" << FairLogger::endl;
98  }
99  return fRStrips.at(componentBranchName).at(stripNb).fPhi;
100 }
101 //--------------------------------------------------------------------------------------------------
102 Double_t ERTelescopeSetup::GetStripR(const TString& componentBranchName, const Int_t stripNb) const {
103  if (fRStrips.find(componentBranchName) == fRStrips.end()) {
104  LOG(FATAL) << "R strips has not been filled for " << componentBranchName << FairLogger::endl;
105  }
106  if (stripNb >= fRStrips.at(componentBranchName).size()) {
107  LOG(FATAL) << "Wrong strip number " << stripNb << " for station " << componentBranchName << FairLogger::endl;
108  }
109  if (fRStrips.at(componentBranchName).at(stripNb).fR == -1) {
110  LOG(FATAL) << "Strip of station " << componentBranchName << " has not attribute R" << FairLogger::endl;
111  }
112  return fRStrips.at(componentBranchName).at(stripNb).fR;
113 }
114 //--------------------------------------------------------------------------------------------------
115 TVector3 ERTelescopeSetup::
116 GetStripLocalPosition(const TString& componentBranchName, const unsigned int stripNb) const {
117  return TVector3(fStrips.at(componentBranchName).at(stripNb).fLocalX,
118  fStrips.at(componentBranchName).at(stripNb).fLocalY,
119  fStrips.at(componentBranchName).at(stripNb).fLocalZ);
120 }
121 //--------------------------------------------------------------------------------------------------
122 TVector3 ERTelescopeSetup::GetStationTranslation(const TString& componentBranchName) const {
123  return fStationGlobalToLocalMatrixies.at(componentBranchName).GetTranslation();
124 }
125 //--------------------------------------------------------------------------------------------------
126 ERTelescopeSetup::StationType ERTelescopeSetup::GetStationType(const TString& componentBranchName) const {
127  if (fStationTypes.find(componentBranchName) == fStationTypes.end())
128  LOG(FATAL) << "Station type has not been filled for " << componentBranchName << FairLogger::endl;
129  return fStationTypes.at(componentBranchName);
130 }
131 //--------------------------------------------------------------------------------------------------
132 TVector3 ERTelescopeSetup::ToStationCoordinateSystem (const TString& componentBranchName,
133  const TVector3& vectorInGlobalCS) const {
134  Double_t global[3], local[3];
135  for (int i(0); i < 3; i++)
136  global[i] = vectorInGlobalCS[i];
137  fStationGlobalToLocalMatrixies.at(componentBranchName).MasterToLocal(global, local);
138  return TVector3(local);
139 }
140 //--------------------------------------------------------------------------------------------------
141 TVector3 ERTelescopeSetup::ToGlobalCoordinateSystem(const TString& componentBranchName,
142  const TVector3& vectorInStationCS) const {
143  Double_t global[3], local[3];
144  for (int i(0); i < 3; i++)
145  local[i] = vectorInStationCS[i];
146  fStationGlobalToLocalMatrixies.at(componentBranchName).LocalToMaster(local, global);
147  return TVector3(global);
148 }
149 //--------------------------------------------------------------------------------------------------
150 TGeoHMatrix GetGlobalToLocalMatrix(const TString& path) {
151  TGeoIterator nextNode(gGeoManager->GetTopVolume());
152  while(nextNode()) {
153  TString nodePath;
154  nextNode.GetPath(nodePath);
155  if (nodePath == path) {
156  return (*static_cast<const TGeoHMatrix*>(nextNode.GetCurrentMatrix()));
157  }
158  }
159  LOG(FATAL) << "Path " << path << " not found in geometry" << FairLogger::endl;
160  return nullptr;
161 }
162 //--------------------------------------------------------------------------------------------------
163 void ERTelescopeSetup::GetTransInMotherNode (TGeoNode const* node, Double_t b[3]) {
164  memcpy(b, node->GetMatrix()->GetTranslation(), sizeof(double)*3);
165 }
166 // Написать несколько методов для работы с объемами:
167 // 1) Функция получения координаты в глобальной СК. (скорее всего функция в п.1)
168 // Проход от самого низкого дочернего объема до родительского внтри одной функции
169 // ------------------------------- -----------------------------------------------------------------
170 void ERTelescopeSetup::ReadGeoParamsFromParContainer() {
171  if (fGeometryInited)
172  return;
173  LOG(DEBUG) << "Loading Telescope setup from parameters database file." << FairLogger::endl;
174  if ( ! gGeoManager ) {
175  std::cerr << "ERTelescopeSetup: cannot initialise without TGeoManager!"<< std::endl;
176  }
177  gGeoManager->CdTop();
178 
179  TGeoNode* cave = gGeoManager->GetCurrentNode();
180  TGeoNode* qtelescope = NULL;
181  TGeoNode* qtelescopeDetector = NULL;
182  TGeoNode* qtelescopeStation = NULL;
183  const auto strip_width = [](TGeoNode* strip_node) {
184  const auto strip_shape = strip_node->GetVolume()->GetShape();
185  double low, high;
186  return TMath::Min(strip_shape->GetAxisRange(1, low, high),
187  strip_shape->GetAxisRange(2, low, high));
188  };
189  for (Int_t iNode = 0; iNode < cave->GetNdaughters(); iNode++) { // cycle by volumes in TOP
190  TString detectorName = cave->GetDaughter(iNode)->GetName();
191  if ( detectorName.Contains("Telescope", TString::kIgnoreCase) ) {
192  qtelescope = cave->GetDaughter(iNode);
193  for (Int_t iDetector = 0; iDetector < qtelescope->GetNdaughters(); iDetector++) { // cycle by subassemblies in QTelescope
194  qtelescopeDetector = qtelescope->GetDaughter(iDetector);
195  for (Int_t iStation = 0; iStation < qtelescopeDetector->GetNdaughters(); iStation++) { // cycle by components in station
196  qtelescopeStation = qtelescopeDetector->GetDaughter(iStation);
197  TString qtelescopeStationName = qtelescopeStation->GetName();
198  Double_t boxInStripTrans[3];
199  Double_t stripInStationTrans[3];
200  Double_t stripInDetectorTrans[3];
201  Double_t stripGlobTrans[3];
202  const bool is_r_station = qtelescopeStation->GetNdaughters() == 1
203  && TString(qtelescopeStation->GetDaughter(0)->GetName()).Contains("r_station");
204  TString digiBranchName = qtelescopeStationName;
205  digiBranchName.Remove(digiBranchName.Last('_'), digiBranchName.Length());
206  if (qtelescopeStationName.Contains("DoubleSi", TString::kIgnoreCase) ) {
207  TString firstStripArrayName, secondStripArrayName;
208  if (qtelescopeStationName.Contains("XY")) {
209  firstStripArrayName = digiBranchName + "_X";
210  secondStripArrayName = digiBranchName + "_Y";
211  } else {
212  firstStripArrayName = digiBranchName + "_Y";
213  secondStripArrayName = digiBranchName + "_X";
214  }
215  if (is_r_station) {
216  LOG(DEBUG) << "Read geometry info for Double R station corresponded to volume " << qtelescopeStationName
217  << " and branch name " << digiBranchName << FairLogger::endl;
218  auto* r_station = qtelescopeStation->GetDaughter(0);
219  FillRStrips(r_station, firstStripArrayName);
220  auto* any_segment = r_station->GetDaughter(0);
221  FillRStrips(any_segment, secondStripArrayName);
222  fStationTypes[firstStripArrayName] = StationType::RStation;
223  fStationTypes[secondStripArrayName] = StationType::RStation;
224  } else {
225  TGeoNode* doubleSiStrip;
226  Bool_t flagFirstStripReaded = kFALSE;
227  Int_t iDoubleSiStrip = 0;
228  for (; iDoubleSiStrip < qtelescopeStation->GetNdaughters(); iDoubleSiStrip++) {
229  doubleSiStrip = qtelescopeStation->GetDaughter(iDoubleSiStrip);
230  GetTransInMotherNode(doubleSiStrip, stripInStationTrans);
231  // Now we assume that QTelescope center in the global zero.
232  // It is neccesary to write clear function to find global coordinates
233  // by all forefathers nodes if possible. Maybe with some FairRoot methods
234  qtelescopeStation->LocalToMaster(stripInStationTrans, stripInDetectorTrans);
235  qtelescopeDetector->LocalToMaster(stripInDetectorTrans, stripGlobTrans);
236  fStrips[firstStripArrayName].emplace_back(stripGlobTrans, stripInStationTrans,
237  strip_width(doubleSiStrip));
238  LOG(DEBUG) << firstStripArrayName << " strip "
239  << fStrips[firstStripArrayName].size()-1 << " global coordinates: "
240  << stripGlobTrans[0] << ", "
241  << stripGlobTrans[1] << ", "
242  << stripGlobTrans[2] ;
243  LOG(DEBUG) << " | local coordinates: "
244  << boxInStripTrans[0] << ", "
245  << boxInStripTrans[1] << ", "
246  << boxInStripTrans[2] << FairLogger::endl;
247  TGeoNode* doubleSiBox;
248  Int_t iDoubleSiBox = 0;
249  if (!flagFirstStripReaded) {
250  for (; iDoubleSiBox < doubleSiStrip->GetNdaughters(); iDoubleSiBox++) {
251  Double_t siBoxLocalTrans[3];
252  doubleSiBox = doubleSiStrip->GetDaughter(iDoubleSiBox);
253  TString siBoxName = doubleSiBox->GetName();
254  GetTransInMotherNode(doubleSiBox, boxInStripTrans);
255  (qtelescopeStationName.Contains("XY")) ? stripInStationTrans[0] = 0
256  : stripInStationTrans[1] = 0;
257  doubleSiStrip->LocalToMaster(boxInStripTrans, stripInStationTrans);
258  qtelescopeStation->LocalToMaster(stripInStationTrans, stripInDetectorTrans);
259  qtelescopeDetector->LocalToMaster(stripInDetectorTrans, stripGlobTrans);
260  fStrips[secondStripArrayName].emplace_back(stripGlobTrans, boxInStripTrans,
261  strip_width(doubleSiBox));
262  LOG(DEBUG) << secondStripArrayName << " strip "
263  << fStrips[secondStripArrayName].size()-1 << " global coordinates: "
264  << stripGlobTrans[0] << ", "
265  << stripGlobTrans[1] << ", "
266  << stripGlobTrans[2];
267  LOG(DEBUG) << " | local coordinates: "
268  << boxInStripTrans[0] << ", "
269  << boxInStripTrans[1] << ", "
270  << boxInStripTrans[2] << FairLogger::endl;;
271 
272  }
273  flagFirstStripReaded = kTRUE;
274  }
275  }
276  fStationTypes[firstStripArrayName] = StationType::QStation;
277  fStationTypes[secondStripArrayName] = StationType::QStation;
278  }
279  TString stationPath;
280  stationPath.Form("cave/%s/%s/%s", qtelescope->GetName(), qtelescopeDetector->GetName(),
281  qtelescopeStationName.Data());
282  fStationGlobalToLocalMatrixies[firstStripArrayName] = GetGlobalToLocalMatrix(stationPath);
283  fStationGlobalToLocalMatrixies[secondStripArrayName] = fStationGlobalToLocalMatrixies[firstStripArrayName];
284  fStationGlobalToLocalMatrixies[firstStripArrayName].Print();
285  }
286  if (qtelescopeStationName.Contains("SingleSi", TString::kIgnoreCase) ) {
287  if (is_r_station) {
288  LOG(DEBUG) << "Read geometry info for R station corresponded to volume " << qtelescopeStationName
289  << " and branch name " << digiBranchName << FairLogger::endl;
290  FillRStrips(qtelescopeStation->GetDaughter(0), digiBranchName);
291  fStationTypes[digiBranchName] = StationType::RStation;
292  } else {
293  TGeoNode* singleSiStrip;
294  Int_t iSingleSiStrip = 0;
295  for (; iSingleSiStrip < qtelescopeStation->GetNdaughters(); iSingleSiStrip++) {
296  singleSiStrip = qtelescopeStation->GetDaughter(iSingleSiStrip);
297  GetTransInMotherNode(singleSiStrip, stripInStationTrans);
298  qtelescopeStation->LocalToMaster(stripInStationTrans, stripInDetectorTrans);
299  qtelescopeDetector->LocalToMaster(stripInDetectorTrans, stripGlobTrans);
300  LOG(DEBUG) << qtelescopeStationName << " strip "
301  << iSingleSiStrip << " global coordinates: "
302  << stripGlobTrans[0] << ", "
303  << stripGlobTrans[1] << ", "
304  << stripGlobTrans[2];
305  LOG(DEBUG) << " | local coordinates: "
306  << stripInStationTrans[0] << ", "
307  << stripInStationTrans[1] << ", "
308  << stripInStationTrans[2] << FairLogger::endl;
309  fStrips[digiBranchName].emplace_back(stripGlobTrans, stripInStationTrans,
310  strip_width(singleSiStrip));
311  }
312  fStationTypes[digiBranchName] = StationType::QStation;
313  }
314  TString stationPath;
315  stationPath.Form("cave/%s/%s/%s", qtelescope->GetName(), qtelescopeDetector->GetName(),
316  qtelescopeStationName.Data());
317  fStationGlobalToLocalMatrixies[digiBranchName] = GetGlobalToLocalMatrix(stationPath);
318  fStationGlobalToLocalMatrixies[digiBranchName].Print();
319  }
320  }
321  }
322  }
323  }
324  gGeoManager->CdTop();
325  fGeometryInited = true;
326 }
327 //--------------------------------------------------------------------------------------------------
328 void ERTelescopeSetup::FillRStrips(TGeoNode* r_station, const TString& branch_name) {
329  const bool is_phi_station = branch_name.EndsWith("_X");
330  std::cerr << 1 << std::endl;
331  for (int i_strip = 0; i_strip < r_station->GetNdaughters(); i_strip++) {
332  auto* strip = r_station->GetDaughter(i_strip);
333  if (is_phi_station) {
334  auto* combi_trans = dynamic_cast<TGeoCombiTrans*>(strip->GetMatrix());
335  if (!combi_trans) {
336  LOG(FATAL) << "Unexpected matrix type in R telescope station" << FairLogger::endl;
337  }
338  const auto phi = combi_trans->GetRotation()->GetPhiRotation();
339  fRStrips[branch_name].emplace_back(phi, -1.);
340  LOG(DEBUG) << branch_name << " strip " << i_strip << " phi = " << phi << FairLogger::endl;
341  } else {
342  auto shape = dynamic_cast<TGeoSphere*>(strip->GetVolume()->GetShape());
343  if (!shape) {
344  LOG(FATAL) << "Unexpected shape type in R telescope station" << FairLogger::endl;
345  }
346  const auto r_min = TMath::Tan(shape->GetTheta1() * TMath::DegToRad()) * shape->GetRmin();
347  const auto r_max = TMath::Tan(shape->GetTheta2() * TMath::DegToRad()) * shape->GetRmax();
348  const auto r = (r_max + r_min) / 2.;
349  fRStrips[branch_name].emplace_back(-1., r);
350  LOG(DEBUG) << branch_name << " strip " << i_strip << " R = " << r << FairLogger::endl;
351  }
352  }
353 }
354 //--------------------------------------------------------------------------------------------------
355 ClassImp(ERTelescopeSetup)
Generated on Tue Dec 7 2021 13:58:46 for er by   doxygen 1.8.11