CaloGeometry.cxx

/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

#include "ISF_FastCaloSimParametrization/CaloGeometry.h"
#include <TTree.h>
#include <TVector2.h>
#include <TRandom.h>
#include <TCanvas.h>
#include <TH2D.h>
#include <TGraphErrors.h>
#include <TVector3.h>
#include <TLegend.h>
#include <fstream>
#include <sstream>


#include "CaloDetDescr/CaloDetDescrElement.h"
//#include "ISF_FastCaloSimParametrization/CaloDetDescrElement.h"
#include "CaloGeoHelpers/CaloSampling.h"
#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
//#include "TMVA/Tools.h"
//#include "TMVA/Factory.h"

using namespace std;

int CaloGeometry_calocol[24]={1,2,3,4, // LAr barrel
                     1,2,3,4, // LAr EM endcap
                     1,2,3,4, // Hadronic end cap cal.
                     1,2,3,   // Tile barrel
                     -42,-28,6, // Tile gap (ITC & scint)
                     1,2,3,   // Tile extended barrel
                     1,2,3    // Forward EM endcap
                    }; 

const int CaloGeometry::MAX_SAMPLING = CaloCell_ID_FCS::MaxSample; //number of calorimeter layers/samplings

Identifier CaloGeometry::m_debug_identify;
bool CaloGeometry::m_debug=false;

CaloGeometry::CaloGeometry() : m_cells_in_sampling(MAX_SAMPLING),m_cells_in_sampling_for_phi0(MAX_SAMPLING),m_cells_in_regions(MAX_SAMPLING),m_isCaloBarrel(MAX_SAMPLING),m_dographs(false),m_FCal_ChannelMap(0)
{
  //TMVA::Tools::Instance();
  for(int i=0;i<2;++i) {
    m_min_eta_sample[i].resize(MAX_SAMPLING); //[side][calosample]
    m_max_eta_sample[i].resize(MAX_SAMPLING); //[side][calosample]
    m_rmid_map[i].resize(MAX_SAMPLING); //[side][calosample]
    m_zmid_map[i].resize(MAX_SAMPLING); //[side][calosample]
    m_rent_map[i].resize(MAX_SAMPLING); //[side][calosample]
    m_zent_map[i].resize(MAX_SAMPLING); //[side][calosample]
    m_rext_map[i].resize(MAX_SAMPLING); //[side][calosample]
    m_zext_map[i].resize(MAX_SAMPLING); //[side][calosample]
  }
  m_graph_layers.resize(MAX_SAMPLING);
//  for(unsigned int i=CaloCell_ID_FCS::FirstSample;i<CaloCell_ID_FCS::MaxSample;++i) {
  for(unsigned int i=CaloSampling::PreSamplerB;i<=CaloSampling::FCAL2;++i) {
    m_graph_layers[i]=0;
    CaloSampling::CaloSample s=static_cast<CaloSampling::CaloSample>(i);
    m_isCaloBarrel[i]=(CaloSampling::barrelPattern() & CaloSampling::getSamplingPattern(s))!=0;
  }
  m_isCaloBarrel[CaloCell_ID_FCS::TileGap3]=false; 
}

CaloGeometry::~CaloGeometry()
{
}

void CaloGeometry::addcell(const CaloDetDescrElement* cell) 
{
  int sampling=cell->getSampling();
  Identifier identify=cell->identify();
  
  m_cells[identify]=cell;
  m_cells_in_sampling[sampling][identify]=cell;
  
  //m_cells[identify]= new CaloDetDescrElement(*cell);
  //m_cells_in_sampling[sampling][identify]= new CaloDetDescrElement(*cell);
  
  CaloGeometryLookup* lookup=0;
  for(unsigned int i=0;i<m_cells_in_regions[sampling].size();++i) {
    if(m_cells_in_regions[sampling][i]->IsCompatible(cell)) {
      lookup=m_cells_in_regions[sampling][i];
      //cout<<sampling<<":  found compatible map"<<endl;
      break;
    }
  }
  if(!lookup) {
    lookup=new CaloGeometryLookup(m_cells_in_regions[sampling].size());
    m_cells_in_regions[sampling].push_back(lookup);
  }
  lookup->add(cell);
}

void CaloGeometry::PrintMapInfo(int i, int j) 
{
  cout<<"  map "<<j<<": "<<m_cells_in_regions[i][j]->size()<<" cells";
  if(i<21) {
    cout<<", "<<m_cells_in_regions[i][j]->cell_grid_eta()<<"*"<<m_cells_in_regions[i][j]->cell_grid_phi();
    cout<<", deta="<<m_cells_in_regions[i][j]->deta().mean()<<"+-"<<m_cells_in_regions[i][j]->deta().rms();
    cout<<", dphi="<<m_cells_in_regions[i][j]->dphi().mean()<<"+-"<<m_cells_in_regions[i][j]->dphi().rms();
    cout<<", mineta="<<m_cells_in_regions[i][j]->mineta()<<", maxeta="<<m_cells_in_regions[i][j]->maxeta();
    cout<<", minphi="<<m_cells_in_regions[i][j]->minphi()<<", maxphi="<<m_cells_in_regions[i][j]->maxphi();
    cout<<endl<<"         ";  
    cout<<", mineta_raw="<<m_cells_in_regions[i][j]->mineta_raw()<<", maxeta_raw="<<m_cells_in_regions[i][j]->maxeta_raw();
    cout<<", minphi_raw="<<m_cells_in_regions[i][j]->minphi_raw()<<", maxphi_raw="<<m_cells_in_regions[i][j]->maxphi_raw();
    cout<<endl;
  } else {
    cout<<", "<<m_cells_in_regions[i][j]->cell_grid_eta()<<"*"<<m_cells_in_regions[i][j]->cell_grid_phi();
    cout<<", dx="<<m_cells_in_regions[i][j]->dx().mean()<<"+-"<<m_cells_in_regions[i][j]->dx().rms();
    cout<<", dy="<<m_cells_in_regions[i][j]->dy().mean()<<"+-"<<m_cells_in_regions[i][j]->dy().rms();
    cout<<", mindx="<<m_cells_in_regions[i][j]->mindx();
    cout<<", mindy="<<m_cells_in_regions[i][j]->mindy();
    cout<<", minx="<<m_cells_in_regions[i][j]->minx()<<", maxx="<<m_cells_in_regions[i][j]->maxx();
    cout<<", miny="<<m_cells_in_regions[i][j]->miny()<<", maxy="<<m_cells_in_regions[i][j]->maxy();
    cout<<endl<<"         ";  
    cout<<", minx_raw="<<m_cells_in_regions[i][j]->minx_raw()<<", maxx_raw="<<m_cells_in_regions[i][j]->maxx_raw();
    cout<<", miny_raw="<<m_cells_in_regions[i][j]->miny_raw()<<", maxy_raw="<<m_cells_in_regions[i][j]->maxy_raw();
    cout<<endl;
  }
}

void CaloGeometry::post_process(int sampling)
{
  //cout<<"post processing sampling "<<sampling<<endl;
  bool found_overlap=false;
  for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
    /*
    cout<<"Sample "<<sampling<<": checking map "<<j<<"/"<<m_cells_in_regions[sampling].size();
    for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
      cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
    }
    cout<<endl;
    */
    for(unsigned int i=j+1;i<m_cells_in_regions[sampling].size();++i) {
      if(m_cells_in_regions[sampling][i]->has_overlap(m_cells_in_regions[sampling][j])) {
        if(!found_overlap) {
          cout<<"Sample "<<sampling<<", starting maps : "<<m_cells_in_regions[sampling].size();
          for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
            cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
          }
          cout<<endl;
        }
        found_overlap=true;
        /*
        cout<<"Sample "<<sampling<<": Found overlap between map "<<j<<" and "<<i<<", "
            <<m_cells_in_regions[sampling].size()<<" total maps"<<endl;
        cout<<"  current configuration map "<<j<<"/"<<m_cells_in_regions[sampling].size();
        for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
          cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
        }
        cout<<endl;

        PrintMapInfo(sampling,j);
        PrintMapInfo(sampling,i);
        */
        
        CaloGeometryLookup* toremove=m_cells_in_regions[sampling][i];
        toremove->merge_into_ref(m_cells_in_regions[sampling][j]);
        
        /*
        cout<<"  New ";
        PrintMapInfo(sampling,j);
        */
        
        for(unsigned int k=i;k<m_cells_in_regions[sampling].size()-1;++k) {
          m_cells_in_regions[sampling][k]=m_cells_in_regions[sampling][k+1];
          m_cells_in_regions[sampling][k]->set_index(k);
        }
        m_cells_in_regions[sampling].resize(m_cells_in_regions[sampling].size()-1);

        /*
        cout<<"  new configuration map "<<j<<"/"<<m_cells_in_regions[sampling].size();
        for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
          cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
        }
        cout<<endl;
        */
        
        --i;
      }
    }
  }
  
  if(found_overlap) {
    cout<<"Sample "<<sampling<<", final maps : "<<m_cells_in_regions[sampling].size();
    for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
      cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
    }
    cout<<endl;
  }

  if(found_overlap) {
    cout<<"Run another round of ";
    post_process(sampling);
  }  
}

void CaloGeometry::InitRZmaps()
{

  int nok=0;
  
  FSmap< double , double > rz_map_eta [2][MAX_SAMPLING];
  FSmap< double , double > rz_map_rmid[2][MAX_SAMPLING];
  FSmap< double , double > rz_map_zmid[2][MAX_SAMPLING];
  FSmap< double , double > rz_map_rent[2][MAX_SAMPLING];
  FSmap< double , double > rz_map_zent[2][MAX_SAMPLING];
  FSmap< double , double > rz_map_rext[2][MAX_SAMPLING];
  FSmap< double , double > rz_map_zext[2][MAX_SAMPLING];
  FSmap< double , int    > rz_map_n   [2][MAX_SAMPLING];


  for(unsigned int side=0;side<=1;++side) for(unsigned int sample=0;sample<MAX_SAMPLING;++sample)
  {
    m_min_eta_sample[side][sample]=+1000;
    m_max_eta_sample[side][sample]=-1000;
  }  
  
  
  for(t_cellmap::iterator calo_iter=m_cells.begin();calo_iter!=m_cells.end();++calo_iter)
  {
    const CaloDetDescrElement* theDDE=(*calo_iter).second;
    if(theDDE)
    {
      ++nok;
      unsigned int sample=theDDE->getSampling();

      int side=0;
      int sign_side=-1;
      double eta_raw=theDDE->eta_raw();
      if(eta_raw>0) {
        side=1;
        sign_side=+1;
      }
      
      if(!m_cells_in_sampling_for_phi0[sample][eta_raw]) {
        m_cells_in_sampling_for_phi0[sample][eta_raw]=theDDE;
      } else {
        if(TMath::Abs(theDDE->phi()) < TMath::Abs(m_cells_in_sampling_for_phi0[sample][eta_raw]->phi())) {
          m_cells_in_sampling_for_phi0[sample][eta_raw]=theDDE;
        }
      }  
      
      double min_eta=theDDE->eta()-theDDE->deta()/2;
      double max_eta=theDDE->eta()+theDDE->deta()/2;
      if(min_eta<m_min_eta_sample[side][sample]) m_min_eta_sample[side][sample]=min_eta;
      if(max_eta>m_max_eta_sample[side][sample]) m_max_eta_sample[side][sample]=max_eta;
      
      if(rz_map_eta[side][sample].find(eta_raw)==rz_map_eta[side][sample].end()) {
        rz_map_eta [side][sample][eta_raw]=0;
        rz_map_rmid[side][sample][eta_raw]=0;
        rz_map_zmid[side][sample][eta_raw]=0;
        rz_map_rent[side][sample][eta_raw]=0;
        rz_map_zent[side][sample][eta_raw]=0;
        rz_map_rext[side][sample][eta_raw]=0;
        rz_map_zext[side][sample][eta_raw]=0;
        rz_map_n   [side][sample][eta_raw]=0;
      }
      rz_map_eta [side][sample][eta_raw]+=theDDE->eta();
      rz_map_rmid[side][sample][eta_raw]+=theDDE->r();
      rz_map_zmid[side][sample][eta_raw]+=theDDE->z();
      double drh=theDDE->dr()/2;
      double dzh=theDDE->dz();
      if(sample<=CaloSampling::EMB3) { // ensure we have a valid sampling
        drh=theDDE->dr();
      } 
      // An `else` would be good here since we can't continue with a *bad* sampling...
      // Should most likely handle the situation nicely rather than with a crash.

      rz_map_rent[side][sample][eta_raw]+=theDDE->r()-drh;
      rz_map_zent[side][sample][eta_raw]+=theDDE->z()-dzh*sign_side;
      rz_map_rext[side][sample][eta_raw]+=theDDE->r()+drh;
      rz_map_zext[side][sample][eta_raw]+=theDDE->z()+dzh*sign_side;
      rz_map_n   [side][sample][eta_raw]++;
      
    }
  }
  
  for(int side=0;side<=1;++side)
  {
   for(int sample=0;sample<MAX_SAMPLING;++sample)
   {

    if(rz_map_n[side][sample].size()>0)
    {
      for(FSmap< double , int >::iterator iter=rz_map_n[side][sample].begin();iter!=rz_map_n[side][sample].end();++iter)
      {
        double eta_raw=iter->first;
        if(iter->second<1)
        {
          //ATH_MSG_WARNING("rz-map for side="<<side<<" sample="<<sample<<" eta_raw="<<eta_raw<<" : #cells="<<iter->second<<" !!!");
        }
        else
        {
          double eta =rz_map_eta[side][sample][eta_raw]/iter->second;
          double rmid=rz_map_rmid[side][sample][eta_raw]/iter->second;
          double zmid=rz_map_zmid[side][sample][eta_raw]/iter->second;
          double rent=rz_map_rent[side][sample][eta_raw]/iter->second;
          double zent=rz_map_zent[side][sample][eta_raw]/iter->second;
          double rext=rz_map_rext[side][sample][eta_raw]/iter->second;
          double zext=rz_map_zext[side][sample][eta_raw]/iter->second;
          
          m_rmid_map[side][sample][eta]=rmid;
          m_zmid_map[side][sample][eta]=zmid;
          m_rent_map[side][sample][eta]=rent;
          m_zent_map[side][sample][eta]=zent;
          m_rext_map[side][sample][eta]=rext;
          m_zext_map[side][sample][eta]=zext;
        }
      }
      //ATH_MSG_DEBUG("rz-map for side="<<side<<" sample="<<sample<<" #etas="<<m_rmid_map[side][sample].size());
     }
    else
    {
     std::cout<<"rz-map for side="<<side<<" sample="<<sample<<" is empty!!!"<<std::endl;
    }
   } //sample
  } //side
  
  if(DoGraphs()) {
    for(int sample=0;sample<MAX_SAMPLING;++sample) {
      m_graph_layers[sample]=new TGraphErrors(rz_map_n[0][sample].size()+rz_map_n[1][sample].size());
      m_graph_layers[sample]->SetMarkerColor(TMath::Abs(CaloGeometry_calocol[sample]));
      m_graph_layers[sample]->SetLineColor(TMath::Abs(CaloGeometry_calocol[sample]));
      int np=0;
      for(int side=0;side<=1;++side) {
        for(FSmap< double , int >::iterator iter=rz_map_n[side][sample].begin();iter!=rz_map_n[side][sample].end();++iter) {
          double eta_raw=iter->first;
          int sign_side=-1;
          if(eta_raw>0) sign_side=+1;
          //double eta =rz_map_eta[side][sample][eta_raw]/iter->second;
          double rmid=rz_map_rmid[side][sample][eta_raw]/iter->second;
          double zmid=rz_map_zmid[side][sample][eta_raw]/iter->second;
          //double rent=rz_map_rent[side][sample][eta_raw]/iter->second;
          //double zent=rz_map_zent[side][sample][eta_raw]/iter->second;
          double rext=rz_map_rext[side][sample][eta_raw]/iter->second;
          double zext=rz_map_zext[side][sample][eta_raw]/iter->second;
          m_graph_layers[sample]->SetPoint(np,zmid,rmid);
          /*
          if(isCaloBarrel(sample)) {
            m_graph_layers[sample]->SetPointError(np,0,rext-rmid);
          } else {
            m_graph_layers[sample]->SetPointError(np,(zext-zent)*sign_side,0);
          }
          */
          m_graph_layers[sample]->SetPointError(np,(zext-zmid)*sign_side,rext-rmid);
          ++np;
        }
      }  
    }
  }
}

TGraph* CaloGeometry::DrawGeoSampleForPhi0(int sample, int calocol, bool print)
{
  TGraph* firstgr=0;
  cout<<"Start sample "<<sample<<" ("<<SamplingName(sample)<<")"<<endl;
  int ngr=0;
  for(t_eta_cellmap::iterator calo_iter=m_cells_in_sampling_for_phi0[sample].begin();calo_iter!=m_cells_in_sampling_for_phi0[sample].end();++calo_iter) {
    const CaloDetDescrElement* theDDE=(*calo_iter).second;
    if(theDDE) {
      TVector3 cv;
      TGraph* gr=new TGraph(5);
      gr->SetLineColor(TMath::Abs(calocol));
      gr->SetFillColor(TMath::Abs(calocol));
      if(calocol<0) {
        gr->SetFillStyle(1001);
      } else {
        gr->SetFillStyle(0);
      }
      gr->SetLineWidth(2);
      double r=theDDE->r();
      double dr=theDDE->dr();
      double x=theDDE->x();
      double dx=theDDE->dx();
      double y=theDDE->y();
      double dy=theDDE->dy();
      double z=theDDE->z();
      double dz=theDDE->dz()*2;
      double eta=theDDE->eta();
      double deta=theDDE->deta();
      
      if(CaloSampling::PreSamplerB<=sample && sample<=CaloSampling::EMB3) {
       dr*=2;
      }
      if(print) {
        cout<<"sample="<<sample<<" r="<<r<<" dr="<<dr<<" eta="<<eta<<" deta="<<deta<<" x="<<x<<" y="<<y<<" z="<<z<<" dz="<<dz<<endl;
      }
      if(isCaloBarrel(sample)) {
        cv.SetPtEtaPhi(r-dr/2,eta-deta/2,0);
        gr->SetPoint(0,cv.Z(),cv.Pt());
        gr->SetPoint(4,cv.Z(),cv.Pt());
        cv.SetPtEtaPhi(r-dr/2,eta+deta/2,0);
        gr->SetPoint(1,cv.Z(),cv.Pt());
        cv.SetPtEtaPhi(r+dr/2,eta+deta/2,0);
        gr->SetPoint(2,cv.Z(),cv.Pt());
        cv.SetPtEtaPhi(r+dr/2,eta-deta/2,0);
        gr->SetPoint(3,cv.Z(),cv.Pt());
      } else {
        if(sample<CaloSampling::FCAL0) {
          cv.SetPtEtaPhi(1,eta-deta/2,0);cv*=(z-dz/2)/cv.Z();
          gr->SetPoint(0,cv.Z(),cv.Pt());
          gr->SetPoint(4,cv.Z(),cv.Pt());
          cv.SetPtEtaPhi(1,eta+deta/2,0);cv*=(z-dz/2)/cv.Z();
          gr->SetPoint(1,cv.Z(),cv.Pt());
          cv.SetPtEtaPhi(1,eta+deta/2,0);cv*=(z+dz/2)/cv.Z();
          gr->SetPoint(2,cv.Z(),cv.Pt());
          cv.SetPtEtaPhi(1,eta-deta/2,0);cv*=(z+dz/2)/cv.Z();
          gr->SetPoint(3,cv.Z(),cv.Pt());
        } else {
          double minr=r;
          double maxr=r;
          for(double px=x-dx/2;px<=x+dx/2;px+=dx) {
            for(double py=y-dy/2;py<=y+dy/2;py+=dy) {
              double pr=TMath::Sqrt(px*px+py*py);
              minr=TMath::Min(minr,pr);
              maxr=TMath::Max(maxr,pr);
            }
          }
          cv.SetXYZ(minr,0,z-dz/2);
          gr->SetPoint(0,cv.Z(),cv.Pt());
          gr->SetPoint(4,cv.Z(),cv.Pt());
          cv.SetXYZ(maxr,0,z-dz/2);
          gr->SetPoint(1,cv.Z(),cv.Pt());
          cv.SetXYZ(maxr,0,z+dz/2);
          gr->SetPoint(2,cv.Z(),cv.Pt());
          cv.SetXYZ(minr,0,z+dz/2);
          gr->SetPoint(3,cv.Z(),cv.Pt());
        }  
      }
      //if(calocol[sample]>0) gr->Draw("Lsame");
      // else gr->Draw("LFsame");
      gr->Draw("LFsame");
      if(ngr==0) firstgr=gr;
      ++ngr;
    }
  }  
  cout<<"Done sample "<<sample<<" ("<<SamplingName(sample)<<")="<<ngr<<endl;
  return firstgr;
}

TCanvas* CaloGeometry::DrawGeoForPhi0()
{
  TCanvas* c=new TCanvas("CaloGeoForPhi0","Calo geometry for #phi~0");
  TH2D* hcalolayout=new TH2D("hcalolayoutPhi0","Reconstruction geometry: calorimeter layout;z [mm];r [mm]",50,-7000,7000,50,0,4000);
  hcalolayout->Draw();
  hcalolayout->SetStats(0);
  hcalolayout->GetYaxis()->SetTitleOffset(1.4);
  
  TLegend* leg=new TLegend(0.30,0.13,0.70,0.37);
  leg->SetFillStyle(0);
  leg->SetFillColor(10);
  leg->SetBorderSize(1);
  leg->SetNColumns(2);

  for(int sample=0;sample<MAX_SAMPLING;++sample) {
    TGraph* gr=DrawGeoSampleForPhi0(sample,CaloGeometry_calocol[sample],(sample==17));
    if(gr) {
      std::string sname=Form("Sampling %2d : ",sample);
      sname+=SamplingName(sample);
      leg->AddEntry(gr,sname.c_str(),"LF");
    }
  }
  leg->Draw();
  return c;
}

const CaloDetDescrElement* CaloGeometry::getDDE(Identifier identify) 
{
  return m_cells[identify];
}
const CaloDetDescrElement* CaloGeometry::getDDE(int sampling,Identifier identify) 
{
  return m_cells_in_sampling[sampling][identify];
}

const CaloDetDescrElement* CaloGeometry::getDDE(int sampling,float eta,float phi,float* distance,int* steps) 
{
  if(sampling<0) return 0;
  if(sampling>=MAX_SAMPLING) return 0;
  if(m_cells_in_regions[sampling].size()==0) return 0;
  
  float dist;
  const CaloDetDescrElement* bestDDE=0;
  if(!distance) distance=&dist;
  *distance=+10000000;
  int intsteps;
  int beststeps;
  if(steps) beststeps=(*steps);
   else beststeps=0;
  
  if(sampling<21) {
    for(int skip_range_check=0;skip_range_check<=1;++skip_range_check) {
      for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
        if(!skip_range_check) {
          if(eta<m_cells_in_regions[sampling][j]->mineta()) continue;
          if(eta>m_cells_in_regions[sampling][j]->maxeta()) continue;
        }  
        if(steps) intsteps=(*steps);
         else intsteps=0;
        if(m_debug) {
          cout<<"CaloGeometry::getDDE : check map"<<j<<" skip_range_check="<<skip_range_check<<endl;
        }
        float newdist;
        const CaloDetDescrElement* newDDE=m_cells_in_regions[sampling][j]->getDDE(eta,phi,&newdist,&intsteps);
        if(m_debug) {
          cout<<"CaloGeometry::getDDE : map"<<j<<" dist="<<newdist<<" best dist="<<*distance<<" steps="<<intsteps<<endl;
        }
        if(newdist<*distance) {
          bestDDE=newDDE;
          *distance=newdist;
          if(steps) beststeps=intsteps;
          if(newdist<-0.1) break; //stop, we are well within the hit cell
        }
      }
      if(bestDDE) break;
    }
  } else {
		return 0;
    //for(int skip_range_check=0;skip_range_check<=1;++skip_range_check) {
      //for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
        //if(!skip_range_check) {
          //if(eta<m_cells_in_regions[sampling][j]->minx()) continue;
          //if(eta>m_cells_in_regions[sampling][j]->maxx()) continue;
          //if(phi<m_cells_in_regions[sampling][j]->miny()) continue;
          //if(phi>m_cells_in_regions[sampling][j]->maxy()) continue;
        //}  
        //if(steps) intsteps=*steps;
         //else intsteps=0;
        //if(m_debug) {
          //cout<<"CaloGeometry::getDDE : check map"<<j<<" skip_range_check="<<skip_range_check<<endl;
        //}
        //float newdist;
        //const CaloGeoDetDescrElement* newDDE=m_cells_in_regions[sampling][j]->getDDE(eta,phi,&newdist,&intsteps);
        //if(m_debug) {
          //cout<<"CaloGeometry::getDDE : map"<<j<<" dist="<<newdist<<" best dist="<<*distance<<" steps="<<intsteps<<endl;
        //}
        //if(newdist<*distance) {
          //bestDDE=newDDE;
          //*distance=newdist;
          //if(steps) beststeps=intsteps;
          //if(newdist<-0.1) break; //stop, we are well within the hit cell
        //}
      //}
      //if(bestDDE) break;
    //}  
  }
  if(steps) *steps=beststeps;
  return bestDDE;
}

const CaloDetDescrElement* CaloGeometry::getFCalDDE(int sampling,float x,float y,float z,float* distance,int* steps){
  int isam = sampling - 20;
  int iphi(-100000),ieta(-100000);
  Long64_t mask1[]{0x34,0x34,0x35};
  Long64_t mask2[]{0x36,0x36,0x37};
  bool found = m_FCal_ChannelMap.getTileID(isam, x, y, ieta, iphi);
  if(steps && found) *steps=0;
  if(!found) {
    //cout << "Warning: Hit is not matched with any FCal cell! Looking for the closest cell" << endl;
    found = getClosestFCalCellIndex(sampling, x, y, ieta, iphi,steps);
  }
  if(!found) {
    cout << "Error: Unable to find the closest FCal cell!" << endl;
    return nullptr;
  }
  
  
  //cout << "CaloGeometry::getFCalDDE: x:" << x << " y: " << y << " ieta: " << ieta << " iphi: " << iphi << " cmap->x(): " << m_FCal_ChannelMap.x(isam,ieta,iphi) << " cmap->y(): " << m_FCal_ChannelMap.y(isam,ieta,iphi) << endl;
  Long64_t id = (ieta << 5) + 2*iphi;
  if(isam==2)id+= (8<<8);
  
  
  
  if(z>0) id+=(mask2[isam-1] << 12);
  else id+=(mask1[isam-1] << 12);
  
  id = id << 44; 
  Identifier identify((unsigned long long)id);
  
  const CaloDetDescrElement* foundcell=m_cells[identify];
  if(distance) {
    *distance=sqrt(pow(foundcell->x() - x,2) +  pow(foundcell->y() - y,2)  );
  }
  
  return foundcell;
}


bool CaloGeometry::getClosestFCalCellIndex(int sampling,float x,float y,int& ieta, int& iphi,int* steps){
  
  double rmin = m_FCal_rmin[sampling-21];
  double rmax = m_FCal_rmax[sampling-21];
  int isam=sampling-20;
  double a=1.;
  const double b=0.01;
  const int nmax=100;
  int i=0;
  
  const double r = sqrt(x*x +y*y);
  if(r==0.) return false;
  const double r_inverse=1./r;
  
  if((r/rmax)>(rmin*r_inverse)){
    x=x*rmax*r_inverse;
    y=y*rmax*r_inverse;
    while((!m_FCal_ChannelMap.getTileID(isam, a*x, a*y, ieta, iphi)) && i<nmax){
      a-=b;
      i++;
    }
  }
  else {
    x=x*rmin*r_inverse;
    y=y*rmin*r_inverse;
    while((!m_FCal_ChannelMap.getTileID(isam, a*x, a*y, ieta, iphi)) && i<nmax){
      a+=b;
      i++;
    }
    
  }
  if(steps)*steps=i+1;
  return i<nmax ? true : false;
}

bool CaloGeometry::PostProcessGeometry()
{
  for(int i=0;i<MAX_SAMPLING;++i) {
    post_process(i);
    for(unsigned int j=0;j<m_cells_in_regions[i].size();++j) {
      m_cells_in_regions[i][j]->post_process();
    }
    //if(i>=21) break;
  }
  
  InitRZmaps(); 
  
  /*
  cout<<"all : "<<m_cells.size()<<endl;
  for(int sampling=0;sampling<MAX_SAMPLING;++sampling) {
    cout<<"cells sampling "<<sampling<<" : "<<m_cells_in_sampling[sampling].size()<<" cells";
    cout<<", "<<m_cells_in_regions[sampling].size()<<" lookup map(s)"<<endl;
    for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
      PrintMapInfo(sampling,j);
      //break;
    }
    //if(i>0) break;
  } 
  */ 

  return true;
}

void CaloGeometry::Validate(int nrnd)
{
  int ntest=0;
  cout<<"start CaloGeometry::Validate()"<<endl;
  for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
    const CaloDetDescrElement* cell=ic->second;
    int sampling=cell->getSampling();
    //if(sampling>=21) continue;

    if(m_debug_identify==cell->identify()) {
      cout<<"CaloGeometry::Validate(), cell "<<ntest<<" id="<<cell->identify()<<endl; 
      m_debug=true;
    }  
    
    for(int irnd=0;irnd<nrnd;++irnd) {
      std::stringstream pos;
      std::stringstream cellpos;
      std::stringstream foundcellpos;
      int steps=0;
      float distance=0;
      bool is_inside;
      bool is_inside_foundcell=false;
      const CaloDetDescrElement* foundcell=0;
      if(sampling<21) {
        float eta=cell->eta()+1.95*(gRandom->Rndm()-0.5)*cell->deta();
        float phi=cell->phi()+1.95*(gRandom->Rndm()-0.5)*cell->dphi();
        foundcell=getDDE(sampling,eta,phi,&distance,&steps);
        
        pos<<"eta="<<eta<<" phi="<<phi;
        cellpos<<"eta="<<cell->eta()<<" eta_raw="<<cell->eta_raw()<<" deta="<<cell->deta()
               <<" ("<<(cell->eta_raw()-cell->eta())/cell->deta()<<") ; "
               <<"phi="<<cell->phi()<<" phi_raw="<<cell->phi_raw()<<" dphi="<<cell->dphi()
               <<" ("<<(cell->phi_raw()-cell->phi())/cell->dphi()<<")";
        if(foundcell) {
          foundcellpos<<"eta="<<foundcell->eta()<<" eta_raw="<<foundcell->eta_raw()<<" deta="<<foundcell->deta()
                      <<" ("<<(foundcell->eta_raw()-foundcell->eta())/foundcell->deta()<<") ; "
                      <<"phi="<<foundcell->phi()<<" phi_raw="<<foundcell->phi_raw()<<" dphi="<<foundcell->dphi()
                      <<" ("<<(foundcell->phi_raw()-foundcell->phi())/cell->dphi()<<")";
          is_inside_foundcell=TMath::Abs( (eta-foundcell->eta())/foundcell->deta() )<0.55 && TMath::Abs( (phi-foundcell->phi())/foundcell->dphi() )<0.55;
        }
        is_inside=TMath::Abs( (eta-cell->eta())/cell->deta() )<0.49 && TMath::Abs( (phi-cell->phi())/cell->dphi() )<0.49;
      } else {
        float x=cell->x()+1.95*(gRandom->Rndm()-0.5)*cell->dx();
        float y=cell->y()+1.95*(gRandom->Rndm()-0.5)*cell->dy();
        float z=cell->z();
        foundcell=getFCalDDE(sampling,x,y,z);
        
        pos<<"x="<<x<<" y="<<y<<" z="<<z;
        cellpos<<"x="<<cell->x()<<" x_raw="<<cell->x_raw()<<" dx="<<cell->dx()
               <<" ("<<(cell->x_raw()-cell->x())/cell->dx()<<") ; "
               <<"y="<<cell->y()<<" y_raw="<<cell->y_raw()<<" dy="<<cell->dy()
               <<" ("<<(cell->y_raw()-cell->y())/cell->dy()<<")";
        if(foundcell) {
          foundcellpos<<"x="<<foundcell->x()<<" x_raw="<<foundcell->x_raw()<<" dx="<<foundcell->dx()
                      <<" ("<<(foundcell->x_raw()-foundcell->x())/foundcell->dx()<<") ; "
                      <<"y="<<foundcell->y()<<" y_raw="<<foundcell->y_raw()<<" dy="<<foundcell->dy()
                      <<" ("<<(foundcell->y_raw()-foundcell->y())/cell->dy()<<")";
          is_inside_foundcell=TMath::Abs( (x-foundcell->x())/foundcell->dx() )<0.75 && TMath::Abs( (y-foundcell->y())/foundcell->dy() )<0.75;
        }       
        is_inside=TMath::Abs( (x-cell->x())/cell->dx() )<0.49 && TMath::Abs( (y-cell->y())/cell->dy() )<0.49;
        //m_debug=true;
      }  

      if(m_debug && foundcell) {
        cout<<"CaloGeometry::Validate(), irnd="<<irnd<<": cell id="<<cell->identify()<<", sampling="<<sampling
                                                     <<", foundcell id="<<foundcell->identify()<<", "<<steps<<" steps"<<endl; 
        cout<<"  "<<cellpos.str()<<endl;
      }  
      if(cell==foundcell) {
        if(steps>3 && distance<-0.01) { 
          cout<<"cell id="<<cell->identify()<<", sampling="<<sampling<<" found in "<<steps<<" steps, dist="<<distance<<" "<<pos.str()<<endl;
          cout<<"  "<<cellpos.str()<<endl;
        }
      } else {
        if(!foundcell) {
          cout<<"cell id="<<cell->identify()<<" not found!";
          cout<<" No cell found in "<<steps<<" steps, dist="<<distance<<" "<<pos.str()<<endl;
          cout<<"  input sampling="<<sampling<<" "<<cellpos.str()<<endl;
          //return;
        }
        if( is_inside && foundcell && !is_inside_foundcell) {
          cout<<"cell id="<<cell->identify()<<" not found, but inside cell area!";
          cout<<" Found instead id="<<foundcell->identify()<<" in "<<steps<<" steps, dist="<<distance<<" "<<pos.str()<<endl;
          cout<<"  input sampling="<<sampling<<" "<<cellpos.str()<<endl;
          cout<<"  output sampling="<<foundcell->getSampling()<<" "<<foundcellpos.str()<<endl;
          return;
        }  
      }
    }  
    m_debug=false;
    if(ntest%25000==0) cout<<"Validate cell "<<ntest<<" with "<<nrnd<<" random hits"<<endl;
    ++ntest;
  }
  cout<<"end CaloGeometry::Validate()"<<endl;
}

double CaloGeometry::deta(int sample,double eta) const
{
  int side=0;
  if(eta>0) side=1;

  double mineta=m_min_eta_sample[side][sample];
  double maxeta=m_max_eta_sample[side][sample];

  if(eta<mineta)
  {
    return fabs(eta-mineta);
  }
  else if(eta>maxeta)
  {
   return fabs(eta-maxeta);
	}
	else
	{
   double d1=fabs(eta-mineta);
   double d2=fabs(eta-maxeta);
   if(d1<d2) return -d1;
   else return -d2;
  }
}


void CaloGeometry::minmaxeta(int sample,double eta,double& mineta,double& maxeta) const 
{
  int side=0;
  if(eta>0) side=1;
  
  mineta=m_min_eta_sample[side][sample];
  maxeta=m_max_eta_sample[side][sample];
}

double CaloGeometry::rmid(int sample,double eta) const 
{
  int side=0;
  if(eta>0) side=1;

  return m_rmid_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::zmid(int sample,double eta) const 
{
  int side=0;
  if(eta>0) side=1;

  return m_zmid_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rzmid(int sample,double eta) const
{
 int side=0;
 if(eta>0) side=1;
	
 if(isCaloBarrel(sample)) return m_rmid_map[side][sample].find_closest(eta)->second;
 else                     return m_zmid_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rent(int sample,double eta) const 
{
  int side=0;
  if(eta>0) side=1;

  return m_rent_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::zent(int sample,double eta) const 
{
  int side=0;
  if(eta>0) side=1;

  return m_zent_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rzent(int sample,double eta) const
{
  int side=0;
  if(eta>0) side=1;

  if(isCaloBarrel(sample)) return m_rent_map[side][sample].find_closest(eta)->second;
   else                    return m_zent_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rext(int sample,double eta) const 
{
  int side=0;
  if(eta>0) side=1;

  return m_rext_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::zext(int sample,double eta) const 
{
  int side=0;
  if(eta>0) side=1;

  return m_zext_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rzext(int sample,double eta) const
{
  int side=0;
  if(eta>0) side=1;

  if(isCaloBarrel(sample)) return m_rext_map[side][sample].find_closest(eta)->second;
   else                    return m_zext_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rpos(int sample,double eta,int subpos) const
{
	
  int side=0;
  if(eta>0) side=1;
   
  if(subpos==CaloSubPos::SUBPOS_ENT) return m_rent_map[side][sample].find_closest(eta)->second;
  if(subpos==CaloSubPos::SUBPOS_EXT) return m_rext_map[side][sample].find_closest(eta)->second;
    
  return m_rmid_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::zpos(int sample,double eta,int subpos) const
{
  int side=0;
  if(eta>0) side=1;

  if(subpos==CaloSubPos::SUBPOS_ENT) return m_zent_map[side][sample].find_closest(eta)->second;
  if(subpos==CaloSubPos::SUBPOS_EXT) return m_zext_map[side][sample].find_closest(eta)->second;
  return m_zmid_map[side][sample].find_closest(eta)->second;
}

double CaloGeometry::rzpos(int sample,double eta,int subpos) const
{
  int side=0;
  if(eta>0) side=1;
 
  if(isCaloBarrel(sample)) {
    if(subpos==CaloSubPos::SUBPOS_ENT) return m_rent_map[side][sample].find_closest(eta)->second;
    if(subpos==CaloSubPos::SUBPOS_EXT) return m_rext_map[side][sample].find_closest(eta)->second;
    return m_rmid_map[side][sample].find_closest(eta)->second;
  } else {
    if(subpos==CaloSubPos::SUBPOS_ENT) return m_zent_map[side][sample].find_closest(eta)->second;
    if(subpos==CaloSubPos::SUBPOS_EXT) return m_zext_map[side][sample].find_closest(eta)->second;
    return m_zmid_map[side][sample].find_closest(eta)->second;
  }  
}

std::string CaloGeometry::SamplingName(int sample)
{
  return CaloSampling::getSamplingName(sample);
}

void  CaloGeometry::calculateFCalRminRmax(){
   
   m_FCal_rmin.resize(3,FLT_MAX);
   m_FCal_rmax.resize(3,0.);

   double x(0.),y(0.),r(0.);
   for(int imap=1;imap<=3;imap++)for(auto it=m_FCal_ChannelMap.begin(imap);it!=m_FCal_ChannelMap.end(imap);it++){
      x=it->second.x();
      y=it->second.y();
      r=sqrt(x*x+y*y);
      if(r<m_FCal_rmin[imap-1])m_FCal_rmin[imap-1]=r;
      if(r>m_FCal_rmax[imap-1])m_FCal_rmax[imap-1]=r;
   }
   
}


bool CaloGeometry::checkFCalGeometryConsistency(){
  
  unsigned long long phi_index,eta_index;
  float x,y,dx,dy;
  long id;
  
  long mask1[]{0x34,0x34,0x35};
  long mask2[]{0x36,0x36,0x37};
  
  m_FCal_rmin.resize(3,FLT_MAX);
  m_FCal_rmax.resize(3,0.);
  
  
  for(int imap=1;imap<=3;imap++){
    
    int sampling = imap+20;
    
    if((int)m_cells_in_sampling[sampling].size() != 2*std::distance(m_FCal_ChannelMap.begin(imap), m_FCal_ChannelMap.end(imap))){
      cout << "Error: Incompatibility between FCalChannel map and GEO file: Different number of cells in m_cells_in_sampling and FCal_ChannelMap" << endl;
      cout << "m_cells_in_sampling: " << m_cells_in_sampling[sampling].size() << endl;
      cout << "FCal_ChannelMap: " << 2*std::distance(m_FCal_ChannelMap.begin(imap), m_FCal_ChannelMap.end(imap)) << endl;
      return false;
    }

    for(auto it=m_FCal_ChannelMap.begin(imap);it!=m_FCal_ChannelMap.end(imap);it++){


      phi_index = it->first & 0xffff;
      eta_index = it->first >> 16;
      x=it->second.x();
      y=it->second.y();
      m_FCal_ChannelMap.tileSize(imap, eta_index, phi_index,dx,dy);

      id=(mask1[imap-1]<<12) + (eta_index << 5) +2*phi_index;

      if(imap==2) id+= (8<<8);

      Identifier id1((unsigned long long)(id<<44));
      const CaloDetDescrElement *DDE1 =getDDE(id1);

      id=(mask2[imap-1]<<12) + (eta_index << 5) +2*phi_index;
      if(imap==2) id+= (8<<8);
  		Identifier id2((unsigned long long)(id<<44));
      const CaloDetDescrElement *DDE2=getDDE(id2);

      if(!TMath::AreEqualRel(x, DDE1->x(),1.E-8) || !TMath::AreEqualRel(y, DDE1->y(),1.E-8) || !TMath::AreEqualRel(x, DDE2->x(),1.E-8) || !TMath::AreEqualRel(y, DDE2->y(),1.E-8) ){
	 cout << "Error: Incompatibility between FCalChannel map and GEO file \n" << x << " " << DDE1->x() << " " << DDE2->x() << y << " " << DDE1->y() << " " << DDE2->y() << endl;
	 return false;   
      }
    }



  }
   
  return true;
}