Particle Flow Cluster Calibration

Complete: 1

This page describes the ECAL and HCAL cluster calibration used in the PFlow framework in the charged hadron hypothesis, for neutral hadron identification. Links to the code are provided.

Responsible Persons

  • Patrick Janot, Chris Silkworth

Introduction

Why is a specific calibration needed ?

Even after the basic Hcal and Ecal calibration the reconstructed energy from a PF charged hadron that leaves deposits in either the Hcal, Ecal or both is still far from the actual energy of the particle. This is especially true at low energy, where much of the time the reconstructed energy of the particle is lost to noise and threshold conditions. With this calibration we give a better estimate of the reconstructed energy, improve resolution and provide a smooth transition from the low energy regime to the high energy regime. This is crucial to the particle flow algorithm since excesses of energies in the Ecal or Hcal lead to the reconstruction of other particles, such as neutral hadrons or photons.

How is the calibration applied in the PFlow framework ?

After performing the procedure described below we are left with a set of calibration parameters that describe a function. This function takes in a "true" energy of the PFBlock and returns a calibrated energy. The function, with the proper parameters, is given to RecoParticleFlow/PFClusterTools/src/PFEnergyCalibration.cc which takes in the uncorrected Hcal and Ecal engeries, then calibrates them. The function and the parameter are hard coded into PFEnergyCalibration.cc.

Deriving PF cluster calibration with the full simulation

Preliminary Stuff

  • Setup area, check out necessary packages and compile
cmsrel CMSSW_3_11_2
cd CMSSW_3_11_2/src
cmsenv
cvs co -r V12-03-02-02 RecoParticleFlow/Configuration
cvs co -r V13-02-00-01 RecoParticleFlow/PFRootEvent
cvs co IOMC/ParticleGuns
cvs co Configuration/Generator #not really necessary just get this to keep my directory neat
scram b -j 4
  • Open a new file in emacs called template.py in the directory $CMSSW_BASE/Configuration/Generator/test/ and paste this into it:

import FWCore.ParameterSet.Config as cms

process = cms.Process("PROD")

process.load("Configuration.StandardSequences.SimulationRandomNumberGeneratorSeeds_cff")

process.load("Configuration.StandardSequences.Services_cff")
process.load("Configuration.StandardSequences.GeometryDB_cff")
process.load("Configuration.StandardSequences.MagneticField_38T_cff")
process.load("FWCore.MessageService.MessageLogger_cfi")
process.load('IOMC.EventVertexGenerators.VtxSmearedRealistic7TeVCollision_cfi')
process.load("Configuration.StandardSequences.Generator_cff")
process.load("Configuration.StandardSequences.SimIdeal_cff")
process.load("SimGeneral.MixingModule.mixNoPU_cfi")
process.load('Configuration.StandardSequences.Digi_cff')
process.load("Configuration.StandardSequences.SimL1Emulator_cff")
process.load("Configuration.StandardSequences.DigiToRaw_cff")
process.load("Configuration.StandardSequences.RawToDigi_cff")
process.load('Configuration.StandardSequences.L1Reco_cff')
process.load("Configuration.StandardSequences.Reconstruction_cff")
# Event output                                                                                                                
process.load("Configuration.EventContent.EventContent_cff")

process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
process.GlobalTag.globaltag = 'START311_V1A::All'

#process.GlobalTag.toGet = cms.VPSet(                      #Add these lines if you want to change the tag.
#   cms.PSet(record = cms.string("HcalRespCorrsRcd"),      #In this case the Hcal response tag was changed
#          tag = cms.string("HcalRespCorrs_v3.00_TEST"),   #to HcalRespCorrs_v3.00_TEST.
#          connect = cms.untracked.string("frontier://FrontierProd/CMS_COND_31X_HCAL")
#          )
#   )


#from Configuration.PyReleaseValidation.autoCond import autoCond                                                              
#process.GlobalTag.globaltag = autoCond['mc']                                                                                 

process.maxEvents = cms.untracked.PSet(
    input = cms.untracked.int32(5000)
)
### RANDOM setting (change last digit(s) to make runs different !)                                                            
process.RandomNumberGeneratorService.generator.initialSeed = 12345XXXX

#--- Special: equi-populater log-ranges of momentum 1-10/10-100/100-1000

process.source = cms.Source("EmptySource")
process.generator = cms.EDProducer("ExpoRandomPGunProducer",
    PGunParameters = cms.PSet(
        # you can request more than 1 particle                                                                                
        PartID = cms.vint32(-211),
        MinEta = cms.double(-3.1),
        MaxEta = cms.double(3.1),
        MinPhi = cms.double(-3.14159265359),
   MaxPhi = cms.double(3.14159265359),
        MinP   = cms.double(1.),
   MaxP   = cms.double(2000.),
    ),
    firstRun = cms.untracked.uint32(1),
    AddAntiParticle = cms.bool(False)
)

process.output = cms.OutputModule("PoolOutputModule",
    outputCommands = process.FEVTSIMEventContent.outputCommands,
    fileName = cms.untracked.string('output.root'),
    dataset = cms.untracked.PSet(
    dataTier = cms.untracked.string('GEN-SIM-DIGI-RECO'),
    filterName = cms.untracked.string('')
    )
)


#--- Neither L1RECO, not HLT included ---------------------------------------                                                 
process.p0 = cms.Path(process.generator*process.pgen)
process.p1 = cms.Path(process.psim)
process.p2 = cms.Path(process.pdigi)
process.p3 = cms.Path(process.SimL1Emulator)
process.p4 = cms.Path(process.DigiToRaw)
process.p5 = cms.Path(process.RawToDigi)
###process.L1Reco_step = cms.Path(process.L1Reco)                                                                             
process.p6 = cms.Path(process.reconstruction)
process.outpath = cms.EndPath(process.output)

#-- Regular schedule definition (replaced)                                                                                    
# process.schedule = cms.Schedule(process.generation_step,process.genfiltersummary_step,process.simulation_step,process.digitisation_step,process.L1simulation_step,process.digi2raw_step)                                                                 
# process.schedule.extend(process.HLTSchedule)                                                                                
# process.schedule.extend([process.raw2digi_step,process.L1Reco_step,process.reconstruction_step,process.endjob_step,process.FEVTSIMoutput_step])                                                                                                          


#-- Re-defined short schedule                                                                                                 
process.schedule = cms.Schedule(process.p0,process.p1,process.p2,process.p3,process.p4,process.p5,process.p6,process.outpath)
  • Open a new file in emacs called ExpoRandomPGunProducer.cc in the directory $CMSSW_BASE/IOMC/ParticleGuns/src/ and paste this into it.
/*                                                                                                                            
 *  $Date: 2011/09/02 20:10:05 $                                                                                              
 *  $Revision: 1.8 $                                                                                                          
 *  \author Jean-Roch Vlimant                                                                                                 
 *  modified by S.Abdullin 04/02/2011                                                                                         
 */

#include <ostream>

#include "IOMC/ParticleGuns/interface/ExpoRandomPGunProducer.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

// #include "FWCore/Utilities/interface/Exception.h"                                                                          

// #include "CLHEP/Random/RandExpo.h"                                                                                         

using namespace edm;
using namespace std;

ExpoRandomPGunProducer::ExpoRandomPGunProducer(const ParameterSet& pset) :
   BaseFlatGunProducer(pset)
{


   ParameterSet defpset ;
   ParameterSet pgun_params =
      pset.getParameter<ParameterSet>("PGunParameters") ;

   fMinP = pgun_params.getParameter<double>("MinP");
   fMaxP = pgun_params.getParameter<double>("MaxP");

   produces<HepMCProduct>();
   produces<GenEventInfoProduct>();

   //the explonential generator                                                                                               
   //   fRandomExpoGenerator = new CLHEP::RandExponential(fRandomEngine,fMeanP);                                              

}

ExpoRandomPGunProducer::~ExpoRandomPGunProducer()
{
   // no need to cleanup GenEvent memory - done in HepMCProduct                                                               
}

void ExpoRandomPGunProducer::produce(Event &e, const EventSetup& es)
{

   if ( fVerbosity > 0 )
   {
     std::cout << " ExpoRandomPGunProducer : Begin New Event Generation"
               << std::endl ;
   }
   // event loop (well, another step in it...)                                                                                

   // no need to clean up GenEvent memory - done in HepMCProduct                                                              
   //                                                                                                                         

   // here re-create fEvt (memory)                                                                                            
   //                                                                                                                         
   fEvt = new HepMC::GenEvent() ;

   // now actualy, cook up the event from PDGTable and gun parameters                                                         
   //                                                                                                                         
   // 1st, primary vertex                                                                                                     
   //                                                                                                                         
   //HepMC::GenVertex* Vtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.,0.,0.));                                         
   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0.,0.,0.));

   // loop over particles                                                                                                     
   //                                                                                                                         
   int barcode = 1 ;
   for (unsigned int ip=0; ip<fPartIDs.size(); ++ip)
   {

     //                                                                                                                       

     double pmom   = fRandomGenerator->fire(fMinP, fMaxP);
     double y      = (1./fMinP) * fRandomGenerator->fire(0.0,1.0);
     double f      = 1./pmom;
     bool   accpt  = ( y < f);

     /*                                                                                                                       
         std::cout << std::endl                    
                 << "pmom " << pmom << "  f " << f << "   y " << y                                                          
                 << "   accept" << accpt << std::endl;                                                                      
     */

     //shoot until in the designated range                                                                                    
     while ((pmom < fMinP || pmom > fMaxP) || !accpt)
       {
         pmom   = fRandomGenerator->fire(fMinP, fMaxP);
         y      = (1./fMinP) * fRandomGenerator->fire(0.0,1.0);
         f      = 1./pmom;
         accpt  = (y < f);

         /*                                                                                                                   
         std::cout << "pmom " << pmom << "  f " << f << "   y " << y                                                          
                   << "   accept" << accpt << std::endl;                                                                      
         */

       }

     //     std::cout << "result = " <<  pmom << std::endl;                                                                   


       double eta    = fRandomGenerator->fire(fMinEta, fMaxEta) ;
       double phi    = fRandomGenerator->fire(fMinPhi, fMaxPhi) ;
       int PartID = fPartIDs[ip] ;
       const HepPDT::ParticleData*
          PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
       double mass   = PData->mass().value() ;
       double theta  = 2.*atan(exp(-eta)) ;

       // p -> pt                                                                                                             
       double mom    = pmom;
       double pt     = mom * sin(theta);
       double px     = pt  * cos(phi) ;
       double py     = pt  * sin(phi) ;
       //                                                                                                                     
       double pz     = mom*cos(theta) ;
       double energy2= mom*mom + mass*mass ;
       double energy = sqrt(energy2) ;
       //CLHEP::Hep3Vector p(px,py,pz) ;                                                                                      
       //HepMC::GenParticle* Part =                                                                                           
       //    new HepMC::GenParticle(CLHEP::HepLorentzVector(p,energy),PartID,1);                                              
       HepMC::FourVector p(px,py,pz,energy) ;
       HepMC::GenParticle* Part =
           new HepMC::GenParticle(p,PartID,1);
       Part->suggest_barcode( barcode ) ;
       barcode++ ;
       Vtx->add_particle_out(Part);

       if ( fAddAntiParticle )
       {
          //CLHEP::Hep3Vector ap(-px,-py,-pz) ;                                                                               
          HepMC::FourVector ap(-px,-py,-pz,energy) ;
          int APartID = -PartID ;
          if ( PartID == 22 || PartID == 23 )
          {
             APartID = PartID ;
          }
          //HepMC::GenParticle* APart =                                                                                       
          //   new HepMC::GenParticle(CLHEP::HepLorentzVector(ap,energy),APartID,1);                                          
          HepMC::GenParticle* APart =
             new HepMC::GenParticle(ap,APartID,1);
          APart->suggest_barcode( barcode ) ;
          barcode++ ;
          Vtx->add_particle_out(APart) ;
       }

   }

   fEvt->add_vertex(Vtx) ;
   fEvt->set_event_number(e.id().event()) ;
   fEvt->set_signal_process_id(20) ;

   if ( fVerbosity > 0 )
   {
      fEvt->print() ;
   }

   auto_ptr<HepMCProduct> BProduct(new HepMCProduct()) ;
   BProduct->addHepMCData( fEvt );
   e.put(BProduct);

   auto_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
   e.put(genEventInfo);

   if ( fVerbosity > 0 )
   {
      // for testing purpose only                                                                                             
      // fEvt->print() ; // prints empty info after it's made into edm::Event                                                 
     std::cout << " FlatRandomPGunProducer : Event Generation Done " << std::endl;
   }
}

  • Open a new file in emacs named ExpoRandomPGunProducer.h in directory $CMSSW_BASE/IOMC/ParticleGuns/interface/ and paste this into it.
#ifndef ExpoRandomPGunProducer_H
#define ExpoRandomPGunProducer_H

/** \class ExpoRandomPGunProducer                                                                                             
 *                                                                                                                            
 * Generates single particle gun in HepMC format                                                                              
 * Jean-Roch Vlimant                                                                                                          
 * modificed by S.Abdulin 04/02/2011                                                                                          
 ***************************************/

#include "IOMC/ParticleGuns/interface/BaseFlatGunProducer.h"
#include "CLHEP/Random/RandExponential.h"
namespace edm
{

  class ExpoRandomPGunProducer : public BaseFlatGunProducer
  {

  public:
    ExpoRandomPGunProducer(const ParameterSet & pset);
    virtual ~ExpoRandomPGunProducer();

  private:

    virtual void produce(Event & e, const EventSetup& es);

  protected :

    // data members                                                                                                           

    double            fMinP   ;
    double            fMaxP   ;
    double            fMeanP ;
    CLHEP::RandExponential * fRandomExpoGenerator;

  };
}

#endif
  • Now just add this line close to the top of the file $CMSSW_BASE/IOMC/ParticleGuns/src/SealModule.cc:
#include "IOMC/ParticleGuns/interface/ExpoRandomPGunProducer.h"

  • And these Lines near the bottom of the same file:
using edm::ExpoRandomPGunProducer;
DEFINE_FWK_MODULE(ExpoRandomPGunProducer);

  • Compile
cd $CMSSW_BASE
scram b -j 4

Create K0L, PiPlus, PiMinus datasets from scratch

  • cd to directory:
cd $CMSSW_BASE/Configuration/Generator/test/
  • Copy template.py to another file name that you wish to make datasets of (K0L, PiMinus, PiPlus). The default is PiMinus since if you look in the file the pdgId is set to -211.
  • Open up the copied file and make the necessary changes:
    • line 33:
      process.RandomNumberGeneratorService.generator.initialSeed = 12345XXXX 
      - change the XXXX to random 4-digit number.
    • line 40:
      PartID
      - change to the pdgId of the particle you want to generate.
    • line 55:
      fileName
      - change output.root to the file name of your choice.
  • Run the generator:
cmsRun yourFile_cfg.py

Create Root NTuple

  • cd to directory:
cd $CMSSW_BASE/RecoParticleFlow/Configuration/test 
  • You can either open up one of the pre-written cfgs (PiMinus_cfg.py, PiPlus_cfg.py, K0L_cfg.py) in order to change the source of what to run over (i.e. to the generated files you created above) or you can simply run as is:
cmsRun K0L_cfg.py # or substitute K0L by PiPlus or PiMinus
  • This will create output files that look like pfcalib_k0L.root ( pfcalib_piminus.root, pfcalib_piplus.root).
  • Copy the file over to $CMSSW_BASE/RecoParcticleFlow/PFRootEvent/test

Deriving the calibration coefficients

The calibration "theory"

See Patrick_Calib2.pdf attachment at the bottom.

Explanation of calib.C

This macro is built out of two main classes:
  • ABCAlphaBeta - This class holds all the calibration coefficients (A, B, C, Alpha, Beta, see calibration "theory" for their meaning) for a Particular bin in ETrue (True energy). It takes in a vector of ETrue values (all within the ETrue bin range), a vector of energy deposits in the ecal, a vector of energy deposits in the hcal, and a vector of eta values. The ith term in each vector correspond to the ith event. This class has the member functions compute*(), where the * corresponds to one or two of the calibration coefficients that are to be calculated (If you have energy in both the ecal and the hcal you must use the computeBC() rather than computeB() and computeC() since in this former case B and C are correlated). The computeA() is really just a place holder and simply sets the A value, this could be changed to actually calculate the best value for A in the future. There are also functions computeETrueAverage() and computeETrueRMS() which must be used since these values will be used later on for fitting.
  • Calibration - This class holds the total calibration information for all ETrue values. It takes in ABCAlphaBeta objects and then plots the coefficients vs. Etrue (with error bars in both x and y). These coefficients are then fit to some function which is defined in the main part of the code. We are left with each coefficient as a function of ETrue, and we can now get the calibrated energy. There is a function that will allowed you to draw the calibration graph that is called drawCoeffGraph(coefficient) where coefficient is a string that takes in which coefficient it is you want to draw. There are two member functions of getCalibratedEnergy. One takes eta in as an argument and one does not. The one without the eta, is the calibration with no eta dependence, while the one with eta has the eta dependence.

There are two global functions defined in this macro:

  • getValuesFromTree() - this function simply takes a root tree apart and puts the ETrue, ecal, hcal and eta values into vectors for easier handling.
  • drawGausFit() - This function takes in a TH2 of (true energy - corrected energy) and true energy. It then splits up the TH2 into a bunch of TH1s (one for each ETrue bin) and then approximates each TH1 as a Gaussian. From this Gaussian you get two values: the average and the error. Each of the values are then plotted vs ETrue and are called response (for the averages) and resolution (for the errors).

The main part of the code is near the bottom and is called calib(). The part of the code simply takes apart the tree and gets the vectors of ETrue, ecal energies, hcal energies , etas and phis (although these are not used). Then we create all the vectors of ABCAlphaBeta objects, for which there are 6 different cases, since they will get different calibrations:

  • barrelWithEcalHcal - calibration for barrel with energy deposits in ecal and hcal.
  • barrelWithEcal - calibration for barrel with energy deposits only in ecal.
  • barrelWithHcal - calibration for barrel with energy deposits only in hcal.
  • endcapWithEcalHcal - calibration for endcap with energy deposits in ecal and hcal.
  • endcapWithEcal - calibration for endcap with energy deposits only in ecal.
  • endcapWithHcal - calibration for endcap with energy deposits only in hcal.
However the barrelWithEcal and endcapWithHcal are not really needed since their calibration is almost identical to that of those with ecal and hcal. After we fill and compute all the coefficients for each ABCAlphaBeta object we create the 6 calibration objects (one for each vector of ABCAlphaBeta objects). From this we pass them functions and then fit them to these functions. We then fill a bunch of different TH2s that can be used to make response and resolution plots. You can see all the way at the bottom where the user can decide which plots he/she would like to draw.

Run calib.C

  • cd to directory:
cd $CMSSW_BASE/RecoParticleFlow/PFRootEvent/test   
  • Open up Macros/calib.C in an editor and search for calib(). A few lines below this you should see inputFile = TFile::Open("fileName"). Replace the file name with one of the root files you created above.
  • Open up root and run the macro:
root -l 
.x Macros/calib.C++

Deriving PF cluster calibration with the data

Putting the calibration functions in the global tag

Validating the calibration functions

Responsible: PatrickJanot
- 7 March 2011

Topic attachments
I Attachment History Action Size Date Who Comment
PDFpdf Patrick_Calib2.pdf r1 manage 918.3 K 2011-04-14 - 16:19 ChristopherSilkworth "Theory" of the calibration
Edit | Attach | Watch | Print version | History: r8 < r7 < r6 < r5 < r4 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r8 - 2011-09-02 - ChristopherSilkworth



 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    CMSPublic All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback