[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()
{
// 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
	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	}
618e1dc0	100
618e1dc0	101	AliGenHerwig::~AliGenHerwig()
	102	{
	103	// Destructor
	104	}
	105
e2054d85	106	void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
	107	{
	108	fEv1Pr = ++eventFirst;
	109	fEv2Pr = ++eventLast;
	110	if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
	111	}
	112
618e1dc0	113	void AliGenHerwig::Init()
	114	{
	115	// Initialisation
	116	fTarget.Resize(8);
	117	fProjectile.Resize(8);
	118	SetMC(new THerwig6());
d274cd60	119	fHerwig=(THerwig6*) fMCEvGen;
618e1dc0	120	// initialize common blocks
	121	fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
	122	// reset parameters according to user needs
	123	InitPDF();
	124	fHerwig->SetPTMIN(fPtHardMin);
	125	fHerwig->SetPTRMS(fPtRMS);
	126	fHerwig->SetMAXPR(fMaxPr);
	127	fHerwig->SetMAXER(fMaxErrors);
	128	fHerwig->SetENSOF(fEnSoft);
7677b281	129
e2054d85	130	fHerwig->SetEV1PR(fEv1Pr);
	131	fHerwig->SetEV2PR(fEv2Pr);
	132
	133	// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
	134	// RMASS(1)=0.32
	135	// RMASS(2)=0.32
	136	// RMASS(3)=0.5
	137	// RMASS(4)=1.55
7677b281	138	// RMASS(5)=4.75
e2054d85	139	// RMASS(6)=174.3
	140	// RMASS(13)=0.75
	141
	142	fHerwig->SetRMASS(4,1.2);
	143	fHerwig->SetRMASS(5,4.75);
7677b281	144
e2054d85	145	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
	146
	147	fHerwig->Hwusta("PI0 ");
	148
618e1dc0	149	// compute parameter dependent constants
	150	fHerwig->PrepareRun();
	151	}
	152
7677b281	153	void AliGenHerwig::InitJimmy()
	154	{
	155	// Initialisation
	156	fTarget.Resize(8);
	157	fProjectile.Resize(8);
	158	SetMC(new THerwig6());
	159	fHerwig=(THerwig6*) fMCEvGen;
	160	// initialize common blocks
	161	fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
	162	// reset parameters according to user needs
	163	InitPDF();
	164	fHerwig->SetPTMIN(fPtHardMin);
	165	fHerwig->SetPTRMS(fPtRMS);
	166	fHerwig->SetMAXPR(fMaxPr);
	167	fHerwig->SetMAXER(fMaxErrors);
	168	fHerwig->SetENSOF(fEnSoft);
	169
	170	fHerwig->SetEV1PR(fEv1Pr);
	171	fHerwig->SetEV2PR(fEv2Pr);
	172
	173	// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
	174	// RMASS(1)=0.32
	175	// RMASS(2)=0.32
	176	// RMASS(3)=0.5
	177	// RMASS(4)=1.55
	178	// RMASS(5)=4.75
	179	// RMASS(6)=174.3
	180	// RMASS(13)=0.75
	181
	182	fHerwig->SetRMASS(4,1.2);
	183	fHerwig->SetRMASS(5,4.75);
	184
	185	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
	186
	187	fHerwig->Hwusta("PI0 ");
	188
	189	// compute parameter dependent constants
	190	fHerwig->PrepareRunJimmy();
	191	}
	192
618e1dc0	193	void AliGenHerwig::InitPDF()
	194	{
	195	switch(fStrucFunc)
	196	{
7677b281	197	// ONLY USES LHAPDF STRUCTURE FUNCTIONS
618e1dc0	198	case kGRVLO98:
7677b281	199	fModPDF=80060;
7677b281	200	fAutPDF="HWLHAPDF";
618e1dc0	201	break;
7677b281	202	case kCTEQ6:
	203	fModPDF=10040;
	204	fAutPDF="HWLHAPDF";
	205	break;
	206	case kCTEQ61:
	207	fModPDF=10100;
	208	fAutPDF="HWLHAPDF";
	209	break;
	210	case kCTEQ6m:
	211	fModPDF=10050;
	212	fAutPDF="HWLHAPDF";
	213	break;
	214	case kCTEQ6l:
	215	fModPDF=10041;
	216	fAutPDF="HWLHAPDF";
	217	break;
	218	case kCTEQ6ll:
	219	fModPDF=10042;
	220	fAutPDF="HWLHAPDF";
	221	break;
	222	case kCTEQ5M:
	223	fModPDF=19050;
	224	fAutPDF="HWLHAPDF";
618e1dc0	225	break;
618e1dc0	226	case kCTEQ5L:
7677b281	227	fModPDF=19070;
	228	fAutPDF="HWLHAPDF";
	229	break;
	230	case kCTEQ4M:
	231	fModPDF=19150;
	232	fAutPDF="HWLHAPDF";
	233	break;
	234	case kCTEQ4L:
	235	fModPDF=19170;
	236	fAutPDF="HWLHAPDF";
	237	break;
56174e25	238	// case kMRST2004nlo:
	239	// fModPDF=20400;
	240	// fAutPDF="HWLHAPDF";
	241	// break;
618e1dc0	242	default:
	243	cerr << "This structure function is not inplemented " << fStrucFunc << endl;
	244	break;
	245	}
7677b281	246	fAutPDF.Resize(20);
618e1dc0	247	fHerwig->SetMODPDF(1,fModPDF);
	248	fHerwig->SetMODPDF(2,fModPDF);
	249	fHerwig->SetAUTPDF(1,fAutPDF);
	250	fHerwig->SetAUTPDF(2,fAutPDF);
	251	}
	252
	253	void AliGenHerwig::Generate()
	254	{
	255	// Generate one event
	256
	257	Float_t polar[3] = {0,0,0};
	258	Float_t origin[3]= {0,0,0};
	259	Float_t origin0[3]= {0,0,0};
	260	Float_t p[4], random[6];
7677b281	261
618e1dc0	262	static TClonesArray *particles;
	263	// converts from mm/c to s
	264	const Float_t kconv=0.001/2.999792458e8;
	265	//
	266	Int_t nt=0;
	267	Int_t jev=0;
	268	Int_t j, kf, ks, imo;
	269	kf=0;
7677b281	270
618e1dc0	271	if(!particles) particles=new TClonesArray("TParticle",10000);
7677b281	272
618e1dc0	273	fTrials=0;
	274	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
	275	if(fVertexSmear==kPerEvent) {
	276	Rndm(random,6);
	277	for (j=0;j<3;j++) {
	278	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	279	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	280	}
	281	}
7677b281	282
618e1dc0	283	while(1)
	284	{
	285	fHerwig->GenerateEvent();
	286	fTrials++;
	287	fHerwig->ImportParticles(particles,"All");
	288	Int_t np = particles->GetEntriesFast()-1;
	289	if (np == 0 ) continue;
7677b281	290
618e1dc0	291	Int_t nc=0;
7677b281	292
618e1dc0	293	Int_t * newPos = new Int_t[np];
618e1dc0	294	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
7677b281	295
618e1dc0	296	for (Int_t i = 0; i<np; i++) {
	297	TParticle * iparticle = (TParticle *) particles->At(i);
	298	imo = iparticle->GetFirstMother();
	299	kf = iparticle->GetPdgCode();
	300	ks = iparticle->GetStatusCode();
7677b281	301	if (ks != 3 &&
7677b281	302	KinematicSelection(iparticle,0))
618e1dc0	303	{
	304	nc++;
	305	p[0]=iparticle->Px();
	306	p[1]=iparticle->Py();
	307	p[2]=iparticle->Pz();
	308	p[3]=iparticle->Energy();
	309	origin[0]=origin0[0]+iparticle->Vx()/10;
	310	origin[1]=origin0[1]+iparticle->Vy()/10;
	311	origin[2]=origin0[2]+iparticle->Vz()/10;
	312	Float_t tof = kconv*iparticle->T();
	313	Int_t iparent = (imo > -1) ? newPos[imo] : -1;
	314	Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf	315	PushTrack(trackIt, iparent, kf,
93831bcf	316	p[0], p[1], p[2], p[3],
7677b281	317	origin[0], origin[1], origin[2],
93831bcf	318	tof,
93831bcf	319	polar[0], polar[1], polar[2],
e2054d85	320	kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0	321	KeepTrack(nt);
	322	newPos[i]=nt;
	323	} // end of if: selection of particle
	324	} // end of for: particle loop
	325	if (newPos) delete[] newPos;
618e1dc0	326	// MakeHeader();
618e1dc0	327	if (nc > 0) {
	328	jev+=nc;
	329	if (jev >= fNpart \|\| fNpart == -1) {
	330	fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0	331	break;
	332	}
	333	}
	334	}
	335	SetHighWaterMark(nt);
	336	// adjust weight due to kinematic selection
	337	AdjustWeights();
	338	// get cross-section
	339	fXsection=fHerwig->GetAVWGT();
	340	}
	341
	342	void AliGenHerwig::AdjustWeights()
	343	{
	344	// Adjust the weights after generation of all events
	345	TParticle *part;
5d12ce38	346	Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0	347	for (Int_t i=0; i<ntrack; i++) {
5d12ce38	348	part= gAlice->GetMCApp()->Particle(i);
618e1dc0	349	part->SetWeight(part->GetWeight()*fKineBias);
	350	}
	351	}
7677b281	352
618e1dc0	353
	354	void AliGenHerwig::KeepFullEvent()
	355	{
	356	fKeep=1;
	357	}
	358
	359	Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
	360	{
	361	//
	362	// Looks recursively if one of the daughters has been selected
	363	//
	364	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	365	Int_t imin=-1;
	366	Int_t imax=-1;
	367	Int_t i;
	368	Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
	369	Bool_t selected=kFALSE;
	370	if (hasDaughters) {
	371	imin=iparticle->GetFirstDaughter();
7677b281	372	imax=iparticle->GetLastDaughter();
618e1dc0	373	for (i=imin; i<= imax; i++){
7677b281	374	TParticle * jparticle = (TParticle *) particles->At(i);
618e1dc0	375	Int_t ip=jparticle->GetPdgCode();
	376	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	377	selected=kTRUE; break;
	378	}
	379	if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	380	}
	381	} else {
	382	return kFALSE;
	383	}
	384
	385	return selected;
	386	}
	387
	388
	389	Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
	390	{
	391	// Select flavor of particle
	392	// 0: all
	393	// 4: charm and beauty
	394	// 5: beauty
	395	if (fFlavor == 0) return kTRUE;
7677b281	396
618e1dc0	397	Int_t ifl=TMath::Abs(pid/100);
	398	if (ifl > 10) ifl/=10;
	399	return (fFlavor == ifl);
	400	}
	401
	402	Bool_t AliGenHerwig::Stable(TParticle* particle)
	403	{
	404	// Return true for a stable particle
	405	//
	406	Int_t kf = TMath::Abs(particle->GetPdgCode());
7677b281	407
618e1dc0	408	if ( (particle->GetFirstDaughter() < 0 ) \|\| (kf == 1000*fFlavor+122))
7677b281	409
618e1dc0	410	{
	411	return kTRUE;
	412	} else {
	413	return kFALSE;
	414	}
	415	}
	416
	417	void AliGenHerwig::FinishRun()
	418	{
	419	fHerwig->Hwefin();
	420	}
	421
7677b281	422	void AliGenHerwig::FinishRunJimmy()
	423	{
	424	fHerwig->Hwefin();
	425	fHerwig->Jmefin();
	426
	427	}
618e1dc0	428