[u/mrichter/AliRoot.git] / EVGEN / 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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.2.2.3  2002/10/11 10:40:46  hristov
Default case added

Revision 1.2.2.2  2002/07/26 18:34:02  alibrary
Updating VirtualMC

Revision 1.3  2002/07/26 15:32:24  hristov
stream.h doesn't exest on Sun, removed from includes

Revision 1.2  2002/07/19 11:43:10  morsch
- Write full stack.
- Use SetTrack passing energy.

Revision 1.1  2002/07/16 11:33:26  morsch
First commit.

*/


// 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 "AliRun.h"

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

#include "Riostream.h"

 ClassImp(AliGenHerwig)

AliGenHerwig::AliGenHerwig()
{
// Constructor
}

AliGenHerwig::AliGenHerwig(Int_t npart)
    :AliGenMC(npart)
{
    SetBeamMomenta();
    SetTarget();
    SetProjectile();
    SetStrucFunc(kGRVLO98);
    fKeep=0;
    fTrigger=0;
    fDecaysOff=1;
    fSelectAll=0;
    fFlavor=0;
    fPtHardMin=10.;
    fPtRMS=0.0;
    fEnSoft=1.0;
    fMaxPr=1;
    fMaxErrors=1000;
//  Set random number
    if (!sRandom) sRandom=fRandom;
}

AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
{
// copy constructor
}


AliGenHerwig::~AliGenHerwig()
{
// Destructor
}

void AliGenHerwig::Init()
{
// Initialisation
  fTarget.Resize(8);
  fProjectile.Resize(8);
  SetMC(new THerwig6());
  fHerwig=(THerwig6*) fgMCEvGen;
  // 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);
  // compute parameter dependent constants
  fHerwig->PrepareRun();
}

void AliGenHerwig::InitPDF()
{
  switch(fStrucFunc)
    {
    case kGRVLO:
      fModPDF=5;
      fAutPDF="GRV";
      break;
    case kGRVHO:
      fModPDF=6;
      fAutPDF="GRV";
      break;
    case kGRVLO98:
      fModPDF=12;
      fAutPDF="GRV";
      break;
    case kMRSDminus:
      fModPDF=31;
      fAutPDF="MRS";
      break;
    case kMRSD0:
      fModPDF=30;
      fAutPDF="MRS";
      break;
    case kMRSG:
      fModPDF=41;
      fAutPDF="MRS";
      break;
    case kMRSTcgLO:
      fModPDF=72;
      fAutPDF="MRS";
      break;
    case kCTEQ4M:
      fModPDF=34;
      fAutPDF="CTEQ";
      break;
    case kCTEQ5L:
      fModPDF=46;
      fAutPDF="CTEQ";
      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;
		gAlice->SetTrack(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, 1., ks);
		KeepTrack(nt);
		newPos[i]=nt;
	    } // end of if: selection of particle
	} // end of for: particle loop
	if (newPos) delete[] newPos;
	printf("\n I've put %i particles on the stack \n",nc);
	//      MakeHeader();
	printf("nc: %d %d\n", nc, fNpart);
	
	if (nc > 0) {
	    jev+=nc;
	    if (jev >= fNpart || fNpart == -1) {
		fKineBias=Float_t(fNpart)/Float_t(fTrials);
		printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
		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->GetNtrack();
    for (Int_t i=0; i<ntrack; i++) {
        part= gAlice->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
    return *this;
}


extern "C" {
  Double_t hwr_() {return sRandom->Rndm();}
}
Commit	Line	Data
b46943fe	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
	16	/*
	17	$Log$
b9d0a01d	18	Revision 1.2.2.3 2002/10/11 10:40:46 hristov
	19	Default case added
	20
	21	Revision 1.2.2.2 2002/07/26 18:34:02 alibrary
	22	Updating VirtualMC
	23
	24	Revision 1.3 2002/07/26 15:32:24 hristov
	25	stream.h doesn't exest on Sun, removed from includes
	26
1784616d	27	Revision 1.2 2002/07/19 11:43:10 morsch
	28	- Write full stack.
	29	- Use SetTrack passing energy.
	30
dc2691e2	31	Revision 1.1 2002/07/16 11:33:26 morsch
	32	First commit.
	33
b46943fe	34	*/
	35
	36
	37
	38	// Generator using Herwig as an external generator
	39	// The main Herwig options are accessable for the user through this interface.
	40	// Uses the THerwig implementation of TGenerator.
	41
b46943fe	42	#include "AliGenHerwig.h"
	43	#include "AliRun.h"
	44
	45	#include <TParticle.h>
	46	#include "THerwig6.h"
	47
b9d0a01d	48	#include "Riostream.h"
b46943fe	49
	50	ClassImp(AliGenHerwig)
	51
	52	AliGenHerwig::AliGenHerwig()
	53	{
	54	// Constructor
	55	}
	56
	57	AliGenHerwig::AliGenHerwig(Int_t npart)
dc2691e2	58	:AliGenMC(npart)
b46943fe	59	{
	60	SetBeamMomenta();
	61	SetTarget();
	62	SetProjectile();
	63	SetStrucFunc(kGRVLO98);
	64	fKeep=0;
	65	fTrigger=0;
	66	fDecaysOff=1;
	67	fSelectAll=0;
	68	fFlavor=0;
	69	fPtHardMin=10.;
	70	fPtRMS=0.0;
	71	fEnSoft=1.0;
	72	fMaxPr=1;
	73	fMaxErrors=1000;
	74	// Set random number
	75	if (!sRandom) sRandom=fRandom;
	76	}
	77
	78	AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
	79	{
	80	// copy constructor
	81	}
	82
	83
	84	AliGenHerwig::~AliGenHerwig()
	85	{
	86	// Destructor
	87	}
	88
	89	void AliGenHerwig::Init()
	90	{
	91	// Initialisation
	92	fTarget.Resize(8);
	93	fProjectile.Resize(8);
	94	SetMC(new THerwig6());
	95	fHerwig=(THerwig6*) fgMCEvGen;
	96	// initialize common blocks
	97	fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
	98	// reset parameters according to user needs
	99	InitPDF();
	100	fHerwig->SetPTMIN(fPtHardMin);
	101	fHerwig->SetPTRMS(fPtRMS);
	102	fHerwig->SetMAXPR(fMaxPr);
	103	fHerwig->SetMAXER(fMaxErrors);
	104	fHerwig->SetENSOF(fEnSoft);
	105	// compute parameter dependent constants
	106	fHerwig->PrepareRun();
	107	}
	108
	109	void AliGenHerwig::InitPDF()
	110	{
	111	switch(fStrucFunc)
	112	{
	113	case kGRVLO:
	114	fModPDF=5;
	115	fAutPDF="GRV";
	116	break;
	117	case kGRVHO:
	118	fModPDF=6;
	119	fAutPDF="GRV";
	120	break;
	121	case kGRVLO98:
	122	fModPDF=12;
123	fAutPDF="GRV";
124	break;
125	case kMRSDminus:
126	fModPDF=31;
127	fAutPDF="MRS";
128	break;
129	case kMRSD0:
130	fModPDF=30;
131	fAutPDF="MRS";
132	break;
133	case kMRSG:
134	fModPDF=41;
135	fAutPDF="MRS";
136	break;
137	case kMRSTcgLO:
138	fModPDF=72;
139	fAutPDF="MRS";
140	break;
141	case kCTEQ4M:
142	fModPDF=34;
143	fAutPDF="CTEQ";
144	break;
145	case kCTEQ5L:
146	fModPDF=46;
147	fAutPDF="CTEQ";
148	break;
b9d0a01d	149	default:
	150	cerr << "This structure function is not inplemented " << fStrucFunc << endl;
	151	break;
b46943fe	152	}
	153	fAutPDF.Resize(20);
	154	fHerwig->SetMODPDF(1,fModPDF);
	155	fHerwig->SetMODPDF(2,fModPDF);
	156	fHerwig->SetAUTPDF(1,fAutPDF);
	157	fHerwig->SetAUTPDF(2,fAutPDF);
	158	}
	159
	160	void AliGenHerwig::Generate()
	161	{
	162	// Generate one event
	163
	164	Float_t polar[3] = {0,0,0};
	165	Float_t origin[3]= {0,0,0};
	166	Float_t origin0[3]= {0,0,0};
dc2691e2	167	Float_t p[4], random[6];
b46943fe	168
	169	static TClonesArray *particles;
	170	// converts from mm/c to s
	171	const Float_t kconv=0.001/2.999792458e8;
	172	//
	173	Int_t nt=0;
	174	Int_t jev=0;
	175	Int_t j, kf, ks, imo;
	176	kf=0;
	177
	178	if(!particles) particles=new TClonesArray("TParticle",10000);
	179
	180	fTrials=0;
	181	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
	182	if(fVertexSmear==kPerEvent) {
	183	Rndm(random,6);
	184	for (j=0;j<3;j++) {
	185	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	186	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	187	}
	188	}
dc2691e2	189
b46943fe	190	while(1)
b46943fe	191	{
dc2691e2	192	fHerwig->GenerateEvent();
	193	fTrials++;
	194	fHerwig->ImportParticles(particles,"All");
	195	Int_t np = particles->GetEntriesFast()-1;
	196	if (np == 0 ) continue;
	197
	198	Int_t nc=0;
	199
	200	Int_t * newPos = new Int_t[np];
	201	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
	202
	203	for (Int_t i = 0; i<np; i++) {
	204	TParticle * iparticle = (TParticle *) particles->At(i);
	205	imo = iparticle->GetFirstMother();
	206	kf = iparticle->GetPdgCode();
	207	ks = iparticle->GetStatusCode();
	208	if (ks != 3 &&
	209	KinematicSelection(iparticle,0))
	210	{
	211	nc++;
	212	p[0]=iparticle->Px();
	213	p[1]=iparticle->Py();
	214	p[2]=iparticle->Pz();
	215	p[3]=iparticle->Energy();
	216	origin[0]=origin0[0]+iparticle->Vx()/10;
	217	origin[1]=origin0[1]+iparticle->Vy()/10;
	218	origin[2]=origin0[2]+iparticle->Vz()/10;
	219	Float_t tof = kconv*iparticle->T();
	220	Int_t iparent = (imo > -1) ? newPos[imo] : -1;
	221	Int_t trackIt = (ks == 1) && fTrackIt;
	222	gAlice->SetTrack(trackIt, iparent, kf,
	223	p[0], p[1], p[2], p[3],
	224	origin[0], origin[1], origin[2],
	225	tof,
	226	polar[0], polar[1], polar[2],
	227	kPPrimary, nt, 1., ks);
	228	KeepTrack(nt);
	229	newPos[i]=nt;
	230	} // end of if: selection of particle
	231	} // end of for: particle loop
	232	if (newPos) delete[] newPos;
	233	printf("\n I've put %i particles on the stack \n",nc);
	234	// MakeHeader();
	235	printf("nc: %d %d\n", nc, fNpart);
	236
	237	if (nc > 0) {
	238	jev+=nc;
	239	if (jev >= fNpart \|\| fNpart == -1) {
	240	fKineBias=Float_t(fNpart)/Float_t(fTrials);
	241	printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
	242	break;
	243	}
b46943fe	244	}
b46943fe	245	}
dc2691e2	246	SetHighWaterMark(nt);
b46943fe	247	// adjust weight due to kinematic selection
dc2691e2	248	AdjustWeights();
b46943fe	249	// get cross-section
dc2691e2	250	fXsection=fHerwig->GetAVWGT();
b46943fe	251	}
	252
	253	void AliGenHerwig::AdjustWeights()
	254	{
	255	// Adjust the weights after generation of all events
	256	TParticle *part;
	257	Int_t ntrack=gAlice->GetNtrack();
	258	for (Int_t i=0; i<ntrack; i++) {
	259	part= gAlice->Particle(i);
	260	part->SetWeight(part->GetWeight()*fKineBias);
	261	}
	262	}
	263
	264
	265	void AliGenHerwig::KeepFullEvent()
	266	{
	267	fKeep=1;
	268	}
	269
	270	Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
	271	{
	272	//
	273	// Looks recursively if one of the daughters has been selected
	274	//
	275	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	276	Int_t imin=-1;
	277	Int_t imax=-1;
	278	Int_t i;
	279	Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
	280	Bool_t selected=kFALSE;
	281	if (hasDaughters) {
	282	imin=iparticle->GetFirstDaughter();
	283	imax=iparticle->GetLastDaughter();
	284	for (i=imin; i<= imax; i++){
	285	TParticle * jparticle = (TParticle *) particles->At(i);
	286	Int_t ip=jparticle->GetPdgCode();
	287	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	288	selected=kTRUE; break;
	289	}
	290	if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	291	}
	292	} else {
	293	return kFALSE;
	294	}
	295
	296	return selected;
	297	}
	298
	299
	300	Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
	301	{
	302	// Select flavor of particle
	303	// 0: all
	304	// 4: charm and beauty
	305	// 5: beauty
	306	if (fFlavor == 0) return kTRUE;
	307
	308	Int_t ifl=TMath::Abs(pid/100);
	309	if (ifl > 10) ifl/=10;
	310	return (fFlavor == ifl);
	311	}
	312
	313	Bool_t AliGenHerwig::Stable(TParticle* particle)
	314	{
315	// Return true for a stable particle
316	//
317	Int_t kf = TMath::Abs(particle->GetPdgCode());
318
319	if ( (particle->GetFirstDaughter() < 0 ) \|\| (kf == 1000*fFlavor+122))
320
321	{
322	return kTRUE;
323	} else {
324	return kFALSE;
325	}
326	}
327
328	void AliGenHerwig::FinishRun()
329	{
330	fHerwig->Hwefin();
331	}
332
dc2691e2	333
b46943fe	334	AliGenHerwig& AliGenHerwig::operator=(const AliGenHerwig& rhs)
	335	{
	336	// Assignment operator
	337	return *this;
	338	}
	339
	340
	341	extern "C" {
	342	Double_t hwr_() {return sRandom->Rndm();}
	343	}