[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 "AliGenHerwigEventHeader.h"
#include "AliHerwigRndm.h"
#include "AliMC.h"
#include "AliRun.h"
#include "driver.h"

using std::cerr;
using std::endl;

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),
    fFileName(0),
    fEtaMinParton(-20.),     
    fEtaMaxParton(20.),     
    fPhiMinParton(0.),     
    fPhiMaxParton(2.* TMath::Pi()),     
    fEtaMinGamma(-20.),      
    fEtaMaxGamma(20.),      
    fPhiMinGamma(0.),      
    fPhiMaxGamma(2. * TMath::Pi()),
    fHeader(0)
{
// 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),
    fFileName(0),
    fEtaMinParton(-20.),     
    fEtaMaxParton(20.),     
    fPhiMinParton(0.),     
    fPhiMaxParton(2.* TMath::Pi()),     
    fEtaMinGamma(-20.),      
    fEtaMaxGamma(20.),      
    fPhiMinGamma(0.),
    fPhiMaxGamma(2. * TMath::Pi()),
    fHeader(0)
{
// Constructor
    fEnergyCMS = 14000;
    SetTarget();
    SetProjectile();
    // Set random number generator
    AliHerwigRndm::SetHerwigRandom(GetRandom());
}

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

void AliGenHerwig::Init()
{
// Initialisation
  fTarget.Resize(8);
  fProjectile.Resize(8);
  SetMC(new THerwig6());
  fHerwig=(THerwig6*) fMCEvGen;
  // initialize common blocks
  fHerwig->Initialize(fProjectile.Data(), fTarget.Data(), fMomentum1, fMomentum2, fProcess); // 
  // reset parameters according to user needs
  InitPDF();
  fHerwig->SetPTMIN(fPtHardMin);
  fHerwig->SetPTMAX(fPtHardMax);
  fHerwig->SetPTRMS(fPtRMS);
  printf("SetMAXPR %15d \n", fMaxPr);
  fHerwig->SetMAXPR(fMaxPr);
  fHerwig->SetMAXER(fMaxErrors);
  fHerwig->SetENSOF(fEnSoft);
// 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(), 49);

  //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);

// 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(), 49);

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

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

void AliGenHerwig::InitPDF()
{
// Initialize PDF
  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 p[4];

  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 kf, ks, imo;
  kf=0;

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

  fTrials=0;

  //  Set collision vertex position 
  if (fVertexSmear == kPerEvent) Vertex();

  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;
	fNprimaries = 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] = fVertex[0] + iparticle->Vx()/10; // [cm]
		origin[1] = fVertex[1] + iparticle->Vy()/10; // [cm]
		origin[2] = fVertex[2] + iparticle->Vz()/10; // [cm]

		Float_t tof = fTime + 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;
		fNprimaries++;
	    } // end of if: selection of particle
	} // end of for: particle loop
	if (newPos) delete[] newPos;
	if (nc > 0) {
	    jev+=nc;
	    if (jev >= fNpart || fNpart == -1) {
		fKineBias=Float_t(fNpart)/Float_t(fTrials);
		break;
	    }
	}
    }
//
  MakeHeader();
//  
  SetHighWaterMark(nt);
//  adjust weight due to kinematic selection
  AdjustWeights();
//  get cross-section
  fXsection=fHerwig->GetAVWGT();
  //printf(">> trials << %d\n",fTrials);
}

Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const 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(const TParticle* iparticle, const 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) const
{
// 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(const TParticle*  particle) const
{
// 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();

}


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