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

#include "AliGenHijing.h"
#include "AliRun.h"

#include <TArrayI.h>
#include <TParticle.h>
#include <THijing.h>

 ClassImp(AliGenHijing)

AliGenHijing::AliGenHijing()
                 :AliGenerator()
{
// Constructor
}

AliGenHijing::AliGenHijing(Int_t npart)
                 :AliGenerator(npart)
{
// Default PbPb collisions at 5. 5 TeV
//
    SetEnergyCMS();
    SetImpactParameterRange();
    SetTarget();
    SetProjectile();
    fKeep=0;
    fQuench=1;
    fShadowing=1;
    fTrigger=0;
    fDecaysOff=1;
    fEvaluate=0;
    fSelectAll=0;
}

AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
{
// copy constructor
}


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

void AliGenHijing::Init()
{
// Initialisation
    SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
		      fAProjectile, fZProjectile, fATarget, fZTarget, 
		      fMinImpactParam, fMaxImpactParam));

    fHijing=(THijing*) fgMCEvGen;
    fHijing->Initialize();
    fHijing->SetIHPR2(3,  fTrigger);
    fHijing->SetIHPR2(4,  fQuench);
    fHijing->SetIHPR2(6,  fShadowing);
    fHijing->SetIHPR2(12, fDecaysOff);    
    fHijing->SetIHPR2(21, fKeep);
//
    if (fEvaluate) EvaluateCrossSections();
}

void AliGenHijing::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];
    Float_t tof;

    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;

    if(!particles) particles=new TClonesArray("TParticle",10000);
    
    fTrials=0;
    for (j=0;j<3;j++) origin0[j]=fOrigin[j];
    if(fVertexSmear==perEvent) {
	gMC->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]));
//	    fHijing->SetMSTP(151,0);
	}
    } else if (fVertexSmear==perTrack) {
//	fHijing->SetMSTP(151,0);
	for (j=0;j<3;j++) {
//	    fHijing->SetPARP(151+j, fOsigma[j]*10.);
	}
    }
    while(1)
    {

	fHijing->GenerateEvent();
	fTrials++;
	fHijing->ImportParticles(particles,"All");
	Int_t np = particles->GetEntriesFast();
	printf("\n **************************************************%d\n",np);
	Int_t nc=0;
	if (np == 0 ) continue;
	Int_t i;
	Int_t * newPos = new Int_t[np];
	for (i = 0; i<np-1; i++) *(newPos+i)=i;
	
	for (i = 0; i<np-1; i++) {
	    TParticle *  iparticle       = (TParticle *) particles->At(i);

	    Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
	    Bool_t  hasDaughter          =  (iparticle->GetFirstDaughter() >=0);
	    Bool_t  selected             =  kTRUE;
	    Bool_t  hasSelectedDaughters =  kFALSE;

	    if (!fSelectAll) selected = KinematicSelection(iparticle);
	    kf        = iparticle->GetPdgCode();
//	    Int_t id1=iparticle->GetFirstDaughter();
//	    Int_t id2=iparticle->GetLastDaughter();	    
//	    printf("\n particle  %d %d %d %d %d %d\n",i, kf, id1, id2, hasDaughter, selected);

	    if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
//
// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
//

	    if (selected || hasSelectedDaughters) {
		nc++;
		ks        = iparticle->GetStatusCode();
		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;
		tof=kconv*iparticle->T();
		imo=-1;
		if (hasMother) {
		    imo=iparticle->GetFirstMother();
		    imo=*(newPos+imo);
		}
		
//		printf("\n selected iparent %d %d %d \n",i, kf, imo);
		if (hasDaughter) {
		    gAlice->SetTrack(0,imo,kf,p,origin,polar,
				     tof,"Primary",nt);
		} else {
		    gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
				     tof,"Secondary",nt);
		}
		*(newPos+i)=nt;
	    } // selected
	} // particle loop 
	delete newPos;
	printf("\n I've put %i particles on the stack \n",nc);
	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;
	    }
	}
    } // event loop
}

Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
{
// Perform kinematic selection
    Float_t px=particle->Px();
    Float_t py=particle->Py();
    Float_t pz=particle->Pz();
    Float_t  e=particle->Energy();

//
//  transverse momentum cut    
    Float_t pt=TMath::Sqrt(px*px+py*py);
    if (pt > fPtMax || pt < fPtMin) 
    {
//	printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
	return kFALSE;
    }
//
// momentum cut
    Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
    if (p > fPMax || p < fPMin) 
    {
//	printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
	return kFALSE;
    }
    
//
// theta cut
    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
    if (theta > fThetaMax || theta < fThetaMin) 
    {
	
//    	printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
	return kFALSE;
    }

//
// rapidity cut
    Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
    if (y > fYMax || y < fYMin)
    {
//	printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
	return kFALSE;
    }

//
// phi cut
    Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
    if (phi > fPhiMax || phi < fPhiMin)
    {
//	printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
	return kFALSE;
    }

    return kTRUE;
}

void AliGenHijing::KeepFullEvent()
{
    fKeep=1;
}

void AliGenHijing::EvaluateCrossSections()
{
//     Glauber Calculation of geometrical x-section
//
    Float_t xTot=0.;          // barn
    Float_t xTotHard=0.;      // barn 
    Float_t xPart=0.;         // barn
    Float_t xPartHard=0.;     // barn 
    Float_t sigmaHard=0.1;    // mbarn
    Float_t bMin=0.;
    Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
    const Float_t kdib=0.2;
    Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;


    printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
    printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
    Int_t i;
    Float_t oldvalue=0.;
    
    for (i=0; i<kMax; i++)
    {
	Float_t xb=bMin+i*kdib;
	Float_t ov;
	ov=fHijing->Profile(xb);
	Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
	Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
	xTot+=gb;
	xTotHard+=gbh;
	if (xb > fMinImpactParam && xb < fMaxImpactParam)
	{
	    xPart+=gb;
	    xPartHard+=gbh;
	}
	
	if ((xTot-oldvalue)/oldvalue<0.0001) break;
	oldvalue=xTot;
	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
	printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
    }
    printf("\n Total cross section (barn): %f \n",xTot);
    printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
    printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
    printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
}

Bool_t AliGenHijing::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();       
//	printf("\n Has daughters %d %d:", imin, imax);
	for (i=imin; i<= imax; i++){
	    TParticle *  jparticle       = (TParticle *) particles->At(i);	
//	    Int_t ip=jparticle->GetPdgCode();
//	    printf("\n consider daughter %d %d", i,ip);
	    
	    if (KinematicSelection(jparticle)) {selected=kTRUE; break;}
	    if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	}
    } else {
	return kFALSE;
    }

    return selected;
}


AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
{
// Assignment operator
    return *this;
}
Commit	Line	Data
110410f9	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	#include "AliGenHijing.h"
	21	#include "AliRun.h"
	22
	23	#include <TArrayI.h>
	24	#include <TParticle.h>
	25	#include <THijing.h>
	26
	27	ClassImp(AliGenHijing)
	28
	29	AliGenHijing::AliGenHijing()
	30	:AliGenerator()
	31	{
	32	// Constructor
	33	}
	34
	35	AliGenHijing::AliGenHijing(Int_t npart)
	36	:AliGenerator(npart)
	37	{
	38	// Default PbPb collisions at 5. 5 TeV
	39	//
	40	SetEnergyCMS();
	41	SetImpactParameterRange();
	42	SetTarget();
	43	SetProjectile();
	44	fKeep=0;
	45	fQuench=1;
	46	fShadowing=1;
	47	fTrigger=0;
	48	fDecaysOff=1;
	49	fEvaluate=0;
	50	fSelectAll=0;
	51	}
	52
	53	AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
	54	{
	55	// copy constructor
	56	}
	57
	58
	59	AliGenHijing::~AliGenHijing()
	60	{
	61	// Destructor
	62	}
	63
	64	void AliGenHijing::Init()
65	{
66	// Initialisation
67	SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
68	fAProjectile, fZProjectile, fATarget, fZTarget,
69	fMinImpactParam, fMaxImpactParam));
70
71	fHijing=(THijing*) fgMCEvGen;
72	fHijing->Initialize();
73	fHijing->SetIHPR2(3, fTrigger);
74	fHijing->SetIHPR2(4, fQuench);
75	fHijing->SetIHPR2(6, fShadowing);
76	fHijing->SetIHPR2(12, fDecaysOff);
77	fHijing->SetIHPR2(21, fKeep);
78	//
79	if (fEvaluate) EvaluateCrossSections();
80	}
81
82	void AliGenHijing::Generate()
83	{
84	// Generate one event
85
86	Float_t polar[3] = {0,0,0};
87	Float_t origin[3]= {0,0,0};
88	Float_t origin0[3]= {0,0,0};
89	Float_t p[3], random[6];
90	Float_t tof;
91
92	static TClonesArray *particles;
93	// converts from mm/c to s
94	const Float_t kconv=0.001/2.999792458e8;
95	//
96	Int_t nt=0;
97	Int_t jev=0;
98	Int_t j, kf, ks, imo;
99
100	if(!particles) particles=new TClonesArray("TParticle",10000);
101
102	fTrials=0;
103	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
104	if(fVertexSmear==perEvent) {
105	gMC->Rndm(random,6);
106	for (j=0;j<3;j++) {
107	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
108	TMath::Sqrt(-2TMath::Log(random[2j+1]));
109	// fHijing->SetMSTP(151,0);
110	}
111	} else if (fVertexSmear==perTrack) {
112	// fHijing->SetMSTP(151,0);
113	for (j=0;j<3;j++) {
114	// fHijing->SetPARP(151+j, fOsigma[j]*10.);
115	}
116	}
117	while(1)
118	{
119
120	fHijing->GenerateEvent();
121	fTrials++;
122	fHijing->ImportParticles(particles,"All");
123	Int_t np = particles->GetEntriesFast();
124	printf("\n **************************************************%d\n",np);
125	Int_t nc=0;
126	if (np == 0 ) continue;
127	Int_t i;
128	Int_t * newPos = new Int_t[np];
129	for (i = 0; i<np-1; i++) *(newPos+i)=i;
130
131	for (i = 0; i<np-1; i++) {
132	TParticle * iparticle = (TParticle *) particles->At(i);
133
134	Bool_t hasMother = (iparticle->GetFirstMother() >=0);
135	Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
136	Bool_t selected = kTRUE;
137	Bool_t hasSelectedDaughters = kFALSE;
138
139	if (!fSelectAll) selected = KinematicSelection(iparticle);
140	kf = iparticle->GetPdgCode();
141	// Int_t id1=iparticle->GetFirstDaughter();
142	// Int_t id2=iparticle->GetLastDaughter();
143	// printf("\n particle %d %d %d %d %d %d\n",i, kf, id1, id2, hasDaughter, selected);
144
145	if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
146	//
147	// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
148	//
149
150	if (selected \|\| hasSelectedDaughters) {
151	nc++;
152	ks = iparticle->GetStatusCode();
153	p[0]=iparticle->Px();
154	p[1]=iparticle->Py();
155	p[2]=iparticle->Pz();
156	origin[0]=origin0[0]+iparticle->Vx()/10;
157	origin[1]=origin0[1]+iparticle->Vy()/10;
158	origin[2]=origin0[2]+iparticle->Vz()/10;
159	tof=kconv*iparticle->T();
160	imo=-1;
161	if (hasMother) {
162	imo=iparticle->GetFirstMother();
163	imo=*(newPos+imo);
164	}
165
166	// printf("\n selected iparent %d %d %d \n",i, kf, imo);
167	if (hasDaughter) {
168	gAlice->SetTrack(0,imo,kf,p,origin,polar,
169	tof,"Primary",nt);
170	} else {
171	gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
172	tof,"Secondary",nt);
173	}
174	*(newPos+i)=nt;
175	} // selected
176	} // particle loop
177	delete newPos;
178	printf("\n I've put %i particles on the stack \n",nc);
179	if (nc > 0) {
180	jev+=nc;
181	if (jev >= fNpart \|\| fNpart == -1) {
182	fKineBias=Float_t(fNpart)/Float_t(fTrials);
183	printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
184	break;
185	}
186	}
187	} // event loop
188	}
189
190	Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
191	{
192	// Perform kinematic selection
193	Float_t px=particle->Px();
194	Float_t py=particle->Py();
195	Float_t pz=particle->Pz();
196	Float_t e=particle->Energy();
197
198	//
199	// transverse momentum cut
200	Float_t pt=TMath::Sqrt(pxpx+pypy);
201	if (pt > fPtMax \|\| pt < fPtMin)
202	{
203	// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
204	return kFALSE;
205	}
206	//
207	// momentum cut
208	Float_t p=TMath::Sqrt(pxpx+pypy+pz*pz);
209	if (p > fPMax \|\| p < fPMin)
210	{
211	// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
212	return kFALSE;
213	}
214
215	//
216	// theta cut
217	Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
218	if (theta > fThetaMax \|\| theta < fThetaMin)
219	{
220
221	// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
222	return kFALSE;
223	}
224
225	//
226	// rapidity cut
227	Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
228	if (y > fYMax \|\| y < fYMin)
229	{
230	// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
231	return kFALSE;
232	}
233
234	//
235	// phi cut
236	Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
237	if (phi > fPhiMax \|\| phi < fPhiMin)
238	{
239	// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
240	return kFALSE;
241	}
242
243	return kTRUE;
244	}
245
246	void AliGenHijing::KeepFullEvent()
247	{
248	fKeep=1;
249	}
250
251	void AliGenHijing::EvaluateCrossSections()
252	{
253	// Glauber Calculation of geometrical x-section
254	//
255	Float_t xTot=0.; // barn
256	Float_t xTotHard=0.; // barn
257	Float_t xPart=0.; // barn
258	Float_t xPartHard=0.; // barn
259	Float_t sigmaHard=0.1; // mbarn
260	Float_t bMin=0.;
261	Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
262	const Float_t kdib=0.2;
263	Int_t kMax=Int_t((bMax-bMin)/kdib)+1;
264
265
266	printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
267	printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
268	Int_t i;
269	Float_t oldvalue=0.;
270
271	for (i=0; i<kMax; i++)
272	{
273	Float_t xb=bMin+i*kdib;
274	Float_t ov;
275	ov=fHijing->Profile(xb);
276	Float_t gb = 2.0.01fHijing->GetHIPR1(40)kdibxb(1.-TMath::Exp(-fHijing->GetHINT1(12)ov));
277	Float_t gbh = 2.0.01fHijing->GetHIPR1(40)kdibxbsigmaHardov;
278	xTot+=gb;
279	xTotHard+=gbh;
280	if (xb > fMinImpactParam && xb < fMaxImpactParam)
281	{
282	xPart+=gb;
283	xPartHard+=gbh;
284	}
285
286	if ((xTot-oldvalue)/oldvalue<0.0001) break;
287	oldvalue=xTot;
288	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
289	printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
290	}
291	printf("\n Total cross section (barn): %f \n",xTot);
292	printf("\n Hard cross section (barn): %f \n \n",xTotHard);
293	printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
294	printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
295	}
296
297	Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
298	{
299	//
300	// Looks recursively if one of the daughters has been selected
301	//
302	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
303	Int_t imin=-1;
304	Int_t imax=-1;
305	Int_t i;
306	Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
307	Bool_t selected=kFALSE;
308	if (hasDaughters) {
309	imin=iparticle->GetFirstDaughter();
310	imax=iparticle->GetLastDaughter();
311	// printf("\n Has daughters %d %d:", imin, imax);
312	for (i=imin; i<= imax; i++){
313	TParticle * jparticle = (TParticle *) particles->At(i);
314	// Int_t ip=jparticle->GetPdgCode();
315	// printf("\n consider daughter %d %d", i,ip);
316
317	if (KinematicSelection(jparticle)) {selected=kTRUE; break;}
318	if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
319	}
320	} else {
321	return kFALSE;
322	}
323
324	return selected;
325	}
326
327
328
329	AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
330	{
331	// Assignment operator
332	return *this;
333	}