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