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"
8using namespace std;
9
10namespace 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
14TauSpinner::EWborn ewbornMu;
15TauSpinner::EWborn ewbornDown;
16TauSpinner::EWborn ewbornUp;
17
18int m_initTables=0; // status of EW tables initialization
19int m_status=-1; // status of other Born parameters initialization
20
21double m_AMZi; // memorized variables of Born initialization
22double m_GAM;
23double m_SWeff;
24double m_alfinv;
25double m_DeltSQ;
26double m_DeltV;
27double m_Gmu;
28int m_keyGSW;
29
30int ID0=-1; // memorized variables for check if Born
31double S0=-5; // kinematics or EW variant changed
32double cost0=-2;
33int key0=-5;
34
35
36double 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
46 return m_initTables;
47 }
48
49// reads in tables with electroweak formfactors.
50int 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.
94int 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;
116cout << "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)
355complex<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
424double 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
438double 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
452double 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
466double 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// Three 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)
526void 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
543// 2) reads in users initialization of parameters to the local storage and change m_status to 0 (parameters re-initialized)
544void ExtraEWparamsSet( double AMZi, double GAM, double SWeff, double alfinv,double DeltSQ, double DeltV, double Gmu,int keyGSW){
545 if(m_AMZi !=AMZi ){ m_AMZi =AMZi; m_status=0;} // AMZi= Amz(ID); // use for initialization from headers of EW tables
546 if(m_GAM !=GAM ){ m_GAM =GAM; m_status=0;} // GAM=Gamz(ID); // use for initialization from headers of EW tables
547 if(m_SWeff !=SWeff ){ m_SWeff =SWeff; m_status=0;} // SWeff=sin2W(ID); // use for initialization from headers of EW tables
548 if(m_alfinv!=alfinv){ m_alfinv=alfinv; m_status=0;}
549 if(m_DeltSQ!=DeltSQ){ m_DeltSQ=DeltSQ; m_status=0;}
550 if(m_DeltV !=DeltV ){ m_DeltV =DeltV; m_status=0;}
551 if(m_Gmu !=Gmu ){ m_Gmu =Gmu; m_status=0;}
552 if(m_keyGSW!=keyGSW){ m_keyGSW=keyGSW; m_status=0;}
553
554}
555
556// 3) returns m_status for the info or for use.
558 return m_status;
559}
560
561// routine initializes parameters and form-factors for t_bornew_
562// whenever it is necessary, otherwise inactive
563int initEWff(int ID,double S,double cost,int key){
564 int keyGSW;
565 //exit(-1);
566 if ( CheckinitTables()==0){
567 cout <<" electroweak tables not initialized. We stop run "<< endl;
568 exit(-1);
569 return 0;}
570
571 if (ID==ID0 && key==key0 && S==S0 && cost==cost0 && ExtraEWparams()!=0){ // all is set no action needed
572 return 1;}
573
574 else if (CheckinitTables()==1){ // initialization itself
575
576 int finID = 11;
577 // if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
578
579 if (ExtraEWparams()==-1){
580 AMZi= Amz(ID); // initialization from headers of EW tables //91.18870000;
581 GAM=Gamz(ID); // initialization from headers of EW tables
582 SWeff=sin2W(ID); // initialization from headers of EW tables // 0.2235200000;// 0.2121517; //0.22352;// 0.22351946; //0.231708; //.231; // dummy
583 DeltSQ=0;
584 DeltV=0;
585 //double Gmu=0.00001166389;// 0.00001166378; //1.16639e-5;
586 //double alfinv=137.0359895;// dummy
587 //int keyGSW=1;
588 }
589 else if (ExtraEWparams()==0 || ID0!=ID || S0!=S || cost0!=cost || key0!=key ){
590 ID0=ID;
591 S0=S;
592 cost0=cost;
593 key0=key;
594 ExtraEWparamsGet(&AMZi,&GAM,&SWeff,&alfinv,&DeltSQ,&DeltV,&Gmu,&keyGSW);
595 }
596
597 double ReGSW1 = 1.0;
598 double ReGSW2 = 1.0;
599 double ReGSW3 = 1.0;
600 double ReGSW4 = 1.0;
601 double ReGSW6 = 1.0;
602
603 double ImGSW1 = 0.0;
604 double ImGSW2 = 0.0;
605 double ImGSW3 = 0.0;
606 double ImGSW4 = 0.0;
607 double ImGSW6 = 0.0;
608
609 double SS=S;
610 double costhe=cost;
611 if(keyGSW==1){ //options
612
613 ReGSW1 = real(EWFACT(ID,0,SS,costhe)); // reGSW[1];
614 ReGSW2 = real(EWFACT(ID,1,SS,costhe)); // reGSW[2];
615 ReGSW3 = real(EWFACT(ID,2,SS,costhe)); // reGSW[3];
616 ReGSW4 = real(EWFACT(ID,3,SS,costhe)); // reGSW[4];
617 // ReGSW6 = real(EWFACT(ID,5,SS,costhe)); // reGSW[6]; //part only
618 ReGSW6 = real(EWFACT(ID,6,SS,costhe)); // reGSW[6];
619
620 ImGSW1 = imag(EWFACT(ID,0,SS,costhe)); //imGSW[1];
621 ImGSW2 = imag(EWFACT(ID,1,SS,costhe)); //imGSW[2];
622 ImGSW3 = imag(EWFACT(ID,2,SS,costhe)); //imGSW[3];
623 ImGSW4 = imag(EWFACT(ID,3,SS,costhe)); //imGSW[4];
624 // ImGSW6 = imag(EWFACT(ID,5,SS,costhe)); //imGSW[6]; //part only
625 ImGSW6 = imag(EWFACT(ID,6,SS,costhe)); //imGSW[6];
626 }
627 else if(keyGSW==3){
628 ReGSW1 = real(EWFACT(ID,0,SS,costhe)); // reGSW[1];
629 ReGSW2 = 1.0;
630 ReGSW3 = 1.0;
631 ReGSW4 = 1.0;
632 ReGSW6 = real(EWFACT(ID,6,SS,costhe)); // reGSW[6];
633
634 ImGSW1 = imag(EWFACT(ID,0,SS,costhe)); //imGSW[1];
635 ImGSW2 = 0.0;
636 ImGSW3 = 0.0;
637 ImGSW4 = 0.0;
638 ImGSW6 = imag(EWFACT(ID,6,SS,costhe)); //imGSW[6];
639
640 }
641 else if(keyGSW==4){
642 ReGSW1 = 1.005000;
643 ReGSW2 = 1.0;
644 ReGSW3 = 1.0;
645 ReGSW4 = 1.0;
646 ReGSW6 = 1.0;
647
648 ImGSW1 = 0.0;
649 ImGSW2 = 0.0;
650 ImGSW3 = 0.0;
651 ImGSW4 = 0.0;
652 ImGSW6 = 0.0;
653
654 }
655
656 if (ID>0) {
657 int IT=2; // ZbW 12.Nov.2019 flipped to 2 , previously it was 1 ?
658 if(ID==1) IT=1;
659 int mode=1; //dummy parameter?
660 initwkswdelt_(&mode, &IT, &finID, &SS, &SWeff, &DeltSQ, &DeltV, &Gmu, &alfinv, &AMZi, &GAM, &keyGSW, &ReGSW1, &ImGSW1, &ReGSW2, &ImGSW2, &ReGSW3, &ImGSW3, &ReGSW4, &ImGSW4, &ReGSW6, &ImGSW6 );
661 }
662 else if (ID==0) { // it is NOT expected to be used with PDFs it is for e+e- clean beams only!
663 int IT=11; // ZbW 12.Nov.2019 flipped to 2 , previously it was 1 ?
664 //if(ID==1) IT=1;
665 int mode=1; //dummy parameter?
666 initwkswdelt_(&mode, &IT, &finID, &SS, &SWeff, &DeltSQ, &DeltV, &Gmu, &alfinv, &AMZi, &GAM, &keyGSW, &ReGSW1, &ImGSW1, &ReGSW2, &ImGSW2, &ReGSW3, &ImGSW3, &ReGSW4, &ImGSW4, &ReGSW6, &ImGSW6 );
667 }
668 else
669 {
670 ID = -ID;
671 int IT=2; // ZbW 12.Nov.2019 flipped to 2 , previously it was 1 ?
672 if(ID==1) IT=1;
673 int mode=1; //dummy?
674 initwkswdelt_(&mode, &IT, &finID, &SS, &SWeff, &DeltSQ, &DeltV, &Gmu, &alfinv, &AMZi, &GAM, &keyGSW, &ReGSW1, &ImGSW1, &ReGSW2, &ImGSW2, &ReGSW3, &ImGSW3, &ReGSW4, &ImGSW4, &ReGSW6, &ImGSW6 );
675 }
676
677
678 }
679 return 1;
680}
681
682
683 //
684 // Calculates Born cross-section summed over final taus spins.
685 // Input parameters:
686 // incoming flavour ID
687 // invariant mass^2 SS
688 // scattering angle costhe
689 // effective weingber SWeff
690 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){
691
692 // BORN x-section.
693 // WARNING: overall sign of costheta may need confirmation
694
695 //if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
696
697 int key=1;
698
699 double AMZi=AMZ0; //91.1876; //91.18870000;
700 double GAM=GAM0; //2.495378;//2.49520000;
701
702 if (mode==1){
703 AMZi= Amz(ID);GAM=Gamz(ID);SWeff=sin2W(ID);// initialization from headers of EW tables
704 }
705 ExtraEWparamsSet(AMZi, GAM, SWeff, alfinv,DeltSQ, DeltV, Gmu,keyGSW);
706 initEWff(ID,SS,costhe, key);
707
708
709 double dOne = 1.0;
710 double dMOne = -1.0;
711
712 // this is for running width
713
714 // sum DY Born over all tau helicity configurations:
715 // EW loop corrections only implemented in t_born
716
717 int Bmode=1; // it is about mass terms
718 return ( t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dOne) + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dMOne)
719 + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dOne) + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dMOne))/alfinv/alfinv;
720
721 // 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.
722 return 1;
723}
724
725 //
726 // Calculates Born cross-section summed over final taus spins.
727 // Input parameters:
728 // incoming flavour ID
729 // invariant mass^2 SS
730 // scattering angle costhe
731 // effective weingber SWeff
732 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){
733
734 // BORN x-section.
735 // WARNING: overall sign of costheta may need confirmation
736
737 // if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
738
739 int key=1;
740
741 double AMZi=AMZ0; // 91.1876;//91.18870000;
742 double GAM=GAM0; //2.495378;// 2.49520000;
743 if (mode==1){
744 AMZi= Amz(ID);GAM=Gamz(ID);SWeff=sin2W(ID);// initialization from headers of EW tables
745 }
746 ExtraEWparamsSet(AMZi, GAM, SWeff, alfinv,DeltSQ, DeltV, Gmu,keyGSW);
747 initEWff(ID,SS,costhe, key);
748
749
750 double dOne = 1.0;
751 double dMOne = -1.0;
752
753 // this is for running width
754
755 // sum DY Born over all tau helicity configurations:
756 // EW loop corrections only implemented in t_born
757
758 int Bmode=1; // it is about mass terms
759 return ( -t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dOne) - t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dOne , &dMOne)
760 + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dOne) + t_bornew_(&Bmode, &keyGSW, &SS, &costhe, &dMOne, &dMOne))/alfinv/alfinv;
761
762 // 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.
763 return 1;
764}
765
766}
void ExtraEWparamsGet(double *AMZi, double *GAM, double *SWeff, double *alfinv, double *DeltSQ, double *DeltV, double *Gmu, int *keyGSW)
Definition: EWtables.cxx:526
double sin2W(int FLAV)
Definition: EWtables.cxx:452
int initTables(char *mumu, char *downdown, char *upup)
Definition: EWtables.cxx:50
int CheckinitTables()
Definition: EWtables.cxx:45
double Amz(int FLAV)
Definition: EWtables.cxx:424
double Gamz(int FLAV)
Definition: EWtables.cxx:438
int ExtraEWparams()
Definition: EWtables.cxx:557
void ExtraEWparamsSet(double AMZi, double GAM, double SWeff, double alfinv, double DeltSQ, double DeltV, double Gmu, int keyGSW)
Definition: EWtables.cxx:544
complex< double > EWFACT(int FLAV, int NO, double s, double costhe)
Definition: EWtables.cxx:355
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:732
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:690
int initEWff(int ID, double S, double cost, int key)
Definition: EWtables.cxx:563
double QCDFACT(int FLAV, int NO, double s)
Definition: EWtables.cxx:466