C++InterfacetoTauola
SANCtable.cxx
1 #include "SANCtable.h"
2 #include "unistd.h"
3 
4 int SANCtable::ns1=0,SANCtable::ns2=0,SANCtable::ns3=0,SANCtable::ncos=0;
5 double SANCtable::smin1=0,SANCtable::smax1=0;
6 double SANCtable::smin2=0,SANCtable::smax2=0;
7 double SANCtable::smin3=0,SANCtable::smax3=0;
8 int SANCtable::iqed=0,SANCtable::iew=0,SANCtable::iborn=0;
9 int SANCtable::gfscheme=0,SANCtable::ifgg=0;
10 double SANCtable::nc=0,SANCtable::fc=0,SANCtable::tlmu2=0;
11 
12 //-----------------------------------------------------------------------------
13 //-Static part-----------------------------------------------------------------
14 //-----------------------------------------------------------------------------
15 void SANCtable::setDimensions(int n1, int n2, int n3, int nc)
16 {
17  ns1=n1;
18  ns2=n2;
19  ns3=n3;
20  ncos=nc;
21 }
22 void 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 }
31 void 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 //-----------------------------------------------------------------------------
42 void 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 
56 bool 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 
71 bool 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 
84 void 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 
122 void SANCtable::open(const char *name)
123 {
124  f.open(name);
125  isOpen = f.is_open();
126 }
127 
128 void SANCtable::close()
129 {
130  f<<"End"<<endl;
131  f.close();
132  isOpen=false;
133 }
134 //-----------------------------------------------------------------------------
135 //-The main computation module-------------------------------------------------
136 //-----------------------------------------------------------------------------
137 double 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 }
Generated on Tue Mar 12 2024 22:57:50 for C++InterfacetoTauola by   doxygen 1.8.7