Updates for mixing.
[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.),
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)
c9a8628a 66{
67// Default constructor
c9a8628a 68 fCollisionGeometry = 0;
69}
70
71AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
1c56e311 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())
c9a8628a 96{
97// Constructor
98 fName = "SlowNucleons";
99 fTitle = "Generator for gray particles in pA collisions";
c9a8628a 100 fCollisionGeometry = 0;
c9a8628a 101}
102
5a9b7f8e 103//____________________________________________________________
c9a8628a 104AliGenSlowNucleons::~AliGenSlowNucleons()
105{
106// Destructor
107 delete fSlowNucleonModel;
108}
109
83b3cd6c 110void 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}
c9a8628a 114
115void 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) {
c8f7f6f9 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.);
5a9b7f8e 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.);
c9a8628a 135 }
136}
137
138void AliGenSlowNucleons::FinishRun()
139{
f687abfe 140// End of run action
141// Show histogram for debugging if requested.
c9a8628a 142 if (fDebug) {
c8f7f6f9 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();
5a9b7f8e 149 c->cd(3);
150 fCosThetaGrayHist->Draw();
c9a8628a 151 }
152}
153
154
155void AliGenSlowNucleons::Generate()
156{
157 //
158 // Generate one event
159 //
160 //
161 // Communication with Gray Particle Model
162 //
1c22a81e 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
99f5114f 170 fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 171 if (fDebug) {
172 printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
c8f7f6f9 173 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
174 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
1c22a81e 175 printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
176 b, nn, nwn, nnw, nwnw);
177 }
178 }
c9a8628a 179
180 //
5a9b7f8e 181 Float_t p[3], theta=0;
c9a8628a 182 Float_t origin[3] = {0., 0., 0.};
183 Float_t polar [3] = {0., 0., 0.};
99f5114f 184 Int_t nt, i, j;
c9a8628a 185 Int_t kf;
99f5114f 186
187 if(fVertexSmear == kPerEvent) {
188 Vertex();
189 for (j=0; j < 3; j++) origin[j] = fVertex[j];
190 } // if kPerEvent
c9a8628a 191//
192// Gray protons
193//
194 fCharge = 1;
195 kf = kProton;
1c22a81e 196 for(i = 0; i < fNgp; i++) {
5a9b7f8e 197 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
198 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 199 PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a 200 0., kPNoProcess, nt, 1.);
201 KeepTrack(nt);
202 }
203//
204// Gray neutrons
205//
206 fCharge = 0;
207 kf = kNeutron;
1c22a81e 208 for(i = 0; i < fNgn; i++) {
5a9b7f8e 209 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
210 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 211 PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a 212 0., kPNoProcess, nt, 1.);
213 KeepTrack(nt);
214 }
215//
216// Black protons
217//
218 fCharge = 1;
219 kf = kProton;
1c22a81e 220 for(i = 0; i < fNbp; i++) {
5a9b7f8e 221 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 222 PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a 223 0., kPNoProcess, nt, 1.);
224 KeepTrack(nt);
225 }
226//
227// Black neutrons
228//
229 fCharge = 0;
230 kf = kNeutron;
1c22a81e 231 for(i = 0; i < fNbn; i++) {
5a9b7f8e 232 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 233 PushTrack(fTrackIt, -1, kf, p, origin, polar,
c9a8628a 234 0., kPNoProcess, nt, 1.);
235 KeepTrack(nt);
236 }
237}
238
239
240
241
5a9b7f8e 242void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
243 Double_t beta, Float_t* q, Float_t &theta)
f687abfe 244
245{
c9a8628a 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
5a9b7f8e 252 Double_t m, pmax, p, f, phi;
c9a8628a 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(2*T*(T+sqrt(T*T+m*m)));
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
5a9b7f8e 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 //
c9a8628a 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);
83b3cd6c 297 q[2] *= fProtonDirection;
c9a8628a 298}
299
300Double_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(p*p+m*m)-m;
305 return (p*p * exp(-ekin/T));
306}
307
308
309void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
310{
311/* Lorentz transform in the direction of q[2] */
312
60e55aee 313 Double_t gamma = 1/sqrt((1.-beta)*(1.+beta));
c9a8628a 314 Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
315 q[2] = gamma * (q[2] + beta*energy);
316}