]>
Commit | Line | Data |
---|---|---|
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 | 38 | ClassImp(AliGenSlowNucleons) |
1c56e311 | 39 | |
c9a8628a | 40 | |
1c56e311 | 41 | AliGenSlowNucleons::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 | ||
69 | AliGenSlowNucleons::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 | ||
99 | //____________________________________________________________ | |
100 | AliGenSlowNucleons::~AliGenSlowNucleons() | |
101 | { | |
102 | // Destructor | |
103 | delete fSlowNucleonModel; | |
104 | } | |
105 | ||
83b3cd6c | 106 | void 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 | |
111 | void 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 | ||
135 | void 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 | ||
152 | void 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, |
88a9f852 | 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, |
88a9f852 | 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, |
88a9f852 | 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, |
88a9f852 | 232 | time, kPNoProcess, nt, 1.,-1); |
c9a8628a | 233 | KeepTrack(nt); |
234 | } | |
235 | } | |
236 | ||
237 | ||
238 | ||
239 | ||
5a9b7f8e | 240 | void 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 | ||
303 | Double_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 | ||
312 | void 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 |