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