[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */


// Generator using Herwig as an external generator
// The main Herwig options are accessable for the user through this interface.
// Uses the THerwig implementation of TGenerator.


#include <Riostream.h>
#include <TClonesArray.h>
#include <TParticle.h>

#include <THerwig6.h>

#include "AliGenHerwig.h"
#include "AliHerwigRndm.h"
#include "AliMC.h"
#include "AliRun.h"
#include "driver.h"

ClassImp(AliGenHerwig)


  AliGenHerwig::AliGenHerwig() :
    AliGenMC(),
    fAutPDF("LHAPDF"),
    fModPDF(19070),
    fStrucFunc(kCTEQ5L),
    fKeep(0),
    fDecaysOff(1),
    fTrigger(0),
    fSelectAll(0),
    fFlavor(0),
    fMomentum1(7000),
    fMomentum2(7000),
    fKineBias(1),
    fTrials(0),
    fXsection(0),
    fHerwig(0x0),
    fProcess(0),
    fPtHardMin(0.),
    fPtHardMax(9999.),
    fPtRMS(0.),
    fMaxPr(10),
    fMaxErrors(1000),
    fEnSoft(1),
    fEv1Pr(0),
    fEv2Pr(0),
    fFileName(0),
    fEtaMinParton(-20.),     
    fEtaMaxParton(20.),     
    fPhiMinParton(0.),     
    fPhiMaxParton(2.* TMath::Pi()),     
    fEtaMinGamma(-20.),      
    fEtaMaxGamma(20.),      
    fPhiMinGamma(0.),      
    fPhiMaxGamma(2. * TMath::Pi())  
{
// Constructor
  fEnergyCMS = 14000;
}

AliGenHerwig::AliGenHerwig(Int_t npart)
    :AliGenMC(npart),
    fAutPDF("LHAPDF"),
    fModPDF(19070),
    fStrucFunc(kCTEQ5L),
    fKeep(0),
    fDecaysOff(1),
    fTrigger(0),
    fSelectAll(0),
    fFlavor(0),
    fMomentum1(7000),
    fMomentum2(7000),
    fKineBias(1),
    fTrials(0),
    fXsection(0),
    fHerwig(0x0),
    fProcess(0),
    fPtHardMin(0.),
    fPtHardMax(9999.),
    fPtRMS(0.),
    fMaxPr(10),
    fMaxErrors(1000),
    fEnSoft(1),
    fEv1Pr(0),
    fEv2Pr(0),
    fFileName(0),
    fEtaMinParton(-20.),     
    fEtaMaxParton(20.),     
    fPhiMinParton(0.),     
    fPhiMaxParton(2.* TMath::Pi()),     
    fEtaMinGamma(-20.),      
    fEtaMaxGamma(20.),      
    fPhiMinGamma(0.),      
    fPhiMaxGamma(2. * TMath::Pi())  
{
    fEnergyCMS = 14000;
    SetTarget();
    SetProjectile();
    // Set random number generator
    AliHerwigRndm::SetHerwigRandom(GetRandom());
}

AliGenHerwig::~AliGenHerwig()
{
// Destructor
}

void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
{
  fEv1Pr = ++eventFirst;
  fEv2Pr = ++eventLast;
  if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
}

void AliGenHerwig::Init()
{
// Initialisation
  fTarget.Resize(8);
  fProjectile.Resize(8);
  SetMC(new THerwig6());
  fHerwig=(THerwig6*) fMCEvGen;
  // initialize common blocks
  fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
  // reset parameters according to user needs
  InitPDF();
  fHerwig->SetPTMIN(fPtHardMin);
  fHerwig->SetPTRMS(fPtRMS);
  fHerwig->SetMAXPR(fMaxPr);
  fHerwig->SetMAXER(fMaxErrors);
  fHerwig->SetENSOF(fEnSoft);

  fHerwig->SetEV1PR(fEv1Pr);
  fHerwig->SetEV2PR(fEv2Pr);

// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
//       RMASS(1)=0.32
//       RMASS(2)=0.32
//       RMASS(3)=0.5
//       RMASS(4)=1.55
//       RMASS(5)=4.75
//       RMASS(6)=174.3
//       RMASS(13)=0.75

  fHerwig->SetRMASS(4,1.2);
  fHerwig->SetRMASS(5,4.75);

  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);

  //fHerwig->Hwusta("PI0     ");

  // compute parameter dependent constants
  fHerwig->PrepareRun();
}

void AliGenHerwig::InitJimmy()
{
// Initialisation
  fTarget.Resize(8);
  fProjectile.Resize(8);
  SetMC(new THerwig6());
  fHerwig=(THerwig6*) fMCEvGen;
  // initialize common blocks
  fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
  // reset parameters according to user needs
  InitPDF();
  fHerwig->SetPTMIN(fPtHardMin);
  fHerwig->SetPTRMS(fPtRMS);
  fHerwig->SetMAXPR(fMaxPr);
  fHerwig->SetMAXER(fMaxErrors);
  fHerwig->SetENSOF(fEnSoft);

  fHerwig->SetEV1PR(fEv1Pr);
  fHerwig->SetEV2PR(fEv2Pr);

// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
//       RMASS(1)=0.32
//       RMASS(2)=0.32
//       RMASS(3)=0.5
//       RMASS(4)=1.55
//       RMASS(5)=4.75
//       RMASS(6)=174.3
//       RMASS(13)=0.75

  fHerwig->SetRMASS(4,1.2);
  fHerwig->SetRMASS(5,4.75);

  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);

  //  fHerwig->Hwusta("PI0     ");

  // compute parameter dependent constants
  fHerwig->PrepareRunJimmy();
}

void AliGenHerwig::InitPDF()
{
  switch(fStrucFunc)
    {
// ONLY USES LHAPDF STRUCTURE FUNCTIONS
    case kGRVLO98:
      fModPDF=80060;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ6:
      fModPDF=10040;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ61:
      fModPDF=10100;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ6m:
      fModPDF=10050;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ6l:
      fModPDF=10041;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ6ll:
      fModPDF=10042;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ5M:
      fModPDF=19050;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ5L:
      fModPDF=19070;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ4M:
      fModPDF=19150;
      fAutPDF="HWLHAPDF";
      break;
    case kCTEQ4L:
      fModPDF=19170;
      fAutPDF="HWLHAPDF";
      break;
//    case kMRST2004nlo:
//      fModPDF=20400;
//      fAutPDF="HWLHAPDF";
//      break;
    default:
      cerr << "This structure function is not inplemented " << fStrucFunc << endl;
      break;
    }
  fAutPDF.Resize(20);
  fHerwig->SetMODPDF(1,fModPDF);
  fHerwig->SetMODPDF(2,fModPDF);
  fHerwig->SetAUTPDF(1,fAutPDF);
  fHerwig->SetAUTPDF(2,fAutPDF);
}

void AliGenHerwig::Generate()
{
  // Generate one event

  Float_t polar[3] =   {0,0,0};
  Float_t origin[3]=   {0,0,0};
  Float_t origin0[3]=  {0,0,0};
  Float_t p[4], random[6];

  static TClonesArray *particles;
  //  converts from mm/c to s
  const Float_t kconv=0.001/2.999792458e8;
  //
  Int_t nt=0;
  Int_t jev=0;
  Int_t j, kf, ks, imo;
  kf=0;

  if(!particles) particles=new TClonesArray("TParticle",10000);

  fTrials=0;
  for (j=0;j<3;j++) origin0[j]=fOrigin[j];
  if(fVertexSmear==kPerEvent) {
    Rndm(random,6);
    for (j=0;j<3;j++) {
      origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
	TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
    }
  }

  while(1)
    {
	fHerwig->GenerateEvent();
	fTrials++;
	fHerwig->ImportParticles(particles,"All");
	Int_t np = particles->GetEntriesFast()-1;
	if (np == 0 ) continue;

	//Check hard partons or direct gamma in kine range

	if (fProcess == kHeJets || fProcess == kHeDirectGamma) {
	  TParticle* parton1 = (TParticle *) particles->At(6);
	  TParticle* parton2 = (TParticle *) particles->At(7);
	  if (!CheckParton(parton1, parton2))  continue ;
	} 

	Int_t nc=0;

	Int_t * newPos = new Int_t[np];
	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;

	for (Int_t i = 0; i<np; i++) {
	    TParticle *  iparticle       = (TParticle *) particles->At(i);
	    imo = iparticle->GetFirstMother();
	    kf        = iparticle->GetPdgCode();
	    ks        = iparticle->GetStatusCode();
	    if (ks != 3 &&
		KinematicSelection(iparticle,0))
	    {
		nc++;
		p[0]=iparticle->Px();
		p[1]=iparticle->Py();
		p[2]=iparticle->Pz();
		p[3]=iparticle->Energy();
		origin[0]=origin0[0]+iparticle->Vx()/10;
		origin[1]=origin0[1]+iparticle->Vy()/10;
		origin[2]=origin0[2]+iparticle->Vz()/10;
		Float_t tof = kconv*iparticle->T();
		Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
		Int_t   trackIt = (ks == 1) && fTrackIt;
		PushTrack(trackIt, iparent, kf,
			  p[0], p[1], p[2], p[3],
			  origin[0], origin[1], origin[2],
			  tof,
			  polar[0], polar[1], polar[2],
			  kPPrimary, nt, fHerwig->GetEVWGT(), ks);
		KeepTrack(nt);
		newPos[i]=nt;
	    } // end of if: selection of particle
	} // end of for: particle loop
	if (newPos) delete[] newPos;
	//      MakeHeader();
	if (nc > 0) {
	    jev+=nc;
	    if (jev >= fNpart || fNpart == -1) {
		fKineBias=Float_t(fNpart)/Float_t(fTrials);
		break;
	    }
	}
    }
  SetHighWaterMark(nt);
//  adjust weight due to kinematic selection
  AdjustWeights();
//  get cross-section
  fXsection=fHerwig->GetAVWGT();
  //printf(">> trials << %d\n",fTrials);
}

Bool_t AliGenHerwig::CheckParton(TParticle* parton1, TParticle* parton2)
{
// Check the kinematic trigger condition
//
//Select events with parton max energy
    if(fPtHardMax < parton1->Pt()) return kFALSE;

// Select events within angular window
    Double_t eta[2];
    eta[0] = parton1->Eta();
    eta[1] = parton2->Eta();
    Double_t phi[2];
    phi[0] = parton1->Phi();
    phi[1] = parton2->Phi();
    Int_t    pdg[2]; 
    pdg[0] = parton1->GetPdgCode();
    pdg[1] = parton2->GetPdgCode();   
    printf("min %f, max %f\n",fPtHardMin, fPtHardMax);
    printf("Parton 1: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton1->GetName(),parton1->Pt(), eta[0],phi[0]*TMath::RadToDeg());
    printf("Parton 2: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton2->GetName(),parton2->Pt(), eta[1],phi[1]*TMath::RadToDeg());
    
    if (fProcess == kHeJets) {
      //Check if one of the 2 outgoing partons are in the eta-phi window
      for(Int_t i = 0; i < 2; i++)
	if ((eta[i] < fEtaMaxParton  && eta[i] > fEtaMinParton) &&
	    (phi[i] < fPhiMaxParton  && phi[i] > fPhiMinParton)) return  kTRUE ;
    }
    
    else {
      //Check if the gamma and the jet  are in the eta-phi window
      Int_t igj = 0;
      Int_t ijj = 0;
      if(pdg[0] == 22) ijj=1;
      else igj=1;
      if ((eta[ijj] < fEtaMaxParton   && eta[ijj] > fEtaMinParton) &&
	  (phi[ijj] < fPhiMaxParton   && phi[ijj] > fPhiMinParton)) {
	
	if ((eta[igj] < fEtaMaxGamma   && eta[igj] > fEtaMinGamma) &&
	    (phi[igj] < fPhiMaxGamma   && phi[igj] > fPhiMinGamma)) return  kTRUE;
	
      }
    }

    return kFALSE ;
}

void AliGenHerwig::AdjustWeights()
{
// Adjust the weights after generation of all events
    TParticle *part;
    Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
    for (Int_t i=0; i<ntrack; i++) {
        part= gAlice->GetMCApp()->Particle(i);
        part->SetWeight(part->GetWeight()*fKineBias);
    }
}


void AliGenHerwig::KeepFullEvent()
{
    fKeep=1;
}

Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
{
//
// Looks recursively if one of the daughters has been selected
//
//    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
    Int_t imin=-1;
    Int_t imax=-1;
    Int_t i;
    Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
    Bool_t selected=kFALSE;
    if (hasDaughters) {
	imin=iparticle->GetFirstDaughter();
	imax=iparticle->GetLastDaughter();
	for (i=imin; i<= imax; i++){
	    TParticle *  jparticle       = (TParticle *) particles->At(i);
	    Int_t ip=jparticle->GetPdgCode();
	    if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
		selected=kTRUE; break;
	    }
	    if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	}
    } else {
	return kFALSE;
    }

    return selected;
}


Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
{
// Select flavor of particle
// 0: all
// 4: charm and beauty
// 5: beauty
    if (fFlavor == 0) return kTRUE;

    Int_t ifl=TMath::Abs(pid/100);
    if (ifl > 10) ifl/=10;
    return (fFlavor == ifl);
}

Bool_t AliGenHerwig::Stable(TParticle*  particle)
{
// Return true for a stable particle
//
    Int_t kf = TMath::Abs(particle->GetPdgCode());

    if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))

    {
	return kTRUE;
    } else {
	return kFALSE;
    }
}

void AliGenHerwig::FinishRun()
{
  fHerwig->Hwefin();
}

void AliGenHerwig::FinishRunJimmy()
{
  fHerwig->Hwefin();
  fHerwig->Jmefin();

}
Commit	Line	Data
	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18
	19	// Generator using Herwig as an external generator
	20	// The main Herwig options are accessable for the user through this interface.
	21	// Uses the THerwig implementation of TGenerator.
	22
	23
	24	#include <Riostream.h>
	25	#include <TClonesArray.h>
	26	#include <TParticle.h>
	27
	28	#include <THerwig6.h>
	29
	30	#include "AliGenHerwig.h"
	31	#include "AliHerwigRndm.h"
	32	#include "AliMC.h"
	33	#include "AliRun.h"
	34	#include "driver.h"
	35
	36	ClassImp(AliGenHerwig)
	37
	38
	39	AliGenHerwig::AliGenHerwig() :
	40	AliGenMC(),
	41	fAutPDF("LHAPDF"),
	42	fModPDF(19070),
	43	fStrucFunc(kCTEQ5L),
	44	fKeep(0),
	45	fDecaysOff(1),
	46	fTrigger(0),
	47	fSelectAll(0),
	48	fFlavor(0),
	49	fMomentum1(7000),
	50	fMomentum2(7000),
	51	fKineBias(1),
	52	fTrials(0),
	53	fXsection(0),
	54	fHerwig(0x0),
	55	fProcess(0),
	56	fPtHardMin(0.),
	57	fPtHardMax(9999.),
	58	fPtRMS(0.),
	59	fMaxPr(10),
	60	fMaxErrors(1000),
	61	fEnSoft(1),
	62	fEv1Pr(0),
	63	fEv2Pr(0),
	64	fFileName(0),
	65	fEtaMinParton(-20.),
	66	fEtaMaxParton(20.),
	67	fPhiMinParton(0.),
	68	fPhiMaxParton(2.* TMath::Pi()),
	69	fEtaMinGamma(-20.),
	70	fEtaMaxGamma(20.),
	71	fPhiMinGamma(0.),
	72	fPhiMaxGamma(2. * TMath::Pi())
	73	{
	74	// Constructor
	75	fEnergyCMS = 14000;
	76	}
	77
	78	AliGenHerwig::AliGenHerwig(Int_t npart)
	79	:AliGenMC(npart),
	80	fAutPDF("LHAPDF"),
	81	fModPDF(19070),
	82	fStrucFunc(kCTEQ5L),
	83	fKeep(0),
	84	fDecaysOff(1),
	85	fTrigger(0),
	86	fSelectAll(0),
	87	fFlavor(0),
	88	fMomentum1(7000),
	89	fMomentum2(7000),
	90	fKineBias(1),
	91	fTrials(0),
	92	fXsection(0),
	93	fHerwig(0x0),
	94	fProcess(0),
	95	fPtHardMin(0.),
	96	fPtHardMax(9999.),
	97	fPtRMS(0.),
	98	fMaxPr(10),
	99	fMaxErrors(1000),
	100	fEnSoft(1),
	101	fEv1Pr(0),
	102	fEv2Pr(0),
	103	fFileName(0),
	104	fEtaMinParton(-20.),
	105	fEtaMaxParton(20.),
	106	fPhiMinParton(0.),
	107	fPhiMaxParton(2.* TMath::Pi()),
	108	fEtaMinGamma(-20.),
	109	fEtaMaxGamma(20.),
	110	fPhiMinGamma(0.),
	111	fPhiMaxGamma(2. * TMath::Pi())
	112	{
	113	fEnergyCMS = 14000;
	114	SetTarget();
	115	SetProjectile();
	116	// Set random number generator
	117	AliHerwigRndm::SetHerwigRandom(GetRandom());
	118	}
	119
	120	AliGenHerwig::~AliGenHerwig()
	121	{
	122	// Destructor
	123	}
	124
	125	void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
	126	{
	127	fEv1Pr = ++eventFirst;
	128	fEv2Pr = ++eventLast;
	129	if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
	130	}
	131
	132	void AliGenHerwig::Init()
	133	{
	134	// Initialisation
	135	fTarget.Resize(8);
	136	fProjectile.Resize(8);
	137	SetMC(new THerwig6());
	138	fHerwig=(THerwig6*) fMCEvGen;
	139	// initialize common blocks
	140	fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
	141	// reset parameters according to user needs
	142	InitPDF();
	143	fHerwig->SetPTMIN(fPtHardMin);
	144	fHerwig->SetPTRMS(fPtRMS);
	145	fHerwig->SetMAXPR(fMaxPr);
	146	fHerwig->SetMAXER(fMaxErrors);
	147	fHerwig->SetENSOF(fEnSoft);
	148
	149	fHerwig->SetEV1PR(fEv1Pr);
	150	fHerwig->SetEV2PR(fEv2Pr);
	151
	152	// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
	153	// RMASS(1)=0.32
	154	// RMASS(2)=0.32
	155	// RMASS(3)=0.5
	156	// RMASS(4)=1.55
	157	// RMASS(5)=4.75
	158	// RMASS(6)=174.3
	159	// RMASS(13)=0.75
	160
	161	fHerwig->SetRMASS(4,1.2);
	162	fHerwig->SetRMASS(5,4.75);
	163
	164	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
	165
	166	//fHerwig->Hwusta("PI0 ");
	167
	168	// compute parameter dependent constants
	169	fHerwig->PrepareRun();
	170	}
	171
	172	void AliGenHerwig::InitJimmy()
	173	{
	174	// Initialisation
	175	fTarget.Resize(8);
	176	fProjectile.Resize(8);
	177	SetMC(new THerwig6());
	178	fHerwig=(THerwig6*) fMCEvGen;
	179	// initialize common blocks
	180	fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
	181	// reset parameters according to user needs
	182	InitPDF();
	183	fHerwig->SetPTMIN(fPtHardMin);
	184	fHerwig->SetPTRMS(fPtRMS);
	185	fHerwig->SetMAXPR(fMaxPr);
	186	fHerwig->SetMAXER(fMaxErrors);
	187	fHerwig->SetENSOF(fEnSoft);
	188
	189	fHerwig->SetEV1PR(fEv1Pr);
	190	fHerwig->SetEV2PR(fEv2Pr);
	191
	192	// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
	193	// RMASS(1)=0.32
	194	// RMASS(2)=0.32
	195	// RMASS(3)=0.5
	196	// RMASS(4)=1.55
	197	// RMASS(5)=4.75
	198	// RMASS(6)=174.3
	199	// RMASS(13)=0.75
	200
	201	fHerwig->SetRMASS(4,1.2);
	202	fHerwig->SetRMASS(5,4.75);
	203
	204	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
	205
	206	// fHerwig->Hwusta("PI0 ");
	207
	208	// compute parameter dependent constants
	209	fHerwig->PrepareRunJimmy();
	210	}
	211
	212	void AliGenHerwig::InitPDF()
	213	{
	214	switch(fStrucFunc)
	215	{
	216	// ONLY USES LHAPDF STRUCTURE FUNCTIONS
	217	case kGRVLO98:
	218	fModPDF=80060;
	219	fAutPDF="HWLHAPDF";
	220	break;
	221	case kCTEQ6:
	222	fModPDF=10040;
	223	fAutPDF="HWLHAPDF";
	224	break;
	225	case kCTEQ61:
	226	fModPDF=10100;
	227	fAutPDF="HWLHAPDF";
	228	break;
	229	case kCTEQ6m:
	230	fModPDF=10050;
	231	fAutPDF="HWLHAPDF";
	232	break;
	233	case kCTEQ6l:
	234	fModPDF=10041;
	235	fAutPDF="HWLHAPDF";
	236	break;
	237	case kCTEQ6ll:
	238	fModPDF=10042;
	239	fAutPDF="HWLHAPDF";
	240	break;
	241	case kCTEQ5M:
	242	fModPDF=19050;
	243	fAutPDF="HWLHAPDF";
	244	break;
	245	case kCTEQ5L:
	246	fModPDF=19070;
	247	fAutPDF="HWLHAPDF";
	248	break;
	249	case kCTEQ4M:
	250	fModPDF=19150;
	251	fAutPDF="HWLHAPDF";
	252	break;
	253	case kCTEQ4L:
	254	fModPDF=19170;
	255	fAutPDF="HWLHAPDF";
	256	break;
	257	// case kMRST2004nlo:
	258	// fModPDF=20400;
	259	// fAutPDF="HWLHAPDF";
	260	// break;
	261	default:
	262	cerr << "This structure function is not inplemented " << fStrucFunc << endl;
	263	break;
	264	}
	265	fAutPDF.Resize(20);
	266	fHerwig->SetMODPDF(1,fModPDF);
	267	fHerwig->SetMODPDF(2,fModPDF);
	268	fHerwig->SetAUTPDF(1,fAutPDF);
	269	fHerwig->SetAUTPDF(2,fAutPDF);
	270	}
	271
	272	void AliGenHerwig::Generate()
	273	{
	274	// Generate one event
	275
	276	Float_t polar[3] = {0,0,0};
	277	Float_t origin[3]= {0,0,0};
	278	Float_t origin0[3]= {0,0,0};
	279	Float_t p[4], random[6];
	280
	281	static TClonesArray *particles;
	282	// converts from mm/c to s
	283	const Float_t kconv=0.001/2.999792458e8;
	284	//
	285	Int_t nt=0;
	286	Int_t jev=0;
	287	Int_t j, kf, ks, imo;
	288	kf=0;
	289
	290	if(!particles) particles=new TClonesArray("TParticle",10000);
	291
	292	fTrials=0;
	293	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
	294	if(fVertexSmear==kPerEvent) {
	295	Rndm(random,6);
	296	for (j=0;j<3;j++) {
	297	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	298	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	299	}
	300	}
	301
	302	while(1)
	303	{
	304	fHerwig->GenerateEvent();
	305	fTrials++;
	306	fHerwig->ImportParticles(particles,"All");
	307	Int_t np = particles->GetEntriesFast()-1;
	308	if (np == 0 ) continue;
	309
	310	//Check hard partons or direct gamma in kine range
	311
	312	if (fProcess == kHeJets \|\| fProcess == kHeDirectGamma) {
	313	TParticle* parton1 = (TParticle *) particles->At(6);
	314	TParticle* parton2 = (TParticle *) particles->At(7);
	315	if (!CheckParton(parton1, parton2)) continue ;
	316	}
	317
	318	Int_t nc=0;
	319
	320	Int_t * newPos = new Int_t[np];
	321	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
	322
	323	for (Int_t i = 0; i<np; i++) {
	324	TParticle * iparticle = (TParticle *) particles->At(i);
	325	imo = iparticle->GetFirstMother();
	326	kf = iparticle->GetPdgCode();
	327	ks = iparticle->GetStatusCode();
	328	if (ks != 3 &&
	329	KinematicSelection(iparticle,0))
	330	{
	331	nc++;
	332	p[0]=iparticle->Px();
	333	p[1]=iparticle->Py();
	334	p[2]=iparticle->Pz();
	335	p[3]=iparticle->Energy();
	336	origin[0]=origin0[0]+iparticle->Vx()/10;
	337	origin[1]=origin0[1]+iparticle->Vy()/10;
	338	origin[2]=origin0[2]+iparticle->Vz()/10;
	339	Float_t tof = kconv*iparticle->T();
	340	Int_t iparent = (imo > -1) ? newPos[imo] : -1;
	341	Int_t trackIt = (ks == 1) && fTrackIt;
	342	PushTrack(trackIt, iparent, kf,
	343	p[0], p[1], p[2], p[3],
	344	origin[0], origin[1], origin[2],
	345	tof,
	346	polar[0], polar[1], polar[2],
	347	kPPrimary, nt, fHerwig->GetEVWGT(), ks);
	348	KeepTrack(nt);
	349	newPos[i]=nt;
	350	} // end of if: selection of particle
	351	} // end of for: particle loop
	352	if (newPos) delete[] newPos;
	353	// MakeHeader();
	354	if (nc > 0) {
	355	jev+=nc;
	356	if (jev >= fNpart \|\| fNpart == -1) {
	357	fKineBias=Float_t(fNpart)/Float_t(fTrials);
	358	break;
	359	}
	360	}
	361	}
	362	SetHighWaterMark(nt);
	363	// adjust weight due to kinematic selection
	364	AdjustWeights();
	365	// get cross-section
	366	fXsection=fHerwig->GetAVWGT();
	367	//printf(">> trials << %d\n",fTrials);
	368	}
	369
	370	Bool_t AliGenHerwig::CheckParton(TParticle* parton1, TParticle* parton2)
	371	{
	372	// Check the kinematic trigger condition
	373	//
	374	//Select events with parton max energy
	375	if(fPtHardMax < parton1->Pt()) return kFALSE;
	376
	377	// Select events within angular window
	378	Double_t eta[2];
	379	eta[0] = parton1->Eta();
	380	eta[1] = parton2->Eta();
	381	Double_t phi[2];
	382	phi[0] = parton1->Phi();
	383	phi[1] = parton2->Phi();
	384	Int_t pdg[2];
	385	pdg[0] = parton1->GetPdgCode();
	386	pdg[1] = parton2->GetPdgCode();
	387	printf("min %f, max %f\n",fPtHardMin, fPtHardMax);
	388	printf("Parton 1: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton1->GetName(),parton1->Pt(), eta[0],phi[0]*TMath::RadToDeg());
	389	printf("Parton 2: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton2->GetName(),parton2->Pt(), eta[1],phi[1]*TMath::RadToDeg());
	390
	391	if (fProcess == kHeJets) {
	392	//Check if one of the 2 outgoing partons are in the eta-phi window
	393	for(Int_t i = 0; i < 2; i++)
	394	if ((eta[i] < fEtaMaxParton && eta[i] > fEtaMinParton) &&
	395	(phi[i] < fPhiMaxParton && phi[i] > fPhiMinParton)) return kTRUE ;
	396	}
	397
	398	else {
	399	//Check if the gamma and the jet are in the eta-phi window
	400	Int_t igj = 0;
	401	Int_t ijj = 0;
	402	if(pdg[0] == 22) ijj=1;
	403	else igj=1;
	404	if ((eta[ijj] < fEtaMaxParton && eta[ijj] > fEtaMinParton) &&
	405	(phi[ijj] < fPhiMaxParton && phi[ijj] > fPhiMinParton)) {
	406
	407	if ((eta[igj] < fEtaMaxGamma && eta[igj] > fEtaMinGamma) &&
	408	(phi[igj] < fPhiMaxGamma && phi[igj] > fPhiMinGamma)) return kTRUE;
	409
	410	}
	411	}
	412
	413	return kFALSE ;
	414	}
	415
	416	void AliGenHerwig::AdjustWeights()
	417	{
	418	// Adjust the weights after generation of all events
	419	TParticle *part;
	420	Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
	421	for (Int_t i=0; i<ntrack; i++) {
	422	part= gAlice->GetMCApp()->Particle(i);
	423	part->SetWeight(part->GetWeight()*fKineBias);
	424	}
	425	}
	426
	427
	428	void AliGenHerwig::KeepFullEvent()
	429	{
	430	fKeep=1;
	431	}
	432
	433	Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
	434	{
	435	//
	436	// Looks recursively if one of the daughters has been selected
	437	//
	438	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	439	Int_t imin=-1;
	440	Int_t imax=-1;
	441	Int_t i;
	442	Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
	443	Bool_t selected=kFALSE;
	444	if (hasDaughters) {
	445	imin=iparticle->GetFirstDaughter();
	446	imax=iparticle->GetLastDaughter();
	447	for (i=imin; i<= imax; i++){
	448	TParticle * jparticle = (TParticle *) particles->At(i);
	449	Int_t ip=jparticle->GetPdgCode();
	450	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	451	selected=kTRUE; break;
	452	}
	453	if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	454	}
	455	} else {
	456	return kFALSE;
	457	}
	458
	459	return selected;
	460	}
	461
	462
	463	Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
	464	{
	465	// Select flavor of particle
	466	// 0: all
	467	// 4: charm and beauty
	468	// 5: beauty
	469	if (fFlavor == 0) return kTRUE;
	470
	471	Int_t ifl=TMath::Abs(pid/100);
	472	if (ifl > 10) ifl/=10;
	473	return (fFlavor == ifl);
	474	}
	475
	476	Bool_t AliGenHerwig::Stable(TParticle* particle)
	477	{
	478	// Return true for a stable particle
	479	//
	480	Int_t kf = TMath::Abs(particle->GetPdgCode());
	481
	482	if ( (particle->GetFirstDaughter() < 0 ) \|\| (kf == 1000*fFlavor+122))
	483
	484	{
	485	return kTRUE;
	486	} else {
	487	return kFALSE;
	488	}
	489	}
	490
	491	void AliGenHerwig::FinishRun()
	492	{
	493	fHerwig->Hwefin();
	494	}
	495
	496	void AliGenHerwig::FinishRunJimmy()
	497	{
	498	fHerwig->Hwefin();
	499	fHerwig->Jmefin();
	500
	501	}
	502