changes from Misha
[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 /* $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
34 #include "AliCollisionGeometry.h"
35 #include "AliGenSlowNucleons.h"
36 #include "AliSlowNucleonModel.h"
37
38 ClassImp(AliGenSlowNucleons)
39
40     
41 AliGenSlowNucleons::AliGenSlowNucleons()
42     :AliGenerator(-1),
43      fCMS(0.),
44      fMomentum(0.),
45      fBeta(0.),
46      fPmax (0.),
47      fCharge(0),
48      fProtonDirection(1.),
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)
64 {
65 // Default constructor
66     fCollisionGeometry = 0;
67 }
68
69 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
70     :AliGenerator(npart),
71      fCMS(14000.),
72      fMomentum(0.),
73      fBeta(0.),
74      fPmax (10.),
75      fCharge(1),
76      fProtonDirection(1.),
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())
92 {
93 // Constructor
94     fName  = "SlowNucleons";
95     fTitle = "Generator for gray particles in pA collisions";
96     fCollisionGeometry = 0;
97 }
98
99 //____________________________________________________________
100 AliGenSlowNucleons::~AliGenSlowNucleons()
101 {
102 // Destructor
103     delete  fSlowNucleonModel;
104 }
105
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 }
110
111 void AliGenSlowNucleons::Init()
112 {
113   //
114   // Initialization
115   //
116     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
117     fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
118     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
119     if (fDebug) {
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.);
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.);
131     }
132 }
133
134 void AliGenSlowNucleons::FinishRun()
135 {
136 // End of run action
137 // Show histogram for debugging if requested.
138     if (fDebug) {
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();
145         c->cd(3);
146         fCosThetaGrayHist->Draw();
147     }
148 }
149
150
151 void AliGenSlowNucleons::Generate()
152 {
153   //
154   // Generate one event
155   //
156   //
157   // Communication with Gray Particle Model 
158   // 
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
166         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
167         if (fDebug) {
168             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
169             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
170             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
171             printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", 
172                    b, nn, nwn, nnw, nwnw);
173         }
174     }     
175
176    //
177     Float_t p[3], theta=0;
178     Float_t origin[3] = {0., 0., 0.};
179     Float_t time = 0.;
180     Float_t polar [3] = {0., 0., 0.};    
181     Int_t nt, i, j;
182     Int_t kf;
183     
184     if(fVertexSmear == kPerEvent) {
185         Vertex();
186         for (j=0; j < 3; j++) origin[j] = fVertex[j];
187         time = fTime;
188     } // if kPerEvent
189 //
190 //  Gray protons
191 //
192     fCharge = 1;
193     kf = kProton;    
194     for(i = 0; i < fNgp; i++) {
195         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
196         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
197         PushTrack(fTrackIt, -1, kf, p, origin, polar,
198                  time, kPNoProcess, nt, 1.);
199         KeepTrack(nt);
200     }
201 //
202 //  Gray neutrons
203 //
204     fCharge = 0;
205     kf = kNeutron;    
206     for(i = 0; i < fNgn; i++) {
207         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
208         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
209         PushTrack(fTrackIt, -1, kf, p, origin, polar,
210                  time, kPNoProcess, nt, 1.);
211         KeepTrack(nt);
212     }
213 //
214 //  Black protons
215 //
216     fCharge = 1;
217     kf = kProton;    
218     for(i = 0; i < fNbp; i++) {
219         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
220         PushTrack(fTrackIt, -1, kf, p, origin, polar,
221                  time, kPNoProcess, nt, 1.);
222         KeepTrack(nt);
223     }
224 //
225 //  Black neutrons
226 //
227     fCharge = 0;
228     kf = kNeutron;    
229     for(i = 0; i < fNbn; i++) {
230         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
231         PushTrack(fTrackIt, -1, kf, p, origin, polar,
232                  time, kPNoProcess, nt, 1.);
233         KeepTrack(nt);
234     }
235 }
236
237
238
239
240 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, 
241         Double_t beta, Float_t* q, Float_t &theta)
242
243 {
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
250  Double_t m, pmax, p, f, phi;
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
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  //
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);
295  q[2] *= fProtonDirection; 
296 }
297
298 Double_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
307 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
308 {
309 /* Lorentz transform in the direction of q[2] */
310  
311     Double_t gamma  = 1./sqrt((1.-beta)*(1.+beta)); 
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 }