C++ Interface to Tauola
vbfdistr.cxx
1#include "TauSpinner/tau_reweight_lib.h"
2#include "TauSpinner/Particle.h"
3#include "TauSpinner/vbfdistr.h"
4#include <cstdlib>
5#include <string.h>
6using namespace Tauolapp;
7
8namespace TauSpinner {
9
10extern double CMSENE;
11extern int nonSM2;
12extern int relWTnonSM;
13extern double WTnonSM;
14extern double Polari;
15extern double WTamplit;
16extern double WTamplitP;
17extern double WTamplitM;
18extern double CMSENE;
19extern double f(double x, int ID, double SS, double cmsene);
20
21const bool DEBUG = 0;
22 int _QCDdefault=1;
23 int _QCDvariant=1;
24
25/** @brief Pointer to a function that modify rezult of alpha_ calc. in vbfdist
26 *
27 * By default this function is not used (pointer is NULL).
28 * It can be changed by the user through TauSpinner::set_vbfdistrModif
29 */
30void (*alphasModif)(double, int, int) = NULL;
31
32/** @brief Set vbfdistrModif function */
33 void set_alphasModif(void (*function)(double, int, int) )
34{
35 alphasModif = function;
36}
37
38/** @brief Pointer to a function that modify rezult of Matrix Element calc. in vbfdist
39 *
40 * By default this function is not used (pointer is NULL).
41 * It can be changed by the user through TauSpinner::set_vbfdistrModif
42 */
43double (*vbfdistrModif)(int, int, int, int, int, int, double[6][4], int, double) = NULL;
44
45/** @brief Set vbfdistrModif function */
46void set_vbfdistrModif(double (*function)(int, int, int, int, int, int, double[6][4], int, double) )
47{
48 vbfdistrModif = function;
49}
50
51void setPDFOpt(int QCDdefault,int QCDvariant){
52 _QCDdefault=QCDdefault;
53 _QCDvariant=QCDvariant;
54}
55
56 int getPDFOpt(int KEY){
57 if (KEY==0 || KEY==1) return _QCDdefault;
58 else return _QCDvariant;
59}
60
61 void alphas(double Q2,int scalePDFOpt, int KEY){
62 if (alphasModif)
63 { alphasModif(Q2,scalePDFOpt,KEY);
64 }
65 else
66 {
67 // any textbook LL calculation
68 const double PI=3.14159265358979324;
69 // number of flavors
70 const double Nf = 5;
71 // alphas_s at mZ
72 const double AlphasMZ = 0.118;
73 const double MZ =91.1876;
74 // alphas_s at scale Q2 (GeV^2)
75 double alfas=AlphasMZ / ( 1 + AlphasMZ/(4*PI) * (11 - 2./3 * Nf) * log(Q2/(MZ*MZ)));
76 // test Alphas ( Q2 = 1000^2 ) = 0.0877445
77 if(scalePDFOpt==0) alfas = 0.118;
78 if(params_r_.as != alfas){// we pass alpha_s to calculation of amplitudes
79 params_r_.as = alfas;
80 // reinitialize constants couplings for amplitude calculations
81 vbf_reinit_(&KEY);
82 }
83 }
84}
85
86/** Wrapper to VBDISTR and place for interface to user provided modification*/
87double vbfdistr(int I1, int I2, int I3, int I4, int H1, int H2, double P[6][4], int KEY)
88{
89 double P_copy[6][4];
90 double original_result = 0.;
91
92 memcpy(P_copy,P,sizeof(double)*6*4);
93
94 if( KEY<2) {
95 // SM mode
96 return vbfdistr_(&I1, &I2, &I3, &I4, &H1, &H2, P_copy, &KEY);
97 }
98 else if( !vbfdistrModif) {
99 printf("TauSpinner::vbfdistr: User function vbfdistrModif not declared. Setting WT_contrib = 0.0 Failed attempt with KEY = %i5. \n",KEY);
100 return original_result;
101 }
102 else if( KEY<4) {
103 // modification mode
104 int KEY_BUF=KEY-2;
105 original_result= vbfdistr_(&I1, &I2, &I3, &I4, &H1, &H2, P_copy, &KEY_BUF);
106 return vbfdistrModif(I1,I2,I3,I4,H1,H2,P_copy,KEY,original_result);
107 }
108 else {
109 // replacement mode
110 return vbfdistrModif(I1,I2,I3,I4,H1,H2,P_copy,KEY,original_result);
111 }
112}
113
114
115/** Get VBF ME2
116 Returns array W[2][2] */
117//---------------------------------------------------------------------------------------------------------------------------------------------------
118void getME2VBF(SimpleParticle &p3i, SimpleParticle &p4i, SimpleParticle &sp_X,SimpleParticle &tau1i, SimpleParticle &tau2i, double (&W)[2][2], int KEY)
119{
120
121 // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation!
122 // FSR photons may need to be treated explicitely or with interpolation procedures.
123 Particle p3(p3i.px(),p3i.py(),p3i.pz(),p3i.e(),0);
124 Particle p4(p4i.px(),p4i.py(),p4i.pz(),p4i.e(),0);
125 double mtaul=sqrt(tau1i.e()*tau1i.e()-tau1i.px()*tau1i.px()-tau1i.py()*tau1i.py()-tau1i.pz()*tau1i.pz());
126 Particle tau1(tau1i.px(),tau1i.py(),tau1i.pz(),tau1i.e(),0);
127 Particle tau2(tau2i.px(),tau2i.py(),tau2i.pz(),tau2i.e(),0);
128
129 // TEST begin
130 /*
131 if(KEY>1){ // this is special arrangement to used alternative calculation to introduce effect of smearing on tau lepton pair (hidden in sp_X) due to showering
132 // in this way we get weights for corelated configs; with and without the shower-kick.
133 Particle P_tautau(tau1.px()+tau2.px(),tau1.py()+tau2.py(),tau1.pz()+tau2.pz(),tau1.e()+tau2.e(),0);
134 tau1.boostToRestFrame(P_tautau);
135 tau2.boostToRestFrame(P_tautau);
136 p3.boostToRestFrame(P_tautau);
137 p4.boostToRestFrame(P_tautau);
138 //v1.boostToRestFrame(P_tautau);
139 //v2.boostToRestFrame(P_tautau);
140
141 Particle P_X(sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(),0);
142 tau1.boostFromRestFrame(P_X);
143 tau2.boostFromRestFrame(P_X);
144 p3.boostFromRestFrame(P_X);
145 p4.boostFromRestFrame(P_X);
146 //v1.boostFromRestFrame(P_X);
147 //v2.boostFromRestFrame(P_X);
148 }
149 */
150 // TEST end
151
152 // we may want to force p3,p4 to be masless.
153 Particle v1(0.,0., 1.,1.,0.);
154 Particle v2(0.,0.,-1.,1.,0.);
155
156 Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px()+1e-8*p3.pz(),
157 p3.py()+p4.py()+tau1.py()+tau2.py()-2e-8*p3.pz(),
158 p3.pz()+p4.pz()+tau1.pz()+tau2.pz(),
159 p3.e() +p4.e() +tau1.e() +tau2.e(), 0 );
160 tau1.boostToRestFrame(P_QQ);
161 tau2.boostToRestFrame(P_QQ);
162 p3.boostToRestFrame(P_QQ);
163 p4.boostToRestFrame(P_QQ);
164 v1.boostToRestFrame(P_QQ);
165 v2.boostToRestFrame(P_QQ);
166
167 // Now we can define vers1 vers2
168 double xn1=sqrt(v1.px()*v1.px()+v1.py()*v1.py()+v1.pz()*v1.pz());
169 double xn2=sqrt(v2.px()*v2.px()+v2.py()*v2.py()+v2.pz()*v2.pz());
170 double xn12=sqrt( (v1.px()-v2.px())*(v1.px()-v2.px()) +(v1.py()-v2.py())*(v1.py()-v2.py())+(v1.pz()-v2.pz())*(v1.pz()-v2.pz()));
171
172 double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
173
174 double x1x2 = SS/CMSENE/CMSENE;
175 double x1Mx2 = P_QQ.pz()/CMSENE*2;
176
177 double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
178 double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
179
180 //---------------------------------------------------------------------------
181 // Construct the matrix for FORTRAN function
182 // NOTE: different order of indices than in FORTRAN!
183 // four options for partonic beams.
184 double P[6][4] = { { sqrt(SS)/2, sqrt(SS)/2/xn12*(v1.px()-v2.px()), sqrt(SS)/2/xn12*(v1.py()-v2.py()), sqrt(SS)/2/xn12*(v1.pz()-v2.pz()) },
185 { sqrt(SS)/2, -sqrt(SS)/2/xn12*(v1.px()-v2.px()), -sqrt(SS)/2/xn12*(v1.py()-v2.py()), -sqrt(SS)/2/xn12*(v1.pz()-v2.pz()) },
186 // double P[6][4] = { { sqrt(SS)/2, -sqrt(SS)/2/xn1*v2.px(), -sqrt(SS)/2/xn1*v2.py(), -sqrt(SS)/2/xn1*v2.pz() },
187 // { sqrt(SS)/2, sqrt(SS)/2/xn1*v2.px(), sqrt(SS)/2/xn1*v2.py(), sqrt(SS)/2/xn1*v2.pz() },
188 // double P[6][4] = { { sqrt(SS)/2, sqrt(SS)/2/xn1*v1.px(), sqrt(SS)/2/xn1*v1.py(), sqrt(SS)/2/xn1*v1.pz() },
189 // { sqrt(SS)/2, -sqrt(SS)/2/xn1*v1.px(), -sqrt(SS)/2/xn1*v1.py(), -sqrt(SS)/2/xn1*v1.pz() },
190 // double P[6][4] = { { sqrt(SS)/2, 0.0, 0.0, sqrt(SS)/2 },
191 // { sqrt(SS)/2, 0.0, 0.0, -sqrt(SS)/2 },
192 { p3.e(), p3.px(), p3.py(), p3.pz() },
193 { p4.e(), p4.px(), p4.py(), p4.pz() },
194 { tau1.e(), tau1.px(), tau1.py(), tau1.pz() },
195 { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } };
196
197 //
198 // Calculate 'f' function for all x1 and all ID1, ID2
199 //
200 double f_x1_ARRAY[11] = { 0.0 };
201 double f_x2_ARRAY[11] = { 0.0 };
202 double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with
203 double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5
204
205 Particle P_tautau( tau1.px()+tau2.px(), tau1.py()+tau2.py(), tau1.pz()+tau2.pz(), tau1.e()+tau2.e(),0 );
206
207 double Q2 = P_tautau.recalculated_mass()*P_tautau.recalculated_mass();
208 double QQ2 = P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
209 double PT2 = P_QQ.px() * P_QQ.px() + P_QQ.py()* P_QQ.py();
210
211 double sumET = p3.e()*sqrt( p3.px()*p3.px()+p3.py()*p3.py())/sqrt( p3.px()*p3.px()+p3.py()*p3.py()+p3.pz()*p3.pz())
212 + p4.e()*sqrt( p4.px()*p4.px()+p4.py()*p4.py())/sqrt( p4.px()*p4.px()+p4.py()*p4.py()+p4.pz()*p4.pz())
213 + tau1.e()*sqrt( tau1.px()*tau1.px()+tau1.py()*tau1.py())/sqrt( tau1.px()*tau1.px()+tau1.py()*tau1.py()+tau1.pz()*tau1.pz() )
214 + tau2.e()*sqrt( tau2.px()*tau2.px()+tau2.py()*tau2.py())/sqrt( tau2.px()*tau2.px()+tau2.py()*tau2.py()+tau2.pz()*tau2.pz() ) ;
215
216 double sumMT = sqrt( p3.recalculated_mass()*p3.recalculated_mass() + p3.px()*p3.px()+p3.py()*p3.py() )
217 + sqrt( p4.recalculated_mass()*p4.recalculated_mass() + p4.px()*p4.px()+p4.py()*p4.py() )
218 + sqrt( tau1.recalculated_mass()*tau1.recalculated_mass() + tau1.px()*tau1.px()+tau1.py()*tau1.py() )
219 + sqrt( tau2.recalculated_mass()*tau2.recalculated_mass() + tau2.px()*tau2.px()+tau2.py()*tau2.py() );
220
221 double fixed_scale = 91.17; // 200.;
222 int scalePDFOpt = 1; // it has to be connected to initialization and flipper
223 scalePDFOpt=getPDFOpt(KEY);
224 double Q2pdf = fixed_scale*fixed_scale;
225 if(scalePDFOpt == 0){
226 //
227 alphas(Q2pdf,scalePDFOpt,KEY);
228 } else if (scalePDFOpt == 1) {
229 Q2pdf = QQ2;
230 alphas(Q2pdf,scalePDFOpt,KEY);
231 } else if (scalePDFOpt == 2) {
232 Q2pdf = sumMT*sumMT;
233 alphas(Q2pdf,scalePDFOpt,KEY);
234 } else if (scalePDFOpt == 3) {
235 Q2pdf = sumET*sumET;
236 alphas(Q2pdf,scalePDFOpt,KEY);
237 } else if (scalePDFOpt == 4) {
238 Q2pdf = Q2;
239 alphas(Q2pdf,scalePDFOpt,KEY);
240 }
241
242 for(int i=-5;i<=5;++i) {
243 f_x1[i] = f(x1,i,Q2pdf,CMSENE);
244 f_x2[i] = f(x2,i,Q2pdf,CMSENE);
245 }
246 // reset to zero
247 W[0][0]=0.0;
248 W[0][1]=0.0;
249 W[1][0]=0.0;
250 W[1][1]=0.0;
251
252 if( DEBUG ){
253 // here you can overwrite kinematical configurations for tests
254 std::cout << " " << std::endl;
255 std::cout << "---podstawiamy 4-pedy na jakies ---" << std::endl;
256 std::cout << "ERW: ----------------------------- " << std::endl;
257 double P1[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
258 { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
259 { 0.8855009E+02, -0.2210038E+02, 0.4007979E+02, -0.7580437E+02 },
260 { 0.3283248E+03, -0.1038482E+03, -0.3019295E+03, 0.7649385E+02 },
261 { 0.1523663E+03, -0.1058795E+03, -0.9770827E+02, 0.4954769E+02 },
262 { 0.4307588E+03, 0.2318280E+03, 0.3595580E+03, -0.5023717E+02 } };
263 std::cout << "ERW: ----------------------------- " << std::endl;
264 double P2[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
265 { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
266 { 0.1177462E+03, -0.6070512E+02, 0.7123011E+02, 0.7145150E+02 },
267 { 0.3509495E+03, -0.3178775E+02, 0.8393832E+02, 0.3392779E+03 },
268 { 0.3493321E+03, 0.1840069E+03, -0.5152712E+02, -0.2924315E+03 },
269 { 0.1819722E+03, -0.9151401E+02, -0.1036413E+03, -0.1182978E+03 } };
270
271 std::cout << "ERW: ----------------------------- " << std::endl;
272 double P3[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
273 { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
274 { 0.2586900E+03, 0.1324670E+03, -0.1696171E+03, -0.1435378E+03 },
275 { 0.1084567E+03, -0.5735712E+02, -0.2162482E+02, -0.8947281E+02 },
276 { 0.4005742E+03, -0.1580760E+03, 0.3563160E+03, 0.9223569E+02 },
277 { 0.2322791E+03, 0.8296613E+02, -0.1650741E+03, 0.1407749E+03 } };
278 std::cout << "ERW: ----------------------------- " << std::endl;
279 double P4[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 },
280 { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 },
281 { 0.1595700E+03, -0.6917808E+02, -0.1395175E+03, -0.3481123E+02 },
282 { 0.2247758E+03, -0.1360140E+03, 0.1650340E+03, -0.6919641E+02 },
283 { 0.2508802E+03, 0.1447863E+01, 0.2499830E+03, -0.2107335E+02 },
284 { 0.3647740E+03, 0.2037442E+03, -0.2754995E+03, 0.1250810E+03 } };
285
286 for (int I1=0;I1<=5;I1++){for (int I2=0;I2<=3;I2++){
287 P[I1][I2]=P1[I1][I2];
288 }}
289
290 printf(" our event : P[0,i]= %16.10f %16.10f %16.10f %16.10f \n",P[0][0],P[0][1],P[0][2],P[0][3]);
291 int II=1;
292 printf(" P[1,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]);
293 II=2;
294 printf(" P[2,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]);
295 II=3;
296 printf(" P[3,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]);
297 II=4;
298 printf(" P[4,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]);
299 II=5;
300 printf(" P[4,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]);
301 printf(" ===========================\n");
302//#################################################################
303 KEY=0; //#################################################################
304//#################################################################
305 printf(" ============== KEY= %2i ==\n",KEY);
306 printf(" ===========================\n");
307 printf(" \n");
308 printf(" ===========================================================\n");
309 printf(" ============== non-zero contributions to <|ME|^2>_spin ==\n");
310 printf(" ===========================================================\n");
311 printf(" \n");
312
313 }
314
315 // these loops need to be cleaned from zero contributions!
316 for(int I1=-5;I1<=5;I1++){ // for test of single production process fix flavour
317 for(int I2=-5;I2<=5;I2++){ // for test of single production process fix flavour
318 for(int I3=-5;I3<=5;I3++){
319 for(int I4=-5;I4<=5;I4++){
320
321 int ID1 = I1; if( ID1 == 0 ) ID1 = 21;
322 int ID2 = I2; if( ID2 == 0 ) ID2 = 21;
323 int ID3 = I3; if( ID3 == 0 ) ID3 = 21;
324 int ID4 = I4; if( ID4 == 0 ) ID4 = 21;
325
326 W[0][0] += f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY); // as in case of nonSM_adopt we change the sign for
327 W[0][1] += f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY); // second tau helicity assuming without checking that
328 W[1][0] += f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4,-1, -1, P, KEY); // for VBF quantization frame conventions are the same.
329 W[1][1] += f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4,-1, 1, P, KEY);
330
331 if( DEBUG ) {
332 if (ID1==51 && ID2==-1 && ID3==4 && ID4==-4 ) /// THIS IS BLOCKED NOW #####################################################
333 { KEY=1;
334 double cosik=vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4,-1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4,-1, 1, P, KEY);
335 std::cout << " ID-s = " <<ID1 <<" "<<ID2 <<" "<<ID3 <<" "<<ID4 <<" "<< " KEY="<<KEY<<std::endl;
336 std::cout << " " << std::endl;
337 std::cout << " ME^2 summed over spin = " << cosik << std::endl;
338 std::cout << "---------" << std::endl;
339 std::cout << " " << std::endl;
340 }
341
342 if( ( W[0][0] > 10) || (W[0][1] > 10) || (W[1][0] > 10) || ( W[1][1] > 10)) {
343 std::cout << "ERWxxx: ID= " <<ID1 << " " << ID2 << " " << ID3 << " " << ID4 << " W[]= "
344 << f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY) << " "
345 << f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY) << " "
346 << f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4,-1, -1, P, KEY) << " "
347 << f_x1[I1]*f_x2[I2]*vbfdistr(ID1,ID2,ID3,ID4,-1, 1, P, KEY)
348 << std::endl;
349 }
350 if( vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY) +vbfdistr(ID1,ID2,ID3,ID4, -1, -1, P, KEY)> 0) {
351 float check = vbfdistr(ID1,ID2,ID4,ID3, 1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY);
352 if( check != 0 || check == 0){
353 int CPSTATUSICP=0;
354 printf(" ID-s= %2i %2i %2i %2i CP used= %1i ## VALUE: %16.10e ## Spin contr.: (+-)= %9.3e (-+)= %9.3e (--)= %9.3e (++)= %9.3e \n", ID1, ID2, ID3,ID4, cpstatus_.icp,
355 vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY)+ vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY)+ vbfdistr(ID1,ID2,ID3,ID4, -1, -1, P, KEY)+ vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY),
356 vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY),
357 vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY),
358 vbfdistr(ID1,ID2,ID3,ID4, -1, -1, P, KEY),
359 vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY) );
360 }
361 }
362 }
363
364 }
365 }
366 }
367 }
368
369
370 // std::cout << "That is it"<< std::endl;
371 //exit(-1);
372
373}
374
375
376
377/*******************************************************************************
378 Calculate weights, case of event record vertex like Z/gamma/H ... -> tau tau decay plus 2 jets.
379
380 Determine taus decay channel, calculates all necessary components from decay HHp HHm and production W[][] for
381 calculation of all weights, later calculates weights.
382
383 Input: X four momentum of Z/Higgs may be larger than sum of tau1 tau2, missing component
384 is assumed to be QED brem four momenta of outgoing partons and taus; vectors of decay products for first and second tau
385
386 Hidden input: relWTnonSM, switch of what is WTnonSM
387 Hidden output: WTamplitM or WTamplitP matrix elements of tau1 tau2 decays
388
389
390 Polari - helicity attributed to taus, 100% correlations between tau+ and tau-
391 accesible with getTauSpin()
392 WTnonSM weight for introduction of matrix element due to nonSM in cross section, or just ME^2*PDF's for relWTnonS==false
393 Explicit output: WT spin correlation weight
394*******************************************************************************/
395double calculateWeightFromParticlesVBF(SimpleParticle &p3, SimpleParticle &p4,SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector<SimpleParticle> &sp_tau1_daughters, vector<SimpleParticle> &sp_tau2_daughters)
396{
397 SimpleParticle sp_tau;
398 SimpleParticle sp_nu_tau;
399 vector<SimpleParticle> sp_tau_daughters;
400 int KEY=0;
401 // here we impose that it is for SM Higgs only
402 if(sp_X.pdgid()==25) KEY=1;
403
404 // First iteration is for tau plus, so the 'nu_tau' is tau minus
405 if (sp_tau1.pdgid() == -15 )
406 {
407 sp_tau = sp_tau1;
408 sp_nu_tau = sp_tau2;
409 sp_tau_daughters = sp_tau1_daughters;
410 }
411 else
412 {
413 sp_tau = sp_tau2;
414 sp_nu_tau = sp_tau1;
415 sp_tau_daughters = sp_tau2_daughters;
416 }
417
418 double *HHp, *HHm;
419
420 // We use this to separate namespace for tau+ and tau-
421 if(true)
422 {
423 // Create Particles from SimpleParticles
424 Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
425 Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
426 Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
427
428 vector<Particle> tau_daughters;
429
430 // tau pdgid
431 int tau_pdgid = sp_tau.pdgid();
432
433 // Create list of tau daughters
434 for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
435 {
436 Particle pp(sp_tau_daughters[i].px(),
437 sp_tau_daughters[i].py(),
438 sp_tau_daughters[i].pz(),
439 sp_tau_daughters[i].e(),
440 sp_tau_daughters[i].pdgid() );
441
442 tau_daughters.push_back(pp);
443 }
444
445 double phi2 = 0.0, theta2 = 0.0;
446
447
448 // Move decay kinematics first to tau rest frame with z axis pointing along nu_tau direction
449 // later rotate again to have neutrino from tau decay along z axis: angles phi2, theta2
450 prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
451
452
453 // Determine decay channel and then calculate polarimetric vector HH
454 HHp = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
455
456 DEBUG
457 (
458 cout<<tau_pdgid<<" -> ";
459 for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
460 cout<<" (HHp: "<<HHp[0]<<" "<<HHp[1]<<" "<<HHp[2]<<" "<<HHp[3]<<") ";
461 cout<<endl;
462 )
463
464 WTamplitP = WTamplit;
465 } // end of tau+
466
467 // Second iteration is for tau minus, so the 'nu_tau' is tau minus
468 if(sp_tau1.pdgid() == 15 )
469 {
470 sp_tau = sp_tau1;
471 sp_nu_tau = sp_tau2;
472 sp_tau_daughters = sp_tau1_daughters;
473 }
474 else
475 {
476 sp_tau = sp_tau2;
477 sp_nu_tau = sp_tau1;
478 sp_tau_daughters = sp_tau2_daughters;
479 }
480
481 // We use this to separate namespace for tau+ and tau-
482 if(true)
483 {
484 // Create Particles from SimpleParticles
485 Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
486 Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
487 Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
488
489 vector<Particle> tau_daughters;
490
491 // tau pdgid
492 int tau_pdgid = sp_tau.pdgid();
493
494 // Create list of tau daughters
495 for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
496 {
497 Particle pp(sp_tau_daughters[i].px(),
498 sp_tau_daughters[i].py(),
499 sp_tau_daughters[i].pz(),
500 sp_tau_daughters[i].e(),
501 sp_tau_daughters[i].pdgid() );
502
503 tau_daughters.push_back(pp);
504 }
505
506 double phi2 = 0.0, theta2 = 0.0;
507
508
509 // Move decay kinematics first to tau rest frame with z axis pointing along nu_tau direction
510 // later rotate again to have neutrino from tau decay along z axis: angles phi2, theta2
511 prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
512
513
514 // Determine decay channel and then calculate polarimetric vector HH
515 HHm = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
516
517 DEBUG
518 (
519 cout<<tau_pdgid<<" -> ";
520 for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
521 cout<<" (HHm: "<<HHm[0]<<" "<<HHm[1]<<" "<<HHm[2]<<" "<<HHm[3]<<") ";
522 cout<<endl;
523 )
524
525 WTamplitM = WTamplit;
526 } // end of tau-
527
528
529 double W[2][2] = { { 0.25, 0.25 },
530 { 0.25, 0.25 } }; // this is trivial W spin (helicity only) density matrix
531 // consist of ME^2*PDF's for production.
532 // pre-initialization all are equal i.e. no spin effects
533
534
535 // calculate ME^2*PDF's for SM
536 // we use: sp_nu_tau, sp_tau; funny names - object order for HH calc.
537 getME2VBF(p3, p4, sp_X, sp_nu_tau, sp_tau, W, KEY);
538 // if(KEY==0 || KEY==1 )getME2VBF(p3, p4, sp_X, sp_tau1, sp_tau2, W, KEY);
539 //else getME2VBF(p3, p4, sp_X, sp_tau1, sp_tau2, W, 0);
540
541
542 double sum=(W[0][0]+W[0][1]+ W[1][0]+W[1][1]); // getME2VBF calculated PDF's times matrix elements squared for each of four helicity configurations
543 // tau+ and tau- using SM process KEY=0 DY, KEY=1 Higgs.
544
545 if(nonSM2==0 ) {
546 WTnonSM=1.0;
547 if(relWTnonSM==0) WTnonSM=sum; // The variable accessible for user stores production weight ME^2 * PDF's summed over flavours and spin states
548 }
549 if(nonSM2==1) {
550 // now we re-calculate W using anomalous variant of matrix elements.
551 // we use: sp_nu_tau, sp_tau; funny names - object order for HH calc.
552 getME2VBF(p3, p4, sp_X, sp_nu_tau, sp_tau, W, KEY+2);
553
554 double sum2=(W[0][0]+W[0][1]+ W[1][0]+W[1][1]);
555
556 WTnonSM=sum2/sum;
557 if(relWTnonSM==0) WTnonSM=sum2;
558 }
559
560
561 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]);
562 WT = WT/(W[0][0]+W[0][1]+ W[1][0]+W[1][1]);
563
564 // we separate cross section into first tau helicity + and - parts. Second tau follow.
565 // This must be tested especially for KEY >=2
566 // case for scalar (H) has flipped sign of second tau spin? Tests needed
567 double RRR = Tauola::randomDouble();
568 Polari=1.0;
569 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;
570
571
572
573 // Print out some info about the channel
574 DEBUG( cout<<" WT: "<<WT<<endl; )
575
576 if (WT<0.0) {
577 printf("WT is: %13.10f. Setting WT = 0.0\n",WT);
578 WT = 0.0;
579 }
580
581 if (WT>4.0) {
582 printf("WT is: %13.10f. Setting WT = 4.0\n",WT);
583 WT = 4.0;
584 }
585
586 delete[] HHp;
587 delete[] HHm;
588
589 return WT;
590}
591
592} // namespace TauSpinner
double vbfdistr_(int *I1, int *I2, int *I3, int *I4, int *H1, int *H2, double P[6][4], int *KEY)
void set_vbfdistrModif(double(*function)(int, int, int, int, int, int, double[6][4], int, double))
Set vbfdistrModif function.
Definition vbfdistr.cxx:46
void setPDFOpt(int QCDdefault, int QCDvariant)
Definition vbfdistr.cxx:51
void set_alphasModif(void(*function)(double, int, int))
Set vbfdistrModif function.
Definition vbfdistr.cxx:33
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
void(* alphasModif)(double, int, int)
Pointer to a function that modify rezult of alpha_ calc. in vbfdist.
Definition vbfdistr.cxx:30
void vbf_reinit_(int *key)
double vbfdistr(int I1, int I2, int I3, int I4, int H1, int H2, SimpleParticle &p1, SimpleParticle &p2, SimpleParticle &tau1, SimpleParticle &tau2, int KEY)
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
void getME2VBF(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X, SimpleParticle &tau1, SimpleParticle &tau2, double(&W)[2][2], int KEY)
Definition vbfdistr.cxx:118
double(* vbfdistrModif)(int, int, int, int, int, int, double[6][4], int, double)
Pointer to a function that modify rezult of Matrix Element calc. in vbfdist.
Definition vbfdistr.cxx:43