//////////////////////////////////////////////////////////////////////////
// ================================================
// TheAnalysis 
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace MyAnalysis
//

// Header include
#include "MyAnalysis/TheAnalysis.h"
#include "GaudiKernel/MsgStream.h"
#include "MyAnalysis/Selectors.h"
// Other stuff
#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"
#include "GeneratorObjects/McEventCollection.h"
#include "StoreGate/DataHandle.h"
#include "GaudiKernel/ObjectVector.h"
namespace MyAnalysis {

  //Constructor
  TheAnalysis:: TheAnalysis(const std::string& name, ISvcLocator* pSvcLocator):
    Algorithm(name,pSvcLocator){
  }
  //__________________________________________________________________________
 //Destructor 
  TheAnalysis::~TheAnalysis()
  {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "TheAnalysis destructor called" << endreq;
  }
  //__________________________________________________________________________
  StatusCode TheAnalysis::initialize()
  {
    
    MsgStream log( messageService(), name() );
    log << MSG::INFO <<"Initializing" << endreq;
    StatusCode sc = service("StoreGateSvc", m_sgSvc);
    if (sc.isFailure()) {
      log << MSG::ERROR << "Could not find StoreGateSvc" << endreq;
      return sc;
    }
    return StatusCode::SUCCESS; 
  }
  //_________________________________________________________________________
  StatusCode TheAnalysis::finalize()
  {
    MsgStream log( messageService(), name() );
    log << MSG::INFO <<"finalizing" << endreq;
    return StatusCode::SUCCESS; 
  }
  //_________________________________________________________________________
  StatusCode TheAnalysis::execute()
  {
    //.............................................
    MsgStream log( messageService(), name() );
    log << MSG::DEBUG << "TheAnalysis execute()" << endreq;

    //=========================================
    // get MC particles datahande from the SG
    //=========================================
    log << MSG::DEBUG << "TheAnalysis execute() stage 1: retrieve McEventCollection" << endreq;
    const DataHandle<McEventCollection> mcCollptr;
    if ( m_sgSvc->retrieve(mcCollptr).isFailure() ) {
      log << MSG::ERROR << "Could not retrieve McEventCollection"
	  << endreq;
      return StatusCode::FAILURE;
    }
    //=========================================
    //copy over some of the particles
    //=========================================
    log << MSG::DEBUG << "TheAnalysis execute() stage 2: Copy particles" << endreq;
    std::vector<const HepMC::GenParticle*>  myMCparticles ;
    McEventCollection::const_iterator itr;
    for (itr = mcCollptr->begin(); itr!=mcCollptr->end(); ++itr) {
      HepMC::GenEvent* genEvt = (*itr)->pGenEvt()  ;
      //---------------------------------------
      //Either use HepMC::copy_if
      //---------------------------------------
      HepMC::copy_if( genEvt->particles_begin(), 
      		      genEvt->particles_end(),
      		      back_inserter( myMCparticles),
      		      is_photon()
      		      );
      
      //---------------------------------------
      //OR use std::copy
      //---------------------------------------
      /**
      std::copy( genEvt->particles_begin(), 
      		      genEvt->particles_end(),
      		      back_inserter( myMCparticles)
       		      );
		      
      **/
    }
    //==========================================
    //print the particle pdg_id
    //==========================================
    log << MSG::DEBUG << "TheAnalysis execute() stage 3: Print ids" << endreq;
    std::vector<const HepMC::GenParticle*>::const_iterator src;
    for(src = myMCparticles.begin() ; src != myMCparticles.end() ; ++src)
      { 
	log << MSG::DEBUG << "Particle x momentum " << (*src)->momentum().px() << endreq;
      }

    return StatusCode::SUCCESS; 
  }
 
} // end of namespace bracket