New reader for the pedestal run and vdrift (Julian) and some bug fixing (Raphaelle)
[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 polar [3] = {0., 0., 0.};    
184     Int_t nt, i, j;
185     Int_t kf;
186     
187     if(fVertexSmear == kPerEvent) {
188         Vertex();
189         for (j=0; j < 3; j++) origin[j] = fVertex[j];
190     } // if kPerEvent
191 //
192 //  Gray protons
193 //
194     fCharge = 1;
195     kf = kProton;    
196     for(i = 0; i < fNgp; i++) {
197         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
198         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
199         PushTrack(fTrackIt, -1, kf, p, origin, polar,
200                  0., kPNoProcess, nt, 1.);
201         KeepTrack(nt);
202     }
203 //
204 //  Gray neutrons
205 //
206     fCharge = 0;
207     kf = kNeutron;    
208     for(i = 0; i < fNgn; i++) {
209         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
210         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
211         PushTrack(fTrackIt, -1, kf, p, origin, polar,
212                  0., kPNoProcess, nt, 1.);
213         KeepTrack(nt);
214     }
215 //
216 //  Black protons
217 //
218     fCharge = 1;
219     kf = kProton;    
220     for(i = 0; i < fNbp; i++) {
221         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
222         PushTrack(fTrackIt, -1, kf, p, origin, polar,
223                  0., kPNoProcess, nt, 1.);
224         KeepTrack(nt);
225     }
226 //
227 //  Black neutrons
228 //
229     fCharge = 0;
230     kf = kNeutron;    
231     for(i = 0; i < fNbn; i++) {
232         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
233         PushTrack(fTrackIt, -1, kf, p, origin, polar,
234                  0., kPNoProcess, nt, 1.);
235         KeepTrack(nt);
236     }
237 }
238
239
240
241
242 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, 
243         Double_t beta, Float_t* q, Float_t &theta)
244
245 {
246 /* 
247    Emit a slow nucleon with "temperature" T [GeV], 
248    from a source moving with velocity beta         
249    Three-momentum [GeV/c] is given back in q[3]    
250 */
251
252  Double_t m, pmax, p, f, phi;
253  TDatabasePDG * pdg = TDatabasePDG::Instance();
254  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
255  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
256  
257  /* Select nucleon type */
258  if(charge == 0) m = kMassNeutron;
259  else m = kMassProton;
260
261  /* Momentum at maximum of Maxwell-distribution */
262
263  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
264
265  /* Try until proper momentum                                  */
266  /* for lack of primitive function of the Maxwell-distribution */
267  /* a brute force trial-accept loop, normalized at pmax        */
268
269  do
270  {
271      p = Rndm() * fPmax;
272      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
273  }
274  while(f < Rndm());
275
276  /* Spherical symmetric emission for black particles (beta=0)*/
277  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
278  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
279  else if(fThetaDistribution!=0){
280    Double_t costheta = fCosTheta->GetRandom();
281    theta = TMath::ACos(costheta);
282  }
283  //
284  phi   = 2. * TMath::Pi() * Rndm();
285
286  /* Determine momentum components in system of the moving source */
287  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
288  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
289  q[2] = p * TMath::Cos(theta);
290
291  /* Transform to system of the target nucleus                             */
292  /* beta is passed as negative, because the gray nucleons are slowed down */
293  Lorentz(m, -beta, q);
294
295  /* Transform to laboratory system */
296  Lorentz(m, fBeta, q);
297  q[2] *= fProtonDirection; 
298 }
299
300 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
301 {
302 /* Relativistic Maxwell-distribution */
303     Double_t ekin;
304     ekin = TMath::Sqrt(p*p+m*m)-m;
305     return (p*p * exp(-ekin/T));
306 }
307
308
309 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
310 {
311 /* Lorentz transform in the direction of q[2] */
312  
313     Double_t gamma  = 1/sqrt((1.-beta)*(1.+beta)); 
314     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
315     q[2] = gamma * (q[2] + beta*energy);
316 }