C++InterfacetoTauola
EWtables.cxx
1 #include <iostream>
2 #include <stdlib.h>
3 #include <fstream>
4 #include <complex>
5 #include <string>
6 #include "TauSpinner/ew_born.h"
7 #include "TauSpinner/EWtables.h"
8 using namespace std;
9 
10 namespace TauSpinner {
11 
12 // for routines migrating from table-parsing-test.cxx. Finally they will resettle to library. May be form own class as well.
13 
14 TauSpinner::EWborn ewbornMu;
15 TauSpinner::EWborn ewbornDown;
16 TauSpinner::EWborn ewbornUp;
17 
18 int m_initTables=0; // status of EW tables initialization
19 int m_status=-1; // status of other Born parameters initialization
20 
21 double m_AMZi; // memorized variables of Born initialization
22 double m_GAM;
23 double m_SWeff;
24 double m_alfinv;
25 double m_DeltSQ;
26 double m_DeltV;
27 double m_Gmu;
28 int m_keyGSW;
29 
30 int ID0=-1; // memorized variables for check if Born
31 double S0=-5; // kinematics or EW variant changed
32 double cost0=-2;
33 int key0=-5;
34 
35 
36 double AMZi= 91.18870000; //memorized variables of initEWff
37  double GAM=2.49520000;// 2.498111432484;
38  double SWeff= 0.2235200000;// 0.2121517; //0.22352;// 0.22351946; //0.231708; //.231; // dummy
39  double DeltSQ=0;
40  double DeltV=0;
41  double Gmu=0.00001166389;// 0.00001166378; //1.16639e-5;
42  double alfinv=137.0359895;// dummy
43 
44 // provides info flag if tables were initialized
45 int CheckinitTables(){
46  return m_initTables;
47  }
48 
49 // reads in tables with electroweak formfactors.
50 int initTables(char* mumu, char* downdown, char* upup) {
51  // char* mumu="table.mu";
52  //char* downdown= "table.down";
53  //char* upup= "table.up";
54 
55  const char *tableLocationMu = mumu;
56  const char *tableLocationDown = downdown;
57  const char *tableLocationUp = upup;
58 
59 
60 
61 
62  cout << "initializing table " << tableLocationMu << ": " << endl;
63  bool resultMu = ewbornMu.FillFromTable(tableLocationMu);
64  if (!resultMu)
65  {
66  cout << "ERROR: could not parse table : " << tableLocationMu << endl;
67  return -1;
68  }
69 
70  cout << "initializing table " << tableLocationDown << ": " << endl;
71  bool resultDown = ewbornDown.FillFromTable(tableLocationDown);
72  if (!resultDown)
73  {
74  cout << "ERROR: could not parse table: " << tableLocationDown << endl;
75  return -1;
76  }
77 
78  cout << "initializing table " << tableLocationUp << ": " << endl;
79  bool resultUp = ewbornUp.FillFromTable(tableLocationUp);
80  if (!resultUp)
81  {
82  cout << "ERROR: could not parse table: " << tableLocationUp << endl;
83  return -1;
84  }
85  m_initTables=1;
86  return 1;
87 }
88 
89  /* Print out some values */
90 
91 /*
92 // ELEMENATY and EARLY test, of little value now
93 //prints stuff on electroweak tables, privides checks if they are properly read.
94 int testit() {
95 
96  //const char *tableLocationMu = "table.mu";
97  //const char *tableLocationDown = "table.down";
98  //const char *tableLocationUp = "table.up";
99  //TauSpinner::EWborn ewbornMu;
100  //TauSpinner::EWborn ewbornDown;
101  //TauSpinner::EWborn ewbornUp;
102 
103 
104  cout.precision(8);
105  cout.setf(std::ios::fixed);
106 
107  cout << " " << endl;
108  cout << " ==Print out some values== " << endl;
109  cout << " " << endl;
110 
111 
112 
113  cout << " " << endl;
114  cout << " ==ewbornMu== " << endl;
115  cout << " " << endl;
116 cout << "HEADER : "
117  << ewbornMu.MZ << " "
118  << ewbornMu.MH << " "
119  << ewbornMu.MT << " "
120  << ewbornMu.SWSQ << " "
121  << ewbornMu.GZ << " "
122  << ewbornMu.MW << " "
123  << ewbornMu.GW << endl << endl;
124  cout << "ranges for mu : "<< endl;
125  cout << " a: = " << ewbornMu.EEa[0] <<" to " << ewbornMu.EEa[ewbornMu.NA-1] << endl;
126  cout << " b: = " << ewbornMu.EEb[0] <<" to " << ewbornMu.EEb[ewbornMu.NB-1] << endl;
127  cout << " c: = " << ewbornMu.EEc[0] <<" to " << ewbornMu.EEc[ewbornMu.NC-1] << endl;
128  cout << " d: = " << ewbornMu.EEd[0] <<" to " << ewbornMu.EEd[ewbornMu.ND-1] << endl;
129  cout << "steps for mu at start : "<< endl;
130  cout << " a: = " << ewbornMu.EEa[1] - ewbornMu.EEa[0] << endl;
131  cout << " b: = " << ewbornMu.EEb[1] - ewbornMu.EEb[0] << endl;
132  cout << " c: = " << ewbornMu.EEc[1] - ewbornMu.EEc[0] << endl;
133  cout << " d: = " << ewbornMu.EEd[1] - ewbornMu.EEd[0] << endl;
134  cout << " " << endl;
135  cout << "steps for mu at higer points : "<< endl;
136  cout << " a: = " << ewbornMu.EEa[5] - ewbornMu.EEa[4] << endl;
137  cout << " b: = " << ewbornMu.EEb[5] - ewbornMu.EEb[4] << endl;
138  cout << " c: = " << ewbornMu.EEc[5] - ewbornMu.EEc[4] << endl;
139  cout << " d: = " << ewbornMu.EEd[5] - ewbornMu.EEd[4] << endl;
140 
141 
142  cout << "random prints : "<< endl;
143  cout << "EEa[100] = " << ewbornMu.EEa[100] << endl;
144  cout << "FFa[100][6] = " << ewbornMu.FFa[100][6] << endl;
145  cout << "FSa[100][1] = " << ewbornMu.FSa[100][1] << endl;
146  cout << endl;
147  cout << "EEb[120][14] = " << ewbornMu.EEb[120] << endl;
148  cout << "FFb[120][14][6] = " << ewbornMu.FFb[120][14][6] << endl;
149  cout << "FSb[120][1] = " << ewbornMu.FSb[120][1] << endl;
150  cout << endl;
151  cout << "EEc[145][30] = " << ewbornMu.EEc[145] << endl;
152  cout << "FFc[145][30][6] = " << ewbornMu.FFc[145][30][6] << endl;
153  cout << "FSc[145][1] = " << ewbornMu.FSc[145][1] << endl;
154  cout << "COSc[30] = " << ewbornMu.COSc[30] << endl;
155  cout << endl;
156  cout << "EEd[ 75][14] = " << ewbornMu.EEd[75] << endl;
157  cout << "FFd[ 75][14][6] = " << ewbornMu.FFd[75][14][6] << endl;
158  cout << "FSd[ 75][1] = " << ewbornMu.FSd[75][1] << endl;
159  cout << endl;
160  cout << "EEd[ 80][14] = " << ewbornMu.EEd[80] << endl;
161  cout << "FFd[ 80][14][6] = " << ewbornMu.FFd[80][14][6] << endl;
162  cout << "FSd[ 80][1] = " << ewbornMu.FSd[80][1] << endl;
163  cout << "COSd[ 0] = " << ewbornMu.COSd[0] << endl;
164  cout << "COSd[14] = " << ewbornMu.COSd[14] << endl;
165 
166 
167 
168  cout << " " << endl;
169  cout << " ==ewbornDown== " << endl;
170  cout << " " << endl;
171  cout << "HEADER : "
172  << ewbornDown.MZ << " "
173  << ewbornDown.MH << " "
174  << ewbornDown.MT << " "
175  << ewbornDown.SWSQ << " "
176  << ewbornDown.GZ << " "
177  << ewbornDown.MW << " "
178  << ewbornDown.GW << endl << endl;
179 
180  cout << "EEa[100] = " << ewbornDown.EEa[100] << endl;
181  cout << "FFa[100][6] = " << ewbornDown.FFa[100][6] << endl;
182  cout << "FSa[100][1] = " << ewbornDown.FSa[100][1] << endl;
183  cout << endl;
184  cout << "EEb[120][14] = " << ewbornDown.EEb[120] << endl;
185  cout << "FFb[120][14][6] = " << ewbornDown.FFb[120][14][6] << endl;
186  cout << "FSb[120][1] = " << ewbornDown.FSb[120][1] << endl;
187  cout << endl;
188  cout << "EEc[145][30] = " << ewbornDown.EEc[145] << endl;
189  cout << "FFc[145][30][6] = " << ewbornDown.FFc[145][30][6] << endl;
190  cout << "FSc[145][1] = " << ewbornDown.FSc[145][1] << endl;
191  cout << "COSc[30] = " << ewbornDown.COSc[30] << endl;
192  cout << endl;
193  cout << "EEd[ 75][14] = " << ewbornDown.EEd[75] << endl;
194  cout << "FFd[ 75][14][6] = " << ewbornDown.FFd[75][14][6] << endl;
195  cout << "FSd[ 75][1] = " << ewbornDown.FSd[75][1] << endl;
196  cout << endl;
197  cout << "EEd[ 80][14] = " << ewbornDown.EEd[80] << endl;
198  cout << "FFd[ 80][14][6] = " << ewbornDown.FFd[80][14][6] << endl;
199  cout << "FSd[ 80][1] = " << ewbornDown.FSd[80][1] << endl;
200  cout << "COSd[ 0] = " << ewbornDown.COSd[0] << endl;
201  cout << "COSd[14] = " << ewbornDown.COSd[14] << endl;
202  cout << " " << endl;
203  cout << " == Table print Sector B for Down tables of formfactors == " << endl;
204  cout << " " << endl;
205  for (int ii=0; ii<5;ii++){
206  int i=60+ii;
207  printf(" * i= %d * \n",i );
208  for (int j=0; j<6;j++){
209  printf("C_%2d:",2*j );
210  for (int k=0; k<4;k++){
211  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornDown.FFb[i][2*j][k]),imag(ewbornDown.FFb[i][2*j][k]));
212  }
213  for (int k=5; k<7;k++){
214  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornDown.FFb[i][2*j][k]),imag(ewbornDown.FFb[i][2*j][k]));
215  }
216 
217  printf(" \n");
218  }
219  }
220  cout << " " << endl;
221  cout << " == Table print Sector C for Down tables of formfactors == " << endl;
222  cout << " " << endl;
223  for (int ii=0; ii<5;ii++){
224  int i=ii*20;
225  printf(" * i= %d * \n",i );
226  for (int j=0; j<6;j++){
227  printf("C_%2d:",2*j );
228  for (int k=0; k<4;k++){
229  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornDown.FFc[i][2*j][k]),imag(ewbornDown.FFc[i][2*j][k]));
230  }
231  for (int k=5; k<7;k++){
232  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornDown.FFc[i][2*j][k]),imag(ewbornDown.FFc[i][2*j][k]));
233  }
234 
235  printf(" \n");
236  }
237  }
238  cout << " " << endl;
239  cout << " == Table print Sector D for Down tables of formfactors == " << endl;
240  cout << " " << endl;
241  for (int ii=0; ii<5;ii++){
242  int i=ii*20;
243  printf(" * i= %d * \n",i );
244  for (int j=0; j<6;j++){
245  printf("C_%2d:",2*j );
246  for (int k=0; k<4;k++){
247  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornDown.FFd[i][2*j][k]),imag(ewbornDown.FFd[i][2*j][k]));
248  }
249  for (int k=5; k<7;k++){
250  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornDown.FFd[i][2*j][k]),imag(ewbornDown.FFd[i][2*j][k]));
251  }
252 
253  printf(" \n");
254  }
255  }
256 
257  cout << " ========== " << endl;
258 
259  cout << " " << endl;
260  cout << " ==ewbornUp== " << endl;
261  cout << " " << endl;
262  cout << "HEADER : "
263  << ewbornUp.MZ << " "
264  << ewbornUp.MH << " "
265  << ewbornUp.MT << " "
266  << ewbornUp.SWSQ << " "
267  << ewbornUp.GZ << " "
268  << ewbornUp.MW << " "
269  << ewbornUp.GW << endl << endl;
270 
271  cout << "EEa[100] = " << ewbornUp.EEa[100] << endl;
272  cout << "FFa[100][6] = " << ewbornUp.FFa[100][6] << endl;
273  cout << "FSa[100][1] = " << ewbornUp.FSa[100][1] << endl;
274  cout << endl;
275  cout << "EEb[120][14] = " << ewbornUp.EEb[120] << endl;
276  cout << "FFb[120][14][6] = " << ewbornUp.FFb[120][14][6] << endl;
277  cout << "FSb[120][1] = " << ewbornUp.FSb[120][1] << endl;
278  cout << endl;
279  cout << "EEc[145][30] = " << ewbornUp.EEc[145] << endl;
280  cout << "FFc[145][30][6] = " << ewbornUp.FFc[145][30][6] << endl;
281  cout << "FSc[145][1] = " << ewbornUp.FSc[145][1] << endl;
282  cout << "COSc[30] = " << ewbornUp.COSc[30] << endl;
283  cout << endl;
284  cout << "EEd[ 75][14] = " << ewbornUp.EEd[75] << endl;
285  cout << "FFd[ 75][14][6] = " << ewbornUp.FFd[75][14][6] << endl;
286  cout << "FSd[ 75][1] = " << ewbornUp.FSd[75][1] << endl;
287  cout << endl;
288  cout << "EEd[ 80][14] = " << ewbornUp.EEd[80] << endl;
289  cout << "FFd[ 80][14][6] = " << ewbornUp.FFd[80][14][6] << endl;
290  cout << "FSd[ 80][1] = " << ewbornUp.FSd[80][1] << endl;
291  cout << "COSd[ 0] = " << ewbornUp.COSd[0] << endl;
292  cout << "COSd[14] = " << ewbornUp.COSd[14] << endl;
293 
294  cout << " " << endl;
295  cout << " == Table print Sector B for Up tables of formfactors == " << endl;
296  cout << " " << endl;
297  for (int ii=0; ii<5;ii++){
298  int i=ii*20;
299  printf(" * i= %d * \n",i );
300  for (int j=0; j<6;j++){
301  printf("C_%2d:",2*j );
302  for (int k=0; k<4;k++){
303  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornUp.FFb[i][2*j][k]),imag(ewbornUp.FFb[i][2*j][k]));
304  }
305  for (int k=5; k<7;k++){
306  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornUp.FFb[i][2*j][k]),imag(ewbornUp.FFb[i][2*j][k]));
307  }
308 
309  printf(" \n");
310  }
311  }
312  cout << " " << endl;
313  cout << " == Table print Sector C for Up tables of formfactors == " << endl;
314  cout << " " << endl;
315  for (int ii=0; ii<5;ii++){
316  int i=ii*20;
317  printf(" * i= %d * \n",i );
318  for (int j=0; j<6;j++){
319  printf("C_%2d:",2*j );
320  for (int k=0; k<4;k++){
321  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornUp.FFc[i][2*j][k]),imag(ewbornUp.FFc[i][2*j][k]));
322  }
323  for (int k=5; k<7;k++){
324  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornUp.FFc[i][2*j][k]),imag(ewbornUp.FFc[i][2*j][k]));
325  }
326 
327  printf(" \n");
328  }
329  }
330  cout << " " << endl;
331  cout << " == Table print Sector D for Up tables of formfactors == " << endl;
332  cout << " " << endl;
333  for (int ii=0; ii<5;ii++){
334  int i=ii*20;
335  printf(" * i= %d * \n",i );
336  for (int j=0; j<6;j++){
337  printf("C_%2d:",2*j );
338  for (int k=0; k<4;k++){
339  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornUp.FFd[i][2*j][k]),imag(ewbornUp.FFd[i][2*j][k]));
340  }
341  for (int k=5; k<7;k++){
342  printf(" F%d=(%12.5e,%12.5e), ", k,real(ewbornUp.FFd[i][2*j][k]),imag(ewbornUp.FFd[i][2*j][k]));
343  }
344 
345  printf(" \n");
346  }
347  }
348 
349 
350  return 1;
351 }
352 */
353 
354 // provides electroweak form-factor for input of (int FLAV, int NO, double s, double costhe)
355 complex<double> EWFACT(int FLAV, int NO, double s, double costhe){
356  complex<double> rezu=complex<double>(0.3,0.4);
357  double Ene=sqrt(s);
358  // cout << "ranges for mu and Ene=: "<< Ene <<endl;
359  TauSpinner::EWborn *ewb = NULL;
360 
361  // choice of flavour, pointer fixing should be better
362  // than copyuing of all struct
363  if (FLAV==0)
364  ewb=&ewbornMu;
365  else if (FLAV==2) //{
366  ewb=&ewbornUp; //cout << "######### up ########## "<< endl;}
367  else if (FLAV==1) //{
368  ewb=&ewbornDown; //cout << "######### down ########## "<< endl;}
369 
370  // ranges for tables
371  double EaMin= ewb->EEa[0]; double EaMax= ewb->EEa[ewb->NA-1];
372  double EbMin= ewb->EEb[0]; double EbMax= ewb->EEb[ewb->NB-1];
373  double EcMin= ewb->EEc[0]; double EcMax= ewb->EEc[ewb->NC-1];
374  double EdMin= ewb->EEd[0]; double EdMax= ewb->EEd[ewb->ND-1];
375 
376  // step for tables
377  double stepA= log(ewb->EEa[1] /ewb->EEa[0]);
378  double stepB= ewb->EEb[1] - ewb->EEb[0];
379  double stepC= ewb->EEc[1] - ewb->EEc[0];
380  double stepD= ewb->EEd[1] - ewb->EEd[0];
381  // cout << "######### stepD "<< stepD<<" "<< ewb->EEd[1]<<" "<< ewb->EEd[0] <<endl;
382  if(Ene>EbMin && Ene<EbMax){ // case B
383  int Index=((Ene-EbMin)/(EbMax-EbMin)*(ewb->NB-1));
384  int Ic=(((costhe+1.0)/2*(ewb->MB-1)));
385  rezu= ewb->FFb[Index][Ic][NO]
386  + (Ene -ewb->EEb[Index]) /stepB * (ewb->FFb[Index+1][Ic][NO] - ewb->FFb[Index][Ic][NO])
387  + ( (costhe+1.)/2. - (1.0*Ic)/(ewb->MB-1) )*(ewb->MB-1)*(
388  (1- (Ene -ewb->EEb[Index])/stepB )*(ewb->FFb[Index ][Ic+1][NO] - ewb->FFb[Index ][Ic][NO])
389  +( (Ene -ewb->EEb[Index])/stepB )*(ewb->FFb[Index+1][Ic+1][NO] - ewb->FFb[Index+1][Ic][NO]) );
390  }
391  else if (Ene>EcMin && Ene<EcMax){ // case C
392  int Index=((Ene-EcMin)/(EcMax-EcMin)*(ewb->NC-1));
393  int Ic=(((costhe+1.0)/2*(ewb->MC-1)));
394  rezu= ewb->FFc[Index][Ic][NO]
395  + (Ene -ewb->EEc[Index]) /stepC * (ewb->FFc[Index+1][Ic][NO] - ewb->FFc[Index][Ic][NO])
396  + ( (costhe+1.)/2. - (1.0*Ic)/(ewb->MC-1) )*(ewb->MC-1)*(
397  (1- (Ene -ewb->EEc[Index])/stepC )*(ewb->FFc[Index ][Ic+1][NO] - ewb->FFc[Index ][Ic][NO])
398  +( (Ene -ewb->EEc[Index])/stepC )*(ewb->FFc[Index+1][Ic+1][NO] - ewb->FFc[Index+1][Ic][NO]) );
399  }
400  else if (Ene>EdMin && Ene<EdMax){ // case D
401  int Index=((Ene-EdMin)/(EdMax-EdMin)*(ewb->ND-1));
402  int Ic=((costhe+1.0)/2.)*(ewb->MD-1);
403  rezu= ewb->FFd[Index][Ic][NO]
404  + (Ene -ewb->EEd[Index]) /stepD * (ewb->FFd[Index+1][Ic][NO] - ewb->FFd[Index][Ic][NO])
405  + ( (costhe+1.)/2. - (1.0*Ic)/(ewb->MD-1) )*(ewb->MD-1)*(
406  (1- (Ene -ewb->EEd[Index])/stepD )*(ewb->FFd[Index ][Ic+1][NO] - ewb->FFd[Index ][Ic][NO])
407  +( (Ene -ewb->EEd[Index])/stepD )*(ewb->FFd[Index+1][Ic+1][NO] - ewb->FFd[Index+1][Ic][NO]) );
408 
409 
410  }
411  else if (Ene>EaMin && Ene<EaMax){ // case A (logarithmic steps)
412  int Index=((log(Ene)-log(EaMin))/(log(EaMax)-log(EaMin))*(ewb->NA-1));
413  rezu= ewb->FFa[Index][NO]
414  + (log(Ene)-log(ewb->EEa[Index])) /stepA * (ewb->FFa[Index+1][NO] - ewb->FFa[Index][NO]);
415  // cout << " strefa A, Ene= "<<Ene<<endl;
416  }
417  else{
418  rezu=complex<double>(1.0,0.0);
419  }
420  return rezu;
421 }
422 
423 // returns Z mass as stored in header of electroweak table for FLAV
424 double Amz(int FLAV){
425  TauSpinner::EWborn *ewb = NULL;
426  if (FLAV==0)
427  ewb=&ewbornMu;
428  else if (FLAV==2)
429  ewb=&ewbornUp;
430  else if (FLAV==1)
431  ewb=&ewbornDown;
432 
433  // Amz= ewb->MZ;
434  return ewb->MZ;
435 }
436 
437 // returns Z widt as stored in header of electroweak table for FLAV
438 double Gamz(int FLAV){
439  TauSpinner::EWborn *ewb = NULL;
440  if (FLAV==0)
441  ewb=&ewbornMu;
442  else if (FLAV==2)
443  ewb=&ewbornUp;
444  else if (FLAV==1)
445  ewb=&ewbornDown;
446 
447  // Amz= ewb->MZ;
448  return ewb->GZ;
449 }
450 
451 // returns sin^2theta_W^eff as stored in header of electroweak table for FLAV
452 double sin2W(int FLAV){
453  TauSpinner::EWborn *ewb = NULL;
454  if (FLAV==0)
455  ewb=&ewbornMu;
456  else if (FLAV==2)
457  ewb=&ewbornUp;
458  else if (FLAV==1)
459  ewb=&ewbornDown;
460 
461  // Amz= ewb->MZ;
462  return ewb->SWSQ;
463 }
464 
465 // returns QCD factor as interpolated from electroweak table for s. For given FLAV and NO
466 double QCDFACT(int FLAV, int NO, double s){
467  double rezu=0.3;
468  double Ene=sqrt(s);
469  // cout << "ranges for mu and Ene=: "<< Ene <<endl;
470  TauSpinner::EWborn *ewb = NULL;
471 
472  // choice of flavour, pointer fixing should be better
473  // than copyuing of all struct
474  if (FLAV==0)
475  ewb=&ewbornMu;
476  else if (FLAV==2)
477  ewb=&ewbornUp;
478  else if (FLAV==1)
479  ewb=&ewbornDown;
480 
481  // ranges for tables
482  double EaMin= ewb->EEa[0]; double EaMax= ewb->EEa[ewb->NA-1];
483  double EbMin= ewb->EEb[0]; double EbMax= ewb->EEb[ewb->NB-1];
484  double EcMin= ewb->EEc[0]; double EcMax= ewb->EEc[ewb->NC-1];
485  double EdMin= ewb->EEd[0]; double EdMax= ewb->EEd[ewb->ND-1];
486 
487  // step for tables
488  double stepA= log(ewb->EEa[1] /ewb->EEa[0]);
489  double stepB= ewb->EEb[1] - ewb->EEb[0];
490  double stepC= ewb->EEc[1] - ewb->EEc[0];
491  double stepD= ewb->EEd[1] - ewb->EEd[0];
492 
493  if(Ene>EbMin && Ene<EbMax){ // case B
494  int Index=((Ene-EbMin)/(EbMax-EbMin)*(ewb->NB-1));
495  rezu= ewb->FSb[Index][NO]
496  + (Ene-ewb->EEb[Index]) /stepB * (ewb->FSb[Index+1][NO] - ewb->FSb[Index][NO]);
497 
498  }
499  else if (Ene>EcMin && Ene<EcMax){ // case C
500  int Index=((Ene-EcMin)/(EcMax-EcMin)*(ewb->NC-1));
501  rezu= ewb->FSc[Index][NO]
502  + (Ene-ewb->EEc[Index]) /stepC * (ewb->FSc[Index+1][NO] - ewb->FSc[Index][NO]);
503 
504 
505  }
506  else if (Ene>EdMin && Ene<EdMax){ // case D
507  int Index=((Ene-EdMin)/(EdMax-EdMin)*(ewb->ND-1));
508  rezu= ewb->FSd[Index][NO]
509  + (Ene-ewb->EEd[Index]) /stepD * (ewb->FSd[Index+1][NO] - ewb->FSd[Index][NO]);
510 
511  }
512  else if (Ene>EaMin && Ene<EaMax){ // case A (logarithmic steps)
513  int Index=((log(Ene)-log(EaMin))/(log(EaMax)-log(EaMin))*(ewb->NA-1));
514  rezu= ewb->FSa[Index][NO]
515  + (log(Ene)-log(ewb->EEa[Index])) /stepA * (ewb->FSa[Index+1][NO] - ewb->FSa[Index][NO]);
516  }
517  else{
518  rezu=0.12;
519  }
520  return rezu;
521 }
522 
523 
524 // Routines to manipulate parameters for user variant for Born.
525 // 1) takes parameters from the local storage and change m_status to 1 (parameters taken)
526 void ExtraEWparamsGet( double *AMZi, double *GAM, double *SWeff, double *alfinv, double *DeltSQ, double* DeltV, double *Gmu,int *keyGSW){
527  if (m_status==-1){
528  cout <<"ERROR Born parameters are not initialized: you can not get them. We stop "<< endl;
529  exit(-1);
530  }
531  *AMZi =m_AMZi;
532  *GAM =m_GAM;
533  *SWeff =m_SWeff;
534  *alfinv=m_alfinv;
535  *DeltSQ=m_DeltSQ;
536  *DeltV =m_DeltV;
537  *Gmu =m_Gmu;
538  *keyGSW=m_keyGSW;
539  m_status=1;
540  // cout << " params get alfinv= "<<*alfinv<<" GAM="<<*GAM<<" keyGSW="<<*keyGSW<< endl;
541 }
542 void keyGSWGet(int *keyGSW){
543  if (m_status==-1){
544  cout <<"WARNING: Born parameters are not initialized you can not adjust keyGSW "<< endl;
545  }
546  else {
547  *keyGSW=m_keyGSW;
548  }
549  // cout << " params get alfinv= "<<*alfinv<<" GAM="<<*GAM<<" keyGSW="<<*keyGSW<< endl;
550 }
551 
552 
553 
554 // 2) reads in users initialization of parameters to the local storage and change m_status to 0 (parameters re-initialized)
555 void ExtraEWparamsSet( double AMZi, double GAM, double SWeff, double alfinv,double DeltSQ, double DeltV, double Gmu,int keyGSW){
556  if(m_AMZi !=AMZi ){ m_AMZi =AMZi; m_status=0;} // AMZi= Amz(ID); // use for initialization from headers of EW tables
557  if(m_GAM !=GAM ){ m_GAM =GAM; m_status=0;} // GAM=Gamz(ID); // use for initialization from headers of EW tables
558  if(m_SWeff !=SWeff ){ m_SWeff =SWeff; m_status=0;} // SWeff=sin2W(ID); // use for initialization from headers of EW tables
559  if(m_alfinv!=alfinv){ m_alfinv=alfinv; m_status=0;}
560  if(m_DeltSQ!=DeltSQ){ m_DeltSQ=DeltSQ; m_status=0;}
561  if(m_DeltV !=DeltV ){ m_DeltV =DeltV; m_status=0;}
562  if(m_Gmu !=Gmu ){ m_Gmu =Gmu; m_status=0;}
563  if(m_keyGSW!=keyGSW){ m_keyGSW=keyGSW; m_status=0;}
564 
565 }
566 
567 // 3) returns m_status for the info or for use.
568 int ExtraEWparams(){
569  return m_status;
570 }
571 
572 // routine initializes parameters and form-factors for t_bornew_
573 // whenever it is necessary, otherwise inactive
574 int initEWff(int ID,double S,double cost,int key){
575  int keyGSW;
576  //exit(-1);
577  if ( CheckinitTables()==0){
578  cout <<" electroweak tables not initialized. We stop run "<< endl;
579  exit(-1);
580  return 0;}
581 
582  if (ID==ID0 && key==key0 && S==S0 && cost==cost0 && ExtraEWparams()!=0){ // all is set no action needed
583  return 1;}
584 
585  else if (CheckinitTables()==1){ // initialization itself
586 
587  int finID = 11;
588  // if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
589 
590  if (ExtraEWparams()==-1){
591  AMZi= Amz(ID); // initialization from headers of EW tables //91.18870000;
592  GAM=Gamz(ID); // initialization from headers of EW tables
593  SWeff=sin2W(ID); // initialization from headers of EW tables // 0.2235200000;// 0.2121517; //0.22352;// 0.22351946; //0.231708; //.231; // dummy
594  DeltSQ=0;
595  DeltV=0;
596  //double Gmu=0.00001166389;// 0.00001166378; //1.16639e-5;
597  //double alfinv=137.0359895;// dummy
598  //int keyGSW=1;
599  }
600  else if (ExtraEWparams()==0 || ID0!=ID || S0!=S || cost0!=cost || key0!=key ){
601  ID0=ID;
602  S0=S;
603  cost0=cost;
604  key0=key;
605  ExtraEWparamsGet(&AMZi,&GAM,&SWeff,&alfinv,&DeltSQ,&DeltV,&Gmu,&keyGSW);
606  }
607  // std::cout << "keyGSW" << keyGSW << std::endl;
608 
609 
610  double ReGSW1 = 1.0;
611  double ReGSW2 = 1.0;
612  double ReGSW3 = 1.0;
613  double ReGSW4 = 1.0;
614  double ReGSW6 = 1.0;
615 
616  double ImGSW1 = 0.0;
617  double ImGSW2 = 0.0;
618  double ImGSW3 = 0.0;
619  double ImGSW4 = 0.0;
620  double ImGSW6 = 0.0;
621 
622  double SS=S;
623  double SSS=S;
624  double costhe=cost;
625  if(keyGSW!=1){costhe=0.0;SSS=AMZi*AMZi;}
626  if(keyGSW==1){ //options
627 
628  ReGSW1 = real(EWFACT(ID,0,SS,costhe)); // reGSW[1];
629  ReGSW2 = real(EWFACT(ID,1,SS,costhe)); // reGSW[2];
630  ReGSW3 = real(EWFACT(ID,2,SS,costhe)); // reGSW[3];
631  ReGSW4 = real(EWFACT(ID,3,SS,costhe)); // reGSW[4];
632  // ReGSW6 = real(EWFACT(ID,5,SS,costhe)); // reGSW[6]; //part only
633  ReGSW6 = real(EWFACT(ID,6,SS,costhe)); // reGSW[6];
634 
635  ImGSW1 = imag(EWFACT(ID,0,SS,costhe)); //imGSW[1];
636  ImGSW2 = imag(EWFACT(ID,1,SS,costhe)); //imGSW[2];
637  ImGSW3 = imag(EWFACT(ID,2,SS,costhe)); //imGSW[3];
638  ImGSW4 = imag(EWFACT(ID,3,SS,costhe)); //imGSW[4];
639  // ImGSW6 = imag(EWFACT(ID,5,SS,costhe)); //imGSW[6]; //part only
640  ImGSW6 = imag(EWFACT(ID,6,SS,costhe)); //imGSW[6];
641  }
642  else if(keyGSW==3){
643  ReGSW1 = 1.0; // reGSW[1];
644  ReGSW2 = 1.0;
645  ReGSW3 = 1.0;
646  ReGSW4 = 1.0;
647  ReGSW6 = real(EWFACT(ID,6,SSS,costhe)); // reGSW[6];
648 
649  ImGSW1 = 0.0; //imGSW[1];
650  ImGSW2 = 0.0;
651  ImGSW3 = 0.0;
652  ImGSW4 = 0.0;
653  ImGSW6 = imag(EWFACT(ID,6,SSS,costhe)); //imGSW[6];
654 
655  }
656  else if(keyGSW==4){
657  ReGSW1 = real(EWFACT(ID,0,SSS,costhe)); // reGSW[1];
658  ReGSW2 = 1.0;
659  ReGSW3 = 1.0;
660  ReGSW4 = 1.0;
661  ReGSW6 = real(EWFACT(ID,6,SSS,costhe)); // reGSW[6];
662 
663  ImGSW1 = imag(EWFACT(ID,0,SSS,costhe)); //imGSW[1];
664  ImGSW2 = 0.0;
665  ImGSW3 = 0.0;
666  ImGSW4 = 0.0;
667  ImGSW6 = imag(EWFACT(ID,6,SSS,costhe)); //imGSW[6];
668 
669  }
670  else if(keyGSW==5){
671  ReGSW1 = real(EWFACT(ID,0,SSS,costhe)); // reGSW[1];
672  ReGSW2 = real(EWFACT(ID,1,SSS,costhe)); // reGSW[2];
673  ReGSW3 = real(EWFACT(ID,2,SSS,costhe)); // reGSW[3];
674  ReGSW4 = real(EWFACT(ID,3,SSS,costhe)); //1.0;
675  ReGSW6 = real(EWFACT(ID,6,SSS,costhe)); // reGSW[6];
676 
677  ImGSW1 = imag(EWFACT(ID,0,SSS,costhe)); //imGSW[1];
678  ImGSW2 = imag(EWFACT(ID,1,SSS,costhe)); //imGSW[2];
679  ImGSW3 = imag(EWFACT(ID,2,SSS,costhe)); //imGSW[3];
680  ImGSW4 = imag(EWFACT(ID,3,SSS,costhe)); //0.0;
681  ImGSW6 = imag(EWFACT(ID,6,SSS,costhe)); //imGSW[6];
682 
683  /*
684  // real v, vv
685  ImGSW2 = 0.0;
686  ImGSW3 = 0.0;
687  ImGSW4 = 0.0;
688  */
689  /*
690  // real v, vv=1
691  ReGSW4 = 1.0;
692  ImGSW2 = 0.0;
693  ImGSW3 = 0.0;
694  ImGSW4 = 0.0;
695  */
696  /*
697  // real pigamgam
698  ReGSW4 = 1.0;
699  ImGSW1 = 0.0;
700  ImGSW2 = 0.0;
701  ImGSW3 = 0.0;
702  ImGSW4 = 0.0;
703  ImGSW6 = 0.0; // make Pi_gammgamma real: formally it is higher order. But
704  // see bardin+grunewald+pasarino-9902452.pdf
705  // why it is inconsistent with the LEP choice
706  */
707  /*
708  // LEP2005 style
709  ReGSW4 = 1.0;
710  ImGSW1 = 0.0;
711  ImGSW2 = 0.0;
712  ImGSW3 = 0.0;
713  ImGSW4 = 0.0;
714  */
715  }
716  else if(keyGSW==12){
717  ReGSW1 = real(EWFACT(0,0,SSS,costhe)); // for (v1)
718  }
719  else if(keyGSW==13){
720  ReGSW1 =real(EWFACT(ID,0,SSS,costhe)); // for (v2)
721  ReGSW2 =real(EWFACT(ID,1,SSS,costhe))/real(EWFACT(0,1,SSS,costhe)); // for (v2)
722  ReGSW3 =real(EWFACT(ID,2,SSS,costhe))/real(EWFACT(0,2,SSS,costhe)); // for (v2)
723  }
724 
725  if (ID>0) {
726  int IT=2; // ZbW 12.Nov.2019 flipped to 2 , previously it was 1 ?
727  if(ID==1) IT=1;
728  int mode=1; //dummy parameter?
729  initwkswdelt_(&mode, &IT, &finID, &SS, &SWeff, &DeltSQ, &DeltV, &Gmu, &alfinv, &AMZi, &GAM, &keyGSW, &ReGSW1, &ImGSW1, &ReGSW2, &ImGSW2, &ReGSW3, &ImGSW3, &ReGSW4, &ImGSW4, &ReGSW6, &ImGSW6 );
730  }
731  else if (ID==0) { // it is NOT expected to be used with PDFs it is for e+e- clean beams only!
732  int IT=11; // ZbW 12.Nov.2019 flipped to 2 , previously it was 1 ?
733  //if(ID==1) IT=1;
734  int mode=1; //dummy parameter?
735  initwkswdelt_(&mode, &IT, &finID, &SS, &SWeff, &DeltSQ, &DeltV, &Gmu, &alfinv, &AMZi, &GAM, &keyGSW, &ReGSW1, &ImGSW1, &ReGSW2, &ImGSW2, &ReGSW3, &ImGSW3, &ReGSW4, &ImGSW4, &ReGSW6, &ImGSW6 );
736  }
737  else
738  {
739  ID = -ID;
740  int IT=2; // ZbW 12.Nov.2019 flipped to 2 , previously it was 1 ?
741  if(ID==1) IT=1;
742  int mode=1; //dummy?
743  initwkswdelt_(&mode, &IT, &finID, &SS, &SWeff, &DeltSQ, &DeltV, &Gmu, &alfinv, &AMZi, &GAM, &keyGSW, &ReGSW1, &ImGSW1, &ReGSW2, &ImGSW2, &ReGSW3, &ImGSW3, &ReGSW4, &ImGSW4, &ReGSW6, &ImGSW6 );
744  }
745 
746 
747  }
748  return 1;
749 }
750 
751 
752  //
753  // Calculates Born cross-section summed over final taus spins.
754  // Input parameters:
755  // incoming flavour ID
756  // invariant mass^2 SS
757  // scattering angle costhe
758  // effective weingber SWeff
759  double sigbornswdelt(int mode, int ID, double SS, double costhe, double SWeff, double DeltSQ, double DeltV, double Gmu, double alfinv, double AMZ0, double GAM0, int keyGSW){
760 
761  // BORN x-section.
762  // WARNING: overall sign of costheta may need confirmation
763 
764  //if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
765 
766  int key=1;
767 
768  double AMZi=AMZ0; //91.1876; //91.18870000;
769  double GAM=GAM0; //2.495378;//2.49520000;
770 
771  if (mode==1){
772  AMZi= Amz(ID);GAM=Gamz(ID);SWeff=sin2W(ID);// initialization from headers of EW tables
773  }
774  ExtraEWparamsSet(AMZi, GAM, SWeff, alfinv,DeltSQ, DeltV, Gmu,keyGSW);
775  initEWff(ID,SS,costhe, key);
776 
777 
778  double dOne = 1.0;
779  double dMOne = -1.0;
780 
781  // this is for running width
782  // sum DY Born over all tau helicity configurations:
783  // EW loop corrections only implemented in t_born
784  int Bmode=1; // it is about mass terms
785  return ( t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dOne) + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dMOne)
786  + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dOne) + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dMOne)); // /alfinv/alfinv;
787 
788  // overall norm. factor .../SS/123231 most probably it is alpha_QED**2/pi/2/SS is from comparison between Born we use and Born used in Warsaw group.
789  return 1;
790 }
791 
792  //
793  // Calculates Born cross-section summed over final taus spins.
794  // Input parameters:
795  // incoming flavour ID
796  // invariant mass^2 SS
797  // scattering angle costhe
798  // effective weingber SWeff
799  double AsNbornswdelt(int mode, int ID, double SS, double costhe, double SWeff, double DeltSQ, double DeltV, double Gmu, double alfinv, double AMZ0, double GAM0, int keyGSW){
800 
801  // BORN x-section.
802  // WARNING: overall sign of costheta may need confirmation
803 
804  // if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
805 
806  int key=1;
807 
808  double AMZi=AMZ0; // 91.1876;//91.18870000;
809  double GAM=GAM0; //2.495378;// 2.49520000;
810  if (mode==1){
811  AMZi= Amz(ID);GAM=Gamz(ID);SWeff=sin2W(ID);// initialization from headers of EW tables
812  }
813  ExtraEWparamsSet(AMZi, GAM, SWeff, alfinv,DeltSQ, DeltV, Gmu,keyGSW);
814  initEWff(ID,SS,costhe, key);
815 
816 
817  double dOne = 1.0;
818  double dMOne = -1.0;
819 
820  // this is for running width
821 
822  // sum DY Born over all tau helicity configurations:
823  // EW loop corrections only implemented in t_born
824 
825  int Bmode=1; // it is about mass terms
826  return ( -t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dOne) - t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dMOne)
827  + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dOne) + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dMOne)); // /alfinv/alfinv;
828 
829  // overall norm. factor .../SS/123231 most probably it is alpha_QED**2/pi/2/SS is from comparison between Born we use and Born used in Warsaw group.
830  return 1;
831 }
832 
833 }
TauSpinner::Amz
double Amz(int FLAV)
Definition: EWtables.cxx:424
TauSpinner::ExtraEWparamsSet
void ExtraEWparamsSet(double AMZi, double GAM, double SWeff, double alfinv, double DeltSQ, double DeltV, double Gmu, int keyGSW)
Definition: EWtables.cxx:555
TauSpinner::ExtraEWparams
int ExtraEWparams()
Definition: EWtables.cxx:568
TauSpinner::EWFACT
complex< double > EWFACT(int FLAV, int NO, double s, double costhe)
Definition: EWtables.cxx:355
TauSpinner::Gamz
double Gamz(int FLAV)
Definition: EWtables.cxx:438
TauSpinner::sigbornswdelt
double sigbornswdelt(int mode, int ID, double SS, double costhe, double SWeff, double DeltSQ, double DeltV, double Gmu, double alfinv, double AMZ0, double GAM0, int keyGSW)
Definition: EWtables.cxx:759
TauSpinner::initTables
int initTables(char *mumu, char *downdown, char *upup)
Definition: EWtables.cxx:50
TauSpinner::initEWff
int initEWff(int ID, double S, double cost, int key)
Definition: EWtables.cxx:574
TauSpinner::AsNbornswdelt
double AsNbornswdelt(int mode, int ID, double SS, double costhe, double SWeff, double DeltSQ, double DeltV, double Gmu, double alfinv, double AMZ0, double GAM0, int keyGSW)
Definition: EWtables.cxx:799
TauSpinner::ExtraEWparamsGet
void ExtraEWparamsGet(double *AMZi, double *GAM, double *SWeff, double *alfinv, double *DeltSQ, double *DeltV, double *Gmu, int *keyGSW)
Definition: EWtables.cxx:526
TauSpinner::keyGSWGet
void keyGSWGet(int *keyGSW)
Definition: EWtables.cxx:542
TauSpinner::CheckinitTables
int CheckinitTables()
Definition: EWtables.cxx:45
TauSpinner::QCDFACT
double QCDFACT(int FLAV, int NO, double s)
Definition: EWtables.cxx:466
TauSpinner::sin2W
double sin2W(int FLAV)
Definition: EWtables.cxx:452
Generated on Tue Mar 12 2024 22:57:50 for C++InterfacetoTauola by   doxygen 1.8.7