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