]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenSlowNucleons.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Correction
[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx
... / ...
CommitLineData
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
16/* $Id$ */
17
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//
26
27#include <TDatabasePDG.h>
28#include <TPDGCode.h>
29#include <TH2F.h>
30#include <TH1F.h>
31#include <TF1.h>
32#include <TCanvas.h>
33#include <TParticle.h>
34
35#include "AliConst.h"
36#include "AliCollisionGeometry.h"
37#include "AliStack.h"
38#include "AliRun.h"
39#include "AliMC.h"
40#include "AliGenSlowNucleons.h"
41#include "AliSlowNucleonModel.h"
42
43ClassImp(AliGenSlowNucleons)
44
45
46AliGenSlowNucleons::AliGenSlowNucleons()
47 :AliGenerator(-1),
48 fCMS(0.),
49 fMomentum(0.),
50 fBeta(0.),
51 fPmax (0.),
52 fCharge(0),
53 fProtonDirection(1.),
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(),
68 fBeamCrossingAngle(0.),
69 fBeamDivergence(0.),
70 fBeamDivEvent(0.),
71 fSlowNucleonModel(0)
72{
73// Default constructor
74 fCollisionGeometry = 0;
75}
76
77AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
78 :AliGenerator(npart),
79 fCMS(14000.),
80 fMomentum(0.),
81 fBeta(0.),
82 fPmax (10.),
83 fCharge(1),
84 fProtonDirection(1.),
85 fTemperatureG(0.05),
86 fBetaSourceG(0.05),
87 fTemperatureB(0.005),
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(),
99 fBeamCrossingAngle(0.),
100 fBeamDivergence(0.),
101 fBeamDivEvent(0.),
102 fSlowNucleonModel(new AliSlowNucleonModel())
103
104{
105// Constructor
106 fName = "SlowNucleons";
107 fTitle = "Generator for gray particles in pA collisions";
108 fCollisionGeometry = 0;
109}
110
111//____________________________________________________________
112AliGenSlowNucleons::~AliGenSlowNucleons()
113{
114// Destructor
115 delete fSlowNucleonModel;
116}
117
118void 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}
122
123void AliGenSlowNucleons::Init()
124{
125 //
126 // Initialization
127 //
128 Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
129 fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
130 fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
131 //printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
132 if (fDebug) {
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.);
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.);
144 }
145
146 printf("\n AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
147}
148
149void AliGenSlowNucleons::FinishRun()
150{
151// End of run action
152// Show histogram for debugging if requested.
153 if (fDebug) {
154 TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
155 c->Divide(2,1);
156 c->cd(1);
157 fDebugHist1->Draw("colz");
158 c->cd(2);
159 fDebugHist2->Draw();
160 c->cd(3);
161 fCosThetaGrayHist->Draw();
162 }
163}
164
165
166void AliGenSlowNucleons::Generate()
167{
168 //
169 // Generate one event
170 //
171 //
172 // Communication with Gray Particle Model
173 //
174 if (fCollisionGeometry) {
175 Float_t b = fCollisionGeometry->ImpactParameter();
176 // Int_t nn = fCollisionGeometry->NN();
177 // Int_t nwn = fCollisionGeometry->NwN();
178 // Int_t nnw = fCollisionGeometry->NNw();
179 // Int_t nwnw = fCollisionGeometry->NwNw();
180
181 //fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
182 fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
183 if (fDebug) {
184 //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
185 printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
186 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
187 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
188
189 }
190 }
191
192 //
193 Float_t p[3] = {0., 0., 0.}, theta=0;
194 Float_t origin[3] = {0., 0., 0.};
195 Float_t time = 0.;
196 Float_t polar [3] = {0., 0., 0.};
197 Int_t nt, i, j;
198 Int_t kf;
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
205 if(fVertexSmear == kPerEvent) {
206 Vertex();
207 for (j=0; j < 3; j++) origin[j] = fVertex[j];
208 time = fTime;
209 } // if kPerEvent
210//
211// Gray protons
212//
213 fCharge = 1;
214 kf = kProton;
215 for(i = 0; i < fNgp; i++) {
216 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
217 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
218 PushTrack(fTrackIt, -1, kf, p, origin, polar,
219 time, kPNoProcess, nt, 1.,-2);
220 KeepTrack(nt);
221 SetProcessID(nt,kGrayProcess);
222 }
223//
224// Gray neutrons
225//
226 fCharge = 0;
227 kf = kNeutron;
228 for(i = 0; i < fNgn; i++) {
229 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
230 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
231 PushTrack(fTrackIt, -1, kf, p, origin, polar,
232 time, kPNoProcess, nt, 1.,-2);
233 KeepTrack(nt);
234 SetProcessID(nt,kGrayProcess);
235 }
236//
237// Black protons
238//
239 fCharge = 1;
240 kf = kProton;
241 for(i = 0; i < fNbp; i++) {
242 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
243 PushTrack(fTrackIt, -1, kf, p, origin, polar,
244 time, kPNoProcess, nt, 1.,-1);
245 KeepTrack(nt);
246 SetProcessID(nt,kBlackProcess);
247 }
248//
249// Black neutrons
250//
251 fCharge = 0;
252 kf = kNeutron;
253 for(i = 0; i < fNbn; i++) {
254 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
255 PushTrack(fTrackIt, -1, kf, p, origin, polar,
256 time, kPNoProcess, nt, 1.,-1);
257 KeepTrack(nt);
258 SetProcessID(nt,kBlackProcess);
259 }
260}
261
262
263
264
265void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
266 Double_t beta, Float_t* q, Float_t &theta)
267
268{
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
275 //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
276
277 Double_t m, pmax, p, f, phi;
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
288 pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
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
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 //
309 phi = 2. * TMath::Pi() * Rndm();
310
311
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);
316 //if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
317
318
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);
322 //if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
323
324 /* Transform to laboratory system */
325 Lorentz(m, fBeta, q);
326 q[2] *= fProtonDirection;
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
332}
333
334Double_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
343//_____________________________________________________________________________
344void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
345{
346/* Lorentz transform in the direction of q[2] */
347
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]);
350 q[2] = gamma * (q[2] + beta*energy);
351 //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
352}
353
354//_____________________________________________________________________________
355void 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//_____________________________________________________________________________
401void 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
434//_____________________________________________________________________________
435void 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}