137746dc677901bd0353cc27144911dd06c4ead2
[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx
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 /*
17 $Log$
18 Revision 1.1  2003/03/24 16:49:23  morsch
19 Slow nucleon generator and model.
20
21 */
22
23 /*
24   Generator for slow nucluons in pA interactions. 
25   Source is modelled by a relativistic Maxwell distributions.
26   Original code by  Ferenc Sikler  <sikler@rmki.kfki.hu>
27  */
28
29 #include <TDatabasePDG.h>
30 #include <TPDGCode.h>
31 #include <TH2F.h>
32
33 #include "AliCollisionGeometry.h"
34 #include "AliGenSlowNucleons.h"
35 #include "AliSlowNucleonModel.h"
36
37  ClassImp(AliGenSlowNucleons)
38     
39  AliGenSlowNucleons::AliGenSlowNucleons():AliGenerator(-1)
40 {
41 // Default constructor
42     fSlowNucleonModel = 0;
43     fCollisionGeometry = 0;
44 }
45
46 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
47     :AliGenerator(npart)
48 {
49 // Constructor
50     fName  = "SlowNucleons";
51     fTitle = "Generator for gray particles in pA collisions";
52     SetPmax();
53     SetTarget();
54     SetNominalCmsEnergy();
55     SetCharge();
56     SetTemperature();
57     SetBetaSource();
58     fSlowNucleonModel = new AliSlowNucleonModel();
59     fCollisionGeometry = 0;
60     fDebug = 0;
61 }
62
63 //____________________________________________________________
64 AliGenSlowNucleons::~AliGenSlowNucleons()
65 {
66 // Destructor
67     delete  fSlowNucleonModel;
68 }
69
70
71 void AliGenSlowNucleons::Init()
72 {
73   //
74   // Initialization
75   //
76     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
77     fMomentum = fCMS/2. * fZTarget / fATarget;
78     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
79     if (fDebug) {
80         fDebugHist = new TH2F("DebugHist", "N_heavy vs nu", 50, 0., 50., 50, 0., 50.);
81     }
82 }
83
84 void AliGenSlowNucleons::FinishRun()
85 {
86     if (fDebug) {
87         fDebugHist->Draw();
88     }
89 }
90
91
92 void AliGenSlowNucleons::Generate()
93 {
94   //
95   // Generate one event
96   //
97   //
98   // Communication with Gray Particle Model 
99   // 
100     Int_t ngp, ngn, nbp, nbn;
101
102     if (fCollisionGeometry) {
103         Float_t b   = fCollisionGeometry->ImpactParameter();
104         Int_t  nn   = fCollisionGeometry->NN();
105         Int_t  nwn  = fCollisionGeometry->NwN();
106         Int_t  nnw  = fCollisionGeometry->NNw();
107         Int_t  nwnw = fCollisionGeometry->NwNw();
108
109         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, ngp, ngn, nbp, nbn);
110         if (fDebug) {
111             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
112             fDebugHist->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
113             printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", 
114                    b, nn, nwn, nnw, nwnw);
115         }
116     }     
117
118    //
119     Float_t p[3];
120     Float_t origin[3] = {0., 0., 0.};
121     Float_t polar [3] = {0., 0., 0.};    
122     Int_t nt, i;
123     Int_t kf;
124 //
125 //  Gray protons
126 //
127     fCharge = 1;
128     kf = kProton;    
129     for(i = 0; i < fNgp; i++) {
130         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
131         SetTrack(fTrackIt, -1, kf, p, origin, polar,
132                  0., kPNoProcess, nt, 1.);
133         KeepTrack(nt);
134     }
135 //
136 //  Gray neutrons
137 //
138     fCharge = 0;
139     kf = kNeutron;    
140     for(i = 0; i < fNgn; i++) {
141         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
142         SetTrack(fTrackIt, -1, kf, p, origin, polar,
143                  0., kPNoProcess, nt, 1.);
144         KeepTrack(nt);
145     }
146 //
147 //  Black protons
148 //
149     fCharge = 1;
150     kf = kProton;    
151     for(i = 0; i < fNbp; i++) {
152         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
153         SetTrack(fTrackIt, -1, kf, p, origin, polar,
154                  0., kPNoProcess, nt, 1.);
155         KeepTrack(nt);
156     }
157 //
158 //  Black neutrons
159 //
160     fCharge = 0;
161     kf = kNeutron;    
162     for(i = 0; i < fNbn; i++) {
163         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
164         SetTrack(fTrackIt, -1, kf, p, origin, polar,
165                  0., kPNoProcess, nt, 1.);
166         KeepTrack(nt);
167     }
168 }
169
170
171
172
173 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, Double_t beta, Float_t* q)
174 /* 
175    Emit a slow nucleon with "temperature" T [GeV], 
176    from a source moving with velocity beta         
177    Three-momentum [GeV/c] is given back in q[3]    
178 */
179
180 {
181  Double_t m, pmax, p, f, theta, phi;
182  
183  TDatabasePDG * pdg = TDatabasePDG::Instance();
184  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
185  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
186  
187  /* Select nucleon type */
188  if(charge == 0) m = kMassNeutron;
189  else m = kMassProton;
190
191  /* Momentum at maximum of Maxwell-distribution */
192
193  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
194
195  /* Try until proper momentum                                  */
196  /* for lack of primitive function of the Maxwell-distribution */
197  /* a brute force trial-accept loop, normalized at pmax        */
198
199  do
200  {
201      p = Rndm() * fPmax;
202      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
203  }
204  while(f < Rndm());
205
206  /* Spherical symmetric emission */
207  theta = TMath::ACos(2. * Rndm() - 1.);
208  phi   = 2. * TMath::Pi() * Rndm();
209
210  /* Determine momentum components in system of the moving source */
211  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
212  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
213  q[2] = p * TMath::Cos(theta);
214
215  /* Transform to system of the target nucleus                             */
216  /* beta is passed as negative, because the gray nucleons are slowed down */
217  Lorentz(m, -beta, q);
218
219  /* Transform to laboratory system */
220  Lorentz(m, fBeta, q);
221 }
222
223 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
224 {
225 /* Relativistic Maxwell-distribution */
226     Double_t ekin;
227     ekin = TMath::Sqrt(p*p+m*m)-m;
228     return (p*p * exp(-ekin/T));
229 }
230
231
232 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
233 {
234 /* Lorentz transform in the direction of q[2] */
235  
236     Double_t gamma  = 1/sqrt(1-beta*beta); 
237     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
238     q[2] = gamma * (q[2] + beta*energy);
239 }
240
241
242
243
244
245
246