new version of trigger patch-cluster matching more detailed histograms
[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx
CommitLineData
c9a8628a 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
88cb7938 16/* $Id$ */
c9a8628a 17
f687abfe 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//
c9a8628a 26
27#include <TDatabasePDG.h>
28#include <TPDGCode.h>
29#include <TH2F.h>
5a9b7f8e 30#include <TH1F.h>
31#include <TF1.h>
c8f7f6f9 32#include <TCanvas.h>
c9a8628a 33
34#include "AliCollisionGeometry.h"
35#include "AliGenSlowNucleons.h"
36#include "AliSlowNucleonModel.h"
37
706938e6 38ClassImp(AliGenSlowNucleons)
1c56e311 39
c9a8628a 40
1c56e311 41AliGenSlowNucleons::AliGenSlowNucleons()
42 :AliGenerator(-1),
43 fCMS(0.),
44 fMomentum(0.),
45 fBeta(0.),
46 fPmax (0.),
1c56e311 47 fCharge(0),
61cf877f 48 fProtonDirection(1.),
1c56e311 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)
c9a8628a 64{
65// Default constructor
c9a8628a 66 fCollisionGeometry = 0;
67}
68
69AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
1c56e311 70 :AliGenerator(npart),
71 fCMS(14000.),
72 fMomentum(0.),
73 fBeta(0.),
74 fPmax (10.),
1c56e311 75 fCharge(1),
61cf877f 76 fProtonDirection(1.),
230d85c6 77 fTemperatureG(0.05),
1c56e311 78 fBetaSourceG(0.05),
230d85c6 79 fTemperatureB(0.005),
1c56e311 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())
c9a8628a 92{
93// Constructor
94 fName = "SlowNucleons";
95 fTitle = "Generator for gray particles in pA collisions";
c9a8628a 96 fCollisionGeometry = 0;
c9a8628a 97}
98
5a9b7f8e 99//____________________________________________________________
c9a8628a 100AliGenSlowNucleons::~AliGenSlowNucleons()
101{
102// Destructor
103 delete fSlowNucleonModel;
104}
105
83b3cd6c 106void 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}
c9a8628a 110
111void AliGenSlowNucleons::Init()
112{
113 //
114 // Initialization
115 //
88a9f852 116 Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
271e6aae 117 fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
c9a8628a 118 fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
88a9f852 119 //printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
c9a8628a 120 if (fDebug) {
c8f7f6f9 121 fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
122 fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
5a9b7f8e 123 fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
124 }
125 //
126 // non-uniform cos(theta) distribution
127 //
128 if(fThetaDistribution != 0) {
129 fCosTheta = new TF1("fCosTheta",
130 "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
131 -1., 1.);
c9a8628a 132 }
133}
134
135void AliGenSlowNucleons::FinishRun()
136{
f687abfe 137// End of run action
138// Show histogram for debugging if requested.
c9a8628a 139 if (fDebug) {
c8f7f6f9 140 TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
141 c->Divide(2,1);
142 c->cd(1);
58776b75 143 fDebugHist1->Draw("colz");
c8f7f6f9 144 c->cd(2);
145 fDebugHist2->Draw();
5a9b7f8e 146 c->cd(3);
147 fCosThetaGrayHist->Draw();
c9a8628a 148 }
149}
150
151
152void AliGenSlowNucleons::Generate()
153{
154 //
155 // Generate one event
156 //
157 //
158 // Communication with Gray Particle Model
159 //
1c22a81e 160 if (fCollisionGeometry) {
161 Float_t b = fCollisionGeometry->ImpactParameter();
230d85c6 162 // Int_t nn = fCollisionGeometry->NN();
163 // Int_t nwn = fCollisionGeometry->NwN();
164 // Int_t nnw = fCollisionGeometry->NNw();
165 // Int_t nwnw = fCollisionGeometry->NwNw();
88a9f852 166
58776b75 167 //fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
168 fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 169 if (fDebug) {
88a9f852 170 //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
230d85c6 171 printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
58776b75 172 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
c8f7f6f9 173 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
1c22a81e 174 }
175 }
c9a8628a 176
177 //
88a9f852 178 Float_t p[3] = {0., 0., 0.}, theta=0;
c9a8628a 179 Float_t origin[3] = {0., 0., 0.};
21391258 180 Float_t time = 0.;
c9a8628a 181 Float_t polar [3] = {0., 0., 0.};
99f5114f 182 Int_t nt, i, j;
c9a8628a 183 Int_t kf;
99f5114f 184
185 if(fVertexSmear == kPerEvent) {
186 Vertex();
187 for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258 188 time = fTime;
99f5114f 189 } // if kPerEvent
c9a8628a 190//
191// Gray protons
192//
193 fCharge = 1;
194 kf = kProton;
1c22a81e 195 for(i = 0; i < fNgp; i++) {
5a9b7f8e 196 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
197 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 198 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 199 time, kPNoProcess, nt, 1.,-2);
c9a8628a 200 KeepTrack(nt);
201 }
202//
203// Gray neutrons
204//
205 fCharge = 0;
206 kf = kNeutron;
1c22a81e 207 for(i = 0; i < fNgn; i++) {
5a9b7f8e 208 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
209 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 210 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 211 time, kPNoProcess, nt, 1.,-2);
c9a8628a 212 KeepTrack(nt);
213 }
214//
215// Black protons
216//
217 fCharge = 1;
218 kf = kProton;
1c22a81e 219 for(i = 0; i < fNbp; i++) {
5a9b7f8e 220 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 221 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 222 time, kPNoProcess, nt, 1.,-1);
c9a8628a 223 KeepTrack(nt);
224 }
225//
226// Black neutrons
227//
228 fCharge = 0;
229 kf = kNeutron;
1c22a81e 230 for(i = 0; i < fNbn; i++) {
5a9b7f8e 231 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 232 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 233 time, kPNoProcess, nt, 1.,-1);
c9a8628a 234 KeepTrack(nt);
235 }
236}
237
238
239
240
5a9b7f8e 241void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
242 Double_t beta, Float_t* q, Float_t &theta)
f687abfe 243
244{
c9a8628a 245/*
246 Emit a slow nucleon with "temperature" T [GeV],
247 from a source moving with velocity beta
248 Three-momentum [GeV/c] is given back in q[3]
249*/
250
88a9f852 251 //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
230d85c6 252
5a9b7f8e 253 Double_t m, pmax, p, f, phi;
c9a8628a 254 TDatabasePDG * pdg = TDatabasePDG::Instance();
255 const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
256 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
257
258 /* Select nucleon type */
259 if(charge == 0) m = kMassNeutron;
260 else m = kMassProton;
261
262 /* Momentum at maximum of Maxwell-distribution */
263
88a9f852 264 pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
c9a8628a 265
266 /* Try until proper momentum */
267 /* for lack of primitive function of the Maxwell-distribution */
268 /* a brute force trial-accept loop, normalized at pmax */
269
270 do
271 {
272 p = Rndm() * fPmax;
273 f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
274 }
275 while(f < Rndm());
276
5a9b7f8e 277 /* Spherical symmetric emission for black particles (beta=0)*/
278 if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
279 /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
280 else if(fThetaDistribution!=0){
281 Double_t costheta = fCosTheta->GetRandom();
282 theta = TMath::ACos(costheta);
283 }
284 //
c9a8628a 285 phi = 2. * TMath::Pi() * Rndm();
286
88a9f852 287
c9a8628a 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
88a9f852 293
c9a8628a 294 /* Transform to system of the target nucleus */
295 /* beta is passed as negative, because the gray nucleons are slowed down */
296 Lorentz(m, -beta, q);
88a9f852 297
c9a8628a 298 /* Transform to laboratory system */
299 Lorentz(m, fBeta, q);
83b3cd6c 300 q[2] *= fProtonDirection;
88a9f852 301
c9a8628a 302}
303
304Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
305{
306/* Relativistic Maxwell-distribution */
307 Double_t ekin;
308 ekin = TMath::Sqrt(p*p+m*m)-m;
309 return (p*p * exp(-ekin/T));
310}
311
312
313void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
314{
315/* Lorentz transform in the direction of q[2] */
316
88a9f852 317 Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
318 Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
c9a8628a 319 q[2] = gamma * (q[2] + beta*energy);
88a9f852 320 //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
c9a8628a 321}
230d85c6 322