[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */


// Generator using Herwig as an external generator
// The main Herwig options are accessable for the user through this interface.
// Uses the THerwig implementation of TGenerator.


#include <Riostream.h>
#include <TClonesArray.h>
#include <TParticle.h>

#include <THerwig6.h>

#include "AliGenHerwig.h"
#include "AliHerwigRndm.h"
#include "AliMC.h"
#include "AliRun.h"
#include "driver.h"

ClassImp(AliGenHerwig)


  AliGenHerwig::AliGenHerwig() :
    AliGenMC(),
    fAutPDF("LHAPDF"),
    fModPDF(19070),
    fStrucFunc(kCTEQ5L),
    fKeep(0),
    fDecaysOff(1),
    fTrigger(0),
    fSelectAll(0),
    fFlavor(0),
    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()
{
// 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();
}

void AliGenHerwig::FinishRunJimmy()
{
  fHerwig->Hwefin();
  fHerwig->Jmefin();

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