CBNT_Atlfast.cxx

Go to the documentation of this file.

#include "AtlfastAlgs/CBNT_Atlfast.h"
#include "AtlfastAlgs/StandardNtupleMaker.h"

#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/HepMC_helper/IsStatusxx.h"
#include "AtlfastUtils/HepMC_helper/All.h"

#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/Jet.h"
#include "AtlfastEvent/EventHeader.h"
#include "AtlfastEvent/MCweightContainerCollection.h"
#include "AtlfastAlgs/GlobalEventData.h"

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



#include <cmath> 

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/ISvcLocator.h"

#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/SmartDataLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
//#include "StoreGate/StoreGateSvc.h"

#include "GaudiKernel/PropertyMgr.h"

#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"


#include <string>

namespace Atlfast{
//<<<<<< CONSTRUCTOR

  CBNT_Atlfast::CBNT_Atlfast(const std::string& name, 
                             ISvcLocator* pSvcLocator) : 
    Algorithm(name, pSvcLocator),
    m_NtupleLocID("/FILE1/CBNT/3333"),
    m_ntPtr(0),
    m_tesIO(0)
  {
    m_atlfastEventLocation    = "/Event/Atlfast";
    m_jetLocation             = "/Event/AtlfastJets";
    m_electronLocation        = "/Event/AtlfastIsolatedElectrons";
    m_isolatedMuonLocation    = "/Event/AtlfastIsolatedMuons";
    m_nonIsolatedMuonLocation = "/Event/AtlfastNonIsolatedMuons";
    m_photonLocation          = "/Event/AtlfastIsolatedPhotons";
    m_mcTruthLocation         = "/Event/McEventCollection";
    m_bPhysicsLocation        = "/Event/AtlfastBPhysics";
    m_eventHeaderLocation     = "/Event/AtlfastEventHeader";

    m_fillIsolatedElectrons = true;
    m_fillIsolatedPhotons = true;
    m_fillIsolatedMuons = true;
    m_fillNonIsolatedMuons = true;
    m_fillJets = true;
    m_fillHistory = false;
    m_fillEventData = true;

    declareProperty( "AtlfastEventLocation",    m_atlfastEventLocation );
    declareProperty( "JetLocation",             m_jetLocation );
    declareProperty( "ElectronLocation",        m_electronLocation );
    declareProperty( "IsolatedMuonLocation",    m_isolatedMuonLocation );
    declareProperty( "NonIsolatedMuonLocation", m_nonIsolatedMuonLocation );
    declareProperty( "PhotonLocation",          m_photonLocation );
    declareProperty( "McTruthLocation",         m_mcTruthLocation );
    declareProperty( "BPhysicsLocation",        m_bPhysicsLocation );
    declareProperty( "EventHeaderLocation",     m_eventHeaderLocation );
    declareProperty( "NtupleLocID",             m_NtupleLocID);
    declareProperty( "FillIsolatedElectrons",   m_fillIsolatedElectrons);
    declareProperty( "FillIsolatedPhotons",     m_fillIsolatedPhotons);
    declareProperty( "FillIsolatedMuons",       m_fillIsolatedMuons);
    declareProperty( "FillNonIsolatedMuons",    m_fillNonIsolatedMuons);
    declareProperty( "FillJets",                m_fillJets);
    declareProperty( "FillHistory",             m_fillHistory);
    declareProperty( "FillEventData",           m_fillEventData);
    
  }
  CBNT_Atlfast::~CBNT_Atlfast(){
    if (m_tesIO) delete m_tesIO;
  }
  
//<<<<<< METHOD initialize

StatusCode CBNT_Atlfast::initialize(){

  StatusCode sc;
  
  MsgStream log(messageService(), name());
  
  // access already booked ntuple m_NtupleLocID
  sc = accessNtuple();
  if (sc.isFailure()) {
    log << MSG::ERROR << "Problem accessing ntuple" << endreq;
    return StatusCode::FAILURE;
  }
  
  NTuple::Tuple* nt=m_ntPtr;
  
  //check the properties
  
  //declare more variables
  
  log << MSG::DEBUG << "declare more variables" << endreq;
  
  // Electrons
  if(m_fillIsolatedElectrons){
    sc = nt->addItem ("ATLFAST/NELE"      , m_nele, 0, 12 );
    sc = nt->addItem ("ATLFAST/KFELE"     , m_nele, m_codeele);
    sc = nt->addItem ("ATLFAST/PXELE"     , m_nele, m_pxele);
    sc = nt->addItem ("ATLFAST/PYELE"     , m_nele, m_pyele);
    sc = nt->addItem ("ATLFAST/PZELE"     , m_nele, m_pzele);
    sc = nt->addItem ("ATLFAST/EEELE"     , m_nele, m_eeele);
  }
  // Photons
  if(m_fillIsolatedPhotons){
    sc = nt->addItem ("ATLFAST/NPHO"      , m_npho, 0, 12 );
    sc = nt->addItem ("ATLFAST/KFPHO"     , m_npho, m_codepho);
    sc = nt->addItem ("ATLFAST/PXPHO"     , m_npho, m_pxpho);
    sc = nt->addItem ("ATLFAST/PYPHO"     , m_npho, m_pypho);
    sc = nt->addItem ("ATLFAST/PZPHO"     , m_npho, m_pzpho);
    sc = nt->addItem ("ATLFAST/EEPHO"     , m_npho, m_eepho);
  }
  // Muons (isolated)
  if(m_fillIsolatedMuons){
    sc = nt->addItem ("ATLFAST/NMUO"      , m_nmuo, 0, 12 );
    sc = nt->addItem ("ATLFAST/KFMUO"     , m_nmuo, m_codemuo);
    sc = nt->addItem ("ATLFAST/PXMUO"     , m_nmuo, m_pxmuo);
    sc = nt->addItem ("ATLFAST/PYMUO"     , m_nmuo, m_pymuo);
    sc = nt->addItem ("ATLFAST/PZMUO"     , m_nmuo, m_pzmuo);
    sc = nt->addItem ("ATLFAST/EEMUO"     , m_nmuo, m_eemuo);
  }
  // MuonsX (non-isolated)
  if(m_fillNonIsolatedMuons){
    sc = nt->addItem ("ATLFAST/NMUX"      , m_nmux, 0, 12 );
    sc = nt->addItem ("ATLFAST/KFMUX"     , m_nmux, m_codemux);
    sc = nt->addItem ("ATLFAST/PXMUX"     , m_nmux, m_pxmux);
    sc = nt->addItem ("ATLFAST/PYMUX"     , m_nmux, m_pymux);
    sc = nt->addItem ("ATLFAST/PZMUX"     , m_nmux, m_pzmux);
    sc = nt->addItem ("ATLFAST/EEMUX"     , m_nmux, m_eemux);
  }
  // Jets
   if(m_fillJets){
     sc = nt->addItem ("ATLFAST/NJET"       , m_njet, 0, 40 );
     sc = nt->addItem ("ATLFAST/KFJET"      , m_njet, m_jetlabel);
     sc = nt->addItem ("ATLFAST/PXJET"      , m_njet, m_pxjet);
     sc = nt->addItem ("ATLFAST/PYJET"      , m_njet, m_pyjet);
     sc = nt->addItem ("ATLFAST/PZJET"      , m_njet, m_pzjet);
     sc = nt->addItem ("ATLFAST/EEJET"      , m_njet, m_eejet);
     sc = nt->addItem ("ATLFAST/BJETDR"     , m_njet, m_bjetdR);
     sc = nt->addItem ("ATLFAST/CJETDR"     , m_njet, m_cjetdR);
     sc = nt->addItem ("ATLFAST/TAUJETDR"   , m_njet, m_taujetdR);
     sc = nt->addItem ("ATLFAST/PTCALO"     , m_njet, m_ptcalo);
     sc = nt->addItem ("ATLFAST/UCORR"      , m_njet, m_ucorr);
     sc = nt->addItem ("ATLFAST/BTAG"       , m_njet, m_btag);
     sc = nt->addItem ("ATLFAST/BCORR"      , m_njet, m_bcorr);
     sc = nt->addItem ("ATLFAST/TAUTAG"     , m_njet, m_tautag);
     sc = nt->addItem ("ATLFAST/TAUCORR"    , m_njet, m_taucorr);
     sc = nt->addItem ("ATLFAST/TAU1PTAG"   , m_njet, m_tau1ptag);
     sc = nt->addItem ("ATLFAST/TAU3PTAG"   , m_njet, m_tau3ptag);
     sc = nt->addItem ("ATLFAST/TAU1P3PCORR", m_njet, m_tau1p3pcorr);
   }
   // History
   if(m_fillHistory){
     sc = nt->addItem ("ATLFAST/NPART"     , m_npart, 0, 40 );
     sc = nt->addItem ("ATLFAST/KPPART"    , m_npart, m_kppart);
     sc = nt->addItem ("ATLFAST/KSPART"    , m_npart, m_kspart);
     sc = nt->addItem ("ATLFAST/KFPART"    , m_npart, m_kfpart);
     sc = nt->addItem ("ATLFAST/KPMOTH"    , m_npart, m_kpmoth);
     sc = nt->addItem ("ATLFAST/KFMOTH"    , m_npart, m_kfmoth);
     sc = nt->addItem ("ATLFAST/PXPART"    , m_npart, m_pxpart);
     sc = nt->addItem ("ATLFAST/PYPART"    , m_npart, m_pypart);
     sc = nt->addItem ("ATLFAST/PZPART"    , m_npart, m_pzpart);
     sc = nt->addItem ("ATLFAST/EEPART"    , m_npart, m_eepart);
   }
   // Missing
   if(m_fillEventData){
     sc = nt->addItem ("ATLFAST/ISUB"      , m_isub);
     sc = nt->addItem ("ATLFAST/JETB"      , m_njetb);
     sc = nt->addItem ("ATLFAST/JETC"      , m_njetc);
     sc = nt->addItem ("ATLFAST/JETTAU"    , m_njettau);
     sc = nt->addItem ("ATLFAST/PXMISS"    , m_pxmiss);
     sc = nt->addItem ("ATLFAST/PYMISS"    , m_pymiss);
     sc = nt->addItem ("ATLFAST/PXNUE"     , m_pxnue);
     sc = nt->addItem ("ATLFAST/PYNUE"     , m_pynue);
     sc = nt->addItem ("ATLFAST/SUMET"     , m_sumET);
     sc = nt->addItem ("ATLFAST/WEIGHT0"   , m_mcWeight0);
     sc = nt->addItem ("ATLFAST/WEIGHT1"   , m_mcWeight1);
     sc = nt->addItem ("ATLFAST/WEIGHT2"   , m_mcWeight2);
   }
   
   if (sc.isFailure()) {
     log << MSG::ERROR << "Problem adding item to ntuple" << endreq;
     return StatusCode::FAILURE;
   }
   
  //get the Global Event Data using singleton pattern
  GlobalEventData* ged = GlobalEventData::Instance();
  // load the location of the MC in StoreGate
  m_mcLocation       = ged -> mcLocation();
  m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
  
  HeaderPrinter hp("Atlfast CBNT NtupleMaker:", log);
  hp.add("MC location                 ", m_mcTruthLocation);
  hp.add("Electron Location           ", m_electronLocation);
  hp.add("Photon   Location           ", m_photonLocation);
  hp.add("Isolated Muon Location      ", m_isolatedMuonLocation);
  hp.add("Non-Isolated Muon Location  ", m_nonIsolatedMuonLocation);
  hp.add("Jet Location                ", m_jetLocation);
  hp.add("EventHeaderLocation         ", m_eventHeaderLocation );
  hp.print();
  return StatusCode::SUCCESS;

}


//<<<<<< METHOD execute

StatusCode CBNT_Atlfast::execute(){

  MsgStream log(messageService(), name());
  log << MSG::DEBUG << "in execute() ..." << endreq;

  log << MSG::DEBUG << " -> retrieving all Atlfast stuff"<< endreq;




  //-------------------
  // Retrieve Electrons 
  //-------------------
  std::vector<ReconstructedParticle*> rp;
  std::vector<ReconstructedParticle*>::const_iterator irp; 
  if(m_fillIsolatedElectrons){
    if(!m_tesIO->copy<ReconstructedParticleCollection>(rp,m_electronLocation)){
      log <<MSG::DEBUG << "Could not find Electrons"<<endreq;
    }
    log << MSG::DEBUG << "processing Electrons: " <<rp.size()<< endreq;
  
    irp = rp.begin();

    m_nele = 0;
    for (; irp != rp.end(); ++irp) {
      m_codeele[m_nele] = (*irp)->pdg_id();
      m_pxele[m_nele]   = (*irp)->momentum().px();
      m_pyele[m_nele]   = (*irp)->momentum().py();
      m_pzele[m_nele]   = (*irp)->momentum().pz();
      m_eeele[m_nele]   = (*irp)->momentum().e();
      if (++m_nele == m_nele->range().distance() ) {
        log << MSG::WARNING << "Electron list truncated in Ntuple" << endreq;
        break;
      }
    }
  }else{
    log << MSG::DEBUG << "Isolated-Electrons switched off!" << endreq;
  }
  //-------------------
  // Retrieve Photons 
  //-------------------
  rp.clear();
  if(m_fillIsolatedPhotons){
    if(!m_tesIO->copy <ReconstructedParticleCollection>(rp, m_photonLocation)){
      log <<MSG::DEBUG << "Could not find Photons"<<endreq;
    }
    log << MSG::DEBUG << "processing Photons: " <<rp.size()<< endreq;
    
    irp  = rp.begin();
    
    m_npho = 0;
    for (; irp != rp.end(); ++irp) {
      m_codepho[m_npho] = (*irp)->pdg_id();
      m_pxpho[m_npho]   = (*irp)->momentum().px();
      m_pypho[m_npho]   = (*irp)->momentum().py();
      m_pzpho[m_npho]   = (*irp)->momentum().pz();
      m_eepho[m_npho]   = (*irp)->momentum().e();
      if (++m_npho == m_npho->range().distance() ) {
        log << MSG::WARNING << "Photon list truncated in Ntuple" << endreq;
        break;
      }
    }
  }else{
    log << MSG::DEBUG << "Isolated-Photons switched off!" << endreq;
  }
  //-------------------
  // Retrieve Isolated Muons
  //-------------------
  rp.clear();
  if(m_fillIsolatedMuons){
    if(!m_tesIO->copy
       <ReconstructedParticleCollection>(rp, m_isolatedMuonLocation)){
      log <<MSG::DEBUG << "Could not find Isolated Muons"<<endreq;
    }
    log << MSG::DEBUG << "Isolated Muons: " <<rp.size()<< endreq;
  
    irp  = rp.begin();
    
    m_nmuo = 0;
    for (; irp != rp.end(); ++irp) {
      m_codemuo[m_nmuo] = (*irp)->pdg_id();
      m_pxmuo[m_nmuo]   = (*irp)->momentum().px();
      m_pymuo[m_nmuo]   = (*irp)->momentum().py();
      m_pzmuo[m_nmuo]   = (*irp)->momentum().pz();
      m_eemuo[m_nmuo]   = (*irp)->momentum().e();
      if (++m_nmuo == m_nmuo->range().distance() ) {
        log << MSG::WARNING 
            << "Isolated Muon list truncated in Ntuple" << endreq;
        break;
      }
    }
  }else{
    log << MSG::DEBUG << "Isolated-Muons switched off!" << endreq;
  }

  //-------------------
  // Retrieve NonIsolated Muons
  //-------------------
  rp.clear();
  if(m_fillNonIsolatedMuons){
    if(!m_tesIO->copy
       <ReconstructedParticleCollection >(rp, m_nonIsolatedMuonLocation)){
      log <<MSG::DEBUG << "Could not find non Isolated muons"<<endreq;
    }
    log << MSG::DEBUG << "NonIsolated Muons: " <<rp.size()<< endreq;
    irp  = rp.begin();
    
    m_nmux = 0;
    
    for (; irp != rp.end(); ++irp) {
      m_codemux[m_nmux] = (*irp)->pdg_id();
      m_pxmux[m_nmux]   = (*irp)->momentum().px();
      m_pymux[m_nmux]   = (*irp)->momentum().py();
      m_pzmux[m_nmux]   = (*irp)->momentum().pz();
      m_eemux[m_nmux]   = (*irp)->momentum().e();
      if (++m_nmux == m_nmux->range().distance() ) {
        log << MSG::WARNING 
            << "Non-isolated Muon list truncated in Ntuple" << endreq;
        break;
      }
    }
  }else{
    log << MSG::DEBUG << "NonIsolated-Muons switched off!" << endreq;
  }
  //-------------------
  // Retrieve Jets
  //-------------------
  
  std::vector<Jet*> jv;
  std::vector<Jet*>::const_iterator ijv; 
  if(m_fillJets){
    if(!m_tesIO->copy<JetCollection >(jv, m_jetLocation)){
      log <<MSG::DEBUG << "Could not find Jets"<<endreq;
    }
    log << MSG::DEBUG << "Processing Jets: " <<jv.size()<< endreq;
    
    ijv = jv.begin();
    
    m_njet = 0;
    
    for (; ijv != jv.end(); ++ijv) {
      m_jetlabel[m_njet]    = (*ijv)->pdg_id();
      m_pxjet[m_njet]       = (*ijv)->px();
      m_pyjet[m_njet]       = (*ijv)->py();
      m_pzjet[m_njet]       = (*ijv)->pz();
      m_eejet[m_njet]       = (*ijv)->e();
      m_bjetdR[m_njet]      = (*ijv)->getdRbquark();
      m_cjetdR[m_njet]      = (*ijv)->getdRcquark();
      m_taujetdR[m_njet]    = (*ijv)->getdRhadtau();
      m_ptcalo[m_njet]      = (*ijv)->pT();
      m_ucorr[m_njet]       = (*ijv)->lTagCorrFactor();
      m_btag[m_njet]        = (*ijv)->isBTagged() ? 1. : 0.;
      m_bcorr[m_njet]       = (*ijv)->bTagCorrFactor();
      m_tautag[m_njet]      = (*ijv)->isTauTagged("Standard") ? 1. : 0.;
      m_taucorr[m_njet]     = (*ijv)->tauTagCorrFactor("Tau");
      m_tau1ptag[m_njet]    = (*ijv)->isTauTagged("Tau1P3P:1prong") ? 1. : 0.;
      m_tau3ptag[m_njet]    = (*ijv)->isTauTagged("Tau1P3P:3prong") ? 1. : 0.;
      m_tau1p3pcorr[m_njet] = (*ijv)->tauTagCorrFactor("Tau1P3P");
      if ( ++m_njet == m_njet->range().distance() ) {break;}
    }
    log << MSG::DEBUG << " Jets in ntuple" << endreq;
  }else{
    log << MSG::DEBUG << "Atlfast Jets switched off!" << endreq;
  }
  
  if(m_fillHistory){
    MCparticleCollection mc;
  
    HepMC_helper::All* all= new  HepMC_helper::All;
    TesIoStat stat = m_tesIO->getMC(mc, all);
    delete all;
  
    if(stat){
      log << MSG::DEBUG << " found MC truth" << endreq;
    }else{
      log << MSG::ERROR << " MC truth not found in TES" << endreq;
    }

    MCparticleCollection::const_iterator imc = mc.begin();
    m_npart=0;
  
    HepMC_helper::IsStatusxx isf(3);
  
    for(;imc != mc.end(); ++imc){
      if(isf(*imc)){
        m_kppart[m_npart]  = (*imc)->barcode();
        m_kspart[m_npart]  = (*imc)->status()%100;
        m_kfpart[m_npart]  = (*imc)->pdg_id();
        std::pair<int, int> parentCodes = getParentCodes(imc);
        m_kfmoth[m_npart]  = parentCodes.first;      
        m_kpmoth[m_npart]  = parentCodes.second;
        m_pxpart[m_npart]  = (*imc)->momentum().x();
        m_pypart[m_npart]  = (*imc)->momentum().y();
        m_pzpart[m_npart]  = (*imc)->momentum().z();
        m_eepart[m_npart]  = (*imc)->momentum().e();
        ++m_npart;
      }
    
      if ( m_npart == m_npart->range().distance() ) break;
    }
  
    log << MSG::DEBUG << " Particles in ntuple" << endreq;      
  }else{
    log << MSG::DEBUG << "Atlfast History switched off!" << endreq;
  }

  if(m_fillEventData){
    const EventHeader* theEventHeader(0);
    if(!(m_tesIO->getDH(theEventHeader))){
      log << MSG::ERROR<<"Could not find the event header in the TES"<<endreq;
    }
  
    if(theEventHeader){
      m_isub    = 0;
      m_njetb   = theEventHeader->nBJets();
      m_njetc   = theEventHeader->nCJets();
      m_njettau = theEventHeader->nTauJets();
      m_pxmiss  = theEventHeader->pMiss().x();
      m_pymiss  = theEventHeader->pMiss().y();
      m_pxnue   = theEventHeader->pEscaped().x();
      m_pynue   = theEventHeader->pEscaped().y();
      m_sumET   = theEventHeader->sumET();
      
      MCweightContainer weights = theEventHeader->firstFewWeights(3);
      m_mcWeight0 = weights[0];
      m_mcWeight1 = weights[1];
      m_mcWeight2 = weights[2];
    } else {
      log << MSG::ERROR << " Invalid EventHeader data handle" << endreq;
    }
    
    log << MSG::DEBUG << " EventHeader quantities are in ntuple" << endreq;
  }else{
    log << MSG::DEBUG << "Atlfast History switched off!" << endreq;
  }
  
  return StatusCode::SUCCESS;
  
}

  std::pair<int, int> 
  CBNT_Atlfast::getParentCodes(MCparticleCollection::const_iterator& 
                                      imc){
    
    HepMC::GenVertex* gv = (*imc)->production_vertex();
    if ( gv == 0 ){ return std::pair<int, int>(0,0);}
    
    if ( gv->particles_in_size() == 0 ) { return std::pair<int, int>(0,0); }

    HepMC::GenVertex::particle_iterator 
      pi = gv->particles_begin(HepMC::parents);
    if( *pi == 0 ){ return std::pair<int, int>(0,0); }
      

    return std::pair<int, int>( (*pi)->pdg_id(), (*pi)->barcode() );
  }

//<<<<<< METHOD finalize

StatusCode CBNT_Atlfast::finalize() {
    //  MsgStream log(messageService(), name());
    //log << MSG::INFO << "in finalize() ..." << endreq;
  return StatusCode::SUCCESS;
}

  
// reaccess Ntuple (should already have been booked)
StatusCode CBNT_Atlfast::accessNtuple() {
  MsgStream log(messageService(), name());
  
  //try to access it  
  NTuplePtr nt(ntupleService(),"/NTUPLES"+m_NtupleLocID);
  
  if (nt) 
    {
      m_ntPtr=nt;
      log << MSG::INFO << "Ntuple " <<m_NtupleLocID 
          << " reaccessed! " << endreq;
      return StatusCode::SUCCESS;
    } else
      {
        log << MSG::FATAL << "Cannot reaccess " << m_NtupleLocID << endreq;
        return StatusCode::FAILURE;
      }
}

}

---