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