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>
6using namespace Tauolapp;
7
8/*This is a group of methods prepared for the testing purpose. The tests include those where fixed kinematics is
9forced. Then external calculation can be used for matrix element calculation and/or PDFs etc.
10Gradually better documentation will be introduced. It will be documented with comment lines and in the draft of the paper.
11To use one has to modify the Makefile the appropriate line of
12objects has to be activated (# removed). Then invoking lines in code have to be commented out as well.
13At 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
16namespace TauSpinner {
17
18extern double CMSENE;
19extern int relWTnonSM;
20extern double WTnonSM;
21extern double Polari;
22extern double WTamplit;
23extern double WTamplitP;
24extern double WTamplitM;
25extern double CMSENE;
26extern double f(double x, int ID, double SS, double cmsene);
27
28const bool DEBUG = 0;
29
30/** Wrapper to VBDISTR */
31double 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
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
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//---------------------------------------------------------------------------------------------------------------------------------------------------
644void 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//---------------------------------------------------------------------------------------------------------------------------------------------------
710void 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//---------------------------------------------------------------------------------------------------------------------------------------------------
795void 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*******************************************************************************/
919double 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
double vbfdistr_(int *I1, int *I2, int *I3, int *I4, int *H1, int *H2, double P[6][4], int *KEY)
void makeSimpleTestME2()
Definition: vbftests.cxx:356
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
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
void makeSimpleTestPDF()
Definition: vbftests.cxx:541
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
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
double vbfdistrWRAP(int I1, int I2, int I3, int I4, int H1, int H2, double P[6][4], int KEY)
Definition: vbftests.cxx:31
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