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