[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 "AliGenHerwig.h"
#include "AliHerwigRndm.h"
#include "AliRun.h"

#include <TParticle.h>
#include "THerwig6.h"

#include "Riostream.h"
#include "AliMC.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),
    fEnergyCMS(14000),
    fMomentum1(7000),
    fMomentum2(7000),
    fKineBias(1),
    fTrials(0),
    fXsection(0),
    fHerwig(0x0),
    fProcess(0),
    fPtHardMin(0.),
    fPtRMS(0.),
    fMaxPr(10),
    fMaxErrors(1000),
    fEnSoft(1),
    fEv1Pr(0),
    fEv2Pr(0),
    fFileName(0)
{
// Constructor
}

AliGenHerwig::AliGenHerwig(Int_t npart)
    :AliGenMC(npart),
    fAutPDF("LHAPDF"),
    fModPDF(19070),
    fStrucFunc(kCTEQ5L),
    fKeep(0),
    fDecaysOff(1),
    fTrigger(0),
    fSelectAll(0),
    fFlavor(0),
    fEnergyCMS(14000),
    fMomentum1(7000),
    fMomentum2(7000),
    fKineBias(1),
    fTrials(0),
    fXsection(0),
    fHerwig(0x0),
    fProcess(0),
    fPtHardMin(0.),
    fPtRMS(0.),
    fMaxPr(10),
    fMaxErrors(1000),
    fEnSoft(1),
    fEv1Pr(0),
    fEv2Pr(0),
    fFileName(0)
{
    SetTarget();
    SetProjectile();
    // Set random number generator
    AliHerwigRndm::SetHerwigRandom(GetRandom());
}

AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
    :AliGenMC(Herwig),
    fAutPDF("LHAPDF"),
    fModPDF(19070),
    fStrucFunc(kCTEQ5L),
    fKeep(0),
    fDecaysOff(1),
    fTrigger(0),
    fSelectAll(0),
    fFlavor(0),
    fEnergyCMS(14000),
    fMomentum1(7000),
    fMomentum2(7000),
    fKineBias(1),
    fTrials(0),
    fXsection(0),
    fHerwig(0x0),
    fProcess(0),
    fPtHardMin(0.),
    fPtRMS(0.),
    fMaxPr(10),
    fMaxErrors(1000),
    fEnSoft(1),
    fEv1Pr(0),
    fEv2Pr(0),
    fFileName(0)
{
// Copy constructor
    Herwig.Copy(*this);
}


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;

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

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


AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return (*this);
}

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