CbmRoot
CbmMuchSegmentSector.cxx
Go to the documentation of this file.
1 
9 #include "CbmMuchSegmentSector.h"
10 
11 #include "CbmGeoMuchPar.h" // for CbmGeoMuchPar
12 #include "CbmMuchAddress.h" // for CbmMuchAddress
13 #include "CbmMuchLayer.h" // for CbmMuchLayer
14 #include "CbmMuchLayerSide.h" // for CbmMuchLayerSide
15 #include "CbmMuchModule.h" // for CbmMuchModule
16 #include "CbmMuchModuleGemRadial.h" // for CbmMuchModuleGemRadial
17 #include "CbmMuchSectorRadial.h" // for CbmMuchSectorRadial
18 #include "CbmMuchStation.h" // for CbmMuchStation
19 
20 #include <FairLogger.h> // for LOG
21 #include <FairRuntimeDb.h> // for FairRuntimeDb
22 #include <FairTask.h> // for FairTask, InitStatus, kSUCCESS
23 
24 #include <TArc.h> // for TArc
25 #include <TCanvas.h> // for TCanvas
26 #include <TFile.h> // for TFile, gFile
27 #include <TMath.h> // for DegToRad, ASin
28 #include <TObjArray.h> // for TObjArray
29 #include <TSystem.h> // for TSystem, gSystem
30 #include <TVector3.h> // for TVector3
31 
32 #include <math.h> // for sqrt
33 #include <stdio.h> // for printf, fprintf, fclose
34 #include <string.h> // for strcat, strlen, strncpy
35 #include <vector> // for vector
36 
37 using std::ifstream;
38 using std::string;
39 using std::vector;
40 
41 // ----- Default constructor -------------------------------------------
43  : FairTask()
44  , fGeoPar(nullptr)
45  , fNStations(0)
46  , fFlag(0)
47  , fStations(nullptr)
48  , fInputFileName()
49  , fDigiFileName()
50  , fNRegions()
51  , fRadii()
52  , fAngles()
53  , fSecLx()
54  , fSecLy()
55  , fNChannels()
56  , fNCols()
57  , fNRows()
58  , fDebug(0) {}
59 // -------------------------------------------------------------------------
60 
61 // ----- Standard constructor ------------------------------------------
63  TString digiFileName,
64  Int_t flag)
65  : FairTask()
66  , fGeoPar(nullptr)
67  , fNStations(0)
68  , fFlag(flag)
69  , fStations(nullptr)
70  , fInputFileName(inputFileName.Data())
71  , fDigiFileName(digiFileName.Data())
72  , fNRegions()
73  , fRadii()
74  , fAngles()
75  , fSecLx()
76  , fSecLy()
77  , fNChannels()
78  , fNCols()
79  , fNRows()
80  , fDebug(0) {}
81 // -------------------------------------------------------------------------
82 
83 // ----- Destructor ----------------------------------------------------
85 // -------------------------------------------------------------------------
86 
87 // ----- Private method SetParContainers --------------------------------
89  // Get runtime database
90  FairRuntimeDb* db = FairRuntimeDb::instance();
91  if (!db) Fatal("Init", "No runtime database");
92  fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
93 }
94 // -------------------------------------------------------------------------
95 
96 // ----- Private method Init ---------------------------------------------
98  printf("\n============================= Inputs segmentation parameters "
99  "================================\n");
100  ReadInputFile();
101  printf("\n==================================================================="
102  "============================\n");
103 
104  // Get MUCH geometry parameter container
106  // LOG(info)<<" Stations = "<<fStations->GetEntries()<<" "<<fNStations;
107  if (!fStations) Fatal("Init", "No input array of MUCH stations.");
108  if (fStations->GetEntries() != fNStations)
109  Fatal("Init", "Incorrect number of stations.");
110 
111  if (fDebug) {
112  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
113  printf("Station %i\n", iStation + 1);
114  for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion)
115  printf(" Region %i: fAngles=%4.1f\n",
116  iRegion + 1,
117  fAngles[iStation][iRegion]);
118  }
119  }
120 
121  // Segment MuCh
122  SegmentMuch();
123  return kSUCCESS;
124 }
125 // -------------------------------------------------------------------------
126 
127 // ----- Public method SegmentMuch --------------------------------------
129  for (Int_t iStation = 0; iStation < fStations->GetEntries(); ++iStation) {
130  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
131 
132  Int_t nLayers = station->GetNLayers();
133  for (Int_t iLayer = 0; iLayer < nLayers; ++iLayer) {
134  CbmMuchLayer* layer = station->GetLayer(iLayer);
135  if (!layer) Fatal("SegmentMuch", "Incomplete layers array.");
136  // Segment layer sides
137  printf("Layer=%d SideF Sectors=%i\n",
138  iLayer,
139  SegmentLayerSide(layer->GetSideF()));
140  printf("Layer=%d SideB Sectors=%i\n",
141  iLayer,
142  SegmentLayerSide(layer->GetSideB()));
143  }
144  printf("Station %i segmented\n", iStation + 1);
145  }
146 
147  // Save parameters
148  TFile* oldfile = gFile;
149  TFile* f = new TFile(fDigiFileName, "RECREATE");
150  fStations->Write("stations", 1);
151 
152  f->Close();
153  gFile = oldfile;
154 
156 }
157 // -------------------------------------------------------------------------
158 
159 // ----- Private method SegmentLayerSide --------------------------------
161  if (!layerSide) Fatal("SegmentLayerSide", "Incomplete layer sides array.");
162  Int_t nModules = layerSide->GetNModules();
163  // LOG(info)<<" Total Modules "<< nModules;
164  Int_t nSectors = 0;
165  for (Int_t iModule = 0; iModule < nModules; iModule++) {
166  CbmMuchModule* module = layerSide->GetModule(iModule);
167  if (module->GetDetectorType() != 3 && module->GetDetectorType() != 4)
168  continue;
169  // if(module->GetDetectorType()!=3) continue;
171  if (nModules > 0)
172  nSectors += SegmentModule(mod, true); // Module design
173  else
174  nSectors += SegmentModule(mod, false); // Monolithic design
175  }
176  return nSectors;
177 }
178 // -------------------------------------------------------------------------
179 
180 // ----- Private method SegmentSector -----------------------------------
182  Bool_t) {
183  Int_t detectorId = module->GetDetectorId();
184  Int_t iStation = CbmMuchAddress::GetStationIndex(detectorId);
185  Int_t iModule = CbmMuchAddress::GetModuleIndex(detectorId);
186  Int_t iLayer = CbmMuchAddress::GetLayerIndex(detectorId);
187  Int_t iSide = CbmMuchAddress::GetLayerSideIndex(detectorId);
188  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
189  Double_t rMin = 0.0, rMax = 0.0, r0 = 0.0;
190  Double_t phi0 = 0.0;
191  if (fFlag == 0) {
192  rMin = station->GetRmin();
193  rMax = station->GetRmax();
194  r0 = module->GetPosition().Perp();
195  phi0 = module->GetPosition().Phi();
196 
197  }
198 
199  else {
200  //============For mini cbm=================================
201  rMin = 18.4654; // Ajit+OS -> rMin for Mv2 module
202  rMax = 97.2187; // Ajit+OS -> rMax for Mv2 module
203  r0 = 58.46; // Ajit+OS -> r0 for Mv2 module
204  phi0 =
205  258.50
206  * TMath::
207  DegToRad(); // Ajit+OS -> rotaion angle for segmented geomtery such that we don't need to rotate point - only translation is needed
208  }
209  Double_t dx1 = module->GetDx1();
210  Double_t dx2 = module->GetDx2();
211  Double_t dy = module->GetDy();
212 
213  //LOG(info)<<" station # "<<iStation<<" Layer No. "<<iLayer<<" Side # "<<iSide<<" Module # "<<iModule<<" rmin "<<rMin<<" rmax "<<rMax<<" phi "<<phi0<<" r0 "<<r0<<" dx1 "<<dx1<<" dx2 "<<dx2<<" dy "<<dy;
214 
215  Double_t t = 2 * dy / (dx2 - dx1);
216  Double_t b = (r0 - dy) - t * dx1;
217  Double_t b2 = b * b;
218  Double_t bt = b * t;
219  Double_t a = t * t + 1;
220 
221  Double_t r1 = rMin;
222  Double_t r2 = rMin;
223  if (fDebug) printf("Debug: %i %i %i %i\n", iStation, iLayer, iSide, iModule);
224  Int_t iSector = 0;
225  for (Int_t i = 0; i < fNRegions[iStation]; i++) {
226  Double_t angle = fAngles[iStation][i] * TMath::DegToRad();
227  while (r2 < fRadii[iStation][i] && r2 < rMax) {
228  r2 = r1 + r1 * angle;
229  Double_t dphi = TMath::ASin((-bt + sqrt(a * r2 * r2 - b2)) / a / r2);
230  // Int_t nPads = 2*Int_t(phiMax/angle)+2;
231  // if (iStation==0 && iLayer==0 && iSide==0 && iModule==0) printf("Sector: %f-%f %i phiMax=%f\n",r1,r2,nPads,phiMax);
232  // LOG(info)<<" sector "<<iSector<< " pad size "<<r2-r1;
233  module->AddSector(new CbmMuchSectorRadial(
234  detectorId, iSector, r1, r2, phi0 - dphi, phi0 + dphi));
235  //LOG(info)<<"r1 "<<r1<<" r2 "<<r2;
236  r1 = r2;
237  iSector++;
238  }
239  }
240  if (fDebug) printf(" Sectors = %i\n", iSector);
241  return iSector;
242 }
243 // -------------------------------------------------------------------------
244 
245 
246 // -------------------------------------------------------------------------
248  ifstream infile;
249  infile.open(fInputFileName);
250  if (!infile) { Fatal("ReadInputFile", "Error: Cannot open the input file."); }
251 
252  Int_t index;
253  string str;
254  vector<string> strs;
255 
256  // Set number of stations
257  OmitDummyLines(infile, str);
258  strs = Split(str, ' ');
259  index = strs.size() - 1;
260  Int_t nStations;
261  StrToNum(strs.at(index), nStations);
262  fNStations = nStations;
263  printf("Number of stations: \t%i", fNStations);
264 
265  // Set number of regions
266  OmitDummyLines(infile, str);
267  strs = Split(str, ' ');
268  printf("\nNumber of regions: ");
269  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
270  index = strs.size() - fNStations + iStation;
271  Int_t nRegions;
272  StrToNum(strs.at(index), nRegions);
273  fNRegions[iStation] = nRegions;
274  fRadii[iStation].resize(nRegions);
275  fAngles[iStation].resize(nRegions);
276  printf("\t%i", fNRegions[iStation]);
277  }
278 
279  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
280  // Set region radii
281  OmitDummyLines(infile, str);
282  strs = Split(str, ' ');
283  printf("\nStation %i", iStation + 1);
284  printf("\n Region radii [cm]: ");
285  for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
286  index = strs.size() - fNRegions[iStation] + iRegion;
287  StrToNum(strs.at(index), fRadii[iStation][iRegion]);
288  printf("\t%4.2f", fRadii[iStation][iRegion]);
289  }
290  // Set pad angles
291  OmitDummyLines(infile, str);
292  strs = Split(str, ' ');
293  printf("\n Pad angles [deg]: ");
294  for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
295  index = strs.size() - fNRegions[iStation] + iRegion;
296  StrToNum(strs.at(index), fAngles[iStation][iRegion]);
297  printf("\t%4.2f", fAngles[iStation][iRegion]);
298  }
299  }
300  printf("\n");
301  infile.close();
302 }
303 
305  // Change file extension
306  char txtfile[100];
307  Int_t length = strlen(fDigiFileName);
308  Int_t iChar;
309  for (iChar = length - 1; iChar >= 0; --iChar) {
310  if (fDigiFileName[iChar] == '.') break;
311  }
312  strncpy(txtfile, fDigiFileName, iChar + 1);
313  txtfile[iChar + 1] = '\0';
314  strcat(txtfile, "txt");
315 
316  FILE* outfile;
317  outfile = fopen(txtfile, "w");
318  /*
319  Int_t colors[] = {kGreen, kMagenta, kCyan, kRed, kBlue, kYellow, kTeal,
320  kPink, kAzure, kOrange, kViolet, kSpring,
321  kGreen+2, kMagenta+2, kCyan+2, kRed+2, kBlue+2, kYellow+2, kTeal+2,
322  kPink+2, kAzure+2, kOrange+2, kViolet+2, kSpring+2,
323  kGreen+4, kMagenta+4, kCyan+4, kRed+4, kBlue+4, kYellow+4, kTeal+4,
324  kPink+4, kAzure+4, kOrange+4, kViolet+4, kSpring+4};
325 */
326  for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
327  fprintf(outfile,
328  "=================================================================="
329  "=========\n");
330  fprintf(outfile, "Station %i\n", iStation + 1);
331  fprintf(outfile,
332  "Sector size, cm Sector position, cm Number of pads Side "
333  "Pad size, cm\n");
334  fprintf(outfile,
335  "------------------------------------------------------------------"
336  "----------\n");
337  TCanvas* c1 = new TCanvas(Form("station%i", iStation + 1),
338  Form("station%i", iStation + 1),
339  1000,
340  1000);
341  c1->SetFillColor(0);
342  c1->Range(-200, -200, 200, 200);
343  c1->Range(-270, -270, 270, 270);
344  c1->Range(-50, -50, 50, 50);
345  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
346  CbmMuchLayer* layer = station->GetLayer(0);
347  for (Int_t iSide = 1; iSide >= 0; iSide--) {
348  CbmMuchLayerSide* layerSide = layer->GetSide(iSide);
349  for (Int_t iModule = 0; iModule < layerSide->GetNModules(); ++iModule) {
350  CbmMuchModule* mod = layerSide->GetModule(iModule);
351  mod->SetFillStyle(0);
352  // mod->Draw();
353  if (mod->GetDetectorType() != 3 && mod->GetDetectorType() != 4)
354  continue;
355  LOG(info) << "Det SEgmentation: " << mod->GetDetectorType()
356  << " Station " << iStation;
357  // if(mod->GetDetectorType() != 3) continue;
359  for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
360  CbmMuchSectorRadial* sector =
361  (CbmMuchSectorRadial*) module->GetSectorByIndex(iSector);
362  sector->AddPads();
363  sector->DrawPads();
364  } // sectors
365  } // modules
366  } // sides
367 
368  // Draw a hole
369  TArc* holeArc = new TArc(0., 0., station->GetRmin());
370  // TArc* holeArc = new TArc(0.,0.,rMin);
371  holeArc->Draw();
372 
373  for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
374  TArc* arc = new TArc(0., 0., fRadii[iStation].at(iRegion));
375  arc->SetLineColor(kBlack);
376  arc->SetLineWidth(2);
377  arc->SetFillStyle(0);
378  arc->Draw();
379  }
380 
381  if (fDebug) {
382  c1->Print(Form(
383  "%s/station%i.eps", gSystem->DirName(fDigiFileName), iStation + 1));
384  c1->Print(Form(
385  "%s/station%i.png", gSystem->DirName(fDigiFileName), iStation + 1));
386  }
387 
388  } //stations
389  fclose(outfile);
390 }
391 
392 
393 //double x0 = -a*c/(a*a+b*b), y0 = -b*c/(a*a+b*b);
394 //if (c*c > r*r*(a*a+b*b)+EPS)
395 // puts ("no points");
396 //else {
397 // double d = r*r - c*c/(a*a+b*b);
398 // double mult = sqrt (d / (a*a+b*b));
399 // double ax,ay,bx,by;
400 // ax = x0 + b * mult;
401 // bx = x0 - b * mult;
402 // ay = y0 - a * mult;
403 // by = y0 + a * mult;
404 //}
405 
CbmMuchSegmentSector::DrawSegmentation
void DrawSegmentation()
Definition: CbmMuchSegmentSector.cxx:304
CbmMuchLayerSide.h
CbmMuchModule.h
CbmMuchSegmentSector::fAngles
std::map< Int_t, std::vector< Double_t > > fAngles
Definition: CbmMuchSegmentSector.h:68
CbmMuchSegmentSector::StrToNum
void StrToNum(std::string &str, T &number)
Definition: CbmMuchSegmentSector.h:150
CbmMuchModuleGem::AddSector
void AddSector(CbmMuchSector *sector)
Definition: CbmMuchModuleGem.h:67
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
CbmMuchSegmentSector::fNStations
Int_t fNStations
Definition: CbmMuchSegmentSector.h:58
CbmMuchStation::GetLayer
CbmMuchLayer * GetLayer(Int_t iLayer) const
Definition: CbmMuchStation.h:58
CbmMuchSectorRadial
Definition: CbmMuchSectorRadial.h:19
CbmMuchStation
Definition: CbmMuchStation.h:22
CbmMuchSegmentSector::fRadii
std::map< Int_t, std::vector< Double_t > > fRadii
Definition: CbmMuchSegmentSector.h:66
CbmMuchSegmentSector
Definition: CbmMuchSegmentSector.h:32
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmMuchSegmentSector::SegmentMuch
void SegmentMuch()
Definition: CbmMuchSegmentSector.cxx:128
CbmMuchSegmentSector::fFlag
Int_t fFlag
Definition: CbmMuchSegmentSector.h:59
CbmMuchSegmentSector::fDigiFileName
TString fDigiFileName
Definition: CbmMuchSegmentSector.h:62
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmMuchSegmentSector::fStations
TObjArray * fStations
Definition: CbmMuchSegmentSector.h:60
CbmGeoMuchPar::GetStations
TObjArray * GetStations()
Definition: CbmGeoMuchPar.h:40
CbmGeoMuchPar
Definition: CbmGeoMuchPar.h:25
CbmMuchLayer::GetSide
CbmMuchLayerSide * GetSide(Bool_t side)
Definition: CbmMuchLayer.h:49
CbmMuchModuleGemRadial.h
CbmMuchModule::GetPosition
TVector3 GetPosition() const
Definition: CbmMuchModule.h:51
CbmMuchAddress::GetModuleIndex
static Int_t GetModuleIndex(Int_t address)
Definition: CbmMuchAddress.h:112
CbmMuchSegmentSector::Split
std::vector< std::string > & Split(const std::string &s, char delim, std::vector< std::string > &elems)
Definition: CbmMuchSegmentSector.h:134
CbmMuchSegmentSector::OmitDummyLines
void OmitDummyLines(std::ifstream &infile, std::string &str)
Definition: CbmMuchSegmentSector.h:127
CbmMuchModuleGemRadial
Definition: CbmMuchModuleGemRadial.h:21
CbmMuchSegmentSector::ReadInputFile
void ReadInputFile()
Definition: CbmMuchSegmentSector.cxx:247
CbmGeoMuchPar.h
CbmMuchStation::GetRmax
Double_t GetRmax() const
Definition: CbmMuchStation.h:52
CbmMuchLayerSide::GetModule
CbmMuchModule * GetModule(Int_t iModule) const
Definition: CbmMuchLayerSide.h:52
CbmMuchSegmentSector::fGeoPar
CbmGeoMuchPar * fGeoPar
Definition: CbmMuchSegmentSector.h:57
CbmMuchSegmentSector::fDebug
Int_t fDebug
Definition: CbmMuchSegmentSector.h:79
CbmMuchSectorRadial::AddPads
void AddPads()
Definition: CbmMuchSectorRadial.cxx:69
CbmMuchSegmentSector::fInputFileName
TString fInputFileName
Definition: CbmMuchSegmentSector.h:61
CbmMuchLayerSide
Definition: CbmMuchLayerSide.h:22
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmMuchModule
Definition: CbmMuchModule.h:24
CbmMuchSectorRadial.h
CbmMuchLayer.h
CbmMuchAddress::GetLayerIndex
static Int_t GetLayerIndex(Int_t address)
Definition: CbmMuchAddress.h:106
CbmMuchSegmentSector::SegmentModule
Int_t SegmentModule(CbmMuchModuleGemRadial *module, Bool_t useModuleDesign)
Definition: CbmMuchSegmentSector.cxx:181
CbmMuchAddress::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchAddress.h:103
CbmMuchModule::GetDetectorId
Int_t GetDetectorId() const
Definition: CbmMuchModule.h:48
CbmMuchSectorRadial::DrawPads
void DrawPads()
Definition: CbmMuchSectorRadial.cxx:87
CbmMuchModule::GetDetectorType
Int_t GetDetectorType() const
Definition: CbmMuchModule.h:52
CbmMuchModuleGemRadial::GetDx1
Double_t GetDx1() const
Definition: CbmMuchModuleGemRadial.h:37
CbmMuchAddress.h
CbmMuchSegmentSector::SetParContainers
virtual void SetParContainers()
Definition: CbmMuchSegmentSector.cxx:88
CbmMuchSegmentSector::SegmentLayerSide
Int_t SegmentLayerSide(CbmMuchLayerSide *layerSide)
Definition: CbmMuchSegmentSector.cxx:160
CbmMuchSegmentSector.h
CbmMuchStation::GetNLayers
Int_t GetNLayers() const
Definition: CbmMuchStation.h:50
CbmMuchLayer::GetSideF
CbmMuchLayerSide * GetSideF()
Definition: CbmMuchLayer.h:47
CbmMuchLayer
Definition: CbmMuchLayer.h:21
CbmMuchLayer::GetSideB
CbmMuchLayerSide * GetSideB()
Definition: CbmMuchLayer.h:48
CbmMuchModuleGemRadial::GetDy
Double_t GetDy() const
Definition: CbmMuchModuleGemRadial.h:39
CbmMuchStation::GetRmin
Double_t GetRmin() const
Definition: CbmMuchStation.h:51
CbmMuchStation.h
CbmMuchModuleGem::GetSectorByIndex
CbmMuchSector * GetSectorByIndex(Int_t iSector)
Definition: CbmMuchModuleGem.h:54
CbmMuchSegmentSector::~CbmMuchSegmentSector
virtual ~CbmMuchSegmentSector()
Definition: CbmMuchSegmentSector.cxx:84
CbmMuchModuleGem::GetNSectors
Int_t GetNSectors() const
Definition: CbmMuchModuleGem.h:57
CbmMuchAddress::GetLayerSideIndex
static Int_t GetLayerSideIndex(Int_t address)
Definition: CbmMuchAddress.h:109
CbmMuchModuleGemRadial::GetDx2
Double_t GetDx2() const
Definition: CbmMuchModuleGemRadial.h:38
CbmMuchLayerSide::GetNModules
Int_t GetNModules() const
Definition: CbmMuchLayerSide.h:47
CbmMuchSegmentSector::Init
virtual InitStatus Init()
Definition: CbmMuchSegmentSector.cxx:97
CbmMuchSegmentSector::fNRegions
std::map< Int_t, Int_t > fNRegions
Definition: CbmMuchSegmentSector.h:64
CbmMuchSegmentSector::CbmMuchSegmentSector
CbmMuchSegmentSector()
Definition: CbmMuchSegmentSector.cxx:42