TuneA settings removed from kPyJets
[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>
6218004b 33#include <TParticle.h>
c9a8628a 34
74cfba91 35#include "AliConst.h"
c9a8628a 36#include "AliCollisionGeometry.h"
6218004b 37#include "AliStack.h"
87f95355 38#include "AliRun.h"
39#include "AliMC.h"
c9a8628a 40#include "AliGenSlowNucleons.h"
41#include "AliSlowNucleonModel.h"
42
706938e6 43ClassImp(AliGenSlowNucleons)
1c56e311 44
c9a8628a 45
1c56e311 46AliGenSlowNucleons::AliGenSlowNucleons()
47 :AliGenerator(-1),
48 fCMS(0.),
49 fMomentum(0.),
50 fBeta(0.),
51 fPmax (0.),
1c56e311 52 fCharge(0),
61cf877f 53 fProtonDirection(1.),
1c56e311 54 fTemperatureG(0.),
55 fBetaSourceG(0.),
56 fTemperatureB(0.),
57 fBetaSourceB(0.),
58 fNgp(0),
59 fNgn(0),
60 fNbp(0),
61 fNbn(0),
62 fDebug(0),
63 fDebugHist1(0),
64 fDebugHist2(0),
65 fThetaDistribution(),
66 fCosThetaGrayHist(),
67 fCosTheta(),
74cfba91 68 fBeamCrossingAngle(0.),
69 fBeamDivergence(0.),
70 fBeamDivEvent(0.),
c3e58d2c 71 fSmearMode(2),
1c56e311 72 fSlowNucleonModel(0)
c9a8628a 73{
74// Default constructor
c9a8628a 75 fCollisionGeometry = 0;
76}
77
78AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
1c56e311 79 :AliGenerator(npart),
80 fCMS(14000.),
81 fMomentum(0.),
82 fBeta(0.),
83 fPmax (10.),
1c56e311 84 fCharge(1),
61cf877f 85 fProtonDirection(1.),
230d85c6 86 fTemperatureG(0.05),
1c56e311 87 fBetaSourceG(0.05),
230d85c6 88 fTemperatureB(0.005),
1c56e311 89 fBetaSourceB(0.),
90 fNgp(0),
91 fNgn(0),
92 fNbp(0),
93 fNbn(0),
94 fDebug(0),
95 fDebugHist1(0),
96 fDebugHist2(0),
97 fThetaDistribution(),
98 fCosThetaGrayHist(),
99 fCosTheta(),
74cfba91 100 fBeamCrossingAngle(0.),
101 fBeamDivergence(0.),
102 fBeamDivEvent(0.),
c3e58d2c 103 fSmearMode(2),
1c56e311 104 fSlowNucleonModel(new AliSlowNucleonModel())
74cfba91 105
c9a8628a 106{
107// Constructor
108 fName = "SlowNucleons";
109 fTitle = "Generator for gray particles in pA collisions";
c9a8628a 110 fCollisionGeometry = 0;
c9a8628a 111}
112
5a9b7f8e 113//____________________________________________________________
c9a8628a 114AliGenSlowNucleons::~AliGenSlowNucleons()
115{
116// Destructor
117 delete fSlowNucleonModel;
118}
119
83b3cd6c 120void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
121// Set direction of the proton to change between pA (1) and Ap (-1)
122 fProtonDirection = dir / TMath::Abs(dir);
123}
c9a8628a 124
125void AliGenSlowNucleons::Init()
126{
127 //
128 // Initialization
129 //
88a9f852 130 Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
271e6aae 131 fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
c9a8628a 132 fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
88a9f852 133 //printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
c9a8628a 134 if (fDebug) {
c8f7f6f9 135 fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
136 fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
5a9b7f8e 137 fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
138 }
139 //
140 // non-uniform cos(theta) distribution
141 //
142 if(fThetaDistribution != 0) {
143 fCosTheta = new TF1("fCosTheta",
144 "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
145 -1., 1.);
c9a8628a 146 }
74cfba91 147
148 printf("\n AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
c9a8628a 149}
150
151void AliGenSlowNucleons::FinishRun()
152{
f687abfe 153// End of run action
154// Show histogram for debugging if requested.
c9a8628a 155 if (fDebug) {
c8f7f6f9 156 TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
157 c->Divide(2,1);
158 c->cd(1);
58776b75 159 fDebugHist1->Draw("colz");
c8f7f6f9 160 c->cd(2);
161 fDebugHist2->Draw();
5a9b7f8e 162 c->cd(3);
163 fCosThetaGrayHist->Draw();
c9a8628a 164 }
165}
166
167
168void AliGenSlowNucleons::Generate()
169{
170 //
171 // Generate one event
172 //
173 //
174 // Communication with Gray Particle Model
175 //
1c22a81e 176 if (fCollisionGeometry) {
177 Float_t b = fCollisionGeometry->ImpactParameter();
230d85c6 178 // Int_t nn = fCollisionGeometry->NN();
179 // Int_t nwn = fCollisionGeometry->NwN();
180 // Int_t nnw = fCollisionGeometry->NNw();
181 // Int_t nwnw = fCollisionGeometry->NwNw();
88a9f852 182
c3e58d2c 183 // (1) Sikler' model
184 if(fSmearMode==0) fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
185 // (2) Model inspired on exp. data at lower energy (Gallio-Oppedisano)
186 // --- smearing the Ncoll fron generator used as input
187 else if(fSmearMode==1) fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
188 // --- smearing directly Nslow
189 else if(fSmearMode==2) fSlowNucleonModel->GetNumberOfSlowNucleons2s(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e 190 if (fDebug) {
88a9f852 191 //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
230d85c6 192 printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
58776b75 193 fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
c8f7f6f9 194 fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
74cfba91 195
1c22a81e 196 }
197 }
c9a8628a 198
199 //
88a9f852 200 Float_t p[3] = {0., 0., 0.}, theta=0;
c9a8628a 201 Float_t origin[3] = {0., 0., 0.};
21391258 202 Float_t time = 0.;
c9a8628a 203 Float_t polar [3] = {0., 0., 0.};
99f5114f 204 Int_t nt, i, j;
c9a8628a 205 Int_t kf;
74cfba91 206
207 // Extracting 1 value per event for the divergence angle
208 Double_t rvec = gRandom->Gaus(0.0, 1.0);
209 fBeamDivEvent = fBeamDivergence * TMath::Abs(rvec);
210 printf("\n AliGenSlowNucleons: applying beam divergence %f mrad to slow nucleons\n",fBeamDivEvent*1000.);
211
99f5114f 212 if(fVertexSmear == kPerEvent) {
213 Vertex();
214 for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258 215 time = fTime;
99f5114f 216 } // if kPerEvent
c9a8628a 217//
218// Gray protons
219//
220 fCharge = 1;
221 kf = kProton;
1c22a81e 222 for(i = 0; i < fNgp; i++) {
5a9b7f8e 223 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
224 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 225 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 226 time, kPNoProcess, nt, 1.,-2);
c9a8628a 227 KeepTrack(nt);
87f95355 228 SetProcessID(nt,kGrayProcess);
c9a8628a 229 }
230//
231// Gray neutrons
232//
233 fCharge = 0;
234 kf = kNeutron;
1c22a81e 235 for(i = 0; i < fNgn; i++) {
5a9b7f8e 236 GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
237 if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf 238 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 239 time, kPNoProcess, nt, 1.,-2);
c9a8628a 240 KeepTrack(nt);
87f95355 241 SetProcessID(nt,kGrayProcess);
c9a8628a 242 }
243//
244// Black protons
245//
246 fCharge = 1;
247 kf = kProton;
1c22a81e 248 for(i = 0; i < fNbp; i++) {
5a9b7f8e 249 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 250 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 251 time, kPNoProcess, nt, 1.,-1);
c9a8628a 252 KeepTrack(nt);
87f95355 253 SetProcessID(nt,kBlackProcess);
c9a8628a 254 }
255//
256// Black neutrons
257//
258 fCharge = 0;
259 kf = kNeutron;
1c22a81e 260 for(i = 0; i < fNbn; i++) {
5a9b7f8e 261 GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf 262 PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75 263 time, kPNoProcess, nt, 1.,-1);
c9a8628a 264 KeepTrack(nt);
87f95355 265 SetProcessID(nt,kBlackProcess);
c9a8628a 266 }
267}
268
269
270
271
5a9b7f8e 272void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
273 Double_t beta, Float_t* q, Float_t &theta)
f687abfe 274
275{
c9a8628a 276/*
277 Emit a slow nucleon with "temperature" T [GeV],
278 from a source moving with velocity beta
279 Three-momentum [GeV/c] is given back in q[3]
280*/
281
88a9f852 282 //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
230d85c6 283
5a9b7f8e 284 Double_t m, pmax, p, f, phi;
c9a8628a 285 TDatabasePDG * pdg = TDatabasePDG::Instance();
286 const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
287 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
288
289 /* Select nucleon type */
290 if(charge == 0) m = kMassNeutron;
291 else m = kMassProton;
292
293 /* Momentum at maximum of Maxwell-distribution */
294
88a9f852 295 pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
c9a8628a 296
297 /* Try until proper momentum */
298 /* for lack of primitive function of the Maxwell-distribution */
299 /* a brute force trial-accept loop, normalized at pmax */
300
301 do
302 {
303 p = Rndm() * fPmax;
304 f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
305 }
306 while(f < Rndm());
307
5a9b7f8e 308 /* Spherical symmetric emission for black particles (beta=0)*/
309 if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
310 /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
311 else if(fThetaDistribution!=0){
312 Double_t costheta = fCosTheta->GetRandom();
313 theta = TMath::ACos(costheta);
314 }
315 //
c9a8628a 316 phi = 2. * TMath::Pi() * Rndm();
317
88a9f852 318
c9a8628a 319 /* Determine momentum components in system of the moving source */
320 q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
321 q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
322 q[2] = p * TMath::Cos(theta);
74cfba91 323 //if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
c9a8628a 324
88a9f852 325
c9a8628a 326 /* Transform to system of the target nucleus */
327 /* beta is passed as negative, because the gray nucleons are slowed down */
328 Lorentz(m, -beta, q);
74cfba91 329 //if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
330
c9a8628a 331 /* Transform to laboratory system */
332 Lorentz(m, fBeta, q);
83b3cd6c 333 q[2] *= fProtonDirection;
74cfba91 334 if(fDebug==1)printf("\n Momentum after LHC boost: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
335
336 if(fBeamCrossingAngle>0.) BeamCrossDivergence(1, q); // applying crossing angle
337 if(fBeamDivergence>0.) BeamCrossDivergence(2, q); // applying divergence
338
c9a8628a 339}
340
341Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
342{
343/* Relativistic Maxwell-distribution */
344 Double_t ekin;
345 ekin = TMath::Sqrt(p*p+m*m)-m;
346 return (p*p * exp(-ekin/T));
347}
348
349
74cfba91 350//_____________________________________________________________________________
c9a8628a 351void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
352{
353/* Lorentz transform in the direction of q[2] */
354
88a9f852 355 Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
356 Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
c9a8628a 357 q[2] = gamma * (q[2] + beta*energy);
88a9f852 358 //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
c9a8628a 359}
230d85c6 360
74cfba91 361//_____________________________________________________________________________
362void AliGenSlowNucleons::BeamCrossDivergence(Int_t iwhat, Float_t *pLab)
363{
364 // Applying beam divergence and crossing angle
365 //
366 Double_t pmod = TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]+pLab[2]*pLab[2]);
367
368 Double_t tetdiv = 0.;
369 Double_t fidiv = 0.;
370 if(iwhat==1){
371 tetdiv = fBeamCrossingAngle;
372 fidiv = k2PI/4.;
373 }
374 else if(iwhat==2){
375 tetdiv = fBeamDivEvent;
376 fidiv = (gRandom->Rndm())*k2PI;
377 }
378
379 Double_t tetpart = TMath::ATan2(TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]), pLab[2]);
380 Double_t fipart=0.;
381 if(TMath::Abs(pLab[1])>0. || TMath::Abs(pLab[0])>0.) fipart = TMath::ATan2(pLab[1], pLab[0]);
382 if(fipart<0.) {fipart = fipart+k2PI;}
383 tetdiv = tetdiv*kRaddeg;
384 fidiv = fidiv*kRaddeg;
385 tetpart = tetpart*kRaddeg;
386 fipart = fipart*kRaddeg;
387
388 Double_t angleSum[2]={0., 0.};
389 AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum);
390
391 Double_t tetsum = angleSum[0];
392 Double_t fisum = angleSum[1];
393 //printf("tetpart %f fipart %f tetdiv %f fidiv %f angleSum %f %f\n",tetpart,fipart,tetdiv,fidiv,angleSum[0],angleSum[1]);
394 tetsum = tetsum*kDegrad;
395 fisum = fisum*kDegrad;
396
397 pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum);
398 pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum);
399 pLab[2] = pmod*TMath::Cos(tetsum);
400 if(fDebug==1){
401 if(iwhat==1) printf(" Beam crossing angle %f mrad ", fBeamCrossingAngle*1000.);
402 else if(iwhat==2) printf(" Beam divergence %f mrad ", fBeamDivEvent*1000.);
403 printf(" p = (%f, %f, %f)\n",pLab[0],pLab[1],pLab[2]);
404 }
405}
406
407//_____________________________________________________________________________
408void AliGenSlowNucleons::AddAngle(Double_t theta1, Double_t phi1,
409 Double_t theta2, Double_t phi2, Double_t *angleSum)
410{
411 // Calculating the sum of 2 angles
412 Double_t temp = -1.;
413 Double_t conv = 180./TMath::ACos(temp);
414
415 Double_t ct1 = TMath::Cos(theta1/conv);
416 Double_t st1 = TMath::Sin(theta1/conv);
417 Double_t cp1 = TMath::Cos(phi1/conv);
418 Double_t sp1 = TMath::Sin(phi1/conv);
419 Double_t ct2 = TMath::Cos(theta2/conv);
420 Double_t st2 = TMath::Sin(theta2/conv);
421 Double_t cp2 = TMath::Cos(phi2/conv);
422 Double_t sp2 = TMath::Sin(phi2/conv);
423 Double_t cx = ct1*cp1*st2*cp2+st1*cp1*ct2-sp1*st2*sp2;
424 Double_t cy = ct1*sp1*st2*cp2+st1*sp1*ct2+cp1*st2*sp2;
425 Double_t cz = ct1*ct2-st1*st2*cp2;
426
427 Double_t rtetsum = TMath::ACos(cz);
428 Double_t tetsum = conv*rtetsum;
429 if(TMath::Abs(tetsum)<1e-4 || tetsum==180.) return;
430
431 temp = cx/TMath::Sin(rtetsum);
432 if(temp>1.) temp=1.;
433 if(temp<-1.) temp=-1.;
434 Double_t fisum = conv*TMath::ACos(temp);
435 if(cy<0) {fisum = 360.-fisum;}
436
437 angleSum[0] = tetsum;
438 angleSum[1] = fisum;
439}
440
87f95355 441//_____________________________________________________________________________
442void AliGenSlowNucleons::SetProcessID(Int_t nt, UInt_t process)
443{
444 // Tag the particle as
445 // gray or black
446 if (fStack)
447 fStack->Particle(nt)->SetUniqueID(process);
448 else
449 gAlice->GetMCApp()->Particle(nt)->SetUniqueID(process);
450}