CellMaker.cxx

Go to the documentation of this file.

// ================================================
// CellMaker 
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace Atlfast
//

#include "AtlfastAlgs/CellMaker.h"
#include "AtlfastAlgs/CellsAboveThreshold.h"
#include "AtlfastEvent/ITransportedParticleFactory.h"
#include "AtlfastAlgs/Calorimeter.h"
#include "AtlfastAlgs/ICellSelector.h"
#include "AtlfastAlgs/ISmearer.h"
#include "AtlfastAlgs/CellSmearer.h"
#include "AtlfastAlgs/GlobalEventData.h"
#include "AtlfastAlgs/EPileupMap.h"
#include "AtlfastAlgs/containerDelete.h"
#include "AtlfastEvent/CollectionDefs.h"

#include "AtlfastUtils/TesIO.h"
#include "AtlfastUtils/HeaderPrinter.h"

#include <cmath>
#include <iomanip>
#include <vector>
// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/Chrono.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"

#include "CLHEP/Units/SystemOfUnits.h"
#include "NavFourMom/INavigable4MomentumCollection.h"

namespace Atlfast {
  CellMaker:: CellMaker(const std::string& name, ISvcLocator* pSvcLocator):
    Algorithm(name,pSvcLocator),
    m_calorimeter(0),
    m_cellSelector(0),
    m_mcSelector(0),
    m_epileupMap(0),
    m_smearer(0),
    m_tesIO(0),
    m_log( messageService(), name )
  {
    
    //default paths for entities in the TES
    m_outputLocation      = "/Event/AtlfastCells" ;
    m_outputLocation_IN4M = "/Event/AtlfastCellsIN4M";

    // default parameters to build the calorimeter
    m_etaCoverage        = 5.0;
    m_minETCell          = 0.0*GeV;
    m_granBarrelEta      = 0.1;
    m_granBarrelPhi      = 0.1;
    m_granForwardEta     = 0.2;
    m_granForwardPhi     = 0.2;
    m_doSmearing         = false;
    m_fastShower         = false;
    
    declareProperty("OutputLocation",     m_outputLocation);
    declareProperty("OutputLocationIN4M", m_outputLocation_IN4M);
    declareProperty("EtaCoverage",        m_etaCoverage);
    declareProperty("MinETCell",          m_minETCell);
    declareProperty("GranBarrelEta",      m_granBarrelEta);
    declareProperty("GranBarrelPhi",      m_granBarrelPhi);
    declareProperty("GranForwardEta",     m_granForwardEta);
    declareProperty("GranForwardPhi",     m_granForwardPhi);
    declareProperty("DoSmearing",     m_doSmearing);
    declareProperty("FastShower",     m_fastShower);
    
  }
  
  //__________________________________________________________________________
  CellMaker::~CellMaker(){

    m_log<<MSG::DEBUG<<"CellMaker destructor starts "<<endreq;

    m_log<<MSG::DEBUG<<"Deleting Calorimeter "<<endreq;
    if (m_calorimeter) delete m_calorimeter;

    m_log<<MSG::DEBUG<<"Deleting TesIO "<<endreq;
    if (m_tesIO) delete m_tesIO;

    m_log<<MSG::DEBUG<<"Deleting Cell Selector "<<endreq;
    if (m_cellSelector) delete m_cellSelector;

    if (m_epileupMap) delete m_epileupMap;
    if (m_smearer) delete m_smearer;

    m_log<<MSG::DEBUG<<"CellMaker destructor ends "<<endreq;
  }
  //__________________________________________________________________________
  StatusCode CellMaker::initialize()
  {
    

    //set up objects used to determine particle charge,
    // and to select input MC 4-vectors.
    m_cellSelector  = new CellsAboveThreshold(m_minETCell);
    

    //get the Global Event Data using singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    m_fieldOn          = ged -> fieldOn();
    m_mcSelector       = ged -> visibleToCal();
    m_barrelForwardEta = ged -> barrelForwardEta();
    m_mcLocation       = ged -> mcLocation();
    m_monopoleIDs      = ged -> monopoleIDs();
 
    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());


    if (fabs(m_etaCoverage) < fabs(m_barrelForwardEta) ) {
      m_log << MSG::ERROR 
          << "asked to makecalo with max extent " 
          << m_etaCoverage 
          << "but barrel-forward boundary at " 
          << m_barrelForwardEta << endreq;
      return StatusCode::FAILURE;
    }
    m_calorimeter = new Calorimeter(
                                    m_log,
                                    m_fastShower,
                                    m_etaCoverage, 
                                    m_barrelForwardEta,
                                    m_granBarrelEta, 
                                    m_granBarrelPhi, 
                                    m_granForwardEta, 
                                    m_granForwardPhi
                                    );
    if(m_doSmearing == true){
      m_epileupMap = new EPileupMap(ged->lumi());
      m_smearer = new CellSmearer(ged->randSeed(),m_barrelForwardEta);
    }else{
      m_epileupMap = NULL;
      m_smearer = NULL;
    }

    HeaderPrinter hp("Atlfast Cell Maker:", m_log);
    
    hp.add("Min Cell ET                         ", m_minETCell);  
    hp.add("Magnetic Field is on?               ", m_fieldOn);
    hp.add("MC Input Location                   ", m_mcLocation);
    hp.add("Output Location                     ", m_outputLocation);
    hp.add("INavigable4Momentum Output Location ", m_outputLocation_IN4M);
    hp.add("DoSmearing                          ", m_doSmearing);    
    hp.add("Run FastShower                      ", m_fastShower);    
    hp.add("Calorimeter Geometry:               ");
    hp.add("Total Eta Range                     ", m_etaCoverage);
    hp.add("Barrel Eta Range                    ", m_barrelForwardEta);
    hp.add("Barrel Eta Granularity              ", m_granBarrelEta);
    hp.add("Barrel Phi Granularity              ", m_granBarrelPhi);
    hp.add("Forward Eta Granularity             ", m_granForwardEta);
    hp.add("Forward Phi Granularity             ", m_granForwardPhi);
    hp.print();


      
    return StatusCode::SUCCESS; 
  }
  
  StatusCode CellMaker::finalize()
  {
    
    m_log<<MSG::INFO<<"Finalizing"<<endreq;
    return StatusCode::SUCCESS; 
  }
  
  //_________________________________________________________________________
  StatusCode CellMaker::execute()
  {
    //.............................................
    
    m_log << MSG::DEBUG << "CellMaker execute()" << endreq;

    ITransportedParticleCollection particlesAtCal;
    deflectParticles(particlesAtCal);
    
    // give to calorimeter (iseq gives the interval [begin,end)
    ITransportedParticleCollectionIter fpart = particlesAtCal.begin();
    ITransportedParticleCollectionIter lpart = particlesAtCal.end();
    m_calorimeter->deposit(fpart, lpart);
    
    //Delete the the container of transported particles and its pointers
    containerDelete(particlesAtCal.begin(), particlesAtCal.end());


    // Add in the energy pileup 
    if(m_doSmearing){

      //Smear the cells resolution first
      m_calorimeter->smearCells(m_smearer);
      m_log << MSG::DEBUG << "Doing Smearing" << endreq;

      m_epileupMap->fillMap();
      m_log << MSG::DEBUG 
          << "Size of EM EPileup vector: " 
          << m_epileupMap->getEMMap()->size() 
          << endreq;
      std::vector<Atlfast::EPileupDeposit*>::iterator fdep = 
        m_epileupMap->getEMMap()->begin();
      std::vector<Atlfast::EPileupDeposit*>::iterator ldep = 
        m_epileupMap->getEMMap()->end();

      m_calorimeter->deposit(fdep, ldep);
      m_log << MSG::DEBUG << "ECAL pileup deposited" << endreq;
      m_log << MSG::DEBUG << "Size of HAD EPileup vector: " 
          << m_epileupMap->getHadMap()->size() 
          << endreq;
      fdep = m_epileupMap->getHadMap()->begin();
      ldep = m_epileupMap->getHadMap()->end();

      m_calorimeter->deposit(fdep, ldep);
      m_log << MSG::DEBUG << "HCAL pileup deposited" << endreq;
    }

    m_log << MSG::DEBUG << "deposit" << endreq;
    
    // Ask for Cells passing selection criteria of m_cellSelector (=Et).

    // For now, we will make, fill and store two collections here
    // The first is to be used iternally by Atlfast Algorithms:
    ITwoCptCellCollection* cells = new ITwoCptCellCollection;
    // The second will be read in by the full jet reconstruction:
    INavigable4MomentumCollection *cells_IN4M = 
      new INavigable4MomentumCollection(SG::VIEW_ELEMENTS);
    

    //************************
    m_calorimeter->giveHitCells(m_cellSelector, cells);
    m_log << MSG::DEBUG << "hit cells retrieved, " << cells->size() 
        << " cells were hit this event " << endreq;

    // Fill the IN4M collection with pointers to the ITwoCptCells
    ITwoCptCellCollection::const_iterator cellItr = cells->begin();
    for (;cellItr != cells->end(); ++cellItr){
      INavigable4Momentum *cell_IN4M = *cellItr;
      cells_IN4M->push_back(cell_IN4M);
    }
    
    // finished with the calorimeter
    m_calorimeter->reset();
    m_log << MSG::DEBUG << "calorimeter reset" << endreq;
    
    // then store and return
    TesIoStat stat = m_tesIO->store(cells, m_outputLocation );  
    std::string mess = stat ? "Stored cells in TES":"Failed store cells in TES";
    m_log << MSG::DEBUG << mess << endreq;

    stat = m_tesIO->store(cells_IN4M, m_outputLocation_IN4M );
    if (stat) m_log << MSG::DEBUG << "Stored cells in TES as INavigable4MomentumCollection" << endreq;
    else m_log << MSG::ERROR << "Failed to store cells in TES as INavigable4MomentumCollection" << endreq;
    
    return stat;

  }
  
  void CellMaker::deflectParticles(ITransportedParticleCollection &itpc){
    
    m_log << MSG::DEBUG << "CellMaker deflectParticles()" << endreq;

    // read MC particles from TES
    MCparticleCollection p;
    TesIoStat stat = m_tesIO->getMC( p, m_mcSelector ) ;
    std::string mess;
    mess = stat? "Retrieved MC from TES ":"Failed MC retrieve from TES ";
    m_log << MSG::DEBUG << mess << p.size()<<endreq;
    
    MCparticleCollectionCIter ip= p.begin();
    for(; ip<p.end(); ++ip){
      
      // Factory comes back with the relevant TransportedParticle
      ITransportedParticle *itp = ITransportedParticleFactory::create(*ip,m_monopoleIDs);
      itp->deflect();
      itpc.push_back(itp);
      
    }   
    
    return;
    
  }
  
  
} // end of namespace bracket

---