EcalHitMaker.cc

Go to the documentation of this file.

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
//#include "DataFormats/EcalDetId/interface/EBDetId.h"
//#include "DataFormats/EcalDetId/interface/EEDetId.h"

//#include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
//#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
//#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"

#include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
#include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
#include "FastSimulation/CaloGeometryTools/interface/CrystalWindowMap.h"
#include "FastSimulation/CaloGeometryTools/interface/CaloDirectionOperations.h"
#include "FastSimulation/Event/interface/FSimVertex.h"
#include "FastSimulation/Event/interface/FSimTrack.h"
#include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer1Properties.h"
#include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer2Properties.h"
#include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
#include "FastSimulation/Utilities/interface/RandomEngine.h"
//#include "FastSimulation/Utilities/interface/Histos.h"

//#include "Math/GenVector/Transform3D.h"
#include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h"

#include <algorithm>
#include <cmath>

typedef ROOT::Math::Plane3D::Vector Vector;
typedef ROOT::Math::Plane3D::Point Point;
typedef ROOT::Math::Transform3DPJ Transform3DR;

EcalHitMaker::EcalHitMaker(CaloGeometryHelper * theCalo,
               const XYZPoint& ecalentrance, 
               const DetId& cell, int onEcal,
               unsigned size, unsigned showertype,
               const RandomEngine* engine):
  CaloHitMaker(theCalo,DetId::Ecal,((onEcal==1)?EcalBarrel:EcalEndcap),onEcal,showertype),
  EcalEntrance_(ecalentrance),
  onEcal_(onEcal),
  myTrack_(NULL),
  random(engine)
{
#ifdef FAMOSDEBUG
  myHistos = Histos::instance();
#endif
  //  myHistos->debug("Constructeur EcalHitMaker");
  simulatePreshower_ = true;
  X0depthoffset_ = 0. ;
  X0PS1_ = 0.;
  X0PS2_ = 0.; 
  X0PS2EE_ = 0.;
  X0ECAL_ = 0.;
  X0EHGAP_ = 0.;
  X0HCAL_ = 0.;
  L0PS1_ = 0.;
  L0PS2_ = 0.;
  L0PS2EE_ = 0.;
  L0ECAL_ = 0.;
  L0EHGAP_ = 0.;
  L0HCAL_ = 0.;
  maxX0_ = 0.;
  totalX0_ = 0;
  totalL0_ = 0.;
  pulledPadProbability_ = 1.;
  outsideWindowEnergy_ = 0.;
  rearleakage_ = 0.;
  bfactor_ = 1.; 
  ncrystals_ = 0;

  doreorg_ = !showertype;

  hitmaphasbeencalculated_ = false;

  if(onEcal) 
    myCalorimeter->buildCrystal(cell,pivot_);
  else
    pivot_=Crystal();
  central_=onEcal==1;
  ecalFirstSegment_=-1;
  
  myCrystalWindowMap_ = 0; 
  // In some cases, a "dummy" grid, not based on a cell, can be built. The previous variables
  // should however be initialized. In such a case onEcal=0
  if(!onEcal) return;

  // Same size in eta-phi
  etasize_ = size;
  phisize_ = size;

  // Build the grid
  // The result is put in CellsWindow and is ordered by distance to the pivot
  myCalorimeter->getWindow(pivot_.getDetId(),size,size,CellsWindow_);

  buildGeometry();
  //  std::cout << " Geometry built " << regionOfInterest_.size() << std::endl;

  truncatedGrid_ = CellsWindow_.size()!=(etasize_*phisize_);

  // A local vector of corners
  mycorners.resize(4);
  corners.resize(4);
  
#ifdef DEBUGGW
  myHistos->fill("h10",EcalEntrance_.eta(),CellsWindow_.size());
  if(onEcal==2) 
    {
      myHistos->fill("h20",EcalEntrance_.perp(),CellsWindow_.size());
      if(EcalEntrance_.perp()>70&&EcalEntrance_.perp()<80&&CellsWindow_.size()<35)
    {
      std::cout << " Truncated grid " << CellsWindow_.size() << " " << EcalEntrance_.perp() << std::endl;
      std::cout << " Pivot " << myCalorimeter->getEcalEndcapGeometry()->getGeometry(pivot_.getDetId())->getPosition().perp();
      std::cout << EEDetId(pivot_.getDetId()) << std::endl;

      std::cout << " Test getClosestCell " << EcalEntrance_ << std::endl;
      DetId testcell = myCalorimeter->getClosestCell(EcalEntrance_, true, false);
      std::cout << " Result "<< EEDetId(testcell) << std::endl;
      std::cout << " Position " << myCalorimeter->getEcalEndcapGeometry()->getGeometry(testcell)->getPosition() << std::endl;
    }
    }

#endif

}

EcalHitMaker::~EcalHitMaker()
{
  if (myCrystalWindowMap_ != 0)
    {
      delete myCrystalWindowMap_;
    }
}

bool 
EcalHitMaker::addHitDepth(double r,double phi,double depth)
{
//  std::cout << " Add hit depth called; Current deph is "  << currentdepth_;
//  std::cout << " Required depth is " << depth << std::endl;
  depth+=X0depthoffset_;
  double sp(1.);
  r*=radiusFactor_;
  CLHEP::Hep2Vector point(r*std::cos(phi),r*std::sin(phi));

  unsigned xtal=fastInsideCell(point,sp);  
  //  if(cellid.isZero()) std::cout << " cell is Zero " << std::endl;
//  if(xtal<1000) 
//    {
//      std::cout << "Result " << regionOfInterest_[xtal].getX0Back() << " " ;
//      std::cout << depth << std::endl;
//    }
//  myHistos->fill("h5000",depth);
  if(xtal<1000)
    {
      //      myHistos->fill("h5002",regionOfInterest_[xtal].getX0Back(),depth);
      //      myHistos->fill("h5003",ecalentrance_.eta(),maxX0_);
      if(regionOfInterest_[xtal].getX0Back()>depth) 
    {
      hits_[xtal]+=spotEnergy;  
      //      myHistos->fill("h5005",r);
      return true;
    }
      else
    {
      rearleakage_+=spotEnergy;
    }
    }
  else
    {
//      std::cout << " Return false " << std::endl;
//      std::cout << " Add hit depth called; Current deph is "  << currentdepth_;
//      std::cout << " Required depth is " << depth << std::endl;
//      std::cout << " R = " << r << " " << radiusFactor_ << std::endl; 
    }

  outsideWindowEnergy_+=spotEnergy;
  return false;
}


//bool EcalHitMaker::addHitDepth(double r,double phi,double depth)
//{
//  std::map<unsigned,float>::iterator itcheck=hitMap_.find(pivot_.getDetId().rawId());
//  if(itcheck==hitMap_.end())
//    {
//      hitMap_.insert(std::pair<uint32_t,float>(pivot_.getDetId().rawId(),spotEnergy));
//    }
//  else
//    {
//      itcheck->second+=spotEnergy;
//    }
//  return true;
//}

bool 
EcalHitMaker::addHit(double r,double phi,unsigned layer)
{
  //  std::cout <<" Addhit " << std::endl;
  //  std::cout << " Before insideCell " << std::endl;
  double sp(1.);
  //  std::cout << " Trying to add " << r << " " << phi << " " << radiusFactor_ << std::endl; 
  r*=radiusFactor_;
  CLHEP::Hep2Vector point(r*std::cos(phi),r*std::sin(phi));
  //  std::cout << "point " << point << std::endl;
  //  CellID cellid=insideCell(point,sp);
  unsigned xtal=fastInsideCell(point,sp);  
  //  if(cellid.isZero()) std::cout << " cell is Zero " << std::endl;
  if(xtal<1000) 
    {
      if(sp==1.)
    hits_[xtal]+=spotEnergy;    
      else      
    hits_[xtal]+=(random->flatShoot()<sp)*spotEnergy;
      return true;
    }

  outsideWindowEnergy_+=spotEnergy;
//  std::cout << " This hit ; r= " << point << " hasn't been added "<<std::endl;
//  std::cout << " Xtal " << xtal << std::endl;
//  for(unsigned ip=0;ip<npadsatdepth_;++ip)
//    {
//      std::cout << padsatdepth_[ip] << std::endl;
//    }
  
  return false;
}

// Temporary solution
//bool EcalHitMaker::addHit(double r,double phi,unsigned layer)
//{
//  std::map<unsigned,float>::iterator itcheck=hitMap_.find(pivot_.getDetId().rawId());
//  if(itcheck==hitMap_.end())
//    {
//      hitMap_.insert(std::pair<uint32_t,float>(pivot_.getDetId().rawId(),spotEnergy));
//    }
//  else
//    {
//      itcheck->second+=spotEnergy;
//    }
//  return true;
//}


unsigned 
EcalHitMaker::fastInsideCell(const CLHEP::Hep2Vector & point,double & sp,bool debug) 
{

  //  debug = true;
  bool found=false;
  unsigned niter=0;
  // something clever has to be implemented here
  unsigned d1,d2;
  convertIntegerCoordinates(point.x(),point.y(),d1,d2);
  //  std::cout << "Fastinside cell " << point.x() << " " <<  point.y() << " " << d1 << " "<< d2 << " " << nx_ << " " << ny_ << std::endl;
  if(d1>=nx_||d2>=ny_) 
    {
      //      std::cout << " Not in the map " <<std::endl;
      return 9999;
    }
  unsigned cell=myCrystalNumberArray_[d1][d2];
  // We are likely to be lucky
  //  std::cout << " Got the cell " << cell << std::endl;
  if (validPads_[cell]&&padsatdepth_[cell].inside(point)) 
    {
      //      std::cout << " We are lucky " << cell << std::endl;
      sp = padsatdepth_[cell].survivalProbability();
      return cell;
    }
  
  //  std::cout << "Starting the loop " << std::endl;
  bool status(true);
  const std::vector<unsigned>& localCellVector(myCrystalWindowMap_->getCrystalWindow(cell,status));
  if(status) 
    {
      unsigned size=localCellVector.size();      
      //      std::cout << " Starting from " << EBDetId(regionOfInterest_[cell].getDetId()) << std::endl;
      //      const std::vector<DetId>& neighbours=myCalorimeter->getNeighbours(regionOfInterest_[cell].getDetId());
//      std::cout << " The neighbours are " << std::endl;
//      for(unsigned ic=0;ic<neighbours.size(); ++ic)
//  {
//    std::cout << EBDetId(neighbours[ic]) << std::endl;
//  }
//      std::cout << " Done " << std::endl;
      for(unsigned ic=0;ic<8&&ic<size;++ic)
    {
      unsigned iq=localCellVector[ic];
//    std::cout << " Testing " << EBDetId(regionOfInterest_[iq].getDetId()) << std::endl; ; 
//    std::cout << " " << iq << std::endl;
//    std::cout << padsatdepth_[iq] ;
      if(validPads_[iq]&&padsatdepth_[iq].inside(point))
        {
          //          std::cout << " Yes " << std::endl;
          //          myHistos->fill("h1000",niter);
          sp = padsatdepth_[iq].survivalProbability();
          //          std::cout << "Finished the loop " << niter << std::endl;
          //          std::cout << "Inside " << std::endl;
          return iq;
        }
      //      std::cout << " Not inside " << std::endl;
      //      std::cout << "No " << std::endl;
      ++niter;
    }
          }
  if(debug) std::cout << " not found in a quad, let's check the " << ncrackpadsatdepth_ << " cracks " << std::endl;
  //  std::cout << "Finished the loop " << niter << std::endl;
  // Let's check the cracks 
  //  std::cout << " Let's check the cracks " << ncrackpadsatdepth_ << " " << crackpadsatdepth_.size() << std::endl;
  unsigned iquad=0;
  unsigned iquadinside=999;
  while(iquad<ncrackpadsatdepth_&&!found)
    {
      //      std::cout << " Inside the while " << std::endl;
      if(crackpadsatdepth_[iquad].inside(point)) 
    {
      iquadinside=iquad;
      found=true;
      sp = crackpadsatdepth_[iquad].survivalProbability();
    }
      ++iquad;
      ++niter;
    }
  //  myHistos->fill("h1002",niter);
  if(!found&&debug) std::cout << " Not found in the cracks " << std::endl;
  return (found) ? crackpadsatdepth_[iquadinside].getNumber(): 9999;
}


void 
EcalHitMaker::setTrackParameters(const XYZNormal& normal,
                 double X0depthoffset,
                 const FSimTrack& theTrack)
{
  //  myHistos->debug("setTrackParameters");
  //  std::cout << " Track " << theTrack << std::endl;
  intersections_.clear();
  // This is certainly enough
  intersections_.reserve(50);
  myTrack_=&theTrack;
  normal_=normal.Unit();
  X0depthoffset_=X0depthoffset;
  cellLine(intersections_);
  buildSegments(intersections_);
//  std::cout << " Segments " << segments_.size() << std::endl;
//  for(unsigned ii=0; ii<segments_.size() ; ++ii)
//    {
//      std::cout << segments_[ii] << std::endl;
//    }


  // This is only needed in case of electromagnetic showers 
  if(EMSHOWER&&onEcal_&&ecalTotalX0()>0.)
    {
      //      std::cout << "Total X0  " << ecalTotalX0() << std::endl;
      for(unsigned ic=0;ic<ncrystals_;++ic)
    {
      for(unsigned idir=0;idir<4;++idir)
        {
          XYZVector norm=regionOfInterest_[ic].exitingNormal(CaloDirectionOperations::Side(idir));
          regionOfInterest_[ic].crystalNeighbour(idir).setToBeGlued((norm.Dot(normal_)<0.));
        }
      // Now calculate the distance in X0 of the back sides of the crystals
      // (only for EM showers)
      if(EMSHOWER)
        {
          XYZVector dir=regionOfInterest_[ic].getBackCenter()-segments_[ecalFirstSegment_].entrance();
          double dist=dir.Dot(normal_);
          double absciss= dist+segments_[ecalFirstSegment_].sEntrance();
          std::vector<CaloSegment>::const_iterator segiterator;
          // First identify the correct segment
          //      std::cout << " Crystal " << ic  << regionOfInterest_[ic].getBackCenter() ;
          //      std::cout << " Entrance : " << segments_[ecalFirstSegment_].entrance()<< std::endl;
          //      std::cout << " Looking for the segment " << dist << std::endl;
          segiterator = find_if(segments_.begin(),segments_.end(),CaloSegment::inSegment(absciss));
          //      std::cout << " Done " << std::endl;
          if(segiterator==segments_.end() )
        {
          // in this case, we won't have any problem. No need to 
          // calculate the real depth. 
          regionOfInterest_[ic].setX0Back(9999);
        }
          else
        {
          DetId::Detector det(segiterator->whichDetector());
          if(det!=DetId::Ecal)
            {
              regionOfInterest_[ic].setX0Back(9999);
            }
          else
            {
              double x0=segiterator->x0FromCm(dist);
              if(x0<maxX0_) maxX0_=x0;
              regionOfInterest_[ic].setX0Back(x0);
            }
        }
          //  myHistos->fill("h4000",ecalentrance_.eta(), regionOfInterest_[ic].getX0Back());
          //}
          //      else
          //{
          //   // in this case, we won't have any problem. No need to 
          //  // calculate the real depth. 
          //  regionOfInterest_[ic].setX0Back(9999);
          //}
        }//EMSHOWER  
    } // ndir
      //      myHistos->fill("h6000",segments_[ecalFirstSegment_].entrance().eta(),maxX0_);
    }
  //  std::cout << "Leaving setTrackParameters" << std::endl
}


void 
EcalHitMaker::cellLine(std::vector<CaloPoint>& cp) 
{
  cp.clear();
  //  if(myTrack->onVFcal()!=2)
  //    {
  if(!central_&&onEcal_&&simulatePreshower_) preshowerCellLine(cp);
  if(onEcal_)ecalCellLine(EcalEntrance_,EcalEntrance_+normal_,cp);
  //    }
  
  XYZPoint vertex(myTrack_->vertex().position().Vect());

  //sort the points by distance (in the ECAL they are not necessarily ordered)
  XYZVector dir(0.,0.,0.);
  if(myTrack_->onLayer1())
    {
      vertex=(myTrack_->layer1Entrance().vertex()).Vect();
      dir=myTrack_->layer1Entrance().Vect().Unit();
    }
  else if(myTrack_->onLayer2())
    {
      vertex=(myTrack_->layer2Entrance().vertex()).Vect();
      dir=myTrack_->layer2Entrance().Vect().Unit();
    }  
  else if(myTrack_->onEcal())
    {
      vertex=(myTrack_->ecalEntrance().vertex()).Vect();
      dir=myTrack_->ecalEntrance().Vect().Unit();
    }
  else if(myTrack_->onHcal())
    {
      vertex=(myTrack_->hcalEntrance().vertex()).Vect();
      dir=myTrack_->hcalEntrance().Vect().Unit();
    }
  else if(myTrack_->onVFcal()==2)
    {
      vertex=(myTrack_->vfcalEntrance().vertex()).Vect();
      dir=myTrack_->vfcalEntrance().Vect().Unit();
    }
  else
    {
      std::cout << " Problem with the grid " << std::endl;
    }
 
  // Move the vertex for distance comparison (5cm)
  vertex -= 5.*dir;
  CaloPoint::DistanceToVertex myDistance(vertex);
  sort(cp.begin(),cp.end(),myDistance);
  
  // The intersections with the HCAL shouldn't need to be sorted
  // with the N.I it is actually a source of problems
  hcalCellLine(cp);

//  std::cout << " Intersections ordered by distance to " << vertex << std::endl;
//
//  for (unsigned ic=0;ic<cp.size();++ic)
//    {
//      XYZVector t=cp[ic]-vertex;
//      std::cout << cp[ic] << " " << t.mag() << std::endl;
//    }
}


void 
EcalHitMaker::preshowerCellLine(std::vector<CaloPoint>& cp) const
{
  //  FSimEvent& mySimEvent = myEventMgr->simSignal();
  //  FSimTrack myTrack = mySimEvent.track(fsimtrack_);
  //  std::cout << "FsimTrack " << fsimtrack_<< std::endl;
  //  std::cout << " On layer 1 " << myTrack.onLayer1() << std::endl;
  // std::cout << " preshowerCellLine " << std::endl; 
  if(myTrack_->onLayer1())
    {
      XYZPoint point1=(myTrack_->layer1Entrance().vertex()).Vect();
      double phys_eta=myTrack_->layer1Entrance().eta();
      double cmthickness = 
    myCalorimeter->layer1Properties(1)->thickness(phys_eta);

      if(cmthickness>0)
    {
      XYZVector dir=myTrack_->layer1Entrance().Vect().Unit();
      XYZPoint point2=point1+dir*cmthickness;
      
      CaloPoint cp1(DetId::Ecal,EcalPreshower,1,point1);
      CaloPoint cp2(DetId::Ecal,EcalPreshower,1,point2);
      cp.push_back(cp1);
      cp.push_back(cp2);
    }
      else
    {
      //      std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1<< std::endl;
    }
    }

  //  std::cout << " On layer 2 " << myTrack.onLayer2() << std::endl;
  if(myTrack_->onLayer2())
    {
      XYZPoint point1=(myTrack_->layer2Entrance().vertex()).Vect();
      double phys_eta=myTrack_->layer2Entrance().eta();
      double cmthickness = 
    myCalorimeter->layer2Properties(1)->thickness(phys_eta);
      if(cmthickness>0)
    {
      XYZVector dir=myTrack_->layer2Entrance().Vect().Unit();
      XYZPoint point2=point1+dir*cmthickness;
      

      CaloPoint cp1(DetId::Ecal,EcalPreshower,2,point1);
      CaloPoint cp2(DetId::Ecal,EcalPreshower,2,point2);

      cp.push_back(cp1);
      cp.push_back(cp2);
    }
      else
    {
      //      std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1 << std::endl;
    }
    }
  //  std::cout << " Exit preshower CellLine " << std::endl;
}

void 
EcalHitMaker::hcalCellLine(std::vector<CaloPoint>& cp) const
{
  //  FSimEvent& mySimEvent = myEventMgr->simSignal();
  //  FSimTrack myTrack = mySimEvent.track(fsimtrack_);
  int onHcal=myTrack_->onHcal();

  if(onHcal<=2&&onHcal>0)
    {
      XYZPoint point1=(myTrack_->hcalEntrance().vertex()).Vect();

      double eta=point1.eta();
      // HCAL thickness in cm (assuming that the particle is coming from 000)
      double thickness= myCalorimeter->hcalProperties(onHcal)->thickness(eta);
      cp.push_back(CaloPoint(DetId::Hcal,point1));
      XYZVector dir=myTrack_->hcalEntrance().Vect().Unit();
      XYZPoint point2=point1+dir*thickness;

      cp.push_back(CaloPoint(DetId::Hcal,point2));
      
    }
  int onVFcal=myTrack_->onVFcal();
  if(onVFcal==2)
    {
      XYZPoint point1=(myTrack_->vfcalEntrance().vertex()).Vect();
      double eta=point1.eta();
      // HCAL thickness in cm (assuming that the particle is coming from 000)
      double thickness= myCalorimeter->hcalProperties(3)->thickness(eta);
      cp.push_back(CaloPoint(DetId::Hcal,point1));
      XYZVector dir=myTrack_->vfcalEntrance().Vect().Unit();
      if(thickness>0)
    {
      XYZPoint point2=point1+dir*thickness;
      cp.push_back(CaloPoint(DetId::Hcal,point2));
    }

    }
}

void 
EcalHitMaker::ecalCellLine(const XYZPoint& a,const XYZPoint& b,std::vector<CaloPoint>& cp) 
{
  //  std::vector<XYZPoint> corners;
  //  corners.resize(4);
  unsigned ic=0;
  double t;
  XYZPoint xp;
  DetId c_entrance,c_exit;
  bool entrancefound(false),exitfound(false);
  //  std::cout << " Look for intersections " << ncrystals_ << std::endl;
  //  std::cout << " regionOfInterest_ " << truncatedGrid_ << " " << regionOfInterest_.size() << std::endl;
  // try to determine the number of crystals to test
  // First determine the incident angle
  double angle=std::acos(normal_.Dot(regionOfInterest_[0].getAxis().Unit()));

  //  std::cout << " Normal " << normal_<< " Axis " << regionOfInterest_[0].getAxis().Unit() << std::endl;
  double backdistance=std::sqrt(regionOfInterest_[0].getAxis().mag2())*std::tan(angle);
  // 1/2.2cm = 0.45
//   std::cout << " Angle " << angle << std::endl;
//   std::cout << " Back distance " << backdistance << std::endl;
  unsigned ncrystals=(unsigned)(backdistance*0.45);
  unsigned highlim=(ncrystals+4);
  highlim*=highlim;
  if(highlim>ncrystals_) highlim=ncrystals_;
  //  unsigned lowlim=(ncrystals>2)? (ncrystals-2):0;
  //  std::cout << " Ncrys " << ncrystals << std::endl;
  
  
  while(ic<ncrystals_&&(ic<highlim||!exitfound))
    {
      // Check front side
      //      if(!entrancefound)
    {
      const Plane3D& plan=regionOfInterest_[ic].getFrontPlane();
//    XYZVector axis1=(plan.Normal());
//    XYZVector axis2=regionOfInterest_[ic].getFirstEdge();
      xp=intersect(plan,a,b,t,false);
      regionOfInterest_[ic].getFrontSide(corners);
      //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(0),axis1,axis2);
      //      if(pad.globalinside(xp)) 
      if(inside3D(corners,xp))
        {
          cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(),UP,xp));
          entrancefound=true;
          c_entrance=regionOfInterest_[ic].getDetId();
          //      myHistos->fill("j12",highlim,ic);
        }
    }
      
      // check rear side
    //  if(!exitfound)
    {
      const Plane3D& plan=regionOfInterest_[ic].getBackPlane();
//    XYZVector axis1=(plan.Normal());
//    XYZVector axis2=regionOfInterest_[ic].getFifthEdge();
      xp=intersect(plan,a,b,t,false);
      regionOfInterest_[ic].getBackSide(corners);
      //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(4),axis1,axis2);
      //      if(pad.globalinside(xp)) 
      if(inside3D(corners,xp))
        {
          cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(),DOWN,xp));
          exitfound=true;
          c_exit=regionOfInterest_[ic].getDetId();
          //          std::cout << " Crystal : " << ic << std::endl;
          //          myHistos->fill("j10",highlim,ic);
        }
    }
      
      if(entrancefound&&exitfound&&c_entrance==c_exit) return;
      // check lateral sides 
      for(unsigned iside=0;iside<4;++iside)
    {
      const Plane3D& plan=regionOfInterest_[ic].getLateralPlane(iside);
      xp=intersect(plan,a,b,t,false);
//    XYZVector axis1=(plan.Normal());
//    XYZVector axis2=regionOfInterest_[ic].getLateralEdge(iside);
      regionOfInterest_[ic].getLateralSide(iside,corners);
      //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(iside),axis1,axis2);
      //      if(pad.globalinside(xp)) 
      if(inside3D(corners,xp))
        {
          cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(),CaloDirectionOperations::Side(iside),xp)); 
          //          std::cout << cp[cp.size()-1] << std::endl;
        }     
    }
      // Go to next crystal 
      ++ic;
    }    
}


void 
EcalHitMaker::buildSegments(const std::vector<CaloPoint>& cp)
{
  //  myHistos->debug();
  //  TimeMe theT("FamosGrid::buildSegments");
  unsigned size=cp.size();
  if(size%2!=0) 
    {
      //      std::cout << " There is a problem " << std::endl;
      return;
    }
  //  myHistos->debug();
  unsigned nsegments=size/2;
  segments_.reserve(nsegments);
  if (size==0) return;
  // curv abs
  double s=0.;
  double sX0=0.;
  double sL0=0.;
  
  unsigned ncrossedxtals = 0;
  unsigned is=0;
  while(is<nsegments)
    {

      if(cp[2*is].getDetId()!=cp[2*is+1].getDetId()&&
     cp[2*is].whichDetector()!=DetId::Hcal&&
     cp[2*is+1].whichDetector()!=DetId::Hcal) 
    {
//    std::cout << " Problem with the segments " << std::endl;
//    std::cout << cp[2*is].whichDetector() << " " << cp[2*is+1].whichDetector() << std::endl;
//    std::cout << is << " " <<cp[2*is].getDetId().rawId() << std::endl; 
//    std::cout << (2*is+1) << " " <<cp[2*is+1].getDetId().rawId() << std::endl; 
      ++is;
      continue;
    }

      // Check if it is a Preshower segment - Layer 1 
      // One segment per layer, nothing between
      //      myHistos->debug("Just avant Preshower"); 
      if(cp[2*is].whichDetector()==DetId::Ecal && cp[2*is].whichSubDetector()==EcalPreshower && cp[2*is].whichLayer()==1)
    {
      if(cp[2*is+1].whichDetector()==DetId::Ecal && cp[2*is+1].whichSubDetector()==EcalPreshower && cp[2*is+1].whichLayer()==1)
        {
          CaloSegment preshsegment(cp[2*is],cp[2*is+1],s,sX0,sL0,CaloSegment::PS,myCalorimeter);
          segments_.push_back(preshsegment);
          //          std::cout << " Added (1-1)" << preshsegment << std::endl;
              s+=preshsegment.length();
          sX0+=preshsegment.X0length();
          sL0+=preshsegment.L0length();
          X0PS1_+=preshsegment.X0length();
          L0PS1_+=preshsegment.L0length();
        }
      else
        {
          std::cout << " Strange segment between Preshower1 and " << cp[2*is+1].whichDetector();
          std::cout << std::endl;
        }
      ++is;
      continue;
    }

      // Check if it is a Preshower segment - Layer 2 
      // One segment per layer, nothing between
      if(cp[2*is].whichDetector()==DetId::Ecal && cp[2*is].whichSubDetector()==EcalPreshower && cp[2*is].whichLayer()==2)
    {
      if(cp[2*is+1].whichDetector()==DetId::Ecal && cp[2*is+1].whichSubDetector()==EcalPreshower && cp[2*is+1].whichLayer()==2)
        {
          CaloSegment preshsegment(cp[2*is],cp[2*is+1],s,sX0,sL0,CaloSegment::PS,myCalorimeter);
          segments_.push_back(preshsegment);
          //          std::cout << " Added (1-2)" << preshsegment << std::endl;
              s+=preshsegment.length();
          sX0+=preshsegment.X0length();
          sL0+=preshsegment.L0length();
          X0PS2_+=preshsegment.X0length();
          L0PS2_+=preshsegment.L0length();

          // material between preshower and EE
          if(is<(nsegments-1) && cp[2*is+2].whichDetector()==DetId::Ecal && cp[2*is+2].whichSubDetector()==EcalEndcap)
                {
          CaloSegment gapsef(cp[2*is+1],cp[2*is+2],s,sX0,sL0,CaloSegment::PSEEGAP,myCalorimeter);
          segments_.push_back(gapsef);
          s+=gapsef.length();     
          sX0+=gapsef.X0length();                                                                                          
          sL0+=gapsef.L0length();                                                                                          
          X0PS2EE_+=gapsef.X0length();                        
          L0PS2EE_+=gapsef.L0length();      
          //          std::cout << " Created  a segment " << gapsef.length()<< " " << gapsef.X0length()<< std::endl;
        }         
        }
      else
        {
          std::cout << " Strange segment between Preshower2 and " << cp[2*is+1].whichDetector();
          std::cout << std::endl;
        }
      ++is;
      continue;
    }
      // Now deal with the ECAL
      // One segment in each crystal. Segment corresponding to cracks/gaps are added
      //      myHistos->debug("Just avant ECAL"); 
      if(cp[2*is].whichDetector()==DetId::Ecal && (cp[2*is].whichSubDetector()==EcalBarrel || cp[2*is].whichSubDetector()==EcalEndcap))
    {
      if(cp[2*is+1].whichDetector()==DetId::Ecal && (cp[2*is+1].whichSubDetector()==EcalBarrel || cp[2*is+1].whichSubDetector()==EcalEndcap) )
        {
          DetId cell2=cp[2*is+1].getDetId();
          // set the real entrance
          if (ecalFirstSegment_<0) ecalFirstSegment_=segments_.size();

          // !! Approximatiom : the first segment is always in a crystal
          if(cp[2*is].getDetId()==cell2)
        {
          CaloSegment segment(cp[2*is],cp[2*is+1],s,sX0,sL0,CaloSegment::PbWO4,myCalorimeter);
         segments_.push_back(segment);
         // std::cout << " Added (2)" << segment << std::endl;
         s+=segment.length();
         sX0+=segment.X0length();
         sL0+=segment.L0length(); 
         X0ECAL_+=segment.X0length();
         L0ECAL_+=segment.L0length();
         ++ncrossedxtals;
         ++is; 
        }
          else
        {
          std::cout << " One more bug in the segment " <<std::endl;
          ++is;
        }
          // Now check if a gap or crack should be added
          if(is>0 && is<nsegments)
        {         
          DetId cell3=cp[2*is].getDetId();
          if(cp[2*is].whichDetector()!=DetId::Hcal) 
            {
              // Crack inside the ECAL
              bool bordercrossing=myCalorimeter->borderCrossing(cell2,cell3);
              CaloSegment cracksegment(cp[2*is-1],cp[2*is],s,sX0,sL0,(bordercrossing)?CaloSegment::CRACK:CaloSegment::GAP,myCalorimeter);
              segments_.push_back(cracksegment);
              s+=cracksegment.length();
              sX0+=cracksegment.X0length();
              sL0+=cracksegment.L0length();
              X0ECAL_+=cracksegment.X0length();
              L0ECAL_+=cracksegment.L0length();
              //   std::cout <<" Added(3) "<< cracksegment << std::endl;
            }
          else
            {
              // a segment corresponding to ECAL/HCAL transition should be
              // added here
              CaloSegment cracksegment(cp[2*is-1],cp[2*is],s,sX0,sL0,CaloSegment::ECALHCALGAP,myCalorimeter);
              segments_.push_back(cracksegment);
              s+=cracksegment.length();
              sX0+=cracksegment.X0length();
              sL0+=cracksegment.L0length();
              X0EHGAP_+=cracksegment.X0length();
              L0EHGAP_+=cracksegment.L0length();
            }
        }
          continue;
        }
      else
        {
          std::cout << " Strange segment between " << cp[2*is].whichDetector();
          std::cout << " and " << cp[2*is+1].whichDetector() << std::endl;
          ++is;
          continue;
        }     
    }
      //      myHistos->debug("Just avant HCAL"); 
      // HCAL
      if(cp[2*is].whichDetector()==DetId::Hcal&&cp[2*is+1].whichDetector()==DetId::Hcal)
    {
      CaloSegment segment(cp[2*is],cp[2*is+1],s,sX0,sL0,CaloSegment::HCAL,myCalorimeter);
      segments_.push_back(segment);
      s+=segment.length();
      sX0+=segment.X0length();
      sL0+=segment.L0length(); 
      X0HCAL_+=segment.X0length();
      L0HCAL_+=segment.L0length();
      //      std::cout <<" Added(4) "<< segment << std::endl;
      ++is; 
    }
    }
//  std::cout << " PS1 " << X0PS1_ << " " << L0PS1_ << std::endl;
//  std::cout << " PS2 " << X0PS2_ << " " << L0PS2_ << std::endl;
//  std::cout << " ECAL " << X0ECAL_ << " " << L0ECAL_ << std::endl;
//  std::cout << " HCAL " << X0HCAL_ << " " << L0HCAL_ << std::endl;
  
  totalX0_ = X0PS1_+X0PS2_+X0PS2EE_+X0ECAL_+X0EHGAP_+X0HCAL_;
  totalL0_ = L0PS1_+L0PS2_+L0PS2EE_+L0ECAL_+L0EHGAP_+L0HCAL_;
  //  myHistos->debug("Just avant le fill"); 

  #ifdef DEBUGCELLLINE
  myHistos->fill("h200",fabs(EcalEntrance_.eta()),X0ECAL_);
  myHistos->fill("h210",EcalEntrance_.phi(),X0ECAL_);
  if(X0ECAL_<20)
    myHistos->fill("h212",EcalEntrance_.phi(),X0ECAL_);
//  if(X0ECAL_<1.) 
//    {
//      for(unsigned ii=0; ii<segments_.size() ; ++ii)
//  {
//    std::cout << segments_[ii] << std::endl;
//  }
//    }
  myHistos->fillByNumber("h30",ncrossedxtals,EcalEntrance_.eta(),X0ECAL_);
  
  double zvertex = myTrack_->vertex().position().z();

  myHistos->fill("h310",EcalEntrance_.eta(),X0ECAL_);
  if(X0ECAL_<22)   myHistos->fill("h410",EcalEntrance_.phi());
  myHistos->fill("h400",zvertex,X0ECAL_);
  #endif
  //  std::cout << " Finished the segments " << std::endl;
    
}

void 
EcalHitMaker::buildGeometry()
{
  configuredGeometry_ = false;
  ncrystals_ = CellsWindow_.size();
  // create the vector with of pads with the appropriate size
  padsatdepth_.resize(ncrystals_);
  
  // This is fully correct in the barrel. 
  ny_= phisize_;
  nx_=ncrystals_/ny_;
  std::vector<unsigned> empty;
  empty.resize(ny_,0);
  myCrystalNumberArray_.reserve((unsigned)nx_);
  for(unsigned inx=0;inx<(unsigned)nx_;++inx)
    {      
      myCrystalNumberArray_.push_back(empty);
    }

  hits_.resize(ncrystals_,0.);
  regionOfInterest_.clear();
  regionOfInterest_.resize(ncrystals_);
  validPads_.resize(ncrystals_);
  for(unsigned ic=0;ic<ncrystals_;++ic)
    {
      myCalorimeter->buildCrystal(CellsWindow_[ic],regionOfInterest_[ic]);
      regionOfInterest_[ic].setNumber(ic);
      DetIdMap_.insert(std::pair<DetId,unsigned>(CellsWindow_[ic],ic));     
    }
  
  // Computes the map of the neighbours
  myCrystalWindowMap_ = new CrystalWindowMap(myCalorimeter,regionOfInterest_);
}



// depth is in X0 , L0 (depending on EMSHOWER/HADSHOWER) or in CM if inCM
bool 
EcalHitMaker::getPads(double depth,bool inCm) 
{
  //std::cout << " New depth " << depth << std::endl;
  // The first time, the relationship between crystals must be calculated
  // but only in the case of EM showers

  if(EMSHOWER && !configuredGeometry_) configureGeometry();


  radiusFactor_ = (EMSHOWER) ? moliereRadius*radiusCorrectionFactor_:interactionLength;
  detailedShowerTail_ = false;
  if(EMSHOWER)
    currentdepth_ = depth+X0depthoffset_;
  else
    currentdepth_ = depth;
  
//  if(currentdepth_>maxX0_+ps1TotalX0()+ps2TotalX0()) 
//    {
//      currentdepth_=maxX0_+ps1TotalX0()+ps2TotalX0()-1.; // the -1 is for safety
//      detailedShowerTail_=true;
//    }
    
//  std::cout << " FamosGrid::getQuads " << currentdepth_ << " " << maxX0_ << std::endl;

  ncrackpadsatdepth_=0;
  
  xmin_=ymin_=999;
  xmax_=ymax_=-999;
  double locxmin,locxmax,locymin,locymax;

  // Get the depth of the pivot
  std::vector<CaloSegment>::const_iterator segiterator;
  // First identify the correct segment

  if(inCm) // centimeter
    {
      segiterator = find_if(segments_.begin(),segments_.end(),CaloSegment::inSegment(currentdepth_));
    }
  else
    {
      // EM shower 
      if(EMSHOWER)
    segiterator = find_if(segments_.begin(),segments_.end(),CaloSegment::inX0Segment(currentdepth_));
      
      //Hadron shower 
      if(HADSHOWER)
    segiterator = find_if(segments_.begin(),segments_.end(),CaloSegment::inL0Segment(currentdepth_));
    }
  if(segiterator==segments_.end()) 
    {
      std::cout << " FamosGrid: Could not go at such depth " << depth << std::endl;
      std::cout << " EMSHOWER " << EMSHOWER << std::endl;
      std::cout << " Track " << *myTrack_ << std::endl;
      std::cout << " Segments " << segments_.size() << std::endl;
      for(unsigned ii=0; ii<segments_.size() ; ++ii)
    {
      std::cout << segments_[ii] << std::endl;
    }
      
      return false;
    }
  //  std::cout << *segiterator << std::endl;
  
  if(segiterator->whichDetector()!=DetId::Ecal)
    {
      std::cout << " In  " << segiterator->whichDetector() << std::endl;
      //      buildPreshower();
      return false;
    }
  
  //  std::cout << *segiterator << std::endl;
  // get the position of the origin
  
  XYZPoint origin;
  if(inCm)
    {
      origin=segiterator->positionAtDepthincm(currentdepth_);
    }
  else
    {
      if(EMSHOWER)
    origin=segiterator->positionAtDepthinX0(currentdepth_);
      if(HADSHOWER)
    origin=segiterator->positionAtDepthinL0(currentdepth_);
    }
  //  std::cout << " currentdepth_ " << currentdepth_ << " " << origin << std::endl;
  XYZVector newaxis=pivot_.getFirstEdge().Cross(normal_);

//  std::cout  << "Normal " << normal_ << std::endl;
//  std::cout << " New axis " << newaxis << std::endl;

//  std::cout << " ncrystals  " << ncrystals << std::endl;
  plan_ = Plane3D((Vector)normal_,(Point)origin);

  unsigned nquads=0;
  double sign=(central_) ? -1.: 1.;
  Transform3DR trans((Point)origin,(Point)(origin+normal_),(Point)(origin+newaxis),
             Point(0,0,0), Point(0.,0.,sign),      Point(0.,1.,0.));
  for(unsigned ic=0;ic<ncrystals_;++ic)
    {
      //      std::cout << " Building geometry for " << regionOfInterest_[ic].getCellID() << std::endl;
      XYZPoint a,b;

      //      std::cout << " Origin " << origin << std::endl;

      //      std::vector<XYZPoint> corners;
      //      corners.reserve(4);
      double dummyt;
      bool hasbeenpulled=false;
      bool behindback=false;
      for(unsigned il=0;il<4;++il)
    {
      // a is the il-th front corner of the crystal. b is the corresponding rear corner
      regionOfInterest_[ic].getLateralEdges(il,a,b);
      
      // pull the surface if necessary (only in the front of the crystals) 
      XYZPoint aprime=a;
      if(pulled(origin,normal_,a)) 
        {
          b=aprime;
          hasbeenpulled=true;
        }
      
      // compute the intersection. 
      // Check that the intersection is in the [a,b] segment  if HADSHOWER
      // if EMSHOWER the intersection is calculated as if the crystals were infinite
      XYZPoint xx=(EMSHOWER)?intersect(plan_,a,b,dummyt,false):intersect(plan_,a,b,dummyt,true);
      
      if(dummyt>1) behindback=true;
      //      std::cout << " Intersect " << il << " " << a << " " << b << " " << plan_ << " " << xx << std::endl;
      // check that the intersection actually exists 
      if(xx.mag2()!=0)
        {
          corners[il] = xx;
        }     
    }    
      //      std::cout << " ncorners " << corners.size() << std::endl;
      if(behindback&&EMSHOWER) detailedShowerTail_=true;
      // If the quad is completly defined. Store it ! 
      if(corners.size()==4)
    {
      padsatdepth_[ic]=CrystalPad(ic,corners,trans,bfactor_,!central_);
      // Parameter to be tuned
      if(hasbeenpulled) padsatdepth_[ic].setSurvivalProbability(pulledPadProbability_);
      validPads_[ic]=true;
      ++nquads; 
      // In principle, this should be done after the quads reorganization. But it would cost one more loop
      //  quadsatdepth_[ic].extrems(locxmin,locxmax,locymin,locymax);
      //  if(locxmin<xmin_) xmin_=locxmin;
      //  if(locymin<ymin_) ymin_=locymin;
      //  if(locxmax>xmax_) xmax_=locxmax;
      //  if(locymax>ymax_) ymax_=locymax;
    }
      else
    {
      padsatdepth_[ic]=CrystalPad();
      validPads_[ic]=false;
    }
    }
  //  std::cout << " Number of quads " << quadsatdepth_.size() << std::endl; 
  if(doreorg_)reorganizePads();
  //  std::cout << "Finished to reorganize " << std::endl;
  npadsatdepth_=nquads;
  //  std::cout << " prepareCellIDMap " << std::endl;
  

  // Resize the Quads to allow for some numerical inaccuracy 
  // in the "inside" function
  for(unsigned ic=0;ic<ncrystals_;++ic) 
    {
    
      if (!validPads_[ic]) continue;

      if(EMSHOWER) padsatdepth_[ic].resetCorners();

      padsatdepth_[ic].extrems(locxmin,locxmax,locymin,locymax);
      if(locxmin<xmin_) xmin_=locxmin;
      if(locymin<ymin_) ymin_=locymin;
      if(locxmax>xmax_) xmax_=locxmax;
      if(locymax>ymax_) ymax_=locymax;
    }
  
  sizex_=(xmax_-xmin_)/nx_;
  sizey_=(ymax_-ymin_)/ny_;


  // Make sure that sizex_ and sizey_ are set before running prepareCellIDMap
  prepareCrystalNumberArray();

  //  std::cout << " Finished  prepareCellIDMap " << std::endl;
  ncrackpadsatdepth_=crackpadsatdepth_.size();

  return true;
}




void 
EcalHitMaker::configureGeometry()
{ 
  configuredGeometry_=true;
  for(unsigned ic=0;ic<ncrystals_;++ic)
    {
      //      std::cout << " Building " << cellids_[ic] << std::endl;
      for(unsigned idir=0;idir<8;++idir)
    {

      unsigned oppdir=CaloDirectionOperations::oppositeDirection(idir);
      // Is there something else to do ? 
      // The relationship with the neighbour may have been set previously.
          if(regionOfInterest_[ic].crystalNeighbour(idir).status()>=0) 
        {
          //          std::cout << " Nothing to do " << std::endl;
          continue ;
        }

      const DetId & oldcell(regionOfInterest_[ic].getDetId());
      CaloDirection dir=CaloDirectionOperations::neighbourDirection(idir);
      DetId newcell(oldcell);
      if(!myCalorimeter->move(newcell,dir))
        {
          // no neighbour in this direction
          regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
          continue;
        }
      // Determine the number of this neighbour
      //      std::cout << " The neighbour is " << newcell << std::endl;
      std::map<DetId,unsigned>::const_iterator niter(DetIdMap_.find(newcell));
      if(niter==DetIdMap_.end())
        {
          //          std::cout << " The neighbour is not in the map " << std::endl;
          regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
          continue;
        }
      // Now there is a neighbour     
      //      std::cout << " The neighbour is " << niter->second << " " << cellids_[niter->second] << std::endl;
      regionOfInterest_[ic].crystalNeighbour(idir).setNumber(niter->second);
      //      std::cout << " Managed to set crystalNeighbour " << ic << " " << idir << std::endl;
      //      std::cout << " Trying  " << niter->second << " " << oppdir << std::endl;
      regionOfInterest_[niter->second].crystalNeighbour(oppdir).setNumber(ic);
      //      std::cout << " Crack/gap " << std::endl;
      if(myCalorimeter->borderCrossing(oldcell,newcell))
        {
          regionOfInterest_[ic].crystalNeighbour(idir).setStatus(1);
          regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(1);
          //          std::cout << " Crack ! " << std::endl;
        }
      else
        {
          regionOfInterest_[ic].crystalNeighbour(idir).setStatus(0);
          regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(0);
          //          std::cout << " Gap" << std::endl;
        }       
    }
    }
    // Magnetic field a la Charlot
  double theta=EcalEntrance_.theta();
  if(theta>M_PI_2) theta=M_PI-theta;
  bfactor_=1./(1.+0.133*theta);
  if(myCalorimeter->magneticField()==0. ) bfactor_=1.;
}

// project fPoint on the plane (original,normal)
bool 
EcalHitMaker::pulled(const XYZPoint & origin, 
             const XYZNormal& normal, 
             XYZPoint & fPoint) const 
{
  // check if fPoint is behind the origin
  double dotproduct=normal.Dot(fPoint-origin);
  if(dotproduct<=0.) return false;

  //norm of normal is 1 
  fPoint -= (1+dotproduct)*normal;
  return true;
}


void 
EcalHitMaker::prepareCrystalNumberArray()
{
  for(unsigned iq=0;iq<npadsatdepth_;++iq)
    {
      if(!validPads_[iq]) continue;
      unsigned d1,d2;
      convertIntegerCoordinates(padsatdepth_[iq].center().x(),padsatdepth_[iq].center().y(),d1,d2);
      myCrystalNumberArray_[d1][d2]=iq;
    }
}

void EcalHitMaker::convertIntegerCoordinates(double x, double y,unsigned &ix,unsigned &iy) const 
{
  int tix=(int)((x-xmin_)/sizex_);
  int tiy=(int)((y-ymin_)/sizey_);
  ix=iy=9999;
  if(tix>=0) ix=(unsigned)tix;
  if(tiy>=0) iy=(unsigned)tiy;
}

const std::map<uint32_t,float>& EcalHitMaker::getHits() 
{
  if (hitmaphasbeencalculated_) return hitMap_;
  for(unsigned ic=0;ic<ncrystals_;++ic)
    {
      hitMap_.insert(std::pair<uint32_t,double>(regionOfInterest_[ic].getDetId().rawId(),hits_[ic]));
    }
  hitmaphasbeencalculated_=true;
  return hitMap_;
}



void 
EcalHitMaker::reorganizePads()
{

  // Some cleaning first 
  //  std::cout << " Starting reorganize " << std::endl;
  crackpadsatdepth_.clear();
  crackpadsatdepth_.reserve(etasize_*phisize_);
  ncrackpadsatdepth_=0;
  std::vector<neighbour> gaps;
  std::vector<std::vector<neighbour> > cracks;  
  //  cracks.clear();
  cracks.resize(ncrystals_);


  for(unsigned iq=0;iq<ncrystals_;++iq)
    {
      if(!validPads_[iq]) continue;

      gaps.clear();
      //   std::cout << padsatdepth_[iq] << std::endl;
      //check all the directions
      for(unsigned iside=0;iside<4;++iside)
    {
      //      std::cout << " To be glued " << iside << " " << regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued() << std::endl;
      CaloDirection thisside=CaloDirectionOperations::Side(iside);
      if(regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued())
        {
          // look for the neighbour and check that it exists
          int neighbourstatus=regionOfInterest_[iq].crystalNeighbour(iside).status();
          if(neighbourstatus<0)
        continue;

          unsigned neighbourNumber=regionOfInterest_[iq].crystalNeighbour(iside).number();
          if(!validPads_[neighbourNumber]) continue;
          // there is a crack between 
          if(neighbourstatus==1)
        {
          //          std::cout << " 1 Crack : " << thisside << " " << cellids_[iq]<< " " << cellids_[neighbourNumber] << std::endl;
          cracks[iq].push_back(neighbour(thisside,neighbourNumber));
        } // else it is a gap 
          else
        {
          gaps.push_back(neighbour(thisside,neighbourNumber));
        }
        }
    }
      // Now lift the gaps
      gapsLifting(gaps,iq);
    }

  unsigned ncracks=cracks.size();
  //  std::cout << " Cracks treatment : " << cracks.size() << std::endl;
  for(unsigned icrack=0;icrack<ncracks;++icrack)
    {
      //      std::cout << " Crack number " << crackiter->first << std::endl;
      cracksPads(cracks[icrack],icrack);
    }

}

//dir 2 = N,E,W,S
CLHEP::Hep2Vector & 
EcalHitMaker::correspondingEdge(neighbour& myneighbour,CaloDirection dir2 ) 
{
  CaloDirection dir=CaloDirectionOperations::oppositeSide(myneighbour.first);
  CaloDirection corner=CaloDirectionOperations::add2d(dir,dir2);
  //  std::cout << "Corresponding Edge " << dir<< " " << dir2 << " " << corner << std::endl;
  return padsatdepth_[myneighbour.second].edge(corner);  
}

bool EcalHitMaker::diagonalEdge(unsigned myPad, CaloDirection dir,CLHEP::Hep2Vector & point)
{
  unsigned idir=CaloDirectionOperations::neighbourDirection(dir);
  if(regionOfInterest_[myPad].crystalNeighbour(idir).status()<0)
    return false;
  unsigned nneighbour=regionOfInterest_[myPad].crystalNeighbour(idir).number();
  if(!validPads_[nneighbour])
    {
      //      std::cout << " Wasn't able to move " << std::endl;
      return false;
    }
  point = padsatdepth_[nneighbour].edge(CaloDirectionOperations::oppositeSide(dir));
  return true;
}

bool EcalHitMaker::unbalancedDirection(const std::vector<neighbour>& dirs,unsigned & unb,unsigned & dir1, unsigned & dir2)
{
  if(dirs.size()==1) return false;
  if(dirs.size()%2==0) return false;
  CaloDirection tmp;
  tmp=CaloDirectionOperations::add2d(dirs[0].first,dirs[1].first);
  if(tmp==NONE) 
    {
      unb=2;
      dir1=0;
      dir2=1;
      return true;
    }
  tmp=CaloDirectionOperations::add2d(dirs[0].first,dirs[2].first);
  if(tmp==NONE)
    { 
      unb=1;
      dir1=0;
      dir2=2;
      return true;
    }
  unb=0;
  dir1=1;
  dir2=2;
  return true;
}


void 
EcalHitMaker::gapsLifting(std::vector<neighbour>& gaps,unsigned iq)
{
  //  std::cout << " Entering gapsLifting "  << std::endl;
  CrystalPad & myPad = padsatdepth_[iq];
  unsigned ngaps=gaps.size();
  static bool debug=false;
  if(ngaps==1)
    {
      if(debug)
    {
      std::cout << " Avant " << ngaps << " " <<gaps[0].first<< std::endl;
      std::cout << myPad << std::endl;
    }
      if(gaps[0].first==NORTH||gaps[0].first==SOUTH)
    {
      CaloDirection dir1=CaloDirectionOperations::add2d(gaps[0].first,EAST);
      CaloDirection dir2=CaloDirectionOperations::add2d(gaps[0].first,WEST);
      myPad.edge(dir1)=correspondingEdge(gaps[0],EAST);
      myPad.edge(dir2)=correspondingEdge(gaps[0],WEST);
    }
      else
    {
      CaloDirection dir1=CaloDirectionOperations::add2d(gaps[0].first,NORTH);
      CaloDirection dir2=CaloDirectionOperations::add2d(gaps[0].first,SOUTH);
      myPad.edge(dir1)=correspondingEdge(gaps[0],NORTH);
      myPad.edge(dir2)=correspondingEdge(gaps[0],SOUTH);
    }
      if(debug)
    {
          std::cout << " Apres " << std::endl;
          std::cout << myPad << std::endl;
    }
    }
  else
    if(ngaps==2)
      {
    if(debug)
      {
        std::cout << " Avant " << ngaps << " " <<gaps[0].first<< " " <<gaps[1].first << std::endl;
        std::cout << myPad << std::endl;
        std::cout << " Voisin 1 " << (gaps[0].second) << std::endl;
        std::cout << " Voisin 2 " << (gaps[1].second) << std::endl;
      }
    CaloDirection corner0=CaloDirectionOperations::add2d(gaps[0].first,gaps[1].first);

    CLHEP::Hep2Vector point(0.,0.);
    if(corner0!=NONE&&diagonalEdge(iq,corner0,point))
      {
        CaloDirection corner1=CaloDirectionOperations::add2d(CaloDirectionOperations::oppositeSide(gaps[0].first),gaps[1].first);
        CaloDirection corner2=CaloDirectionOperations::add2d(gaps[0].first,CaloDirectionOperations::oppositeSide(gaps[1].first));
        myPad.edge(corner0) = point;
        myPad.edge(corner1) = correspondingEdge(gaps[1],CaloDirectionOperations::oppositeSide(gaps[0].first));
        myPad.edge(corner2) = correspondingEdge(gaps[0],CaloDirectionOperations::oppositeSide(gaps[1].first));
      }
    else
      if(corner0==NONE)
        {
          if(gaps[0].first==EAST||gaps[0].first==WEST)
        {
          CaloDirection corner1=CaloDirectionOperations::add2d(gaps[0].first,NORTH);
          CaloDirection corner2=CaloDirectionOperations::add2d(gaps[0].first,SOUTH);
          myPad.edge(corner1)=correspondingEdge(gaps[0],NORTH);
          myPad.edge(corner2)=correspondingEdge(gaps[0],SOUTH);
          
          corner1=CaloDirectionOperations::add2d(gaps[1].first,NORTH);
          corner2=CaloDirectionOperations::add2d(gaps[1].first,SOUTH);
          myPad.edge(corner1)=correspondingEdge(gaps[1],NORTH);
          myPad.edge(corner2)=correspondingEdge(gaps[1],SOUTH);
        }
          else
        {
          CaloDirection corner1=CaloDirectionOperations::add2d(gaps[0].first,EAST);
          CaloDirection corner2=CaloDirectionOperations::add2d(gaps[0].first,WEST);
          myPad.edge(corner1)=correspondingEdge(gaps[0],EAST);
          myPad.edge(corner2)=correspondingEdge(gaps[0],WEST);
          
          corner1=CaloDirectionOperations::add2d(gaps[1].first,EAST);
          corner2=CaloDirectionOperations::add2d(gaps[1].first,WEST);
          myPad.edge(corner1)=correspondingEdge(gaps[1],EAST);
          myPad.edge(corner2)=correspondingEdge(gaps[1],WEST);
        }
        }
    if(debug)
      {
        std::cout << " Apres " << std::endl;
        std::cout << myPad << std::endl;
      }
      }
    else
      if(ngaps==3)
    {
      // in this case the four corners have to be changed 
      unsigned iubd,idir1,idir2;
      CaloDirection diag;
      CLHEP::Hep2Vector point(0.,0.);
      //      std::cout << " Yes : 3 gaps" << std::endl;
      if(unbalancedDirection(gaps,iubd,idir1,idir2))        
        {
          CaloDirection ubd(gaps[iubd].first),dir1(gaps[idir1].first);
          CaloDirection dir2(gaps[idir2].first);
          
          //          std::cout << " Avant " << std::endl << myPad << std::endl;
          //          std::cout << ubd << " " << dir1 << " " << dir2 << std::endl;
          diag=CaloDirectionOperations::add2d(ubd,dir1);
          if(diagonalEdge(iq,diag,point))
        myPad.edge(diag)=point;
          diag=CaloDirectionOperations::add2d(ubd,dir2);
          if(diagonalEdge(iq,diag,point))
        myPad.edge(diag)=point;
          CaloDirection oppside=CaloDirectionOperations::oppositeSide(ubd);
          myPad.edge(CaloDirectionOperations::add2d(oppside,dir1))=correspondingEdge(gaps[idir1],oppside);
          myPad.edge(CaloDirectionOperations::add2d(oppside,dir2))=correspondingEdge(gaps[idir2],oppside);
          //          std::cout << " Apres " << std::endl << myPad << std::endl;
        }                 
    }
      else
    if(ngaps==4)
      {
        //      std::cout << " Waouh :4 gaps" << std::endl;
        //      std::cout << " Avant " << std::endl;
        //      std::cout << myPad<< std::endl;
        CLHEP::Hep2Vector point(0.,0.);
        if(diagonalEdge(iq,NORTHEAST,point)) 
          myPad.edge(NORTHEAST)=point;
        if(diagonalEdge(iq,NORTHWEST,point)) 
          myPad.edge(NORTHWEST)=point;
        if(diagonalEdge(iq,SOUTHWEST,point)) 
          myPad.edge(SOUTHWEST)=point;
        if(diagonalEdge(iq,SOUTHEAST,point)) 
          myPad.edge(SOUTHEAST)=point;  
        //      std::cout << " Apres " << std::endl;
        //      std::cout << myPad<< std::endl;
      }   
}

void 
EcalHitMaker::cracksPads(std::vector<neighbour> & cracks, unsigned iq)
{
  //  std::cout << " myPad " << &myPad << std::endl;
  unsigned ncracks=cracks.size();
  CrystalPad & myPad = padsatdepth_[iq];
  for(unsigned ic=0;ic<ncracks;++ic)
    {
      //      std::vector<CLHEP::Hep2Vector> mycorners;
      //      mycorners.reserve(4);
      switch(cracks[ic].first)
    {
    case NORTH:
      {
        mycorners[0] = (padsatdepth_[cracks[ic].second].edge(SOUTHWEST));
        mycorners[1] = (padsatdepth_[cracks[ic].second].edge(SOUTHEAST));
        mycorners[2] = (myPad.edge(NORTHEAST));
        mycorners[3] = (myPad.edge(NORTHWEST));
      }
      break;
    case SOUTH:
      {
        mycorners[0] = (myPad.edge(SOUTHWEST));
        mycorners[1] = (myPad.edge(SOUTHEAST));
        mycorners[2] = (padsatdepth_[cracks[ic].second].edge(NORTHEAST));
        mycorners[3] = (padsatdepth_[cracks[ic].second].edge(NORTHWEST));
      }
      break;
    case EAST:
      {
        mycorners[0] = (myPad.edge(NORTHEAST));
        mycorners[1] = (padsatdepth_[cracks[ic].second].edge(NORTHWEST));
        mycorners[2] = (padsatdepth_[cracks[ic].second].edge(SOUTHWEST));
        mycorners[3] = (myPad.edge(SOUTHEAST));
      }
      break;
    case WEST:
      {
        mycorners[0] = (padsatdepth_[cracks[ic].second].edge(NORTHEAST));
        mycorners[1] = (myPad.edge(NORTHWEST));
        mycorners[2] = (myPad.edge(SOUTHWEST));
        mycorners[3] = (padsatdepth_[cracks[ic].second].edge(SOUTHEAST));
      }
      break;
    default:
      {
      }
    }  
      CrystalPad crackpad(ic,mycorners);
      // to be tuned. A simpleconfigurable should be used
      crackpad.setSurvivalProbability(crackPadProbability_);   
      crackpadsatdepth_.push_back(crackpad);
    }
  //  std::cout << " Finished cracksPads " << std::endl;
}


bool EcalHitMaker::inside3D(const std::vector<XYZPoint>& corners, 
                const XYZPoint& p) const
{
  // corners and p are in the same plane
  // p is inside "corners" if the four crossproducts (corners[i]xcorners[i+1]) 
  // are in the same direction

  XYZVector crossproduct(0.,0.,0.),previouscrossproduct(0.,0.,0.);

  for(unsigned ip=0;ip<4 ; ++ip) {
    crossproduct = (corners[ip]-p).Cross(corners[(ip+1)%4]-p);
    if(ip==0)
    previouscrossproduct=crossproduct;
    else
      if (crossproduct.Dot(previouscrossproduct)<0.) return false;       
  }

  return true;
}