[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();
	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("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->NwN(), 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.,-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.,-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.,-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
	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	Double_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	//printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
	120	if (fDebug) {
	121	fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
	122	fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
	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.);
	132	}
	133	}
	134
	135	void AliGenSlowNucleons::FinishRun()
	136	{
	137	// End of run action
	138	// Show histogram for debugging if requested.
	139	if (fDebug) {
	140	TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
	141	c->Divide(2,1);
	142	c->cd(1);
	143	fDebugHist1->Draw();
	144	c->cd(2);
	145	fDebugHist2->Draw();
	146	c->cd(3);
	147	fCosThetaGrayHist->Draw();
	148	}
	149	}
	150
	151
	152	void AliGenSlowNucleons::Generate()
	153	{
	154	//
	155	// Generate one event
	156	//
	157	//
	158	// Communication with Gray Particle Model
	159	//
	160	if (fCollisionGeometry) {
	161	Float_t b = fCollisionGeometry->ImpactParameter();
	162	// Int_t nn = fCollisionGeometry->NN();
	163	// Int_t nwn = fCollisionGeometry->NwN();
	164	// Int_t nnw = fCollisionGeometry->NNw();
	165	// Int_t nwnw = fCollisionGeometry->NwNw();
	166
	167	fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	168	if (fDebug) {
	169	//printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
	170	printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
	171	fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
	172	fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
	173	}
	174	}
	175
	176	//
	177	Float_t p[3] = {0., 0., 0.}, theta=0;
	178	Float_t origin[3] = {0., 0., 0.};
	179	Float_t time = 0.;
	180	Float_t polar [3] = {0., 0., 0.};
	181	Int_t nt, i, j;
	182	Int_t kf;
	183
	184	if(fVertexSmear == kPerEvent) {
	185	Vertex();
	186	for (j=0; j < 3; j++) origin[j] = fVertex[j];
	187	time = fTime;
	188	} // if kPerEvent
	189	//
	190	// Gray protons
	191	//
	192	fCharge = 1;
	193	kf = kProton;
	194	for(i = 0; i < fNgp; i++) {
	195	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
	196	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
	197	PushTrack(fTrackIt, -1, kf, p, origin, polar,
	198	time, kPNoProcess, nt, 1.,-1);
	199	KeepTrack(nt);
	200	}
	201	//
	202	// Gray neutrons
	203	//
	204	fCharge = 0;
	205	kf = kNeutron;
	206	for(i = 0; i < fNgn; i++) {
	207	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
	208	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
	209	PushTrack(fTrackIt, -1, kf, p, origin, polar,
	210	time, kPNoProcess, nt, 1.,-1);
	211	KeepTrack(nt);
	212	}
	213	//
	214	// Black protons
	215	//
	216	fCharge = 1;
	217	kf = kProton;
	218	for(i = 0; i < fNbp; i++) {
	219	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
	220	PushTrack(fTrackIt, -1, kf, p, origin, polar,
	221	time, kPNoProcess, nt, 1.,-1);
	222	KeepTrack(nt);
	223	}
	224	//
	225	// Black neutrons
	226	//
	227	fCharge = 0;
	228	kf = kNeutron;
	229	for(i = 0; i < fNbn; i++) {
	230	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
	231	PushTrack(fTrackIt, -1, kf, p, origin, polar,
	232	time, kPNoProcess, nt, 1.,-1);
	233	KeepTrack(nt);
	234	}
	235	}
	236
	237
	238
	239
	240	void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
	241	Double_t beta, Float_t* q, Float_t &theta)
	242
	243	{
	244	/*
	245	Emit a slow nucleon with "temperature" T [GeV],
	246	from a source moving with velocity beta
	247	Three-momentum [GeV/c] is given back in q[3]
	248	*/
	249
	250	//printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
	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+TMath::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
	287	/* Determine momentum components in system of the moving source */
	288	q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
	289	q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
	290	q[2] = p * TMath::Cos(theta);
	291
	292
	293	/* Transform to system of the target nucleus */
	294	/* beta is passed as negative, because the gray nucleons are slowed down */
	295	Lorentz(m, -beta, q);
	296
	297	/* Transform to laboratory system */
	298	Lorentz(m, fBeta, q);
	299	q[2] *= fProtonDirection;
	300
	301	}
	302
	303	Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
	304	{
	305	/* Relativistic Maxwell-distribution */
	306	Double_t ekin;
	307	ekin = TMath::Sqrt(pp+mm)-m;
	308	return (pp exp(-ekin/T));
	309	}
	310
	311
	312	void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
	313	{
	314	/* Lorentz transform in the direction of q[2] */
	315
	316	Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
	317	Double_t energy = TMath::Sqrt(mm + q[0]q[0] + q[1]q[1] + q[2]q[2]);
	318	q[2] = gamma * (q[2] + beta*energy);
	319	//printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
	320	}
	321