doxy: install THtml converter
[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx
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 ClassImp(AliGenSlowNucleons)
44
45     
46 AliGenSlowNucleons::AliGenSlowNucleons()
47     :AliGenerator(-1),
48      fCMS(0.),
49      fMomentum(0.),
50      fBeta(0.),
51      fPmax (0.),
52      fCharge(0),
53      fProtonDirection(1.),
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(),
68      fBeamCrossingAngle(0.),
69      fBeamDivergence(0.),
70      fBeamDivEvent(0.),
71      fSmearMode(2),
72      fSlowNucleonModel(0)
73 {
74 // Default constructor
75     fCollisionGeometry = 0;
76 }
77
78 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
79     :AliGenerator(npart),
80      fCMS(14000.),
81      fMomentum(0.),
82      fBeta(0.),
83      fPmax (10.),
84      fCharge(1),
85      fProtonDirection(1.),
86      fTemperatureG(0.05), 
87      fBetaSourceG(0.05),
88      fTemperatureB(0.005),
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(),
100      fBeamCrossingAngle(0.),
101      fBeamDivergence(0.),
102      fBeamDivEvent(0.),
103      fSmearMode(2),
104      fSlowNucleonModel(new AliSlowNucleonModel())
105
106 {
107 // Constructor
108     fName  = "SlowNucleons";
109     fTitle = "Generator for gray particles in pA collisions";
110     fCollisionGeometry = 0;
111 }
112
113 //____________________________________________________________
114 AliGenSlowNucleons::~AliGenSlowNucleons()
115 {
116 // Destructor
117     delete  fSlowNucleonModel;
118 }
119
120 void 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 }
124
125 void AliGenSlowNucleons::Init()
126 {
127   //
128   // Initialization
129   //
130     Double_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
131     fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
132     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
133     //printf("  fMomentum %f    fBeta %1.10f\n",fMomentum, fBeta);
134     if (fDebug) {
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.);
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.);
146     }
147    
148     printf("\n  AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
149 }
150
151 void AliGenSlowNucleons::FinishRun()
152 {
153 // End of run action
154 // Show histogram for debugging if requested.
155     if (fDebug) {
156         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
157         c->Divide(2,1);
158         c->cd(1);
159         fDebugHist1->Draw("colz");
160         c->cd(2);
161         fDebugHist2->Draw();
162         c->cd(3);
163         fCosThetaGrayHist->Draw();
164     }
165 }
166
167
168 void AliGenSlowNucleons::Generate()
169 {
170   //
171   // Generate one event
172   //
173   //
174   // Communication with Gray Particle Model 
175   // 
176     if (fCollisionGeometry) {
177         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         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);
190         if (fDebug) {
191             //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
192             printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
193             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
194             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
195
196         }
197     }     
198
199    //
200     Float_t p[3] = {0., 0., 0.}, theta=0;
201     Float_t origin[3] = {0., 0., 0.};
202     Float_t time = 0.;
203     Float_t polar [3] = {0., 0., 0.};    
204     Int_t nt, i, j;
205     Int_t kf;
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
212     if(fVertexSmear == kPerEvent) {
213         Vertex();
214         for (j=0; j < 3; j++) origin[j] = fVertex[j];
215         time = fTime;
216     } // if kPerEvent
217 //
218 //  Gray protons
219 //
220     fCharge = 1;
221     kf = kProton;    
222     for(i = 0; i < fNgp; i++) {
223         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
224         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
225         PushTrack(fTrackIt, -1, kf, p, origin, polar,
226                  time, kPNoProcess, nt, 1.,-2);
227         KeepTrack(nt);
228         SetProcessID(nt,kGrayProcess);
229     }
230 //
231 //  Gray neutrons
232 //
233     fCharge = 0;
234     kf = kNeutron;    
235     for(i = 0; i < fNgn; i++) {
236         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
237         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
238         PushTrack(fTrackIt, -1, kf, p, origin, polar,
239                  time, kPNoProcess, nt, 1.,-2);
240         KeepTrack(nt);
241         SetProcessID(nt,kGrayProcess);
242     }
243 //
244 //  Black protons
245 //
246     fCharge = 1;
247     kf = kProton;    
248     for(i = 0; i < fNbp; i++) {
249         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
250         PushTrack(fTrackIt, -1, kf, p, origin, polar,
251                  time, kPNoProcess, nt, 1.,-1);
252         KeepTrack(nt);
253         SetProcessID(nt,kBlackProcess);
254     }
255 //
256 //  Black neutrons
257 //
258     fCharge = 0;
259     kf = kNeutron;    
260     for(i = 0; i < fNbn; i++) {
261         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
262         PushTrack(fTrackIt, -1, kf, p, origin, polar,
263                  time, kPNoProcess, nt, 1.,-1);
264         KeepTrack(nt);
265         SetProcessID(nt,kBlackProcess);
266     }
267 }
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  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
295  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  do
302  {
303      p = Rndm() * fPmax;
304      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
305  }
306  while(f < Rndm());
307
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  //
316  phi   = 2. * TMath::Pi() * Rndm();
317
318
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);
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  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  Lorentz(m, fBeta, q);
333  q[2] *= fProtonDirection; 
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  
339 }
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     ekin = TMath::Sqrt(p*p+m*m)-m;
346     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     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]);
357     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 }
360
361 //_____________________________________________________________________________
362 void 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 //_____________________________________________________________________________
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   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
441 //_____________________________________________________________________________
442 void 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 }