StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SANCinterface.c
1 #include <math.h>
2 #include <stdio.h>
3 #include <time.h>
4 
5 extern void upup_(int *L1,int *L2, int *L3,int *L4,double *s,double *t,double *u,int *iz,double *har,double *hai);
6 extern void downdown_(int *L1,int *L2, int *L3,int *L4,double *s,double *t,double *u,int *iz,double *har,double *hai);
7 extern void flagset_(int *iqedx,int *iewx,int *ibornx,int *gfschemex,int *ifggx,double *ncx,double *fcx,double *tlmu2x);
8 extern void paraget_(double *mtax, double *conhcx, double *pix);
9 extern void printconsts_(int *pmg);
10 
11 
12 int iqed=0,iew=1,iborn=0,gfscheme=0,ifgg=1;
13 double nc=1.0,fc=3.0,tlmu2=1e-5;
14 int buf=0;
15 double R[4][4];
16 double divv;
17 double borndivv=123.456;
18 
19 struct compl {double im,re;};
20 
21 struct 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 
29 struct 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 
37 double abso(struct compl x)
38 {
39  return sqrt(x.re*x.re+x.im*x.im);
40 }
41 
42 
43 void 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 
199 int 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 }
Generated by   doxygen 1.8.5