CbmRoot
Loading...
Searching...
No Matches
CbmMuchSegmentManual.cxx
Go to the documentation of this file.
1/* Copyright (C) 2009-2021 St. Petersburg Polytechnic University, St. Petersburg
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Mikhail Ryzhinskiy [committer] */
4
14
15#include "CbmGeoMuchPar.h" // for CbmGeoMuchPar
16#include "CbmMuchAddress.h" // for CbmMuchAddress
17#include "CbmMuchLayer.h" // for CbmMuchLayer
18#include "CbmMuchLayerSide.h" // for CbmMuchLayerSide
19#include "CbmMuchModule.h" // for CbmMuchModule
20#include "CbmMuchModuleGemRectangular.h" // for CbmMuchModuleGemRectangular
21#include "CbmMuchSectorRectangular.h" // for CbmMuchSectorRectangular
22#include "CbmMuchStation.h" // for CbmMuchStation
23
24#include <FairRuntimeDb.h> // for FairRuntimeDb
25#include <FairTask.h> // for FairTask, InitStatus, kSUCCESS
26
27#include <TArc.h> // for TArc
28#include <TCanvas.h> // for TCanvas
29#include <TFile.h> // for TFile, gFile
30#include <TMath.h> // for Sqrt, Power, Log2
31#include <TMathBase.h> // for Abs, Max, Min
32#include <TObjArray.h> // for TObjArray
33#include <TSystem.h> // for TSystem, gSystem
34#include <TVector3.h> // for TVector3
35
36#include <cstdio> // for printf, fprintf, fclose
37#include <cstring> // for strcat, strlen, strncpy
38#include <limits> // for numeric_limits
39
40using std::ifstream;
41using std::string;
42using std::vector;
43
44// ----- Default constructor -------------------------------------------
46 : FairTask()
47 , fGeoPar(nullptr)
48 , fNStations(0)
49 , fStations(nullptr)
50 , fInputFileName((char*) "much.digi.root")
52 , fNRegions()
53 , fRadii()
54 , fSecLx()
55 , fSecLy()
56 , fNChannels()
57 , fNCols()
58 , fNRows()
59 , fDebug(0)
60{
61}
62// -------------------------------------------------------------------------
63
64// ----- Standard constructor ------------------------------------------
65CbmMuchSegmentManual::CbmMuchSegmentManual(char* inputFileName, char* digiFileName)
66 : FairTask()
67 , fGeoPar(nullptr)
68 , fNStations(0)
69 , fStations(nullptr)
70 , fInputFileName(inputFileName)
71 , fDigiFileName(digiFileName)
72 , fNRegions()
73 , fRadii()
74 , fSecLx()
75 , fSecLy()
76 , fNChannels()
77 , fNCols()
78 , fNRows()
79 , fDebug(0)
80{
81}
82// -------------------------------------------------------------------------
83
84// ----- Destructor ----------------------------------------------------
86// -------------------------------------------------------------------------
87
88// ----- Public method SetNRegions --------------------------------------
90{
91 if (iStation < 0 || iStation >= fNStations) Fatal("SetNRegions", "iStation is out of range.");
92 fNRegions[iStation] = nRegions;
93 fRadii[iStation].resize(nRegions);
94 fSecLx[iStation].resize(nRegions);
95 fSecLy[iStation].resize(nRegions);
96 fNCols[iStation].resize(nRegions);
97 fNRows[iStation].resize(nRegions);
98
99 if (fDebug) { printf("Station %i has %i regions\n", iStation + 1, nRegions); }
100
101 // Deal with channels more universally
102 Int_t n = Int_t(TMath::Log2(fNChannels[iStation]) + 1e-2);
103 Int_t nChans = Int_t(TMath::Power(2, n) + 1e-2);
104 if (nChans != fNChannels[iStation]) Fatal("Init", "Number of channels should be equal to two with integer power.");
105 Int_t nPower = n / 2;
106 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
107 fNCols[iStation].at(iRegion) = (Int_t) TMath::Power(2, nPower);
108 fNRows[iStation].at(iRegion) = n % 2 != 0 ? (Int_t) TMath::Power(2, nPower + 1) : fNCols[iStation].at(iRegion);
109
110 if (fDebug) {
111 printf("Region %i has %i columns and %i rows per sector\n", iRegion + 1, fNCols[iStation].at(iRegion),
112 fNRows[iStation].at(iRegion));
113 }
114 }
115}
116// -------------------------------------------------------------------------
117
118// ----- Public method SetNChannels --------------------------------------
120{
121 if (iStation < 0 || iStation >= fNStations) Fatal("SetNChannels", "iStation is out of range.");
122
123 fNChannels[iStation] = nChannels;
124
125 if (fDebug) { printf("Station %i has %i channels per sector\n", iStation + 1, nChannels); }
126}
127// -------------------------------------------------------------------------
128
129// ----- Public method SetRegionRadius ----------------------------------
130void CbmMuchSegmentManual::SetRegionRadius(Int_t iStation, Int_t iRegion, Double_t radius)
131{
132 if (iStation < 0 || iStation >= fNStations) Fatal("SetRegionRadius", "iStation is out of range.");
133 if (iRegion < 0 || iRegion >= fNRegions[iStation]) Fatal("SetRegionRadius", "iRegion is out of range.");
134 fRadii[iStation].at(iRegion) = radius;
135 if (fDebug) {
136 printf("Radius of the Region %i of station %i is %4.2f cm\n", iRegion + 1, iStation + 1,
137 fRadii[iStation].at(iRegion));
138 }
139}
140// -------------------------------------------------------------------------
141
142// ----- Public method SetSigma -----------------------------------------
143void CbmMuchSegmentManual::SetSigma(Int_t iStation, Int_t iRegion, Double_t sigmaX, Double_t sigmaY)
144{
145 if (iStation < 0 || iStation >= fNStations) Fatal("SetSigma", "iStation is out of range.");
146 if (iRegion < 0 || iRegion >= fNRegions[iStation]) Fatal("SetSigma", "iRegion is out of range.");
147
148 Double_t secLx = fSecLx[iStation].at(iRegion) = fNCols[iStation].at(iRegion) * TMath::Sqrt(12.) * sigmaX;
149 Double_t secLy = fSecLy[iStation].at(iRegion) = fNRows[iStation].at(iRegion) * TMath::Sqrt(12.) * sigmaY;
150 if (TMath::Abs(secLy - secLx) < 1e-5) return;
151 if (TMath::Abs(secLy / secLx - 2) > 1e-5) {
152 if (secLy > 2 * secLx) {
153 fNCols[iStation].at(iRegion) *= 2;
154 fNRows[iStation].at(iRegion) /= 2;
155 }
156 else {
157 fNCols[iStation].at(iRegion) /= 2;
158 fNRows[iStation].at(iRegion) *= 2;
159 }
160 SetSigma(iStation, iRegion, sigmaX, sigmaY);
161 }
162}
163// -------------------------------------------------------------------------
164
165// ----- Public method SetPadSize ----------------------------------------
166void CbmMuchSegmentManual::SetPadSize(Int_t iStation, Int_t iRegion, Double_t padLx, Double_t padLy)
167{
168 if (iStation < 0 || iStation >= fNStations) Fatal("SetPadSize", "iStation is out of range.");
169 if (iRegion < 0 || iRegion >= fNRegions[iStation]) Fatal("SetPadSize", "iRegion is out of range.");
170
171 Double_t secLx = fSecLx[iStation].at(iRegion) = fNCols[iStation].at(iRegion) * padLx;
172 Double_t secLy = fSecLy[iStation].at(iRegion) = fNRows[iStation].at(iRegion) * padLy;
173 if (TMath::Abs(secLy - secLx) < 1e-5) return;
174 if (TMath::Abs(secLy / secLx - 2) > 1e-5) {
175 if (secLy > 2 * secLx) {
176 fNCols[iStation].at(iRegion) *= 2;
177 fNRows[iStation].at(iRegion) /= 2;
178 }
179 else {
180 fNCols[iStation].at(iRegion) /= 2;
181 fNRows[iStation].at(iRegion) *= 2;
182 }
183 SetPadSize(iStation, iRegion, padLx, padLy);
184 }
185}
186// -------------------------------------------------------------------------
187
188// ----- Private method SetParContainers --------------------------------
190{
191 // Get runtime database
192 FairRuntimeDb* db = FairRuntimeDb::instance();
193 if (!db) Fatal("Init", "No runtime database");
194 fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
195}
196// -------------------------------------------------------------------------
197
198// ----- Private method Init ---------------------------------------------
200{
201 printf("\n============================= Inputs segmentation parameters "
202 "================================\n");
204 printf("\n==================================================================="
205 "============================\n");
206
207
208 // Get MUCH geometry parameter container
209 fStations = fGeoPar->GetStations();
210 if (!fStations) Fatal("Init", "No input array of MUCH stations.");
211 if (fStations->GetEntriesFast() != fNStations) Fatal("Init", "Incorrect number of stations.");
212
213 if (fDebug) {
214 for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
215 printf("Station %i\n", iStation + 1);
216 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
217 printf("Region %i: fSecLx = %f fSecLy = %f\n", iRegion + 1, fSecLx[iStation].at(iRegion),
218 fSecLy[iStation].at(iRegion));
219 printf("fNCols = %i fNRows = %i\n", fNCols[iStation].at(iRegion), fNRows[iStation].at(iRegion));
220 }
221 }
222 }
223 // Segment MuCh
224 SegmentMuch();
225 return kSUCCESS;
226}
227// -------------------------------------------------------------------------
228
229// ----- Public method SegmentMuch --------------------------------------
231{
232 for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
233 CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
234
235 Int_t nLayers = station->GetNLayers();
236 for (Int_t iLayer = 0; iLayer < nLayers; ++iLayer) {
237 CbmMuchLayer* layer = station->GetLayer(iLayer);
238 if (!layer) Fatal("SegmentMuch", "Incomplete layers array.");
239 // Segment layer sides
240 SegmentLayerSide(layer->GetSideF());
241 SegmentLayerSide(layer->GetSideB());
242 }
243 printf("Station %i segmented\n", iStation + 1);
244 }
245
247 TFile* oldFile = gFile;
248 TDirectory* oldDir = gDirectory;
249
250 // Save parameters
251 TFile* f = new TFile(fDigiFileName, "RECREATE");
252 fStations->Write("stations", 1);
253
254 f->Close();
256 gFile = oldFile;
257 gDirectory = oldDir;
258
259 // Draw colored stations
261
262 // Print some output info
263 Print();
264}
265// -------------------------------------------------------------------------
266
267// ----- Private method SegmentLayerSide --------------------------------
269{
270 if (!layerSide) Fatal("SegmentLayerSide", "Incomplete layer sides array.");
271 Int_t nModules = layerSide->GetNModules();
272 for (Int_t iModule = 0; iModule < nModules; iModule++) {
273 CbmMuchModule* module = layerSide->GetModule(iModule);
274 if (module->GetDetectorType() != 1) continue;
276 if (nModules > 1) SegmentModule(mod, true); // Module design
277 else
278 SegmentModule(mod, false); // Monolithic design
279 }
280}
281// -------------------------------------------------------------------------
282
283// ----- Private method SegmentSector -----------------------------------
285{
286 Int_t detectorId = module->GetDetectorId();
287 Int_t iStation = CbmMuchAddress::GetStationIndex(detectorId);
288 // CbmMuchStation* station = (CbmMuchStation*)fStations->At(iStation);
289 Int_t iRegion = -1;
290 Double_t secMaxLx = GetSectorMaxSize(module, "Width", iRegion);
291 Double_t secMaxLy = GetSectorMaxSize(module, "Length", iRegion);
292 Double_t padMaxLx = GetPadMaxSize(module, "Width");
293 Double_t padMaxLy = GetPadMaxSize(module, "Length");
294 TVector3 size = module->GetSize();
295 Double_t modLx = size.X();
296 Double_t modLy = size.Y();
297 Double_t modLz = size.Z();
298 TVector3 position = module->GetPosition();
299 Double_t modX = position.X();
300 Double_t modY = position.Y();
301 Double_t modZ = position.Z();
302
303 Int_t nCols = Int_t(modLx / secMaxLx);
304 Int_t nRows = Int_t(modLy / secMaxLy);
305 Int_t nX = modX < 0 ? -1 : 1;
306 Int_t nY = modY < 0 ? -1 : 1;
307
308 Double_t secX, secY; //, secLx, secLy;
309 TVector3 secSize, secPosition;
310 Int_t iSector;
311 for (Int_t iCol = 0; iCol < nCols; ++iCol) {
312 secX = useModuleDesign ? modX + nX * ((iCol + 0.5) * secMaxLx - modLx / 2.)
313 : secMaxLx * nCols / 2. - (0.5 + iCol) * secMaxLx;
314 for (Int_t iRow = 0; iRow < nRows; ++iRow) {
315 secY = useModuleDesign ? modY + nY * ((iRow + 0.5) * secMaxLy - modLy / 2.)
316 : secMaxLy * nRows / 2. - (0.5 + iRow) * secMaxLy;
317 iSector = module->GetNSectors();
318 secSize.SetXYZ(secMaxLx, secMaxLy, modLz);
319 secPosition.SetXYZ(secX, secY, modZ);
320 SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, secPosition, secSize,
321 fNCols[iStation].at(iRegion), fNRows[iStation].at(iRegion)));
322 }
323 }
324
325 // Process incomplete sectors
326 Int_t nPadRows, nPadCols;
327 Int_t n = useModuleDesign ? 1 : 2;
328 Double_t ly = useModuleDesign ? modLy - secMaxLy * nRows : (modLy - secMaxLy * nRows) / 2.;
329 if (ly > padMaxLy) {
330 nPadRows = Int_t(ly / padMaxLy);
331 for (Int_t i = 0; i < n; ++i) {
332 secY = useModuleDesign ? modY + nY * (modLy / 2. - ly / 2.) : TMath::Power(-1, i) * (modLy / 2. - ly / 2.);
333 for (Int_t iCol = 0; iCol < nCols; ++iCol) {
334 secX = useModuleDesign ? modX + nX * ((iCol + 0.5) * secMaxLx - modLx / 2.)
335 : secMaxLx * nCols / 2. - (0.5 + iCol) * secMaxLx;
336 iSector = module->GetNSectors();
337 secSize.SetXYZ(secMaxLx, ly, modLz);
338 secPosition.SetXYZ(secX, secY, modZ);
339 SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, secPosition, secSize,
340 fNCols[iStation].at(iRegion), nPadRows));
341 }
342 }
343 }
344
345 Double_t lx = useModuleDesign ? modLx - secMaxLx * nCols : (modLx - secMaxLx * nCols) / 2.;
346 if (lx > padMaxLx) {
347 nPadCols = (Int_t)(lx / padMaxLx);
348 for (Int_t i = 0; i < n; ++i) {
349 secX = useModuleDesign ? modX + nX * (modLx / 2. - lx / 2.) : TMath::Power(-1, i) * (modLx / 2. - lx / 2.);
350 for (Int_t iRow = 0; iRow < nRows; ++iRow) {
351 secY = useModuleDesign ? modY + nY * ((iRow + 0.5) * secMaxLy - modLy / 2.)
352 : secMaxLy * nRows / 2. - (0.5 + iRow) * secMaxLy;
353 iSector = module->GetNSectors();
354 secSize.SetXYZ(lx, secMaxLy, modLz);
355 secPosition.SetXYZ(secX, secY, modZ);
356 SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, secPosition, secSize, nPadCols,
357 fNRows[iStation].at(iRegion)));
358 }
359 }
360 }
361
362 if (lx > padMaxLx && ly > padMaxLy) {
363 nPadCols = (Int_t)(lx / padMaxLx);
364 nPadRows = (Int_t)(ly / padMaxLy);
365 secX = modX + nX * (modLx / 2. - lx / 2.);
366 secY = modY + nY * (modLy / 2. - ly / 2.);
367 iSector = module->GetNSectors();
368 secSize.SetXYZ(lx, ly, modLz);
369 secPosition.SetXYZ(secX, secY, modZ);
370 SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, secPosition, secSize, nPadCols, nPadRows));
371 }
372}
373// -------------------------------------------------------------------------
374
375// ----- Private method SegmentSector -----------------------------------
377{
378 TVector3 secSize = sector->GetSize();
379 TVector3 secPosition = sector->GetPosition();
380 Int_t detectorId = module->GetDetectorId();
381 Int_t iStation = CbmMuchAddress::GetStationIndex(detectorId);
382 Int_t iSector = module->GetNSectors();
383 Double_t secLx = secSize.X();
384 Double_t secLy = secSize.Y();
385 Double_t secLz = secSize.Z();
386 Double_t padLx = sector->GetPadDx();
387 Double_t padLy = sector->GetPadDy();
388 Double_t secX = secPosition.X();
389 Double_t secY = secPosition.Y();
390 Double_t secZ = secPosition.Z();
391 Bool_t isIncomplete = IsIncompleteSector(sector);
392 Int_t nCols = sector->GetPadNx();
393 Int_t nRows = sector->GetPadNy();
394 Int_t nX = secX < 0 ? -1 : 1;
395 Int_t nY = secY < 0 ? -1 : 1;
396
397 Int_t iRegion = -1;
398 Bool_t resultX = ShouldSegment(sector, "X", iRegion);
399 Bool_t resultY = ShouldSegment(sector, "Y", iRegion);
400
401 if (!resultX && !resultY) {
402 delete sector;
403
404 CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
405 Double_t rMax = station->GetRmax();
406 if (IntersectsRad(sector, module->GetCutRadius()) == 2 || !IntersectsRad(sector, rMax)) return;
407 nCols = isIncomplete ? nCols : fNCols[iStation].at(iRegion);
408 nRows = isIncomplete ? nRows : fNRows[iStation].at(iRegion);
409 module->AddSector(new CbmMuchSectorRectangular(detectorId, iSector, secPosition, secSize, nCols, nRows));
410 return;
411 }
412
413 // Needed for the case of incomplete sectors
414 Int_t iReg = -1;
415 Int_t nC = Int_t(padLx / GetPadMaxSize(module, "Width"));
416 Int_t nR = Int_t(padLy / GetPadMaxSize(module, "Length"));
417 Double_t pLx = nC == 0 ? padLx : nC * GetPadMaxSize(module, "Width");
418 Double_t pLy = nR == 0 ? padLy : nR * GetPadMaxSize(module, "Length");
419 Double_t sLx = nC == 0 ? secLx : nC * GetSectorMaxSize(module, "Width", iReg);
420 Double_t sLy = nR == 0 ? secLy : nR * GetSectorMaxSize(module, "Length", iReg);
421 nCols = Int_t(sLx / pLx);
422 nRows = Int_t(sLy / pLy);
423 delete sector;
424
425 TVector3 position, size;
426 Double_t newSecLx = 0., newSecLy = 0., newSecX = 0., newSecY = 0.;
427 for (Int_t i = 0; i < 2; ++i) {
428 newSecLx = resultX ? i * (secLx - newSecLx) + (1 - i) * pLx / 2. * nCols : secLx;
429 newSecLy = resultY ? i * (secLy - newSecLy) + (1 - i) * pLy / 2. * nRows : secLy;
430 newSecX = resultX ? secX - TMath::Power(-1, i) * nX * (secLx / 2. - newSecLx / 2.) : secX;
431 newSecY = resultY ? secY - TMath::Power(-1, i) * nY * (secLy / 2. - newSecLy / 2.) : secY;
432
433 Double_t newPadLx = resultX ? pLx / 2. : padLx;
434 Double_t newPadLy = resultY ? pLy / 2. : padLy;
435 nCols = Int_t(newSecLx / newPadLx);
436 nRows = Int_t(newSecLy / newPadLy);
437
438 position.SetXYZ(newSecX, newSecY, secZ);
439 size.SetXYZ(newSecLx, newSecLy, secLz);
440 iSector = module->GetNSectors();
441 SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, position, size, nCols, nRows));
442 }
443}
444// -------------------------------------------------------------------------
445
446// ----- Private method IntersectsRad -----------------------------------
448{
449 if (radius < 0) return 0;
450
451 Int_t intersects = 0;
452
453 Double_t modLx = module->GetSize()[0];
454 Double_t modLy = module->GetSize()[1];
455 Double_t modX = module->GetPosition()[0];
456 Double_t modY = module->GetPosition()[1];
457
458 Double_t ulR = TMath::Sqrt((modX - modLx / 2.) * (modX - modLx / 2.) + (modY + modLy / 2.) * (modY + modLy / 2.));
459 Double_t urR = TMath::Sqrt((modX + modLx / 2.) * (modX + modLx / 2.) + (modY + modLy / 2.) * (modY + modLy / 2.));
460 Double_t blR = TMath::Sqrt((modX - modLx / 2.) * (modX - modLx / 2.) + (modY - modLy / 2.) * (modY - modLy / 2.));
461 Double_t brR = TMath::Sqrt((modX + modLx / 2.) * (modX + modLx / 2.) + (modY - modLy / 2.) * (modY - modLy / 2.));
462
463 if (ulR < radius) intersects++;
464 if (urR < radius) intersects++;
465 if (blR < radius) intersects++;
466 if (brR < radius) intersects++;
467
468 if (intersects == 4) return 2;
469 if (intersects) return 1;
470 else
471 return 0;
472}
473// -------------------------------------------------------------------------
474
475// ----- Private method IntersectsRad -----------------------------------
477{
478 if (radius < 0) return 0;
479
480 Int_t intersects = 0;
481
482 Double_t secLx = sector->GetSize()[0];
483 Double_t secLy = sector->GetSize()[1];
484 Double_t secX = sector->GetPosition()[0];
485 Double_t secY = sector->GetPosition()[1];
486
487 Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
488 Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
489 Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
490 Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
491
492 if (ulR < radius) intersects++;
493 if (urR < radius) intersects++;
494 if (blR < radius) intersects++;
495 if (brR < radius) intersects++;
496
497 if (intersects == 4) return 2;
498 if (intersects) return 1;
499 else
500 return 0;
501}
502// -------------------------------------------------------------------------
503
504// ----- Private method ShouldSegment -----------------------------------
506{
507 Double_t secLx = sector->GetSize()[0];
508 Double_t secLy = sector->GetSize()[1];
509 Double_t secArea = secLx * secLy;
510 Double_t secL = direction == "X" ? secLx : secLy;
511 // Bool_t isIncomplete = IsIncompleteSector(sector);
512
514
515 // Get region index for the sector
516 iRegion = GetRegionIndex(sector);
517 Double_t secRegL = direction == "X" ? fSecLx[iStation].at(iRegion) : fSecLy[iStation].at(iRegion);
518 Double_t secRegArea = fSecLx[iStation].at(iRegion) * fSecLy[iStation].at(iRegion);
519
520 if (secArea > secRegArea) {
521 // If sector length is larger than it's width
522 if (secLy > secLx && direction == "X" && ShouldSegment(sector, "Y", iRegion)) return false;
523 // If sector width is larger than or equal to it's length
524 if (secLy <= secLx && direction == "Y") return false;
525
526 // If sector size is larger than that corresponding to the region
527 if (secL > secRegL) return true;
528 }
529 return false;
530}
531// -------------------------------------------------------------------------
532
533// ----- Private method GetRegionIndex ----------------------------------
535{
537 Double_t secLx = sector->GetSize()[0];
538 Double_t secLy = sector->GetSize()[1];
539 Double_t secArea = secLx * secLy;
540 Double_t sX = TMath::Abs(sector->GetPosition()[0]) - secLx / 2.;
541 Double_t sY = TMath::Abs(sector->GetPosition()[1]) - secLy / 2.;
542 Double_t secRad = TMath::Sqrt(sX * sX + sY * sY);
543
544 Int_t iRegion = fNRegions[iStation] - 1;
545
546 for (Int_t iReg = 0; iReg < fNRegions[iStation]; ++iReg) {
547 Double_t secRegArea = fSecLx[iStation].at(iReg) * fSecLy[iStation].at(iReg);
548 Double_t regionRad = fRadii[iStation].at(iReg);
549 if (secRad > regionRad) continue;
550
551 iRegion = iReg;
552 if (iReg > 0 && !IsIncompleteSector(sector)) {
553 Double_t secPrevRegArea = fSecLx[iStation].at(iReg - 1) * fSecLy[iStation].at(iReg - 1);
554 if (secArea < secRegArea && secArea >= secPrevRegArea) iRegion--;
555 }
556 break;
557 }
558 return iRegion;
559}
560// -------------------------------------------------------------------------
561
562// ----- Private method GetSectorMaxSize --------------------------------
563Double_t CbmMuchSegmentManual::GetSectorMaxSize(CbmMuchModuleGemRectangular* module, const TString side, Int_t& iRegion)
564{
566 Int_t nRegions = fNRegions[iStation];
567 for (iRegion = 0; iRegion < nRegions; ++iRegion) {
568 Double_t rad = fRadii[iStation].at(iRegion);
569 Int_t result = IntersectsRad(module, rad);
570 if (result == 2) return side == "Width" ? fSecLx[iStation].at(iRegion) : fSecLy[iStation].at(iRegion);
571 }
572 iRegion = nRegions - 1;
573 return side == "Width" ? fSecLx[iStation].at(iRegion) : fSecLy[iStation].at(iRegion);
574}
575// -------------------------------------------------------------------------
576
577// ----- Private method GetPadMaxSize -----------------------------------
579{
581 Int_t iRegion = -1;
582 Double_t sectorSize = GetSectorMaxSize(module, side, iRegion);
583 return side == "Width" ? sectorSize / fNCols[iStation].at(iRegion) : sectorSize / fNRows[iStation].at(iRegion);
584}
585// -------------------------------------------------------------------------
586
587// ----- Private method IsIncompleteSector ------------------------------
589{
590 Bool_t result = false;
592 Double_t secLx = sector->GetSize()[0];
593 Double_t secLy = sector->GetSize()[1];
594 Double_t minL = TMath::Min(secLx, secLy);
595 Double_t maxL = TMath::Max(secLx, secLy);
596 Int_t nFrac = Int_t((maxL + 1e-5) / minL);
597 Int_t nPower = Int_t(TMath::Log2(nFrac) + 1e-2);
598 Double_t maxL1 = minL * TMath::Power(2, nPower);
599
600 if (TMath::Abs(maxL - maxL1) > 1e-5 || sector->GetNChannels() < fNChannels[iStation]) result = true;
601 return result;
602}
603// -------------------------------------------------------------------------
604
605// ----- Private method Print -------------------------------------------
606void CbmMuchSegmentManual::Print(Option_t*) const
607{
608 printf("Segmentation written to file %s\n", fDigiFileName);
609 Int_t nTotSectors = 0;
610 Int_t nTotChannels = 0;
611 Int_t nTotGems = 0;
612 Int_t nTotStraws = 0;
613 printf("====================================================================="
614 "============================\n");
615 for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
616 CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
617 Int_t nGems = 0;
618 Int_t nStraws = 0;
619 Int_t nSectors = 0;
620 Int_t nChannels = 0;
621 Double_t padMaxLx = std::numeric_limits<Double_t>::min();
622 Double_t padMinLx = std::numeric_limits<Double_t>::max();
623 Double_t padMaxLy = std::numeric_limits<Double_t>::min();
624 Double_t padMinLy = std::numeric_limits<Double_t>::max();
625 Double_t secMaxLx = std::numeric_limits<Double_t>::min();
626 Double_t secMinLx = std::numeric_limits<Double_t>::max();
627 Double_t secMaxLy = std::numeric_limits<Double_t>::min();
628 Double_t secMinLy = std::numeric_limits<Double_t>::max();
629 if (!station) continue;
630 for (Int_t iLayer = 0; iLayer < station->GetNLayers(); ++iLayer) {
631 CbmMuchLayer* layer = station->GetLayer(iLayer);
632 if (!layer) continue;
633 for (Int_t iSide = 0; iSide < 2; ++iSide) {
634 CbmMuchLayerSide* side = layer->GetSide(iSide);
635 if (!side) continue;
636 for (Int_t iModule = 0; iModule < side->GetNModules(); ++iModule) {
637 CbmMuchModule* mod = side->GetModule(iModule);
638 if (!mod) continue;
639 switch (mod->GetDetectorType()) {
640 case 1: { // GEMs
642 if (!module) continue;
643 nGems++;
644 nSectors += module->GetNSectors();
645 for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
647 if (!sector) continue;
648 Double_t padLx = sector->GetPadDx();
649 Double_t padLy = sector->GetPadDy();
650 if (padLx > padMaxLx) {
651 padMaxLx = padLx;
652 secMaxLx = sector->GetPadNx() * padLx;
653 }
654 if (padLx < padMinLx) {
655 padMinLx = padLx;
656 secMinLx = sector->GetPadNx() * padLx;
657 }
658 if (padLy > padMaxLy) {
659 padMaxLy = padLy;
660 secMaxLy = sector->GetPadNy() * padLy;
661 }
662 if (padLy < padMinLy) {
663 padMinLy = padLy;
664 secMinLy = sector->GetPadNy() * padLy;
665 }
666 nChannels += sector->GetNChannels();
667 }
668 break;
669 }
670 case 2: // Straw tubes
671 nStraws++;
672 break;
673 }
674 }
675 }
676 }
677 printf("Station %i:\n", iStation + 1);
678 printf(" GEM modules: %i\n", nGems);
679 if (nGems) {
680 printf(" Sectors: %i, Min.Sector size:%3.2fx%3.2f, Max.Sector "
681 "size:%3.2fx%3.2f\n",
682 nSectors, secMinLx, secMinLy, secMaxLx, secMaxLy);
683 for (Int_t iReg = 0; iReg < fNRegions.at(iStation); ++iReg) {
684 printf("Region %i: size %fx%f\n", iReg, fSecLx.at(iStation).at(iReg), fSecLy.at(iStation).at(iReg));
685 }
686 printf(" Pads: %i, Min.Pad size:%3.2fx%3.2f, Max.Pad size:%3.2fx%3.2f\n", nChannels, padMinLx, padMinLy,
687 padMaxLx, padMaxLy);
688 }
689 printf(" Straw modules: %i\n", nStraws);
690 nTotSectors += nSectors;
691 nTotChannels += nChannels;
692 nTotGems += nGems;
693 nTotStraws += nStraws;
694 }
695 printf("---------------------------------------------------------------------"
696 "----------------------------\n");
697 printf(" Summary: \n GEM modules: %i\n Sectors: %i, Pads: %i\n "
698 "Straw modules: %i\n",
699 nTotGems, nTotSectors, nTotChannels, nTotStraws);
700 printf("====================================================================="
701 "============================\n");
702}
703
705{
706 // Change file extension
707 char txtfile[100];
708 Int_t length = strlen(fDigiFileName);
709 Int_t iChar;
710 for (iChar = length - 1; iChar >= 0; --iChar) {
711 if (fDigiFileName[iChar] == '.') break;
712 }
713 strncpy(txtfile, fDigiFileName, iChar + 1);
714 txtfile[iChar + 1] = '\0';
715 strcat(txtfile, "txt");
716
717 FILE* outfile;
718 outfile = fopen(txtfile, "w");
719 /*
720 Int_t colors[] = {kGreen, kMagenta, kCyan, kRed, kBlue, kYellow, kTeal,
721 kPink, kAzure, kOrange, kViolet, kSpring,
722 kGreen+2, kMagenta+2, kCyan+2, kRed+2, kBlue+2, kYellow+2, kTeal+2,
723 kPink+2, kAzure+2, kOrange+2, kViolet+2, kSpring+2,
724 kGreen+4, kMagenta+4, kCyan+4, kRed+4, kBlue+4, kYellow+4, kTeal+4,
725 kPink+4, kAzure+4, kOrange+4, kViolet+4, kSpring+4};
726*/
727 for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
728 fprintf(outfile, "=================================================================="
729 "=========\n");
730 fprintf(outfile, "Station %i\n", iStation + 1);
731 fprintf(outfile, "Sector size, cm Sector position, cm Number of pads Side "
732 "Pad size, cm\n");
733 fprintf(outfile, "------------------------------------------------------------------"
734 "----------\n");
735 TCanvas* c1 = new TCanvas(Form("station%i", iStation + 1), Form("station%i", iStation + 1), 800, 800);
736 c1->SetFillColor(0);
737 c1->Range(-250, -250, 250, 250);
738 CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
739 CbmMuchLayer* layer = station->GetLayer(0);
740 for (Int_t iSide = 1; iSide >= 0; iSide--) {
741 CbmMuchLayerSide* layerSide = layer->GetSide(iSide);
742 for (Int_t iModule = 0; iModule < layerSide->GetNModules(); ++iModule) {
743 CbmMuchModule* mod = layerSide->GetModule(iModule);
744 if (mod->GetDetectorType() != 1) continue;
746 module->InitModule();
747 for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
749 // Reject incomplete sectors by size
750 Double_t secLx = sector->GetSize()[0];
751 Double_t secLy = sector->GetSize()[1];
752 // Bool_t isIncomplete = IsIncompleteSector(sector);
753 // TODO
754 // if(isIncomplete){
755 // iSide ? sector->SetFillColor(TColor::GetColorDark(kGray+1)) : sector->SetFillColor(kGray + 1);
756 // }
757 // else{
758 // Int_t i = Int_t((secLx+1e-5)/fSecLx[iStation].at(0)) - 1;
759 // Int_t j = Int_t((secLy+1e-5)/fSecLy[iStation].at(0)) - 1;
760 // sector->SetFillColor(iSide ? TColor::GetColorDark(colors[i+j]) : colors[i+j]);
761 // }
762 sector->DrawPads();
763 const char* side = iSide ? "Back" : "Front";
764 fprintf(outfile, "%-4.2fx%-10.2f %-6.2fx%-12.2f %-14i %-5s ", secLx, secLy, sector->GetPosition()[0],
765 sector->GetPosition()[1], sector->GetNChannels(), side);
766 fprintf(outfile, "%-4.2fx%-4.2f\n", sector->GetPadDx(), sector->GetPadDy());
767 } // sectors
768 } // modules
769 } // sides
770
771 // Draw a hole
772 TArc* holeArc = new TArc(0., 0., station->GetRmin());
773 holeArc->Draw();
774
775 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
776 TArc* arc = new TArc(0., 0., fRadii[iStation].at(iRegion));
777 arc->SetLineColor(kBlack);
778 arc->SetLineWidth(2);
779 arc->SetFillStyle(0);
780 arc->Draw();
781 }
782
783 c1->Print(Form("%s/station%i.eps", gSystem->DirName(fDigiFileName), iStation + 1));
784 c1->Print(Form("%s/station%i.png", gSystem->DirName(fDigiFileName), iStation + 1));
785 } //stations
786 fclose(outfile);
787}
788
790{
791 ifstream infile;
792 infile.open(fInputFileName);
793 if (!infile) { Fatal("ReadInputFile", "Error: Cannot open the input file."); }
794
795 Int_t index;
796 string str;
797 vector<string> strs;
798
799 // Set number of stations
800 OmitDummyLines(infile, str);
801 strs = Split(str, ' ');
802 index = strs.size() - 1;
803 Int_t nStations;
804 StrToNum(strs.at(index), nStations);
805 SetNStations(nStations);
806 printf("Number of stations: \t%i", fNStations);
807
808 // Set number of channels
809 OmitDummyLines(infile, str);
810 strs = Split(str, ' ');
811 printf("\nNumber of channels: ");
812 for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
813 index = strs.size() - fNStations + iStation;
814 Int_t nChannels;
815 StrToNum(strs.at(index), nChannels);
816 SetNChannels(iStation, nChannels);
817 printf("\t%i", fNChannels[iStation]);
818 }
819
820 // Set number of regions
821 OmitDummyLines(infile, str);
822 strs = Split(str, ' ');
823 printf("\nNumber of regions: ");
824 for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
825 index = strs.size() - fNStations + iStation;
826 Int_t nRegions;
827 StrToNum(strs.at(index), nRegions);
828 SetNRegions(iStation, nRegions);
829 printf("\t%i", fNRegions[iStation]);
830 }
831
832 for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
833 vector<Double_t> padWidth(fNRegions[iStation], 0);
834 vector<Double_t> padLength(fNRegions[iStation], 0);
835
836 // Set region radii
837 OmitDummyLines(infile, str);
838 strs = Split(str, ' ');
839 printf("\nStation %i", iStation + 1);
840 printf("\n Region radii [cm]: ");
841 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
842 index = strs.size() - fNRegions[iStation] + iRegion;
843 Double_t radius;
844 StrToNum(strs.at(index), radius);
845 SetRegionRadius(iStation, iRegion, radius);
846 printf("\t%4.2f", fRadii[iStation].at(iRegion));
847 }
848
849 // Set pad width
850 OmitDummyLines(infile, str);
851 strs = Split(str, ' ');
852 printf("\n Pad width [cm]: ");
853 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
854 index = strs.size() - fNRegions[iStation] + iRegion;
855 StrToNum(strs.at(index), padWidth.at(iRegion));
856 printf("\t%4.2f", padWidth.at(iRegion));
857 }
858
859 // Set pad length
860 OmitDummyLines(infile, str);
861 strs = Split(str, ' ');
862 printf("\n Pad length [cm]: ");
863 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion) {
864 index = strs.size() - fNRegions[iStation] + iRegion;
865 StrToNum(strs.at(index), padLength.at(iRegion));
866 printf("\t%4.2f", padLength.at(iRegion));
867 }
868
869 for (Int_t iRegion = 0; iRegion < fNRegions[iStation]; ++iRegion)
870 SetPadSize(iStation, iRegion, padWidth.at(iRegion), padLength.at(iRegion));
871 }
872 infile.close();
873}
874
ClassImp(CbmConverterManager)
static constexpr size_t size()
Definition KfSimdPseudo.h:2
int Int_t
bool Bool_t
static int32_t GetStationIndex(int32_t address)
Int_t GetNModules() const
CbmMuchModule * GetModule(Int_t iModule) const
CbmMuchLayerSide * GetSide(Bool_t side)
CbmMuchLayerSide * GetSideF()
CbmMuchLayerSide * GetSideB()
CbmMuchSector * GetSectorByIndex(Int_t iSector)
Int_t GetNSectors() const
Int_t GetDetectorId() const
Double_t GetCutRadius() const
Int_t GetDetectorType() const
UInt_t GetAddress() const
Int_t GetNChannels() const
void SetNRegions(Int_t iStation, Int_t nRegions)
std::vector< std::string > & Split(const std::string &s, char delim, std::vector< std::string > &elems)
void SetRegionRadius(Int_t iStation, Int_t iRegion, Double_t radius)
void Print(Option_t *="") const
Bool_t ShouldSegment(CbmMuchSectorRectangular *sector, const TString direction, Int_t &iRegion)
std::map< Int_t, std::vector< Double_t > > fSecLy
void SetNStations(Int_t nStations)
Bool_t IsIncompleteSector(CbmMuchSectorRectangular *sector)
void StrToNum(std::string &str, T &number)
std::map< Int_t, std::vector< Int_t > > fNCols
std::map< Int_t, Int_t > fNChannels
void SetSigma(Int_t iStation, Int_t iRegion, Double_t sigmaX, Double_t sigmaY)
Int_t IntersectsRad(CbmMuchModuleGemRectangular *module, Double_t radius)
void SegmentSector(CbmMuchModuleGemRectangular *module, CbmMuchSectorRectangular *sector)
void SetNChannels(Int_t iStation, Int_t nChannels)
std::map< Int_t, std::vector< Int_t > > fNRows
std::map< Int_t, std::vector< Double_t > > fRadii
void SetPadSize(Int_t iStation, Int_t iRegion, Double_t padLx, Double_t padLy)
void OmitDummyLines(std::ifstream &infile, std::string &str)
void SegmentLayerSide(CbmMuchLayerSide *layerSide)
void SegmentModule(CbmMuchModuleGemRectangular *module, Bool_t useModuleDesign)
std::map< Int_t, std::vector< Double_t > > fSecLx
Int_t GetRegionIndex(CbmMuchSectorRectangular *sector)
Double_t GetSectorMaxSize(CbmMuchModuleGemRectangular *module, const TString side, Int_t &iRegion)
std::map< Int_t, Int_t > fNRegions
Double_t GetPadMaxSize(CbmMuchModuleGemRectangular *module, const TString side)
virtual InitStatus Init()
Double_t GetRmin() const
CbmMuchLayer * GetLayer(Int_t iLayer) const
Double_t GetRmax() const
Int_t GetNLayers() const