//*LastUpdate :  jsf-1-14  31-January-2000  By Akiya Miyamoto
//*LastUpdate :  jsf-1-7-2  16-April-1999  By Akiya Miyamoto
//*LastUpdate :  jsf-1-7  6-April-1999  By Akiya Miyamoto
//*-- Author  : Akiya Miyamoto  6-April-1999

///////////////////////////////////////////////////////////////////
//
// JSFReadGenerator
//
// Read a generator data file, and stored into JSFGeneratorParticle class
//
// Currently, HEPEVT format file is supported. 
//
// HEPEVT format:  This file is accessed by subroutines in readgenutil.F.
// The file is a Fortran binary file.  and should be created as follows.
//
// INTEGER*4 NENDIAN
// INTEGER*4 NVERS
// INTEGER*4 NEVHEP
// INTEGER*4 NHEP
// INTEGER*4 ISTHEP(NMXHEP), IDHEP(NMXHEP)
// INTEGER*4 JMOHEP(2,NMXHEP), JDAHEP(2,NMXHEP)
// REAL*8    PHEP(5,NMXHEP), VHEP(4,NMXHEP)
//
//   IU=1
//   NENDIAN=1296651082
//   WRITE(IU) NENDIAN, NVERS, NEVHEP, NHEP, (ITHEP(K),K=1,NHEP),
//  .          (IDHEP(K),K=1,NHEP),
//  .          (JMOHEP(1,K),JMOHEP(2,K),K=1,NHEP),
//  .          (JDAHEP(1,K),JDAHEP(2,K),K=1,NHEP),
//  .          ((PHEP(J,K),J=1,5),K=1,NHEP),
//  .          ((VHEP(J,K),J=1,5),K=1,NHEP)
//
// (Update)
//   16-Apr-2000 A.Miyamoto  Put bug fixes written by I.Nakamura in ReadOneRecord()
//
//   31-Jan-2000 A.Miyamoto  Mange jmohep[i][0]=0 and isthep[i]=1 | 2 case in ReadOneRecord()
//   21-Jan-2003 A.Miyamoto  Create ReadHepEvent method.
//
//$Id: JSFReadGenerator.cxx,v 1.7 2003/01/22 00:17:20 miyamoto Exp $
//
//////////////////////////////////////////////////////////////////

#ifdef __USEISOCXX__
#include <iostream>
#else
#include <iostream.h>
#endif

#include "TPythia6.h"
#include "JSFSteer.h"
#include "JSFLCFULL.h"
#include "JSFGenerator.h"
#include "JSFReadGenerator.h"
#include <TMath.h>

ClassImp(JSFReadGenerator)
ClassImp(JSFReadGeneratorBuf)

extern "C" {
extern void readgenopen_(Int_t *unit, Char_t *name, Int_t len);
extern void readgenclose_(Int_t *unit);
extern void readgenhepevt_(Int_t *iunit, 
	Int_t *endian, Int_t *nvers, Int_t *nevhep, Int_t *nhep,
        Int_t isthep[], Int_t idhep[],
        Int_t jmohep[][2], Int_t jdahep[][2],
	Double_t phep[][5], Double_t vhep[][4], Int_t *nret);

#if __PYTHIA_VERSION__ >= 6 
  extern Int_t pychge_(Int_t *kf);
  extern int pycomp_(int *kf);
#else
  extern Int_t luchge_(Int_t *kf);
  extern Float_t ulctau_(Int_t *kf);
#endif

};

#if __PYTHIA_VERSION__ >= 6 
struct COMMON_PYDAT2_t {
  int KCHG[4][500];
  double PMAS[4][500];
  double PARF[2000];
  double VCKM[4][4];
};

extern COMMON_PYDAT2_t pydat2_;
#endif

using namespace std;

//_____________________________________________________________________________
JSFReadGenerator::JSFReadGenerator(const char *name, const char *title)
       : JSFGenerator(name,title)
{
  fEventBuf = new JSFReadGeneratorBuf("JSFReadGeneratorBuf", 
	       "JSFReadGenerator event buffer",this);
 
  sscanf(gJSF->Env()->GetValue("JSFReadGenerator.DataFile","genevent.dat"),
	 "%s",fDataFileName);
  sscanf(gJSF->Env()->GetValue("JSFReadGenerator.Unit","10"),"%d",&fUnit);
  sscanf(gJSF->Env()->GetValue("JSFReadGenerator.Format","HEPEVT"),
	 "%s",fFormat);

}



// ---------------------------------------------------------------
Bool_t JSFReadGenerator::BeginRun(Int_t nrun)
{
  if( !IsWritable() ) return kTRUE; // Quit when not output mode.
  Bool_t rc=kTRUE;
// Open a file 
  Int_t lenf=strlen(fDataFileName);
  readgenopen_(&fUnit, fDataFileName, lenf);

  return rc;
}

// ---------------------------------------------------------------
Bool_t JSFReadGenerator::EndRun()
{
  if( !IsWritable() ) return kTRUE;
  readgenclose_(&fUnit);
  return kTRUE;
}


// ---------------------------------------------------------------
Bool_t JSFReadGenerator::Process(Int_t nev)
{
//

  JSFReadGeneratorBuf *simb=(JSFReadGeneratorBuf*)EventBuf();
  return simb->ReadOneRecord();
}

// ---------------------------------------------------------------
 JSFReadGeneratorBuf::JSFReadGeneratorBuf(const char *name,
     const char *title, JSFReadGenerator *module)
  :JSFGeneratorBuf(name, title, (JSFGenerator*)module)
{
  SetStartSeed(12345);
}



// ---------------------------------------------------------------
 Bool_t JSFReadGeneratorBuf::ReadHepEvent(const Int_t maxhep, Int_t &nevhep, Int_t &nhep, 
	 Int_t isthep[], Int_t idhep[], Int_t jmohep[][2], Int_t jdahep[][2],
	 Double_t phep[][5], Double_t vhep[][4])
{
// Read Generator data and saved to the class.

// Due to imcompatibility in JSFGeneratorParticle class and HEPEVT,
// second mother information is not saved.
//
// 

  Int_t nvers=1;
  Int_t ivers, endian ;
  Int_t nret;
  Int_t unit=((JSFReadGenerator*)Module())->GetUnit();
  readgenhepevt_(&unit,&endian, &ivers, &nevhep,  &nhep,
	isthep, idhep, jmohep, jdahep, phep, vhep, &nret);

  if( nvers != ivers ) {
    printf("Errror at JSFReadGenerator.. Version number of ");
    printf("data is %d, while program can process only ",ivers);
    printf("version %dn",nvers);
    return kFALSE;
  }
  else if(nret == -1) { 
    printf("In JSFReadGeneratorBuf::UnpackDST .. ");
    printf("read end-of-file of the input file.n");
    return kFALSE;
  }
  else if( nret == -2 ) {
    printf("In JSFReadGeneratorBuf::UnpackDST .. ");
    printf("read error detected.n");
    return kFALSE;
  }

  return kTRUE;
}

// ---------------------------------------------------------------
 Bool_t JSFReadGeneratorBuf::ReadOneRecord()
{
// Read Generator data and saved to the class.

// Due to imcompatibility in JSFGeneratorParticle class and HEPEVT,
// second mother information is not saved.
//
// 

  const Int_t kNMXHEP=4000;
  Int_t nevhep, nhep;
  Int_t isthep[kNMXHEP], idhep[kNMXHEP];
  Int_t jmohep[kNMXHEP][2], jdahep[kNMXHEP][2];
  Double_t phep[kNMXHEP][5], vhep[kNMXHEP][4];


  if( !ReadHepEvent(kNMXHEP, nevhep, nhep, isthep, idhep, jmohep, jdahep,
		    phep, vhep) ) return kFALSE;

  // ***************************************
  // Scan data and reassign mother-daughter map
  // ***************************************

  Int_t i;
  struct modau {
    Int_t mother;
    Int_t ndau;
    Int_t dau1st;
    Int_t nser;
  } map[kNMXHEP];
  Int_t nser=0;
  Int_t n1stfinal=0;
  Int_t n1stdoc=0;
  Int_t jlist[kNMXHEP];

  for(i=0;i<nhep;i++){
    map[i].nser=0;
    if( isthep[i] == 0 ) continue;

    nser++;
    jlist[nser]=i;
    /*
    if( jmohep[i][0] < 1 && (isthep[i]==1 || isthep[i]==2 )) {
      printf("Fatal error in JSFReadGenerator::ReadOneRecord().n");
      printf("jmohep should not be 0 when isthep=1 or 2 but n");
      printf("jmohep[%d][0]=%d, isthep[%d]=%dn",i,jmohep[i][0],i,isthep[i]);
      return kFALSE;
    }
    */

    switch ( isthep[i] ) {
      case 1:
        map[i].ndau=0;
        map[i].dau1st=0;
	if( jmohep[i][0] <= 0 ) map[i].mother=0;
        else map[i].mother=map[jmohep[i][0]-1].nser; // modified by I.Nakamura
        map[i].nser=nser;
	if( n1stfinal == 0 ) n1stfinal=nser;
	if( map[i].mother == 0 ) map[i].mother=n1stdoc;
	break;
      case 2:
        map[i].ndau=jdahep[i][1]-jdahep[i][0]+1;
        map[i].dau1st=-2;
	if( jmohep[i][0] <= 0 ) map[i].mother=0;
        else map[i].mother=map[jmohep[i][0]-1].nser; // modified by I.Nakamura
        map[i].nser=nser;
	if( n1stfinal == 0 ) n1stfinal=nser;
	if( map[i].mother == 0 ) map[i].mother=n1stdoc;
	break;
      case 3:
	if( n1stdoc == 0 ) n1stdoc=nser;
        map[i].ndau=1;
        map[i].dau1st=-3;
        map[i].mother=-3;
        map[i].nser=nser;
	break;
      default:
        map[i].ndau=1;
        map[i].dau1st=-4;
        map[i].mother=-isthep[i];
        map[i].nser=nser;
    }
  }

  for(i=0;i<nhep;i++){
    if( map[i].nser == 0 ) continue;
    switch (map[i].dau1st){
      case -2:
	map[i].dau1st=map[jdahep[i][0]-1].nser; // modified by I.Nakamura
	break;
      case -3:
      case -4:
	map[i].dau1st=n1stfinal;
	break;
    }
  }

  // ***************************************
  // Fill GeneratorParticle Array
  // ***************************************
  
  TClonesArray &particles=*fParticles;

  TVector p(4);
  TVector v(4);
  Int_t iser=0;
  
  Int_t is;
  for(is=1;is<=nser;is++){
    i=jlist[is];
    if( isthep[i] == 0 ) continue;

    Int_t id=idhep[i];
    Float_t mass=phep[i][4];
    Int_t mother=map[i].mother;
    Int_t firstdaughter=map[i].dau1st;
    Int_t ndaughter=map[i].ndau;

#if __PYTHIA_VERSION__ >= 6
    Float_t charge=((Float_t)pychge_(&id))/3.0;
    Int_t kc=pycomp_(&id);
    Float_t xctau=pydat2_.PMAS[3][kc-1]*0.1;
#else
    Float_t charge=((Float_t)luchge_(&id))/3.0;
    Float_t xctau=((Float_t)ulctau_(&id));
#endif
    Float_t dl=0.0;
    //    if( mother < 0 ) xctau=0.0;

    p(0)=phep[i][3]; p(1)=phep[i][0]; p(2)=phep[i][1]; p(3)=phep[i][2];
    v(0)=vhep[i][3]; v(1)=vhep[i][0]; v(2)=vhep[i][1]; v(3)=vhep[i][2];

    iser++;
    new(particles[is-1]) JSFGeneratorParticle(iser, id, mass, charge,
	              p, v, ndaughter, firstdaughter, mother,
		      xctau, dl);
  }

  SetNparticles(iser);

  return kTRUE;
}

#if __ROOT_FULLVERSION >= 30000
//______________________________________________________________________________
void JSFReadGenerator::Streamer(TBuffer &R__b)
{
   // Stream an object of class JSFReadGenerator.

   if (R__b.IsReading()) {
     UInt_t R__s, R__c;
     Version_t R__v=R__b.ReadVersion(&R__s, &R__c);
     if( R__v > 1 ) {
       JSFReadGenerator::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
       return;
     }
     JSFGenerator::Streamer(R__b);
     R__b >> fUnit;
     R__b.ReadStaticArray(fFormat);
     R__b.CheckByteCount(R__s, R__c, JSFReadGenerator::IsA());

   } else {
     JSFReadGenerator::Class()->WriteBuffer(R__b, this);
   }

}

#else
//______________________________________________________________________________
void JSFReadGenerator::Streamer(TBuffer &R__b)
{
   // Stream an object of class JSFReadGenerator.

   if (R__b.IsReading()) {
      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
      JSFGenerator::Streamer(R__b);
      R__b >> fUnit;
      R__b.ReadStaticArray(fFormat);
   } else {
      R__b.WriteVersion(JSFReadGenerator::IsA());
      JSFGenerator::Streamer(R__b);
      R__b << fUnit;
      R__b.WriteArray(fFormat, 32);
   }
}

#endif