[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$
*/


// 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 <stream.h>
#include "AliGenHerwig.h"
#include "AliRun.h"

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


 ClassImp(AliGenHerwig)

AliGenHerwig::AliGenHerwig()
{
// Constructor
}

AliGenHerwig::AliGenHerwig(Int_t 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;
    }
  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[3], 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)=i;
      
      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 == 1 && kf != 0 && KinematicSelection(iparticle,0)) {
	  nc++;
	  p[0]=iparticle->Px();
	  p[1]=iparticle->Py();
	  p[2]=iparticle->Pz();
	  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();
	  SetTrack(fTrackIt,-1,kf,p,origin,polar,tof,kPPrimary,nt);
	  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();

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