7.3 ECAL Clustering

Complete:
Detailed Review status

Goals of this page

Reconstruct BasicClusters (local maxima) and SuperClusters (clusters of BasicClusters to collect radiated energy) in the ECAL starting from Digis.

Contents

• Introduction
  • Hybrid superclusters
  • Island basic clusters and superclusters
• Procedure
• Set up your Environment
• Run the Reconstruction
  • Clusters from EcalRecHits
  • Clusters directly from Digis
  • Contents of output root file
• Fill Histograms with EgammaSimpleAnalyzer
  • Structure of analyzer
  • Configuration and parameters of analyzer
  • Book histograms in ::beginJob()
  • Do analysis in ::analyze()
  • Store output in ::endJob()
  • Output of EgammaSimpleAnalyzer
• Comparision of SuperCluster η with η of generated particle
• Access to cluster shape objects for basic clusters
• Information Sources
• Review status

Introduction

Electron and photon showers deposit their energy in several crystals in the ECAL. Approximately 94% of the incident energy of a single electron or photon is contained in 3x3 crystals, and 97% in 5x5 crystals. Summing the energy measured in such fixed arrays gives the best performance for unconverted photons, or for electrons in the test beam. The presence in CMS of material in front of the calorimeter results in bremsstrahlung and photon conversions. Because of the strong magnetic field the energy reaching the calorimeter is spread in φ. The spread energy is clustered by building a cluster of clusters, a "supercluster", which is extended in φ.

Hybrid superclusters

The η-φ geometry of the barrel crystals allows the the Hybrid algorithm to exploit our knowledge of the lateral shower shape in the η direction (by taking a fixed bar of 3 or 5 crystals in η), while searching dynamically for separated energy in the φ direction. Unlike the Island algorithm, the Hybrid algorithm cannot be used without making superclusters. The algorithm makes superclusters which are then decomposed into sub-clusters ("basic clusters").

Island basic clusters and superclusters

The Island algorithm starts with a seed crystal, local energy maximum above a defined threshold. Starting from the seed position, adjacent crystals are examined, scanning first in φ and then in η. Along each scan line, crystals are added to the cluster until a rise in energy (meaning that the area with energies at the noise level or a crystal belonging to another cluster is reached) or crystal that has not been read out (i.e. without recorded energy) is encountered.

For a crystal to be added to a cluster the following must be true:

• The crystal contains a rechit with positive energy.
• The crystal has not been assigned to another cluster already.
• The previous crystal added (in the same direction) has higher energy.

In much the same way as energy is clustered at the level of the crystals, non-overlapping Island clusters can in turn be clustered into superclusters. The procedure is seeded by searching for the most energetic cluster and then collecting all the other nearby clusters in a very narrow η-window, and much wider φ-window.

Procedure

The sequence for ECAL standalone reconstruction is illustrated below. The first step is to perform local reconstruction of EcalRecHits. This is followed by hybrid clustering (barrel only) and island clustering (barrel and endcaps). The next step is to apply an energy scale correction to set the peak of the reconstructed distribution of ET/ET_true at 1.0. The correction is a parametric function of the number of crystals in the seed BasicCluster of the SuperCluster. Addition of energy deposited in the preshower is the final step in the standalone ECAL reconstruction procedure.

Set up your Environment

Set your runtime environment (shown for release %WBRELEASENEW%):
cd CMSSW_%WBRELEASENEW%/src
cmsenv

• Full instructions
• Summary for every login session
• Access modules from the CVS repository.

Run the Reconstruction

Checkout the following code into the CMSSW_7_4_15/src directory:

cvs co -r CMSSW_1_6_0 RecoEcal/EgammaClusterProducers

The reconstruction procedure is defined configuration files under RecoEcal/EgammaClusterProducers/test. Go to this directory and update the config files:

cd RecoEcal/EgammaClusterProducers/test
cvs update -r egammaClustersTutorial160 egammaDigiToClusters.cfg egammaRecHitsToClusters.cfg

Clusters from EcalRecHits

The production team provides event collections that contain EcalRecHits, and which can be used to reconstruct BasicClusters and SuperClusters. The configuration file egammaRecHitsToClusters.cfg defines the procedure. The configuration file contains the top-level process block named eg which includes several other blocks. The event number set to -1 in the source block means that the program runs over all the events in the input data file. The components (the C++ classes to be instantiated) of the cluster reconstruction are defined in module blocks for Hybrid clustering (barrel) Island clustering (endcaps) and preshower clustering, called within the cff file:

# all clustering in one sequence
include "RecoEcal/EgammaClusterProducers/data/ecalClusteringSequence.cff"

which in turn includes the following:

# create sequence for island clustering
include "RecoEcal/EgammaClusterProducers/data/islandClusteringSequence.cff"

# hybrid clustering sequence
include "RecoEcal/EgammaClusterProducers/data/hybridClusteringSequence.cff"
 
# clusters in preshower
include "RecoEcal/EgammaClusterProducers/data/preshowerClusteringSequence.cff"

and defines:

sequence ecalClusteringSequence = {
                                     islandClusteringSequence,
                                     hybridClusteringSequence,
                                     preshowerClusteringSequence
                                   }

In the component cff files, the clustering algorithm modules are put together (sequenced for inclusion in a path) into sequence blocks:

sequence islandClusteringSequence =  {
             islandBasicClusterProducer,
             islandSuperClusterProducer,
             correctedIslandSuperClusterProducer }

sequence hybridClusteringSequence = {
           hybridSuperClusterProducer,
           correctedHybridSuperClusterProducer }

sequence preshowerClusteringSequence = { correctedEndcapSuperClustersWithPreshower }

All of these modules are run as a result of including ecalClusteringSequence in the cfg path statement:

path p = { ecalClusteringSequence }

Edit egammaRecHitsToClusters.cfg to run over your favorite collection (must contain RecHits produced with CMSSW_1_6_X).

To run it go to RecoEcal/EgammaClusterProducers/test directory and type:

cmsRun egammaRecHitsToClusters.cfg

The output of your job is the egamma_clusters.root root file which contains the reconstruced basic and super clusters (see below for more details on the content of the file).

Clusters directly from Digis

Alternatively, the configuration file egammaDigiToClusters.cfg for going from Digis to SuperClusters without storing the intermediate RecHits is also provided.

However, keep in mind that in order to use this configuration file, the input collection MUST contain the digis. Otherwise no rechits will be made and therefore no clusters. In general you don't need to use this configuration file, since the official production will provide collections that contain directly the EcalRecHits.

In addition to the blocks needed for clustering, egammaDigiToClusters.cfg contains also the block ecalLocalRecoSequence for local reconstruction in the ECAL to make EcalRecHits from digis. This sequence is defined in Reconstruction.cff, which must be included in the cfg:

  #Standard sequences
  include "Configuration/StandardSequences/data/FakeConditions.cff"
  include "Configuration/StandardSequences/data/Reconstruction.cff"

  # all clustering in one sequence
  include "RecoEcal/EgammaClusterProducers/data/ecalClusteringSequence.cff"

  # all clustering in one sequence
  include "RecoEcal/EgammaClusterProducers/data/ecalClusteringSequence.cff"
  
  # path to run everything
  path p = { ecalLocalRecoSequence, ecalClusteringSequence } 

To run the reconstruction, type in RecoEgamma/EgammaClusterProducers/test:

cvs update -r 1.14 egammaDigiToClusters.cfg
cmsRun egammaDigiToClusters.cfg

Contents of output root file:

You can browse the content of the outout root file (egamma_clusters.root in this tutorial) with TBrowser from within root. As shown in this figure, you will see branches in the Events tree whose name include data type (recoBasicClusters and recoSuperClusters), name of the producer(e.g. islandBasicClusterProducer), the secondary optional label (e.g. islandBarrelBasicClustercollection) and the process name (e.g. Rec). Collections ending with __Rec were already in the events that were read in. The collections reconstructed in the example end with __eg, the name of the process in the first line of egammaDigiToClusters.cfg

You can find the complete list of collections made by the producers in RecoEgamma/EgammaClusterProducers, including module names and secondary labels needed to get these collections in framework modules, at SWGuideDataFormatRecoEcal

Note that there are several collections for reco::SuperCluster. In general, the user should only ever need to access two, one for barrel and one for endcaps:

   recoSuperClusters_correctedHybridSuperClusters__eg
   recoSuperClusters_correctedEndcapSuperClustersWithPreshower__eg

Note that the following methods can be used to return the raw SuperCluster energy (i.e. sum of BasicCluster energies, without any corrections):

double reco::SuperCluster::rawEnergy() const;

and the following method will return the preshowerEnergy (returns zero for barrel SuperClusters):

double reco::SuperCluster::preshowerEnergy() const;

Note that it is not recommended to use Island SuperClusters in the barrel. Conversely, Hybrid BasicClusters are not considered useful as standalone objects (i.e. other than components of a hybrid SuperCluster). In particular, Island BasicClusters should be used for ECAL isolation in both barrel and endcaps.

Fill Histograms with EgammaSimpleAnalyzer

The EDAnalyzer class EgammaSimpleAnalyzer is set up in the RecoEcal/EgammaClusterProducers package. The purpose of this simple analyzer module is to illustrate how users can

• fetch collections of objects from Event
• perform simple analysis
• book and store historgrams in output

A dedicated configuration file runEgammaSimpleAnalyzer.cfg has been added to read a collection of events containing the reconstructed BasicClusters and SuperClusters (for example the one produced by running cmsRun egammRecHitsToClusters.cfg).

You need to edit the config file to use your favorite input collection

  source = PoolSource {
    untracked vstring fileNames = {'file:egamma_clusters.root'}
  }

You can run it with:

cmsRun runEgammaSimpleAnalyzer.cfg

The structure and logic of EgammaSimpleAnalyzer is the following (Note: the EgammaSimpleAnalyzer code has advanced relative to what is below, but the following code gives a simple illustration of the concepts):

• all histograms and the root file to store them are data members of EgammaSimpleAnalyzer (by pointer);

class EgammaSimpleAnalyzer : public edm::EDAnalyzer {

 private:
   std::string islandBarrelBasicClusterCollection_;
   std::string islandBarrelBasicClusterProducer_;
   // root file to store histograms
   TFile*  rootFile_;
   std::string outputFile_; // output file
   // per event quantities
   TH1F* h1_nIslandEBBC_;
   TH1F* h1_nIslandEEBC_;
}

• ParameterSet mechanism is used in EgammaSimpleAnalyzer::EgammaSimpleAnalyzer() to provide runtime configuration from user to the module, such as name of the root file to store histograms;

EgammaSimpleAnalyzer::EgammaSimpleAnalyzer( const edm::ParameterSet& ps ) {

  nbinHist_ = ps.getParameter<int>("nbinHist");
  correctedHybridSuperClusterCollection_ = 
                  ps.getParameter<std::string>("correctedHybridSuperClusterCollection");
  correctedHybridSuperClusterProducer_   = 
                  ps.getParameter<std::string>("correctedHybridSuperClusterProducer");
  
  outputFile_   = ps.getParameter<std::string>("outputFile");

}

• histograms are booked and root file is opened (dynamic allocation) in EgammaSimpleAnalyzer::beginJob();

void EgammaSimpleAnalyzer::beginJob(edm::EventSetup const&) {

  // go to *OUR* rootfile and book histograms
  rootFile_->cd();
  h1_nIslandEBBC_ = new TH1F("nIslandEBBC","# basic clusters with island in barrel",
                                                11,-0.5,10.5);
  h1_nIslandEEBC_ = new TH1F("nIslandEEBC","# basic clusters with island in endcap",
                                                11,-0.5,10.5);
}

• analysis is done in EgammaSimpleAnalyzer::analyze() by fetching collections of objects from Event, to iterate over them and fill histograms;

void EgammaSimpleAnalyzer::analyze( const edm::Event& evt, const edm::EventSetup& es ) {
  // Get island basic clusters
  Handle<reco::BasicClusterCollection> pIslandBarrelBasicClusters;
  evt.getByLabel(islandBarrelBasicClusterProducer_, 
  islandBarrelBasicClusterCollection_,pIslandBarrelBasicClusters);
 
  const reco::BasicClusterCollection* islandBarrelBasicClusters = 
            pIslandBarrelBasicClusters.product();
  h1_nIslandEBBC_->Fill(islandBarrelBasicClusters->size());
  
  // loop over the Basic clusters and fill the histogram 
  for(reco::BasicClusterCollection::const_iterator aClus = islandBarrelBasicClusters->begin(); 
              aClus != islandBarrelBasicClusters->end(); aClus++) {
    h1_islandEBBCEnergy_->Fill( aClus->energy() );
    h1_islandEBBCXtals_->Fill(  aClus->getHitsByDetId().size() );
  }
}

• histograms are written into root file in EgammaSimpleAnalyzer::endJob();

void EgammaSimpleAnalyzer::endJob() {

  // go to *OUR* root file and store histograms
  rootFile_->cd();

  h1_nIslandEBBC_->Write();
  h1_nIslandEEBC_->Write();
}

• the dynamically allocated root file is deleted in EgammaSimpleAnalyzer::~EgammaSimpleAnalyzer().
  NB: No need to delete the histograms...

EgammaSimpleAnalyzer::~EgammaSimpleAnalyzer() {
  delete rootFile_;
}

Output of EgammaSimpleAnalyzer

As usual, you can inspect the histograms stored with EgammaSimpleAnalyzer using TBrowser:

Below are a few plots obtained with a sample of single electrons with transverse energy = 35GeV:

Number of BasicClusters contained in Hybrid SuperClusters:

Corrected reconstructed SuperCluster energy, for Hybrid SuperClusters (barrel) and for Island Superclusters (endcaps).

Note on SuperCluster Position

The method SuperCluster::position() returns the position of the SuperCluster as an math::XYZPoint. It is important to note that the pseudorapidity obtained through position().eta() is a detector position coordinate, defined relative to the centre of CMS at (0.,0.,0.). This is not the same as the particle direction, measured relative to the event vertex. Therefore if one wants to compare the position().eta() of a SuperCluster with the η of a generated electron or photon, it is necessary to convert the η of the generated particle to a detector η. The following function can be used to perform this transformation, given the η of the generated particle, and the z and transverse component ρ of the vertex position:

float ecalEta(float EtaParticle , float Zvertex, float RhoVertex)
{
 const float R_ECAL           = 136.5;
 const float Z_Endcap         = 328.0;
 const float etaBarrelEndcap  = 1.479;
 
 if (EtaParticle!= 0.)
    {
      float Theta = 0.0  ;
      float ZEcal = (R_ECAL-RhoVertex)*sinh(EtaParticle)+Zvertex;
      
      if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
      if(Theta<0.0) Theta = Theta+Geom::pi() ;

      float ETA = - log(tan(0.5*Theta));
      
      if( fabs(ETA) > etaBarrelEndcap )
   {
     float Zend = Z_Endcap ;
     if(EtaParticle<0.0 )  Zend = -Zend ;
     float Zlen = Zend - Zvertex ;
     float RR = Zlen/sinh(EtaParticle);
     Theta = atan((RR+RhoVertex)/Zend);
     if(Theta<0.0) Theta = Theta+Geom::pi() ;
     ETA = - log(tan(0.5*Theta));
   }
      return ETA;
    }
  else
    {
      edm::LogWarning("")  << "[EcalPositionFromTrack::etaTransformation] 
      Warning: Eta equals to zero, not correcting" ;
#edit the above 2 lines to be a single line
      return EtaParticle;
    }
}

Cluster shape information

The collection cluster shape variables have a number of uses. The energy in nxn windows is available. The fixed windows are centred on the maximum energy (seed) crystal for odd values of n, and centred on the maximum energy 2x2 which inclues the seed for even values of n. These windows can be used for alternative estimates of the shower energy - for example, 5x5 provides the best energy resolution for unconverted photons in the barrel; and 3x3 is used in the studies of pizero calibration. They can also be used to construct useful variables for cuts in further analysis. For example: R9, defined as the ratio of the 3x3 energy divided by the supercluster energy, provides a very effective photon conversion detector.

Shape variables calculated from the covariance matrix are also available, covEtaEta, covEtaPhi, covPhiPhi, and have been found to be useful in discriminating between signal and background in electron identification.

In order to access ClusterShape information corresponding to a given BasicCluster, it is necessary to use an AssociationMap. The map can be first read in from the event as follows. Note that there are separate maps for the barrel and endcaps:

    // Get association maps linking BasicClusters to ClusterShape
    edm::Handle<reco::BasicClusterShapeAssociationCollection> barrelClShpHandle;
    iEvent.getByLabel("hybridSuperClusters","hybridShapeAssoc", barrelClShpHandle);
    edm::Handle<reco::BasicClusterShapeAssociationCollection> endcapClShpHandle;
    iEvent.getByLabel("islandBasicClusters","islandEndcapShapeAssoc", endcapClShpHandle);

A typical use case is to access the shape of the seed BasicCluster of the SuperCluster. The entry in the map is identified as follows:

    reco::BasicClusterShapeAssociationCollection::const_iterator seedShpItr;

    // Find the entry in the map corresponding to the seed BasicCluster of the SuperCluster
    DetId id = mySuperCluster->seed()->getHitsByDetId()[0];
    if (id.subdetId() == EcalBarrel) {
      seedShpItr = barrelClShpHandle->find(mySuperCluster->seed());
    } else {
      seedShpItr = endcapClShpHandle->find(mySuperCluster->seed());
    }

Note that it is necessary to know whether the BasicCluster concerned is in the barrel or endcap in order to search in the correct AssociationMap. Here we haev used the DetId. One could alternatively require: fabs(mySuperCluster->seed()->eta())<1.479. Finally, the ClusterShapeRef corresponding to the seed BasicCluster is obtained from the map:

    // Get the ClusterShapeRef corresponding to the BasicCluster
    const reco::ClusterShapeRef& seedShapeRef = seedShpItr->val;
    double S1overS9 = seedShapeRef->eMax()/seedShapeRef->e3x3();
  }

Information Sources

https://twiki.cern.ch/twiki/bin/view/CMS/EgammaTutorialSegments

Review status

Reviewer/Editor and Date	Comments
KatiLassilaPerini - 12 Jun 2006	added intro from PTDR vol1 and more info on config files

Responsible: DavidFutyan
Last reviewed by: Chris Seez - 18 Dec 2006