ClusterKtStrategy.cxx

Go to the documentation of this file.

// ================================================
//        ClusterKtStrategy class Implementation
//
//  + PreCluster helper class
//
// ================================================
//
// Namespace Atlfast
//
#include "AtlfastAlgs/ClusterKtStrategy.h"
#include "AtlfastEvent/Cell.h"
#include "AtlfastEvent/KtCluster.h"
#include "AtlfastEvent/ICluster.h"
#include "AtlfastEvent/CollectionDefs.h"

#include "KtJet/KtLorentzVector.h"
#include "KtJet/KtEvent.h"

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

//Other
#include "CLHEP/Vector/LorentzVector.h"
#include <map>
#include <vector>


namespace Atlfast {

  ClusterKtStrategy::ClusterKtStrategy(double minet, double rp, std::string angle, std::string recom, int subjet) : 
    m_minClusterET(minet), m_rParameter(rp), m_subjet(subjet){
    //set kt angle scheme
    if(angle == "angular") {m_angle = 1;}
    else if (angle == "deltaR") {m_angle = 2;}
    else if (angle == "QCD") {m_angle = 3;}
    else{
      std::cerr << "Warning, unknown Kt angle scheme. Scheme set to \"angular\"" << std::endl;
      m_angle = 1;
    }
    //set kt rcombination scheme
    if(recom == "E") {m_recom = 1;}
    else if (recom == "pt") {m_recom = 2;}
    else if (recom == "pt2") {m_recom = 3;}
    else if (recom == "et") {m_recom = 4;}
    else if (recom == "et2") {m_recom = 5;}
    else{
      std::cerr << "Warning, unknown Kt recombination scheme. Scheme set to \"E\"" << std::endl;
      m_recom = 1;
    }
  }

  void 
  ClusterKtStrategy::makeClusters( 
                                  MsgStream& log, 
                                  const std::vector<IKinematic*>& storedCells,
                                  IKinematicVector& unusedCells,
                                  IClusterCollection* clusters) const {
    log << MSG::DEBUG << storedCells.size() << " cells to start with"<<endreq;
    
    std::vector<IKinematic*>::const_iterator itr = storedCells.begin();
    std::vector<KtJet::KtLorentzVector> jetvec;
    std::map<KtJet::KtLorentzVector,IKinematic*> CellMap;
    for(; itr != storedCells.end() ; ++itr){
      KtJet::KtLorentzVector tempVec = KtJet::KtLorentzVector((*itr)->momentum());
      jetvec.push_back(tempVec);
      CellMap[tempVec] = (*itr);
    }
    log << MSG::DEBUG << "Made jetvec"<<endreq;
    KtJet::KtEvent ev(jetvec,4,m_angle,m_recom,m_rParameter);
    log << MSG::DEBUG << "Made KtEvent"<<endreq;

    vector<KtJet::KtLorentzVector> jets0 = ev.getJets();
    log << MSG::DEBUG << "Got events jets "<< jets0.size() << endreq;
    std::vector<KtJet::KtLorentzVector>::const_iterator iter = jets0.begin();
    for (; iter != jets0.end() ; ++iter){
      if((*iter).et() > m_minClusterET) { 
        std::vector<KtJet::KtLorentzVector> tompVec;
        tompVec = (*iter).copyConstituents();
        //Make a vector of associated cells
        std::vector<IKinematic*> tempCellVec;
        std::vector<KtJet::KtLorentzVector>::const_iterator iter2  = tompVec.begin();
        for (; iter2 != tompVec.end() ; ++iter2){
          if(CellMap.find((*iter2)) == CellMap.end()){
            log << MSG::ERROR <<  "Cannot find this key" << endreq;
          }else {
            tempCellVec.push_back(CellMap[(*iter2)]);
            CellMap.erase((*iter2));
            //log << MSG::DEBUG << "Removed cell from map" << endreq;
          }
        }
        //now do sub-jet analysis if required
        double yCut = 0.;
        if(m_subjet){
          KtJet::KtEvent subev(*iter,m_angle,m_recom);
          yCut = subev.getYMerge(m_subjet); 
        }
        //Make a Atlfast::KtCluster
        ICluster* clu = new KtCluster(*iter,tempCellVec.begin(),tempCellVec.end(),yCut);
        log << MSG::DEBUG << "Made a KtCluster" << endreq;
        clusters->push_back(clu);
      }
    }
    //return;
    log << MSG::DEBUG 
        <<  "Number of clusters Made by Kt algo " << clusters->size() 
        << endreq;
    std::map<KtJet::KtLorentzVector,IKinematic*>::iterator celit = 
      CellMap.begin();

    for (; celit != CellMap.end() ; ++celit){

      unusedCells.push_back( ( (*celit).second ) );

      log << MSG::DEBUG 
          << "Size of unusedCells vector " << unusedCells.size() 
          << endreq;
    }

    return;
  }
  
} // end of namespace bracket

---