[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
1c44e4d8	1
618e1dc0	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
7cdba479	17	/* $Id$ */
618e1dc0	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
618e1dc0	24
37b09b91	25	#include <Riostream.h>
37b09b91	26	#include <TClonesArray.h>
618e1dc0	27	#include <TParticle.h>
618e1dc0	28
37b09b91	29	#include <THerwig6.h>
618e1dc0	30
37b09b91	31	#include "AliGenHerwig.h"
aa9da500	32	#include "AliGenHerwigEventHeader.h"
37b09b91	33	#include "AliHerwigRndm.h"
	34	#include "AliMC.h"
	35	#include "AliRun.h"
e2054d85	36	#include "driver.h"
e2054d85	37
c64cb1f6	38	using std::cerr;
	39	using std::endl;
	40
014a9521	41	ClassImp(AliGenHerwig)
618e1dc0	42
7cdba479	43
a1d47d99	44	AliGenHerwig::AliGenHerwig() :
a1d47d99	45	AliGenMC(),
7677b281	46	fAutPDF("LHAPDF"),
7677b281	47	fModPDF(19070),
2ad1e1c2	48	fStrucFunc(kCTEQ5L),
a1d47d99	49	fKeep(0),
	50	fDecaysOff(1),
	51	fTrigger(0),
	52	fSelectAll(0),
	53	fFlavor(0),
a1d47d99	54	fMomentum1(7000),
	55	fMomentum2(7000),
	56	fKineBias(1),
	57	fTrials(0),
	58	fXsection(0),
	59	fHerwig(0x0),
	60	fProcess(0),
2d654757	61	fPtHardMin(0.),
9305fc24	62	fPtHardMax(9999.),
2d654757	63	fPtRMS(0.),
a1d47d99	64	fMaxPr(10),
a1d47d99	65	fMaxErrors(1000),
e2054d85	66	fEnSoft(1),
9305fc24	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.),
aa9da500	75	fPhiMaxGamma(2. * TMath::Pi()),
aa9da500	76	fHeader(0)
618e1dc0	77	{
618e1dc0	78	// Constructor
fc7e1b1c	79	fEnergyCMS = 14000;
618e1dc0	80	}
	81
	82	AliGenHerwig::AliGenHerwig(Int_t npart)
2d654757	83	:AliGenMC(npart),
7677b281	84	fAutPDF("LHAPDF"),
7677b281	85	fModPDF(19070),
2d654757	86	fStrucFunc(kCTEQ5L),
	87	fKeep(0),
	88	fDecaysOff(1),
	89	fTrigger(0),
	90	fSelectAll(0),
	91	fFlavor(0),
2d654757	92	fMomentum1(7000),
	93	fMomentum2(7000),
	94	fKineBias(1),
	95	fTrials(0),
	96	fXsection(0),
	97	fHerwig(0x0),
	98	fProcess(0),
	99	fPtHardMin(0.),
9305fc24	100	fPtHardMax(9999.),
2d654757	101	fPtRMS(0.),
	102	fMaxPr(10),
	103	fMaxErrors(1000),
	104	fEnSoft(1),
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
	129	void AliGenHerwig::Init()
	130	{
	131	// Initialisation
	132	fTarget.Resize(8);
	133	fProjectile.Resize(8);
	134	SetMC(new THerwig6());
d274cd60	135	fHerwig=(THerwig6*) fMCEvGen;
618e1dc0	136	// initialize common blocks
455f128d	137	fHerwig->Initialize(fProjectile.Data(), fTarget.Data(), fMomentum1, fMomentum2, fProcess); //
618e1dc0	138	// reset parameters according to user needs
	139	InitPDF();
	140	fHerwig->SetPTMIN(fPtHardMin);
7a5495ad	141	fHerwig->SetPTMAX(fPtHardMax);
618e1dc0	142	fHerwig->SetPTRMS(fPtRMS);
1c44e4d8	143	printf("SetMAXPR %15d \n", fMaxPr);
618e1dc0	144	fHerwig->SetMAXPR(fMaxPr);
	145	fHerwig->SetMAXER(fMaxErrors);
	146	fHerwig->SetENSOF(fEnSoft);
e2054d85	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
7677b281	152	// RMASS(5)=4.75
e2054d85	153	// RMASS(6)=174.3
	154	// RMASS(13)=0.75
	155
	156	fHerwig->SetRMASS(4,1.2);
	157	fHerwig->SetRMASS(5,4.75);
7677b281	158
09bba509	159	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
e2054d85	160
9305fc24	161	//fHerwig->Hwusta("PI0 ");
e2054d85	162
618e1dc0	163	// compute parameter dependent constants
	164	fHerwig->PrepareRun();
	165	}
	166
7677b281	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
7677b281	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
09bba509	196	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
7677b281	197
9305fc24	198	// fHerwig->Hwusta("PI0 ");
7677b281	199
	200	// compute parameter dependent constants
	201	fHerwig->PrepareRunJimmy();
	202	}
	203
618e1dc0	204	void AliGenHerwig::InitPDF()
618e1dc0	205	{
552f205f	206	// Initialize PDF
618e1dc0	207	switch(fStrucFunc)
618e1dc0	208	{
7677b281	209	// ONLY USES LHAPDF STRUCTURE FUNCTIONS
618e1dc0	210	case kGRVLO98:
7677b281	211	fModPDF=80060;
7677b281	212	fAutPDF="HWLHAPDF";
618e1dc0	213	break;
7677b281	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";
618e1dc0	237	break;
618e1dc0	238	case kCTEQ5L:
7677b281	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;
56174e25	250	// case kMRST2004nlo:
	251	// fModPDF=20400;
	252	// fAutPDF="HWLHAPDF";
	253	// break;
618e1dc0	254	default:
	255	cerr << "This structure function is not inplemented " << fStrucFunc << endl;
	256	break;
	257	}
7677b281	258	fAutPDF.Resize(20);
618e1dc0	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
aa9da500	269	Float_t polar[3] = {0,0,0};
	270	Float_t origin[3] = {0,0,0};
	271	Float_t p[4];
7677b281	272
618e1dc0	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;
aa9da500	279	Int_t kf, ks, imo;
618e1dc0	280	kf=0;
7677b281	281
618e1dc0	282	if(!particles) particles=new TClonesArray("TParticle",10000);
7677b281	283
618e1dc0	284	fTrials=0;
aa9da500	285
	286	// Set collision vertex position
	287	if (fVertexSmear == kPerEvent) Vertex();
7677b281	288
618e1dc0	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;
7677b281	296
9305fc24	297	//Check hard partons or direct gamma in kine range
	298
	299	if (fProcess == kHeJets \|\| fProcess == kHeDirectGamma) {
aa9da500	300	TParticle* parton1 = (TParticle *) particles->At(6);
	301	TParticle* parton2 = (TParticle *) particles->At(7);
	302	if (!CheckParton(parton1, parton2)) continue ;
9305fc24	303	}
9305fc24	304
aa9da500	305	Int_t nc = 0;
	306	fNprimaries = 0;
	307
618e1dc0	308	Int_t * newPos = new Int_t[np];
618e1dc0	309	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
7677b281	310
618e1dc0	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();
7677b281	316	if (ks != 3 &&
7677b281	317	KinematicSelection(iparticle,0))
618e1dc0	318	{
	319	nc++;
	320	p[0]=iparticle->Px();
	321	p[1]=iparticle->Py();
	322	p[2]=iparticle->Pz();
	323	p[3]=iparticle->Energy();
aa9da500	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
21391258	329	Float_t tof = fTime + kconv*iparticle->T();
618e1dc0	330	Int_t iparent = (imo > -1) ? newPos[imo] : -1;
618e1dc0	331	Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf	332	PushTrack(trackIt, iparent, kf,
93831bcf	333	p[0], p[1], p[2], p[3],
7677b281	334	origin[0], origin[1], origin[2],
93831bcf	335	tof,
93831bcf	336	polar[0], polar[1], polar[2],
e2054d85	337	kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0	338	KeepTrack(nt);
618e1dc0	339	newPos[i]=nt;
aa9da500	340	fNprimaries++;
618e1dc0	341	} // end of if: selection of particle
	342	} // end of for: particle loop
	343	if (newPos) delete[] newPos;
618e1dc0	344	if (nc > 0) {
	345	jev+=nc;
	346	if (jev >= fNpart \|\| fNpart == -1) {
	347	fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0	348	break;
	349	}
	350	}
	351	}
aa9da500	352	//
	353	MakeHeader();
	354	//
618e1dc0	355	SetHighWaterMark(nt);
	356	// adjust weight due to kinematic selection
	357	AdjustWeights();
	358	// get cross-section
	359	fXsection=fHerwig->GetAVWGT();
9305fc24	360	//printf(">> trials << %d\n",fTrials);
	361	}
	362
552f205f	363	Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const TParticle* parton2)
9305fc24	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 ;
618e1dc0	407	}
	408
	409	void AliGenHerwig::AdjustWeights()
	410	{
	411	// Adjust the weights after generation of all events
	412	TParticle *part;
5d12ce38	413	Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0	414	for (Int_t i=0; i<ntrack; i++) {
5d12ce38	415	part= gAlice->GetMCApp()->Particle(i);
618e1dc0	416	part->SetWeight(part->GetWeight()*fKineBias);
	417	}
	418	}
7677b281	419
618e1dc0	420
	421	void AliGenHerwig::KeepFullEvent()
	422	{
	423	fKeep=1;
	424	}
	425
552f205f	426	Bool_t AliGenHerwig::DaughtersSelection(const TParticle* iparticle, const TClonesArray* particles)
618e1dc0	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();
7677b281	439	imax=iparticle->GetLastDaughter();
618e1dc0	440	for (i=imin; i<= imax; i++){
7677b281	441	TParticle * jparticle = (TParticle *) particles->At(i);
618e1dc0	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
552f205f	456	Bool_t AliGenHerwig::SelectFlavor(Int_t pid) const
618e1dc0	457	{
	458	// Select flavor of particle
	459	// 0: all
	460	// 4: charm and beauty
	461	// 5: beauty
	462	if (fFlavor == 0) return kTRUE;
7677b281	463
618e1dc0	464	Int_t ifl=TMath::Abs(pid/100);
	465	if (ifl > 10) ifl/=10;
	466	return (fFlavor == ifl);
	467	}
	468
552f205f	469	Bool_t AliGenHerwig::Stable(const TParticle* particle) const
618e1dc0	470	{
	471	// Return true for a stable particle
	472	//
	473	Int_t kf = TMath::Abs(particle->GetPdgCode());
7677b281	474
618e1dc0	475	if ( (particle->GetFirstDaughter() < 0 ) \|\| (kf == 1000*fFlavor+122))
7677b281	476
618e1dc0	477	{
	478	return kTRUE;
	479	} else {
	480	return kFALSE;
	481	}
	482	}
	483
	484	void AliGenHerwig::FinishRun()
	485	{
	486	fHerwig->Hwefin();
	487	}
	488
7677b281	489	void AliGenHerwig::FinishRunJimmy()
	490	{
	491	fHerwig->Hwefin();
	492	fHerwig->Jmefin();
	493
	494	}
618e1dc0	495
aa9da500	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	//
616d81ca	511	// Event weight (cross section)
	512	((AliGenHerwigEventHeader*) fHeader)->SetWeight(fHerwig->GetEVWGT());
	513	//
aa9da500	514	// Event Vertex
aa9da500	515	fHeader->SetPrimaryVertex(fVertex);
21391258	516	fHeader->SetInteractionTime(fTime);
aa9da500	517	//
	518	// Number of primaries
	519	fHeader->SetNProduced(fNprimaries);
	520	// Pass header
	521	//
	522	AddHeader(fHeader);
	523	fHeader = 0x0;
	524	}