//_____________________________________________________________________________
// Mar 10 2001 P.Murat: initialization of GENP block starting from the
//                      TMcGenEvent
//                      GENP block is missing event weight
//_____________________________________________________________________________
#include <cmath>
#include "TROOT.h"
#include "TFolder.h"
#include "TParticle.h"
#include "TMcGenEvent.hh"

#include "Stntuple/obj/TStnNode.hh"
#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/obj/TGenpBlock.hh"

//_____________________________________________________________________________
int StnGenInitGenpBlock(TStnDataBlock* Block, 
			TStnEvent*     Event, 
			Int_t          Mode ) 
{
  // init GENP block starting from the output of MC generators

  static TMcGenEvent* mcev = 0;

  int        pdg_code, passed, ncodes, np, nint;
  Double_t   eta, pt;
  TParticle* part;

  if (! mcev) {
    mcev = (TMcGenEvent*) 
      gROOT->GetRootFolder()->FindObject("Ana/GeneratorModule/McGenEvent");
  }

  TGenpBlock* genp = (TGenpBlock*) Block;
  genp->Clear();
//-----------------------------------------------------------------------------
// this is just a convention - MC event should be named "McGenEvent"
//-----------------------------------------------------------------------------
  if (! mcev ) return 0;
  nint = mcev->NPrimaryInteractions();

  for (int i=0; i<nint; i++) {
    TPrimaryInteraction* pint = mcev->PrimaryInteraction(i);
    genp->AddPrimaryVertex();
    np = pint->NParticles();
    for (int ip=0; ip<np; ip++) {

      part = pint->Particle(ip);

      passed = 1;
      if (genp->UseCuts()) {
	passed    = 0;
	eta = part->Eta();
	pt  = part->Pt();
	
	if ( (pt > genp->PtMin()) && 
	     (eta > genp->EtaMin()) && (eta <= genp->EtaMax()) ) {
	  passed = 1;
	}
      }
      if (! passed) goto NEXT_PARTICLE;

      pdg_code = part->GetPdgCode();
      ncodes = genp->NPdgCodes();
      if (ncodes == 0) {
					// if there is no codes defined, pass
					// on all the particles
	passed = 1;
      }
      else {
	passed = 0;
					// store only the particles with 
					// PDG codes defined by the user

	for (int i=0; i<genp->NPdgCodes(); i++) {
	  if (TMath::Abs(pdg_code) == genp->PdgCode(i)) {
	    passed = 1;
	    break;
	  }
	}
      }
      
      if (passed) {
	genp->NewParticle(ip,part);
      }
    NEXT_PARTICLE:;
    }
  }

  return 0;
}