ClusterMaker.cxx

Go to the documentation of this file.

// ================================================
//        ClusterMaker class Implementation
//
//  + PreCluster helper class
//
// ================================================
//
// Namespace Atlfast
//

// This package includes
#include "AtlfastAlgs/ClusterMaker.h"
#include "AtlfastAlgs/ClusterConeStrategy.h"
#include "AtlfastAlgs/SharedConeStrategy.h"
#include "AtlfastAlgs/ClusterKtStrategy.h"
#include "AtlfastUtils/FunctionObjects.h"
#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastEvent/Cell.h"
#include "AtlfastEvent/TwoCptCell.h"
#include "AtlfastEvent/ITwoCptCell.h"
#include "AtlfastEvent/Cluster.h"
#include "AtlfastUtils/IClusterStrategy.h"
#include "AtlfastUtils/TesIO.h"
#include "AtlfastAlgs/GlobalEventData.h"
//#include "AtlfastAlgs/TypeRecoverer.h"
#include "AtlfastAlgs/ContainerDispatcher.h"
#include "AtlfastEvent/TypeVisitor.h"

#include <cmath> 
#include <algorithm>
#include <vector>

// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"

//Other
#include "CLHEP/Vector/LorentzVector.h"

namespace Atlfast {
  //--------------------------------
  // Constructors and destructors
  //--------------------------------
  
  ClusterMaker::ClusterMaker
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ),
      m_clusterStrategy(0),
      m_tesIO(0)
  {
    
    // Setting the parameter defaults.
    m_rConeBarrel         = 0.401; //+.001 avoids grid problem
    m_rConeForward        = 0.401;
    m_minInitiatorET      = 1.5*GeV;
    m_minClusterET        = 5.0*GeV;
    m_processCells        = true;
    m_processTracks       = false;
    m_inputCellLocation   = "/Event/AtlfastCells" ;
    m_inputTrackLocation  = "/Event/AtlfastTracks" ;
    m_outputLocation      = "/Event/AtlfastClusters";
    m_unusedCellLocation  = "/Event/AtlfastUnusedCells";
    m_unusedTrackLocation = "/Event/AtlfastUnusedTracks";
    m_masslessJets        = true;
    m_strategy            = "Cone";
    //=====================
    m_ktRParameter         = 1.0;
    m_ktAngle             = "deltaR";
    m_ktRecomScheme        = "E";
    m_ktYCut              = 0;
    // Declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    
    declareProperty( "RConeBarrel",         m_rConeBarrel ) ;
    declareProperty( "RConeForward",        m_rConeForward ) ;
    declareProperty( "minInitiatorET",      m_minInitiatorET ) ;
    declareProperty( "minClusterET",        m_minClusterET ) ;
    declareProperty( "Strategy",            m_strategy ) ;
    declareProperty( "ProcessCells",        m_processCells ) ;
    declareProperty( "ProcessTracks",       m_processTracks ) ;
    declareProperty( "InputCellLocation",   m_inputCellLocation ) ;
    declareProperty( "InputTrackLocation",  m_inputTrackLocation ) ;
    declareProperty( "OutputLocation",      m_outputLocation ) ;
    declareProperty( "UnusedCellLocation",  m_unusedCellLocation ) ;
    declareProperty( "UnusedTrackLocation", m_unusedTrackLocation ) ;
    declareProperty( "MasslessJets",        m_masslessJets ) ;
    //==============================================
    declareProperty( "KtRParameter",         m_ktRParameter ) ;
    declareProperty( "KtAngle",            m_ktAngle ) ;
    declareProperty( "KtRecomScheme",      m_ktRecomScheme ) ;
    declareProperty( "KtYCut",             m_ktYCut ) ;
    
  }

  // Destructor
  ClusterMaker::~ClusterMaker() {
    if(m_tesIO){ 
      delete m_tesIO;
    }
    if(m_clusterStrategy){ 
      delete m_clusterStrategy;
    }
  } 

  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode ClusterMaker::initialize(){
 
   MsgStream log( messageService(), name() ) ;
   log << MSG::DEBUG<< "Cluster Maker initialising " << endreq;
 
   //moved this before TesIO
   GlobalEventData* ged = GlobalEventData::Instance();
   m_mcLocation  = ged -> mcLocation();

   //    m_tesIO = new TesIO( eventDataService() );
   m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());

    if(m_strategy == "Cone"){
      m_clusterStrategy = new ClusterConeStrategy(
                                                  m_rConeBarrel,   
                                                  m_rConeForward,
                                                  m_minInitiatorET,
                                                  m_minClusterET,
                                                  m_masslessJets);

    } else if(m_strategy == "Kt"){
      m_clusterStrategy = new ClusterKtStrategy(m_minClusterET,m_ktRParameter,m_ktAngle,m_ktRecomScheme,m_ktYCut);
    } else if(m_strategy == "Shared"){
      double barrelForwardEta = ged -> barrelForwardEta();
      m_clusterStrategy = new SharedConeStrategy(m_rConeBarrel,   
                                                 m_rConeForward,
                                                 m_minInitiatorET,
                                                 m_minClusterET,
                                                 barrelForwardEta);
    }else{
      m_clusterStrategy = new ClusterConeStrategy(
                                                  m_rConeBarrel,   
                                                  m_rConeForward,
                                                  m_minInitiatorET,
                                                  m_minClusterET,
                                                  m_masslessJets);
      
      m_strategy = "Cone";
    }

    HeaderPrinter hp("Atlfast Cluster Maker:", log);
    hp.add("Cluster Strategy            ", m_strategy);  
    if(m_strategy == "Kt"){
      hp.add( "Kt R-Parameter              ", m_ktRParameter ) ;
      hp.add( "Kt Angle                    ", m_ktAngle ) ;
      hp.add( "Kt Recombination Scheme     ", m_ktRecomScheme ) ;
      hp.add( "Kt Y-Merge value            ", m_ktYCut ) ;
    }else{
      hp.add("Endcap Cone Size             ", m_rConeForward);
      hp.add("Barrel Cone Size             ", m_rConeBarrel);
      hp.add("Min ET for Cell initiator    ", m_minInitiatorET);
    }
    hp.add("Min ET for cluster          ", m_minClusterET);
    hp.add("Process Cells               ", m_processCells);
    hp.add("Process Tracks              ", m_processTracks);
    if(m_processCells){
      hp.add("Input CellLocation          ", m_inputCellLocation);
    }
    if(m_processTracks) {
      hp.add("Input TrackLocation         ", m_inputTrackLocation);
    }
    hp.add("Output Location             ", m_outputLocation);    
    hp.add("Unused Cell Location        ", m_unusedCellLocation);    
    hp.add("Unused Track Location       ", m_unusedTrackLocation);    
    hp.add("Massless Jets               ", m_masslessJets);    
    hp.print();

    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode ClusterMaker::finalize(){
    return StatusCode::SUCCESS ;
  }
  

  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  // This execute method is written to take care of the administration
  // only (accessing TES, etc).
  //
  // The actual clustering strategy is delegated to a further method
  // (to become an object at some point)
  //
  // All Clusters created are output to the TES
  
  StatusCode ClusterMaker::execute( ){
    MsgStream log( messageService(), name() ) ;
    std::string message;
    //.........................................................
    // Extract the elements (cells, tracks..) from the TES 
    //
    //copy to vector -cannot use DataVector fo r type change!
    std::vector<IKinematic*> elements;  

    TesIoStat stat;
    if(m_processCells){
      stat = m_tesIO->copy<ITwoCptCellCollection>(elements, 
                                                  m_inputCellLocation);
      message = stat?  "Found Cells in TES":"No Cells found in TES";
      log<<MSG::DEBUG << message <<" "<<elements.size()<<endreq;
    }
    if(m_processTracks){
      stat = m_tesIO->copy<TrackCollection>(elements, m_inputTrackLocation);
      message = stat?  "Found Tracks in TES":"No Tracks found in TES";
      log<<MSG::DEBUG << message <<" "<<elements.size()<<endreq;
    }
    // ......................
    // Make a containers to store the new Cluster pointers,
    // and unused elements
    
    IClusterCollection* clusters = new IClusterCollection ;
    IKinematicVector unusedLocalElements;

    log << MSG::DEBUG << "Starting Clustering Strategy"  << endreq ;
    
    m_clusterStrategy->makeClusters( 
                                    log, 
                                    elements, 
                                    unusedLocalElements, 
                                    clusters
                                    ) ;
    
    //......................................
    // Set up for conversions of unused IKinematics back to concrete types
    //
    log<<MSG::DEBUG
       <<"Number Of Unused Local Elements "
       << unusedLocalElements.size()<<endreq;
    //    TypeRecoverer atr(unusedLocalElements);
    TypeVisitor types = ContainerDispatcher(
                                          unusedLocalElements.begin(),
                                          unusedLocalElements.end(),
                                          TypeVisitor()
                                          );
    //log<<MSG::DEBUG <<"Recovered Types "<<endreq;


    //Register the newly made clusters
    stat = m_tesIO -> store(clusters, m_outputLocation);
    message = stat ?  "Clusters stored":"Failed to store Clusters ";
    log<<MSG::DEBUG<<message<<endreq;
    //  Register the any unused IKs if they are Cells
    stat = m_tesIO -> store(new TwoCptCellVector(types.typeVector(TwoCptCell())), 
                            m_unusedCellLocation);
    message =stat? "unused cells stored":"unused cells store failed";
    log<<MSG::DEBUG<<message<<endreq;
    //log<<MSG::DEBUG<<"Number Of Unused Cells "<< atr.cells()->size()<<endreq;
    //  Register the any unused IKs if they are Tracks
    stat = m_tesIO -> store(new TrackVector(types.typeVector(Track())), 
                            m_unusedTrackLocation);
    message =stat? "unused tracks stored":"unused cells store failed";
    log<<MSG::DEBUG<<message<<endreq;
    //log<<MSG::DEBUG<<"Number Of Unused Tracks "<< atr.tracks()->size()<<endreq;
    return stat ;
    
  }
  
} // end of namespace bracket

---