]>
Commit | Line | Data |
---|---|---|
1 | /************************************************************************* | |
2 | * Copyright(c) 1998-2048, 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 | //Author: Dariusz Miskowiec | |
17 | //Date: 2010 | |
18 | ||
19 | //============================================================================= | |
20 | // Coulomb correlation function | |
21 | //============================================================================= | |
22 | ||
23 | #include <TRandom2.h> | |
24 | #include <TComplex.h> | |
25 | #include <TFile.h> | |
26 | #include <TH2.h> | |
27 | #include "AliUnicorCoulomb.h" | |
28 | ||
29 | ClassImp(AliUnicorCoulomb) | |
30 | ||
31 | //============================================================================= | |
32 | AliUnicorCoulomb::AliUnicorCoulomb(int zz, double mass, double R) : TGraph(1001) { | |
33 | ||
34 | // constructor | |
35 | // zz = +1 for pi+pi+, -1 for pi-pi-, etc. | |
36 | // mass is the reduced mass: m1*m2/(m1+m2) | |
37 | ||
38 | SetMarkerStyle(20); | |
39 | int nstep = 10000; | |
40 | ||
41 | double rx = R; | |
42 | double ry = R; | |
43 | double rz = R; | |
44 | double r0 = 0; | |
45 | TRandom2 ran; | |
46 | ||
47 | if (zz>0) SetPoint(0,0,0); | |
48 | else SetPoint(0,0,1e9); | |
49 | for (int i=1; i<1000; i++) { | |
50 | double qinv = i/1000.0; | |
51 | double k = qinv/2.0; | |
52 | double sum = 0; | |
53 | for (int j=0; j<nstep; j++) { | |
54 | double x0 = ran.Gaus(0.0,rx); | |
55 | double y0 = ran.Gaus(0.0,ry); | |
56 | double z0 = ran.Gaus(0.0,rz); | |
57 | double t0 = ran.Gaus(0.0,r0); | |
58 | double x1 = ran.Gaus(0.0,rx); | |
59 | double y1 = ran.Gaus(0.0,ry); | |
60 | double z1 = ran.Gaus(0.0,rz); | |
61 | double t1 = ran.Gaus(0.0,r0); | |
62 | double dx = x1-x0; | |
63 | double dy = y1-y0; | |
64 | double dz = z1-z0; | |
65 | double dt = t1-t0; | |
66 | dz += dt*k/mass; | |
67 | sum += WaveFunction2(zz,mass,k,dx,dy,dz); | |
68 | } | |
69 | SetPoint(i,qinv,sum/nstep); | |
70 | } | |
71 | SetPoint(1000,1000,1); | |
72 | } | |
73 | //============================================================================= | |
74 | double AliUnicorCoulomb::WaveFunction2(int zz, double mass, double k, double x, double y, double z) { | |
75 | ||
76 | // Square of normalized Coulomb wave function. | |
77 | // Hypergeometrical function diverges for numerical reasons for large Qinv, z and x. | |
78 | // Out of the convergence limit I will return 1.0 | |
79 | ||
80 | double qinv12 = 1.2 * 2.0 * k; | |
81 | double p0 = -9.1569/qinv12; | |
82 | double p2 = 2.735e-2*qinv12; | |
83 | double wf2 = 1.0; | |
84 | int converg = (z>p0+p2*(x*x+y*y)); | |
85 | if (converg) { | |
86 | TComplex co = WaveFunction(zz,mass,k,x,y,z); | |
87 | wf2 = co.Rho2()*Gamow(zz,mass,k); | |
88 | } | |
89 | return wf2; | |
90 | } | |
91 | //============================================================================= | |
92 | TComplex AliUnicorCoulomb::WaveFunction(int zz, double mass, double k, double x, double y, double z) { | |
93 | ||
94 | // Coulomb wave function in parabolic coordinates | |
95 | // Merzbacher, Quantum Mechanics, p.248 | |
96 | // Landau-Lifshitz, Quantum Mechanics (Non-Relativistic Theory), p.518 | |
97 | ||
98 | double r = sqrt(x*x+y*y+z*z); | |
99 | // double ksi = r+z; | |
100 | double eta = r-z; | |
101 | TComplex jjj(0,1); | |
102 | TComplex a = -zz*mass/k/137.036*jjj; | |
103 | TComplex b(1.0,0.0); | |
104 | TComplex c = jjj*k*eta/0.197327; | |
105 | TComplex wf = TComplex::Exp(jjj*k*z/0.197327)*F1(a,b,c); | |
106 | //printf("%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f \n",a.Re(),a.Im(),b.Re(),b.Im(),c.Re(),c.Im(),F1(a,b,c).Rho()); | |
107 | return wf; | |
108 | } | |
109 | //============================================================================= | |
110 | double AliUnicorCoulomb::Gamow(int zz, double m, double k) { | |
111 | ||
112 | // Squared Coulomb wave function in infinity = Gamow | |
113 | // zz - product of charges | |
114 | // m - reduced mass | |
115 | // k = Qinv/2 for equal masses | |
116 | ||
117 | double n = -zz*m/k/137.036; | |
118 | double a = 2*TMath::Pi()*n; | |
119 | double ga = a/(1.0-exp(-a)); | |
120 | return ga; | |
121 | } | |
122 | //============================================================================= | |
123 | TComplex AliUnicorCoulomb::F1(TComplex alpha, TComplex gamma, TComplex z) { | |
124 | ||
125 | // confluent hypergeometric function 1F1 | |
126 | ||
127 | double prec = 0.0000001; | |
128 | TComplex term(1.0,0.0); | |
129 | TComplex sum(1.0,0.0); | |
130 | for (int n=1; n<1000; n++) { | |
131 | double u = (double) n; | |
132 | TComplex v = (alpha+u-1.0)/(gamma+u-1.0); | |
133 | TComplex w = z/u; | |
134 | term *= v; | |
135 | term *= w; | |
136 | sum += term; | |
137 | // printf("%10d %10.3f %10.3f %10.3f %10.3f\n", n, v.Rho(), w.Rho(), term.Rho(), sum.Rho()); | |
138 | if (TComplex::Abs(term/sum) < prec) return sum; | |
139 | } | |
140 | printf("F1 Maximum number of iterations reached\n"); | |
141 | return sum; | |
142 | } | |
143 | //============================================================================= | |
144 | void AliUnicorCoulomb::Makehist(int zz, double m, const char *outfil) { | |
145 | ||
146 | // Make a two-dim histogram coulomb(R,Q) and store it on outfil. | |
147 | // Later to be used via TH2::Interpolate(). | |
148 | // zz and m are the product of charges and the reduced mass. | |
149 | ||
150 | TH2D *hi = new TH2D("co","coulomb",11,-0.5,10.5,999,0.0005,0.9995); | |
151 | hi->SetXTitle("R (fm)"); | |
152 | hi->SetYTitle("Q (GeV/c)"); | |
153 | AliUnicorCoulomb *co = 0; | |
154 | for (int i=1; i<=hi->GetNbinsX(); i++) { | |
155 | double R = hi->GetXaxis()->GetBinCenter(i); | |
156 | printf("R = %.1f fm\n",R); | |
157 | co = new AliUnicorCoulomb(zz, m, R); | |
158 | for (int j=1; j<=hi->GetNbinsY(); j++) { | |
159 | double Q = hi->GetYaxis()->GetBinCenter(j); | |
160 | hi->SetBinContent(i,j,co->Eval(Q)); | |
161 | } | |
162 | delete co; | |
163 | } | |
164 | TFile::Open(outfil,"recreate"); | |
165 | hi->Write(); | |
166 | gFile->Close(); | |
167 | } | |
168 | //============================================================================= |