/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \-------------  arXiv:0903.2225v1  ------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
**   FROG: [hep-ex/0901.2718v1]                                       **
**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   Main authors:                                                    **
**   -------------                                                    **
**                                                                    **
**                Severine Ovyn                Xavier Rouby           **
**          severine.ovyn@uclouvain.be      xavier.rouby@cern         **
**                                                                    **
**         Center for Particle Physics and Phenomenology (CP3)        **
**                Universite catholique de Louvain (UCL)              **
**                     Louvain-la-Neuve, Belgium                      **
**                                                                    **
**                      Copyright (C) 2008-2009,                      **
**                        All rights reserved.                        **
**                                                                    **
***********************************************************************/

/// \file D_CaloUtil.cc
/// \brief Implementation of D_CaloElement, D_CaloResolution and D_CaloList classes

#include "CaloUtil.h"
#include <cmath>
#include <algorithm>
#include <iostream>
#include "TRandom.h"

D_CaloElement& D_CaloElement::operator=(const D_CaloElement& calo) {
  if(this == &calo) return *this;
   _name   = calo._name;
   _etamin = calo._etamin;
   _etamax = calo._etamax;
   _resol_em = calo._resol_em;
   _resol_had= calo._resol_had;
  return *this;
}

D_CaloResolution& D_CaloResolution::operator=(const D_CaloResolution& resolution) {
  if(this == &resolution) return *this;
  _c = resolution._c;
  _n = resolution._n;
  _s = resolution._s;
  return *this;
}

void D_CaloResolution::print() const {
  std::cout << "C: " << _c << "  N:" << _n << "  S:" << _s << std::endl;
}

const D_CaloElement& D_CaloList::getElement(const float eta) const {
   for (unsigned int i=0; i<_calorimeters.size(); i++) {
      if(_calorimeters[i].getEtamin() < eta && eta < _calorimeters[i].getEtamax()) return _calorimeters[i];
   } 
   // if not found
   return dummyCalo;
}

void D_CaloList::addElement(const D_CaloElement & calo) {
  _sorted = false;
  _calorimeters.push_back(calo);
}

void D_CaloList::sortElements() {
  // here: it sorts the elements with respect to their eta
  sort(_calorimeters.begin(),_calorimeters.end(),ordering());
 
   for (unsigned int i=0; i<_calorimeters.size()-1; i++) {
	if(_calorimeters[i].getEtamax() > _calorimeters[i+1].getEtamin() ) {
		_etamin = UNDEFINED;
		_etamax = UNDEFINED;
		cerr << "** Error: the calorimeter list should never contain overlapping    **"<< endl;
		cerr << "**    calorimeters.                                                **"<< endl;
		return;
	}		
   }
  _etamin = _calorimeters[0].getEtamin();
  _etamax = _calorimeters[_calorimeters.size()-1].getEtamax();
  _sorted = true;
}

D_CaloList& D_CaloList::operator=(const D_CaloList& list) {
  if(this == &list) return *this;
    _calorimeters = list._calorimeters;
    _etamin       = list._etamin;
    _etamax       = list._etamax;
    _sorted       = list._sorted;
  return *this;
}

void D_CaloList::print() const {
  if(_sorted) cout << "The list has been sorted\n"; 
  else cout << "The list has not been sorted\n";

  for (unsigned int i=0; i<_calorimeters.size(); i++) {
	std::cout << _calorimeters[i].getName() << " " << _calorimeters[i].getEtamin() << 
                                                   " " << _calorimeters[i].getEtamax() << std::endl;	  
        _calorimeters[i].getElectromagneticResolution().print();
        _calorimeters[i].getHadronicResolution().print();
  }

}








D_CaloTower::D_CaloTower(const float iEta, const float iPhi) : 
	_eta(iEta), _phi(iPhi), _Eem(UNDEFINED), _Ehad(UNDEFINED), _E(UNDEFINED), _ET(UNDEFINED) {}

D_CaloTower::D_CaloTower(const D_CaloTower& tower) : 
	_eta(tower._eta), _phi(tower._phi), _Eem(tower._Eem), 
        _Ehad(tower._Ehad), _E(tower._E), _ET(tower._ET) {}

D_CaloTower& D_CaloTower::operator=(const D_CaloTower& tower) {
   if(this==&tower) return *this;
	_eta  = tower._eta; 	_phi  = tower._phi; 
        _Eem  = tower._Eem; 	_Ehad = tower._Ehad;
        _E    = tower._E; 	_ET   = tower._ET;
   return *this;
}


const float D_CaloResolution::Smear(const float energy) const {
    float sE = gRandom->Gaus(energy,
           sqrt( pow(_n,2) +
                 pow(_c*energy,2) +
                 pow(_s*sqrt(energy),2) ));
    return (sE>0)? sE : 0;
}

const bool D_CaloTower::operator==(const D_CaloTower& el) const { if (this == &el) return true; return this->_eta==el._eta && this->_phi == el._phi;}

void D_CaloTowerList::addTower(const D_CaloTower& tower) {
  _calotowers.push_back(tower);
  _sorted = false;
  _merged = false;
}

D_CaloTowerList& D_CaloTowerList::operator=(const D_CaloTowerList& list) {
  if(this==&list) return *this;
     _calotowers = list._calotowers;
     _sorted = list._sorted;
     _merged = list._merged;
  return *this;
}

void D_CaloTowerList::sortElements() {
  _sorted = true;
  if(_calotowers.size()==0) return;
  sort(_calotowers.begin(),_calotowers.end(),ordering());
}

void D_CaloTowerList::mergeDuplicates() {
  if(!_sorted) sortElements();
  _merged = true;
  if (_calotowers.size() ==0) { return;}
 
 vector<D_CaloTower> temp;
 temp.push_back(_calotowers[0]);
 for (unsigned int i=1; i<_calotowers.size(); i++) {

   if(temp[temp.size()-1] == _calotowers[i]) { 
        //i.e. if towers have the same eta/phi
	float em = temp[temp.size()-1].getEem() + _calotowers[i].getEem();
        float had= temp[temp.size()-1].getEhad()+ _calotowers[i].getEhad();
        temp[temp.size()-1].Set_Eem_Ehad_E_ET(em,had);
   }
   else temp.push_back(_calotowers[i]); 

 } 
   _calotowers = temp;
}

void D_CaloTowerList::print() const {
  for (unsigned int i=0; i<_calotowers.size(); i++) {
      std::cout << _calotowers[i].getEta() << " " <<
              _calotowers[i].getPhi() << " " <<
              _calotowers[i].getEem() << " " <<
              _calotowers[i].getEhad() << std::endl;
  }
}

void D_CaloTowerList::smearTowers(const D_CaloList& list_of_calorimeters) {
   for (unsigned int i=0; i<_calotowers.size(); i++) {
         // first, finds in which calo the tower is, in order to have its resolution em /had
         D_CaloElement currentCalo = list_of_calorimeters.getElement(_calotowers[i].getEta() );
         // then, smears separately the em and had parts
         float sEem  = currentCalo.getElectromagneticResolution().Smear( _calotowers[i].getEem() );
         float sEhad = currentCalo.getHadronicResolution().Smear( _calotowers[i].getEhad() );
         // finally, sets the tower smeared energy
         _calotowers[i].Set_Eem_Ehad_E_ET(sEem,sEhad);
   }
}


const bool D_CaloTower::operator>(const D_CaloTower& tocomp) const {
        if(_eta == tocomp._eta) {
                if(_phi > tocomp._phi) return true;
                else return false;
        }
        else if(_eta > tocomp._eta)  return true;
             else  return false;
};

const bool D_CaloTower::operator<(const D_CaloTower& tocomp) const {
        if(_eta == tocomp._eta) {
                if(_phi < tocomp._phi) return true;
                else return false;
        }
        else if(_eta < tocomp._eta)  return true;
             else  return false;
};

const D_CaloTower& D_CaloTowerList::getElement(const float eta, const float phi) {
  if(!_sorted) sortElements();
  for (unsigned int i=0; i< _calotowers.size(); i++) {
	if(_calotowers[i].getEta() == eta &&
	   _calotowers[i].getPhi() == phi )
		return _calotowers[i];
  }
  return dummyTower; // if not found
}

const D_CaloTower& D_CaloTowerList::getElement(const float eta, const float phi) const {
  if(!_sorted) cout << "Warning: unsorted list!";
  for (unsigned int i=0; i< _calotowers.size(); i++) {
        if(_calotowers[i].getEta() == eta &&
           _calotowers[i].getPhi() == phi )
                return _calotowers[i];
  }
  return dummyTower; // if not found
}