Introducing the interaction time into the aliroot generators. In case of gaussian...
[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      fATarget (0.),
48      fZTarget (0.),
49      fCharge(0),
50      fProtonDirection(0.),
51      fTemperatureG(0.), 
52      fBetaSourceG(0.),
53      fTemperatureB(0.),
54      fBetaSourceB(0.),
55      fNgp(0),
56      fNgn(0),
57      fNbp(0),
58      fNbn(0),
59      fDebug(0),
60      fDebugHist1(0),
61      fDebugHist2(0),
62      fThetaDistribution(),
63      fCosThetaGrayHist(),
64      fCosTheta(),
65      fSlowNucleonModel(0)
66 {
67 // Default constructor
68     fCollisionGeometry = 0;
69 }
70
71 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
72     :AliGenerator(npart),
73      fCMS(14000.),
74      fMomentum(0.),
75      fBeta(0.),
76      fPmax (10.),
77      fATarget (208.),
78      fZTarget (82.),
79      fCharge(1),
80      fProtonDirection(0.),
81      fTemperatureG(0.04), 
82      fBetaSourceG(0.05),
83      fTemperatureB(0.004),
84      fBetaSourceB(0.),
85      fNgp(0),
86      fNgn(0),
87      fNbp(0),
88      fNbn(0),
89      fDebug(0),
90      fDebugHist1(0),
91      fDebugHist2(0),
92      fThetaDistribution(),
93      fCosThetaGrayHist(),
94      fCosTheta(),
95      fSlowNucleonModel(new AliSlowNucleonModel())
96 {
97 // Constructor
98     fName  = "SlowNucleons";
99     fTitle = "Generator for gray particles in pA collisions";
100     fCollisionGeometry = 0;
101 }
102
103 //____________________________________________________________
104 AliGenSlowNucleons::~AliGenSlowNucleons()
105 {
106 // Destructor
107     delete  fSlowNucleonModel;
108 }
109
110 void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
111 // Set direction of the proton to change between pA (1) and Ap (-1)
112   fProtonDirection = dir / TMath::Abs(dir);
113 }
114
115 void AliGenSlowNucleons::Init()
116 {
117   //
118   // Initialization
119   //
120     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
121     fMomentum = fCMS/2. * fZTarget / fATarget;
122     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
123     if (fDebug) {
124         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
125         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
126         fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
127     }
128     //
129     // non-uniform cos(theta) distribution
130     //
131     if(fThetaDistribution != 0) {
132         fCosTheta = new TF1("fCosTheta",
133                             "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
134                             -1., 1.);
135     }
136 }
137
138 void AliGenSlowNucleons::FinishRun()
139 {
140 // End of run action
141 // Show histogram for debugging if requested.
142     if (fDebug) {
143         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
144         c->Divide(2,1);
145         c->cd(1);
146         fDebugHist1->Draw();
147         c->cd(2);
148         fDebugHist2->Draw();
149         c->cd(3);
150         fCosThetaGrayHist->Draw();
151     }
152 }
153
154
155 void AliGenSlowNucleons::Generate()
156 {
157   //
158   // Generate one event
159   //
160   //
161   // Communication with Gray Particle Model 
162   // 
163     if (fCollisionGeometry) {
164         Float_t b   = fCollisionGeometry->ImpactParameter();
165         Int_t  nn   = fCollisionGeometry->NN();
166         Int_t  nwn  = fCollisionGeometry->NwN();
167         Int_t  nnw  = fCollisionGeometry->NNw();
168         Int_t  nwnw = fCollisionGeometry->NwNw();
169
170         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
171         if (fDebug) {
172             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
173             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
174             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
175             printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", 
176                    b, nn, nwn, nnw, nwnw);
177         }
178     }     
179
180    //
181     Float_t p[3], theta=0;
182     Float_t origin[3] = {0., 0., 0.};
183     Float_t time = 0.;
184     Float_t polar [3] = {0., 0., 0.};    
185     Int_t nt, i, j;
186     Int_t kf;
187     
188     if(fVertexSmear == kPerEvent) {
189         Vertex();
190         for (j=0; j < 3; j++) origin[j] = fVertex[j];
191         time = fTime;
192     } // if kPerEvent
193 //
194 //  Gray protons
195 //
196     fCharge = 1;
197     kf = kProton;    
198     for(i = 0; i < fNgp; i++) {
199         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
200         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
201         PushTrack(fTrackIt, -1, kf, p, origin, polar,
202                  time, kPNoProcess, nt, 1.);
203         KeepTrack(nt);
204     }
205 //
206 //  Gray neutrons
207 //
208     fCharge = 0;
209     kf = kNeutron;    
210     for(i = 0; i < fNgn; i++) {
211         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
212         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
213         PushTrack(fTrackIt, -1, kf, p, origin, polar,
214                  time, kPNoProcess, nt, 1.);
215         KeepTrack(nt);
216     }
217 //
218 //  Black protons
219 //
220     fCharge = 1;
221     kf = kProton;    
222     for(i = 0; i < fNbp; i++) {
223         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
224         PushTrack(fTrackIt, -1, kf, p, origin, polar,
225                  time, kPNoProcess, nt, 1.);
226         KeepTrack(nt);
227     }
228 //
229 //  Black neutrons
230 //
231     fCharge = 0;
232     kf = kNeutron;    
233     for(i = 0; i < fNbn; i++) {
234         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
235         PushTrack(fTrackIt, -1, kf, p, origin, polar,
236                  time, kPNoProcess, nt, 1.);
237         KeepTrack(nt);
238     }
239 }
240
241
242
243
244 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, 
245         Double_t beta, Float_t* q, Float_t &theta)
246
247 {
248 /* 
249    Emit a slow nucleon with "temperature" T [GeV], 
250    from a source moving with velocity beta         
251    Three-momentum [GeV/c] is given back in q[3]    
252 */
253
254  Double_t m, pmax, p, f, phi;
255  TDatabasePDG * pdg = TDatabasePDG::Instance();
256  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
257  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
258  
259  /* Select nucleon type */
260  if(charge == 0) m = kMassNeutron;
261  else m = kMassProton;
262
263  /* Momentum at maximum of Maxwell-distribution */
264
265  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
266
267  /* Try until proper momentum                                  */
268  /* for lack of primitive function of the Maxwell-distribution */
269  /* a brute force trial-accept loop, normalized at pmax        */
270
271  do
272  {
273      p = Rndm() * fPmax;
274      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
275  }
276  while(f < Rndm());
277
278  /* Spherical symmetric emission for black particles (beta=0)*/
279  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
280  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
281  else if(fThetaDistribution!=0){
282    Double_t costheta = fCosTheta->GetRandom();
283    theta = TMath::ACos(costheta);
284  }
285  //
286  phi   = 2. * TMath::Pi() * Rndm();
287
288  /* Determine momentum components in system of the moving source */
289  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
290  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
291  q[2] = p * TMath::Cos(theta);
292
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);
296
297  /* Transform to laboratory system */
298  Lorentz(m, fBeta, q);
299  q[2] *= fProtonDirection; 
300 }
301
302 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
303 {
304 /* Relativistic Maxwell-distribution */
305     Double_t ekin;
306     ekin = TMath::Sqrt(p*p+m*m)-m;
307     return (p*p * exp(-ekin/T));
308 }
309
310
311 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
312 {
313 /* Lorentz transform in the direction of q[2] */
314  
315     Double_t gamma  = 1/sqrt((1.-beta)*(1.+beta)); 
316     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
317     q[2] = gamma * (q[2] + beta*energy);
318 }