C++ Interface to Tauola
SANCinterface.c
1#include <math.h>
2#include <stdio.h>
3#include <time.h>
4
5extern void upup_(int *L1,int *L2, int *L3,int *L4,double *s,double *t,double *u,int *iz,double *har,double *hai);
6extern void downdown_(int *L1,int *L2, int *L3,int *L4,double *s,double *t,double *u,int *iz,double *har,double *hai);
7extern void flagset_(int *iqedx,int *iewx,int *ibornx,int *gfschemex,int *ifggx,double *ncx,double *fcx,double *tlmu2x);
8extern void paraget_(double *mtax, double *conhcx, double *pix);
9extern void printconsts_(int *pmg);
10
11
12int iqed=0,iew=1,iborn=0,gfscheme=0,ifgg=1;
13double nc=1.0,fc=3.0,tlmu2=1e-5;
14int buf=0;
15double R[4][4];
16double divv;
17double borndivv=123.456;
18
19struct compl {double im,re;};
20
21struct compl mul(struct compl x, struct compl y)
22{
23 struct compl c ;
24 c.re=x.re*y.re-x.im*y.im;
25 c.im=x.re*y.im+x.im*y.re;
26 return c;
27}
28
29struct compl conju(struct compl x)
30{
31 struct compl c ;
32 c.re=x.re;
33 c.im=-x.im;
34 return c;
35}
36
37double abso(struct compl x)
38{
39 return sqrt(x.re*x.re+x.im*x.im);
40}
41
42
43void Rcalc(int flav, double sloop,double costhetloop)
44{
45 int L1,L2,L3,L4,iz,i,j,M3,M4;
46 double s,t,u;
47 double cosf,betaf,har,hai,MM,sinf,xnorm,xx;
48
49 double mta,conhc,pi;
50
51 double sum,sig[4][4][2], Bench[4][4]={0};
52 double ha[2];
53 struct compl amp[2][2][2][2][2]={0}, sigma[4][2][2]={0}, c;
54 sigma[0][0][0].re=1;
55 sigma[0][1][1].re=1;
56
57 sigma[1][0][1].re=1;
58 sigma[1][1][0].re=1;
59
60 sigma[2][0][1].im=-1;
61 sigma[2][1][0].im= 1;
62
63 sigma[3][0][0].re=1;
64 sigma[3][1][1].re=-1;
65
66 s=sloop;
67 cosf=costhetloop;
68// --------------------------
69 paraget_(&mta,&conhc,&pi);
70 betaf = sqrt(1e0-4e0*mta*mta/s);
71 t = mta*mta - s/2*(1e0-betaf*cosf);
72 u = mta*mta - s/2*(1e0+betaf*cosf);
73 for(iz = 0;iz<2;iz++)
74 {
75 for(L1=1;L1<=2;L1++)
76 {
77 for(L2=1;L2<=2;L2++)
78 {
79 for(L3=1;L3<=2;L3++)
80 {
81 for(L4=1;L4<=2;L4++)
82 {
83 if(flav==1) downdown_(&L1,&L2,&L3,&L4,&s,&t,&u,&iz,&har,&hai);
84 else upup_(&L1,&L2,&L3,&L4,&s,&t,&u,&iz,&har,&hai);
85 amp[L1-1][L2-1][L3-1][L4-1][iz].im=hai;
86 amp[L1-1][L2-1][L3-1][L4-1][iz].re=har;
87 }
88 }
89 }
90 }
91 }
92 for(i=0;i<4;i++){
93 for(j=0;j<4;j++){
94 for(iz = 0;iz<2;iz++)
95 {
96 sum = 0e0;
97 for(L1=1;L1<=2;L1++)
98 {
99 for(L2=1;L2<=2;L2++)
100 {
101 for(L3=1;L3<=2;L3++)
102 {
103 for(L4=1;L4<=2;L4++)
104 for(M3=1;M3<=2;M3++)
105 for(M4=1;M4<=2;M4++)
106 {
107 c=mul(mul(mul(amp[L1-1][L2-1][L3-1][L4-1][iz],conju(amp[L1-1][L2-1][M3-1][M4-1][iz])),
108 sigma[i][L3-1][M3-1]),(sigma[j][M4-1][L4-1]));
109 sum=sum+ c.re;
110 }
111 }
112 }
113 }
114 sig[i][j][iz]=sum;
115 }
116 R[i][j] = conhc * // to pbarn
117 nc/fc*1e0/2.0/s *
118 1e0/4 * // spin sum
119 (sig[i][j][1] - sig[i][j][0]) * // |Amp|^2 - linearized
120 betaf/16/pi; // phase_space/dcos{theta}
121 if(i==0 && j==0)
122 divv=R[0][0];
123
124 R[i][j]=R[i][j]/divv;
125 if(i!=0) R[i][j] =- R[i][j]; // necessary for tau- but further investigation needed
126
127 }
128 }
129
130 /* // there should be rotation like this one, however we do not need that.
131 // Wrong definition of Pauli matrix for antiparticle, missing conjugation?
132 for(j=0;j<4;j++)
133 {
134 xx= R[1][j];
135 R[3][j] =- R[3][j];
136 R[1][j] = -xx;
137 }
138
139 */
140
141
142 for(i=0;i<4;i++) // rotations around z axis must be checked
143 {
144 xx= R[i][1];
145 R[i][1] =- R[i][2];
146 R[i][2] = xx;
147 }
148
149
150 for(j=0;j<4;j++) // rotations around z axis must be checked
151 {
152 xx= R[1][j];
153 R[1][j] = R[2][j];
154 R[2][j] =- xx;
155 }
156
157
158 /*
159 printf("\n");
160 printf("d_sigma/d_cos{theta} = %.9f \n",divv);
161 printf("\n");
162 printf("R(i,j)= \n");
163 printf("\n");
164 for(i=0;i<4;i++)
165 {
166 for(j=0;j<4;j++)
167 printf("%.5f ",R[i][j]);
168 printf("\n");
169 }
170
171 sinf=sqrt(1-cosf*cosf);
172 MM=sqrt(4e0*mta*mta/s);
173 xnorm=1+cosf*cosf + MM*MM*sinf*sinf;
174 Bench[0][0]=(1+cosf*cosf + MM*MM*sinf*sinf)/xnorm;
175 Bench[1][1]=(-(1- MM*MM)*sinf*sinf)/xnorm;
176 Bench[2][2]=( (1+ MM*MM)*sinf*sinf)/xnorm;
177 Bench[3][3]=(1+cosf*cosf - MM*MM*sinf*sinf)/xnorm;
178 Bench[2][3]=(2*MM*sinf*cosf)/xnorm;
179 Bench[3][2]=(2*MM*sinf*cosf)/xnorm;
180 printf("\n");
181 printf("Bench(i,j)= \n");
182 printf("\n");
183
184 for(i=0;i<4;i++)
185 {
186 for(j=0;j<4;j++)
187 printf("%.5f ",Bench[i][j]);
188 printf("\n");
189 }
190
191 printf("\n");
192 printf(" \n");
193 printf("Bench(i,j)=R(i,j) at low energies, where ew corr and Z contribute little \n");
194 printf("\n");
195 */
196
197}
198
199int main(int argc, char **argv){
200 double sloop=100;
201 double costhetloop=0.0;
202 int i,j,k,l;
203 int NS1=100, NS2=100, NS3=100, NCOS=21;
204 double smin1=log(6*6);
205 double smax1=log(17000*17000);
206
207 double smin2=85*85;
208 double smax2=110*110;
209
210 double smin3=160*160;
211 double smax3=220*220;
212 int flav;
213 double step;
214 FILE *f=NULL,*d=NULL,*dd=NULL;
215 unsigned char c;
216 char ftime[255],path[255];
217 time_t utime;
218 if (argc > 1) flav=atoi(argv[1]);
219 else flav=1;
220 if (flav==1) f=fopen("table1-1.txt","w");
221 else f=fopen("table2-2.txt","w");
222
223 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
224 printconsts_(&buf);
225 dd=fopen("SancLib_v1_02/lib.txt","r");
226 if(!dd) { printf("No info file on electrowek library SANC (file SancLib_v1_02/lib.txt missing), we better stop \n"); return -1; }
227 d=fopen("var.dump","r");
228 if(!d) { printf("No initialization variables of SANC (file var.dump missing), we better stop \n"); return -1; }
229
230 fprintf(f,"Dimensions: %i %i %i %i\n",NS1,NS2,NS3,NCOS);
231 fprintf(f,"Ranges: %.5f %.5f %.5f %.5f %.5f %.5f \n",smin1,smax1,smin2,smax2,smin3,smax3);
232 getcwd(path,255);
233 time(&utime);
234 strftime(ftime,255,"%d %b %Y, %H:%M:%S, GMT%z",localtime(&utime));
235 fprintf(f,"Timestamp: %s\n",ftime);
236 fprintf(f,"Path: %s\n\n",path);
237
238 while(1)
239 {
240 c=0;
241 fscanf(dd,"%c",&c);
242 if(c==0) break;
243 fprintf(f,"%c",c);
244 }
245
246 while(1)
247 {
248 c=0;
249 fscanf(d,"%c",&c);
250 if(c==0) break;
251 fprintf(f,"%c",c);
252 }
253 fprintf(f,"\nBeginRange1\n");
254 printf("\nBeginRange1\n");
255
256 step=(smax1-smin1)/(NS1-1);
257 for(i=0;i<NS1;i++)
258 {
259 sloop=exp(smin1+i*step);
260 printf("%.5f \n",sloop);
261 for(j=0;j<NCOS;j++)
262 {
263 costhetloop=-1.+1.0/NCOS+j*2./(NCOS);
264
265 iew=0;
266 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
267 Rcalc(flav,sloop,costhetloop);
268 borndivv=divv;
269
270 iew=1;
271 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
272 Rcalc(flav,sloop,costhetloop);
273
274 for(k=0;k<4;k++)
275 for(l=0;l<4;l++)
276 fprintf(f,"%.5f ",R[k][l]);
277 fprintf(f,"%.8f ",divv);
278 fprintf(f,"%.8f\n",borndivv);
279 // fprintf(f,"%.5f \n",sloop);
280 }
281 }
282
283
284 fprintf(f,"\nBeginRange2\n");
285 printf("\nBeginRange2\n");
286
287 step=(smax2-smin2)/(NS2-1);
288 for(i=0;i<NS2;i++)
289 {
290 sloop=(smin2+i*step);
291 printf("%.5f \n",sloop);
292 for(j=0;j<NCOS;j++)
293 {
294 costhetloop=-1.+1.0/NCOS+j*2./(NCOS);
295
296 iew=0;
297 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
298 Rcalc(flav,sloop,costhetloop);
299 borndivv=divv;
300
301 iew=1;
302 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
303 Rcalc(flav,sloop,costhetloop);
304
305 for(k=0;k<4;k++)
306 for(l=0;l<4;l++)
307 fprintf(f,"%.5f ",R[k][l]);
308 fprintf(f,"%.5f ",divv);
309 fprintf(f,"%.5f \n",borndivv);
310 }
311 }
312
313
314 fprintf(f,"\nBeginRange3\n");
315 printf("\nBeginRange3\n");
316
317 step=(smax3-smin3)/(NS3-1);
318 for(i=0;i<NS3;i++)
319 {
320 sloop=(smin3+i*step);
321 printf("%.5f \n",sloop);
322 for(j=0;j<NCOS;j++)
323 {
324 costhetloop=-1.+1.0/NCOS+j*2./(NCOS);
325
326 iew=0;
327 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
328 Rcalc(flav,sloop,costhetloop);
329 borndivv=divv;
330
331 iew=1;
332 flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
333 Rcalc(flav,sloop,costhetloop);
334
335 for(k=0;k<4;k++)
336 for(l=0;l<4;l++)
337 fprintf(f,"%.5f ",R[k][l]);
338 fprintf(f,"%.5f ",divv);
339 fprintf(f,"%.5f \n",borndivv);
340 }
341 }
342 fprintf(f,"End\n");
343 fclose(f);
344 return 0;
345}