JetCalibrator.cxx

Go to the documentation of this file.

// ================================================
// JetCalibrator class Implementation
// ================================================
//
#include "AtlfastAlgs/JetCalibrator.h"
#include "AtlfastAlgs/GlobalEventData.h"
#include "AtlfastUtils/IsAssociated.h"

// visitors stuff
#include "AtlfastUtils/JetVisitor.h"
#include "AtlfastUtils/ClusEcalEtSumVisitor.h"
#include "AtlfastUtils/ClusHcalEtSumVisitor.h"


#include "AtlfastEvent/Jet.h"
#include "AtlfastEvent/ICell.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/ParticleCodes.h"
#include "AtlfastEvent/MsgStreamDefs.h"
#include "AtlfastEvent/Cluster.h"
#include "AtlfastEvent/MCparticleCollection.h"

#include "AtlfastUtils/TesIO.h"
#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/FunctionObjects.h"
#include "AtlfastEvent/SimpleKinematic.h"

#include <cmath> 
#include <algorithm>
#include <numeric>
#include <assert.h>

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



// CLHEP,HepMC
#include "GeneratorObjects/McEventCollection.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"


namespace Atlfast{

    

//--------------------------------
// Constructors and destructors
//--------------------------------
  JetCalibrator::JetCalibrator 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator )
    , m_tesIO( 0 ) {
    
    // Setting the parameter defaults.
    
    // Default paths for entities in the TES
    m_inputLocation             = "/Event/AtlfastJets";
    m_outputLocation            = "/Event/AtlfastCalibratedJets";
    
    // This is how you declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    declareProperty( "InputLocation",           m_inputLocation ) ;
    declareProperty( "OutputLocation",          m_outputLocation ) ;
  }
  
  JetCalibrator::~JetCalibrator() {
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "Destructor Called" << endreq;
    
    if (m_tesIO) {
      delete m_tesIO;
    }
  } 
 

 
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode JetCalibrator::initialize(){
    MsgStream log( messageService(), name() ) ;
    
    GlobalEventData* ged = GlobalEventData::Instance();
    m_mcLocation       = ged -> mcLocation();

    //    m_tesIO = new TesIO(eventDataService());  
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());  
    
    HeaderPrinter hp("Atlfast JetCalibrator:", log);
    
    hp.add("TES Locations:              ");
    hp.add(" Uncalibrated Jets from     ", m_inputLocation); 
    hp.add(" Calibrated Jets to         ", m_outputLocation);
    hp.print();
    
    
    return StatusCode::SUCCESS ;
  }
  
  
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode JetCalibrator::finalize(){

    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "finalizing" << endreq;    
    return StatusCode::SUCCESS ;
  }
  
  
  
  //---------------------------------
  // execute() 
  //---------------------------------
  
  StatusCode JetCalibrator::execute( ){
    
    MsgStream log( messageService(), name() ) ;    
    
    std::string mess;
    // make a message logging stream
    
      
    JetCollection* calibratedJets=new JetCollection;
    JetVector uncalibratedJets;
    
    if( ! m_tesIO->copy<JetCollection> ( uncalibratedJets, m_inputLocation ) ) {
      log << MSG::INFO << "No Jets in TES " << endreq;
    } else {
      log << MSG::DEBUG << uncalibratedJets.size()<<" Uncalibrated Jets from TES " << endreq;
    }
    
    JetVector::iterator src;
    for(src = uncalibratedJets.begin(); src != uncalibratedJets.end(); ++src) { 
      
      // retrieve cluster/muons associated to jets
      m_jetVisitor = SimpleAssocsDispatcher(*src, JetVisitor());
      
      Jet* jet = new Jet(**src);
      
      // check if any muons associated to the jet
      // NOTE! currently only muons are associated to the jet
      bool hasAssociatedMuons = IsAssociated<ReconstructedParticle>(*src);
      ReconstructedParticleVector  muons;
      if ( hasAssociatedMuons ){
        // get muons associated to jet
        ReconstructedParticle rp;
        muons = m_jetVisitor.typeVector(rp);
        
        // undo muons
        this->undoMuons(jet, muons);
      }
      
      
      // calibrate the jet
      this->calibrate(jet);
      
      
      if ( hasAssociatedMuons ){
        // add muons
        this->addMuons(jet, muons);
      }
      
      
      log << MSG::DEBUG << "Pushing back a calibrated Jet" << endreq;
      calibratedJets->push_back( jet );
      log<<MSG::DEBUG<<"jet "<<*jet<<endreq;
      
    }
    
    //......................................
    // Register any Jets which have been successfully created in the 
    // transient event store. 
    TesIoStat stat;
    log<<MSG::DEBUG<<"Storing "<< calibratedJets->size() <<" Jets"<<endreq;
    stat = m_tesIO->store( calibratedJets, m_outputLocation ) ;
    if(!stat){
      log<<MSG::ERROR<<"Could not store jets"<<endreq;
      return stat;
    }
    return StatusCode::SUCCESS;
  }
  
  
  
  //----------------------------------
  // calibrate()
  //----------------------------------
  
  void JetCalibrator::calibrate ( Jet* jet )
  {
    
    MsgStream log( messageService(), name() ) ;
    
    Cluster* cp(0);
    // get clusters associated to jet
    ClusterVector clusters = m_jetVisitor.typeVector(*cp);
    // currently only one cluster associated to jet
    ClusterVector::const_iterator cItr = clusters.begin();
    

    // set the cone radius to 0.1 (default = 0.2)
    //double coreRadius = 0.1; 

    TwoCptCell tCC;
    // get ecal deposited transverse energy
    ClusEcalEtSumVisitor 
      eEtVisitor = SimpleAssocsDispatcher((*cItr), ClusEcalEtSumVisitor(*cItr));
    // get hcal deposited transverse energy
    ClusHcalEtSumVisitor 
      hEtVisitor = SimpleAssocsDispatcher((*cItr), ClusHcalEtSumVisitor(*cItr));
    
    double eClusEt = eEtVisitor.getSum(tCC);
    double hClusEt = hEtVisitor.getSum(tCC);
    

    //     
    log << MSG::DEBUG << "Cluster : " << (*cItr)->eT() << endreq;
    log << MSG::DEBUG << "E/H cal : " << eClusEt << " / " <<hClusEt<< endreq;
    log << MSG::DEBUG << "Delta Et: " << (*cItr)->eT()-eClusEt-hClusEt << endreq;
    log << MSG::DEBUG << "alpha/beta: " << alpha((*cItr)->eta(),(eClusEt+hClusEt))
                               << " / " << beta((*cItr)->eta(),(eClusEt+hClusEt)) << endreq;
    
    
    log << MSG::DEBUG << "Ecal: " << endreq;
    log << MSG::DEBUG << "   Cluster Core : " << eEtVisitor.getCore(tCC) << endreq;
    log << MSG::DEBUG << "   Cluster Isolation Ring : " << eEtVisitor.getIsolRing(tCC) << endreq;
    log << MSG::DEBUG << "Hcal: " << endreq;
    log << MSG::DEBUG << "   Cluster Core : " << hEtVisitor.getCore(tCC) << endreq;
    log << MSG::DEBUG << "   Cluster Isolation Ring : " << hEtVisitor.getIsolRing(tCC) << endreq;

    double clusEt = alpha( (*cItr)->eta(), (eClusEt+hClusEt) ) * eClusEt
                   + beta( (*cItr)->eta(), (eClusEt+hClusEt) ) * hClusEt
                   + ( (*cItr)->eT() - (eClusEt+hClusEt) );
    
    log << MSG::DEBUG << " --> Calibrated Cluster Et : " << clusEt << endreq;
    log << MSG::DEBUG << " --> Scale -> CalibratedEt/UncalibrateEt : " 
        << clusEt/(*cItr)->eT() << endreq;
    
    // calibrate jet momentum
    jet->setMomentum( ( clusEt/(jet->eT()) ) * jet->momentum() );
    
    log << MSG::DEBUG << "Calibrated  Jet: " << jet << endreq ;
    
  }
  
  
  
  
  
  double JetCalibrator::alpha(double eta, double et)
  {

    // Make sure et is in units of GeV
    et = et/GeV;

    double a1[21]={-0.0235573, -0.0149172, -0.0112624,  -0.00620699, -0.0184942, 
                   -0.0102938, -0.0180827, -0.00621978, -0.0297793,  -0.0278862, 
                   -0.0117897, -0.0232403, -0.0379402,  -0.0250858,  -0.0303011, 
                   -0.0347847, -0.0374433, -0.121275,   -0.0740044,  -0.156905,  -0.0339867};

    double a2[21]={ 0.0525405,  0.0511558,  0.0431571,   0.0529257,   0.0509096, 
                    0.0491337,  0.0759449,  0.0426821,   0.0712299,   0.0776627, 
                    0.054591,   0.0640166,  0.0533651,   0.0723055,   0.0992998, 
                    0.0918459,  0.0838721,  0.138749,    0.136141,    0.261209,   0.0605451};

    double a3[21]={ 1.01634,    1.01021,    1.01848,     1.00097,     1.02683, 
                    1.04378,    1.02553,    1.03694,     1.03106,     1.01113, 
                    1.02839,    1.05985,    1.14132,     1.11199,     1.03202, 
                    1.00291,    1.02473,    1.05478,     1.00299,     1.04378,    0.986938};

    double y;
    if (std::abs(eta)<=1.0){
      int i = int( std::abs(eta)/0.1 );
      double x = (et>10)? et : 10;
      y = a2[i]*std::exp(a1[i]*x)+a3[i];
    } else if (std::abs(eta)<=3.2){
      int i = 10 + int( (std::abs(eta)-1.0)/0.2 );
      double x = (et>10)? et : 10;
      y = a2[i]*std::exp(a1[i]*x)+a3[i];
    } else {
      y = 1;
    }

    //Scale return value according to units used.
    
    return (y>1)? y/GeV : 1.0/GeV;
  }

  double JetCalibrator::beta(double eta, double et)
  {

    // Make sure et is in units of GeV
    et = et/GeV;

    double b1[21]={-0.0163571, -0.0419745, -0.0172721, -0.0635324, -0.0137478, 
                   -0.0272831, -0.0164739, -0.0226414, -0.0407372, -0.0303132, 
                   -0.0274939, -0.0230969, -0.0112391, -0.0295145, -0.0802241, 
                   -0.0119718, -0.0192763, -0.0325638, -0.01964,   -0.0187867, -0.0341681};

    double b2[21]={ 0.116963,   0.184845,   0.108203,   0.258461,   0.0990949, 
                    0.11335,    0.0541237,  0.109091,   0.144866,   0.106829, 
                    0.0942077,  0.0960772,  0.133047,   0.0779446,  0.151293, 
                    0.0855836,  0.07601,    0.121707,   0.11542,    0.111608,   0.133449};

    double b3[21]={ 1.05823,    1.06982,    1.04692,    1.08131,    1.04635, 
                    1.05313,    1.06387,    1.04447,    1.05965,    1.06855, 
                    1.05192,    1.04673,    1.08641,    1.29738,    1.33545, 
                    1.30941,    1.3002,     1.30029,    1.24078,    1.27179,    1.239};

    double y;
    if (std::abs(eta)<=1.0){
      int i = int( std::abs(eta)/0.1 );
      double x = (et>10)? et : 10;
      y = b2[i]*std::exp(b1[i]*x)+b3[i];
    } else if (std::abs(eta)<=3.2){
      int i = 10 + int( (std::abs(eta)-1.0)/0.2 );
      double x = (et>10)? et : 10;
      y = b2[i]*std::exp(b1[i]*x)+b3[i];
    } else {
      y = 1;
    }

    //Scale return value according to units used.
    
    return (y>1)? y/GeV : 1.0/GeV;
  }





  //-------------------------------------------
  // undomuons
  //-------------------------------------------
  
  // undo muons from the jet
  
  void JetCalibrator::undoMuons(Jet* jet, ReconstructedParticleVector& muons )
  {
    
    if (muons.size() == 0)  return;
    
    MsgStream log( messageService(), name() ) ;
    
    ReconstructedParticleVector::const_iterator first = muons.begin() ;
    ReconstructedParticleVector::const_iterator last  = muons.end() ;
    ReconstructedParticleVector::const_iterator itr;
    for (itr=first; itr!=last; ++itr) {
      log << MSG::DEBUG 
          << "undoing muon with momentum: "
          << (*itr)->momentum() 
          << endreq;
      
      log << MSG::DEBUG << "Old jet momentum: " << *jet << endreq;
      jet->setMomentum(jet->momentum() - (*itr)->momentum()); 
      log << MSG::DEBUG << "New jet momentum: " << *jet << endreq ;
    }

  }
  

 

  //-------------------------------------------
  // addmuons
  //-------------------------------------------
  
  // add muons to the jet
  
  void JetCalibrator::addMuons(Jet* jet, ReconstructedParticleVector& muons )
  {
    
    if (muons.size() == 0)  return;
    
    MsgStream log( messageService(), name() ) ;
    
    // association
    IAOO* ia = jet;
    
    ReconstructedParticleVector::const_iterator first = muons.begin() ;
    ReconstructedParticleVector::const_iterator last  = muons.end() ;
    ReconstructedParticleVector::const_iterator itr;
    for (itr=first; itr!=last; ++itr) {
      log << MSG::DEBUG 
          << "Adding muon with momentum: "
          << (*itr)->momentum() 
          << endreq;
      
      log << MSG::DEBUG << "Old jet momentum: " << *jet << endreq;
      jet->setMomentum(jet->momentum() + (*itr)->momentum()); 
      log << MSG::DEBUG << "New jet momentum: " << *jet << endreq ;
      
      // create association
      ia->associate(*itr);
    }

  }
  

  


} // end namespace bracket

---