#ifdef USE_CDFEDM2
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Component: SourceHEPGHadronsWjets.cc (based on SourceHEPGHerwigHadrons.cc)
//
// Purpose:   This is a concrete class that uses the HEPG bank to fill the EnergyData object. This class inherits from InputSource.
//            It selects stable hadrons before hadronization from a HERWIG 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)
//            09/07/2004  Modified bdc
//
// History:  July 4 2003 (J.-F. Arguin): Minor bug fix noticed by Mario
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



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

SourceHEPGPythiaHadronsWjets::SourceHEPGPythiaHadronsWjets():_sourceName("HEPGHerwigHadronsWjets"){}

Id SourceHEPGPythiaHadronsWjets::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,"SourceHEPGPythiaHadronsWjets: ") << "Could not find HEPG_StorableBank"  << endmsg;
    return 0;
  }
  
  //HEPG_StorableBank::particleIter thisParticle(*_dataHandle, 0);
  
  _iter = 0;
  _end = _dataHandle->n_particles();
  _ietaCount = 0;
  _iphiCount = 0;
  _statusCount = 0;

  return _dataHandle->object_id();  // return the object ID
  
}

Id SourceHEPGPythiaHadronsWjets::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,"SourceHEPGPythiaHadronsWjets: ") << "Could not find HEPG StorableBank with requested object ID"  << endmsg;
    return 0;
  }
  _iter = 0;
  _end = _dataHandle->n_particles();
  _ietaCount = 0;
  _iphiCount = 0;
  _statusCount = 0;

  return 1;
  
}

EnergyData* SourceHEPGPythiaHadronsWjets::next(){ 

  EnergyData* eData = 0;

  while(_iter < _end && _dataHandle->istdhep(_iter) != 1)
    _iter++;

  if(_iter < _end){

    int status_mother = -1;
    int mother  = -1;
    int mother_pdg = -1;
    if(_dataHandle->jmo1hep(_iter)) {
      mother = _dataHandle->jmo1hep(_iter) - 1;
      status_mother = _dataHandle->istdhep(mother);
      mother_pdg = _dataHandle->idhep(mother);
    }
    // debugging output
    
    /*    std::cout << "Particle = " << _iter << std::endl;    
    std::cout << "Particle PDG = " << _dataHandle->idhep(_iter) << std::endl;
    std::cout << "Particle status code = " << _dataHandle->istdhep(_iter) << std::endl;
    std::cout << "Mother  = " << mother << std::endl;
    std::cout << "Mother status code = " << status_mother << std::endl;
    std::cout << "Mother PDG = " << mother_pdg << std::endl;
    */ 
    // in W or Z case, the first two stable (statuscode=1)
    //particles with parents with status 3 are the W decay products .
    //Here we mask all neutrinos and electrons and muons. Taus and their
    //products are left alone to form jets
    // These should not enter as jets.
    // All other stable particles are the final state hadrons and their decay 
    //products, which we do want to regard as jets

    if(!((mother_pdg==11 || mother_pdg==-11
          || mother_pdg==12 || mother_pdg==-12
          || mother_pdg==13 || mother_pdg==-13
          || mother_pdg==14 || mother_pdg==-14
          || mother_pdg==16 || mother_pdg==-16) && status_mother==3)) {
      
      HepLorentzVector vec(_dataHandle->Px(_iter),_dataHandle->Py(_iter),_dataHandle->Pz(_iter),_dataHandle->E(_iter));
      
      if(vec.perp() > 0 && fabs(vec.e()) > fabs(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,"SourceHEPGPythiaHadronsWjets::next()") << "No space to store next hadron." << 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