[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"

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()),
    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),
    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()),
    fHeader(0)
{
// Constructor
    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 = fEv1Pr;
}

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

	// 
	if (gAlice) {
	    if (gAlice->GetEvNumber()>=fEv1Pr &&
		gAlice->GetEvNumber()<=fEv2Pr) fHerwig->PrintEvt();

	}

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