New volumes in the thermal shield (M. Sitta)
[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx
CommitLineData
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>
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 38ClassImp(AliGenSlowNucleons)
1c56e311 39
c9a8628a 40
1c56e311 41AliGenSlowNucleons::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{
65// Default constructor
c9a8628a 66 fCollisionGeometry = 0;
67}
68
69AliGenSlowNucleons::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
5a9b7f8e 99//____________________________________________________________
c9a8628a 100AliGenSlowNucleons::~AliGenSlowNucleons()
101{
102// Destructor
103 delete fSlowNucleonModel;
104}
105
83b3cd6c 106void 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
111void AliGenSlowNucleons::Init()
112{
113 //
114 // Initialization
115 //
116 Float_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);
119 if (fDebug) {
c8f7f6f9 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.);
5a9b7f8e 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.);
c9a8628a 131 }
132}
133
134void AliGenSlowNucleons::FinishRun()
135{
f687abfe 136// End of run action
137// Show histogram for debugging if requested.
c9a8628a 138 if (fDebug) {
c8f7f6f9 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();
5a9b7f8e 145 c->cd(3);
146 fCosThetaGrayHist->Draw();
c9a8628a 147 }
148}
149
150
151void AliGenSlowNucleons::Generate()
152{
153 //
154 // Generate one event
155 //
156 //
157 // Communication with Gray Particle Model
158 //
1c22a81e 159 if (fCollisionGeometry) {
160 Float_t b = fCollisionGeometry->ImpactParameter();
230d85c6 161 // Int_t nn = fCollisionGeometry->NN();
162 // Int_t nwn = fCollisionGeometry->NwN();
163 // Int_t nnw = fCollisionGeometry->NNw();
164 // Int_t nwnw = fCollisionGeometry->NwNw();
1c22a81e 165
99f5114f 166 fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 167 if (fDebug) {
230d85c6 168 printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
c8f7f6f9 169 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
170 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
230d85c6 171 //printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
1c22a81e 172 }
173 }
c9a8628a 174
175 //
5a9b7f8e 176 Float_t p[3], theta=0;
c9a8628a 177 Float_t origin[3] = {0., 0., 0.};
21391258 178 Float_t time = 0.;
c9a8628a 179 Float_t polar [3] = {0., 0., 0.};
99f5114f 180 Int_t nt, i, j;
c9a8628a 181 Int_t kf;
99f5114f 182
183 if(fVertexSmear == kPerEvent) {
184 Vertex();
185 for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258 186 time = fTime;
99f5114f 187 } // if kPerEvent
c9a8628a 188//
189// Gray protons
190//
191 fCharge = 1;
192 kf = kProton;
1c22a81e 193 for(i = 0; i < fNgp; i++) {
5a9b7f8e 194 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
195 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 196 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 197 time, kPNoProcess, nt, 1.);
c9a8628a 198 KeepTrack(nt);
199 }
200//
201// Gray neutrons
202//
203 fCharge = 0;
204 kf = kNeutron;
1c22a81e 205 for(i = 0; i < fNgn; i++) {
5a9b7f8e 206 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
207 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 208 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 209 time, kPNoProcess, nt, 1.);
c9a8628a 210 KeepTrack(nt);
211 }
212//
213// Black protons
214//
215 fCharge = 1;
216 kf = kProton;
1c22a81e 217 for(i = 0; i < fNbp; i++) {
5a9b7f8e 218 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 219 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 220 time, kPNoProcess, nt, 1.);
c9a8628a 221 KeepTrack(nt);
222 }
223//
224// Black neutrons
225//
226 fCharge = 0;
227 kf = kNeutron;
1c22a81e 228 for(i = 0; i < fNbn; i++) {
5a9b7f8e 229 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 230 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 231 time, kPNoProcess, nt, 1.);
c9a8628a 232 KeepTrack(nt);
233 }
234}
235
236
237
238
5a9b7f8e 239void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
240 Double_t beta, Float_t* q, Float_t &theta)
f687abfe 241
242{
c9a8628a 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
230d85c6 249 //printf("Generating slow nuc. with: charge %d. temp. %1.3f, beta %1.2f\n",charge,T,beta);
250
5a9b7f8e 251 Double_t m, pmax, p, f, phi;
c9a8628a 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(2*T*(T+sqrt(T*T+m*m)));
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
5a9b7f8e 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 //
c9a8628a 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);
83b3cd6c 296 q[2] *= fProtonDirection;
c9a8628a 297}
298
299Double_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(p*p+m*m)-m;
304 return (p*p * exp(-ekin/T));
305}
306
307
308void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
309{
310/* Lorentz transform in the direction of q[2] */
311
61cf877f 312 Double_t gamma = 1./sqrt((1.-beta)*(1.+beta));
c9a8628a 313 Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
314 q[2] = gamma * (q[2] + beta*energy);
315}
230d85c6 316