Update to subtract fixed values of v2,3 (Redmer)
[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 //
88a9f852 116 Double_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);
88a9f852 119 //printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
c9a8628a 120 if (fDebug) {
c8f7f6f9 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.);
5a9b7f8e 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.);
c9a8628a 132 }
133}
134
135void AliGenSlowNucleons::FinishRun()
136{
f687abfe 137// End of run action
138// Show histogram for debugging if requested.
c9a8628a 139 if (fDebug) {
c8f7f6f9 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();
5a9b7f8e 146 c->cd(3);
147 fCosThetaGrayHist->Draw();
c9a8628a 148 }
149}
150
151
152void AliGenSlowNucleons::Generate()
153{
154 //
155 // Generate one event
156 //
157 //
158 // Communication with Gray Particle Model
159 //
1c22a81e 160 if (fCollisionGeometry) {
161 Float_t b = fCollisionGeometry->ImpactParameter();
230d85c6 162 // Int_t nn = fCollisionGeometry->NN();
163 // Int_t nwn = fCollisionGeometry->NwN();
164 // Int_t nnw = fCollisionGeometry->NNw();
165 // Int_t nwnw = fCollisionGeometry->NwNw();
88a9f852 166
99f5114f 167 fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 168 if (fDebug) {
88a9f852 169 //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
230d85c6 170 printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
c8f7f6f9 171 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
172 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
1c22a81e 173 }
174 }
c9a8628a 175
176 //
88a9f852 177 Float_t p[3] = {0., 0., 0.}, theta=0;
c9a8628a 178 Float_t origin[3] = {0., 0., 0.};
21391258 179 Float_t time = 0.;
c9a8628a 180 Float_t polar [3] = {0., 0., 0.};
99f5114f 181 Int_t nt, i, j;
c9a8628a 182 Int_t kf;
99f5114f 183
184 if(fVertexSmear == kPerEvent) {
185 Vertex();
186 for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258 187 time = fTime;
99f5114f 188 } // if kPerEvent
c9a8628a 189//
190// Gray protons
191//
192 fCharge = 1;
193 kf = kProton;
1c22a81e 194 for(i = 0; i < fNgp; i++) {
5a9b7f8e 195 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
196 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 197 PushTrack(fTrackIt, -1, kf, p, origin, polar,
1ffc7bfa 198 time, kPNoProcess, nt, 1.,1);
c9a8628a 199 KeepTrack(nt);
200 }
201//
202// Gray neutrons
203//
204 fCharge = 0;
205 kf = kNeutron;
1c22a81e 206 for(i = 0; i < fNgn; i++) {
5a9b7f8e 207 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
208 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 209 PushTrack(fTrackIt, -1, kf, p, origin, polar,
1ffc7bfa 210 time, kPNoProcess, nt, 1.,1);
c9a8628a 211 KeepTrack(nt);
212 }
213//
214// Black protons
215//
216 fCharge = 1;
217 kf = kProton;
1c22a81e 218 for(i = 0; i < fNbp; i++) {
5a9b7f8e 219 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 220 PushTrack(fTrackIt, -1, kf, p, origin, polar,
1ffc7bfa 221 time, kPNoProcess, nt, 1.,1);
c9a8628a 222 KeepTrack(nt);
223 }
224//
225// Black neutrons
226//
227 fCharge = 0;
228 kf = kNeutron;
1c22a81e 229 for(i = 0; i < fNbn; i++) {
5a9b7f8e 230 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 231 PushTrack(fTrackIt, -1, kf, p, origin, polar,
1ffc7bfa 232 time, kPNoProcess, nt, 1.,1);
c9a8628a 233 KeepTrack(nt);
234 }
235}
236
237
238
239
5a9b7f8e 240void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
241 Double_t beta, Float_t* q, Float_t &theta)
f687abfe 242
243{
c9a8628a 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
88a9f852 250 //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
230d85c6 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
88a9f852 263 pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
c9a8628a 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
88a9f852 286
c9a8628a 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
88a9f852 292
c9a8628a 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);
88a9f852 296
c9a8628a 297 /* Transform to laboratory system */
298 Lorentz(m, fBeta, q);
83b3cd6c 299 q[2] *= fProtonDirection;
88a9f852 300
c9a8628a 301}
302
303Double_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(p*p+m*m)-m;
308 return (p*p * exp(-ekin/T));
309}
310
311
312void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
313{
314/* Lorentz transform in the direction of q[2] */
315
88a9f852 316 Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
317 Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
c9a8628a 318 q[2] = gamma * (q[2] + beta*energy);
88a9f852 319 //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
c9a8628a 320}
230d85c6 321