CbmRoot
CbmMuchSegmentAuto.cxx
Go to the documentation of this file.
1 
9 #include "CbmMuchSegmentAuto.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 "CbmMuchModuleGem.h" // for CbmMuchModuleGem
17 #include "CbmMuchPoint.h" // for CbmMuchPoint
18 #include "CbmMuchSectorRectangular.h" // for CbmMuchSectorRectangular
19 #include "CbmMuchStation.h" // for CbmMuchStation
20 
21 #include <FairRootManager.h> // for FairRootManager
22 #include <FairRuntimeDb.h> // for FairRuntimeDb
23 #include <FairTask.h> // for FairTask, InitStatus, kSUCCESS
24 
25 #include <TArc.h> // for TArc
26 #include <TAxis.h> // for TAxis
27 #include <TCanvas.h> // for TCanvas
28 #include <TClonesArray.h> // for TClonesArray
29 #include <TF1.h> // for TF1
30 #include <TFile.h> // for TFile, gFile
31 #include <TH1.h> // for TH1D
32 #include <TMath.h> // for Sqrt, Log2, Pi
33 #include <TMathBase.h> // for Abs
34 #include <TObjArray.h> // for TObjArray
35 #include <TStyle.h> // for TStyle, gStyle
36 #include <TSystem.h> // for TSystem, gSystem
37 #include <TVector3.h> // for TVector3
38 
39 #include <cassert> // for assert
40 #include <iosfwd> // for string
41 #include <limits> // for numeric_limits
42 #include <math.h> // for exp
43 #include <stdio.h> // for printf, fprintf, fclose
44 
45 using std::string;
46 
47 // ----- Default constructor -------------------------------------------
48 CbmMuchSegmentAuto::CbmMuchSegmentAuto()
49  : FairTask()
50  , fEvents(0)
51  , fPoints(nullptr)
52  , fHistHitDensity()
53  , fNStations(0)
54  , fStations(nullptr)
55  , fDigiFileName()
56  , fGeoPar(nullptr)
57  , fExp0()
58  , fExp1()
59  , fSigmaXmin()
60  , fSigmaYmin()
61  , fSigmaXmax()
62  , fSigmaYmax()
63  , fOccupancyMax() {}
64 // -------------------------------------------------------------------------
65 
66 // ----- Standard constructor ------------------------------------------
67 CbmMuchSegmentAuto::CbmMuchSegmentAuto(const char* digiFileName)
68  : FairTask()
69  , fEvents(0)
70  , fPoints(nullptr)
71  , fHistHitDensity()
72  , fNStations(0)
73  , fStations(nullptr)
74  , fDigiFileName(digiFileName)
75  , fGeoPar(nullptr)
76  , fExp0()
77  , fExp1()
78  , fSigmaXmin()
79  , fSigmaYmin()
80  , fSigmaXmax()
81  , fSigmaYmax()
82  , fOccupancyMax() {}
83 // -------------------------------------------------------------------------
84 
85 // ----- Destructor ----------------------------------------------------
86 CbmMuchSegmentAuto::~CbmMuchSegmentAuto() {}
87 // -------------------------------------------------------------------------
88 
89 void CbmMuchSegmentAuto::SetNStations(Int_t nStations) {
90  fNStations = nStations;
91  fSigmaXmin.resize(fNStations);
92  fSigmaYmin.resize(fNStations);
93  fSigmaXmax.resize(fNStations);
94  fSigmaYmax.resize(fNStations);
95  fOccupancyMax.resize(fNStations);
96 }
97 
98 void CbmMuchSegmentAuto::SetSigmaMin(Double_t* sigmaXmin, Double_t* sigmaYmin) {
99  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
100  fSigmaXmin.at(iStation) = sigmaXmin[iStation];
101  fSigmaYmin.at(iStation) = sigmaYmin[iStation];
102  }
103 }
104 void CbmMuchSegmentAuto::SetSigmaMax(Double_t* sigmaXmax, Double_t* sigmaYmax) {
105  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
106  fSigmaXmax.at(iStation) = sigmaXmax[iStation];
107  fSigmaYmax.at(iStation) = sigmaYmax[iStation];
108  }
109 }
110 
111 void CbmMuchSegmentAuto::SetOccupancyMax(Double_t* occupancyMax) {
112  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
113  fOccupancyMax.at(iStation) = occupancyMax[iStation];
114  }
115 }
116 
117 // ----- Private method SetParContainers --------------------------------
118 void CbmMuchSegmentAuto::SetParContainers() {
119  // Get runtime database
120  FairRuntimeDb* db = FairRuntimeDb::instance();
121  if (!db) Fatal("SetParContainers", "No runtime database");
122  // Get MUCH geometry param
123  fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
124 }
125 // -------------------------------------------------------------------------
126 
127 
128 // ----- Private method Init --------------------------------------------
129 InitStatus CbmMuchSegmentAuto::Init() {
130  FairRootManager* fManager = FairRootManager::Instance();
131  fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
132  fEvents = 0;
133 
134  fStations = fGeoPar->GetStations();
135  if (!fStations) Fatal("Init", "No input array of MUCH stations.");
136  if (fNStations != fStations->GetEntries())
137  Fatal("Init", "Incorrect number of stations set.");
138  //fNStations = fStations->GetEntries();
139  fHistHitDensity = new TH1D*[fNStations];
140 
141  for (Int_t i = 0; i < fNStations; i++) {
142  // CbmMuchStation* station = (CbmMuchStation*) fStations->At(i);
143  fHistHitDensity[i] = new TH1D(
144  Form("hStation%i", i + 1), Form("Station %i", i + 1), 110, 0, 220);
145  fHistHitDensity[i]->GetXaxis()->SetTitle("r, cm");
146  fHistHitDensity[i]->GetYaxis()->SetTitle("hits/(event#timescm^{2})");
147  fHistHitDensity[i]->SetLineColor(4);
148  fHistHitDensity[i]->SetLineWidth(2);
149  fHistHitDensity[i]->GetYaxis()->SetTitleOffset(1.27);
150  }
151  return kSUCCESS;
152 }
153 // -------------------------------------------------------------------------
154 
155 // -------------------------------------------------------------------------
156 void CbmMuchSegmentAuto::Exec(Option_t*) {
157  fEvents++;
158  printf("Event: %i\n", fEvents);
159 
160  gStyle->SetOptStat(0);
161  for (Int_t iPoint = 0; iPoint < fPoints->GetEntriesFast(); iPoint++) {
162  CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(iPoint);
163  if (!point) continue;
164  Int_t iStation = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
165  // CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
166  // if (station->GetDetectorType()==2) continue;
167 
168  Int_t iLayer = CbmMuchAddress::GetLayerIndex(point->GetDetectorID());
169  // printf("iStation = %i\n", iStation);
170  // printf("detId = %qd\n", point->GetDetectorID());
171  // printf("iPoint = %i\n", iPoint);
172  assert(iStation >= 0 && iStation < fNStations);
173  if (iLayer) continue;
174  TVector3 pos;
175  point->Position(pos);
176  ((TH1D*) fHistHitDensity[iStation])->Fill(pos.Pt());
177  }
178 }
179 // -------------------------------------------------------------------------
180 
181 // -------------------------------------------------------------------------
182 void CbmMuchSegmentAuto::FinishTask() {
183  // Create normalization histo
184  TH1D* hNorm = new TH1D("hNorm", "", 110, 0, 220);
185  Double_t binSize = 2.;
186  for (Int_t l = 0; l < 100; l++) {
187  Double_t R1 = l * binSize;
188  Double_t R2 = (l + 1) * binSize;
189  Double_t s = TMath::Pi() * (R2 * R2 - R1 * R1);
190  hNorm->SetBinContent(l + 1, s * fEvents);
191  }
192 
193  for (Int_t i = 0; i < fNStations; i++) {
194  CbmMuchStation* station = (CbmMuchStation*) fStations->At(i);
195  TH1D* h = fHistHitDensity[i];
196  h->Divide(hNorm);
197  TCanvas* c1 = new TCanvas(
198  Form("cStation%i", i + 1), Form("Station %i", i + 1), 10, 10, 500, 500);
199  c1->SetLogy();
200  c1->SetLeftMargin(0.12);
201  c1->SetRightMargin(0.08);
202  TF1* f1 = new TF1("f1", "exp([0] + [1]*x)");
203  f1->SetParameter(0, 0.5);
204  f1->SetParameter(1, -0.1);
205  h->Fit("f1", "QLL", "", station->GetRmin(), station->GetRmax());
206  Double_t exp0 = h->GetFunction("f1")->GetParameter(0);
207  Double_t exp1 = h->GetFunction("f1")->GetParameter(1);
208  fExp0.push_back(exp0);
209  fExp1.push_back(exp1);
210 
211  h->Draw();
212  c1->Print(
213  Form("%s/hd_station%i.eps", gSystem->DirName(fDigiFileName), i + 1));
214  c1->Print(
215  Form("%s/hd_station%i.png", gSystem->DirName(fDigiFileName), i + 1));
216  }
217 
218  for (Int_t i = 0; i < fNStations; i++) {
219  CbmMuchStation* station = (CbmMuchStation*) fStations->At(i);
220  //if (station->GetDetectorType()==2) continue;
221 
222  Int_t nLayers = station->GetNLayers();
223  for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
224  CbmMuchLayer* layer = station->GetLayer(iLayer);
225  if (!layer) Fatal("Init", "Incomplete layers array.");
226  // Get layer sides
227  InitLayerSide(layer->GetSideF());
228  InitLayerSide(layer->GetSideB());
229  }
230  printf("Station%i segmented\n", i + 1);
231  }
232 
233  // Save parameters
234  TFile* oldfile = gFile;
235  TFile* f = new TFile(fDigiFileName, "RECREATE");
236  fStations->Write("stations", 1);
237 
238  f->Close();
239  gFile = oldfile;
240 
241  DrawSegmentation();
242 
243  Print();
244 }
245 // -------------------------------------------------------------------------
246 
247 
248 // ----- Private method InitLayerSide -----------------------------------
249 void CbmMuchSegmentAuto::InitLayerSide(CbmMuchLayerSide* layerSide) {
250  if (!layerSide) Fatal("Init", "Incomplete layer sides array.");
251  Int_t nModules = layerSide->GetNModules();
252  for (Int_t iModule = 0; iModule < nModules; iModule++) {
253  CbmMuchModule* mod = layerSide->GetModule(iModule);
254  if (mod->GetDetectorType() != 1) continue;
255  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
256  SegmentModule(module);
257  }
258 }
259 // -------------------------------------------------------------------------
260 
261 
262 // ----- Private method SegmentModule -----------------------------------
263 void CbmMuchSegmentAuto::SegmentModule(CbmMuchModuleGem* module) {
264  TVector3 modSize = module->GetSize();
265  Double_t modLx = modSize.X();
266  Double_t modLy = modSize.Y();
267  Double_t modLz = modSize.Z();
268  TVector3 modPosition = module->GetPosition();
269  Double_t modX = modPosition.X();
270  Double_t modY = modPosition.Y();
271  Double_t modZ = modPosition.Z();
272 
273  Bool_t result = modLx > modLy;
274  Int_t iRatio = (result) ? (Int_t)((modLx + 1e-3) / modLy)
275  : (Int_t)((modLy + 1e-3) / modLx);
276  Int_t detectorId = module->GetDetectorId();
277  Double_t secLx = (result) ? modLx / iRatio : modLx;
278  Double_t secLy = (result) ? modLy : modLy / iRatio;
279  for (Int_t i = 0; i < iRatio; i++) {
280  Double_t secX = (result) ? modX - modLx / 2. + (i + 0.5) * secLx : modX;
281  Double_t secY = (result) ? modY : modY - modLy / 2. + (i + 0.5) * secLy;
282  Int_t iSector = module->GetNSectors();
283 
284  TVector3 sectorPosition(secX, secY, modZ);
285  TVector3 sectorSize(secLx, secLy, modLz);
286  CbmMuchSectorRectangular* sector = new CbmMuchSectorRectangular(
287  detectorId, iSector, sectorPosition, sectorSize, 8, 16);
288  SegmentSector(module, sector);
289  }
290 }
291 // -------------------------------------------------------------------------
292 
293 // ----- Private method SegmentSector ------------------------------------
294 void CbmMuchSegmentAuto::SegmentSector(CbmMuchModuleGem* module,
295  CbmMuchSectorRectangular* sector) {
296  TVector3 secSize = sector->GetSize();
297  TVector3 secPosition = sector->GetPosition();
298  Int_t detectorId = module->GetDetectorId();
299  Int_t iStation = CbmMuchAddress::GetStationIndex(detectorId);
300  Int_t iLayer = CbmMuchAddress::GetLayerIndex(detectorId);
301  Int_t iSide = CbmMuchAddress::GetLayerSideIndex(detectorId);
302  Int_t iModule = CbmMuchAddress::GetModuleIndex(detectorId);
303  Double_t secLx = secSize.X();
304  Double_t secLy = secSize.Y();
305  Double_t secLz = secSize.Z();
306  Double_t secX = secPosition.X();
307  Double_t secY = secPosition.Y();
308  Double_t secZ = secPosition.Z();
309 
310  assert(CbmMuchAddress::GetStationIndex(sector->GetAddress()) == iStation);
311  assert(CbmMuchAddress::GetLayerIndex(sector->GetAddress()) == iLayer);
312  assert(CbmMuchAddress::GetLayerSideIndex(sector->GetAddress()) == iSide);
313  assert(CbmMuchAddress::GetModuleIndex(sector->GetAddress()) == iModule);
314 
315  Bool_t result1 = ShouldSegmentByX(sector);
316  Bool_t result2 = ShouldSegmentByY(sector);
317 
318  if (result1 || result2) {
319  delete sector;
320  } else {
321  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
322  Double_t rMax = station->GetRmax();
323  if ((IntersectsRad(sector, module->GetCutRadius()) == 2)
324  || !IntersectsRad(sector, rMax)) {
325  delete sector;
326  return;
327  }
328  module->AddSector(sector);
329  return;
330  }
331 
332  Int_t iSector;
333  Double_t newSecX, newSecY, newSecLx, newSecLy;
334  Bool_t equal = TMath::Abs(secLx - secLy) < 1e-5;
335  Bool_t res = secLx > secLy;
336  TVector3 position, size;
337  if (result1 && result2) {
338  if (equal) { res = kTRUE; }
339  for (Int_t i = 0; i < 2; ++i) {
340  iSector =
341  module
342  ->GetNSectors(); //CbmMuchGeoScheme::GetDetIdFromModule(moduleId, module->GetNSectors());
343  newSecLx = res ? secLx / 2. : secLx;
344  newSecLy = res ? secLy : secLy / 2.;
345  newSecX = res ? secX + (i - 0.5) * newSecLx : secX;
346  newSecY = res ? secY : secY + (i - 0.5) * newSecLy;
347  position.SetXYZ(newSecX, newSecY, secZ);
348  size.SetXYZ(newSecLx, newSecLy, secLz);
349  SegmentSector(module,
350  new CbmMuchSectorRectangular(
351  detectorId, iSector, position, size, 8, 16));
352  }
353  } else if (result1 || result2) {
354  for (Int_t i = 0; i < 2; i++) {
355  iSector =
356  module
357  ->GetNSectors(); //CbmMuchGeoScheme::GetDetIdFromModule(moduleId, module->GetNSectors());
358  newSecLx = result1 ? secLx / 2. : secLx;
359  newSecLy = result1 ? secLy : secLy / 2;
360  newSecX = result1 ? secX + (i - 0.5) * newSecLx : secX;
361  newSecY = result1 ? secY : secY + (i - 0.5) * newSecLy;
362  position.SetXYZ(newSecX, newSecY, secZ);
363  size.SetXYZ(newSecLx, newSecLy, secLz);
364  SegmentSector(module,
365  new CbmMuchSectorRectangular(
366  detectorId, iSector, position, size, 8, 16));
367  }
368  }
369 }
370 // -------------------------------------------------------------------------
371 
372 // ----- Private method ShouldSegmentByX ---------------------------------
373 Bool_t CbmMuchSegmentAuto::ShouldSegmentByX(CbmMuchSectorRectangular* sector) {
374  Double_t secLx = sector->GetSize()[0];
375  Double_t secLy = sector->GetSize()[1];
376  Double_t secX = sector->GetPosition()[0];
377  Double_t secY = sector->GetPosition()[1];
378 
379  Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.)
380  + (secY + secLy / 2.) * (secY + secLy / 2.));
381  Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.)
382  + (secY + secLy / 2.) * (secY + secLy / 2.));
383  Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.)
384  + (secY - secLy / 2.) * (secY - secLy / 2.));
385  Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.)
386  + (secY - secLy / 2.) * (secY - secLy / 2.));
387 
388  Double_t uR = (ulR < urR) ? ulR : urR;
389  Double_t bR = (blR < brR) ? blR : brR;
390  Double_t R = (uR < bR) ? uR : bR;
391 
392  Int_t iStation = CbmMuchAddress::GetStationIndex(sector->GetAddress());
393  //CbmMuchStationGem* station = (CbmMuchStationGem*)fStations->At(iStation);
394  // Check minimum and maximum allowed resolution
395  Double_t sigmaMax =
396  fSigmaXmax.at(iStation); //station->GetSigmaXmax(); //[cm]
397  Double_t sigmaMin =
398  fSigmaXmin.at(iStation); //station->GetSigmaXmin(); //[cm]
399  Double_t sigma = sector->GetSigmaX();
400  if (sigma > sigmaMin && sigma / 2. < sigmaMin) return false;
401  if (sigma > sigmaMax) return true;
402  // Check for occupancy
403  Double_t hitDensity = exp(fExp0.at(iStation) + fExp1.at(iStation) * R);
404  Double_t occupancyMax =
405  fOccupancyMax.at(iStation); //station->GetOccupancyMax();
406  Double_t sPad = secLx * secLy / 128.;
407  Double_t nPads =
408  (1.354
409  - 0.304
410  * TMath::Log2(sPad)); // number of pads fired by one track in average
411  //printf("nPads:%f\n",nPads);
412  Double_t occupancy = hitDensity * sPad * nPads;
413  if (occupancy > occupancyMax) return true;
414  return false;
415 }
416 // -------------------------------------------------------------------------
417 
418 // ----- Private method ShouldSegmentByY ---------------------------------
419 Bool_t CbmMuchSegmentAuto::ShouldSegmentByY(CbmMuchSectorRectangular* sector) {
420  Double_t secLx = sector->GetSize()[0];
421  Double_t secLy = sector->GetSize()[1];
422  Double_t secX = sector->GetPosition()[0];
423  Double_t secY = sector->GetPosition()[1];
424 
425  Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.)
426  + (secY + secLy / 2.) * (secY + secLy / 2.));
427  Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.)
428  + (secY + secLy / 2.) * (secY + secLy / 2.));
429  Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.)
430  + (secY - secLy / 2.) * (secY - secLy / 2.));
431  Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.)
432  + (secY - secLy / 2.) * (secY - secLy / 2.));
433 
434  Double_t uR = (ulR < urR) ? ulR : urR;
435  Double_t bR = (blR < brR) ? blR : brR;
436  Double_t R = (uR < bR) ? uR : bR;
437 
438  Int_t iStation = CbmMuchAddress::GetStationIndex(sector->GetAddress());
439  //CbmMuchStationGem* station = (CbmMuchStationGem*)fStations->At(iStation);
440  // Check minimum and maximum allowed resolution
441  Double_t sigmaMax =
442  fSigmaYmax.at(iStation); //station->GetSigmaYmax(); //[cm]
443  Double_t sigmaMin =
444  fSigmaYmin.at(iStation); //station->GetSigmaYmin(); //[cm]
445  Double_t sigma = sector->GetSigmaY();
446  if (sigma > sigmaMin && sigma / 2. < sigmaMin) return false;
447  if (sigma > sigmaMax) return true;
448  // Check for occupancy
449  Double_t hitDensity = exp(fExp0.at(iStation) + fExp1.at(iStation) * R);
450  Double_t occupancyMax =
451  fOccupancyMax.at(iStation); //station->GetOccupancyMax();
452  Double_t sPad = secLx * secLy / 128.;
453  Double_t nPads =
454  (1.354
455  - 0.304
456  * TMath::Log2(sPad)); // number of pads fired by one track in average
457  Double_t occupancy = hitDensity * sPad * nPads;
458  if (occupancy > occupancyMax) return true;
459  return false;
460 }
461 // -------------------------------------------------------------------------
462 
463 
464 // ----- Private method IntersectsHole -----------------------------------
465 Int_t CbmMuchSegmentAuto::IntersectsRad(CbmMuchSectorRectangular* sector,
466  Double_t radius) {
467  if (radius < 0) return 0;
468 
469  Int_t intersects = 0;
470 
471  Double_t secLx = sector->GetSize()[0];
472  Double_t secLy = sector->GetSize()[1];
473  Double_t secX = sector->GetPosition()[0];
474  Double_t secY = sector->GetPosition()[1];
475 
476  Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.)
477  + (secY + secLy / 2.) * (secY + secLy / 2.));
478  Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.)
479  + (secY + secLy / 2.) * (secY + secLy / 2.));
480  Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.)
481  + (secY - secLy / 2.) * (secY - secLy / 2.));
482  Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.)
483  + (secY - secLy / 2.) * (secY - secLy / 2.));
484 
485  if (ulR < radius) intersects++;
486  if (urR < radius) intersects++;
487  if (blR < radius) intersects++;
488  if (brR < radius) intersects++;
489 
490  if (intersects == 4) return 2;
491  if (intersects)
492  return 1;
493  else
494  return 0;
495 }
496 // -------------------------------------------------------------------------
497 
498 void CbmMuchSegmentAuto::Print(Option_t*) const {
499  printf("Segmentation written to file %s\n", fDigiFileName.Data());
500  Int_t nTotSectors = 0;
501  Int_t nTotChannels = 0;
502  Int_t nTotGems = 0;
503  Int_t nTotStraws = 0;
504  printf("====================================================================="
505  "============================\n");
506  for (Int_t iStation = 0; iStation < fStations->GetEntries(); ++iStation) {
507  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
508  Int_t nGems = 0;
509  Int_t nStraws = 0;
510  Int_t nSectors = 0;
511  Int_t nChannels = 0;
512  Double_t padMaxLx = std::numeric_limits<Double_t>::min();
513  Double_t padMinLx = std::numeric_limits<Double_t>::max();
514  Double_t padMaxLy = std::numeric_limits<Double_t>::min();
515  Double_t padMinLy = std::numeric_limits<Double_t>::max();
516  if (!station) continue;
517  for (Int_t iLayer = 0; iLayer < station->GetNLayers(); ++iLayer) {
518  CbmMuchLayer* layer = station->GetLayer(iLayer);
519  if (!layer) continue;
520  for (Int_t iSide = 0; iSide < 2; ++iSide) {
521  CbmMuchLayerSide* side = layer->GetSide(iSide);
522  if (!side) continue;
523  for (Int_t iModule = 0; iModule < side->GetNModules(); ++iModule) {
524  CbmMuchModule* mod = side->GetModule(iModule);
525  if (!mod) continue;
526  switch (mod->GetDetectorType()) {
527  case 1: { // GEMs
528  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
529  if (!module) continue;
530  nGems++;
531  nSectors += module->GetNSectors();
532  for (Int_t iSector = 0; iSector < module->GetNSectors();
533  ++iSector) {
534  CbmMuchSectorRectangular* sector =
535  (CbmMuchSectorRectangular*) module->GetSector(iSector);
536  if (!sector) continue;
537  Double_t padLx = sector->GetPadDx();
538  Double_t padLy = sector->GetPadDy();
539  if (padLx > padMaxLx) padMaxLx = padLx;
540  if (padLx < padMinLx) padMinLx = padLx;
541  if (padLy > padMaxLy) padMaxLy = padLy;
542  if (padLy < padMinLy) padMinLy = padLy;
543  nChannels += sector->GetNChannels();
544  }
545  break;
546  }
547  case 2: // Straw tubes
548  nStraws++;
549  break;
550  }
551  }
552  }
553  }
554  printf("Station %i:\n", iStation + 1);
555  printf(" GEM modules: %i\n", nGems);
556  if (nGems)
557  printf(" Sectors: %i, Pads: %i, Min.Pad size:%3.2fx%3.2f, Max.Pad "
558  "size:%3.2fx%3.2f\n",
559  nSectors,
560  nChannels,
561  padMinLx,
562  padMinLy,
563  padMaxLx,
564  padMaxLy);
565  printf(" Straw modules: %i\n", nStraws);
566  nTotSectors += nSectors;
567  nTotChannels += nChannels;
568  nTotGems += nGems;
569  nTotStraws += nStraws;
570  }
571  printf("---------------------------------------------------------------------"
572  "----------------------------\n");
573  printf(" Summary: \n GEM modules: %i\n Sectors: %i, Pads: %i\n "
574  "Straw modules: %i\n",
575  nTotGems,
576  nTotSectors,
577  nTotChannels,
578  nTotStraws);
579  printf("====================================================================="
580  "============================\n");
581 }
582 
583 void CbmMuchSegmentAuto::DrawSegmentation() {
584  string digifile(fDigiFileName);
585  Int_t startIndex = digifile.size() - 4;
586  string txtfile = digifile.erase(startIndex, 4);
587  txtfile.append("txt");
588  FILE* outfile;
589  outfile = fopen(txtfile.c_str(), "w");
590 
591  // Int_t colors[] = {kGreen, kBlue, kViolet, kRed, kYellow, kOrange, kMagenta, kCyan, kSpring, kPink, kAzure, kTeal,
592  // kGreen+10, kBlue+10, kViolet+10, kRed+10, kYellow+10, kOrange+10, kMagenta+10, kCyan+10, kSpring+10, kPink+10, kAzure+10, kTeal+10};
593 
594  Double_t secMinLx = std::numeric_limits<Double_t>::max();
595  Double_t secMinLy = std::numeric_limits<Double_t>::max();
596  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
597  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
598  CbmMuchLayer* layer = station->GetLayer(0);
599  for (Int_t iSide = 1; iSide >= 0; iSide--) {
600  CbmMuchLayerSide* side = layer->GetSide(iSide);
601  for (Int_t iModule = 0; iModule < side->GetNModules(); ++iModule) {
602  CbmMuchModule* mod = side->GetModule(iModule);
603  if (mod->GetDetectorType() != 1) continue;
604  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
605  for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
606  CbmMuchSectorRectangular* sector =
607  (CbmMuchSectorRectangular*) module->GetSector(iSector);
608  if (sector->GetSize()[0] < secMinLx) secMinLx = sector->GetSize()[0];
609  if (sector->GetSize()[1] < secMinLy) secMinLy = sector->GetSize()[1];
610  }
611  } // modules
612  } // sides
613  }
614 
615  for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
616  fprintf(outfile,
617  "=================================================================="
618  "==========\n");
619  fprintf(outfile, "Station %i\n", iStation + 1);
620  fprintf(outfile,
621  "Sector size, cm Sector position, cm Number of pads Side "
622  "Pad size, cm\n");
623  fprintf(outfile,
624  "------------------------------------------------------------------"
625  "----------\n");
626  TCanvas* c1 = new TCanvas(Form("station%i", iStation + 1),
627  Form("station%i", iStation + 1),
628  800,
629  800);
630  c1->SetFillColor(0);
631  c1->Range(-250, -250, 250, 250);
632  CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
633  CbmMuchLayer* layer = station->GetLayer(0);
634  for (Int_t iSide = 1; iSide >= 0; iSide--) {
635  CbmMuchLayerSide* layerSide = layer->GetSide(iSide);
636  for (Int_t iModule = 0; iModule < layerSide->GetNModules(); ++iModule) {
637  CbmMuchModule* mod = layerSide->GetModule(iModule);
638  mod->SetFillStyle(0);
639  mod->Draw();
640  if (mod->GetDetectorType() != 1) continue;
641  CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
642  for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
643  CbmMuchSectorRectangular* sector =
644  (CbmMuchSectorRectangular*) module->GetSector(iSector);
645 
646  // Int_t i = Int_t((sector->GetSize()[0]+1e-3)/secMinLx) - 1;
647  // Int_t j = Int_t((sector->GetSize()[1]+1e-3)/secMinLy) - 1;
648  // TODO
649  // sector->SetFillColor(iSide ? TColor::GetColorDark(colors[i+j]) : colors[i+j]);
650  // sector->Draw("f");
651  // sector->Draw();
652  const char* side = iSide ? "Back" : "Front";
653  fprintf(outfile,
654  "%-4.2fx%-10.2f %-6.2fx%-12.2f %-14i %-5s ",
655  sector->GetSize()[0],
656  sector->GetSize()[1],
657  sector->GetPosition()[0],
658  sector->GetPosition()[1],
659  sector->GetNChannels(),
660  side);
661  //TODO fprintf(outfile, "%-4.2fx%-4.2f\n", sector->GetDx(), sector->GetDy());
662  } // sectors
663  } // modules
664  } // sides
665 
666  // Draw a hole
667  TArc* arc = new TArc(0., 0., station->GetRmin());
668  arc->Draw();
669 
670  c1->Print(
671  Form("%s/station%i.eps", gSystem->DirName(fDigiFileName), iStation + 1));
672  c1->Print(
673  Form("%s/station%i.png", gSystem->DirName(fDigiFileName), iStation + 1));
674  } //stations
675  fclose(outfile);
676 }
677 
678 ClassImp(CbmMuchSegmentAuto)
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:134
CbmMuchLayerSide.h
CbmMuchModule.h
CbmMuchSegmentAuto::ShouldSegmentByX
Bool_t ShouldSegmentByX(CbmMuchSectorRectangular *sector)
Definition: CbmMuchSegmentAuto.cxx:373
CbmMuchSegmentAuto::SetSigmaMax
void SetSigmaMax(Double_t *sigmaXmax, Double_t *sigmaYmax)
Definition: CbmMuchSegmentAuto.cxx:104
CbmMuchSegmentAuto::FinishTask
virtual void FinishTask()
Definition: CbmMuchSegmentAuto.cxx:182
CbmMuchSegmentAuto::SetParContainers
virtual void SetParContainers()
Definition: CbmMuchSegmentAuto.cxx:118
CbmMuchModuleGem::AddSector
void AddSector(CbmMuchSector *sector)
Definition: CbmMuchModuleGem.h:67
CbmMuchPoint
Definition: CbmMuchPoint.h:21
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
CbmMuchSegmentAuto::fHistHitDensity
TH1D ** fHistHitDensity
Definition: CbmMuchSegmentAuto.h:53
CbmMuchSegmentAuto
Definition: CbmMuchSegmentAuto.h:29
CbmMuchStation::GetLayer
CbmMuchLayer * GetLayer(Int_t iLayer) const
Definition: CbmMuchStation.h:58
CbmMuchStation
Definition: CbmMuchStation.h:22
CbmMuchSegmentAuto::fSigmaXmin
std::vector< Double_t > fSigmaXmin
Definition: CbmMuchSegmentAuto.h:64
CbmMuchSectorRectangular::GetPadDy
Double_t GetPadDy() const
Definition: CbmMuchSectorRectangular.h:40
CbmMuchSectorRectangular::GetPadDx
Double_t GetPadDx() const
Definition: CbmMuchSectorRectangular.h:39
CbmMuchSectorRectangular::GetSigmaX
Double_t GetSigmaX() const
Definition: CbmMuchSectorRectangular.h:41
CbmMuchSectorRectangular::GetSize
TVector3 GetSize() const
Definition: CbmMuchSectorRectangular.h:36
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmMuchSegmentAuto::SetOccupancyMax
void SetOccupancyMax(Double_t *occupancyMax)
Definition: CbmMuchSegmentAuto.cxx:111
CbmGeoMuchPar::GetStations
TObjArray * GetStations()
Definition: CbmGeoMuchPar.h:40
CbmMuchModuleGem::GetSector
CbmMuchSector * GetSector(Int_t address)
Definition: CbmMuchModuleGem.cxx:71
CbmMuchSegmentAuto::CbmMuchSegmentAuto
CbmMuchSegmentAuto()
Definition: CbmMuchSegmentAuto.cxx:48
CbmMuchSector::GetNChannels
Int_t GetNChannels() const
Definition: CbmMuchSector.h:31
CbmMuchSegmentAuto::InitLayerSide
void InitLayerSide(CbmMuchLayerSide *layerSide)
Definition: CbmMuchSegmentAuto.cxx:249
CbmGeoMuchPar
Definition: CbmGeoMuchPar.h:25
CbmMuchSegmentAuto::fStations
TObjArray * fStations
Definition: CbmMuchSegmentAuto.h:55
CbmMuchLayer::GetSide
CbmMuchLayerSide * GetSide(Bool_t side)
Definition: CbmMuchLayer.h:49
CbmMuchSegmentAuto::SetSigmaMin
void SetSigmaMin(Double_t *sigmaXmin, Double_t *sigmaYmin)
Definition: CbmMuchSegmentAuto.cxx:98
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:33
CbmMuchModule::GetPosition
TVector3 GetPosition() const
Definition: CbmMuchModule.h:51
CbmMuchSegmentAuto::fExp1
std::vector< Double_t > fExp1
Definition: CbmMuchSegmentAuto.h:62
CbmMuchAddress::GetModuleIndex
static Int_t GetModuleIndex(Int_t address)
Definition: CbmMuchAddress.h:112
CbmMuchModule::GetSize
TVector3 GetSize() const
Definition: CbmMuchModule.h:50
CbmMuchSegmentAuto::fGeoPar
CbmGeoMuchPar * fGeoPar
Definition: CbmMuchSegmentAuto.h:58
CbmMuchSectorRectangular.h
CbmMuchSegmentAuto::Print
void Print(Option_t *="") const
Definition: CbmMuchSegmentAuto.cxx:498
CbmMuchSegmentAuto.h
CbmMuchSegmentAuto::fSigmaXmax
std::vector< Double_t > fSigmaXmax
Definition: CbmMuchSegmentAuto.h:66
CbmMuchPoint.h
h
Data class with information on a STS local track.
CbmGeoMuchPar.h
CbmMuchSegmentAuto::DrawSegmentation
void DrawSegmentation()
Definition: CbmMuchSegmentAuto.cxx:583
CbmMuchStation::GetRmax
Double_t GetRmax() const
Definition: CbmMuchStation.h:52
CbmMuchLayerSide::GetModule
CbmMuchModule * GetModule(Int_t iModule) const
Definition: CbmMuchLayerSide.h:52
CbmMuchModule::GetCutRadius
Double_t GetCutRadius() const
Definition: CbmMuchModule.h:49
CbmMuchSegmentAuto::~CbmMuchSegmentAuto
virtual ~CbmMuchSegmentAuto()
Definition: CbmMuchSegmentAuto.cxx:86
CbmMuchModuleGem
Definition: CbmMuchModuleGem.h:24
CbmMuchSectorRectangular
Definition: CbmMuchSectorRectangular.h:24
CbmMuchSegmentAuto::fPoints
TClonesArray * fPoints
Definition: CbmMuchSegmentAuto.h:52
CbmMuchSegmentAuto::ShouldSegmentByY
Bool_t ShouldSegmentByY(CbmMuchSectorRectangular *sector)
Definition: CbmMuchSegmentAuto.cxx:419
CbmMuchSegmentAuto::fOccupancyMax
std::vector< Double_t > fOccupancyMax
Definition: CbmMuchSegmentAuto.h:68
CbmMuchLayerSide
Definition: CbmMuchLayerSide.h:22
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmMuchModule
Definition: CbmMuchModule.h:24
CbmMuchLayer.h
CbmMuchAddress::GetLayerIndex
static Int_t GetLayerIndex(Int_t address)
Definition: CbmMuchAddress.h:106
CbmMuchAddress::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchAddress.h:103
CbmMuchModule::GetDetectorId
Int_t GetDetectorId() const
Definition: CbmMuchModule.h:48
CbmMuchSectorRectangular::GetPosition
TVector3 GetPosition() const
Definition: CbmMuchSectorRectangular.h:35
CbmMuchSegmentAuto::fSigmaYmax
std::vector< Double_t > fSigmaYmax
Definition: CbmMuchSegmentAuto.h:67
CbmMuchModule::GetDetectorType
Int_t GetDetectorType() const
Definition: CbmMuchModule.h:52
CbmMuchSectorRectangular::GetSigmaY
Double_t GetSigmaY() const
Definition: CbmMuchSectorRectangular.h:42
CbmMuchSegmentAuto::SegmentModule
void SegmentModule(CbmMuchModuleGem *module)
Definition: CbmMuchSegmentAuto.cxx:263
CbmMuchAddress.h
CbmMuchSegmentAuto::fExp0
std::vector< Double_t > fExp0
Definition: CbmMuchSegmentAuto.h:60
CbmMuchSegmentAuto::IntersectsRad
Int_t IntersectsRad(CbmMuchSectorRectangular *sector, Double_t radius)
Definition: CbmMuchSegmentAuto.cxx:465
CbmMuchSegmentAuto::Init
virtual InitStatus Init()
Definition: CbmMuchSegmentAuto.cxx:129
CbmMuchSegmentAuto::fNStations
Int_t fNStations
Definition: CbmMuchSegmentAuto.h:54
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmMuchSegmentAuto::fEvents
Int_t fEvents
Definition: CbmMuchSegmentAuto.h:51
CbmMuchStation::GetNLayers
Int_t GetNLayers() const
Definition: CbmMuchStation.h:50
CbmMuchSegmentAuto::SetNStations
void SetNStations(Int_t nStations)
Definition: CbmMuchSegmentAuto.cxx:89
CbmMuchSegmentAuto::fSigmaYmin
std::vector< Double_t > fSigmaYmin
Definition: CbmMuchSegmentAuto.h:65
CbmMuchLayer::GetSideF
CbmMuchLayerSide * GetSideF()
Definition: CbmMuchLayer.h:47
CbmMuchSector::GetAddress
UInt_t GetAddress() const
Definition: CbmMuchSector.h:27
CbmMuchLayer
Definition: CbmMuchLayer.h:21
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
CbmMuchLayer::GetSideB
CbmMuchLayerSide * GetSideB()
Definition: CbmMuchLayer.h:48
CbmMuchStation::GetRmin
Double_t GetRmin() const
Definition: CbmMuchStation.h:51
CbmMuchSegmentAuto::fDigiFileName
TString fDigiFileName
Definition: CbmMuchSegmentAuto.h:57
CbmMuchStation.h
CbmMuchSegmentAuto::Exec
virtual void Exec(Option_t *option)
Definition: CbmMuchSegmentAuto.cxx:156
CbmMuchSegmentAuto::SegmentSector
void SegmentSector(CbmMuchModuleGem *module, CbmMuchSectorRectangular *sector)
Definition: CbmMuchSegmentAuto.cxx:294
CbmMuchModuleGem::GetNSectors
Int_t GetNSectors() const
Definition: CbmMuchModuleGem.h:57
CbmMuchModuleGem.h
CbmMuchAddress::GetLayerSideIndex
static Int_t GetLayerSideIndex(Int_t address)
Definition: CbmMuchAddress.h:109
CbmMuchLayerSide::GetNModules
Int_t GetNModules() const
Definition: CbmMuchLayerSide.h:47