Possibility to use multiplcity as weight in the correction framework task, needed...
[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.),
1c56e311 77 fTemperatureG(0.04),
78 fBetaSourceG(0.05),
79 fTemperatureB(0.004),
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 //
116 Float_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);
119 if (fDebug) {
c8f7f6f9 120 fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
121 fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
5a9b7f8e 122 fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
123 }
124 //
125 // non-uniform cos(theta) distribution
126 //
127 if(fThetaDistribution != 0) {
128 fCosTheta = new TF1("fCosTheta",
129 "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
130 -1., 1.);
c9a8628a 131 }
132}
133
134void AliGenSlowNucleons::FinishRun()
135{
f687abfe 136// End of run action
137// Show histogram for debugging if requested.
c9a8628a 138 if (fDebug) {
c8f7f6f9 139 TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
140 c->Divide(2,1);
141 c->cd(1);
142 fDebugHist1->Draw();
143 c->cd(2);
144 fDebugHist2->Draw();
5a9b7f8e 145 c->cd(3);
146 fCosThetaGrayHist->Draw();
c9a8628a 147 }
148}
149
150
151void AliGenSlowNucleons::Generate()
152{
153 //
154 // Generate one event
155 //
156 //
157 // Communication with Gray Particle Model
158 //
1c22a81e 159 if (fCollisionGeometry) {
160 Float_t b = fCollisionGeometry->ImpactParameter();
161 Int_t nn = fCollisionGeometry->NN();
162 Int_t nwn = fCollisionGeometry->NwN();
163 Int_t nnw = fCollisionGeometry->NNw();
164 Int_t nwnw = fCollisionGeometry->NwNw();
165
99f5114f 166 fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 167 if (fDebug) {
168 printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
c8f7f6f9 169 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
170 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
1c22a81e 171 printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
172 b, nn, nwn, nnw, nwnw);
173 }
174 }
c9a8628a 175
176 //
5a9b7f8e 177 Float_t p[3], 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,
21391258 198 time, kPNoProcess, nt, 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,
21391258 210 time, kPNoProcess, nt, 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,
21391258 221 time, kPNoProcess, nt, 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,
21391258 232 time, kPNoProcess, nt, 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
5a9b7f8e 250 Double_t m, pmax, p, f, phi;
c9a8628a 251 TDatabasePDG * pdg = TDatabasePDG::Instance();
252 const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
253 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
254
255 /* Select nucleon type */
256 if(charge == 0) m = kMassNeutron;
257 else m = kMassProton;
258
259 /* Momentum at maximum of Maxwell-distribution */
260
261 pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
262
263 /* Try until proper momentum */
264 /* for lack of primitive function of the Maxwell-distribution */
265 /* a brute force trial-accept loop, normalized at pmax */
266
267 do
268 {
269 p = Rndm() * fPmax;
270 f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
271 }
272 while(f < Rndm());
273
5a9b7f8e 274 /* Spherical symmetric emission for black particles (beta=0)*/
275 if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
276 /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
277 else if(fThetaDistribution!=0){
278 Double_t costheta = fCosTheta->GetRandom();
279 theta = TMath::ACos(costheta);
280 }
281 //
c9a8628a 282 phi = 2. * TMath::Pi() * Rndm();
283
284 /* Determine momentum components in system of the moving source */
285 q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
286 q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
287 q[2] = p * TMath::Cos(theta);
288
289 /* Transform to system of the target nucleus */
290 /* beta is passed as negative, because the gray nucleons are slowed down */
291 Lorentz(m, -beta, q);
292
293 /* Transform to laboratory system */
294 Lorentz(m, fBeta, q);
83b3cd6c 295 q[2] *= fProtonDirection;
c9a8628a 296}
297
298Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
299{
300/* Relativistic Maxwell-distribution */
301 Double_t ekin;
302 ekin = TMath::Sqrt(p*p+m*m)-m;
303 return (p*p * exp(-ekin/T));
304}
305
306
307void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
308{
309/* Lorentz transform in the direction of q[2] */
310
61cf877f 311 Double_t gamma = 1./sqrt((1.-beta)*(1.+beta));
c9a8628a 312 Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
313 q[2] = gamma * (q[2] + beta*energy);
314}