// ================================================
 //        ClusterFastJetStrategy class Implementation
 //
 //
 // ================================================
 //
 // Namespace Atlfast
 //
 #include "AtlfastAlgs/ClusterFastJetStrategy.h"
 #include "AtlfastEvent/Cell.h"
 #include "AtlfastEvent/KtCluster.h"
 #include "AtlfastEvent/ICluster.h"
 
 #include "fastjet/SISConePlugin.hh"
 #include "JetRec/YSplitter.h"
 
 // Gaudi includes
 #include "GaudiKernel/MsgStream.h"
 
 //Other
 #include "CLHEP/Vector/LorentzVector.h"
 
 namespace Atlfast {
 
   /**
    * This Class builds Atlfast::KtClusters from Atlfast::Cells
    * It uses the FastKt Algorithm as well as the SISCone Algorithm via the
    * FastJet Interface.
    * It stores the Atlfast::KtClusters as Atlfast::Cluster* in a std::vector.
    * All unused cells are put in the a std::vector of Atlfast::Cells
    *
    * @author Jon Couchman (Original)
    * @author Alexander Richards (FastJet update)
    */
 
   ClusterFastJetStrategy::ClusterFastJetStrategy(std::string strategy, 
                                                  double coneRadius,
                                                  double overlapThreshold,
                                                  double minet, 
                                                  double rp, 
                                                  std::string angle, 
                                                  std::string recom, 
                                                  int subjet) : 
     m_minClusterET(minet), m_rParameter(rp), m_subjet(subjet), m_plugin(0), m_jetDef(0) {
 
     
     if(strategy == "Kt" || strategy == "Cam" || strategy == "AntiKt"){
       
       //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 \"deltaR\"" << std::endl;
         m_angle = 2;
       }
       
       fastjet::RecombinationScheme recombScheme;
       
       //set kt rcombination scheme
       if(recom == "E")         {m_recom = 1; recombScheme = fastjet::E_scheme;}
       else if (recom == "pt")  {m_recom = 2; recombScheme = fastjet::pt_scheme;}
       else if (recom == "pt2") {m_recom = 3; recombScheme = fastjet::pt2_scheme;}
       else if (recom == "et")  {m_recom = 4; recombScheme = fastjet::Et_scheme;}
       else if (recom == "et2") {m_recom = 5; recombScheme = fastjet::Et2_scheme;}
       else{
         std::cerr << "Warning, unknown Kt recombination scheme. Scheme set to \"E\"" << std::endl;
         m_recom = 1;
         recombScheme = fastjet::E_scheme;
       }
       
       // What is this?
       fastjet::Strategy ktstrategy = fastjet::Best;
       fastjet::JetAlgorithm jetalgorithm = fastjet::kt_algorithm;
       if ( strategy == "Cam" )     jetalgorithm = fastjet::cambridge_algorithm;
       if ( strategy == "AntiKt" )  jetalgorithm = fastjet::antikt_algorithm;
       
       m_jetDef = new fastjet::JetDefinition(jetalgorithm, m_rParameter, recombScheme, ktstrategy);
       
     }
     else if (strategy == "Cone"){
       // Does this have to be a new statement?
       m_plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold);
       m_jetDef = new fastjet::JetDefinition(m_plugin);
       
     } 
   }
 
   std::string ClusterFastJetStrategy::Details(const char type[], const fastjet::PseudoJet &pj)const{
     char outputDetails[256];
     std::string outString;
     sprintf(outputDetails,"%10g %10g %10g %10g  %s\0",pj.px(),pj.py(),pj.pz(),pj.E(),type);
     outString.append(outputDetails);
     return outString;
   }
   std::string ClusterFastJetStrategy::Details(const char type[])const{
     char outputDetails[256];
     std::string outString;
     sprintf(outputDetails,"%10s %10s %10s %10s  %s\0","px","py","pz","E",type);
     outString.append(outputDetails);
     return outString;
   }
 
    // Destructor
 
   ClusterFastJetStrategy::~ClusterFastJetStrategy()
   {
    delete m_plugin;
    delete m_jetDef;
   } 
 
   void ClusterFastJetStrategy::makeClusters(MsgStream& log, 
                                             const std::vector<IKinematic*>& storedCells,
                                             IKinematicVector& unusedCells,
                                             IClusterCollection* clusters) const {
     
     log << MSG::DEBUG << storedCells.size() << " cells to start with"<<endreq;
 
     /** First we make a std::vector<PseudoJet> from the Atlfast::Cells
         We need to fill a map relating PseudoJets to Cells so that we
         can associate the Atlfast::Cells with the Atlfast::KtClusters    */
     
     std::vector<IKinematic*>::const_iterator itr = storedCells.begin();
     std::vector<fastjet::PseudoJet> input_particles;
     std::map<fastjet::PseudoJet,IKinematic*,JetComp> CellMap;
 
 
     log << MSG::VERBOSE <<"Event Summary..."<<endreq;
     log << MSG::VERBOSE <<Details("Input_Particles")<<endreq;
 
     fastjet::PseudoJet tempVec;
     input_particles.reserve(storedCells.size());
     for(; itr != storedCells.end() ; ++itr){
       tempVec = fastjet::PseudoJet((*itr)->momentum());
       log << MSG::VERBOSE <<Details("Input_Particles",tempVec)<<endreq;
       input_particles.push_back(tempVec);
       CellMap[tempVec] = (*itr);
     }
     log << MSG::DEBUG << "Made vector of input particles"<<endreq;
     log << MSG::VERBOSE << "Number of input Cells "<<CellMap.size()<<". Number of input particles "<<input_particles.size()<<endreq;
     int iRemoved=0; //count the removed cells
     /** Now make a ClusterSequence object from the std::vector<PsuedoJet>
         This Object calculates all the Clusters using the chosen (Kt/SISCone) clustering scheme
         defined by the input variables.
     */
     if(!m_jetDef){ 
       printf("ClusterFastJetStrategy  ERROR  Error no jet definition exists, maybe you didn't type the clustering strategy correctly.\n");
       printf("ClusterFastJetStrategy  ERROR  Returning having made no clusters.\n");
       return; 
     }
     fastjet::ClusterSequence clust_seq(input_particles, *m_jetDef);
     log << MSG::DEBUG << "Made ClusterSequence"<<endreq;
     
     /** Now we retrieve all the Clusters made by the ClusterSequence as inclusive jet vector<PseudoJet> */
     
     
     std::vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets();
     
     
     /** We can construct Atlfast::KtClusters from the returned PseudoJets
      as they are HepLorentzVectors. Also use the map to get list of Clusters associated
     Cells. Also make a second ClusterSequence from each PseudoJets constituents to calculate the
     y-cut */
     log << MSG::DEBUG << "Got events jets "<< inclusive_jets.size() << endreq;
     
     std::vector<fastjet::PseudoJet>::const_iterator ijet = inclusive_jets.begin();
     double et=0;
     double modP=0;
     
     log << MSG::VERBOSE <<Details("Output_Jets")<<endreq;
     log << MSG::VERBOSE <<Details("Output_Constituents")<<endreq;
     for (; ijet != inclusive_jets.end() ; ++ijet){
       log << MSG::VERBOSE <<Details("Output_Jets",(*ijet))<<endreq;
       modP=sqrt( (*ijet).px()*(*ijet).px() + (*ijet).py()*(*ijet).py() + (*ijet).pz()*(*ijet).pz() );
       et=(*ijet).E() * (*ijet).perp()/modP;
       
       if(et > m_minClusterET) { 
         std::vector<fastjet::PseudoJet> jet_constituents;
         jet_constituents = clust_seq.constituents(*ijet);
         //Make a vector of associated cells
         std::vector<IKinematic*> tempCellVec;
         std::vector<fastjet::PseudoJet>::const_iterator iconstituent  = jet_constituents.begin();
         std::vector<HepLorentzVector> ysplit_constituents;
         log << MSG::VERBOSE <<"Number of constituents in jet: "<<jet_constituents.size()<<endreq;
         for (; iconstituent != jet_constituents.end() ; ++iconstituent){
           log << MSG::VERBOSE <<Details("Output_Constituents",(*iconstituent))<<endreq;
           ysplit_constituents.push_back(HepLorentzVector((*iconstituent).px(),(*iconstituent).py(),(*iconstituent).pz(),(*iconstituent).E()));
           if(CellMap.find(*iconstituent) == CellMap.end()){
             log << MSG::WARNING <<  "Cannot find this constituent in map of cells" << endreq;
           }
           else {
             tempCellVec.push_back(CellMap[(*iconstituent)]);
             log << MSG::VERBOSE <<"Removing cell as matched with constituent"<<endreq;
             CellMap.erase(*iconstituent);
             ++iRemoved;
             //log << MSG::DEBUG << "Removed cell from map" << endreq;
           }
         }
         //now do sub-jet analysis if required
         double yCut = 0.;
         //int type = 4 // pp collisions
         if(m_subjet){
           YSplitter ys(ysplit_constituents,/*type*/4, m_angle, m_recom, m_rParameter);//type = 4 - pp collisions
           yCut = ys.getYMergeVals(1).at(0);
         }
 
         //Make a Atlfast::KtCluster
         HepLorentzVector tmp((*ijet).px(), (*ijet).py(), (*ijet).pz(), (*ijet).E());
         ICluster* clu = new KtCluster(tmp,tempCellVec.begin(),tempCellVec.end(),yCut);
         log << MSG::DEBUG << "Made a KtCluster" << endreq;
         /** put cluster into main vector */
         clusters->push_back(clu);
       }
       else{ log << MSG::VERBOSE <<"Jet has not enough et to cluster: JetEt = "<<et<<", min cluster Et = "<<m_minClusterET<<endreq; }
     }
     //return;
     log << MSG::DEBUG <<  "Number of clusters Made by Kt algo " << clusters->size() << endreq;
     /** Put copies of all the unassociated cells into the unusedCells vector */
     std::map<fastjet::PseudoJet,IKinematic*,JetComp>::iterator celit = CellMap.begin();
 
     for (; celit != CellMap.end() ; ++celit){
 
       unusedCells.push_back( ( (*celit).second ) );
     }
     
     log << MSG::VERBOSE << "Number of output clusters "<<inclusive_jets.size()<<".Number of Removed cells "<<iRemoved<<". Number of saved KtClusters "<<clusters->size()<< endreq;
     log << MSG::DEBUG << "Size of unusedCells vector " << unusedCells.size() << endreq;
     return;
   }
   
 } // end of namespace bracket
 
 
 
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.