namespace Geometry{

  const int    nPhiCEM = 24;
  const double dPhiCEM = 2.0*TMath::Pi()/nPhiCEM;

  const double maxEtaCEM = 1.1;
  const int    nEtaCEM   = 20;
  const double dEtaCEM   = 2.*maxEtaCEM/nEtaCEM;

  const double R_CEM     = 175.26;

  const double xTowWidth = 48.48;

  // Muon is 2-tower if closer than 1.58 cm from edge of tower
  const double MuonDZ  = 1.58;
  //  Deta
  //  1
  //  0       Dphi -1 0 1
  // -1
  // val of 0 means knock it out
  const int knock_region_el[][3]={{1,0,0},
				  {0,0,0},
				  {1,0,0}};

  const int knock_region_mu[][3]={{1,0,1},
				  {1,0,1},
				  {1,0,1}};

}

namespace Mip{

  const double mip_avg  = -0.533;
  const double mip_sig  =  0.117;
  const double mip_zero =  0.01;

  const double UE_avg   =  0.149;
}


//----------------------------------------------------------
double getPhi(const TLorentzVector& p){

  return TVector2::Phi_0_2pi(p.Phi());
}

//----------------------------------------------------------
double getEta(const TLorentzVector& p){

  return p.Eta();
}
//----------------------------------------------------------
double getTheta(const TLorentzVector& p){

  return p.Theta();
}
//----------------------------------------------------------
double getZ(const double eta, const double R){

  double z   = R*TMath::SinH(eta);

  return z;
}
//----------------------------------------------------------
int getiPhi(const double phi){

  int iPhi  = (int)(phi/Geometry::dPhiCEM);
  if(iPhi > Geometry::nPhiCEM-1)  iPhi = Geometry::nPhiCEM-1;
  if(iPhi < 0)                    iPhi = 0;

  return iPhi;
}
//----------------------------------------------------------
int getiEta(const double eta){

  int   ieta = -100;

  if (TMath::Abs(eta) < Geometry::maxEtaCEM){
    ieta=10+TMath::Floor(eta/Geometry::dEtaCEM);
  }

  return ieta;
}
//----------------------------------------------------------
double minZTow(int ieta){
  
  // min etas:
  //  0 ... -1.1 
  // 19 ...  1.1 - dEtaCEM

  double mineta = -Geometry::maxEtaCEM+ieta*Geometry::dEtaCEM;
  double z      = Geometry::R_CEM*TMath::SinH(mineta);
  
  return z;
}
double maxZTow(int ieta){
  
  double maxeta = -Geometry::maxEtaCEM+(ieta+1)*Geometry::dEtaCEM;
  double z      = Geometry::R_CEM*TMath::SinH(maxeta);
  
  return z;
}
//----------------------------------------------------------
double getLocalX(const double phi){

  double phiMod = TMath::Abs(fmod(phi,Geometry::dPhiCEM));
  double locx   = Geometry::xTowWidth * (0.5 - phiMod / Geometry::dPhiCEM);

  return locx;
}
//----------------------------------------------------------
double getLocalZ(const double eta){

  double locz = -99999.0;
  int    ieta = getiEta(eta);

  if (ieta>=0 && ieta < 20){

    double minztow = minZTow(ieta);
    double maxztow = maxZTow(ieta);
    double       z = getZ(eta , Geometry::R_CEM);
    
    double dz1     = z - minztow;
    double dz2     = maxztow - z;
    
    // dz is positive if the nearest tower edge is to the left
    locz = (dz1 < dz2) ? dz1 : -dz2 ;
  }
  return locz;
}
//----------------------------------------------------------
// returns 0 if the photon in a crack between towers

double cracks(double localx, double localz) {

  double cracks = 1.0;
  if (TMath::Abs(localx) > 23.1) cracks = 0.0;
  if (TMath::Abs(localz) <  4.2) cracks = 0.0;

  return cracks;
}

//----------------------------------------------------------
int MuonPhoton(const TLorentzVector& pm, 
	       const TLorentzVector& pg){

  // 0 into the recoil
  // 1 merge with the muon
  // 2 knock out
  // 3 in a crack

  int retval=0;
  
  double meta=getEta(pm);
  double mphi=getPhi(pm);
  double geta=getEta(pg);
  double gphi=getPhi(pg);

  int iphi1=getiPhi(mphi);
  int iphi2=getiPhi(gphi);

  int ieta1=getiEta(meta);
  int ieta2=getiEta(geta);
  
  int mu_dz=getLocalZ(meta);
  int ph_dz=getLocalZ(geta);

  int ph_x =getLocalX(gphi);
  

  if (ieta1==-100 || ieta2==-100) retval=0;
  else {
    
    if (cracks(ph_x,ph_dz)) retval=3;
    else{
    
      int dieta = ieta1-ieta2;
      int diphi = iphi1-iphi2;
      if (diphi < -12) {diphi += 24;}
      if (diphi >  12) {diphi -= 24;}
      
      
      //merging
      //-------------------------------------------------------------
      if(diphi == 0) {
	// if em energy is in the same tower as muon, add its energy
	if(dieta == 0) retval = 1;  // one-tower muon
	
	else if (abs(dieta)==1 && 
		 (dieta*mu_dz) < Geometry::MuonDZ && 
		 (dieta*mu_dz) > 0.)
	  retval = 1;  // two-tower muon
      }
      //Knocking out
      //-------------------------------------------------------------
      
      if (TMath::Abs(diphi)<2 && 
	  TMath::Abs(dieta)<2 && retval != 1)
	
	if (Geometry::knock_region_mu[1-dieta][1+diphi]==0) retval = 2;
      
    }
  }
  return retval;

}

//----------------------------------------------------------

int ElectronPhoton(const TLorentzVector& pe, 
		   const TLorentzVector& pg){

  // 0 into the recoil
  // 1 merge with the muon
  // 2 knock out
  // 3 in the crack

  int retval=0;
  
  double eeta=getEta(pe);
  double ephi=getPhi(pe);
  double geta=getEta(pg);
  double gphi=getPhi(pg);

  int iphi1=getiPhi(ephi);
  int iphi2=getiPhi(gphi);

  int ieta1=getiEta(eeta);
  int ieta2=getiEta(geta);

  double eledz=getLocalZ(eeta);
  double phodz=getLocalZ(geta);

  double pho_x=getLocalX(gphi);



  if (ieta1==-100 || ieta2==-100) retval=0;
  else {
    
    if (cracks(pho_x,phodz)) retval=3;
    else{


      // check the sign vs ces x etc
      
      int dieta = ieta2-ieta1;
      int diphi = iphi2-iphi1;
      if (diphi < -12) {diphi+=24;}
      if (diphi >  12) {diphi-=24;}
      
      //merging
      //-------------------------------------------------------------
      // Algorithm as implemented here :
      // Merge photons with d-eta = -1,0,+1 in the same phi wedge. 
      // Exceptions are tower 0 (clusters cannot extend across crack) 
      // and tower 9 (PEM towers cannot enter CEM clusters).
      
      if ( (ieta1>9 && ieta2>9) || (ieta1<10 && ieta2<10) ){
	
	if (diphi == 0){
	  
	  if(dieta == 0) retval = 1;
	  
	  else if ( abs(dieta)==1     && 
		    (-dieta*eledz) > 0.) retval = 1;
	}
      }
      
      //Knocking out
      //-------------------------------------------------------------    
      if (retval != 1){
	
	// the sign of local x, 
	// which is falling towards increasing iphi, so we have a -sign
	// lxces goes from -24 to 24cm when phi increases
	double lxces = -getLocalX(ephi);
	lxces/=TMath::Abs(lxces);
	
	// if dphi>0 means photon on the side of higher phi, 
	// so electron is closer for positive lxces
	// if dphi<0 lxces is closer at negative lxces
	diphi*=(int)lxces;
	
	if (TMath::Abs(diphi)<2 && 
	    TMath::Abs(dieta)<2    )
	  if (Geometry::knock_region_el[1-dieta][1+diphi]==0) retval = 2;
      }
      
    }
  } 
  return retval;
}
  
//----------------------------------------------------------
// Muon energy measured in the EM towers
// MIP contribution + underlying event + photons
// This routine returns MIP + UE
// for photons: check if they need to be added 
// using the return value (=1) of the MuonPhoton routine

double MipE(){

  double mipe = Mip::UE_avg;
  
  double ifzero = gRandom->Rndm();

  if (ifzero > Mip::mip_zero) // 0.01% have zero MIP contribution
    mipe += TMath::Power(10.,gRandom->Gaus(Mip::mip_avg,Mip::mip_sig));
  
  return mipe;
}
//----------------------------------------------------------

bool MipCutFail(const double e_em, const double pt){

  bool fail=false;

  if (pt <= 100) {if (e_em  > 2.0) fail=true;}
  if (pt >  100) {if (e_em  >(2.0 + (pt -100.0 )*0.0115)) fail=true;}


  return fail;
}

//----------------------------------------------------------


void simplesim(){


  TLorentzVector pm = TLorentzVector(1.,0.0,1.00,TMath::Sqrt(2));
  TLorentzVector pg = TLorentzVector(1.,-0.2,1.02,TMath::Sqrt(2));

  cout << ElectronPhoton(pm,pg) << endl;
  cout << MuonPhoton    (pm,pg) << endl;


  // This generates a histogram of MIP energies using the
  // gaussian with the parameters from Mip:: ...

  TH1F* mip=new TH1F("mip","MIP energies without photons",100,0.,3.);

  for (int i=0; i<20000; i++) mip->Fill(MipE());
  
  mip->Draw();

}