#include "TLorentzVector.h"

#include <BaBar/Experiment.hh>
#include <Edm/EventRecord.hh>
#include <Edm/ConstHandle.hh>

#include <AbsEnv/AbsEnv.hh>

#include <evt/Event.hh>

#include <SimulationObjects/HEPG_StorableBank.hh>
#include <SimulationObjects/HEP4_StorableBank.hh>

#include <Stntuple/obj/TGenpBlock.hh>
#include <Stntuple/mod/InitStntupleDataBlocks.hh>

//_____________________________________________________________________________
Int_t StntupleInitGenpBlock(TStnDataBlock* block, AbsEvent* event, int mode) 
{
  // fill generator particle data block (GENP - GENerator Particles - 
  // is the name from Run I)

  TLorentzVector v;
  TGenpBlock* genp = (TGenpBlock*) block;
  genp->Clear();
//-----------------------------------------------------------------------------
// initialization from HEPG, loop over HEPG particles and fill the list
// loop over the existing HEPG banks, add non-initialized interaction
//-----------------------------------------------------------------------------
  HEPG_StorableBank::particleIter hepg_particle;

  int id, ncodes, np(0);

  for (AbsEvent::ConstIterator ib(event,"HEPG_StorableBank"); 
                                                         ib.is_valid(); ib++) {
    ConstHandle<HEPG_StorableBank> h(ib);
    const HEPG_StorableBank& hepg = *h;

    genp->AddPrimaryVertex();

    for (HEPG_StorableBank::particleIter it(hepg); it.is_valid(); ++it) {
      int ip = it.particle();
					// make sure particle has the right
					// code and passes the kinematic cuts

      v.SetXYZT(hepg.Px(it),hepg.Py(it),hepg.Pz(it),hepg.E(it));

      int passed = 1;
      if (genp->UseCuts()) {
	passed    = 0;
	float eta = v.Eta();
	float pt  = v.Pt();

	if ( (pt  >= genp->PtMin ()) && 
	     (eta >= genp->EtaMin()) && (eta < genp->EtaMax()) ) {
	  passed = 1;
	}
      }
      if (! passed) goto NEXT_PARTICLE;
      id     = hepg.idhep(it);
//-----------------------------------------------------------------------------
// if there is no codes defined, pass on all the particles, otherwise store 
// only particles with PDG codes defined by the user
//-----------------------------------------------------------------------------
      ncodes = genp->NPdgCodes();

      if (ncodes == 0) {
	passed = 1;
      }
      else {
	passed = 0;
	for (int i=0; i<genp->NPdgCodes(); i++) {
	  if (abs(id) == genp->PdgCode(i)) {
	    passed = 1;
	    break;
	  }
	}
      }
	
      if (! passed) goto NEXT_PARTICLE;

      genp->NewParticle(np,
			hepg.idhep(it),
			hepg.istdhep(it),
			hepg.jmo1hep(it)-1,
			hepg.jmo2hep(it)-1,
			hepg.jda1hep(it)-1,
			hepg.jda2hep(it)-1,
			hepg.Px(it),
			hepg.Py(it),
			hepg.Pz(it),
			hepg.E(it),
			hepg.Vx(it),
			hepg.Vy(it),
			hepg.Vz(it),
			hepg.Vt(it));
      np++;
    NEXT_PARTICLE:;
    }
  }
//-----------------------------------------------------------------------------
//dt: fill event weight and pt_hat stored in HEP4 bank 
//-----------------------------------------------------------------------------
  for (AbsEvent::ConstIterator ib(event,"HEP4_StorableBank"); 
                                                         ib.is_valid(); ib++) {
    ConstHandle<HEP4_StorableBank> h(ib);
    const HEP4_StorableBank& hep4 = *h;
    genp->SetEventWeight(hep4.eventweightlh());
    genp->SetPthat(hep4.scalelh_6());
  }
  return 0;
}