[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
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);
	144	fHerwig->SetPTRMS(fPtRMS);
	145	fHerwig->SetMAXPR(fMaxPr);
	146	fHerwig->SetMAXER(fMaxErrors);
	147	fHerwig->SetENSOF(fEnSoft);
7677b281	148
e2054d85	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
7677b281	157	// RMASS(5)=4.75
e2054d85	158	// RMASS(6)=174.3
	159	// RMASS(13)=0.75
	160
	161	fHerwig->SetRMASS(4,1.2);
	162	fHerwig->SetRMASS(5,4.75);
7677b281	163
e2054d85	164	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
e2054d85	165
9305fc24	166	//fHerwig->Hwusta("PI0 ");
e2054d85	167
618e1dc0	168	// compute parameter dependent constants
	169	fHerwig->PrepareRun();
	170	}
	171
7677b281	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
9305fc24	206	// fHerwig->Hwusta("PI0 ");
7677b281	207
	208	// compute parameter dependent constants
	209	fHerwig->PrepareRunJimmy();
	210	}
	211
618e1dc0	212	void AliGenHerwig::InitPDF()
	213	{
	214	switch(fStrucFunc)
	215	{
7677b281	216	// ONLY USES LHAPDF STRUCTURE FUNCTIONS
618e1dc0	217	case kGRVLO98:
7677b281	218	fModPDF=80060;
7677b281	219	fAutPDF="HWLHAPDF";
618e1dc0	220	break;
7677b281	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";
618e1dc0	244	break;
618e1dc0	245	case kCTEQ5L:
7677b281	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;
56174e25	257	// case kMRST2004nlo:
	258	// fModPDF=20400;
	259	// fAutPDF="HWLHAPDF";
	260	// break;
618e1dc0	261	default:
	262	cerr << "This structure function is not inplemented " << fStrucFunc << endl;
	263	break;
	264	}
7677b281	265	fAutPDF.Resize(20);
618e1dc0	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];
7677b281	280
618e1dc0	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;
7677b281	289
618e1dc0	290	if(!particles) particles=new TClonesArray("TParticle",10000);
7677b281	291
618e1dc0	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	}
7677b281	301
618e1dc0	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;
7677b281	309
9305fc24	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
618e1dc0	318	Int_t nc=0;
7677b281	319
618e1dc0	320	Int_t * newPos = new Int_t[np];
618e1dc0	321	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
7677b281	322
618e1dc0	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();
7677b281	328	if (ks != 3 &&
7677b281	329	KinematicSelection(iparticle,0))
618e1dc0	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;
93831bcf	342	PushTrack(trackIt, iparent, kf,
93831bcf	343	p[0], p[1], p[2], p[3],
7677b281	344	origin[0], origin[1], origin[2],
93831bcf	345	tof,
93831bcf	346	polar[0], polar[1], polar[2],
e2054d85	347	kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0	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;
618e1dc0	353	// MakeHeader();
618e1dc0	354	if (nc > 0) {
	355	jev+=nc;
	356	if (jev >= fNpart \|\| fNpart == -1) {
	357	fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0	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();
9305fc24	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 ;
618e1dc0	414	}
	415
	416	void AliGenHerwig::AdjustWeights()
	417	{
	418	// Adjust the weights after generation of all events
	419	TParticle *part;
5d12ce38	420	Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0	421	for (Int_t i=0; i<ntrack; i++) {
5d12ce38	422	part= gAlice->GetMCApp()->Particle(i);
618e1dc0	423	part->SetWeight(part->GetWeight()*fKineBias);
	424	}
	425	}
7677b281	426
618e1dc0	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();
7677b281	446	imax=iparticle->GetLastDaughter();
618e1dc0	447	for (i=imin; i<= imax; i++){
7677b281	448	TParticle * jparticle = (TParticle *) particles->At(i);
618e1dc0	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;
7677b281	470
618e1dc0	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());
7677b281	481
618e1dc0	482	if ( (particle->GetFirstDaughter() < 0 ) \|\| (kf == 1000*fFlavor+122))
7677b281	483
618e1dc0	484	{
	485	return kTRUE;
	486	} else {
	487	return kFALSE;
	488	}
	489	}
	490
	491	void AliGenHerwig::FinishRun()
	492	{
	493	fHerwig->Hwefin();
	494	}
	495
7677b281	496	void AliGenHerwig::FinishRunJimmy()
	497	{
	498	fHerwig->Hwefin();
	499	fHerwig->Jmefin();
	500
	501	}
618e1dc0	502