source: svn/trunk/src/CaloUtil.cc@ 366

/*
  ---- Delphes ----
  A Fast Simulator for general purpose LHC detector
  S. Ovyn ~~~~ severine.ovyn@uclouvain.be
  
  Center for Particle Physics and Phenomenology (CP3)
  Universite Catholique de Louvain (UCL)
  Louvain-la-Neuve, Belgium
*/

/// \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
}