---

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/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"
// 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 <pair.h>

namespace Atlfast {

  StandardNtupleMaker::StandardNtupleMaker(const std::string& name, ISvcLocator* pSvcLocator) 
    :  Algorithm(name, pSvcLocator)      {
    
    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_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( "JetBLocation",            m_jetBLocation );
    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(){
    delete m_tesIO;
    delete m_jetCal;
  }
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  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_codejet);
        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/PTcalo"    , m_njet, m_ptcalo);
        status = nt1->addItem ("PPJETS/PTbjet"    , m_njet, m_ptbjet);
        status = nt1->addItem ("PPJETS/PTujet"    , m_njet, m_ptujet);
        // Jets calibrated by AtlfastB
        status = nt1->addItem ("PPJETS/NJETB"      , m_njetB, 0, 40 );
        status = nt1->addItem ("PPJETS/KFJETB"     , m_njetB, m_codejetB);
        status = nt1->addItem ("PPJETS/PXJETB"     , m_njetB, m_pxjetB);
        status = nt1->addItem ("PPJETS/PYJETB"     , m_njetB, m_pyjetB);
        status = nt1->addItem ("PPJETS/PZJETB"     , m_njetB, m_pzjetB);
        status = nt1->addItem ("PPJETS/EEJETB"     , m_njetB, m_eejetB);
        // History
        // 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);
        
        //status = nt1->addItem ("PHISTORY/NPART"     , m_npart, 0, 40 );
        //status = nt1->addItem ("PHISTORY/KFPART"    , m_npart, m_kfpart);
        //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);
        //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;
      }
    }

    m_jetCal = new JetRecalibrator();
    log<<MSG::DEBUG<<"returning sith status"<< status<<endreq;
    //    m_tesIO = new TesIO(eventDataService());
    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();
    //    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_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;

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

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

---