#ifdef USE_CDFEDM2
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Component: SourceHEPGPythiaPartons.cc
//
// Purpose:   This is a concrete class that uses the HEPG bank to fill the EnergyData object. This class inherits from InputSource.
//            It selects partons before hadronization from a PYTHIA 2->2 QCD event.
//
//            The following variable will be filled in EnergyData:
//            - emPhi: phi of particle
//            - hadPhi: samething as above
//            - emEnergy: 60% of particle energy
//            - hadEnergy: 40% of particle energy
//            - particle's 4-vector
//
//            Note: The ratio em/had is based on an average for jets.
//
//            iEta and iPhi are dummy, and are incremented by one at each new particle. They need to be filled to store in a
//            PhysicsTowerData.
//
// Created:   25/11/2002  Matthias Toennesmann
//
// History:
//            July 4 (J.F. Arguin): Minor bug fix noticed by Mario
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include "Edm/GenericConstHandle.hh"
#include "Calor/SourceHEPGPythiaPartons.hh"
#include "CalorGeometry/CalConstants.hh"
#include "ErrorLogger_i/gERRLOG.hh"
// namespace calor {
  

SourceHEPGPythiaPartons::SourceHEPGPythiaPartons():_sourceName("HEPGPythiaPartons"){}

Id SourceHEPGPythiaPartons::initEvent(AbsEvent* anEvent){
  
  StorableObject::SelectByClassName cName( "HEPG_StorableBank" );
  EventRecord::ConstIterator hepgIterator(anEvent, cName);
  if ( hepgIterator.is_valid() ) 
    _dataHandle = ConstHandle<HEPG_StorableBank>(hepgIterator);
  else {
    _dataHandle.set_null();
    ERRLOG(ELerror,"SourceHEPGPythiaPartons: ") << "Could not find HEPG_StorableBank"  << endmsg;
    return 0;
  }
  
  //HEPG_StorableBank::particleIter thisParticle(*_dataHandle, 0);
  
  _iter = 0;
  _end = _dataHandle->n_particles();
  for(int iParticle = 0; iParticle < _end; iParticle++){
    if(_iter == 0 && _dataHandle->istdhep(iParticle) != 3)
      _iter = iParticle;
    if(_iter != 0 && _end == _dataHandle->n_particles() && _dataHandle->jmo1hep(iParticle) > _iter)
      _end = iParticle;
  }

  _ietaCount = 0;
  _iphiCount = 0;
  
  return _dataHandle->object_id();  // return the object ID
  
}

Id SourceHEPGPythiaPartons::initEvent(AbsEvent* anEvent, Id objectId){

  StorableObject::SelectByObjectId selector( objectId);
  
  EventRecord::ConstIterator hepgIterator( anEvent, selector );
  if ( hepgIterator.is_valid() ) 
    _dataHandle = ConstHandle<HEPG_StorableBank>(hepgIterator);
  else {
    _dataHandle.set_null();
    ERRLOG(ELerror,"SourceHEPGPythiaPartons: ") << "Could not find HEPG StorableBank with requested object ID"  << endmsg;
    return 0;
  }
  _iter = 0;
  _end = _dataHandle->n_particles();
  for(int iParticle = 0; iParticle < _end; iParticle++){
    if(_iter == 0 && _dataHandle->istdhep(iParticle) != 3)
      _iter = iParticle;
    if(_iter != 0 && _end == _dataHandle->n_particles() && _dataHandle->jmo1hep(iParticle) > _iter)
      _end = iParticle;
  }

  _ietaCount = 0;
  _iphiCount = 0;
  
  return 1;
  
}

EnergyData* SourceHEPGPythiaPartons::next(){ 

  EnergyData* eData = 0;

  while(_iter < _end &&
	(_dataHandle->jda1hep(_iter) != 0 && _dataHandle->jda1hep(_iter) < _end ||
	 abs(_dataHandle->idhep(_iter)) >= 11 && abs(_dataHandle->idhep(_iter)) <= 16))
    _iter++;

  if(_iter < _end){

    HepLorentzVector vec(_dataHandle->Px(_iter),_dataHandle->Py(_iter),_dataHandle->Pz(_iter),_dataHandle->E(_iter));

    if(vec.perp() > 0 && abs(vec.e()) > abs(vec.pz())){

      float eta = -log(tan(asin(vec.perp()/vec.rho())/2));

      if(fabs(eta) < 7.5){

	float thePhi = vec.phi();
	if(thePhi < 0) thePhi += TWOPI;

	if(_ietaCount >= TOWER_NETA)
	  ERRLOG(ELerror,"SourceHEPGPythiaPartons::next()") << "No space to store next parton." << endmsg;
	else
	  eData = new EnergyData(thePhi,
				 thePhi,
				 0.60 * _dataHandle->E(_iter),
				 0.40 * _dataHandle->E(_iter),
				 _ietaCount,
				 _iphiCount,
				 vec);

	// Here, we increment iphi by one, and ieta if iphi is at the end of the calorimeter grid
	_iphiCount++;
	if(_iphiCount == TOWER_PHI_SEG[TOWER_TYPE[_ietaCount]]){
	  _ietaCount++;
	  _iphiCount = 0;
	}
      }
    }

    _iter++;
  }

  return eData;
}

#endif  // USE_CDFEDM2