C++ Interface to Tauola
vbftests.cxx
1 #include "TauSpinner/tau_reweight_lib.h"
2 #include "TauSpinner/Particle.h"
3 #include "TauSpinner/vbftests.h"
4 #include <cstdlib>
5 #include <string.h>
6 using namespace Tauolapp;
7 
8 /*This is a group of methods prepared for the testing purpose. The tests include those where fixed kinematics is
9 forced. Then external calculation can be used for matrix element calculation and/or PDFs etc.
10 Gradually better documentation will be introduced. It will be documented with comment lines and in the draft of the paper.
11 To use one has to modify the Makefile the appropriate line of
12 objects has to be activated (# removed). Then invoking lines in code have to be commented out as well.
13 At a certain moment such commented out lines will be removed from the code or will remain. We will have to make up the mind.
14 -*/
15 
16 namespace TauSpinner {
17 
18 extern double CMSENE;
19 extern int relWTnonSM;
20 extern double WTnonSM;
21 extern double Polari;
22 extern double WTamplit;
23 extern double WTamplitP;
24 extern double WTamplitM;
25 extern double CMSENE;
26 extern double f(double x, int ID, double SS, double cmsene);
27 
28 const bool DEBUG = 0;
29 
30 /** Wrapper to VBDISTR */
31 double vbfdistrWRAP(int I1, int I2, int I3, int I4, int H1, int H2, double P[6][4], int KEY)
32 {
33  double P_copy[6][4];
34  memcpy(P_copy,P,sizeof(double)*6*4);
35 
36  return vbfdistr_(&I1, &I2, &I3, &I4, &H1, &H2, P_copy, &KEY);
37 
38 }
39 
40 //---------------------------------------------------------------------------------------------------------------------------------------------------
41  void calcXsect(int IDPROD, SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY)
42 {
43  /*
44  // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation!
45  // FSR photons may need to be treated explicitely or with interpolation procedures.
46  Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(),
47  p3.py()+p4.py()+tau1.py()+tau2.py(),
48  p3.pz()+p4.pz()+tau1.pz()+tau2.pz(),
49  p3.e() +p4.e() +tau1.e() +tau2.e(), 0 );
50  */
51  Particle P_QQ( p3.px()+p4.px()+sp_X.px(),
52  p3.py()+p4.py()+sp_X.py(),
53  p3.pz()+p4.pz()+sp_X.pz(),
54  p3.e() +p4.e() +sp_X.e() , 0 );
55 
56 
57  double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
58 
59  double x1x2 = SS/CMSENE/CMSENE;
60  double x1Mx2 = P_QQ.pz()/CMSENE*2;
61 
62  double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
63  double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
64 
65  //---------------------------------------------------------------------------
66  // Construct the matrix for FORTRAN function
67  // NOTE: different order of indices than in FORTRAN!
68  double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 },
69  { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 },
70  { p3.e(), p3.px(), p3.py(), p3.pz() },
71  { p4.e(), p4.px(), p4.py(), p4.pz() },
72  { tau1.e(), tau1.px(), tau1.py(), tau1.pz() },
73  { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } };
74 
75 
76  W[0][0]=0.0;
77  W[0][1]=0.0;
78  W[1][0]=0.0;
79  W[1][1]=0.0;
80 
81  //
82  // Calculate 'f' function for all x1 and all ID1, ID2
83  //
84  double f_x1_ARRAY[11] = { 0.0 };
85  double f_x2_ARRAY[11] = { 0.0 };
86  double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with
87  double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5
88 
89  double SSfix = 91.188*91.188;
90  double CMSENEfix = 91.188;
91  for(int i=-5;i<=5;++i) {
92  f_x1[i] = f(x1,i,SSfix,CMSENEfix);
93  f_x2[i] = f(x2,i,SSfix,CMSENEfix);
94  // std::cout << "ERW: x, Q2, f(x,Q2), f(x,Q)=" << x1 << " " << SS << " " << f(x1,i,SS,CMSENE) << " " << f(x1,i,sqrt(SS),CMSENE) << std::endl;
95  }
96 
97  W[0][0]=0.0;
98  W[0][1]=0.0;
99  W[1][0]=0.0;
100  W[1][1]=0.0;
101 
102  // these loops need to be cleaned from zero contributions!
103  for(int I1=-4;I1<=4;I1++){ // for test of single production process fix flavour
104  for(int I2=-4;I2<=4;I2++){ // for test of single production process fix flavour
105  for(int I3=-4;I3<=4;I3++){
106  for(int I4=-4;I4<=4;I4++){
107 
108  int ID1 = I1; if( ID1 == 0 ) ID1 = 21;
109  int ID2 = I2; if( ID2 == 0 ) ID2 = 21;
110  int ID3 = I3; if( ID3 == 0 ) ID3 = 21;
111  int ID4 = I4; if( ID4 == 0 ) ID4 = 21;
112 
113  //preselect production group
114  bool accept = false;
115  // any process
116  if( IDPROD == 0) accept = true;
117  // only GG processes
118  if( IDPROD == 1 && ID1 == 21 && ID2 == 21 && abs (ID3) < 5 && abs (ID4) < 5 ) accept = true;
119  // only GG processes ... details
120  if( IDPROD == 110 && ID1 == 21 && ID2 == 21 ) {
121  accept = true;
122  if( abs (ID3) != 2 || abs (ID4) != 2 ) accept = false;
123  }
124  if( IDPROD == 111 && ID1 == 21 && ID2 == 21 ) {
125  accept = true;
126  if( abs (ID3) != 1 || abs (ID4) != 1 ) accept = false;
127  }
128  if( IDPROD == 112 && ID1 == 21 && ID2 == 21 ) {
129  accept = true;
130  if( abs (ID3) != 3 || abs (ID4) != 3 ) accept = false;
131  }
132  if( IDPROD == 113 && ID1 == 21 && ID2 == 21 ) {
133  accept = true;
134  if( abs (ID3) != 4 || abs (ID4) != 4 ) accept = false;
135  }
136  // only GQ processes
137  if( IDPROD == 2 && ID1 == 21 && abs(ID2) < 5 ) accept = true ;
138  if( IDPROD == 2 && ID2 == 21 && abs(ID1) < 5 ) accept = true ;
139  // only GQ processes ... details
140  if( (IDPROD == 210) && ( ID1 == 21) && ( ID2 < 5 ) && (ID2 > 0) ) accept = true;
141  if( (IDPROD == 211) && ( ID1 == 21) && ( ID2 < 5 ) && (ID2 < 0) ) accept = true;
142  if( (IDPROD == 212) && ( ID2 == 21) && ( ID1 < 5 ) && (ID1 > 0) ) accept = true;
143  if( (IDPROD == 213) && ( ID2 == 21) && ( ID1 < 5 ) && (ID1 < 0) ) accept = true;
144  // only QQ, QXQX processes
145  if( IDPROD == 3 && (ID1 * ID2 > 0 ) && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true;
146  // only QQ processes .. details
147  if( IDPROD == 31 && ID1 > 0 && ID2 > 0 && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true;
148  if( IDPROD == 32 && ID1 < 0 && ID2 < 0 && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true;
149  // only QQ processes .. details
150  if( IDPROD == 310 && ID1 > 0 && ID2 > 0 && ID1 < 5 && ID2 < 5 && ID1 == ID2 ) accept = true;
151  if( IDPROD == 311 && ID1 > 0 && ID2 > 0 && ID1 < 5 && ID2 < 5 && ID1 != ID2 ) accept = true;
152  if( IDPROD == 312 && ID1 < 0 && ID2 < 0 && ID1 == ID2 ) accept = true;
153  if( IDPROD == 313 && ID1 < 0 && ID2 < 0 && ID1 != ID2 ) accept = true;
154  // only QQ processes .. details
155  if( IDPROD == 320 && ID1 * ID2 > 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) == abs(ID2) ) accept = true;
156  if( IDPROD == 321 && ID1 * ID2 > 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) != abs(ID2) ) accept = true;
157  // only QQ processes .. details
158  if( IDPROD == 3101 && ID1 == 2 && ID2 == 2 && ID3 == 2 && ID4 == 2 ) accept = true;
159  if( IDPROD == 3102 && ID1 == -2 && ID2 == -2 && ID3 == -2 && ID4 == -2 ) accept = true;
160  if( IDPROD == 3103 && ID1 == 2 && ID2 == 2 ) accept = true;
161  if( IDPROD == 3104 && ID1 == -2 && ID2 == -2 ) accept = true;
162  // only QQX processes
163  if( IDPROD == 4 && (ID1 * ID2 < 0 ) && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true;
164  // only QQX processes .. details
165  if( IDPROD == 420 && ID1 * ID2 < 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) == abs(ID2) ) accept = true;
166  if( IDPROD == 421 && ID1 * ID2 < 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) != abs(ID2) ) accept = true;
167  // only QQX processes, first family in initial/final
168  if( IDPROD == 41 && (ID1 * ID2 < 0 ) && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID3) < 5 && abs(ID4) < 5 ) {
169  accept = true;
170  if( abs (ID1) > 2 || abs (ID2) > 2 || abs (ID3) >2 || abs (ID4) > 2 ) accept = false;
171  }
172  // only QQX processes, second family in initial/final state
173  if( IDPROD == 42 && (ID1 * ID2 < 0 ) && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID3) < 5 && abs(ID4) < 5 ) {
174  accept = true;
175  if( abs (ID1) < 3 || abs (ID2) < 3 || abs (ID3) < 3 || abs (ID4) < 3 ) accept = false;
176  }
177  // only QQX processes, mixed families in initial/final state
178  if( IDPROD == 43 && ( ID1 * ID2 < 0) && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID3) < 5 && abs(ID4) < 5 ) {
179  accept = true;
180  if( abs (ID1) < 3 && abs (ID2) < 3 && abs (ID3) < 3 && abs (ID4) < 3 ) accept = false;
181  if( abs (ID1) > 2 && abs (ID2) > 2 && abs (ID3) > 2 && abs (ID4) > 2 ) accept = false;
182  }
183  // only QQX processes, first family in initial state
184  if( IDPROD == 44 && (ID1 * ID2 < 0 ) && abs(ID1) < 5 && abs(ID2) < 5 ) {
185  accept = true;
186  if( abs (ID1) > 2 || abs (ID2) > 2 ) accept = false;
187  }
188  // only QQX processes ... details
189  if( IDPROD == 4101 && (ID1 * ID2 < 0 ) ) {
190  accept = true;
191  // only U UX --> U UX
192  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 2 || abs (ID4) != 2 ) accept = false;
193  if( ID1 != 2 ) accept = false;
194  }
195  if( IDPROD == 4102 && (ID1 * ID2 < 0 ) ) {
196  accept = true;
197  // only U UX --> U UX
198  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 2 || abs (ID4) != 2 ) accept = false;
199  if( ID1 != -2 ) accept = false;
200  }
201  if( IDPROD == 4103&& (ID1 * ID2 < 0 ) ) {
202  accept = true;
203  // only U UX --> Q QX
204  if( abs (ID1) != 2 || abs (ID2) != 2 ) accept = false;
205  if( ID1 != 2 ) accept = false;
206  }
207  if( IDPROD == 4104&& (ID1 * ID2 < 0 ) ) {
208  accept = true;
209  // only UX U --> Q QX
210  if( abs (ID1) != 2 || abs (ID2) != 2 ) accept = false;
211  if( ID1 != -2 ) accept = false;
212  }
213  if( IDPROD == 4105&& (ID1 * ID2 < 0 ) ) {
214  accept = true;
215  // only D DX --> Q QX
216  if( abs (ID1) != 1 || abs (ID2) != 1 ) accept = false;
217  if( ID1 != 1 ) accept = false;
218  }
219  if( IDPROD == 4106&& (ID1 * ID2 < 0 ) ) {
220  accept = true;
221  // only DX D --> Q QX
222  if( abs (ID1) != 1 || abs (ID2) != 1 ) accept = false;
223  if( ID1 != -1 ) accept = false;
224  }
225  if( IDPROD == 4107&& (ID1 * ID2 < 0 ) ) {
226  accept = true;
227  // only Q QX --> Q QX
228  if( ID1 < 0 || ID2 > 0) accept = false;
229  }
230  if( IDPROD == 4108&& (ID1 * ID2 < 0 ) ) {
231  accept = true;
232  // only QX Q --> Q QX
233  if( ID1 > 0 || ID2 < 0) accept = false;
234  }
235  // only QQX processes ... details
236  if( IDPROD == 410 && (ID1 * ID2 < 0 ) ) {
237  accept = true;
238  // only U UX --> U UX
239  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 2 || abs (ID4) != 2 ) accept = false;
240  }
241  if( IDPROD == 411 && (ID1 * ID2 < 0 ) ) {
242  accept = true;
243  // only U UX --> C CX
244  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 4 || abs (ID4) != 4 ) accept = false;
245  }
246  if( IDPROD == 412 && (ID1 * ID2 < 0 ) ) {
247  accept = true;
248  // only U UX --> S DX
249  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 3 || abs (ID4) != 1 ) accept = false;
250  }
251  if( IDPROD == 413 && (ID1 * ID2 < 0 ) ) {
252  accept = true;
253  // only U UX --> D SX
254  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 1 || abs (ID4) != 3 ) accept = false;
255  }
256  if( IDPROD == 414 && (ID1 * ID2 < 0 ) ) {
257  accept = true;
258  // only U UX --> D DX
259  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 1 || abs (ID4) != 1 ) accept = false;
260  }
261  if( IDPROD == 415 && (ID1 * ID2 < 0 ) ) {
262  accept = true;
263  // only U UX --> S SX
264  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 3 || abs (ID4) != 3 ) accept = false;
265  }
266  if( IDPROD == 416 && (ID1 * ID2 < 0 ) ) {
267  accept = true;
268  // only U UX --> Q QX
269  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
270  }
271  if( IDPROD == 417 && (ID1 * ID2 < 0 ) ) {
272  accept = true;
273  // only U UX --> G G
274  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) == 21 || abs (ID4) == 21 ) accept = false;
275  }
276  if( IDPROD == 418 && (ID1 * ID2 < 0 ) ) {
277  accept = true;
278  // only U UX --> J J
279  if( abs (ID1) != 2 || abs (ID2) != 2 ) accept = false;
280  }
281  if( IDPROD == 419 && (ID1 * ID2 < 0 ) ) {
282  accept = true;
283  // only D DX --> D DX
284  if( abs (ID1) != 1 || abs (ID2) != 1 || abs (ID3) != 1 || abs (ID4) != 1 ) accept = false;
285  }
286  if( IDPROD == 511 && (ID1 * ID2 < 0 ) ) {
287  accept = true;
288  // only D DX, DX D --> G G
289  if( abs (ID1) != 1 || abs (ID2) != 1 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
290  }
291  if( IDPROD == 512 && (ID1 * ID2 < 0 ) ) {
292  accept = true;
293  // only D DX --> G G
294  if( ID1 != 1 || ID2 != -1 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
295  }
296  if( IDPROD == 513 && (ID1 * ID2 < 0 ) ) {
297  accept = true;
298  // only DX D --> G G
299  if( ID1 != -1 || ID2 != 1 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
300  }
301  if( IDPROD == 521 && (ID1 * ID2 < 0 ) ) {
302  accept = true;
303  // only U UX, UX U --> G G
304  if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
305  }
306  if( IDPROD == 522 && (ID1 * ID2 < 0 ) ) {
307  accept = true;
308  // only U UX --> G G
309  if( ID1 != 2 || ID2 != -2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
310  }
311  if( IDPROD == 523 && (ID1 * ID2 < 0 ) ) {
312  accept = true;
313  // only UX U --> G G
314  if( ID1 != -2 || ID2 != 2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
315  }
316  if( IDPROD == 531 && (ID1 * ID2 < 0 ) ) {
317  accept = true;
318  // only Q QX, QX Q --> G G
319  if( abs(ID1) > 4 || abs(ID2) > 4 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
320  }
321  if( IDPROD == 532 && (ID1 * ID2 < 0 ) ) {
322  accept = true;
323  // only Q QX --> G G
324  if( abs(ID1) > 4 || abs(ID2) > 4 || ID1 < 0 || abs(ID3) != 21 || abs(ID4) != 21 ) accept = false;
325  }
326  if( IDPROD == 533 && (ID1 * ID2 < 0 ) ) {
327  accept = true;
328  // only QX Q --> G G
329  if( abs(ID1) > 4 || abs(ID2) > 4 || ID2 < 0 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false;
330  }
331  if( IDPROD == 541 ) {
332  accept = false;
333  // only U SX --> D CX
334  if( ID1 == 2 && ID2 == -3 && ID3 == 1 && ID4 == -4 ) accept = true;
335  if( ID1 == -3 && ID2 == 2 && ID3 == 1 && ID4 == -4 ) accept = true;
336  if( ID1 == 2 && ID2 == -3 && ID3 ==-4 && ID4 == 1 ) accept = true;
337  if( ID1 == -3 && ID2 == 2 && ID3 ==-4 && ID4 == 1 ) accept = true;
338  }
339 
340  if(accept ){
341 
342  W[0][0] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P, KEY); // as in case of nonSM_adopt we change the sign for
343  W[0][1] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P, KEY); // second tau helicity assuming without checking that
344  W[1][0] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, -1, P, KEY); // for VBF quantization frame conventions are the same.
345  W[1][1] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, 1, P, KEY);
346 
347  }
348  }
349  }
350  }
351  }
352 }
353 
354 
355 
356 void makeSimpleTestME2()
357 {
358 
359  //---------------------------------------------------------------------------
360  // Construct the matrix for FORTRAN function
361  // NOTE: different order of indices than in FORTRAN!
362 
363  std::cout << "ERW: ----------------------------- " << std::endl;
364  double P1[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
365  { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
366  { 0.8855009E+02, -0.2210038E+02, 0.4007979E+02, -0.7580437E+02 },
367  { 0.3283248E+03, -0.1038482E+03, -0.3019295E+03, 0.7649385E+02 },
368  { 0.1523663E+03, -0.1058795E+03, -0.9770827E+02, 0.4954769E+02 },
369  { 0.4307588E+03, 0.2318280E+03, 0.3595580E+03, -0.5023717E+02 } };
370  calcTestME2(1, P1);
371 
372  Particle p1 ( 0.0, 0.0, 0.5000000E+03, 0.5000000E+03, 1 );
373  Particle p2 ( 0.0, 0.0, -0.5000000E+03, 0.5000000E+03, -1 );
374  Particle p3 ( -0.2210038E+02, 0.4007979E+02, -0.7580437E+02, 0.8855009E+02, 2 );
375  Particle p4 ( -0.1038482E+03, -0.3019295E+03, 0.7649385E+02, 0.3283248E+03, -2 );
376  Particle tau1 ( -0.1058795E+03, -0.9770827E+02, 0.4954769E+02, 0.1523663E+03, 15 );
377  Particle tau2 ( 0.2318280E+03, 0.3595580E+03, -0.5023717E+02, 0.4307588E+03, -15 );
378 
379  std::cout << "ERW: 4-momenta before boost ----------------------------- " << std::endl;
380 
381  p1.print();
382  p2.print();
383  p3.print();
384  p4.print();
385  tau1.print();
386  tau2.print();
387  double x1=0.03; // x1 for the first parton
388  double xm = 1000.0;
389  double cmsene = 13000.0;
390  double x2=xm/cmsene*xm/cmsene/x1; // x2 adjusted from cms energy (13000), and virtuality of the whole system 1000
391  std::cout << "x1, x2= " << x1 << " " << x2 << std::endl;
392  if (x2>=1.0) std::cout << "ERW: wrong x1, x2 cannot be defined (is above 1); x2=" << x2 << std::endl;
393 
394  Particle P_Z0( tau1.px()+tau2.px(), tau1.py()+tau2.py(), tau1.pz()+tau2.pz(), tau1.e()+tau2.e(), 23 );
395  P_Z0.print();
396  double Q2 = P_Z0.recalculated_mass() * P_Z0.recalculated_mass();
397  std::cout << " pdfs = " << f(x1,0,Q2,cmsene) * f(x2,1,Q2,cmsene) << std::endl;
398 
399  Particle P_QQ( 0.0, 1.0e-10, (x1-x2)*cmsene/2.0, (x1+x2)*cmsene/2.0, 10 );
400  P_QQ.print();
401 
402  p1.boostFromRestFrame(P_QQ);
403  p2.boostFromRestFrame(P_QQ);
404  p3.boostFromRestFrame(P_QQ);
405  p4.boostFromRestFrame(P_QQ);
406  P_Z0.boostFromRestFrame(P_QQ);
407  tau1.boostFromRestFrame(P_QQ);
408  tau2.boostFromRestFrame(P_QQ);
409 
410  std::cout << "ERW: 4-momenta after boost ----------------------------- " << std::endl;
411 
412  p1.print();
413  p2.print();
414  p3.print();
415  p4.print();
416  P_Z0.print();
417  tau1.print();
418  tau2.print();
419 
420  std::cout << "ERW: ----------------------------- " << std::endl;
421  double P2[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
422  { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
423  { 0.1177462E+03, -0.6070512E+02, 0.7123011E+02, 0.7145150E+02 },
424  { 0.3509495E+03, -0.3178775E+02, 0.8393832E+02, 0.3392779E+03 },
425  { 0.3493321E+03, 0.1840069E+03, -0.5152712E+02, -0.2924315E+03 },
426  { 0.1819722E+03, -0.9151401E+02, -0.1036413E+03, -0.1182978E+03 } };
427  calcTestME2(2, P2);
428  std::cout << "ERW: ----------------------------- " << std::endl;
429  double P3[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
430  { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
431  { 0.2586900E+03, 0.1324670E+03, -0.1696171E+03, -0.1435378E+03 },
432  { 0.1084567E+03, -0.5735712E+02, -0.2162482E+02, -0.8947281E+02 },
433  { 0.4005742E+03, -0.1580760E+03, 0.3563160E+03, 0.9223569E+02 },
434  { 0.2322791E+03, 0.8296613E+02, -0.1650741E+03, 0.1407749E+03 } };
435  calcTestME2(3, P3);
436  std::cout << "ERW: ----------------------------- " << std::endl;
437  double P4[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
438  { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
439  { 0.1595700E+03, -0.6917808E+02, -0.1395175E+03, -0.3481123E+02 },
440  { 0.2247758E+03, -0.1360140E+03, 0.1650340E+03, -0.6919641E+02 },
441  { 0.2508802E+03, 0.1447863E+01, 0.2499830E+03, -0.2107335E+02 },
442  { 0.3647740E+03, 0.2037442E+03, -0.2754995E+03, 0.1250810E+03 } };
443  calcTestME2(4, P4);
444 
445 }
446 
447 /*******************************************************************************/
448  void calcTestME2(int iter, double P1[6][4])
449 {
450 
451  double ME2ref, ME2;
452  int ID1, ID2, ID3, ID4;
453 
454  int KEY=0;
455 
456  // case GD -> GD
457  ID1 = 21; ID2 = 1; ID3 = 21; ID4 = 1;
458 
459  if(iter == 1) ME2ref = 1.5961579933867344E-010;
460  if(iter == 2) ME2ref = 3.8749742050329834E-010;
461  if(iter == 3) ME2ref = 5.0434636937545397E-012;
462  if(iter == 4) ME2ref = 7.9077184257060427E-012;
463 
464  ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY)
465  + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY)
466  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY)
467  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY);
468 
469  std::cout << "ERW: ME2 GD->GD (ref) = " << ME2ref << std::endl;
470  std::cout << "ERW: ME2 GD->GD = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl;
471 
472  // case GU -> GU
473  ID1 = 21; ID2 = 2; ID3 = 21; ID4 = 2;
474 
475  if(iter == 1) ME2ref = 2.9195503763051040E-010;
476 
477  ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY)
478  + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY)
479  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY)
480  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY);
481 
482  if(iter == 1)std::cout << "ERW: ME2 GU->GU (ref) = " << ME2ref << std::endl;
483  if(iter == 1)std::cout << "ERW: ME2 GU->GU = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl;
484 
485  // case DD -> DD
486  ID1 = 1; ID2 = 1; ID3 = 1; ID4 = 1;
487 
488  if(iter == 1) ME2ref = 3.3953129762581284E-017;
489  if(iter == 2) ME2ref = 6.0054072781075002E-017;
490  if(iter == 3) ME2ref = 3.9707833415682912E-018;
491  if(iter == 4) ME2ref = 1.5342177192807347E-018;
492 
493  ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY)
494  + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY)
495  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY)
496  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY);
497 
498  std::cout << "ERW: ME2 DD->DD (ref) = " << ME2ref << std::endl;
499  std::cout << "ERW: ME2 DD->DD = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl;
500 
501  // case UU -> UU
502  ID1 = 2; ID2 = 2; ID3 = 2; ID4 = 2;
503 
504  if(iter == 1) ME2ref = 1.9412924824120248E-017;
505  if(iter == 2) ME2ref = 4.0830132078559964E-017;
506  if(iter == 3) ME2ref = 2.1297931556738857E-018;
507  if(iter == 4) ME2ref = 7.9215386777281075E-019;
508 
509  ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY)
510  + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY)
511  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY)
512  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY);
513 
514  std::cout << "ERW: ME2 UU->UU (ref) = " << ME2ref << std::endl;
515  std::cout << "ERW: ME2 UU->UU = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl;
516 
517  // case DDX -> CCX
518  ID1 = 1; ID2 = -1; ID3 = 4; ID4 = -4;
519 
520  if(iter == 1) ME2ref = 3.6780119739685137E-020;
521  if(iter == 2) ME2ref = 5.2280923900274694E-018;
522  if(iter == 3) ME2ref = 2.9001589049209953E-019;
523  if(iter == 4) ME2ref = 5.8509026904432882E-020;
524 
525  ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY)
526  + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY)
527  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY)
528  + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY);
529 
530  std::cout << "ERW: ME2 DDX->CCX (ref) = " << ME2ref << std::endl;
531  std::cout << "ERW: ME2 DDX->CCX = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl;
532 
533 }
534 
535  //
536  // -------------------------------------------------------------------------------------------------------
537  // -------------------------------------------------------------------------------------------------------
538  // -------------------------------------------------------------------------------------------------------
539  //
540 
541 void makeSimpleTestPDF()
542 {
543 
544  std::cout << "------------------------------------" << std::endl;
545  std::cout << "ERW: TEST PDF only = " << std::endl;
546  std::cout << "ERW: PDF:LHAPDFset = cteq6ll.LHpdf " << std::endl;
547  std::cout << "ERW: reference from PYTHIA printout " << std::endl;
548  std::cout << "------------------------------------" << std::endl;
549  std::cout << "ERW: ID1,ID2 =2 -2 " << std::endl;
550  std::cout << "ERW: Q,Q2 = 97.868, 9578.118 " << std::endl;
551  double pdf1= f(2.51297e-01,2,9578.118,97.868);
552  std::cout << "ERW: x1, xpdf1= " << 2.51297e-01 << " " << 2.51297e-01 *pdf1 << " ref pdf = " << 0.411 << std::endl;
553  double pdf2= f(1.94463e-04,-2,9578.118,97.868);
554  std::cout << "ERW: x2, xpdf2= " << 1.94463e-04 << " " << 1.94463e-04 *pdf2 << " ref pdf = " << 2.989 << std::endl;
555  std::cout << "------------------------------------" << std::endl;
556  std::cout << "ERW: ID1,ID2 =2 -2 " << std::endl;
557  std::cout << "ERW: Q,Q2 = 38.747, 1501.314 " << std::endl;
558  pdf1= f(4.17758e-04,-2,1501.314,38.747);
559  std::cout << "ERW: x1, xpdf1 = " << 4.17758e-04 << " " << 4.17758e-04 *pdf1 << " ref pdf = " << 1.768 << std::endl;
560  pdf2= f(1.83354e-02, 2,1501.314,38.747);
561  std::cout << "ERW: x2, xpdf2 = " << 1.83354e-02 << " " << 1.83354e-02 *pdf2 << " ref pdf = " << 0.610 << std::endl;
562  std::cout << "------------------------------------" << std::endl;
563  std::cout << "ERW: ID1,ID2 =1 -1 " << std::endl;
564  std::cout << "ERW: Q,Q2 = 91.585, 8387.831 " << std::endl;
565  pdf1= f(6.29049e-02, 1,8387.831,91.585);
566  std::cout << "ERW: x1, xpdf1 = " << 6.29049e-02 << " " << 6.29049e-02 *pdf1 << " ref pdf = " << 0.386 << std::endl;
567  pdf2= f(6.80314e-04, -1,8387.831,91.585);
568  std::cout << "ERW: x2, xpdf2 = " << 6.80314e-04 << " " << 6.80314e-04 *pdf2 << " ref pdf = " << 1.751 << std::endl;
569  std::cout << "------------------------------------" << std::endl;
570  std::cout << "ERW: ID1,ID2 =-2 2 " << std::endl;
571  std::cout << "ERW: Q,Q2 = 12.979, 168.453 " << std::endl;
572  pdf1= f( 1.29739e-01, -2,168.453,12.979);
573  std::cout << "ERW: x1, xpdf1 = " << 1.29739e-01 << " " << 1.29739e-01*pdf1 << " ref pdf = " << 0.066 << std::endl;
574  pdf2= f( 6.62449e-06, 2,168.453,12.979);
575  std::cout << "ERW: x2, xpdf2 = " << 6.62449e-06 << " " << 6.62449e-06 *pdf2 << " ref pdf = " << 4.559 << std::endl;
576  std::cout << "------------------------------------" << std::endl;
577  std::cout << "ERW: ID1,ID2 = 2 -2 " << std::endl;
578  std::cout << "ERW: Q,Q2 = 91.707, 8410.127 " << std::endl;
579  pdf1= f(1.73131e-02 , 2,8410.127,91.707);
580  std::cout << "ERW: x1, xpdf1 = " << 1.73131e-02 << " " << 1.73131e-02 *pdf1 << " ref pdf = " << 0.643 << std::endl;
581  pdf2= f( 2.47840e-03, -2,8410.127,91.707);
582  std::cout << "ERW: x2, xpdf2 = " << 2.47840e-03 << " " << 2.47840e-03 *pdf2 << " ref pdf = " << 0.986 << std::endl;
583  std::cout << "------------------------------------" << std::endl;
584  std::cout << "ERW: ID1,ID2 = 1 -1 " << std::endl;
585  std::cout << "ERW: Q,Q2 = 91.203 8317.936 " << std::endl;
586  pdf1= f( 3.87999e-04 , 1, 8317.936, 91.203);
587  std::cout << "ERW: x1, xpdf1 = " << 3.87999e-04 << " " << 3.87999e-04 *pdf1 << " ref pdf = " << 2.256 << std::endl;
588  pdf2= f( 1.09378e-01, -1,8317.936, 91.203);
589  std::cout << "ERW: x2, xpdf2 = " << 1.09378e-01 << " " << 1.09378e-01 *pdf2 << " ref pdf = " << 0.106 << std::endl;
590  std::cout << "------------------------------------" << std::endl;
591  std::cout << "ERW: ID1,ID2 = 1 0 " << std::endl;
592  std::cout << "ERW: Q,Q2 = 51.963, 2700.201 " << std::endl;
593  pdf1= f( 9.22157e-02 , 1, 2700.201, 51.963);
594  std::cout << "ERW: x1, xpdf1 = " << 9.22157e-02 << " " << 9.22157e-02 *pdf1 << " ref pdf = " << 0.353 << std::endl;
595  pdf2= f( 2.50786e-03 , 0,2700.201, 51.963);
596  std::cout << "ERW: x2, xpdf2 = " << 2.50786e-03 << " " << 2.50786e-03 *pdf2 << " ref pdf = " << 20.362 << std::endl;
597  std::cout << "------------------------------------" << std::endl;
598  std::cout << "ERW: ID1,ID2 = 1 0 " << std::endl;
599  std::cout << "ERW: Q,Q2 = 55.719, 3104.648 " << std::endl;
600  pdf1= f( 1.28045e-01 , 1, 3104.648, 55.719);
601  std::cout << "ERW: x1, xpdf1 = " << 1.28045e-01 << " " << 1.28045e-01 *pdf1 << " ref pdf = " << 0.312 << std::endl;
602  pdf2= f( 3.12973e-03 , 0, 3104.648, 55.719);
603  std::cout << "ERW: x2, xpdf2 = " << 3.12973e-03 << " " << 3.12973e-03 *pdf2 << " ref pdf = " << 17.938 << std::endl;
604  std::cout << "------------------------------------" << std::endl;
605  std::cout << "ERW: ID1,ID2 = 3 0 " << std::endl;
606  std::cout << "ERW: Q,Q2 = 106.126, 11262.651 " << std::endl;
607  pdf1= f( 8.02832e-02 , 3, 11262.651, 106.126);
608  std::cout << "ERW: x1, xpdf1 = " << 8.02832e-02 << " " << 8.02832e-02 *pdf1 << " ref pdf = " << 0.077 << std::endl;
609  pdf2= f( 3.85011e-03 , 0, 11262.651, 106.126);
610  std::cout << "ERW: x2, xpdf2 = " << 3.85011e-03 << " " << 3.85011e-03 *pdf2 << " ref pdf = " << 16.542 << std::endl;
611  std::cout << "------------------------------------" << std::endl;
612 
613  std::cout << "ERW: Q= 98.069051344553060, x =1.36851248162658170E-002 " << std::endl;
614  double x = 1.36851248162658170E-002;
615  std::cout << "ERW: f(-5) = " << f( 1.36851248162658170E-002 , -5, 9617.53, 98.069 ) << " ref=12.649157008859417" << std::endl;
616  std::cout << "ERW: f(-4) = " << f( 1.36851248162658170E-002 , -4, 9617.53, 98.069 ) << " ref=19.196982264477899" << std::endl;
617  std::cout << "ERW: f(-3) = " << f( 1.36851248162658170E-002 , -3, 9617.53, 98.069 ) << " ref=23.971647011421041" << std::endl;
618  std::cout << "ERW: f(-2) = " << f( 1.36851248162658170E-002 , -2, 9617.53, 98.069 ) << " ref=30.602439987620933" << std::endl;
619  std::cout << "ERW: f(-1) = " << f( 1.36851248162658170E-002 , -1, 9617.53, 98.069 ) << " ref=31.664848276050574" << std::endl;
620  std::cout << "ERW: f(0) = " << f( 1.36851248162658170E-002 , 0, 9617.53, 98.069 ) << " ref=483.13004227606962" << std::endl;
621  std::cout << "ERW: f( 1) = " << f( 1.36851248162658170E-002 , 1, 9617.53, 98.069 ) << " ref=41.868651977841559" << std::endl;
622  std::cout << "ERW: f( 2) = " << f( 1.36851248162658170E-002 , 2, 9617.53, 98.069 ) << " ref=49.256024133250129" << std::endl;
623  std::cout << "ERW: f( 3) = " << f( 1.36851248162658170E-002 , 3, 9617.53, 98.069 ) << " ref=23.971647011421041" << std::endl;
624  std::cout << "ERW: f( 4) = " << f( 1.36851248162658170E-002 , 4, 9617.53, 98.069 ) << " ref=19.196982264477899" << std::endl;
625  std::cout << "ERW: f( 5) = " << f( 1.36851248162658170E-002 , 5, 9617.53, 98.069 ) << " ref=12.649157008859417" << std::endl;
626 
627  std::cout << "ERW: f(-5) = " << f( 1.36851248162658170E-002 , -5, 98.069, 98.069 ) << " ref=12.649157008859417" << std::endl;
628  std::cout << "ERW: f(-4) = " << f( 1.36851248162658170E-002 , -4, 98.069, 98.069 ) << " ref=19.196982264477899" << std::endl;
629  std::cout << "ERW: f(-3) = " << f( 1.36851248162658170E-002 , -3, 98.069, 98.069 ) << " ref=23.971647011421041" << std::endl;
630  std::cout << "ERW: f(-2) = " << f( 1.36851248162658170E-002 , -2, 98.069, 98.069 ) << " ref=30.602439987620933" << std::endl;
631  std::cout << "ERW: f(-1) = " << f( 1.36851248162658170E-002 , -1, 98.069, 98.069 ) << " ref=31.664848276050574" << std::endl;
632  std::cout << "ERW: f(0) = " << f( 1.36851248162658170E-002 , 0, 98.069, 98.069 ) << " ref=483.13004227606962" << std::endl;
633  std::cout << "ERW: f( 1) = " << f( 1.36851248162658170E-002 , 1, 98.069, 98.069 ) << " ref=41.868651977841559" << std::endl;
634  std::cout << "ERW: f( 2) = " << f( 1.36851248162658170E-002 , 2, 98.069, 98.069 ) << " ref=49.256024133250129" << std::endl;
635  std::cout << "ERW: f( 3) = " << f( 1.36851248162658170E-002 , 3, 98.069, 98.069 ) << " ref=23.971647011421041" << std::endl;
636  std::cout << "ERW: f( 4) = " << f( 1.36851248162658170E-002 , 4, 98.069, 98.069 ) << " ref=19.196982264477899" << std::endl;
637  std::cout << "ERW: f( 5) = " << f( 1.36851248162658170E-002 , 5, 98.069, 98.069 ) << " ref=12.649157008859417" << std::endl;
638 
639 
640 }
641 
642 
643 //---------------------------------------------------------------------------------------------------------------------------------------------------
644 void calcSumME2(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY,
645  int ID1, int ID2, int ID3, int ID4)
646 {
647  /*
648  // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation!
649  // FSR photons may need to be treated explicitely or with interpolation procedures.
650  Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(),
651  p3.py()+p4.py()+tau1.py()+tau2.py(),
652  p3.pz()+p4.pz()+tau1.pz()+tau2.pz(),
653  p3.e() +p4.e() +tau1.e() +tau2.e(), 0 );
654  */
655  Particle P_QQ( p3.px()+p4.px()+sp_X.px(),
656  p3.py()+p4.py()+sp_X.py(),
657  p3.pz()+p4.pz()+sp_X.pz(),
658  p3.e() +p4.e() +sp_X.e() , 0 );
659 
660 
661  double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
662 
663  double x1x2 = SS/CMSENE/CMSENE;
664  double x1Mx2 = P_QQ.pz()/CMSENE*2;
665 
666  double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
667  double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
668 
669  //---------------------------------------------------------------------------
670  // Construct the matrix for FORTRAN function
671  // NOTE: different order of indices than in FORTRAN!
672  double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 },
673  { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 },
674  { p3.e(), p3.px(), p3.py(), p3.pz() },
675  { p4.e(), p4.px(), p4.py(), p4.pz() },
676  { tau1.e(), tau1.px(), tau1.py(), tau1.pz() },
677  { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } };
678 
679 
680  W[0][0] = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P, KEY);
681  W[0][1] = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P, KEY);
682  W[1][0] = vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, -1, P, KEY);
683  W[1][1] = vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, 1, P, KEY);
684 
685  double sumW =(W[0][0]+W[0][1]+ W[1][0]+W[1][1]);
686  if( sumW < 0 ){
687  std::cout << "ERW: ==== " << std::endl;
688  std::cout << "ERW: W[] = " << W[0][0] << " " << W[0][1] << " " << W[1][0] << " " << W[1][1] << std::endl;
689  std::cout << "ERW: ==== " << std::endl;
690  std::cout << "ERW: p1 = " << CMSENE/2*x1 << " " << 0.0 << " " << 0.0 << " " << CMSENE/2*x1 << std::endl;
691  std::cout << "ERW: p2 = " << CMSENE/2*x2 << " " << 0.0 << " " << 0.0 << " " << -CMSENE/2*x2 << std::endl;
692  std::cout << "ERW: X = " << sp_X.e() << " " << sp_X.px() << " " << sp_X.py() << " " << sp_X.pz() << std::endl;
693  std::cout << "ERW: p3 = " << p3.e() << " " << p3.px() << " " << p3.py() << " " << p3.pz() << std::endl;
694  std::cout << "ERW: p4 = " << p4.e() << " " << p4.px() << " " << p4.py() << " " << p4.pz() << std::endl;
695  std::cout << "ERW: tau1 = " << tau1.e() << " " << tau1.px() << " " << tau1.py() << " " << tau1.pz() << std::endl;
696  std::cout << "ERW: tau2 = " << tau2.e() << " " << tau2.px() << " " << tau2.py() << " " << tau2.pz() << std::endl;
697  std::cout << "ERW: ==== " << std::endl;
698  std::cout << "ERW: p1+p2= " << P[0][0]+P[1][0] << " " << P[0][1]+P[1][1] << " " << P[0][2]+P[1][2] << " " << P[0][3]+P[1][3] << std::endl;
699  std::cout << "ERW: X+p3+p4= " << P[2][0]+P[3][0]+P[4][0]+P[5][0] << " " << P[2][1]+P[3][1]+P[4][1]+P[5][1] << " "
700  << P[2][2]+P[3][2]+P[4][2]+P[5][2] << " " << P[2][3]+P[3][3]+P[4][3]+P[5][3] << std::endl;
701 
702  double p3sqr = p3.e()*p3.e()-p3.px()*p3.px()-p3.py()*p3.py()-p3.pz()*p3.pz();
703  double p4sqr = p4.e()*p4.e()-p4.px()*p4.px()-p4.py()*p4.py()-p4.pz()*p4.pz();
704  std::cout << "ERW: p3^2, p4^2 =" << p3sqr << " " << p4sqr << std::endl;
705  }
706 
707 }
708 
709 //---------------------------------------------------------------------------------------------------------------------------------------------------
710 void calcPDFs(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY,
711  int ID1, int ID2, int ID3, int ID4, int pdfOpt)
712 {
713  /*
714  // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation!
715  // FSR photons may need to be treated explicitely or with interpolation procedures.
716  Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(),
717  p3.py()+p4.py()+tau1.py()+tau2.py(),
718  p3.pz()+p4.pz()+tau1.pz()+tau2.pz(),
719  p3.e() +p4.e() +tau1.e() +tau2.e(), 0 );
720  */
721  Particle P_QQ( p3.px()+p4.px()+sp_X.px(),
722  p3.py()+p4.py()+sp_X.py(),
723  p3.pz()+p4.pz()+sp_X.pz(),
724  p3.e() +p4.e() +sp_X.e() , 0 );
725 
726 
727  double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
728 
729  double x1x2 = SS/CMSENE/CMSENE;
730  double x1Mx2 = P_QQ.pz()/CMSENE*2;
731 
732  double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
733  double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
734 
735  //---------------------------------------------------------------------------
736  // Construct the matrix for FORTRAN function
737  // NOTE: different order of indices than in FORTRAN!
738  double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 },
739  { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 },
740  { p3.e(), p3.px(), p3.py(), p3.pz() },
741  { p4.e(), p4.px(), p4.py(), p4.pz() },
742  { tau1.e(), tau1.px(), tau1.py(), tau1.pz() },
743  { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } };
744 
745 
746  //
747  // Calculate 'f' function for all x1 and all ID1, ID2
748  //
749  double f_x1_ARRAY[11] = { 0.0 };
750  double f_x2_ARRAY[11] = { 0.0 };
751  double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with
752  double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5
753 
754 
755  Particle P_X( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e() , 0 );
756 
757 
758  double Q2 = P_X.recalculated_mass()*P_X.recalculated_mass();
759  double PT2 = sp_X.px()* sp_X.px() + sp_X.py()* sp_X.py();
760 
761  if(pdfOpt == 0){
762  for(int i=-5;i<=5;++i) {
763  f_x1[i] = f(x1,i,Q2,CMSENE);
764  f_x2[i] = f(x2,i,Q2,CMSENE);
765  }
766  } else if (pdfOpt == 1) {
767  for(int i=-5;i<=5;++i) {
768  f_x1[i] = f(x1,i,PT2,CMSENE);
769  f_x2[i] = f(x2,i,PT2,CMSENE);
770  }
771  } else if (pdfOpt == 2) {
772  double PT24 = PT2/4.;
773  for(int i=-5;i<=5;++i) {
774  f_x1[i] = f(x1,i,PT24,CMSENE);
775  f_x2[i] = f(x2,i,PT24,CMSENE);
776  }
777  }
778 
779 
780  int I1=ID1;
781  int I2=ID2;
782  if( I1 == 21) I1=0;
783  if( I2 == 21) I2=0;
784 
785  W[0][0] = f_x1[I1]*f_x2[I2];
786  W[0][1] = f_x1[I1]*f_x2[I2];
787  W[1][0] = f_x1[I1]*f_x2[I2];
788  W[1][1] = f_x1[I1]*f_x2[I2];
789 
790 
791 }
792 
793 
794 //---------------------------------------------------------------------------------------------------------------------------------------------------
795 void calcProdMatrix(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY,
796  int ID1, int ID2, int ID3, int ID4, int pdfOpt)
797 {
798  /*
799  // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation!
800  // FSR photons may need to be treated explicitely or with interpolation procedures.
801  Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(),
802  p3.py()+p4.py()+tau1.py()+tau2.py(),
803  p3.pz()+p4.pz()+tau1.pz()+tau2.pz(),
804  p3.e() +p4.e() +tau1.e() +tau2.e(), 0 );
805  */
806  Particle P_QQ( p3.px()+p4.px()+sp_X.px(),
807  p3.py()+p4.py()+sp_X.py(),
808  p3.pz()+p4.pz()+sp_X.pz(),
809  p3.e() +p4.e() +sp_X.e() , 0 );
810 
811 
812  double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
813 
814  double x1x2 = SS/CMSENE/CMSENE;
815  double x1Mx2 = P_QQ.pz()/CMSENE*2;
816 
817  double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
818  double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
819 
820  //---------------------------------------------------------------------------
821  // Construct the matrix for FORTRAN function
822  // NOTE: different order of indices than in FORTRAN!
823  double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 },
824  { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 },
825  { p3.e(), p3.px(), p3.py(), p3.pz() },
826  { p4.e(), p4.px(), p4.py(), p4.pz() },
827  { tau1.e(), tau1.px(), tau1.py(), tau1.pz() },
828  { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } };
829 
830 
831  //
832  // Calculate 'f' function for all x1 and all ID1, ID2
833  //
834  double f_x1_ARRAY[11] = { 0.0 };
835  double f_x2_ARRAY[11] = { 0.0 };
836  double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with
837  double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5
838 
839 
840  Particle P_X( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e() , 0 );
841 
842 
843  double Q2 = P_X.recalculated_mass()*P_X.recalculated_mass();
844  double PT2 = sp_X.px()* sp_X.px() + sp_X.py()* sp_X.py();
845 
846  if(pdfOpt == 0){
847  for(int i=-5;i<=5;++i) {
848  f_x1[i] = f(x1,i,Q2,CMSENE);
849  f_x2[i] = f(x2,i,Q2,CMSENE);
850  }
851  } else if (pdfOpt == 1) {
852  for(int i=-5;i<=5;++i) {
853  f_x1[i] = f(x1,i,PT2,CMSENE);
854  f_x2[i] = f(x2,i,PT2,CMSENE);
855  }
856  } else if (pdfOpt == 2) {
857  double PT24 = PT2/4.;
858  for(int i=-5;i<=5;++i) {
859  f_x1[i] = f(x1,i,PT24,CMSENE);
860  f_x2[i] = f(x2,i,PT24,CMSENE);
861  }
862  }
863 
864 
865  int I1=ID1;
866  int I2=ID2;
867  if( I1 == 21) I1=0;
868  if( I2 == 21) I2=0;
869 
870  W[0][0] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P, KEY);
871  W[0][1] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P, KEY);
872  W[1][0] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, -1, P, KEY);
873  W[1][1] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, 1, P, KEY);
874 
875  double sumW =(W[0][0]+W[0][1]+ W[1][0]+W[1][1]);
876  if( sumW < 0 ){
877  std::cout << "ERW: f1*f2= " << f_x1[I1]*f_x2[I2] << " x1=" << x1 << " x2=" << x2 << " f1=" << f_x1[I1] << " f2=" << f_x2[I1] << std::endl;
878  std::cout << "ERW: ==== " << std::endl;
879  std::cout << "ERW: W[] = " << W[0][0] << " " << W[0][1] << " " << W[1][0] << " " << W[1][1] << std::endl;
880  std::cout << "ERW: ==== " << std::endl;
881  std::cout << "ERW: p1 = " << CMSENE/2*x1 << " " << 0.0 << " " << 0.0 << " " << CMSENE/2*x1 << std::endl;
882  std::cout << "ERW: p2 = " << CMSENE/2*x2 << " " << 0.0 << " " << 0.0 << " " << -CMSENE/2*x2 << std::endl;
883  std::cout << "ERW: X = " << sp_X.e() << " " << sp_X.px() << " " << sp_X.py() << " " << sp_X.pz() << std::endl;
884  std::cout << "ERW: p3 = " << p3.e() << " " << p3.px() << " " << p3.py() << " " << p3.pz() << std::endl;
885  std::cout << "ERW: p4 = " << p4.e() << " " << p4.px() << " " << p4.py() << " " << p4.pz() << std::endl;
886  std::cout << "ERW: tau1 = " << tau1.e() << " " << tau1.px() << " " << tau1.py() << " " << tau1.pz() << std::endl;
887  std::cout << "ERW: tau2 = " << tau2.e() << " " << tau2.px() << " " << tau2.py() << " " << tau2.pz() << std::endl;
888  std::cout << "ERW: ==== " << std::endl;
889  std::cout << "ERW: p1+p2= " << P[0][0]+P[1][0] << " " << P[0][1]+P[1][1] << " " << P[0][2]+P[1][2] << " " << P[0][3]+P[1][3] << std::endl;
890  std::cout << "ERW: X+p3+p4= " << P[2][0]+P[3][0]+P[4][0]+P[5][0] << " " << P[2][1]+P[3][1]+P[4][1]+P[5][1] << " "
891  << P[2][2]+P[3][2]+P[4][2]+P[5][2] << " " << P[2][3]+P[3][3]+P[4][3]+P[5][3] << std::endl;
892 
893  double p3sqr = p3.e()*p3.e()-p3.px()*p3.px()-p3.py()*p3.py()-p3.pz()*p3.pz();
894  double p4sqr = p4.e()*p4.e()-p4.px()*p4.px()-p4.py()*p4.py()-p4.pz()*p4.pz();
895  std::cout << "ERW: p3^2, p4^2 =" << p3sqr << " " << p4sqr << std::endl;
896  }
897 
898 }
899 
900 
901 /*******************************************************************************
902  Calculate weights, case of event record vertex like Z/gamma/H ... -> tau tau decay plus 2 jets.
903 
904  Determine taus decay channel, calculates all necessary components from decay HHp HHm and production W[][] for
905  calculation of all weights, later calculates weights.
906 
907  Input: X four momentum of Z/Higgs may be larger than sum of tau1 tau2, missing component
908  is assumed to be QED brem four momenta of outgoing partons and taus; vectors of decay products for first and second tau
909 
910  Hidden input: relWTnonSM, switch of what is WTnonSM
911  Hidden output: WTamplitM or WTamplitP matrix elements of tau1 tau2 decays
912 
913 
914  Polari - helicity attributed to taus, 100% correlations between tau+ and tau-
915  accesible with getTauSpin()
916  WTnonSM weight for introduction of matrix element due to nonSM in cross section, or just ME^2*PDF's for relWTnonS==false
917  Explicit output: WT spin correlation weight
918 *******************************************************************************/
919 double calculateWeightFromParticlesVBFPROD(int IDPROD, SimpleParticle &p3, SimpleParticle &p4,SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector<SimpleParticle> &sp_tau1_daughters, vector<SimpleParticle> &sp_tau2_daughters, int KEY)
920 {
921  SimpleParticle sp_tau;
922  SimpleParticle sp_nu_tau;
923  vector<SimpleParticle> sp_tau_daughters;
924 
925  // First iteration is for tau plus, so the 'nu_tau' is tau minus
926  if (sp_tau1.pdgid() == -15 )
927  {
928  sp_tau = sp_tau1;
929  sp_nu_tau = sp_tau2;
930  sp_tau_daughters = sp_tau1_daughters;
931  }
932  else
933  {
934  sp_tau = sp_tau2;
935  sp_nu_tau = sp_tau1;
936  sp_tau_daughters = sp_tau2_daughters;
937  }
938 
939  double *HHp, *HHm;
940 
941  // We use this to separate namespace for tau+ and tau-
942  if(true)
943  {
944  // Create Particles from SimpleParticles
945  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
946  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
947  Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
948 
949  vector<Particle> tau_daughters;
950 
951  // tau pdgid
952  int tau_pdgid = sp_tau.pdgid();
953 
954  // Create list of tau daughters
955  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
956  {
957  Particle pp(sp_tau_daughters[i].px(),
958  sp_tau_daughters[i].py(),
959  sp_tau_daughters[i].pz(),
960  sp_tau_daughters[i].e(),
961  sp_tau_daughters[i].pdgid() );
962 
963  tau_daughters.push_back(pp);
964  }
965 
966  double phi2 = 0.0, theta2 = 0.0;
967 
968 
969  // Move decay kinematics first to tau rest frame with z axis pointing along nu_tau direction
970  // later rotate again to have neutrino from tau decay along z axis: angles phi2, theta2
971  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
972 
973 
974  // Determine decay channel and then calculate polarimetric vector HH
975  HHp = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
976 
977  DEBUG
978  (
979  cout<<tau_pdgid<<" -> ";
980  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
981  cout<<" (HHp: "<<HHp[0]<<" "<<HHp[1]<<" "<<HHp[2]<<" "<<HHp[3]<<") ";
982  cout<<endl;
983  )
984 
985  WTamplitP = WTamplit;
986  } // end of tau+
987 
988  // Second iteration is for tau minus, so the 'nu_tau' is tau minus
989  if(sp_tau1.pdgid() == 15 )
990  {
991  sp_tau = sp_tau1;
992  sp_nu_tau = sp_tau2;
993  sp_tau_daughters = sp_tau1_daughters;
994  }
995  else
996  {
997  sp_tau = sp_tau2;
998  sp_nu_tau = sp_tau1;
999  sp_tau_daughters = sp_tau2_daughters;
1000  }
1001 
1002  // We use this to separate namespace for tau+ and tau-
1003  if(true)
1004  {
1005  // Create Particles from SimpleParticles
1006  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
1007  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
1008  Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
1009 
1010  vector<Particle> tau_daughters;
1011 
1012  // tau pdgid
1013  int tau_pdgid = sp_tau.pdgid();
1014 
1015  // Create list of tau daughters
1016  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
1017  {
1018  Particle pp(sp_tau_daughters[i].px(),
1019  sp_tau_daughters[i].py(),
1020  sp_tau_daughters[i].pz(),
1021  sp_tau_daughters[i].e(),
1022  sp_tau_daughters[i].pdgid() );
1023 
1024  tau_daughters.push_back(pp);
1025  }
1026 
1027  double phi2 = 0.0, theta2 = 0.0;
1028 
1029 
1030  // Move decay kinematics first to tau rest frame with z axis pointing along nu_tau direction
1031  // later rotate again to have neutrino from tau decay along z axis: angles phi2, theta2
1032  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
1033 
1034 
1035  // Determine decay channel and then calculate polarimetric vector HH
1036  HHm = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
1037 
1038  DEBUG
1039  (
1040  cout<<tau_pdgid<<" -> ";
1041  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
1042  cout<<" (HHm: "<<HHm[0]<<" "<<HHm[1]<<" "<<HHm[2]<<" "<<HHm[3]<<") ";
1043  cout<<endl;
1044  )
1045 
1046  WTamplitM = WTamplit;
1047  } // end of tau-
1048 
1049 
1050  double W[2][2] = { { 0.25, 0.25 },
1051  { 0.25, 0.25 } }; // this is trivial W spin (helicity only) density matrix
1052  // consist of ME^2*PDF's for production.
1053  // pre-initialization all equal i.e. no spin effects
1054 
1055 
1056  // calculate ME^2*PDF's for SM
1057  if(KEY==0 || KEY==1 )calcXsect(IDPROD, p3, p4, sp_X, sp_tau1, sp_tau2, W, KEY);
1058  else calcXsect(IDPROD, p3, p4, sp_X, sp_tau1, sp_tau2, W, 0);
1059 
1060 
1061  double sum=(W[0][0]+W[0][1]+ W[1][0]+W[1][1]); // getVBFspins calculated PDF's times matrix elements squared for each of four helicity configurations
1062  // tau+ and tau- using SM process KEY=0 DY, KEY=1 Higgs.
1063 
1064  if(KEY==0 || KEY==1 ) {
1065  WTnonSM=1.0;
1066  if(relWTnonSM==0) WTnonSM=sum; // The variable accessible for user stores production weight ME^2 * PDF's summed over flavours and spin states
1067  }
1068  if(KEY>=2) {
1069  // now we re-calculate W using anomalous variant of matrix elements.
1070  calcXsect(IDPROD, p3, p4, sp_X, sp_tau1, sp_tau2, W, KEY);
1071 
1072  double sum2=(W[0][0]+W[0][1]+ W[1][0]+W[1][1]);
1073 
1074  WTnonSM=sum2/sum;
1075  if(relWTnonSM==0) WTnonSM=sum2;
1076  }
1077 
1078  double WT = W[0][0]*(1+HHp[2])*(1+HHm[2])+W[0][1]*(1+HHp[2])*(1-HHm[2])+ W[1][0]*(1-HHp[2])*(1+HHm[2])+W[1][1]*(1-HHp[2])*(1-HHm[2]);
1079  WT = WT/(W[0][0]+W[0][1]+ W[1][0]+W[1][1]);
1080 
1081  // we separate cross section into first tau helicity + and - parts. Second tau follow.
1082  // This must be tested especially for KEY >=2
1083  // case for scalar (H) has flipped sign of second tau spin? Tests needed
1084  double RRR = Tauola::randomDouble();
1085  Polari=1.0;
1086  if (RRR<(W[0][0]*(1+HHp[2])*(1+HHm[2])+W[0][1]*(1+HHp[2])*(1-HHm[2]))/(W[0][0]*(1+HHp[2])*(1+HHm[2])+W[0][1]*(1+HHp[2])*(1-HHm[2])+W[1][0]*(1-HHp[2])*(1+HHm[2])+W[1][1]*(1-HHp[2])*(1-HHm[2]))) Polari=-1.0;
1087 
1088 
1089 
1090  // Print out some info about the channel
1091  DEBUG( cout<<" WT: "<<WT<<endl; )
1092 
1093  if (WT<0.0) {
1094  printf("WT is: %13.10f. Setting WT = 0.0\n",WT);
1095  WT = 0.0;
1096  }
1097 
1098  if (WT>4.0) {
1099  printf("WT is: %13.10f. Setting WT = 4.0\n",WT);
1100  WT = 4.0;
1101  }
1102 
1103  delete[] HHp;
1104  delete[] HHm;
1105 
1106  return WT;
1107 }
1108 
1109 } // namespace TauSpinner
1110 
TauSpinner::vbfdistr_
double vbfdistr_(int *I1, int *I2, int *I3, int *I4, int *H1, int *H2, double P[6][4], int *KEY)
TauSpinner::makeSimpleTestME2
void makeSimpleTestME2()
Definition: vbftests.cxx:356
TauSpinner::vbfdistrWRAP
double vbfdistrWRAP(int I1, int I2, int I3, int I4, int H1, int H2, double P[6][4], int KEY)
Definition: vbftests.cxx:31
TauSpinner::calcSumME2
void calcSumME2(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X, SimpleParticle &tau1, SimpleParticle &tau2, double(&W)[2][2], int KEY, int ID1, int ID2, int ID3, int ID4)
Definition: vbftests.cxx:644
TauSpinner::calcPDFs
void calcPDFs(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X, SimpleParticle &tau1, SimpleParticle &tau2, double(&W)[2][2], int KEY, int ID1, int ID2, int ID3, int ID4, int pdfOpt)
Definition: vbftests.cxx:710
TauSpinner::makeSimpleTestPDF
void makeSimpleTestPDF()
Definition: vbftests.cxx:541
TauSpinner::calcXsect
void calcXsect(int IDPROD, SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X, SimpleParticle &tau1, SimpleParticle &tau2, double(&W)[2][2], int KEY)
Definition: vbftests.cxx:41
TauSpinner::calculateHH
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
Definition: tau_reweight_lib.cxx:1082
TauSpinner::prepareKinematicForHH
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
Definition: tau_reweight_lib.cxx:985
TauSpinner::calcProdMatrix
void calcProdMatrix(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X, SimpleParticle &tau1, SimpleParticle &tau2, double(&W)[2][2], int KEY, int ID1, int ID2, int ID3, int ID4, int pdfOpt)
Definition: vbftests.cxx:795
Generated by doxygen 1.9.1