---

XANAClusters.cpp

Go to the documentation of this file.

#include <Utilities/Configuration/interface/Architecture.h>

#include <XANADOO/XANAUtilities/interface/XANAToolkit.h>
#include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
#include <XANADOO/XANAClusters/interface/XANAHadCluster.h>
#include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
#include <XANADOO/XANAClusters/interface/XANACellID.h>

#include <ElectronPhoton/EgammaClusters/interface/EgammaBasicCluster.h>
#include <ElectronPhoton/EgammaClusters/interface/EgammaBasicDynamic.h>
#include <ElectronPhoton/EgammaClusters/interface/EgammaBasicIsland.h>
#include <ElectronPhoton/EgammaClusters/interface/EgammaBasicHybrid.h>
#include <ElectronPhoton/EgammaClusters/interface/EGBLogPosCorr.h>
#include <ElectronPhoton/ClusterTools/interface/IsoCalculator.h>
#include <ElectronPhoton/EgammaAnalysis/interface/EgammaHAnalyzer.h>
#include <ElectronPhoton/EgammaAnalysis/interface/EgammaAlgorithm.h>

#include <Calorimetry/CaloRecHit/interface/CaloRecHit.h>

#include <CARF/G3Event/interface/G3EventProxy.h>
#include <CARF/Reco/interface/AutoRecCollection.h>

#include <Utilities/Notification/interface/Observer.h>
#include <Utilities/Notification/interface/PackageInitializer.h>
#include <Utilities/UI/interface/SimpleConfigurable.h> 
#include <Utilities/Notification/interface/Singleton.h>

#include <vector>

struct ClusterObserver : public Observer<G3EventProxy *> {

  int nbClusters_ ;
  double energySum_ ;
  SimpleConfigurable<bool> writeCaloRecHits_;
  SimpleConfigurable<bool> produceReference_;
  SimpleConfigurable<string> clusterizer_;

  ClusterObserver()
   : nbClusters_(0), energySum_(0)
   { 
    // initialize observer
    Observer<G3EventProxy *>::init() ;      
    
    // steering
    writeCaloRecHits_=SimpleConfigurable<bool>(false,"XANAClusters:ESDWithCaloRecHits");
    produceReference_=SimpleConfigurable<bool>(false,"XANAClusters:produceReference");
    clusterizer_=SimpleConfigurable<string>("EGBCluster","XANAClusters:Clusterizer");
   }
   
  ~ClusterObserver()
   {
    std::cout<<"[OVAL] number of clusters : "<<nbClusters_<<endl ;
    std::cout<<"[OVAL] mean energy : "<<energySum_/nbClusters_<<endl ;
   }

   
  void upDate( G3EventProxy * ev )
   { 
    
    if (!ev) return ;    
    
    RecCollection<EgammaBasicCluster> clusters ( RecQuery(clusterizer_.value()) );
    RecCollection<EgammaBasicCluster> ::const_iterator cluster ;
    
    // to produce reference
    if (produceReference_) {
      // BasicClusters
      for (cluster=clusters.begin(); cluster!=clusters.end(); cluster++)  
       {
        // remove double reconstruction in case of EGBCluster
        if (clusterizer_.value() == "EGBCluster") {
          // hybrid is algo>=100 and is used in barrel only
          // island is algo 1 and is used everywhere
          if ( (cluster->whichAlgo()>99 && fabs(cluster->eta())>1.4791) || 
             (cluster->whichAlgo()<99 && fabs(cluster->eta())<1.4791) ) continue;
        }
        nbClusters_++ ;
        energySum_ += cluster->energy() ;
        std::cout<<"[OVAL] new cluster with energy "<<cluster->energy()<<flush ;
        std::cout<<" and "<<cluster->numberOfDigis()<<" crystal(s)"<<endl ;
       }
       return ;
    } 
    
    // construct XANAEmClusters from BasicClusters
    vector<XANACluster> xclusters;
    XANACluster *xcluster;
    for (cluster=clusters.begin(); cluster!=clusters.end(); cluster++)  
     {
      // remove double reconstruction in case of EGBCluster
      if (clusterizer_.value() == "EGBCluster") {
        // hybrid is algo>=100 and is used in barrel only
        // island is algo 1 and is used everywhere
        if ( (cluster->whichAlgo()>99 && fabs(cluster->eta())>1.4791) || 
           (cluster->whichAlgo()<99 && fabs(cluster->eta())<1.4791) ) continue;
      }
      xana::out(2) << "[update] new cluster with energy " << cluster->energy() << " from algo " << cluster->whichAlgo() << std::endl;
      nbClusters_++ ;
      // for naming, stick to the RecAlgo names
      std::string name = "";
      if (cluster->whichAlgo() == 1) name="EGBCluster";
      else if (cluster->whichAlgo() == 2) name="EGBDynamicCluster";
      else if (cluster->whichAlgo() == 100) name="EGBCluster";
      xcluster = new XANAEmCluster(cluster->energy(), cluster->position(), cluster->chi2(), name); 
      if (writeCaloRecHits_) {
        std::vector<const CaloRecHit *> rhits=cluster->rHits();
        for (vector<const CaloRecHit *>::const_iterator
         it=rhits.begin(); it != rhits.end(); it++) {
          xcluster->addHit(new XANACaloRecHit((*it)->energy(), (*it)->time(),
          (*it)->chi2(), XANACellID((*it)->getMyCell().position())));
//        (*it)->chi2(), new XANACellID((*it)->getMyCell().position())));
        } 
      } else 
        xcluster->setNumberOfRecHits(cluster->numberOfDigis());
      xclusters.push_back(*xcluster);
     }
     
    // read XANAEmClusters 
    vector<XANACluster>::iterator it;
    for (it=xclusters.begin(); it!=xclusters.end(); it++)  
     {
      std::cout<<"[OVAL] new cluster with energy "<<it->getEnergy()<<flush ;
      std::cout<<" and position "<<it->getPosition()<<endl ;
     }
   }
      
} ;
 
static PKBuilder<ClusterObserver> clusterObserver("ClusterObserver") ;

---