[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.),
     fATarget (0.),
     fZTarget (0.),
     fCharge(0),
     fProtonDirection(0.),
     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.),
     fATarget (208.),
     fZTarget (82.),
     fCharge(1),
     fProtonDirection(0.),
     fTemperatureG(0.04), 
     fBetaSourceG(0.05),
     fTemperatureB(0.004),
     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(const AliGenSlowNucleons & sn):
    AliGenerator(sn),
    fCMS(0.),
    fMomentum(0.),
    fBeta(0.),
    fPmax (0.),
    fATarget (0.),
    fZTarget (0.),
    fCharge(0),
    fProtonDirection(0.),
    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)
{
// Copy constructor
    sn.Copy(*this);
}

//____________________________________________________________
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
  //
    Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
    fMomentum = fCMS/2. * fZTarget / fATarget;
    fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
    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();
	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);
	if (fDebug) {
	    printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
	    fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
	    fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
	    printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", 
		   b, nn, nwn, nnw, nwnw);
	}
    }     

   //
    Float_t p[3], theta=0;
    Float_t origin[3] = {0., 0., 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];
    } // 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,
		 0., kPNoProcess, nt, 1.);
	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,
		 0., kPNoProcess, nt, 1.);
	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,
		 0., kPNoProcess, nt, 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,
		 0., kPNoProcess, nt, 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]    
*/

 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+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/sqrt(1-beta*beta); 
    Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    q[2] = gamma * (q[2] + beta*energy);
}

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

void AliGenSlowNucleons::Copy(TObject&) const
{
    //
    // Copy 
    //
    Fatal("Copy","Not implemented!\n");
}
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.),
	47	fATarget (0.),
	48	fZTarget (0.),
	49	fCharge(0),
	50	fProtonDirection(0.),
	51	fTemperatureG(0.),
	52	fBetaSourceG(0.),
	53	fTemperatureB(0.),
	54	fBetaSourceB(0.),
	55	fNgp(0),
	56	fNgn(0),
	57	fNbp(0),
	58	fNbn(0),
	59	fDebug(0),
	60	fDebugHist1(0),
	61	fDebugHist2(0),
	62	fThetaDistribution(),
	63	fCosThetaGrayHist(),
	64	fCosTheta(),
	65	fSlowNucleonModel(0)
c9a8628a	66	{
c9a8628a	67	// Default constructor
c9a8628a	68	fCollisionGeometry = 0;
	69	}
	70
	71	AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
1c56e311	72	:AliGenerator(npart),
	73	fCMS(14000.),
	74	fMomentum(0.),
	75	fBeta(0.),
	76	fPmax (10.),
	77	fATarget (208.),
	78	fZTarget (82.),
	79	fCharge(1),
	80	fProtonDirection(0.),
	81	fTemperatureG(0.04),
	82	fBetaSourceG(0.05),
	83	fTemperatureB(0.004),
	84	fBetaSourceB(0.),
	85	fNgp(0),
	86	fNgn(0),
	87	fNbp(0),
	88	fNbn(0),
	89	fDebug(0),
	90	fDebugHist1(0),
	91	fDebugHist2(0),
	92	fThetaDistribution(),
	93	fCosThetaGrayHist(),
	94	fCosTheta(),
	95	fSlowNucleonModel(new AliSlowNucleonModel())
c9a8628a	96	{
	97	// Constructor
	98	fName = "SlowNucleons";
	99	fTitle = "Generator for gray particles in pA collisions";
c9a8628a	100	fCollisionGeometry = 0;
c9a8628a	101	}
c9a8628a	102
5a9b7f8e	103	//____________________________________________________________
f687abfe	104	AliGenSlowNucleons::AliGenSlowNucleons(const AliGenSlowNucleons & sn):
1c56e311	105	AliGenerator(sn),
	106	fCMS(0.),
	107	fMomentum(0.),
	108	fBeta(0.),
	109	fPmax (0.),
	110	fATarget (0.),
	111	fZTarget (0.),
	112	fCharge(0),
	113	fProtonDirection(0.),
	114	fTemperatureG(0.),
	115	fBetaSourceG(0.),
	116	fTemperatureB(0.),
	117	fBetaSourceB(0.),
	118	fNgp(0),
	119	fNgn(0),
	120	fNbp(0),
	121	fNbn(0),
	122	fDebug(0),
	123	fDebugHist1(0),
	124	fDebugHist2(0),
	125	fThetaDistribution(),
	126	fCosThetaGrayHist(),
	127	fCosTheta(),
	128	fSlowNucleonModel(0)
f687abfe	129	{
	130	// Copy constructor
	131	sn.Copy(*this);
	132	}
	133
c9a8628a	134	//____________________________________________________________
	135	AliGenSlowNucleons::~AliGenSlowNucleons()
	136	{
	137	// Destructor
	138	delete fSlowNucleonModel;
	139	}
	140
83b3cd6c	141	void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
	142	// Set direction of the proton to change between pA (1) and Ap (-1)
	143	fProtonDirection = dir / TMath::Abs(dir);
	144	}
c9a8628a	145
	146	void AliGenSlowNucleons::Init()
	147	{
	148	//
	149	// Initialization
	150	//
	151	Float_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
	152	fMomentum = fCMS/2. * fZTarget / fATarget;
	153	fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
	154	if (fDebug) {
c8f7f6f9	155	fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
c8f7f6f9	156	fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
5a9b7f8e	157	fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
	158	}
	159	//
	160	// non-uniform cos(theta) distribution
	161	//
	162	if(fThetaDistribution != 0) {
	163	fCosTheta = new TF1("fCosTheta",
	164	"(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))exp(2.x/3.14159265358979312)",
	165	-1., 1.);
c9a8628a	166	}
	167	}
	168
	169	void AliGenSlowNucleons::FinishRun()
	170	{
f687abfe	171	// End of run action
f687abfe	172	// Show histogram for debugging if requested.
c9a8628a	173	if (fDebug) {
c8f7f6f9	174	TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
	175	c->Divide(2,1);
	176	c->cd(1);
	177	fDebugHist1->Draw();
	178	c->cd(2);
	179	fDebugHist2->Draw();
5a9b7f8e	180	c->cd(3);
5a9b7f8e	181	fCosThetaGrayHist->Draw();
c9a8628a	182	}
	183	}
	184
	185
	186	void AliGenSlowNucleons::Generate()
	187	{
	188	//
	189	// Generate one event
	190	//
	191	//
	192	// Communication with Gray Particle Model
	193	//
1c22a81e	194	if (fCollisionGeometry) {
	195	Float_t b = fCollisionGeometry->ImpactParameter();
	196	Int_t nn = fCollisionGeometry->NN();
	197	Int_t nwn = fCollisionGeometry->NwN();
	198	Int_t nnw = fCollisionGeometry->NNw();
	199	Int_t nwnw = fCollisionGeometry->NwNw();
	200
99f5114f	201	fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e	202	if (fDebug) {
1c22a81e	203	printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
c8f7f6f9	204	fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
c8f7f6f9	205	fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
1c22a81e	206	printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
	207	b, nn, nwn, nnw, nwnw);
	208	}
	209	}
c9a8628a	210
c9a8628a	211	//
5a9b7f8e	212	Float_t p[3], theta=0;
c9a8628a	213	Float_t origin[3] = {0., 0., 0.};
c9a8628a	214	Float_t polar [3] = {0., 0., 0.};
99f5114f	215	Int_t nt, i, j;
c9a8628a	216	Int_t kf;
99f5114f	217
	218	if(fVertexSmear == kPerEvent) {
	219	Vertex();
	220	for (j=0; j < 3; j++) origin[j] = fVertex[j];
	221	} // if kPerEvent
c9a8628a	222	//
	223	// Gray protons
	224	//
	225	fCharge = 1;
	226	kf = kProton;
1c22a81e	227	for(i = 0; i < fNgp; i++) {
5a9b7f8e	228	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
5a9b7f8e	229	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf	230	PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a	231	0., kPNoProcess, nt, 1.);
	232	KeepTrack(nt);
	233	}
	234	//
	235	// Gray neutrons
	236	//
	237	fCharge = 0;
	238	kf = kNeutron;
1c22a81e	239	for(i = 0; i < fNgn; i++) {
5a9b7f8e	240	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
5a9b7f8e	241	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf	242	PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a	243	0., kPNoProcess, nt, 1.);
	244	KeepTrack(nt);
	245	}
	246	//
	247	// Black protons
	248	//
	249	fCharge = 1;
	250	kf = kProton;
1c22a81e	251	for(i = 0; i < fNbp; i++) {
5a9b7f8e	252	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf	253	PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a	254	0., kPNoProcess, nt, 1.);
	255	KeepTrack(nt);
	256	}
	257	//
	258	// Black neutrons
	259	//
	260	fCharge = 0;
	261	kf = kNeutron;
1c22a81e	262	for(i = 0; i < fNbn; i++) {
5a9b7f8e	263	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf	264	PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a	265	0., kPNoProcess, nt, 1.);
	266	KeepTrack(nt);
	267	}
	268	}
	269
	270
	271
	272
5a9b7f8e	273	void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
5a9b7f8e	274	Double_t beta, Float_t* q, Float_t &theta)
f687abfe	275
f687abfe	276	{
c9a8628a	277	/*
	278	Emit a slow nucleon with "temperature" T [GeV],
	279	from a source moving with velocity beta
	280	Three-momentum [GeV/c] is given back in q[3]
	281	*/
	282
5a9b7f8e	283	Double_t m, pmax, p, f, phi;
c9a8628a	284	TDatabasePDG * pdg = TDatabasePDG::Instance();
	285	const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
	286	const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
	287
	288	/* Select nucleon type */
	289	if(charge == 0) m = kMassNeutron;
	290	else m = kMassProton;
	291
	292	/* Momentum at maximum of Maxwell-distribution */
	293
	294	pmax = TMath::Sqrt(2T(T+sqrt(TT+mm)));
	295
	296	/* Try until proper momentum */
	297	/* for lack of primitive function of the Maxwell-distribution */
	298	/* a brute force trial-accept loop, normalized at pmax */
	299
	300	do
	301	{
	302	p = Rndm() * fPmax;
	303	f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
	304	}
	305	while(f < Rndm());
	306
5a9b7f8e	307	/* Spherical symmetric emission for black particles (beta=0)*/
	308	if(beta==0 \|\| fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
	309	/* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
	310	else if(fThetaDistribution!=0){
	311	Double_t costheta = fCosTheta->GetRandom();
	312	theta = TMath::ACos(costheta);
	313	}
	314	//
c9a8628a	315	phi = 2. * TMath::Pi() * Rndm();
	316
	317	/* Determine momentum components in system of the moving source */
	318	q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
	319	q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
	320	q[2] = p * TMath::Cos(theta);
	321
	322	/* Transform to system of the target nucleus */
	323	/* beta is passed as negative, because the gray nucleons are slowed down */
	324	Lorentz(m, -beta, q);
	325
	326	/* Transform to laboratory system */
	327	Lorentz(m, fBeta, q);
83b3cd6c	328	q[2] *= fProtonDirection;
c9a8628a	329	}
	330
	331	Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
	332	{
	333	/* Relativistic Maxwell-distribution */
	334	Double_t ekin;
	335	ekin = TMath::Sqrt(pp+mm)-m;
	336	return (pp exp(-ekin/T));
	337	}
	338
	339
	340	void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
	341	{
	342	/* Lorentz transform in the direction of q[2] */
	343
	344	Double_t gamma = 1/sqrt(1-beta*beta);
	345	Double_t energy = sqrt(mm + q[0]q[0] + q[1]q[1] + q[2]q[2]);
	346	q[2] = gamma * (q[2] + beta*energy);
	347	}
	348
f687abfe	349
	350	AliGenSlowNucleons& AliGenSlowNucleons::operator=(const AliGenSlowNucleons& rhs)
	351	{
	352	// Assignment operator
	353	rhs.Copy(*this);
	354	return *this;
	355	}
	356
dc1d768c	357	void AliGenSlowNucleons::Copy(TObject&) const
f687abfe	358	{
	359	//
	360	// Copy
	361	//
	362	Fatal("Copy","Not implemented!\n");
	363	}
	364
c9a8628a	365
	366
	367
	368
	369
	370