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