CbmMuchSegmentAuto.cxx

Go to the documentation of this file.

/* Copyright (C) 2007-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Mikhail Ryzhinskiy [committer], Florian Uhlig */
 
#include "CbmMuchSegmentAuto.h"
 
#include "CbmGeoMuchPar.h"             // for CbmGeoMuchPar
#include "CbmMuchAddress.h"            // for CbmMuchAddress
#include "CbmMuchLayer.h"              // for CbmMuchLayer
#include "CbmMuchLayerSide.h"          // for CbmMuchLayerSide
#include "CbmMuchModule.h"             // for CbmMuchModule
#include "CbmMuchModuleGem.h"          // for CbmMuchModuleGem
#include "CbmMuchPoint.h"              // for CbmMuchPoint
#include "CbmMuchSectorRectangular.h"  // for CbmMuchSectorRectangular
#include "CbmMuchStation.h"            // for CbmMuchStation
 
#include <FairRootManager.h>  // for FairRootManager
#include <FairRuntimeDb.h>    // for FairRuntimeDb
#include <FairTask.h>         // for FairTask, InitStatus, kSUCCESS
 
#include <TArc.h>          // for TArc
#include <TAxis.h>         // for TAxis
#include <TCanvas.h>       // for TCanvas
#include <TClonesArray.h>  // for TClonesArray
#include <TF1.h>           // for TF1
#include <TFile.h>         // for TFile, gFile
#include <TH1.h>           // for TH1D
#include <TMath.h>         // for Sqrt, Log2, Pi
#include <TMathBase.h>     // for Abs
#include <TObjArray.h>     // for TObjArray
#include <TStyle.h>        // for TStyle, gStyle
#include <TSystem.h>       // for TSystem, gSystem
#include <TVector3.h>      // for TVector3
 
#include <cassert>  // for assert
#include <iosfwd>   // for string
#include <limits>   // for numeric_limits
 
#include <math.h>   // for exp
#include <stdio.h>  // for printf, fprintf, fclose
 
using std::string;
 
// -----   Default constructor   -------------------------------------------
CbmMuchSegmentAuto::CbmMuchSegmentAuto()
  : FairTask()
  , fEvents(0)
  , fPoints(nullptr)
  , fHistHitDensity()
  , fNStations(0)
  , fStations(nullptr)
  , fDigiFileName()
  , fGeoPar(nullptr)
  , fExp0()
  , fExp1()
  , fSigmaXmin()
  , fSigmaYmin()
  , fSigmaXmax()
  , fSigmaYmax()
  , fOccupancyMax()
{
}
// -------------------------------------------------------------------------
 
// -----   Standard constructor   ------------------------------------------
CbmMuchSegmentAuto::CbmMuchSegmentAuto(const char* digiFileName)
  : FairTask()
  , fEvents(0)
  , fPoints(nullptr)
  , fHistHitDensity()
  , fNStations(0)
  , fStations(nullptr)
  , fDigiFileName(digiFileName)
  , fGeoPar(nullptr)
  , fExp0()
  , fExp1()
  , fSigmaXmin()
  , fSigmaYmin()
  , fSigmaXmax()
  , fSigmaYmax()
  , fOccupancyMax()
{
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
CbmMuchSegmentAuto::~CbmMuchSegmentAuto() {}
// -------------------------------------------------------------------------
 
void CbmMuchSegmentAuto::SetNStations(Int_t nStations)
{
  fNStations = nStations;
  fSigmaXmin.resize(fNStations);
  fSigmaYmin.resize(fNStations);
  fSigmaXmax.resize(fNStations);
  fSigmaYmax.resize(fNStations);
  fOccupancyMax.resize(fNStations);
}
 
void CbmMuchSegmentAuto::SetSigmaMin(Double_t* sigmaXmin, Double_t* sigmaYmin)
{
  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
    fSigmaXmin.at(iStation) = sigmaXmin[iStation];
    fSigmaYmin.at(iStation) = sigmaYmin[iStation];
  }
}
void CbmMuchSegmentAuto::SetSigmaMax(Double_t* sigmaXmax, Double_t* sigmaYmax)
{
  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
    fSigmaXmax.at(iStation) = sigmaXmax[iStation];
    fSigmaYmax.at(iStation) = sigmaYmax[iStation];
  }
}
 
void CbmMuchSegmentAuto::SetOccupancyMax(Double_t* occupancyMax)
{
  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
    fOccupancyMax.at(iStation) = occupancyMax[iStation];
  }
}
 
// -----   Private method SetParContainers  --------------------------------
void CbmMuchSegmentAuto::SetParContainers()
{
  // Get runtime database
  FairRuntimeDb* db = FairRuntimeDb::instance();
  if (!db) Fatal("SetParContainers", "No runtime database");
  // Get MUCH geometry param
  fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
}
// -------------------------------------------------------------------------
 
 
// -----   Private method Init  --------------------------------------------
InitStatus CbmMuchSegmentAuto::Init()
{
  FairRootManager* fManager = FairRootManager::Instance();
  fPoints                   = (TClonesArray*) fManager->GetObject("MuchPoint");
  fEvents                   = 0;
 
  fStations = fGeoPar->GetStations();
  if (!fStations) Fatal("Init", "No input array of MUCH stations.");
  if (fNStations != fStations->GetEntriesFast()) Fatal("Init", "Incorrect number of stations set.");
  //fNStations = fStations->GetEntriesFast();
  fHistHitDensity = new TH1D*[fNStations];
 
  for (Int_t i = 0; i < fNStations; i++) {
    //    CbmMuchStation* station = (CbmMuchStation*) fStations->At(i);
    fHistHitDensity[i] = new TH1D(Form("hStation%i", i + 1), Form("Station %i", i + 1), 110, 0, 220);
    fHistHitDensity[i]->GetXaxis()->SetTitle("r, cm");
    fHistHitDensity[i]->GetYaxis()->SetTitle("hits/(event#timescm^{2})");
    fHistHitDensity[i]->SetLineColor(4);
    fHistHitDensity[i]->SetLineWidth(2);
    fHistHitDensity[i]->GetYaxis()->SetTitleOffset(1.27);
  }
  return kSUCCESS;
}
// -------------------------------------------------------------------------
 
// -------------------------------------------------------------------------
void CbmMuchSegmentAuto::Exec(Option_t*)
{
  fEvents++;
  printf("Event: %i\n", fEvents);
 
  gStyle->SetOptStat(0);
  for (Int_t iPoint = 0; iPoint < fPoints->GetEntriesFast(); iPoint++) {
    CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(iPoint);
    if (!point) continue;
    Int_t iStation = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
    //        CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
    //        if (station->GetDetectorType()==2) continue;
 
    Int_t iLayer = CbmMuchAddress::GetLayerIndex(point->GetDetectorID());
    //    printf("iStation = %i\n", iStation);
    //    printf("detId = %qd\n", point->GetDetectorID());
    //    printf("iPoint = %i\n", iPoint);
    assert(iStation >= 0 && iStation < fNStations);
    if (iLayer) continue;
    TVector3 pos;
    point->Position(pos);
    ((TH1D*) fHistHitDensity[iStation])->Fill(pos.Pt());
  }
}
// -------------------------------------------------------------------------
 
// -------------------------------------------------------------------------
void CbmMuchSegmentAuto::FinishTask()
{
  // Create normalization histo
  TH1D* hNorm      = new TH1D("hNorm", "", 110, 0, 220);
  Double_t binSize = 2.;
  for (Int_t l = 0; l < 100; l++) {
    Double_t R1 = l * binSize;
    Double_t R2 = (l + 1) * binSize;
    Double_t s  = TMath::Pi() * (R2 * R2 - R1 * R1);
    hNorm->SetBinContent(l + 1, s * fEvents);
  }
 
  for (Int_t i = 0; i < fNStations; i++) {
    CbmMuchStation* station = (CbmMuchStation*) fStations->At(i);
    TH1D* h                 = fHistHitDensity[i];
    h->Divide(hNorm);
    TCanvas* c1 = new TCanvas(Form("cStation%i", i + 1), Form("Station %i", i + 1), 10, 10, 500, 500);
    c1->SetLogy();
    c1->SetLeftMargin(0.12);
    c1->SetRightMargin(0.08);
    TF1* f1 = new TF1("f1", "exp([0] + [1]*x)");
    f1->SetParameter(0, 0.5);
    f1->SetParameter(1, -0.1);
    h->Fit("f1", "QLL", "", station->GetRmin(), station->GetRmax());
    Double_t exp0 = h->GetFunction("f1")->GetParameter(0);
    Double_t exp1 = h->GetFunction("f1")->GetParameter(1);
    fExp0.push_back(exp0);
    fExp1.push_back(exp1);
 
    h->Draw();
    c1->Print(Form("%s/hd_station%i.eps", gSystem->DirName(fDigiFileName), i + 1));
    c1->Print(Form("%s/hd_station%i.png", gSystem->DirName(fDigiFileName), i + 1));
  }
 
  for (Int_t i = 0; i < fNStations; i++) {
    CbmMuchStation* station = (CbmMuchStation*) fStations->At(i);
    //if (station->GetDetectorType()==2) continue;
 
    Int_t nLayers = station->GetNLayers();
    for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
      CbmMuchLayer* layer = station->GetLayer(iLayer);
      if (!layer) Fatal("Init", "Incomplete layers array.");
      // Get layer sides
      InitLayerSide(layer->GetSideF());
      InitLayerSide(layer->GetSideB());
    }
    printf("Station%i segmented\n", i + 1);
  }
 
  TFile* oldFile     = gFile;
  TDirectory* oldDir = gDirectory;
 
  // Save parameters
  TFile* f = new TFile(fDigiFileName, "RECREATE");
  fStations->Write("stations", 1);
 
  f->Close();
  gFile      = oldFile;
  gDirectory = oldDir;
 
  DrawSegmentation();
 
  Print();
}
// -------------------------------------------------------------------------
 
 
// -----   Private method InitLayerSide  -----------------------------------
void CbmMuchSegmentAuto::InitLayerSide(CbmMuchLayerSide* layerSide)
{
  if (!layerSide) Fatal("Init", "Incomplete layer sides array.");
  Int_t nModules = layerSide->GetNModules();
  for (Int_t iModule = 0; iModule < nModules; iModule++) {
    CbmMuchModule* mod = layerSide->GetModule(iModule);
    if (mod->GetDetectorType() != 1) continue;
    CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
    SegmentModule(module);
  }
}
// -------------------------------------------------------------------------
 
 
// -----   Private method SegmentModule  -----------------------------------
void CbmMuchSegmentAuto::SegmentModule(CbmMuchModuleGem* module)
{
  TVector3 modSize     = module->GetSize();
  Double_t modLx       = modSize.X();
  Double_t modLy       = modSize.Y();
  Double_t modLz       = modSize.Z();
  TVector3 modPosition = module->GetPosition();
  Double_t modX        = modPosition.X();
  Double_t modY        = modPosition.Y();
  Double_t modZ        = modPosition.Z();
 
  Bool_t result    = modLx > modLy;
  Int_t iRatio     = (result) ? (Int_t)((modLx + 1e-3) / modLy) : (Int_t)((modLy + 1e-3) / modLx);
  Int_t detectorId = module->GetDetectorId();
  Double_t secLx   = (result) ? modLx / iRatio : modLx;
  Double_t secLy   = (result) ? modLy : modLy / iRatio;
  for (Int_t i = 0; i < iRatio; i++) {
    Double_t secX = (result) ? modX - modLx / 2. + (i + 0.5) * secLx : modX;
    Double_t secY = (result) ? modY : modY - modLy / 2. + (i + 0.5) * secLy;
    Int_t iSector = module->GetNSectors();
 
    TVector3 sectorPosition(secX, secY, modZ);
    TVector3 sectorSize(secLx, secLy, modLz);
    CbmMuchSectorRectangular* sector =
      new CbmMuchSectorRectangular(detectorId, iSector, sectorPosition, sectorSize, 8, 16);
    SegmentSector(module, sector);
  }
}
// -------------------------------------------------------------------------
 
// -----  Private method SegmentSector  ------------------------------------
void CbmMuchSegmentAuto::SegmentSector(CbmMuchModuleGem* module, CbmMuchSectorRectangular* sector)
{
  TVector3 secSize     = sector->GetSize();
  TVector3 secPosition = sector->GetPosition();
  Int_t detectorId     = module->GetDetectorId();
  Int_t iStation       = CbmMuchAddress::GetStationIndex(detectorId);
  Int_t iLayer         = CbmMuchAddress::GetLayerIndex(detectorId);
  Int_t iSide          = CbmMuchAddress::GetLayerSideIndex(detectorId);
  Int_t iModule        = CbmMuchAddress::GetModuleIndex(detectorId);
  Double_t secLx       = secSize.X();
  Double_t secLy       = secSize.Y();
  Double_t secLz       = secSize.Z();
  Double_t secX        = secPosition.X();
  Double_t secY        = secPosition.Y();
  Double_t secZ        = secPosition.Z();
 
  assert(CbmMuchAddress::GetStationIndex(sector->GetAddress()) == iStation);
  assert(CbmMuchAddress::GetLayerIndex(sector->GetAddress()) == iLayer);
  assert(CbmMuchAddress::GetLayerSideIndex(sector->GetAddress()) == iSide);
  assert(CbmMuchAddress::GetModuleIndex(sector->GetAddress()) == iModule);
 
  Bool_t result1 = ShouldSegmentByX(sector);
  Bool_t result2 = ShouldSegmentByY(sector);
 
  if (result1 || result2) { delete sector; }
  else {
    CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
    Double_t rMax           = station->GetRmax();
    if ((IntersectsRad(sector, module->GetCutRadius()) == 2) || !IntersectsRad(sector, rMax)) {
      delete sector;
      return;
    }
    module->AddSector(sector);
    return;
  }
 
  Int_t iSector;
  Double_t newSecX, newSecY, newSecLx, newSecLy;
  Bool_t equal = TMath::Abs(secLx - secLy) < 1e-5;
  Bool_t res   = secLx > secLy;
  TVector3 position, size;
  if (result1 && result2) {
    if (equal) { res = kTRUE; }
    for (Int_t i = 0; i < 2; ++i) {
      iSector  = module->GetNSectors();  //CbmMuchGeoScheme::GetDetIdFromModule(moduleId, module->GetNSectors());
      newSecLx = res ? secLx / 2. : secLx;
      newSecLy = res ? secLy : secLy / 2.;
      newSecX  = res ? secX + (i - 0.5) * newSecLx : secX;
      newSecY  = res ? secY : secY + (i - 0.5) * newSecLy;
      position.SetXYZ(newSecX, newSecY, secZ);
      size.SetXYZ(newSecLx, newSecLy, secLz);
      SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, position, size, 8, 16));
    }
  }
  else if (result1 || result2) {
    for (Int_t i = 0; i < 2; i++) {
      iSector  = module->GetNSectors();  //CbmMuchGeoScheme::GetDetIdFromModule(moduleId, module->GetNSectors());
      newSecLx = result1 ? secLx / 2. : secLx;
      newSecLy = result1 ? secLy : secLy / 2;
      newSecX  = result1 ? secX + (i - 0.5) * newSecLx : secX;
      newSecY  = result1 ? secY : secY + (i - 0.5) * newSecLy;
      position.SetXYZ(newSecX, newSecY, secZ);
      size.SetXYZ(newSecLx, newSecLy, secLz);
      SegmentSector(module, new CbmMuchSectorRectangular(detectorId, iSector, position, size, 8, 16));
    }
  }
}
// -------------------------------------------------------------------------
 
// -----  Private method ShouldSegmentByX  ---------------------------------
Bool_t CbmMuchSegmentAuto::ShouldSegmentByX(CbmMuchSectorRectangular* sector)
{
  Double_t secLx = sector->GetSize()[0];
  Double_t secLy = sector->GetSize()[1];
  Double_t secX  = sector->GetPosition()[0];
  Double_t secY  = sector->GetPosition()[1];
 
  Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
  Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
  Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
  Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
 
  Double_t uR = (ulR < urR) ? ulR : urR;
  Double_t bR = (blR < brR) ? blR : brR;
  Double_t R  = (uR < bR) ? uR : bR;
 
  Int_t iStation = CbmMuchAddress::GetStationIndex(sector->GetAddress());
  //CbmMuchStationGem* station = (CbmMuchStationGem*)fStations->At(iStation);
  // Check minimum and maximum allowed resolution
  Double_t sigmaMax = fSigmaXmax.at(iStation);  //station->GetSigmaXmax(); //[cm]
  Double_t sigmaMin = fSigmaXmin.at(iStation);  //station->GetSigmaXmin(); //[cm]
  Double_t sigma    = sector->GetSigmaX();
  if (sigma > sigmaMin && sigma / 2. < sigmaMin) return false;
  if (sigma > sigmaMax) return true;
  // Check for occupancy
  Double_t hitDensity   = exp(fExp0.at(iStation) + fExp1.at(iStation) * R);
  Double_t occupancyMax = fOccupancyMax.at(iStation);  //station->GetOccupancyMax();
  Double_t sPad         = secLx * secLy / 128.;
  Double_t nPads        = (1.354 - 0.304 * TMath::Log2(sPad));  // number of pads fired by one track in average
  //printf("nPads:%f\n",nPads);
  Double_t occupancy = hitDensity * sPad * nPads;
  if (occupancy > occupancyMax) return true;
  return false;
}
// -------------------------------------------------------------------------
 
// -----  Private method ShouldSegmentByY  ---------------------------------
Bool_t CbmMuchSegmentAuto::ShouldSegmentByY(CbmMuchSectorRectangular* sector)
{
  Double_t secLx = sector->GetSize()[0];
  Double_t secLy = sector->GetSize()[1];
  Double_t secX  = sector->GetPosition()[0];
  Double_t secY  = sector->GetPosition()[1];
 
  Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
  Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
  Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
  Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
 
  Double_t uR = (ulR < urR) ? ulR : urR;
  Double_t bR = (blR < brR) ? blR : brR;
  Double_t R  = (uR < bR) ? uR : bR;
 
  Int_t iStation = CbmMuchAddress::GetStationIndex(sector->GetAddress());
  //CbmMuchStationGem* station = (CbmMuchStationGem*)fStations->At(iStation);
  // Check minimum and maximum allowed resolution
  Double_t sigmaMax = fSigmaYmax.at(iStation);  //station->GetSigmaYmax(); //[cm]
  Double_t sigmaMin = fSigmaYmin.at(iStation);  //station->GetSigmaYmin(); //[cm]
  Double_t sigma    = sector->GetSigmaY();
  if (sigma > sigmaMin && sigma / 2. < sigmaMin) return false;
  if (sigma > sigmaMax) return true;
  // Check for occupancy
  Double_t hitDensity   = exp(fExp0.at(iStation) + fExp1.at(iStation) * R);
  Double_t occupancyMax = fOccupancyMax.at(iStation);  //station->GetOccupancyMax();
  Double_t sPad         = secLx * secLy / 128.;
  Double_t nPads        = (1.354 - 0.304 * TMath::Log2(sPad));  // number of pads fired by one track in average
  Double_t occupancy    = hitDensity * sPad * nPads;
  if (occupancy > occupancyMax) return true;
  return false;
}
// -------------------------------------------------------------------------
 
 
// -----  Private method IntersectsHole  -----------------------------------
Int_t CbmMuchSegmentAuto::IntersectsRad(CbmMuchSectorRectangular* sector, Double_t radius)
{
  if (radius < 0) return 0;
 
  Int_t intersects = 0;
 
  Double_t secLx = sector->GetSize()[0];
  Double_t secLy = sector->GetSize()[1];
  Double_t secX  = sector->GetPosition()[0];
  Double_t secY  = sector->GetPosition()[1];
 
  Double_t ulR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
  Double_t urR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY + secLy / 2.) * (secY + secLy / 2.));
  Double_t blR = TMath::Sqrt((secX - secLx / 2.) * (secX - secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
  Double_t brR = TMath::Sqrt((secX + secLx / 2.) * (secX + secLx / 2.) + (secY - secLy / 2.) * (secY - secLy / 2.));
 
  if (ulR < radius) intersects++;
  if (urR < radius) intersects++;
  if (blR < radius) intersects++;
  if (brR < radius) intersects++;
 
  if (intersects == 4) return 2;
  if (intersects) return 1;
  else
    return 0;
}
// -------------------------------------------------------------------------
 
void CbmMuchSegmentAuto::Print(Option_t*) const
{
  printf("Segmentation written to file %s\n", fDigiFileName.Data());
  Int_t nTotSectors  = 0;
  Int_t nTotChannels = 0;
  Int_t nTotGems     = 0;
  Int_t nTotStraws   = 0;
  printf("====================================================================="
         "============================\n");
  for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
    CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
    Int_t nGems             = 0;
    Int_t nStraws           = 0;
    Int_t nSectors          = 0;
    Int_t nChannels         = 0;
    Double_t padMaxLx       = std::numeric_limits<Double_t>::min();
    Double_t padMinLx       = std::numeric_limits<Double_t>::max();
    Double_t padMaxLy       = std::numeric_limits<Double_t>::min();
    Double_t padMinLy       = std::numeric_limits<Double_t>::max();
    if (!station) continue;
    for (Int_t iLayer = 0; iLayer < station->GetNLayers(); ++iLayer) {
      CbmMuchLayer* layer = station->GetLayer(iLayer);
      if (!layer) continue;
      for (Int_t iSide = 0; iSide < 2; ++iSide) {
        CbmMuchLayerSide* side = layer->GetSide(iSide);
        if (!side) continue;
        for (Int_t iModule = 0; iModule < side->GetNModules(); ++iModule) {
          CbmMuchModule* mod = side->GetModule(iModule);
          if (!mod) continue;
          switch (mod->GetDetectorType()) {
            case 1: {  // GEMs
              CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
              if (!module) continue;
              nGems++;
              nSectors += module->GetNSectors();
              for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
                CbmMuchSectorRectangular* sector = (CbmMuchSectorRectangular*) module->GetSector(iSector);
                if (!sector) continue;
                Double_t padLx = sector->GetPadDx();
                Double_t padLy = sector->GetPadDy();
                if (padLx > padMaxLx) padMaxLx = padLx;
                if (padLx < padMinLx) padMinLx = padLx;
                if (padLy > padMaxLy) padMaxLy = padLy;
                if (padLy < padMinLy) padMinLy = padLy;
                nChannels += sector->GetNChannels();
              }
              break;
            }
            case 2:  // Straw tubes
              nStraws++;
              break;
          }
        }
      }
    }
    printf("Station %i:\n", iStation + 1);
    printf("   GEM modules: %i\n", nGems);
    if (nGems)
      printf("      Sectors: %i,  Pads: %i, Min.Pad size:%3.2fx%3.2f, Max.Pad "
             "size:%3.2fx%3.2f\n",
             nSectors, nChannels, padMinLx, padMinLy, padMaxLx, padMaxLy);
    printf("   Straw modules: %i\n", nStraws);
    nTotSectors += nSectors;
    nTotChannels += nChannels;
    nTotGems += nGems;
    nTotStraws += nStraws;
  }
  printf("---------------------------------------------------------------------"
         "----------------------------\n");
  printf(" Summary: \n   GEM modules: %i\n      Sectors: %i, Pads: %i\n   "
         "Straw modules: %i\n",
         nTotGems, nTotSectors, nTotChannels, nTotStraws);
  printf("====================================================================="
         "============================\n");
}
 
void CbmMuchSegmentAuto::DrawSegmentation()
{
  string digifile(fDigiFileName.Data());
  Int_t startIndex = digifile.size() - 4;
  string txtfile   = digifile.erase(startIndex, 4);
  txtfile.append("txt");
  FILE* outfile;
  outfile = fopen(txtfile.c_str(), "w");
 
  //  Int_t colors[] = {kGreen, kBlue, kViolet, kRed, kYellow, kOrange, kMagenta, kCyan,  kSpring, kPink, kAzure, kTeal,
  //      kGreen+10, kBlue+10, kViolet+10, kRed+10, kYellow+10, kOrange+10, kMagenta+10, kCyan+10,  kSpring+10, kPink+10, kAzure+10, kTeal+10};
 
  Double_t secMinLx = std::numeric_limits<Double_t>::max();
  Double_t secMinLy = std::numeric_limits<Double_t>::max();
  for (Int_t iStation = 0; iStation < fNStations; ++iStation) {
    CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
    CbmMuchLayer* layer     = station->GetLayer(0);
    for (Int_t iSide = 1; iSide >= 0; iSide--) {
      CbmMuchLayerSide* side = layer->GetSide(iSide);
      for (Int_t iModule = 0; iModule < side->GetNModules(); ++iModule) {
        CbmMuchModule* mod = side->GetModule(iModule);
        if (mod->GetDetectorType() != 1) continue;
        CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
        for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
          CbmMuchSectorRectangular* sector = (CbmMuchSectorRectangular*) module->GetSector(iSector);
          if (sector->GetSize()[0] < secMinLx) secMinLx = sector->GetSize()[0];
          if (sector->GetSize()[1] < secMinLy) secMinLy = sector->GetSize()[1];
        }
      }  // modules
    }    // sides
  }
 
  for (Int_t iStation = 0; iStation < fStations->GetEntriesFast(); ++iStation) {
    fprintf(outfile, "=================================================================="
                     "==========\n");
    fprintf(outfile, "Station %i\n", iStation + 1);
    fprintf(outfile, "Sector size, cm   Sector position, cm   Number of pads   Side   "
                     "Pad size, cm\n");
    fprintf(outfile, "------------------------------------------------------------------"
                     "----------\n");
    TCanvas* c1 = new TCanvas(Form("station%i", iStation + 1), Form("station%i", iStation + 1), 800, 800);
    c1->SetFillColor(0);
    c1->Range(-250, -250, 250, 250);
    CbmMuchStation* station = (CbmMuchStation*) fStations->At(iStation);
    CbmMuchLayer* layer     = station->GetLayer(0);
    for (Int_t iSide = 1; iSide >= 0; iSide--) {
      CbmMuchLayerSide* layerSide = layer->GetSide(iSide);
      for (Int_t iModule = 0; iModule < layerSide->GetNModules(); ++iModule) {
        CbmMuchModule* mod = layerSide->GetModule(iModule);
        mod->SetFillStyle(0);
        mod->Draw();
        if (mod->GetDetectorType() != 1) continue;
        CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
        for (Int_t iSector = 0; iSector < module->GetNSectors(); ++iSector) {
          CbmMuchSectorRectangular* sector = (CbmMuchSectorRectangular*) module->GetSector(iSector);
 
          //          Int_t i = Int_t((sector->GetSize()[0]+1e-3)/secMinLx) - 1;
          //          Int_t j = Int_t((sector->GetSize()[1]+1e-3)/secMinLy) - 1;
          // TODO
          //          sector->SetFillColor(iSide ? TColor::GetColorDark(colors[i+j]) : colors[i+j]);
          //          sector->Draw("f");
          //          sector->Draw();
          const char* side = iSide ? "Back" : "Front";
          fprintf(outfile, "%-4.2fx%-10.2f   %-6.2fx%-12.2f   %-14i   %-5s   ", sector->GetSize()[0],
                  sector->GetSize()[1], sector->GetPosition()[0], sector->GetPosition()[1], sector->GetNChannels(),
                  side);
          //TODO          fprintf(outfile, "%-4.2fx%-4.2f\n", sector->GetDx(), sector->GetDy());
        }  // sectors
      }    // modules
    }      // sides
 
    // Draw a hole
    TArc* arc = new TArc(0., 0., station->GetRmin());
    arc->Draw();
 
    c1->Print(Form("%s/station%i.eps", gSystem->DirName(fDigiFileName), iStation + 1));
    c1->Print(Form("%s/station%i.png", gSystem->DirName(fDigiFileName), iStation + 1));
  }  //stations
  fclose(outfile);
}
 
ClassImp(CbmMuchSegmentAuto)