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