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