---

CBNT_Atlfast.cxx

Go to the documentation of this file.

#include "AtlfastAlgs/CBNT_Atlfast.h"

#include "AtlfastAlgs/StandardNtupleMaker.h"
#include "AtlfastAlgs/JetRecalibrator.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"

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



#include <cmath> 
#include <pair.h>

#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 "StoreGate/DataHandle.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/51") {
    m_atlfastEventLocation    = "/Event/Atlfast";
    m_jetLocation             = "/Event/AtlfastJets";
    m_jetBLocation            = "/Event/AtlfastBJets";
    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";
    
    declareProperty( "AtlfastEventLocation",    m_atlfastEventLocation );
    declareProperty( "JetLocation",             m_jetLocation );
    declareProperty( "JetBLocation",            m_jetBLocation );
    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);
    
  }
  CBNT_Atlfast::~CBNT_Atlfast(){
    delete m_tesIO;
    delete m_jetCal;
  }
  
//<<<<<< METHOD initialize

StatusCode CBNT_Atlfast::initialize(){

  StatusCode sc;
  
  MsgStream log(messageService(), name());
  
  // access already booked ntuple m_NtupleLocID
  sc = accessNtuple();
  if (!sc) {
    return StatusCode::FAILURE;
  }
  
  NTuple::Tuple* nt=p_nt;
  
  //check the properties
  
  //declare more variables
  
  log << MSG::DEBUG << "declare more variables" << endreq;
  
  // Electrons
  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
  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)
  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)
  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
  sc = nt->addItem ("ATLFAST/NJET"      , m_njet, 0, 40 );
  sc = nt->addItem ("ATLFAST/KFJET"     , m_njet, m_codejet);
  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/PTcalo"    , m_njet, m_ptcalo);
  sc = nt->addItem ("ATLFAST/PTbjet"    , m_njet, m_ptbjet);
  sc = nt->addItem ("ATLFAST/PTujet"    , m_njet, m_ptujet);
  // Jets calibrated by AtlfastB
  sc = nt->addItem ("ATLFAST/NJETB"      , m_njetB, 0, 40 );
  sc = nt->addItem ("ATLFAST/KFJETB"     , m_njetB, m_codejetB);
  sc = nt->addItem ("ATLFAST/PXJETB"     , m_njetB, m_pxjetB);
  sc = nt->addItem ("ATLFAST/PYJETB"     , m_njetB, m_pyjetB);
  sc = nt->addItem ("ATLFAST/PZJETB"     , m_njetB, m_pzjetB);
  sc = nt->addItem ("ATLFAST/EEJETB"     , m_njetB, m_eejetB);
  // History
  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
  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);
  
  
  m_jetCal = new JetRecalibrator();
  m_tesIO = new TesIO();
  
  HeaderPrinter hp("Atlfast StandardNtuleMaker:", 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;
  if(!m_tesIO->copy<ReconstructedParticleCollection>(rp,m_electronLocation)){
    log <<MSG::DEBUG << "Could not find Electrons"<<endreq;
  }
  log << MSG::DEBUG << "processing Electrons: " <<rp.size()<< endreq;
  
  std::vector<ReconstructedParticle*>::const_iterator 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;
    }
  }

  //-------------------
  // Retrieve Photons 
  //-------------------

  rp.clear();
  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;
    }
  }

  //-------------------
  // Retrieve Isolated Muons
  //-------------------
    
  rp.clear();
  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;
    }
  }

  //-------------------
  // Retrieve NonIsolated Muons
  //-------------------
  
  rp.clear();
  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;
    }
  }

  //-------------------
  // Retrieve Jets
  //-------------------
  
  std::vector<Jet*> jv;
  if(!m_tesIO->copy<JetCollection >(jv, m_jetLocation)){
    log <<MSG::DEBUG << "Could not find Jets"<<endreq;
  }
  log << MSG::DEBUG << "Processing Jets: " <<jv.size()<< endreq;
  
  std::vector<Jet*>::const_iterator ijv  = jv.begin();
  
  m_njet = 0;
  
  for (; ijv != jv.end(); ++ijv) {
    m_codejet[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_ptcalo[m_njet] = (*ijv)->pT();
    m_ptbjet[m_njet] = m_jetCal->bCalibrate( (*ijv)->pT() );
    m_ptujet[m_njet] = m_jetCal->uCalibrate( (*ijv)->pT() );
    if ( ++m_njet == m_njet->range().distance() ) break;
  }
  log << MSG::DEBUG << " Jets in ntuple" << endreq;
  
  
  //-------------------
  // Retrieve AtlfastB Jets
  //-------------------
    
  std::vector<Jet*> jvb;
  if(!m_tesIO->copy<JetCollection >(jvb, m_jetBLocation)){
    log <<MSG::DEBUG << "Could not find AtlfastB Jets"<<endreq;
  }
  log << MSG::DEBUG << "Processing AtlfastB Jets: " <<jvb.size()<< endreq;
  
  ijv  = jvb.begin();
  
  m_njetB = 0;
  
  for (; ijv != jvb.end(); ++ijv) {
    m_codejetB[m_njetB]= (*ijv)->pdg_id();
    m_pxjetB[m_njetB]  = (*ijv)->px();
    m_pyjetB[m_njetB]  = (*ijv)->py();
    m_pzjetB[m_njetB]  = (*ijv)->pz();
    m_eejetB[m_njetB]  = (*ijv)->e();
    if ( ++m_njetB == m_njetB->range().distance() ) break;
  }
  log << MSG::DEBUG << " AtlfastB Jets in ntuple" << endreq;
  
  
  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;      
    
  const DataHandle<EventHeader> theEventHeader;
  if(!(m_tesIO->getDH(theEventHeader))){
    log << MSG::ERROR<<"Could not find the event header in the TES"<<endreq;
  }
  
  if(theEventHeader.isValid()){
    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();
  } else {
    log << MSG::ERROR << " Invalid EventHeader data handle" << endreq;
  }
  
  log << MSG::DEBUG << " EventHeader quantities are in ntuple" << 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) 
  {
      p_nt=nt;
      log << MSG::INFO << "Ntuple " <<m_NtupleLocID 
          << " reaccessed! " << endreq;
      return StatusCode::SUCCESS;
  } else
  {
      log << MSG::FATAL << "Cannot reaccess " << m_NtupleLocID << endreq;
      return StatusCode::FAILURE;
  }
}
 




}

---