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