// ================================================
// MyAnalysis class Implementation
// ================================================
//
//
// Namespace ATLFast
//

// Utilities
#include <cmath> 
#include "AtlfastCode/Phi.h"
#include "AtlfastCode/HeaderPrinter.h"

// This is part of the STL and allows you to use
// generic sort algorithms (amongst other things)
#include <algorithm>
#include "AtlfastCode/FunctionObjects.h"

// Athena/Gaudi includes
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/SmartRefVector.h"


// This algorithm includes
#include "AtlfastCode/MyAnalysis.h"



//-------------------------------
// Athena magic lines.
// These are so that it can make instances of this algorithm
// You never need to know why or how !!!
//-------------------------------

static const AlgFactory<Atlfast::MyAnalysis> s_factory ;
const IAlgFactory& MyAnalysisFactory = s_factory ;

namespace Atlfast {




//--------------------------------
// Constructors and destructors
//--------------------------------

MyAnalysis::MyAnalysis 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator )
{

  // Setting the parameter defaults.
  m_particleLocation   = DEFAULT_particleLocation ;
  m_histStart          = DEFAULT_histStart ;
  m_hist_pt        = NULL ;
  m_hist_eta        = NULL ;
  m_hist_2d        = NULL ;



  // This is how you declare the paramemters to Gaudi so that
  // they can be over-written via the job options file
  declareProperty( "ParticleLocation", m_particleLocation ) ;
  declareProperty( "HistStart",        m_histStart ) ;

}

// Destructor
MyAnalysis::~MyAnalysis() {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "Destructor called" << endreq;
} 


//---------------------------------
// initialise() 
//---------------------------------

StatusCode MyAnalysis::initialize()
{

  // We must here instantiate items which can only be made after
  // any job options have been set

   MsgStream log( messageService(), name() ) ;
   log << MSG::DEBUG << "Initialising" << endreq;
   HeaderPrinter hp("MyAnalysis:",log);
   hp.add( "ParticleLocation", m_particleLocation ) ;
   hp.add( "HistStart       ",        m_histStart ) ;

   m_hist_pt =
     histoSvc()->book("/stat/simple1d/",
		      m_histStart++,
		      "pT",
		      160,
		      0.0,
		      160.0
		      ); 

   m_hist_eta =
     histoSvc()->book("/stat/simple1d/",
		      m_histStart++,
		      "Eta",
		      100,
		      -5.0,
		      5.0
		      ); 

   m_hist_2d =
     histoSvc()->book("/stat/simple2d/",
		      m_histStart++,
		      "Eta vs pT",
		      100,
		      -5.0,
		      5.0,
		      160,
		      0.0,
		      160.0
		      );
 
   return StatusCode::SUCCESS ;
}


//---------------------------------
// finalise() 
//---------------------------------

StatusCode MyAnalysis::finalize()
{
   MsgStream log( messageService(), name() ) ;
   log << MSG::DEBUG << "Finalizing" << endreq;
   return StatusCode::SUCCESS ;
}


//----------------------------------------------
// execute() method called once per event
//----------------------------------------------
//


StatusCode MyAnalysis::execute( )
{

  //................................
  // make a message logging stream 

    MsgStream log( messageService(), name() ) ;

  //.........................................................
  // This is how you extract the things you want to look at. 
  // Look at the userguide for details of ReconstructedParticle, Cell
  // and Cluster classes etc.

    // The ReconstructedParticles are stored in a collection object which
    // is typedefined by ReconstructedParticleCollection. Here we make
    // a pointer to this collection.

    SmartDataPtr<ReconstructedParticleCollection> 
       particles( this->eventDataService(), m_particleLocation ) ;

    // Test to see if particles successfully retrieved
    if( ! particles ) { 
      log << MSG::DEBUG 
	  << "No ReconstructedParticles found in TES at " 
	  << m_particleLocation 
	  << endreq ;
    }

    // Same for Cells and Clusters

    // Iterate over all particles and dump them out
    if( particles )
      this->dumpParticles( log, particles ) ;

  // .................................
  // go home  

    // No tidying to do

    return StatusCode::SUCCESS ;

}


//-------------------------------------------
// private: dumpParticles
//------------------------------------------


void MyAnalysis::dumpParticles
    ( MsgStream& log, 
      ReconstructedParticleCollection* particles 
    )
{

  log << MSG::DEBUG << "Starting particle dump " << endreq ;

  // Make an iterator to a particle in the collection which 
  // has been supplied
  ReconstructedParticleCollection::iterator particle ;

  // traverse collection of particles
      int noPart = 0;
  for( particle = particles->begin(); 
       particle < particles->end(); 
       ++particle 
     )
    {
      // If you just want to print its pT say, then:
      double pT = (*particle)->pT() ;
      log << MSG::DEBUG << " Particle pT was " << pT << endreq;
      ++noPart;
      m_hist_pt->fill(pT,1.0);
      m_hist_eta->fill((*particle)->eta(),1.0);
      m_hist_2d->fill(pT,(*particle)->eta(),1.0);

   }
  log << MSG::DEBUG << "No. Particles " << noPart << endreq;
  return;
}



} // end of namespace bracket