[u/mrichter/AliRoot.git] / EVGEN / AliGenHIJINGpara.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.13  2002/10/14 14:55:35  hristov
Merging the VirtualMC branch to the main development branch (HEAD)

Revision 1.11.4.1  2002/07/24 08:56:28  alibrary
Updating EVGEN on TVirtulaMC

Revision 1.12  2002/06/19 06:56:34  hristov
Memory leak corrected

Revision 1.11  2002/03/20 10:21:13  hristov
Set fPtMax to 15 GeV in order to avoid some numerical problems

Revision 1.10  2001/10/15 16:44:46  morsch
- Possibility for vertex distribution truncation.
- Write mc header with vertex position.

Revision 1.9  2001/07/27 17:09:36  morsch
Use local SetTrack, KeepTrack and SetHighWaterMark methods
to delegate either to local stack or to stack owned by AliRun.
(Piotr Skowronski, A.M.)

Revision 1.8  2001/07/20 11:03:58  morsch
Issue warning message if used outside allowed eta range (-8 to 8).

Revision 1.7  2001/07/17 12:41:01  morsch
- Calculation of fraction of event corresponding to selected pt-range corrected
(R. Turrisi)
- Parent weight corrected.

Revision 1.6  2001/05/16 14:57:10  alibrary
New files for folders and Stack

Revision 1.5  2000/12/21 16:24:06  morsch
Coding convention clean-up

Revision 1.4  2000/11/30 07:12:50  alibrary
Introducing new Rndm and QA classes

Revision 1.3  2000/10/02 21:28:06  fca
Removal of useless dependecies via forward declarations

Revision 1.2  2000/07/11 18:24:55  fca
Coding convention corrections + few minor bug fixes

Revision 1.1  2000/06/09 20:20:30  morsch
Same class as previously in AliSimpleGen.cxx
All coding rule violations except RS3 corrected (AM)

*/

// Parameterisation of pi and K, eta and pt distributions
// used for the ALICE TDRs.
// eta: according to HIJING (shadowing + quenching)
// pT : according to CDF measurement at 1.8 TeV
// Author: andreas.morsch@cern.ch


//Begin_Html
/*
<img src="picts/AliGeneratorClass.gif">
</pre>
<br clear=left>
<font size=+2 color=red>
<p>The responsible person for this module is
<a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
</font>
<pre>
*/
//End_Html
//                                                               //
///////////////////////////////////////////////////////////////////

#include "AliGenHIJINGpara.h"
#include "AliGenEventHeader.h"
#include "AliRun.h"
#include "AliConst.h"
#include "AliDecayer.h"
#include "AliDecayerPythia.h"
#include "AliPDG.h"

#include <TF1.h>
#include <TArrayF.h>
#include <TDatabasePDG.h>
#include <TClonesArray.h>
#include <TParticle.h>

ClassImp(AliGenHIJINGpara)

AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
{
// copy constructor
}

//_____________________________________________________________________________
static Double_t ptpi(Double_t *px, Double_t *)
{
  //
  //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  //     POWER LAW FOR PT > 500 MEV
  //     MT SCALING BELOW (T=160 MEV)
  //
  const Double_t kp0 = 1.3;
  const Double_t kxn = 8.28;
  const Double_t kxlim=0.5;
  const Double_t kt=0.160;
  const Double_t kxmpi=0.139;
  const Double_t kb=1.;
  Double_t y, y1, xmpi2, ynorm, a;
  Double_t x=*px;
  //
  y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  xmpi2=kxmpi*kxmpi;
  ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  a=ynorm/y1;
  if (x > kxlim)
    y=a*TMath::Power(kp0/(kp0+x),kxn);
  else
    y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  return y*x;
}

//_____________________________________________________________________________
static Double_t ptscal(Double_t pt, Int_t np)
{
    //    SCALING EN MASSE PAR RAPPORT A PTPI
    //     MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
    const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
    //     VALUE MESON/PI AT 5 GEV
    const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
    np--;
    Double_t f5=TMath::Power(((
	sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
    Double_t fmax2=f5/kfmax[np];
    // PIONS
    Double_t ptpion=100.*ptpi(&pt, (Double_t*) 0);
    Double_t fmtscal=TMath::Power(((
	sqrt(pt*pt+0.018215)+2.)/ (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ 
	fmax2;
    return fmtscal*ptpion;
}

//_____________________________________________________________________________
static Double_t ptka( Double_t *px, Double_t *)
{
    //
    // pt parametrisation for k
    //
    return ptscal(*px,2);
}


//_____________________________________________________________________________
static Double_t etapic( Double_t *py, Double_t *)
{
  //
  // eta parametrisation for pi
  //
    const Double_t ka1    = 4913.;
    const Double_t ka2    = 1819.;
    const Double_t keta1  = 0.22;
    const Double_t keta2  = 3.66;
    const Double_t kdeta1 = 1.47;
    const Double_t kdeta2 = 1.51;
    Double_t y=TMath::Abs(*py);
    //
    Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
    Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
    return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
}

//_____________________________________________________________________________
static Double_t etakac( Double_t *py, Double_t *)
{
    //
    // eta parametrisation for ka
    //
    const Double_t ka1    = 497.6;
    const Double_t ka2    = 215.6;
    const Double_t keta1  = 0.79;
    const Double_t keta2  = 4.09;
    const Double_t kdeta1 = 1.54;
    const Double_t kdeta2 = 1.40;
    Double_t y=TMath::Abs(*py);
    //
    Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
    Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
    return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
}

//_____________________________________________________________________________
AliGenHIJINGpara::AliGenHIJINGpara()
  :AliGenerator()
{
    //
    // Default constructor
    //
    fPtpi    = 0;
    fPtka    = 0;
    fETApic  = 0;
    fETAkac  = 0;
    fDecayer = 0;
    SetCutVertexZ();
    SetPtRange();
    SetPi0Decays();
}

//_____________________________________________________________________________
AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
  :AliGenerator(npart)
{
  // 
  // Standard constructor
  //
    fName="HIGINGpara";
    fTitle="HIJING Parametrisation Particle Generator";
    fPtpi    = 0;
    fPtka    = 0;
    fETApic  = 0;
    fETAkac  = 0;
    fDecayer = 0;
    SetCutVertexZ();
    SetPtRange();
    SetPi0Decays();
}

//_____________________________________________________________________________
AliGenHIJINGpara::~AliGenHIJINGpara()
{
  //
  // Standard destructor
  //
    delete fPtpi;
    delete fPtka;
    delete fETApic;
    delete fETAkac;
}

//_____________________________________________________________________________
void AliGenHIJINGpara::Init()
{
  //
  // Initialise the HIJING parametrisation
  //
    Float_t etaMin =-TMath::Log(TMath::Tan(
	TMath::Min((Double_t)fThetaMax/2,TMath::Pi()/2-1.e-10)));
    Float_t etaMax = -TMath::Log(TMath::Tan(
	TMath::Max((Double_t)fThetaMin/2,1.e-10)));
    fPtpi   = new TF1("ptpi",&ptpi,0,20,0);
    fPtka   = new TF1("ptka",&ptka,0,20,0);
    fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
    fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);

    TF1 etaPic0("etapic",&etapic,-7,7,0);
    TF1 etaKac0("etakac",&etakac,-7,7,0);

    TF1 ptPic0("ptpi",&ptpi,0.,15.,0);
    TF1 ptKac0("ptka",&ptka,0.,15.,0);

    Float_t intETApi  = etaPic0.Integral(-0.5, 0.5);
    Float_t intETAka  = etaKac0.Integral(-0.5, 0.5);
    Float_t scalePi   = 7316/(intETApi/1.5);
    Float_t scaleKa   =  684/(intETAka/2.0);

//  Fraction of events corresponding to the selected pt-range    
    Float_t intPt    = (0.877*ptPic0.Integral(0, 15)+
			0.123*ptKac0.Integral(0, 15));
    Float_t intPtSel = (0.877*ptPic0.Integral(fPtMin, fPtMax)+
			0.123*ptKac0.Integral(fPtMin, fPtMax));
    Float_t ptFrac   = intPtSel/intPt;

//  Fraction of events corresponding to the selected eta-range    
    Float_t intETASel  = (scalePi*etaPic0.Integral(etaMin, etaMax)+
			  scaleKa*etaKac0.Integral(etaMin, etaMax));
//  Fraction of events corresponding to the selected phi-range    
    Float_t phiFrac    = (fPhiMax-fPhiMin)/2/TMath::Pi();

    fParentWeight = Float_t(fNpart)/(intETASel*ptFrac*phiFrac);
    
    printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", 
	   ClassName(),100.*fParentWeight);

// Issue warning message if etaMin or etaMax are outside the alowed range 
// of the parametrization
    if (etaMin < -8.001 || etaMax > 8.001) {
	printf("\n \n WARNING FROM AliGenHIJINGPara !");
	printf("\n YOU ARE USING THE PARAMETERISATION OUTSIDE ");	
	printf("\n THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)");	    
	printf("\n YOUR LIMITS: %f %f \n \n ", etaMin, etaMax);
    }
//
//
    if (fPi0Decays) fDecayer = new AliDecayerPythia();
}


//_____________________________________________________________________________
void AliGenHIJINGpara::Generate()
{
  //
  // Generate one trigger
  //

  
    const Float_t kRaKpic=0.14;
    const Float_t kBorne=1/(1+kRaKpic);
    Float_t polar[3]= {0,0,0};
    //
    const Int_t kPions[3] = {kPi0, kPiPlus, kPiMinus};
    const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
    //
    Float_t origin[3];
    Float_t pt, pl, ptot;
    Float_t phi, theta;
    Float_t p[3];
    Int_t i, part, nt, j;
    //
    TF1 *ptf;
    TF1 *etaf;
    //
    Float_t random[6];
    //
    for (j=0;j<3;j++) origin[j]=fOrigin[j];

    if(fVertexSmear == kPerEvent) {
	Float_t dv[3];
	dv[2] = 1.e10;
	while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
	    Rndm(random,6);
	    for (j=0; j < 3; j++) {
		dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
		    TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
	    }
	}
	for (j=0; j < 3; j++) origin[j] += dv[j];
    } // if kPerEvent
    TArrayF eventVertex;
    eventVertex.Set(3);
    eventVertex[0] = origin[0];
    eventVertex[1] = origin[1];
    eventVertex[2] = origin[2];

    for(i=0;i<fNpart;i++) {
	while(1) {
	    Rndm(random,3);
	    if(random[0]<kBorne) {
		part=kPions[Int_t (random[1]*3)];
		ptf=fPtpi;
		etaf=fETApic;
	    } else {
		part=kKaons[Int_t (random[1]*4)];
		ptf=fPtka;
		etaf=fETAkac;
	    }
	    phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
	    theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
	    if(theta<fThetaMin || theta>fThetaMax) continue;
	    pt=ptf->GetRandom();
	    pl=pt/TMath::Tan(theta);
	    ptot=TMath::Sqrt(pt*pt+pl*pl);
	    if(ptot<fPMin || ptot>fPMax) continue;
	    p[0]=pt*TMath::Cos(phi);
	    p[1]=pt*TMath::Sin(phi);
	    p[2]=pl;
	    if(fVertexSmear==kPerTrack) {
		Rndm(random,6);
		for (j=0;j<3;j++) {
		    origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
			TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
		}
	    }
	    SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,nt,fParentWeight);
//
//          Decay pi0 if requested
	    if (part == kPi0 && fPi0Decays) DecayPi0(origin, p);
	    break;
	}
    }
// Header
    AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
// Event Vertex
    header->SetPrimaryVertex(eventVertex);
    gAlice->SetGenEventHeader(header); 
}

AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
{
// Assignment operator
    return *this;
}

void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
    AliGenerator::SetPtRange(ptmin, ptmax);
}

void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p) 
{
//
//    Decay the pi0
//    and put decay products on the stack
//
    static TClonesArray *particles;
    if(!particles) particles = new TClonesArray("TParticle",1000);
//    
    const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
    Float_t       e     = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] +p[2] * p[2] * kMass * kMass);
//
//  Decay the pi0    
    TLorentzVector pmom(p[0], p[1], p[2], e);
    fDecayer->Decay(kPi0, &pmom);
//
// Put decay particles on the stack
//
    Int_t nt = 0;
    Float_t polar[3] = {0., 0., 0.};
    Int_t np = fDecayer->ImportParticles(particles);
    
    for (Int_t i = 1; i < np; i++)
    {
	TParticle* iParticle =  (TParticle *) particles->At(1);
	p[0]=iParticle->Px();
	p[1]=iParticle->Py();
	p[2]=iParticle->Pz();
	Int_t part = iParticle->GetPdgCode();
	SetTrack(fTrackIt, 0, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
    }
}
Commit	Line	Data
790bbabf	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$
2067f62c	18	Revision 1.13 2002/10/14 14:55:35 hristov
	19	Merging the VirtualMC branch to the main development branch (HEAD)
	20
b9d0a01d	21	Revision 1.11.4.1 2002/07/24 08:56:28 alibrary
	22	Updating EVGEN on TVirtulaMC
	23
	24	Revision 1.12 2002/06/19 06:56:34 hristov
	25	Memory leak corrected
	26
aecf26a3	27	Revision 1.11 2002/03/20 10:21:13 hristov
	28	Set fPtMax to 15 GeV in order to avoid some numerical problems
	29
6b236b6f	30	Revision 1.10 2001/10/15 16:44:46 morsch
	31	- Possibility for vertex distribution truncation.
	32	- Write mc header with vertex position.
	33
971816d4	34	Revision 1.9 2001/07/27 17:09:36 morsch
	35	Use local SetTrack, KeepTrack and SetHighWaterMark methods
	36	to delegate either to local stack or to stack owned by AliRun.
	37	(Piotr Skowronski, A.M.)
	38
a99cf51f	39	Revision 1.8 2001/07/20 11:03:58 morsch
	40	Issue warning message if used outside allowed eta range (-8 to 8).
	41
5b0e4c0f	42	Revision 1.7 2001/07/17 12:41:01 morsch
	43	- Calculation of fraction of event corresponding to selected pt-range corrected
	44	(R. Turrisi)
	45	- Parent weight corrected.
	46
f25521b2	47	Revision 1.6 2001/05/16 14:57:10 alibrary
	48	New files for folders and Stack
	49
9e1a0ddb	50	Revision 1.5 2000/12/21 16:24:06 morsch
	51	Coding convention clean-up
	52
675e9664	53	Revision 1.4 2000/11/30 07:12:50 alibrary
	54	Introducing new Rndm and QA classes
	55
65fb704d	56	Revision 1.3 2000/10/02 21:28:06 fca
	57	Removal of useless dependecies via forward declarations
	58
94de3818	59	Revision 1.2 2000/07/11 18:24:55 fca
	60	Coding convention corrections + few minor bug fixes
	61
aee8290b	62	Revision 1.1 2000/06/09 20:20:30 morsch
	63	Same class as previously in AliSimpleGen.cxx
	64	All coding rule violations except RS3 corrected (AM)
	65
790bbabf	66	*/
675e9664	67
	68	// Parameterisation of pi and K, eta and pt distributions
	69	// used for the ALICE TDRs.
	70	// eta: according to HIJING (shadowing + quenching)
	71	// pT : according to CDF measurement at 1.8 TeV
	72	// Author: andreas.morsch@cern.ch
	73
	74
790bbabf	75	//Begin_Html
	76	/*
	77	<img src="picts/AliGeneratorClass.gif">
	78	</pre>
	79	<br clear=left>
	80	<font size=+2 color=red>
	81	<p>The responsible person for this module is
	82	<a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
	83	</font>
	84	<pre>
	85	*/
	86	//End_Html
	87	// //
	88	///////////////////////////////////////////////////////////////////
	89
	90	#include "AliGenHIJINGpara.h"
971816d4	91	#include "AliGenEventHeader.h"
790bbabf	92	#include "AliRun.h"
790bbabf	93	#include "AliConst.h"
2067f62c	94	#include "AliDecayer.h"
2067f62c	95	#include "AliDecayerPythia.h"
790bbabf	96	#include "AliPDG.h"
790bbabf	97
971816d4	98	#include <TF1.h>
971816d4	99	#include <TArrayF.h>
2067f62c	100	#include <TDatabasePDG.h>
	101	#include <TClonesArray.h>
	102	#include <TParticle.h>
971816d4	103
790bbabf	104	ClassImp(AliGenHIJINGpara)
	105
	106	AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
	107	{
	108	// copy constructor
	109	}
	110
	111	//_____________________________________________________________________________
	112	static Double_t ptpi(Double_t px, Double_t )
	113	{
	114	//
	115	// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
	116	// POWER LAW FOR PT > 500 MEV
	117	// MT SCALING BELOW (T=160 MEV)
	118	//
	119	const Double_t kp0 = 1.3;
	120	const Double_t kxn = 8.28;
	121	const Double_t kxlim=0.5;
	122	const Double_t kt=0.160;
	123	const Double_t kxmpi=0.139;
	124	const Double_t kb=1.;
	125	Double_t y, y1, xmpi2, ynorm, a;
	126	Double_t x=*px;
	127	//
	128	y1=TMath::Power(kp0/(kp0+kxlim),kxn);
	129	xmpi2=kxmpi*kxmpi;
	130	ynorm=kb(TMath::Exp(-sqrt(kxlimkxlim+xmpi2)/kt));
	131	a=ynorm/y1;
	132	if (x > kxlim)
	133	y=a*TMath::Power(kp0/(kp0+x),kxn);
	134	else
	135	y=kbTMath::Exp(-sqrt(xx+xmpi2)/kt);
	136	return y*x;
	137	}
	138
	139	//_____________________________________________________________________________
	140	static Double_t ptscal(Double_t pt, Int_t np)
	141	{
	142	// SCALING EN MASSE PAR RAPPORT A PTPI
	143	// MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
	144	const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
	145	// VALUE MESON/PI AT 5 GEV
	146	const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
	147	np--;
	148	Double_t f5=TMath::Power(((
	149	sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
	150	Double_t fmax2=f5/kfmax[np];
	151	// PIONS
	152	Double_t ptpion=100.ptpi(&pt, (Double_t) 0);
	153	Double_t fmtscal=TMath::Power(((
	154	sqrt(ptpt+0.018215)+2.)/ (sqrt(ptpt+khm[np]*khm[np])+2.0)),12.3)/
	155	fmax2;
	156	return fmtscal*ptpion;
	157	}
	158
	159	//_____________________________________________________________________________
	160	static Double_t ptka( Double_t px, Double_t )
	161	{
	162	//
	163	// pt parametrisation for k
	164	//
	165	return ptscal(*px,2);
	166	}
	167
168
169	//_____________________________________________________________________________
170	static Double_t etapic( Double_t py, Double_t )
171	{
172	//
173	// eta parametrisation for pi
174	//
175	const Double_t ka1 = 4913.;
176	const Double_t ka2 = 1819.;
177	const Double_t keta1 = 0.22;
178	const Double_t keta2 = 3.66;
179	const Double_t kdeta1 = 1.47;
180	const Double_t kdeta2 = 1.51;
181	Double_t y=TMath::Abs(*py);
182	//
183	Double_t ex1 = (y-keta1)(y-keta1)/(2kdeta1*kdeta1);
184	Double_t ex2 = (y-keta2)(y-keta2)/(2kdeta2*kdeta2);
185	return ka1TMath::Exp(-ex1)+ka2TMath::Exp(-ex2);
186	}
187
188	//_____________________________________________________________________________
189	static Double_t etakac( Double_t py, Double_t )
190	{
191	//
192	// eta parametrisation for ka
193	//
194	const Double_t ka1 = 497.6;
195	const Double_t ka2 = 215.6;
196	const Double_t keta1 = 0.79;
197	const Double_t keta2 = 4.09;
198	const Double_t kdeta1 = 1.54;
199	const Double_t kdeta2 = 1.40;
200	Double_t y=TMath::Abs(*py);
201	//
202	Double_t ex1 = (y-keta1)(y-keta1)/(2kdeta1*kdeta1);
203	Double_t ex2 = (y-keta2)(y-keta2)/(2kdeta2*kdeta2);
204	return ka1TMath::Exp(-ex1)+ka2TMath::Exp(-ex2);
205	}
206
207	//_____________________________________________________________________________
208	AliGenHIJINGpara::AliGenHIJINGpara()
209	:AliGenerator()
210	{
211	//
212	// Default constructor
213	//
2067f62c	214	fPtpi = 0;
	215	fPtka = 0;
	216	fETApic = 0;
	217	fETAkac = 0;
	218	fDecayer = 0;
971816d4	219	SetCutVertexZ();
6b236b6f	220	SetPtRange();
2067f62c	221	SetPi0Decays();
790bbabf	222	}
	223
	224	//_____________________________________________________________________________
	225	AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
	226	:AliGenerator(npart)
	227	{
	228	//
	229	// Standard constructor
	230	//
	231	fName="HIGINGpara";
	232	fTitle="HIJING Parametrisation Particle Generator";
2067f62c	233	fPtpi = 0;
	234	fPtka = 0;
	235	fETApic = 0;
	236	fETAkac = 0;
	237	fDecayer = 0;
971816d4	238	SetCutVertexZ();
6b236b6f	239	SetPtRange();
2067f62c	240	SetPi0Decays();
790bbabf	241	}
	242
	243	//_____________________________________________________________________________
	244	AliGenHIJINGpara::~AliGenHIJINGpara()
	245	{
	246	//
	247	// Standard destructor
	248	//
	249	delete fPtpi;
	250	delete fPtka;
	251	delete fETApic;
	252	delete fETAkac;
	253	}
	254
	255	//_____________________________________________________________________________
	256	void AliGenHIJINGpara::Init()
	257	{
	258	//
	259	// Initialise the HIJING parametrisation
	260	//
	261	Float_t etaMin =-TMath::Log(TMath::Tan(
	262	TMath::Min((Double_t)fThetaMax/2,TMath::Pi()/2-1.e-10)));
	263	Float_t etaMax = -TMath::Log(TMath::Tan(
	264	TMath::Max((Double_t)fThetaMin/2,1.e-10)));
f25521b2	265	fPtpi = new TF1("ptpi",&ptpi,0,20,0);
f25521b2	266	fPtka = new TF1("ptka",&ptka,0,20,0);
790bbabf	267	fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
790bbabf	268	fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
f25521b2	269
aecf26a3	270	TF1 etaPic0("etapic",&etapic,-7,7,0);
aecf26a3	271	TF1 etaKac0("etakac",&etakac,-7,7,0);
f25521b2	272
aecf26a3	273	TF1 ptPic0("ptpi",&ptpi,0.,15.,0);
aecf26a3	274	TF1 ptKac0("ptka",&ptka,0.,15.,0);
f25521b2	275
aecf26a3	276	Float_t intETApi = etaPic0.Integral(-0.5, 0.5);
aecf26a3	277	Float_t intETAka = etaKac0.Integral(-0.5, 0.5);
f25521b2	278	Float_t scalePi = 7316/(intETApi/1.5);
	279	Float_t scaleKa = 684/(intETAka/2.0);
	280
	281	// Fraction of events corresponding to the selected pt-range
aecf26a3	282	Float_t intPt = (0.877*ptPic0.Integral(0, 15)+
	283	0.123*ptKac0.Integral(0, 15));
	284	Float_t intPtSel = (0.877*ptPic0.Integral(fPtMin, fPtMax)+
	285	0.123*ptKac0.Integral(fPtMin, fPtMax));
f25521b2	286	Float_t ptFrac = intPtSel/intPt;
	287
	288	// Fraction of events corresponding to the selected eta-range
aecf26a3	289	Float_t intETASel = (scalePi*etaPic0.Integral(etaMin, etaMax)+
aecf26a3	290	scaleKa*etaKac0.Integral(etaMin, etaMax));
f25521b2	291	// Fraction of events corresponding to the selected phi-range
	292	Float_t phiFrac = (fPhiMax-fPhiMin)/2/TMath::Pi();
	293
	294	fParentWeight = Float_t(fNpart)/(intETASelptFracphiFrac);
790bbabf	295
9e1a0ddb	296	printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ",
9e1a0ddb	297	ClassName(),100.*fParentWeight);
5b0e4c0f	298
	299	// Issue warning message if etaMin or etaMax are outside the alowed range
	300	// of the parametrization
	301	if (etaMin < -8.001 \|\| etaMax > 8.001) {
	302	printf("\n \n WARNING FROM AliGenHIJINGPara !");
	303	printf("\n YOU ARE USING THE PARAMETERISATION OUTSIDE ");
	304	printf("\n THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)");
	305	printf("\n YOUR LIMITS: %f %f \n \n ", etaMin, etaMax);
	306	}
2067f62c	307	//
	308	//
	309	if (fPi0Decays) fDecayer = new AliDecayerPythia();
790bbabf	310	}
790bbabf	311
2067f62c	312
790bbabf	313	//_____________________________________________________________________________
	314	void AliGenHIJINGpara::Generate()
	315	{
	316	//
	317	// Generate one trigger
	318	//
	319
	320
	321	const Float_t kRaKpic=0.14;
	322	const Float_t kBorne=1/(1+kRaKpic);
	323	Float_t polar[3]= {0,0,0};
	324	//
	325	const Int_t kPions[3] = {kPi0, kPiPlus, kPiMinus};
	326	const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
	327	//
	328	Float_t origin[3];
	329	Float_t pt, pl, ptot;
	330	Float_t phi, theta;
	331	Float_t p[3];
	332	Int_t i, part, nt, j;
	333	//
	334	TF1 *ptf;
	335	TF1 *etaf;
	336	//
	337	Float_t random[6];
	338	//
	339	for (j=0;j<3;j++) origin[j]=fOrigin[j];
971816d4	340
	341	if(fVertexSmear == kPerEvent) {
	342	Float_t dv[3];
	343	dv[2] = 1.e10;
	344	while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
	345	Rndm(random,6);
	346	for (j=0; j < 3; j++) {
	347	dv[j] = fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	348	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	349	}
790bbabf	350	}
971816d4	351	for (j=0; j < 3; j++) origin[j] += dv[j];
	352	} // if kPerEvent
	353	TArrayF eventVertex;
	354	eventVertex.Set(3);
	355	eventVertex[0] = origin[0];
	356	eventVertex[1] = origin[1];
	357	eventVertex[2] = origin[2];
	358
790bbabf	359	for(i=0;i<fNpart;i++) {
790bbabf	360	while(1) {
65fb704d	361	Rndm(random,3);
790bbabf	362	if(random[0]<kBorne) {
	363	part=kPions[Int_t (random[1]*3)];
	364	ptf=fPtpi;
2067f62c	365	etaf=fETApic;
790bbabf	366	} else {
	367	part=kKaons[Int_t (random[1]*4)];
	368	ptf=fPtka;
	369	etaf=fETAkac;
	370	}
	371	phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
	372	theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
	373	if(theta<fThetaMin \|\| theta>fThetaMax) continue;
	374	pt=ptf->GetRandom();
	375	pl=pt/TMath::Tan(theta);
	376	ptot=TMath::Sqrt(ptpt+plpl);
	377	if(ptot<fPMin \|\| ptot>fPMax) continue;
	378	p[0]=pt*TMath::Cos(phi);
	379	p[1]=pt*TMath::Sin(phi);
	380	p[2]=pl;
aee8290b	381	if(fVertexSmear==kPerTrack) {
65fb704d	382	Rndm(random,6);
790bbabf	383	for (j=0;j<3;j++) {
	384	origin[j]=fOrigin[j]+fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	385	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	386	}
	387	}
a99cf51f	388	SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,nt,fParentWeight);
2067f62c	389	//
	390	// Decay pi0 if requested
	391	if (part == kPi0 && fPi0Decays) DecayPi0(origin, p);
790bbabf	392	break;
	393	}
	394	}
971816d4	395	// Header
	396	AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
	397	// Event Vertex
	398	header->SetPrimaryVertex(eventVertex);
	399	gAlice->SetGenEventHeader(header);
790bbabf	400	}
	401
	402	AliGenHIJINGpara& AliGenHIJINGpara::operator=(const AliGenHIJINGpara& rhs)
	403	{
	404	// Assignment operator
	405	return *this;
	406	}
	407
6b236b6f	408	void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
2067f62c	409	AliGenerator::SetPtRange(ptmin, ptmax);
6b236b6f	410	}
790bbabf	411
2067f62c	412	void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p)
	413	{
	414	//
	415	// Decay the pi0
	416	// and put decay products on the stack
	417	//
	418	static TClonesArray *particles;
	419	if(!particles) particles = new TClonesArray("TParticle",1000);
	420	//
	421	const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
	422	Float_t e = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] +p[2] * p[2] * kMass * kMass);
	423	//
	424	// Decay the pi0
	425	TLorentzVector pmom(p[0], p[1], p[2], e);
	426	fDecayer->Decay(kPi0, &pmom);
	427	//
	428	// Put decay particles on the stack
	429	//
	430	Int_t nt = 0;
	431	Float_t polar[3] = {0., 0., 0.};
	432	Int_t np = fDecayer->ImportParticles(particles);
	433
	434	for (Int_t i = 1; i < np; i++)
	435	{
	436	TParticle* iParticle = (TParticle *) particles->At(1);
	437	p[0]=iParticle->Px();
	438	p[1]=iParticle->Py();
	439	p[2]=iParticle->Pz();
	440	Int_t part = iParticle->GetPdgCode();
	441	SetTrack(fTrackIt, 0, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
	442	}
	443	}