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