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