StandardNtupleMaker.cxx

Go to the documentation of this file.

//  ====================================================================
//  StandardNtupleMaker.cxx
//  -----------------r---------------------------------------------------
//
//  This is the implementation of Standard Ntuple Maker which 
//  creates the combined Ntuple based upon entities in the TES
//  This code is based upon AtlfastSTLNtupAlg by E. Richter-Was 
//  
//
//  Authors   : H.T. Phillips, P. Clarke, P. Sherwood, 
//                 R. Steward, E. Richter-Was
//
//  ====================================================================
#include "AtlfastAlgs/StandardNtupleMaker.h"
#include "AtlfastAlgs/GlobalEventData.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"
// Framework include files
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/NTuple.h"
#include "GaudiKernel/SmartDataPtr.h" //need for ntuple core code

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


#include "GaudiKernel/ObjectList.h"

// Example related include files
 //
#include <cmath> 
#include <utility>

namespace Atlfast {

  StandardNtupleMaker::StandardNtupleMaker(const std::string& name, ISvcLocator* pSvcLocator) 
    :  Algorithm(name, pSvcLocator),
       p_nt1(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_trackLocation           = DEFAULT_trackLocation;
    m_mcTruthLocation         = "/Event/McEventCollection";
    m_bPhysicsLocation        = "/Event/AtlfastBPhysics";
    //    m_triggerLocation         = DEFAULT_triggerLocation;
    m_eventHeaderLocation     = "/Event/AtlfastEventHeader";
    
    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( "TrackLocation", m_trackLocation );
    declareProperty( "McTruthLocation",         m_mcTruthLocation );
    declareProperty( "BPhysicsLocation",        m_bPhysicsLocation );
    //    declareProperty( "TriggerLocation",m_triggerLocation );
    declareProperty( "EventHeaderLocation",     m_eventHeaderLocation );
    
    
  }
  //********************************************************
  StandardNtupleMaker::~StandardNtupleMaker(){
    if (m_tesIO) delete m_tesIO;
  }
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  StatusCode StandardNtupleMaker::initialize() {
    MsgStream log(msgSvc(), name());
    log << MSG::DEBUG << "Initialising SNM" << endreq;
    StatusCode status;
    
    m_nEvents = 0;   
    //
    log << MSG::DEBUG << "Ntuple call 51" << endreq;
    NTuple::Directory* col = 0;
    log << MSG::DEBUG << "About to open ntuple file" << endreq;
    NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");
    if ( file1 )    {
      log << MSG::DEBUG << "Opened ntuple file" << endreq;
      if ( !(col=ntupleSvc()->createDirectory("/NTUPLES/FILE1/Atlfast")) )   {
        log << MSG::ERROR << "Cannot create directory /NTUPLES/FILE1/Atlfast" 
            << endreq;  
        return StatusCode::FAILURE;
      }else{
        log << MSG::DEBUG << "Created ntuple directory" << endreq;
      }
    }
    else    {
      log << MSG::ERROR << "Cannot access file /NTUPLES/FILE1" << endreq;  
      return StatusCode::FAILURE;
    }
    // First: A column wise N tuple
    NTuplePtr nt1(ntupleSvc(), "/NTUPLES/FILE1/Atlfast/51");
    if ( !nt1 )    {    // Check if already booked
      nt1 = ntupleSvc()->book (col, 51, CLID_ColumnWiseTuple, "Atlfast51");
      if ( nt1 )    {
        log << MSG::DEBUG << "booked ntuple " << endreq;
        // Add an index column
        // Electrons
        status = nt1->addItem ("PLEPTONS/NELE"      , m_nele, 0, 12 );
        status = nt1->addItem ("PLEPTONS/KFELE"     , m_nele, m_codeele);
        status = nt1->addItem ("PLEPTONS/PXELE"     , m_nele, m_pxele);
        status = nt1->addItem ("PLEPTONS/PYELE"     , m_nele, m_pyele);
        status = nt1->addItem ("PLEPTONS/PZELE"     , m_nele, m_pzele);
        status = nt1->addItem ("PLEPTONS/EEELE"     , m_nele, m_eeele);
        // Photons
        status = nt1->addItem ("PPHOTONS/NPHO"      , m_npho, 0, 12 );
        status = nt1->addItem ("PPHOTONS/KFPHO"     , m_npho, m_codepho);
        status = nt1->addItem ("PPHOTONS/PXPHO"     , m_npho, m_pxpho);
        status = nt1->addItem ("PPHOTONS/PYPHO"     , m_npho, m_pypho);
        status = nt1->addItem ("PPHOTONS/PZPHO"     , m_npho, m_pzpho);
        status = nt1->addItem ("PPHOTONS/EEPHO"     , m_npho, m_eepho);
        // Muons (isolated)
        status = nt1->addItem ("PLEPTONS/NMUO"      , m_nmuo, 0, 12 );
        status = nt1->addItem ("PLEPTONS/KFMUO"     , m_nmuo, m_codemuo);
        status = nt1->addItem ("PLEPTONS/PXMUO"     , m_nmuo, m_pxmuo);
        status = nt1->addItem ("PLEPTONS/PYMUO"     , m_nmuo, m_pymuo);
        status = nt1->addItem ("PLEPTONS/PZMUO"     , m_nmuo, m_pzmuo);
        status = nt1->addItem ("PLEPTONS/EEMUO"     , m_nmuo, m_eemuo);
        // MuonsX (non-isolated)
        status = nt1->addItem ("PMUXS/NMUX"      , m_nmux, 0, 12 );
        status = nt1->addItem ("PMUXS/KFMUX"     , m_nmux, m_codemux);
        status = nt1->addItem ("PMUXS/PXMUX"     , m_nmux, m_pxmux);
        status = nt1->addItem ("PMUXS/PYMUX"     , m_nmux, m_pymux);
        status = nt1->addItem ("PMUXS/PZMUX"     , m_nmux, m_pzmux);
        status = nt1->addItem ("PMUXS/EEMUX"     , m_nmux, m_eemux);
        // Jets
        status = nt1->addItem ("PPJETS/NJET"       , m_njet, 0, 40 );
        status = nt1->addItem ("PPJETS/KFJET"      , m_njet, m_jetlabel);
        status = nt1->addItem ("PPJETS/PXJET"      , m_njet, m_pxjet);
        status = nt1->addItem ("PPJETS/PYJET"      , m_njet, m_pyjet);
        status = nt1->addItem ("PPJETS/PZJET"      , m_njet, m_pzjet);
        status = nt1->addItem ("PPJETS/EEJET"      , m_njet, m_eejet);
        status = nt1->addItem ("PPJETS/BJETDR"     , m_njet, m_bjetdR);
        status = nt1->addItem ("PPJETS/CJETDR"     , m_njet, m_cjetdR);
        status = nt1->addItem ("PPJETS/TAUJETDR"   , m_njet, m_taujetdR);
        status = nt1->addItem ("PPJETS/PTCALO"     , m_njet, m_ptcalo);
        status = nt1->addItem ("PPJETS/UCORR"      , m_njet, m_ucorr);
        status = nt1->addItem ("PPJETS/BTAG"       , m_njet, m_btag);
        status = nt1->addItem ("PPJETS/BCORR"      , m_njet, m_bcorr);
        status = nt1->addItem ("PPJETS/TAUTAG"     , m_njet, m_tautag);
        status = nt1->addItem ("PPJETS/TAUCORR"    , m_njet, m_taucorr);
        status = nt1->addItem ("PPJETS/TAU1PTAG"   , m_njet, m_tau1ptag);
        status = nt1->addItem ("PPJETS/TAU3PTAG"   , m_njet, m_tau3ptag);
        status = nt1->addItem ("PPJETS/TAU1P3PCORR", m_njet, m_tau1p3pcorr);
        // History
        status = nt1->addItem ("PHISTORY/NPART"     , m_npart, 0, 40 );
        status = nt1->addItem ("PHISTORY/KPPART"    , m_npart, m_kppart);
        status = nt1->addItem ("PHISTORY/KSPART"    , m_npart, m_kspart);
        status = nt1->addItem ("PHISTORY/KFPART"    , m_npart, m_kfpart);
        status = nt1->addItem ("PHISTORY/KPMOTH"    , m_npart, m_kpmoth);
        status = nt1->addItem ("PHISTORY/KFMOTH"    , m_npart, m_kfmoth);
        status = nt1->addItem ("PHISTORY/PXPART"    , m_npart, m_pxpart);
        status = nt1->addItem ("PHISTORY/PYPART"    , m_npart, m_pypart);
        status = nt1->addItem ("PHISTORY/PZPART"    , m_npart, m_pzpart);
        status = nt1->addItem ("PHISTORY/EEPART"    , m_npart, m_eepart);
        // Missing
        status = nt1->addItem ("INFO/ISUB"      , m_isub);
        status = nt1->addItem ("INFO/JETB"      , m_njetb);
        status = nt1->addItem ("INFO/JETC"      , m_njetc);
        status = nt1->addItem ("INFO/JETTAU"    , m_njettau);
        status = nt1->addItem ("PMISSING/PXMISS"    , m_pxmiss);
        status = nt1->addItem ("PMISSING/PYMISS"    , m_pymiss);
        status = nt1->addItem ("PMISSING/PXNUE"     , m_pxnue);
        status = nt1->addItem ("PMISSING/PYNUE"     , m_pynue);
        status = nt1->addItem ("MCWEIGHT/W0"        , m_mcWeight0);
        status = nt1->addItem ("MCWEIGHT/W1"        , m_mcWeight1);
        status = nt1->addItem ("MCWEIGHT/W2"        , m_mcWeight2);
        //tracks
        //status = nt1->addItem ("PTRACKS/NTRA"      , m_ntra, 0, 500 );
        //status = nt1->addItem ("PTRACKS/KPTRA"     , m_ntra, m_kpTruth);
        //status = nt1->addItem ("PTRACKS/KFTRA"     , m_ntra, m_kfTruth);
        //
        //status = nt1->addItem ("PTRACKS/KPM1TRA"   , m_ntra, m_kpm1tra);
        //status = nt1->addItem ("PTRACKS/KFM1TRA"   , m_ntra, m_kfm1tra);
        //status = nt1->addItem ("PTRACKS/KPM2TRA"   , m_ntra, m_kpm2tra);
        //status = nt1->addItem ("PTRACKS/KFM2TRA"   , m_ntra, m_kfm2tra);
        //status = nt1->addItem ("PTRACKS/KPM3TRA"   , m_ntra, m_kpm3tra);
        //status = nt1->addItem ("PTRACKS/KFM3TRA"   , m_ntra, m_kfm3tra);
        //status = nt1->addItem ("PTRACKS/KPM4TRA"   , m_ntra, m_kpm4tra);
        //status = nt1->addItem ("PTRACKS/KFM4TRA"   , m_ntra, m_kfm4tra);
        //status = nt1->addItem ("PTRACKS/KPM5TRA"   , m_ntra, m_kpm5tra);
        //status = nt1->addItem ("PTRACKS/KFM5TRA"   , m_ntra, m_kfm5tra);
        //status = nt1->addItem ("PTRACKS/KPM6TRA"   , m_ntra, m_kpm6tra);
        //status = nt1->addItem ("PTRACKS/KFM6TRA"   , m_ntra, m_kfm6tra);
          //
        //status = nt1->addItem ("PTRACKS/D0TRA  "      , m_ntra, m_d0Truth);
        //status = nt1->addItem ("PTRACKS/Z0TRA  "      , m_ntra, m_z0Truth);
        //status = nt1->addItem ("PTRACKS/PHITRA "      , m_ntra, m_phiTruth);
        //status = nt1->addItem ("PTRACKS/COTTRA "      , m_ntra, m_cotTruth);
        //status = nt1->addItem ("PTRACKS/PTINVTRA "    , m_ntra, m_ptInvTruth);
        //
        //status = nt1->addItem ("PTRACKS/D0TRACRU  "   , m_ntra, m_d0Track);
        //status = nt1->addItem ("PTRACKS/Z0TRACRU  "   , m_ntra, m_z0Track);
        //status = nt1->addItem ("PTRACKS/PHITRACRU "   , m_ntra, m_phiTrack);
        //status = nt1->addItem ("PTRACKS/COTTRACRU "   , m_ntra, m_cotTrack);
        //status = nt1->addItem ("PTRACKS/PTINVTRACRU " , m_ntra, m_ptInvTrack);
        //
        //status = nt1->addItem ("PTRACKS/CORR11TRA "   , m_ntra, m_corr11tra);
        //status = nt1->addItem ("PTRACKS/CORR21TRA "   , m_ntra, m_corr12tra);
        //status = nt1->addItem ("PTRACKS/CORR31TRA "   , m_ntra, m_corr13tra);
        //status = nt1->addItem ("PTRACKS/CORR41TRA "   , m_ntra, m_corr14tra);
        //status = nt1->addItem ("PTRACKS/CORR51TRA "   , m_ntra, m_corr15tra);
        //status = nt1->addItem ("PTRACKS/CORR22TRA "   , m_ntra, m_corr22tra);
        //status = nt1->addItem ("PTRACKS/CORR23TRA "   , m_ntra, m_corr23tra);
        //status = nt1->addItem ("PTRACKS/CORR24TRA "   , m_ntra, m_corr24tra);
        //status = nt1->addItem ("PTRACKS/CORR25TRA "   , m_ntra, m_corr25tra);
        //status = nt1->addItem ("PTRACKS/CORR33TRA "   , m_ntra, m_corr33tra);
        //status = nt1->addItem ("PTRACKS/CORR34TRA "   , m_ntra, m_corr34tra);
        //status = nt1->addItem ("PTRACKS/CORR34TRA "   , m_ntra, m_corr34tra);
        //status = nt1->addItem ("PTRACKS/CORR35TRA "   , m_ntra, m_corr35tra);
        //status = nt1->addItem ("PTRACKS/CORR44TRA "   , m_ntra, m_corr44tra);
        //status = nt1->addItem ("PTRACKS/CORR45TRA "   , m_ntra, m_corr45tra);
        //status = nt1->addItem ("PTRACKS/CORR55TRA "   , m_ntra, m_corr55tra);
      }
      else    {   // did not manage to book the N tuple....
        return StatusCode::FAILURE;
      }
    }

    log<<MSG::DEBUG<<"returning sith status"<< status<<endreq;


    GlobalEventData* ged = GlobalEventData::Instance();
    // load the location of the MC in StoreGate
    m_mcLocation       = ged -> mcLocation();

    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());

    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();
    //    declareProperty( "TriggerLocation",m_triggerLocation );


    return StatusCode::SUCCESS;
  }
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  StatusCode StandardNtupleMaker::execute() {
    // This just makes the code below a bit easier to read (and type)
    MsgStream log(msgSvc(), name());
    log << MSG::DEBUG << "****Start execute()****" << endreq;
    
    // Counter of events processed
    if( m_nEvents % 50 == 0) {
      log << MSG::INFO << "    processing event number " 
          << m_nEvents << endreq;
    }     
    
    m_nEvents++;  
    
    StatusCode status;

    //-------------------
    // 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_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;


    //-------------------
    // Retrieve Status 21 quarks (fudged with final state collector FIXME)
    //-------------------
    //***PS    std::vector<HepMC::GenParticle*> mc;
    //***PS    m_tesIO->getMC( mc, m_mcTruthLocation, &isf) ;
    //***PS
    //***PS    log << MSG::DEBUG << "Processing partons" <<mc.size()<< endreq;
    //***PS
    //***PS    std::vector<HepMC::GenParticle*>::const_iterator imc  = mc.begin();
    
    MCparticleCollection mc;
    //    TesIoStat stat = m_tesIO->getMC(mc, &isf);

    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;
        //}else{
        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;

      //      m_kfpart[m_npart]  = (*imc)->pdg_id();
      //      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();
    }
    
    log << MSG::DEBUG << " Particles in ntuple" << endreq;      

    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();


      MCweightContainer weights = theEventHeader->firstFewWeights(3);
      // MCweightContainer weights(3, 0.);
      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;

    // Write out ntuples
    
    log << MSG::DEBUG << "Writing ntuple part 51" << endreq;
    status = ntupleSvc()->writeRecord("/NTUPLES/FILE1/Atlfast/51");
    if ( !status.isSuccess() ) log <<MSG::ERROR 
                                   <<"Cannot fill id 51"
                                   <<endreq;
    log << MSG::DEBUG << "Done all the ntuple writing" << endreq;
    return StatusCode::SUCCESS;
  }
  std::pair<int, int> 
  StandardNtupleMaker::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() );
  }

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  StatusCode StandardNtupleMaker::finalize() {
    
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "    Number of events processed      = " << m_nEvents << endreq;
    
    return StatusCode::SUCCESS;
  }
  
  
} // end of namesapce bracket

---