]>
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> |
6218004b | 33 | #include <TParticle.h> |
c9a8628a | 34 | |
74cfba91 | 35 | #include "AliConst.h" |
c9a8628a | 36 | #include "AliCollisionGeometry.h" |
6218004b | 37 | #include "AliStack.h" |
87f95355 | 38 | #include "AliRun.h" |
39 | #include "AliMC.h" | |
c9a8628a | 40 | #include "AliGenSlowNucleons.h" |
41 | #include "AliSlowNucleonModel.h" | |
42 | ||
706938e6 | 43 | ClassImp(AliGenSlowNucleons) |
1c56e311 | 44 | |
c9a8628a | 45 | |
1c56e311 | 46 | AliGenSlowNucleons::AliGenSlowNucleons() |
47 | :AliGenerator(-1), | |
48 | fCMS(0.), | |
49 | fMomentum(0.), | |
50 | fBeta(0.), | |
51 | fPmax (0.), | |
1c56e311 | 52 | fCharge(0), |
61cf877f | 53 | fProtonDirection(1.), |
1c56e311 | 54 | fTemperatureG(0.), |
55 | fBetaSourceG(0.), | |
56 | fTemperatureB(0.), | |
57 | fBetaSourceB(0.), | |
58 | fNgp(0), | |
59 | fNgn(0), | |
60 | fNbp(0), | |
61 | fNbn(0), | |
62 | fDebug(0), | |
63 | fDebugHist1(0), | |
64 | fDebugHist2(0), | |
65 | fThetaDistribution(), | |
66 | fCosThetaGrayHist(), | |
67 | fCosTheta(), | |
74cfba91 | 68 | fBeamCrossingAngle(0.), |
69 | fBeamDivergence(0.), | |
70 | fBeamDivEvent(0.), | |
1c56e311 | 71 | fSlowNucleonModel(0) |
c9a8628a | 72 | { |
73 | // Default constructor | |
c9a8628a | 74 | fCollisionGeometry = 0; |
75 | } | |
76 | ||
77 | AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart) | |
1c56e311 | 78 | :AliGenerator(npart), |
79 | fCMS(14000.), | |
80 | fMomentum(0.), | |
81 | fBeta(0.), | |
82 | fPmax (10.), | |
1c56e311 | 83 | fCharge(1), |
61cf877f | 84 | fProtonDirection(1.), |
230d85c6 | 85 | fTemperatureG(0.05), |
1c56e311 | 86 | fBetaSourceG(0.05), |
230d85c6 | 87 | fTemperatureB(0.005), |
1c56e311 | 88 | fBetaSourceB(0.), |
89 | fNgp(0), | |
90 | fNgn(0), | |
91 | fNbp(0), | |
92 | fNbn(0), | |
93 | fDebug(0), | |
94 | fDebugHist1(0), | |
95 | fDebugHist2(0), | |
96 | fThetaDistribution(), | |
97 | fCosThetaGrayHist(), | |
98 | fCosTheta(), | |
74cfba91 | 99 | fBeamCrossingAngle(0.), |
100 | fBeamDivergence(0.), | |
101 | fBeamDivEvent(0.), | |
1c56e311 | 102 | fSlowNucleonModel(new AliSlowNucleonModel()) |
74cfba91 | 103 | |
c9a8628a | 104 | { |
105 | // Constructor | |
106 | fName = "SlowNucleons"; | |
107 | fTitle = "Generator for gray particles in pA collisions"; | |
c9a8628a | 108 | fCollisionGeometry = 0; |
c9a8628a | 109 | } |
110 | ||
111 | //____________________________________________________________ | |
112 | AliGenSlowNucleons::~AliGenSlowNucleons() | |
113 | { | |
114 | // Destructor | |
115 | delete fSlowNucleonModel; | |
116 | } | |
117 | ||
83b3cd6c | 118 | void AliGenSlowNucleons::SetProtonDirection(Float_t dir) { |
119 | // Set direction of the proton to change between pA (1) and Ap (-1) | |
120 | fProtonDirection = dir / TMath::Abs(dir); | |
121 | } | |
c9a8628a | 122 | |
123 | void AliGenSlowNucleons::Init() | |
124 | { | |
125 | // | |
126 | // Initialization | |
127 | // | |
88a9f852 | 128 | Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass(); |
271e6aae | 129 | fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget); |
c9a8628a | 130 | fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum); |
88a9f852 | 131 | //printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta); |
c9a8628a | 132 | if (fDebug) { |
c8f7f6f9 | 133 | fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.); |
134 | fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.); | |
5a9b7f8e | 135 | fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.); |
136 | } | |
137 | // | |
138 | // non-uniform cos(theta) distribution | |
139 | // | |
140 | if(fThetaDistribution != 0) { | |
141 | fCosTheta = new TF1("fCosTheta", | |
142 | "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)", | |
143 | -1., 1.); | |
c9a8628a | 144 | } |
74cfba91 | 145 | |
146 | printf("\n AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.); | |
c9a8628a | 147 | } |
148 | ||
149 | void AliGenSlowNucleons::FinishRun() | |
150 | { | |
f687abfe | 151 | // End of run action |
152 | // Show histogram for debugging if requested. | |
c9a8628a | 153 | if (fDebug) { |
c8f7f6f9 | 154 | TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700); |
155 | c->Divide(2,1); | |
156 | c->cd(1); | |
58776b75 | 157 | fDebugHist1->Draw("colz"); |
c8f7f6f9 | 158 | c->cd(2); |
159 | fDebugHist2->Draw(); | |
5a9b7f8e | 160 | c->cd(3); |
161 | fCosThetaGrayHist->Draw(); | |
c9a8628a | 162 | } |
163 | } | |
164 | ||
165 | ||
166 | void AliGenSlowNucleons::Generate() | |
167 | { | |
168 | // | |
169 | // Generate one event | |
170 | // | |
171 | // | |
172 | // Communication with Gray Particle Model | |
173 | // | |
1c22a81e | 174 | if (fCollisionGeometry) { |
175 | Float_t b = fCollisionGeometry->ImpactParameter(); | |
230d85c6 | 176 | // Int_t nn = fCollisionGeometry->NN(); |
177 | // Int_t nwn = fCollisionGeometry->NwN(); | |
178 | // Int_t nnw = fCollisionGeometry->NNw(); | |
179 | // Int_t nwnw = fCollisionGeometry->NwNw(); | |
88a9f852 | 180 | |
58776b75 | 181 | //fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn); |
182 | fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn); | |
1c22a81e | 183 | if (fDebug) { |
88a9f852 | 184 | //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw); |
230d85c6 | 185 | printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn); |
58776b75 | 186 | fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.); |
c8f7f6f9 | 187 | fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.); |
74cfba91 | 188 | |
1c22a81e | 189 | } |
190 | } | |
c9a8628a | 191 | |
192 | // | |
88a9f852 | 193 | Float_t p[3] = {0., 0., 0.}, theta=0; |
c9a8628a | 194 | Float_t origin[3] = {0., 0., 0.}; |
21391258 | 195 | Float_t time = 0.; |
c9a8628a | 196 | Float_t polar [3] = {0., 0., 0.}; |
99f5114f | 197 | Int_t nt, i, j; |
c9a8628a | 198 | Int_t kf; |
74cfba91 | 199 | |
200 | // Extracting 1 value per event for the divergence angle | |
201 | Double_t rvec = gRandom->Gaus(0.0, 1.0); | |
202 | fBeamDivEvent = fBeamDivergence * TMath::Abs(rvec); | |
203 | printf("\n AliGenSlowNucleons: applying beam divergence %f mrad to slow nucleons\n",fBeamDivEvent*1000.); | |
204 | ||
99f5114f | 205 | if(fVertexSmear == kPerEvent) { |
206 | Vertex(); | |
207 | for (j=0; j < 3; j++) origin[j] = fVertex[j]; | |
21391258 | 208 | time = fTime; |
99f5114f | 209 | } // if kPerEvent |
c9a8628a | 210 | // |
211 | // Gray protons | |
212 | // | |
213 | fCharge = 1; | |
214 | kf = kProton; | |
1c22a81e | 215 | for(i = 0; i < fNgp; i++) { |
5a9b7f8e | 216 | GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta); |
217 | if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta)); | |
642f15cf | 218 | PushTrack(fTrackIt, -1, kf, p, origin, polar, |
58776b75 | 219 | time, kPNoProcess, nt, 1.,-2); |
c9a8628a | 220 | KeepTrack(nt); |
87f95355 | 221 | SetProcessID(nt,kGrayProcess); |
c9a8628a | 222 | } |
223 | // | |
224 | // Gray neutrons | |
225 | // | |
226 | fCharge = 0; | |
227 | kf = kNeutron; | |
1c22a81e | 228 | for(i = 0; i < fNgn; i++) { |
5a9b7f8e | 229 | GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta); |
230 | if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta)); | |
642f15cf | 231 | PushTrack(fTrackIt, -1, kf, p, origin, polar, |
58776b75 | 232 | time, kPNoProcess, nt, 1.,-2); |
c9a8628a | 233 | KeepTrack(nt); |
87f95355 | 234 | SetProcessID(nt,kGrayProcess); |
c9a8628a | 235 | } |
236 | // | |
237 | // Black protons | |
238 | // | |
239 | fCharge = 1; | |
240 | kf = kProton; | |
1c22a81e | 241 | for(i = 0; i < fNbp; i++) { |
5a9b7f8e | 242 | GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta); |
642f15cf | 243 | PushTrack(fTrackIt, -1, kf, p, origin, polar, |
58776b75 | 244 | time, kPNoProcess, nt, 1.,-1); |
c9a8628a | 245 | KeepTrack(nt); |
87f95355 | 246 | SetProcessID(nt,kBlackProcess); |
c9a8628a | 247 | } |
248 | // | |
249 | // Black neutrons | |
250 | // | |
251 | fCharge = 0; | |
252 | kf = kNeutron; | |
1c22a81e | 253 | for(i = 0; i < fNbn; i++) { |
5a9b7f8e | 254 | GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta); |
642f15cf | 255 | PushTrack(fTrackIt, -1, kf, p, origin, polar, |
58776b75 | 256 | time, kPNoProcess, nt, 1.,-1); |
c9a8628a | 257 | KeepTrack(nt); |
87f95355 | 258 | SetProcessID(nt,kBlackProcess); |
c9a8628a | 259 | } |
260 | } | |
261 | ||
262 | ||
263 | ||
264 | ||
5a9b7f8e | 265 | void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, |
266 | Double_t beta, Float_t* q, Float_t &theta) | |
f687abfe | 267 | |
268 | { | |
c9a8628a | 269 | /* |
270 | Emit a slow nucleon with "temperature" T [GeV], | |
271 | from a source moving with velocity beta | |
272 | Three-momentum [GeV/c] is given back in q[3] | |
273 | */ | |
274 | ||
88a9f852 | 275 | //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta); |
230d85c6 | 276 | |
5a9b7f8e | 277 | Double_t m, pmax, p, f, phi; |
c9a8628a | 278 | TDatabasePDG * pdg = TDatabasePDG::Instance(); |
279 | const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass(); | |
280 | const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass(); | |
281 | ||
282 | /* Select nucleon type */ | |
283 | if(charge == 0) m = kMassNeutron; | |
284 | else m = kMassProton; | |
285 | ||
286 | /* Momentum at maximum of Maxwell-distribution */ | |
287 | ||
88a9f852 | 288 | pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m))); |
c9a8628a | 289 | |
290 | /* Try until proper momentum */ | |
291 | /* for lack of primitive function of the Maxwell-distribution */ | |
292 | /* a brute force trial-accept loop, normalized at pmax */ | |
293 | ||
294 | do | |
295 | { | |
296 | p = Rndm() * fPmax; | |
297 | f = Maxwell(m, p, T) / Maxwell(m , pmax, T); | |
298 | } | |
299 | while(f < Rndm()); | |
300 | ||
5a9b7f8e | 301 | /* Spherical symmetric emission for black particles (beta=0)*/ |
302 | if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.); | |
303 | /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/ | |
304 | else if(fThetaDistribution!=0){ | |
305 | Double_t costheta = fCosTheta->GetRandom(); | |
306 | theta = TMath::ACos(costheta); | |
307 | } | |
308 | // | |
c9a8628a | 309 | phi = 2. * TMath::Pi() * Rndm(); |
310 | ||
88a9f852 | 311 | |
c9a8628a | 312 | /* Determine momentum components in system of the moving source */ |
313 | q[0] = p * TMath::Sin(theta) * TMath::Cos(phi); | |
314 | q[1] = p * TMath::Sin(theta) * TMath::Sin(phi); | |
315 | q[2] = p * TMath::Cos(theta); | |
74cfba91 | 316 | //if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]); |
c9a8628a | 317 | |
88a9f852 | 318 | |
c9a8628a | 319 | /* Transform to system of the target nucleus */ |
320 | /* beta is passed as negative, because the gray nucleons are slowed down */ | |
321 | Lorentz(m, -beta, q); | |
74cfba91 | 322 | //if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]); |
323 | ||
c9a8628a | 324 | /* Transform to laboratory system */ |
325 | Lorentz(m, fBeta, q); | |
83b3cd6c | 326 | q[2] *= fProtonDirection; |
74cfba91 | 327 | if(fDebug==1)printf("\n Momentum after LHC boost: p = (%f, %f, %f)\n",q[0],q[1],q[2]); |
328 | ||
329 | if(fBeamCrossingAngle>0.) BeamCrossDivergence(1, q); // applying crossing angle | |
330 | if(fBeamDivergence>0.) BeamCrossDivergence(2, q); // applying divergence | |
331 | ||
c9a8628a | 332 | } |
333 | ||
334 | Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T) | |
335 | { | |
336 | /* Relativistic Maxwell-distribution */ | |
337 | Double_t ekin; | |
338 | ekin = TMath::Sqrt(p*p+m*m)-m; | |
339 | return (p*p * exp(-ekin/T)); | |
340 | } | |
341 | ||
342 | ||
74cfba91 | 343 | //_____________________________________________________________________________ |
c9a8628a | 344 | void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q) |
345 | { | |
346 | /* Lorentz transform in the direction of q[2] */ | |
347 | ||
88a9f852 | 348 | Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta)); |
349 | Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]); | |
c9a8628a | 350 | q[2] = gamma * (q[2] + beta*energy); |
88a9f852 | 351 | //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]); |
c9a8628a | 352 | } |
230d85c6 | 353 | |
74cfba91 | 354 | //_____________________________________________________________________________ |
355 | void AliGenSlowNucleons::BeamCrossDivergence(Int_t iwhat, Float_t *pLab) | |
356 | { | |
357 | // Applying beam divergence and crossing angle | |
358 | // | |
359 | Double_t pmod = TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]+pLab[2]*pLab[2]); | |
360 | ||
361 | Double_t tetdiv = 0.; | |
362 | Double_t fidiv = 0.; | |
363 | if(iwhat==1){ | |
364 | tetdiv = fBeamCrossingAngle; | |
365 | fidiv = k2PI/4.; | |
366 | } | |
367 | else if(iwhat==2){ | |
368 | tetdiv = fBeamDivEvent; | |
369 | fidiv = (gRandom->Rndm())*k2PI; | |
370 | } | |
371 | ||
372 | Double_t tetpart = TMath::ATan2(TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]), pLab[2]); | |
373 | Double_t fipart=0.; | |
374 | if(TMath::Abs(pLab[1])>0. || TMath::Abs(pLab[0])>0.) fipart = TMath::ATan2(pLab[1], pLab[0]); | |
375 | if(fipart<0.) {fipart = fipart+k2PI;} | |
376 | tetdiv = tetdiv*kRaddeg; | |
377 | fidiv = fidiv*kRaddeg; | |
378 | tetpart = tetpart*kRaddeg; | |
379 | fipart = fipart*kRaddeg; | |
380 | ||
381 | Double_t angleSum[2]={0., 0.}; | |
382 | AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum); | |
383 | ||
384 | Double_t tetsum = angleSum[0]; | |
385 | Double_t fisum = angleSum[1]; | |
386 | //printf("tetpart %f fipart %f tetdiv %f fidiv %f angleSum %f %f\n",tetpart,fipart,tetdiv,fidiv,angleSum[0],angleSum[1]); | |
387 | tetsum = tetsum*kDegrad; | |
388 | fisum = fisum*kDegrad; | |
389 | ||
390 | pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum); | |
391 | pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum); | |
392 | pLab[2] = pmod*TMath::Cos(tetsum); | |
393 | if(fDebug==1){ | |
394 | if(iwhat==1) printf(" Beam crossing angle %f mrad ", fBeamCrossingAngle*1000.); | |
395 | else if(iwhat==2) printf(" Beam divergence %f mrad ", fBeamDivEvent*1000.); | |
396 | printf(" p = (%f, %f, %f)\n",pLab[0],pLab[1],pLab[2]); | |
397 | } | |
398 | } | |
399 | ||
400 | //_____________________________________________________________________________ | |
401 | void AliGenSlowNucleons::AddAngle(Double_t theta1, Double_t phi1, | |
402 | Double_t theta2, Double_t phi2, Double_t *angleSum) | |
403 | { | |
404 | // Calculating the sum of 2 angles | |
405 | Double_t temp = -1.; | |
406 | Double_t conv = 180./TMath::ACos(temp); | |
407 | ||
408 | Double_t ct1 = TMath::Cos(theta1/conv); | |
409 | Double_t st1 = TMath::Sin(theta1/conv); | |
410 | Double_t cp1 = TMath::Cos(phi1/conv); | |
411 | Double_t sp1 = TMath::Sin(phi1/conv); | |
412 | Double_t ct2 = TMath::Cos(theta2/conv); | |
413 | Double_t st2 = TMath::Sin(theta2/conv); | |
414 | Double_t cp2 = TMath::Cos(phi2/conv); | |
415 | Double_t sp2 = TMath::Sin(phi2/conv); | |
416 | Double_t cx = ct1*cp1*st2*cp2+st1*cp1*ct2-sp1*st2*sp2; | |
417 | Double_t cy = ct1*sp1*st2*cp2+st1*sp1*ct2+cp1*st2*sp2; | |
418 | Double_t cz = ct1*ct2-st1*st2*cp2; | |
419 | ||
420 | Double_t rtetsum = TMath::ACos(cz); | |
421 | Double_t tetsum = conv*rtetsum; | |
422 | if(TMath::Abs(tetsum)<1e-4 || tetsum==180.) return; | |
423 | ||
424 | temp = cx/TMath::Sin(rtetsum); | |
425 | if(temp>1.) temp=1.; | |
426 | if(temp<-1.) temp=-1.; | |
427 | Double_t fisum = conv*TMath::ACos(temp); | |
428 | if(cy<0) {fisum = 360.-fisum;} | |
429 | ||
430 | angleSum[0] = tetsum; | |
431 | angleSum[1] = fisum; | |
432 | } | |
433 | ||
87f95355 | 434 | //_____________________________________________________________________________ |
435 | void AliGenSlowNucleons::SetProcessID(Int_t nt, UInt_t process) | |
436 | { | |
437 | // Tag the particle as | |
438 | // gray or black | |
439 | if (fStack) | |
440 | fStack->Particle(nt)->SetUniqueID(process); | |
441 | else | |
442 | gAlice->GetMCApp()->Particle(nt)->SetUniqueID(process); | |
443 | } |