Use fVertex instead of fEventVertex
[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
16/*
17$Log$
99f5114f 18Revision 1.2 2003/03/25 09:55:24 morsch
19Numbers of slow nucleons either from model or user set.
20
1c22a81e 21Revision 1.1 2003/03/24 16:49:23 morsch
22Slow nucleon generator and model.
23
c9a8628a 24*/
25
26/*
27 Generator for slow nucluons in pA interactions.
28 Source is modelled by a relativistic Maxwell distributions.
29 Original code by Ferenc Sikler <sikler@rmki.kfki.hu>
30 */
31
32#include <TDatabasePDG.h>
33#include <TPDGCode.h>
34#include <TH2F.h>
35
36#include "AliCollisionGeometry.h"
37#include "AliGenSlowNucleons.h"
38#include "AliSlowNucleonModel.h"
39
40 ClassImp(AliGenSlowNucleons)
41
42 AliGenSlowNucleons::AliGenSlowNucleons():AliGenerator(-1)
43{
44// Default constructor
45 fSlowNucleonModel = 0;
46 fCollisionGeometry = 0;
47}
48
49AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
50 :AliGenerator(npart)
51{
52// Constructor
53 fName = "SlowNucleons";
54 fTitle = "Generator for gray particles in pA collisions";
55 SetPmax();
56 SetTarget();
57 SetNominalCmsEnergy();
58 SetCharge();
59 SetTemperature();
60 SetBetaSource();
61 fSlowNucleonModel = new AliSlowNucleonModel();
62 fCollisionGeometry = 0;
63 fDebug = 0;
64}
65
66//____________________________________________________________
67AliGenSlowNucleons::~AliGenSlowNucleons()
68{
69// Destructor
70 delete fSlowNucleonModel;
71}
72
73
74void AliGenSlowNucleons::Init()
75{
76 //
77 // Initialization
78 //
79 Float_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
80 fMomentum = fCMS/2. * fZTarget / fATarget;
81 fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
82 if (fDebug) {
83 fDebugHist = new TH2F("DebugHist", "N_heavy vs nu", 50, 0., 50., 50, 0., 50.);
84 }
85}
86
87void AliGenSlowNucleons::FinishRun()
88{
89 if (fDebug) {
90 fDebugHist->Draw();
91 }
92}
93
94
95void AliGenSlowNucleons::Generate()
96{
97 //
98 // Generate one event
99 //
100 //
101 // Communication with Gray Particle Model
102 //
1c22a81e 103 if (fCollisionGeometry) {
104 Float_t b = fCollisionGeometry->ImpactParameter();
105 Int_t nn = fCollisionGeometry->NN();
106 Int_t nwn = fCollisionGeometry->NwN();
107 Int_t nnw = fCollisionGeometry->NNw();
108 Int_t nwnw = fCollisionGeometry->NwNw();
109
99f5114f 110 fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 111 if (fDebug) {
112 printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
113 fDebugHist->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
114 printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
115 b, nn, nwn, nnw, nwnw);
116 }
117 }
c9a8628a 118
119 //
120 Float_t p[3];
121 Float_t origin[3] = {0., 0., 0.};
122 Float_t polar [3] = {0., 0., 0.};
99f5114f 123 Int_t nt, i, j;
c9a8628a 124 Int_t kf;
99f5114f 125
126 if(fVertexSmear == kPerEvent) {
127 Vertex();
128 for (j=0; j < 3; j++) origin[j] = fVertex[j];
129 } // if kPerEvent
c9a8628a 130//
131// Gray protons
132//
133 fCharge = 1;
134 kf = kProton;
1c22a81e 135 for(i = 0; i < fNgp; i++) {
c9a8628a 136 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
137 SetTrack(fTrackIt, -1, kf, p, origin, polar,
138 0., kPNoProcess, nt, 1.);
139 KeepTrack(nt);
140 }
141//
142// Gray neutrons
143//
144 fCharge = 0;
145 kf = kNeutron;
1c22a81e 146 for(i = 0; i < fNgn; i++) {
c9a8628a 147 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
148 SetTrack(fTrackIt, -1, kf, p, origin, polar,
149 0., kPNoProcess, nt, 1.);
150 KeepTrack(nt);
151 }
152//
153// Black protons
154//
155 fCharge = 1;
156 kf = kProton;
1c22a81e 157 for(i = 0; i < fNbp; i++) {
c9a8628a 158 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
159 SetTrack(fTrackIt, -1, kf, p, origin, polar,
160 0., kPNoProcess, nt, 1.);
161 KeepTrack(nt);
162 }
163//
164// Black neutrons
165//
166 fCharge = 0;
167 kf = kNeutron;
1c22a81e 168 for(i = 0; i < fNbn; i++) {
c9a8628a 169 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
170 SetTrack(fTrackIt, -1, kf, p, origin, polar,
171 0., kPNoProcess, nt, 1.);
172 KeepTrack(nt);
173 }
174}
175
176
177
178
179void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, Double_t beta, Float_t* q)
180/*
181 Emit a slow nucleon with "temperature" T [GeV],
182 from a source moving with velocity beta
183 Three-momentum [GeV/c] is given back in q[3]
184*/
185
186{
187 Double_t m, pmax, p, f, theta, phi;
188
189 TDatabasePDG * pdg = TDatabasePDG::Instance();
190 const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
191 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
192
193 /* Select nucleon type */
194 if(charge == 0) m = kMassNeutron;
195 else m = kMassProton;
196
197 /* Momentum at maximum of Maxwell-distribution */
198
199 pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
200
201 /* Try until proper momentum */
202 /* for lack of primitive function of the Maxwell-distribution */
203 /* a brute force trial-accept loop, normalized at pmax */
204
205 do
206 {
207 p = Rndm() * fPmax;
208 f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
209 }
210 while(f < Rndm());
211
212 /* Spherical symmetric emission */
213 theta = TMath::ACos(2. * Rndm() - 1.);
214 phi = 2. * TMath::Pi() * Rndm();
215
216 /* Determine momentum components in system of the moving source */
217 q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
218 q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
219 q[2] = p * TMath::Cos(theta);
220
221 /* Transform to system of the target nucleus */
222 /* beta is passed as negative, because the gray nucleons are slowed down */
223 Lorentz(m, -beta, q);
224
225 /* Transform to laboratory system */
226 Lorentz(m, fBeta, q);
227}
228
229Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
230{
231/* Relativistic Maxwell-distribution */
232 Double_t ekin;
233 ekin = TMath::Sqrt(p*p+m*m)-m;
234 return (p*p * exp(-ekin/T));
235}
236
237
238void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
239{
240/* Lorentz transform in the direction of q[2] */
241
242 Double_t gamma = 1/sqrt(1-beta*beta);
243 Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
244 q[2] = gamma * (q[2] + beta*energy);
245}
246
247
248
249
250
251
252