CalSection.cxx

Go to the documentation of this file.

//                                                                      //
// PS   CalSection.                                                     //
// 00/09/06                                                             //
//                                                                      //
// build an array of cells spanning the eta and phi ranges passed in the//
// constructor. The granularities used are the closest to the requested //
// values, but which divide exactly into the requested ranges.          //
//                                                                      //
//                                                                      //
//                                                                      //
//                                                                      //
//                                                                      //


#include "AtlfastAlgs/CalSection.h"
#include <iostream>
#include <assert.h>
#include <algorithm>
#include "AtlfastEvent/TwoCptCell.h"
#include "AtlfastAlgs/ICellSelector.h"
#include "AtlfastAlgs/ISmearer.h"
#include "AtlfastEvent/ITransportedParticle.h"
#include "AtlfastAlgs/DeleteObject.h"
#include "GaudiKernel/MsgStream.h"
#include "AtlfastEvent/EPileupDeposit.h"
#include "FastShowerUtils/Gridlet.h"
#include "CLHEP/Vector/LorentzVector.h"
namespace Atlfast {
  using FastShower::GridletElement;

//____________________________________________________________________
  CalSection::CalSection(
                         MsgStream& log,
                         double minEta,
                         double maxEta,
                         double granEta,
                         double granPhi,
                         double minPhi,
                         double maxPhi
                         ):
    m_log(log),
    m_minEta(minEta), m_maxEta(maxEta),
    m_minPhi(minPhi), m_maxPhi(maxPhi),
    m_calSectionReject(minEta, maxEta, minPhi, maxPhi){
    
    log <<MSG::DEBUG<< "CalSection: construction. Eta range: "
        <<minEta<<" - "<<maxEta<<endreq; 
    log <<MSG::DEBUG<< "                          Phi range: "
        <<minPhi<<" - "<<maxPhi<<endreq; 
    log <<MSG::DEBUG<< " starting granularity Eta: "<<granEta
        << " phi: "<<granPhi
        <<endreq;

    log <<MSG::DEBUG<<"CalSectionReject parameters:"<<endreq;
    log <<MSG::DEBUG<<m_calSectionReject<<endreq;

    int nEta    = int(((maxEta-minEta)/granEta)+0.5);
    m_granEta   = (maxEta-minEta)/int(nEta);
    m_nPhi      = 1+int(((m_maxPhi-m_minPhi)/granPhi)+0.5);
    m_granPhi   = (m_maxPhi-m_minPhi)/double(m_nPhi); 
    log <<MSG::DEBUG<< " final granularity Eta: "<<m_granEta
        << " phi: "<<m_granPhi
        <<endreq; 
    log <<MSG::DEBUG<< " final Bins in Eta: "<<nEta
        << " phi: "<<m_nPhi
        <<endreq; 

    for (int ieta=0; ieta<nEta; ieta++){
      for (int iphi=0; iphi<m_nPhi; iphi++){
        
        double eta        = minEta+((double(ieta)+0.5)*double(m_granEta));
        double phi        = m_minPhi + ((double(iphi)+0.5)*double(m_granPhi));
        CellDescriptor id = CellDescriptor(eta,phi);
        ITwoCptCell* cell = new TwoCptCell(id);
        
        m_cells.push_back(cell);
      }
    }
  }
  CalSection::~CalSection(){
    m_log<<MSG::DEBUG<<"CalSection destructor starts "<<endreq;
    std::for_each(m_cells.begin(), m_cells.end(), DeleteObject());
    m_log<<MSG::DEBUG<<"CalSection destructor ends "<<endreq;
  }
  CalSection::CalSection(const CalSection& c):
    m_log(c.m_log),
    m_minEta(c.m_minEta),
    m_maxEta(c.m_maxEta),
    m_minPhi(c.m_minPhi),
    m_maxPhi(c.m_maxPhi),
    m_granEta(c.m_granEta),
    m_granPhi(c.m_granPhi),
    m_nPhi(c.m_nPhi),
    m_calSectionReject(c.m_calSectionReject){
    std::vector<ITwoCptCell*>::const_iterator i=c.m_cells.begin();
    for(;i<c.m_cells.end();++i) m_cells.push_back((*i)->cloneITCC());
  }
  //----------------------------------------------------------------
  void CalSection::deposit(ITransportedParticleCollectionIter& f, 
                           ITransportedParticleCollectionIter& l){
    // remove particles not in acceptance
    ITransportedParticleCollectionIter divider;
    divider = std::partition(f, l, m_calSectionReject);
    
    // fill cells 
    for(ITransportedParticleCollectionCIter i=divider; i<l; ++i){newHit(*i);}

    // assume sections do not overlap: do not reprocess.
    l=divider;
  }
  //----------------------------------------------------------------
  void 
  CalSection::depositEcal(std::vector<const GridletElement*>::iterator& f, 
                          std::vector<const GridletElement*>::iterator& l
                          ){
    //    cerr<<"Calsection depositEcal begin"<<endl;
  
    //    cerr<<"  no of elements on entry "<<l-f<<endl;
    std::vector<const GridletElement*>::iterator divider;
    divider = std::partition(f, l, m_calSectionReject);
    //    cerr<<"Calsection depositEcal after partition"<<endl;
  
    std::vector<const GridletElement*>::iterator i;
    for( i=divider; i<l; ++i){newEHit(*i);}
    //    cerr<<"Calsection depositEcal after newhit call"<<endl;
  

    // assume sections do not overlap: do not reprocess.
    l=divider;
    //    cerr<<"  no of elements on exit "<<l-f<<endl;
    //    cerr<<"Calsection depositEcal end"<<endl;
  }
  //----------------------------------------------------------------
  void CalSection::depositHcal(std::vector<const GridletElement*>::iterator& f, 
                               std::vector<const GridletElement*>::iterator& l
                               ){
    
    std::vector<const GridletElement*>::iterator divider;
    divider = std::partition(f, l, m_calSectionReject);
  
    std::vector<const GridletElement*>::iterator i;
    for( i=divider; i<l; ++i){newHHit(*i);}

    // assume sections do not overlap: do not reprocess.
    l=divider;
  }
  //----------------------------------------------------------------
  void CalSection::depositEgen(std::vector<const Gridlet*>::iterator& f, 
                               std::vector<const Gridlet*>::iterator& l
                               ){
    std::vector<const Gridlet*>::iterator divider;
    divider = std::partition(f, l, m_calSectionReject);
    
    
    std::vector<const Gridlet*>::iterator i;
    for( i=divider; i<l; ++i){setEgen(*i);}
    
    // assume sections do not overlap: do not reprocess.
    l=divider;
  }
  
  //----------------------------------------------------------------
  void CalSection::newHit( const ITransportedParticle* tp){
    //check cell is reasonable
    const HepMC::GenParticle* particle = tp->particle();
    //const double eta = particle->momentum().pseudoRapidity();
    const double eta = tp->eta();
    const double phi = tp->phi();
    int ind = iindex(phi, eta);
    ITwoCptCell* cell = m_cells[ind];
    assert(std::abs( eta-(cell->eta()) ) <m_granEta);
    assert(std::abs( phi-(cell->phi()) ) <m_granPhi);
    cell->newHit(particle);
    //use a map to store pointers to hit cells to avoid
    //duplicate entries if the cell is hit more than once.
    m_hitCells[ind] = cell;

    return;
  }
  //----------------------------------------------------------------
  void CalSection::newEHit( const GridletElement* ge){
    int ind = index(ge);
    ITwoCptCell* cell = m_cells[ind];
    //    cerr<<"gridletEl, phi, eta, eT:     "
    //    <<ge->phi()<<" "<<ge->eta()<<" "<<ge->et()<<endl;
    //    cerr<<"cell, phi, eta, granularity: "<<cell->phi()<<" "<<cell->eta()
    //  <<" "<<m_granPhi<<" "<<m_granEta<<endl;
    assert(std::abs( (ge->eta()-cell->eta()) )<m_granEta);
    assert(std::abs( (ge->phi()-cell->phi()) )<m_granPhi);
    cell->depositEcal(ge->et());
    //use a map to store pointers to hit cells to avoid
    //duplicate entries if the cell is hit more than once.
    m_hitCells[ind] = cell;

    return;
  }
  //----------------------------------------------------------------
  void CalSection::newHHit( const GridletElement* ge){
    int ind = index(ge);
    ITwoCptCell* cell = m_cells[ind];

    assert(std::abs( (ge->eta())-(cell->eta()) )<m_granEta);
    assert(std::abs( (ge->phi())-(cell->phi()) )<m_granPhi);

    cell->depositHcal(ge->et());
    //use a map to store pointers to hit cells to avoid
    //duplicate entries if the cell is hit more than once.
    m_hitCells[ind] = cell;

    return;
  }
  //----------------------------------------------------------------
  void CalSection::setEgen( const Gridlet* gr){
    int ind = index(gr);
    ITwoCptCell* cell = m_cells[ind];
    cell->addEgen(gr->eGen());
    m_hitCells[ind] = cell;
    return;
  }
  //----------------------------------------------------------------
  int  CalSection::index(const Gridlet* gr) const{
    int ind = this->iindex( gr->phi0(), gr->eta0() );
    return ind;
  }
  //----------------------------------------------------------------
  int  CalSection::index(const GridletElement* ge) const{
    int ind = this->iindex( ge->phi(), ge->eta() );
    return ind;
  }
  //----------------------------------------------------------------
  int  CalSection::index(const ITransportedParticle* tp) const {
    const HepMC::GenParticle* particle = tp->particle();
    const double eta = particle->momentum().pseudoRapidity();
    const double phi = tp->phi();
    int ind = iindex(phi, eta);
    return ind;
  }
  //----------------------------------------------------------------
  int  CalSection::index(const Atlfast::EPileupDeposit* ep) const{
    int ind = iindex(ep->phi(),ep->eta());
    return ind;
  }
  //----------------------------------------------------------------
  int CalSection::iindex(double phi, double eta) const{
    int iphi = int( (phi-m_minPhi)/m_granPhi );
    int ieta = int( (eta-m_minEta)/m_granEta );
    
    return  (m_nPhi*ieta  + iphi);
  }
  //----------------------------------------------------------------
  // FOR ENERGY PILEUP METHODS
  //----------------------------------------------------------------
  void CalSection::deposit(std::vector<Atlfast::EPileupDeposit*>::iterator& f, 
                           std::vector<Atlfast::EPileupDeposit*>::iterator& l){
    // remove pileup deposits not in acceptance
    std::vector<Atlfast::EPileupDeposit*>::iterator divider;
    divider = std::partition(f, l, m_calSectionReject);

    // fill cells 
    std::vector<Atlfast::EPileupDeposit*>::iterator  i=divider;
    for(; i<l; ++i) newHit(*i);

    // assume sections do not overlap: do not reprocess.
    // hits processed here: set last to end of unprocessed hits
    l=divider;
  }

  //----------------------------------------------------------------
  void CalSection::newHit( const Atlfast::EPileupDeposit* ep){
    //check cell is reasonable
    int ind = index(ep);
    ITwoCptCell* cell = m_cells[ind];
    assert(std::abs( (ep->eta()) - (cell->eta()) )<m_granEta);
    assert(std::abs( (ep->phi()) - (cell->phi()) )<m_granPhi);
    cell->newHit(ep);
     //use a map to store pointers to hit cells to avoid
    //duplicate entries if the cell is hit more than once.
    m_hitCells[ind] = cell;
    return;
  }
  //===========================================================================
  void CalSection::giveHits
    (const ICellSelector* p_cellselect, ITwoCptCellCollection* cells) const{
    std::map<int, ITwoCptCell*, std::less<int> >::const_iterator i;
    for(i=m_hitCells.begin(); i!=m_hitCells.end(); ++i){
      //test cell before newing and filling vector
      if( (*p_cellselect) (((*i).second))){ 
        ITwoCptCell* cell = ((*i).second)->cloneITCC();
      //Need to copy cells because Athena deletes container AND POINTERS
      // AND POINTEES
        cells->push_back(cell);
      }
    }

  }
  //===========================================================================
  void CalSection::smearCells
    (ISmearer* p_smearer) {
    std::map<int, ITwoCptCell*, std::less<int> >::const_iterator i;
    for(i=m_hitCells.begin(); i!=m_hitCells.end(); ++i){
      if((*((*i).second)).momentum().e() > 0){
        HepLorentzVector temp = p_smearer->smear((*((*i).second)).momentum());
          (*((*i).second)).setPt(temp);
      }
    }
  }



  void CalSection::reset(){
    //reset the cells of this CalSection

    std::map<int, ITwoCptCell*>::const_iterator i = m_hitCells.begin();
    for(;i!=m_hitCells.end();++i) (*i).second->resetCell();
    m_hitCells.erase(m_hitCells.begin(), m_hitCells.end());
  }
}

    
    
  

---