[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.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 for slow nucleons in pA interactions. 
//  Source is modelled by a relativistic Maxwell distributions.
//  This class cooparates with AliCollisionGeometry if used inside AliGenCocktail.
//  In this case the number of slow nucleons is determined from the number of wounded nuclei
//  using a realisation of AliSlowNucleonModel.
//  Original code by  Ferenc Sikler  <sikler@rmki.kfki.hu>
// 

#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TF1.h>
#include <TCanvas.h>

#include "AliCollisionGeometry.h"
#include "AliGenSlowNucleons.h"
#include "AliSlowNucleonModel.h"

ClassImp(AliGenSlowNucleons)

    
AliGenSlowNucleons::AliGenSlowNucleons()
    :AliGenerator(-1),
     fCMS(0.),
     fMomentum(0.),
     fBeta(0.),
     fPmax (0.),
     fCharge(0),
     fProtonDirection(1.),
     fTemperatureG(0.), 
     fBetaSourceG(0.),
     fTemperatureB(0.),
     fBetaSourceB(0.),
     fNgp(0),
     fNgn(0),
     fNbp(0),
     fNbn(0),
     fDebug(0),
     fDebugHist1(0),
     fDebugHist2(0),
     fThetaDistribution(),
     fCosThetaGrayHist(),
     fCosTheta(),
     fSlowNucleonModel(0)
{
// Default constructor
    fCollisionGeometry = 0;
}

AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
    :AliGenerator(npart),
     fCMS(14000.),
     fMomentum(0.),
     fBeta(0.),
     fPmax (10.),
     fCharge(1),
     fProtonDirection(1.),
     fTemperatureG(0.05), 
     fBetaSourceG(0.05),
     fTemperatureB(0.005),
     fBetaSourceB(0.),
     fNgp(0),
     fNgn(0),
     fNbp(0),
     fNbn(0),
     fDebug(0),
     fDebugHist1(0),
     fDebugHist2(0),
     fThetaDistribution(),
     fCosThetaGrayHist(),
     fCosTheta(),
     fSlowNucleonModel(new AliSlowNucleonModel())
{
// Constructor
    fName  = "SlowNucleons";
    fTitle = "Generator for gray particles in pA collisions";
    fCollisionGeometry = 0;
}

//____________________________________________________________
AliGenSlowNucleons::~AliGenSlowNucleons()
{
// Destructor
    delete  fSlowNucleonModel;
}

void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
// Set direction of the proton to change between pA (1) and Ap (-1)
  fProtonDirection = dir / TMath::Abs(dir);
}

void AliGenSlowNucleons::Init()
{
  //
  // Initialization
  //
    Double_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
    fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
    fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
    //printf("  fMomentum %f    fBeta %1.10f\n",fMomentum, fBeta);
    if (fDebug) {
	fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
	fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
    	fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
    }
    //
    // non-uniform cos(theta) distribution
    //
    if(fThetaDistribution != 0) {
	fCosTheta = new TF1("fCosTheta",
			    "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
			    -1., 1.);
    }
}

void AliGenSlowNucleons::FinishRun()
{
// End of run action
// Show histogram for debugging if requested.
    if (fDebug) {
	TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
	c->Divide(2,1);
	c->cd(1);
	fDebugHist1->Draw("colz");
	c->cd(2);
	fDebugHist2->Draw();
	c->cd(3);
	fCosThetaGrayHist->Draw();
    }
}


void AliGenSlowNucleons::Generate()
{
  //
  // Generate one event
  //
  //
  // Communication with Gray Particle Model 
  // 
    if (fCollisionGeometry) {
	Float_t b   = fCollisionGeometry->ImpactParameter();
	//	Int_t  nn   = fCollisionGeometry->NN();
	//      Int_t  nwn  = fCollisionGeometry->NwN();
	//      Int_t  nnw  = fCollisionGeometry->NNw();
	//      Int_t  nwnw = fCollisionGeometry->NwNw();
	
	//fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	if (fDebug) {
	    //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
	    printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
	    fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
	    fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
	}
    }     

   //
    Float_t p[3] = {0., 0., 0.}, theta=0;
    Float_t origin[3] = {0., 0., 0.};
    Float_t time = 0.;
    Float_t polar [3] = {0., 0., 0.};    
    Int_t nt, i, j;
    Int_t kf;
    
    if(fVertexSmear == kPerEvent) {
	Vertex();
	for (j=0; j < 3; j++) origin[j] = fVertex[j];
	time = fTime;
    } // if kPerEvent
//
//  Gray protons
//
    fCharge = 1;
    kf = kProton;    
    for(i = 0; i < fNgp; i++) {
	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-2);
	KeepTrack(nt);
    }
//
//  Gray neutrons
//
    fCharge = 0;
    kf = kNeutron;    
    for(i = 0; i < fNgn; i++) {
	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-2);
	KeepTrack(nt);
    }
//
//  Black protons
//
    fCharge = 1;
    kf = kProton;    
    for(i = 0; i < fNbp; i++) {
	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-1);
	KeepTrack(nt);
    }
//
//  Black neutrons
//
    fCharge = 0;
    kf = kNeutron;    
    for(i = 0; i < fNbn; i++) {
	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-1);
	KeepTrack(nt);
    }
}


void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, 
	Double_t beta, Float_t* q, Float_t &theta)

{
/* 
   Emit a slow nucleon with "temperature" T [GeV], 
   from a source moving with velocity beta         
   Three-momentum [GeV/c] is given back in q[3]    
*/

 //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
 
 Double_t m, pmax, p, f, phi;
 TDatabasePDG * pdg = TDatabasePDG::Instance();
 const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 
 /* Select nucleon type */
 if(charge == 0) m = kMassNeutron;
 else m = kMassProton;

 /* Momentum at maximum of Maxwell-distribution */

 pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));

 /* Try until proper momentum                                  */
 /* for lack of primitive function of the Maxwell-distribution */
 /* a brute force trial-accept loop, normalized at pmax        */

 do
 {
     p = Rndm() * fPmax;
     f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 }
 while(f < Rndm());

 /* Spherical symmetric emission for black particles (beta=0)*/
 if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 else if(fThetaDistribution!=0){
   Double_t costheta = fCosTheta->GetRandom();
   theta = TMath::ACos(costheta);
 }
 //
 phi   = 2. * TMath::Pi() * Rndm();


 /* Determine momentum components in system of the moving source */
 q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 q[2] = p * TMath::Cos(theta);


 /* Transform to system of the target nucleus                             */
 /* beta is passed as negative, because the gray nucleons are slowed down */
 Lorentz(m, -beta, q);
 
 /* Transform to laboratory system */
 Lorentz(m, fBeta, q);
 q[2] *= fProtonDirection; 

}

Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
{
/* Relativistic Maxwell-distribution */
    Double_t ekin;
    ekin = TMath::Sqrt(p*p+m*m)-m;
    return (p*p * exp(-ekin/T));
}


void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
{
/* Lorentz transform in the direction of q[2] */
 
    Double_t gamma  = 1./TMath::Sqrt((1.-beta)*(1.+beta)); 
    Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    q[2] = gamma * (q[2] + beta*energy);
    //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
}
Commit	Line	Data
c9a8628a	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
88cb7938	16	/* $Id$ */
c9a8628a	17
f687abfe	18	//
	19	// Generator for slow nucleons in pA interactions.
	20	// Source is modelled by a relativistic Maxwell distributions.
	21	// This class cooparates with AliCollisionGeometry if used inside AliGenCocktail.
	22	// In this case the number of slow nucleons is determined from the number of wounded nuclei
	23	// using a realisation of AliSlowNucleonModel.
	24	// Original code by Ferenc Sikler <sikler@rmki.kfki.hu>
	25	//
c9a8628a	26
	27	#include <TDatabasePDG.h>
	28	#include <TPDGCode.h>
	29	#include <TH2F.h>
5a9b7f8e	30	#include <TH1F.h>
5a9b7f8e	31	#include <TF1.h>
c8f7f6f9	32	#include <TCanvas.h>
c9a8628a	33
	34	#include "AliCollisionGeometry.h"
	35	#include "AliGenSlowNucleons.h"
	36	#include "AliSlowNucleonModel.h"
	37
706938e6	38	ClassImp(AliGenSlowNucleons)
1c56e311	39
c9a8628a	40
1c56e311	41	AliGenSlowNucleons::AliGenSlowNucleons()
	42	:AliGenerator(-1),
	43	fCMS(0.),
	44	fMomentum(0.),
	45	fBeta(0.),
	46	fPmax (0.),
1c56e311	47	fCharge(0),
61cf877f	48	fProtonDirection(1.),
1c56e311	49	fTemperatureG(0.),
	50	fBetaSourceG(0.),
	51	fTemperatureB(0.),
	52	fBetaSourceB(0.),
	53	fNgp(0),
	54	fNgn(0),
	55	fNbp(0),
	56	fNbn(0),
	57	fDebug(0),
	58	fDebugHist1(0),
	59	fDebugHist2(0),
	60	fThetaDistribution(),
	61	fCosThetaGrayHist(),
	62	fCosTheta(),
	63	fSlowNucleonModel(0)
c9a8628a	64	{
c9a8628a	65	// Default constructor
c9a8628a	66	fCollisionGeometry = 0;
	67	}
	68
	69	AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
1c56e311	70	:AliGenerator(npart),
	71	fCMS(14000.),
	72	fMomentum(0.),
	73	fBeta(0.),
	74	fPmax (10.),
1c56e311	75	fCharge(1),
61cf877f	76	fProtonDirection(1.),
230d85c6	77	fTemperatureG(0.05),
1c56e311	78	fBetaSourceG(0.05),
230d85c6	79	fTemperatureB(0.005),
1c56e311	80	fBetaSourceB(0.),
	81	fNgp(0),
	82	fNgn(0),
	83	fNbp(0),
	84	fNbn(0),
	85	fDebug(0),
	86	fDebugHist1(0),
	87	fDebugHist2(0),
	88	fThetaDistribution(),
	89	fCosThetaGrayHist(),
	90	fCosTheta(),
	91	fSlowNucleonModel(new AliSlowNucleonModel())
c9a8628a	92	{
	93	// Constructor
	94	fName = "SlowNucleons";
	95	fTitle = "Generator for gray particles in pA collisions";
c9a8628a	96	fCollisionGeometry = 0;
c9a8628a	97	}
	98
	99	//____________________________________________________________
	100	AliGenSlowNucleons::~AliGenSlowNucleons()
	101	{
	102	// Destructor
	103	delete fSlowNucleonModel;
	104	}
	105
83b3cd6c	106	void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
	107	// Set direction of the proton to change between pA (1) and Ap (-1)
	108	fProtonDirection = dir / TMath::Abs(dir);
	109	}
c9a8628a	110
	111	void AliGenSlowNucleons::Init()
	112	{
	113	//
	114	// Initialization
	115	//
88a9f852	116	Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
271e6aae	117	fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
c9a8628a	118	fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
88a9f852	119	//printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
c9a8628a	120	if (fDebug) {
c8f7f6f9	121	fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
c8f7f6f9	122	fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
5a9b7f8e	123	fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
	124	}
	125	//
	126	// non-uniform cos(theta) distribution
	127	//
	128	if(fThetaDistribution != 0) {
	129	fCosTheta = new TF1("fCosTheta",
	130	"(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))exp(2.x/3.14159265358979312)",
	131	-1., 1.);
c9a8628a	132	}
	133	}
	134
	135	void AliGenSlowNucleons::FinishRun()
	136	{
f687abfe	137	// End of run action
f687abfe	138	// Show histogram for debugging if requested.
c9a8628a	139	if (fDebug) {
c8f7f6f9	140	TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
	141	c->Divide(2,1);
	142	c->cd(1);
58776b75	143	fDebugHist1->Draw("colz");
c8f7f6f9	144	c->cd(2);
c8f7f6f9	145	fDebugHist2->Draw();
5a9b7f8e	146	c->cd(3);
5a9b7f8e	147	fCosThetaGrayHist->Draw();
c9a8628a	148	}
	149	}
	150
	151
	152	void AliGenSlowNucleons::Generate()
	153	{
	154	//
	155	// Generate one event
	156	//
	157	//
	158	// Communication with Gray Particle Model
	159	//
1c22a81e	160	if (fCollisionGeometry) {
1c22a81e	161	Float_t b = fCollisionGeometry->ImpactParameter();
230d85c6	162	// Int_t nn = fCollisionGeometry->NN();
	163	// Int_t nwn = fCollisionGeometry->NwN();
	164	// Int_t nnw = fCollisionGeometry->NNw();
	165	// Int_t nwnw = fCollisionGeometry->NwNw();
88a9f852	166
58776b75	167	//fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
58776b75	168	fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e	169	if (fDebug) {
88a9f852	170	//printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
230d85c6	171	printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
58776b75	172	fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
c8f7f6f9	173	fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
1c22a81e	174	}
1c22a81e	175	}
c9a8628a	176
c9a8628a	177	//
88a9f852	178	Float_t p[3] = {0., 0., 0.}, theta=0;
c9a8628a	179	Float_t origin[3] = {0., 0., 0.};
21391258	180	Float_t time = 0.;
c9a8628a	181	Float_t polar [3] = {0., 0., 0.};
99f5114f	182	Int_t nt, i, j;
c9a8628a	183	Int_t kf;
99f5114f	184
	185	if(fVertexSmear == kPerEvent) {
	186	Vertex();
	187	for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258	188	time = fTime;
99f5114f	189	} // if kPerEvent
c9a8628a	190	//
	191	// Gray protons
	192	//
	193	fCharge = 1;
	194	kf = kProton;
1c22a81e	195	for(i = 0; i < fNgp; i++) {
5a9b7f8e	196	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
5a9b7f8e	197	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf	198	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	199	time, kPNoProcess, nt, 1.,-2);
c9a8628a	200	KeepTrack(nt);
	201	}
	202	//
	203	// Gray neutrons
	204	//
	205	fCharge = 0;
	206	kf = kNeutron;
1c22a81e	207	for(i = 0; i < fNgn; i++) {
5a9b7f8e	208	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
5a9b7f8e	209	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf	210	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	211	time, kPNoProcess, nt, 1.,-2);
c9a8628a	212	KeepTrack(nt);
	213	}
	214	//
	215	// Black protons
	216	//
	217	fCharge = 1;
	218	kf = kProton;
1c22a81e	219	for(i = 0; i < fNbp; i++) {
5a9b7f8e	220	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf	221	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	222	time, kPNoProcess, nt, 1.,-1);
c9a8628a	223	KeepTrack(nt);
	224	}
	225	//
	226	// Black neutrons
	227	//
	228	fCharge = 0;
	229	kf = kNeutron;
1c22a81e	230	for(i = 0; i < fNbn; i++) {
5a9b7f8e	231	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf	232	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	233	time, kPNoProcess, nt, 1.,-1);
c9a8628a	234	KeepTrack(nt);
	235	}
	236	}
	237
	238
	239
	240
5a9b7f8e	241	void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
5a9b7f8e	242	Double_t beta, Float_t* q, Float_t &theta)
f687abfe	243
f687abfe	244	{
c9a8628a	245	/*
	246	Emit a slow nucleon with "temperature" T [GeV],
	247	from a source moving with velocity beta
	248	Three-momentum [GeV/c] is given back in q[3]
	249	*/
	250
88a9f852	251	//printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
230d85c6	252
5a9b7f8e	253	Double_t m, pmax, p, f, phi;
c9a8628a	254	TDatabasePDG * pdg = TDatabasePDG::Instance();
	255	const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
	256	const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
	257
	258	/* Select nucleon type */
	259	if(charge == 0) m = kMassNeutron;
	260	else m = kMassProton;
	261
	262	/* Momentum at maximum of Maxwell-distribution */
	263
88a9f852	264	pmax = TMath::Sqrt(2T(T+TMath::Sqrt(TT+mm)));
c9a8628a	265
	266	/* Try until proper momentum */
	267	/* for lack of primitive function of the Maxwell-distribution */
	268	/* a brute force trial-accept loop, normalized at pmax */
	269
	270	do
	271	{
	272	p = Rndm() * fPmax;
	273	f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
	274	}
	275	while(f < Rndm());
	276
5a9b7f8e	277	/* Spherical symmetric emission for black particles (beta=0)*/
	278	if(beta==0 \|\| fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
	279	/* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
	280	else if(fThetaDistribution!=0){
	281	Double_t costheta = fCosTheta->GetRandom();
	282	theta = TMath::ACos(costheta);
	283	}
	284	//
c9a8628a	285	phi = 2. * TMath::Pi() * Rndm();
c9a8628a	286
88a9f852	287
c9a8628a	288	/* Determine momentum components in system of the moving source */
	289	q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
	290	q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
	291	q[2] = p * TMath::Cos(theta);
	292
88a9f852	293
c9a8628a	294	/* Transform to system of the target nucleus */
	295	/* beta is passed as negative, because the gray nucleons are slowed down */
	296	Lorentz(m, -beta, q);
88a9f852	297
c9a8628a	298	/* Transform to laboratory system */
c9a8628a	299	Lorentz(m, fBeta, q);
83b3cd6c	300	q[2] *= fProtonDirection;
88a9f852	301
c9a8628a	302	}
	303
	304	Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
	305	{
	306	/* Relativistic Maxwell-distribution */
	307	Double_t ekin;
	308	ekin = TMath::Sqrt(pp+mm)-m;
	309	return (pp exp(-ekin/T));
	310	}
	311
	312
	313	void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
	314	{
	315	/* Lorentz transform in the direction of q[2] */
	316
88a9f852	317	Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
88a9f852	318	Double_t energy = TMath::Sqrt(mm + q[0]q[0] + q[1]q[1] + q[2]q[2]);
c9a8628a	319	q[2] = gamma * (q[2] + beta*energy);
88a9f852	320	//printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
c9a8628a	321	}
230d85c6	322