C++ Interface to Tauola
SANCtable.cxx
1#include "SANCtable.h"
2#include "unistd.h"
3
4int SANCtable::ns1=0,SANCtable::ns2=0,SANCtable::ns3=0,SANCtable::ncos=0;
5double SANCtable::smin1=0,SANCtable::smax1=0;
6double SANCtable::smin2=0,SANCtable::smax2=0;
7double SANCtable::smin3=0,SANCtable::smax3=0;
8int SANCtable::iqed=0,SANCtable::iew=0,SANCtable::iborn=0;
9int SANCtable::gfscheme=0,SANCtable::ifgg=0;
10double SANCtable::nc=0,SANCtable::fc=0,SANCtable::tlmu2=0;
11
12//-----------------------------------------------------------------------------
13//-Static part-----------------------------------------------------------------
14//-----------------------------------------------------------------------------
15void SANCtable::setDimensions(int n1, int n2, int n3, int nc)
16{
17 ns1=n1;
18 ns2=n2;
19 ns3=n3;
20 ncos=nc;
21}
22void SANCtable::setRanges(double sn1, double sx1, double sn2, double sx2, double sn3, double sx3)
23{
24 smin1=sn1*sn1;
25 smax1=sx1*sx1;
26 smin2=sn2*sn2;
27 smax2=sx2*sx2;
28 smin3=sn3*sn3;
29 smax3=sx3*sx3;
30}
31void SANCtable::setFlags()
32{
33 iqed=0;iew=1;iborn=0;gfscheme=0;ifgg=1;
34 nc=1.0;fc=3.0;tlmu2=1e-5;
35 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
36 int buf=0;
37 printconsts_(&buf);
38}
39//-----------------------------------------------------------------------------
40//-----------------------------------------------------------------------------
41//-----------------------------------------------------------------------------
42void SANCtable::setFixedLength(int prc)
43{
44 if(prc==0)
45 {
46 f.precision(6);
47 f.unsetf(ios::fixed);
48 }
49 else
50 {
51 f.precision(prc);
52 f.setf(ios::fixed);
53 }
54}
55
56bool SANCtable::addHeader()
57{
58 char path[255],ftime[255];
59 time_t utime;
60 if(ns1==0 || smin1==0) return false;
61 f<<"Dimensions: "<<ns1<<" "<<ns2<<" "<<ns3<<" "<<ncos<<endl;
62 f<<"Ranges: "<<log(smin1)<<" "<<log(smax1)<<" "<<smin2<<" "<<smax2<<" "<<smin3<<" "<<smax3<<endl;
63 getcwd(path,255);
64 time(&utime);
65 strftime(ftime,255,"%d %b %Y, %H:%M:%S, GMT%z",localtime(&utime));
66 f<<"Timestamp: "<<ftime<<endl;
67 f<<"Path: "<<path<<endl<<endl;
68 return true;
69}
70
71bool SANCtable::addFile(const char *name)
72{
73 ifstream d(name);
74 unsigned char c;
75 if(!d.is_open()) return false;
76 while(!d.eof())
77 {
78 d>>noskipws>>c;
79 f.put(c);
80 }
81 return true;
82}
83
84void SANCtable::addRange(int rangeNo,bool isLog)
85{
86 double min=0,max=0,steps=0;
87 ofstream::fmtflags last = cout.setf(ios::fixed);
88 f<<"\nBeginRange"<<rangeNo<<endl;
89 cout<<"\nBeginRange"<<rangeNo<<endl;
90 switch(rangeNo)
91 {
92 case 1: min=smin1; max=smax1; steps=ns1; break;
93 case 2: min=smin2; max=smax2; steps=ns2; break;
94 case 3: min=smin3; max=smax3; steps=ns3; break;
95 }
96 double step= (isLog) ? (log(max)-log(min))/(steps-1) : (max-min)/(steps-1);
97 double sloop=0;
98 for(int i=0;i<steps;i++)
99 {
100 sloop= (isLog) ? exp(log(min)+i*step) : min+i*step;
101 cout<<sloop<<endl;
102 for(int j=0;j<ncos;j++)
103 {
104 double costhetloop=-1.+1.0/ncos+j*2./ncos;
105 iew=0;
106 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
107 double borndivv = Rcalc(flav,sloop,costhetloop);
108
109 iew = (born) ? 0 : 1;
110 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
111 double divv = Rcalc(flav,sloop,costhetloop);
112
113 for(int k=0;k<4;k++)
114 for(int l=0;l<4;l++)
115 f<<R[k][l]<<" ";
116 f<<divv<<" "<<borndivv<<endl;
117 }
118 }
119 cout.flags(last);
120}
121
122void SANCtable::open(const char *name)
123{
124 f.open(name);
125 isOpen = f.is_open();
126}
127
128void SANCtable::close()
129{
130 f<<"End"<<endl;
131 f.close();
132 isOpen=false;
133}
134//-----------------------------------------------------------------------------
135//-The main computation module-------------------------------------------------
136//-----------------------------------------------------------------------------
137double SANCtable::Rcalc(int flav, double s,double cosf)
138{
139 double sig[4][4][2],sum=0.,divv=0.;
140 complex<double> amp[2][2][2][2][2]={0}, sigma[4][2][2]={0}, c;
141
142 sigma[0][0][0] = complex<double>(1,0);
143 sigma[0][1][1] = complex<double>(1,0);
144
145 sigma[1][0][1] = complex<double>(1,0);
146 sigma[1][1][0] = complex<double>(1,0);
147
148 sigma[2][0][1] = complex<double>(0,-1);
149 sigma[2][1][0] = complex<double>(0,1);
150
151 sigma[3][0][0] = complex<double>(1,0);
152 sigma[3][1][1] = complex<double>(-1,0);
153
154 double mta,conhc,pi;
155 paraget_(&mta,&conhc,&pi);
156
157 double betaf = sqrt(1e0-4e0*mta*mta/s);
158 double t = mta*mta - s/2*(1e0-betaf*cosf);
159 double u = mta*mta - s/2*(1e0+betaf*cosf);
160 for(int iz = 0;iz<2;iz++)
161 {
162 for(int L1=1;L1<=2;L1++)
163 {
164 for(int L2=1;L2<=2;L2++)
165 {
166 for(int L3=1;L3<=2;L3++)
167 {
168 for(int L4=1;L4<=2;L4++)
169 {
170 double har=0,hai=0;
171 if(flav==1) downdown_(&L1,&L2,&L3,&L4,&s,&t,&u,&iz,&har,&hai);
172 else upup_(&L1,&L2,&L3,&L4,&s,&t,&u,&iz,&har,&hai);
173 amp[L1-1][L2-1][L3-1][L4-1][iz] = complex<double>(har,hai);
174 }
175 }
176 }
177 }
178 }
179
180 for(int i=0;i<4;i++){
181 for(int j=0;j<4;j++){
182 for(int iz = 0;iz<2;iz++)
183 {
184 sum = 0.0;
185 for(int L1=1;L1<=2;L1++)
186 {
187 for(int L2=1;L2<=2;L2++)
188 {
189 for(int L3=1;L3<=2;L3++)
190 {
191 for(int L4=1;L4<=2;L4++)
192 for(int M3=1;M3<=2;M3++)
193 for(int M4=1;M4<=2;M4++)
194 {
195 c = amp[L1-1][L2-1][L3-1][L4-1][iz] *
196 conj(amp[L1-1][L2-1][M3-1][M4-1][iz]) *
197 sigma[i][M3-1][L3-1] *
198 sigma[j][M4-1][L4-1];
199 sum+=real(c);
200 }
201 }
202 }
203 }
204 sig[i][j][iz]=sum;
205 }
206
207 R[i][j] = conhc * // to pbarn
208 nc/fc*1e0/2.0/s *
209 1e0/4 * // spin sum
210 (sig[i][j][1] - sig[i][j][0]) * // |Amp|^2 - linearized
211 betaf/16/pi; // phase_space/dcos{theta}
212
213 if(i==0 && j==0) divv=R[0][0];
214
215 R[i][j]=R[i][j]/divv;
216
217 }
218 }
219
220
221 for(int i=0;i<4;i++) // rotations for tau+ see tau+ KW and SANC frames
222 {
223 double xx= R[i][1];
224 R[i][1] = R[i][2];
225 R[i][2] = -xx;
226 }
227
228 for(int j=0;j<4;j++) // rotations for tau- see tau- KW and SANC frames
229 {
230 double xx= R[1][j];
231 R[1][j] = -R[2][j];
232 R[2][j] = -xx;
233 R[3][j] = -R[3][j];
234 }
235 for(int i=0;i<4;i++) // rotations for tau+/tau- see reaction KW and SANC frames
236 {
237 R[i][1] = -R[i][1];
238 R[i][2] = -R[i][2];
239 }
240
241 for(int j=0;j<4;j++) // rotations for tau+/tau- see reaction KW and SANC frames
242 {
243 R[1][j] = -R[1][j];
244 R[2][j] = -R[2][j];
245 }
246 /*
247 printf("\n");
248 printf("d_sigma/d_cos{theta} = %.9f \n",divv);
249 printf("s = %.9f \n" ,s);
250 printf("cosf = %.9f \n",cosf);
251 printf("\n");
252 printf("R(i,j)= \n");
253 printf("\n");
254 for(int i=0;i<4;i++)
255 {
256 for(int j=0;j<4;j++)
257 printf("%.5f ",R[i][j]);
258 printf("\n");
259 }
260 double sinf;
261 double MM;
262 double xnorm;
263 double Bench[4][4]={0.};
264 sinf=sqrt(1-cosf*cosf);
265 MM=sqrt(4e0*mta*mta/s);
266 xnorm=1+cosf*cosf + MM*MM*sinf*sinf;
267 Bench[0][0]=(1+cosf*cosf + MM*MM*sinf*sinf)/xnorm;
268 Bench[1][1]=(-(1- MM*MM)*sinf*sinf)/xnorm;
269 Bench[2][2]=( (1+ MM*MM)*sinf*sinf)/xnorm;
270 Bench[3][3]=(1+cosf*cosf - MM*MM*sinf*sinf)/xnorm;
271 Bench[2][3]=(2*MM*sinf*cosf)/xnorm;
272 Bench[3][2]=(2*MM*sinf*cosf)/xnorm;
273 printf("\n");
274 printf("Bench(i,j)= \n");
275 printf("\n");
276
277 for(int i=0;i<4;i++)
278 {
279 for(int j=0;j<4;j++)
280 printf("%.5f ",Bench[i][j]);
281 printf("\n");
282 }
283
284 printf("\n");
285 printf(" \n");
286 printf("Bench(i,j)=R(i,j) at low energies, where ew corr and Z contribute little \n");
287 printf("\n");
288
289 */
290
291 return divv;
292}