Beam crossing angle and divergence (Chiara)
[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
34 #include "AliConst.h"
35 #include "AliCollisionGeometry.h"
36 #include "AliGenSlowNucleons.h"
37 #include "AliSlowNucleonModel.h"
38
39 ClassImp(AliGenSlowNucleons)
40
41     
42 AliGenSlowNucleons::AliGenSlowNucleons()
43     :AliGenerator(-1),
44      fCMS(0.),
45      fMomentum(0.),
46      fBeta(0.),
47      fPmax (0.),
48      fCharge(0),
49      fProtonDirection(1.),
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(),
64      fBeamCrossingAngle(0.),
65      fBeamDivergence(0.),
66      fBeamDivEvent(0.),
67      fSlowNucleonModel(0)
68 {
69 // Default constructor
70     fCollisionGeometry = 0;
71 }
72
73 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
74     :AliGenerator(npart),
75      fCMS(14000.),
76      fMomentum(0.),
77      fBeta(0.),
78      fPmax (10.),
79      fCharge(1),
80      fProtonDirection(1.),
81      fTemperatureG(0.05), 
82      fBetaSourceG(0.05),
83      fTemperatureB(0.005),
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(),
95      fBeamCrossingAngle(0.),
96      fBeamDivergence(0.),
97      fBeamDivEvent(0.),
98      fSlowNucleonModel(new AliSlowNucleonModel())
99
100 {
101 // Constructor
102     fName  = "SlowNucleons";
103     fTitle = "Generator for gray particles in pA collisions";
104     fCollisionGeometry = 0;
105 }
106
107 //____________________________________________________________
108 AliGenSlowNucleons::~AliGenSlowNucleons()
109 {
110 // Destructor
111     delete  fSlowNucleonModel;
112 }
113
114 void 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 }
118
119 void AliGenSlowNucleons::Init()
120 {
121   //
122   // Initialization
123   //
124     Double_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
125     fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
126     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
127     //printf("  fMomentum %f    fBeta %1.10f\n",fMomentum, fBeta);
128     if (fDebug) {
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.);
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.);
140     }
141    
142     printf("\n  AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
143 }
144
145 void AliGenSlowNucleons::FinishRun()
146 {
147 // End of run action
148 // Show histogram for debugging if requested.
149     if (fDebug) {
150         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
151         c->Divide(2,1);
152         c->cd(1);
153         fDebugHist1->Draw("colz");
154         c->cd(2);
155         fDebugHist2->Draw();
156         c->cd(3);
157         fCosThetaGrayHist->Draw();
158     }
159 }
160
161
162 void AliGenSlowNucleons::Generate()
163 {
164   //
165   // Generate one event
166   //
167   //
168   // Communication with Gray Particle Model 
169   // 
170     if (fCollisionGeometry) {
171         Float_t b   = fCollisionGeometry->ImpactParameter();
172         //      Int_t  nn   = fCollisionGeometry->NN();
173         //      Int_t  nwn  = fCollisionGeometry->NwN();
174         //      Int_t  nnw  = fCollisionGeometry->NNw();
175         //      Int_t  nwnw = fCollisionGeometry->NwNw();
176         
177         //fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
178         fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
179         if (fDebug) {
180             //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
181             printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
182             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
183             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
184
185         }
186     }     
187
188    //
189     Float_t p[3] = {0., 0., 0.}, theta=0;
190     Float_t origin[3] = {0., 0., 0.};
191     Float_t time = 0.;
192     Float_t polar [3] = {0., 0., 0.};    
193     Int_t nt, i, j;
194     Int_t kf;
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
201     if(fVertexSmear == kPerEvent) {
202         Vertex();
203         for (j=0; j < 3; j++) origin[j] = fVertex[j];
204         time = fTime;
205     } // if kPerEvent
206 //
207 //  Gray protons
208 //
209     fCharge = 1;
210     kf = kProton;    
211     for(i = 0; i < fNgp; i++) {
212         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
213         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
214         PushTrack(fTrackIt, -1, kf, p, origin, polar,
215                  time, kPNoProcess, nt, 1.,-2);
216         KeepTrack(nt);
217     }
218 //
219 //  Gray neutrons
220 //
221     fCharge = 0;
222     kf = kNeutron;    
223     for(i = 0; i < fNgn; i++) {
224         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
225         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
226         PushTrack(fTrackIt, -1, kf, p, origin, polar,
227                  time, kPNoProcess, nt, 1.,-2);
228         KeepTrack(nt);
229     }
230 //
231 //  Black protons
232 //
233     fCharge = 1;
234     kf = kProton;    
235     for(i = 0; i < fNbp; i++) {
236         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
237         PushTrack(fTrackIt, -1, kf, p, origin, polar,
238                  time, kPNoProcess, nt, 1.,-1);
239         KeepTrack(nt);
240     }
241 //
242 //  Black neutrons
243 //
244     fCharge = 0;
245     kf = kNeutron;    
246     for(i = 0; i < fNbn; i++) {
247         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
248         PushTrack(fTrackIt, -1, kf, p, origin, polar,
249                  time, kPNoProcess, nt, 1.,-1);
250         KeepTrack(nt);
251     }
252 }
253
254
255
256
257 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, 
258         Double_t beta, Float_t* q, Float_t &theta)
259
260 {
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
267  //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
268  
269  Double_t m, pmax, p, f, phi;
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
280  pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
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
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  //
301  phi   = 2. * TMath::Pi() * Rndm();
302
303
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);
308  //if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
309
310
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);
314  //if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
315
316  /* Transform to laboratory system */
317  Lorentz(m, fBeta, q);
318  q[2] *= fProtonDirection; 
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  
324 }
325
326 Double_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
335 //_____________________________________________________________________________
336 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
337 {
338 /* Lorentz transform in the direction of q[2] */
339  
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]);
342     q[2] = gamma * (q[2] + beta*energy);
343     //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
344 }
345
346 //_____________________________________________________________________________
347 void 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 //_____________________________________________________________________________
393 void  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