Introducing the interaction time into the aliroot generators. In case of gaussian...
[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.};
21391258 183 Float_t time = 0.;
c9a8628a 184 Float_t polar [3] = {0., 0., 0.};
99f5114f 185 Int_t nt, i, j;
c9a8628a 186 Int_t kf;
99f5114f 187
188 if(fVertexSmear == kPerEvent) {
189 Vertex();
190 for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258 191 time = fTime;
99f5114f 192 } // if kPerEvent
c9a8628a 193//
194// Gray protons
195//
196 fCharge = 1;
197 kf = kProton;
1c22a81e 198 for(i = 0; i < fNgp; i++) {
5a9b7f8e 199 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
200 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 201 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 202 time, kPNoProcess, nt, 1.);
c9a8628a 203 KeepTrack(nt);
204 }
205//
206// Gray neutrons
207//
208 fCharge = 0;
209 kf = kNeutron;
1c22a81e 210 for(i = 0; i < fNgn; i++) {
5a9b7f8e 211 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
212 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 213 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 214 time, kPNoProcess, nt, 1.);
c9a8628a 215 KeepTrack(nt);
216 }
217//
218// Black protons
219//
220 fCharge = 1;
221 kf = kProton;
1c22a81e 222 for(i = 0; i < fNbp; i++) {
5a9b7f8e 223 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 224 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 225 time, kPNoProcess, nt, 1.);
c9a8628a 226 KeepTrack(nt);
227 }
228//
229// Black neutrons
230//
231 fCharge = 0;
232 kf = kNeutron;
1c22a81e 233 for(i = 0; i < fNbn; i++) {
5a9b7f8e 234 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 235 PushTrack(fTrackIt, -1, kf, p, origin, polar,
21391258 236 time, kPNoProcess, nt, 1.);
c9a8628a 237 KeepTrack(nt);
238 }
239}
240
241
242
243
5a9b7f8e 244void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
245 Double_t beta, Float_t* q, Float_t &theta)
f687abfe 246
247{
c9a8628a 248/*
249 Emit a slow nucleon with "temperature" T [GeV],
250 from a source moving with velocity beta
251 Three-momentum [GeV/c] is given back in q[3]
252*/
253
5a9b7f8e 254 Double_t m, pmax, p, f, phi;
c9a8628a 255 TDatabasePDG * pdg = TDatabasePDG::Instance();
256 const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
257 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
258
259 /* Select nucleon type */
260 if(charge == 0) m = kMassNeutron;
261 else m = kMassProton;
262
263 /* Momentum at maximum of Maxwell-distribution */
264
265 pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
266
267 /* Try until proper momentum */
268 /* for lack of primitive function of the Maxwell-distribution */
269 /* a brute force trial-accept loop, normalized at pmax */
270
271 do
272 {
273 p = Rndm() * fPmax;
274 f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
275 }
276 while(f < Rndm());
277
5a9b7f8e 278 /* Spherical symmetric emission for black particles (beta=0)*/
279 if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
280 /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
281 else if(fThetaDistribution!=0){
282 Double_t costheta = fCosTheta->GetRandom();
283 theta = TMath::ACos(costheta);
284 }
285 //
c9a8628a 286 phi = 2. * TMath::Pi() * Rndm();
287
288 /* Determine momentum components in system of the moving source */
289 q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
290 q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
291 q[2] = p * TMath::Cos(theta);
292
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);
296
297 /* Transform to laboratory system */
298 Lorentz(m, fBeta, q);
83b3cd6c 299 q[2] *= fProtonDirection;
c9a8628a 300}
301
302Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
303{
304/* Relativistic Maxwell-distribution */
305 Double_t ekin;
306 ekin = TMath::Sqrt(p*p+m*m)-m;
307 return (p*p * exp(-ekin/T));
308}
309
310
311void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
312{
313/* Lorentz transform in the direction of q[2] */
314
60e55aee 315 Double_t gamma = 1/sqrt((1.-beta)*(1.+beta));
c9a8628a 316 Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
317 q[2] = gamma * (q[2] + beta*energy);
318}