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>
6 using namespace Tauolapp;
7 
8 namespace TauSpinner {
9 
10 extern double CMSENE;
11 extern int nonSM2;
12 extern int relWTnonSM;
13 extern double WTnonSM;
14 extern double Polari;
15 extern double WTamplit;
16 extern double WTamplitP;
17 extern double WTamplitM;
18 extern double CMSENE;
19 extern double f(double x, int ID, double SS, double cmsene);
20 
21 const 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  */
30 void (*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  */
43 double (*vbfdistrModif)(int, int, int, int, int, int, double[6][4], int, double) = NULL;
44 
45 /** @brief Set vbfdistrModif function */
46 void set_vbfdistrModif(double (*function)(int, int, int, int, int, int, double[6][4], int, double) )
47 {
48  vbfdistrModif = function;
49 }
50 
51 void 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*/
87 double 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 //---------------------------------------------------------------------------------------------------------------------------------------------------
118 void 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 *******************************************************************************/
395 double 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
TauSpinner::vbfdistr_
double vbfdistr_(int *I1, int *I2, int *I3, int *I4, int *H1, int *H2, double P[6][4], int *KEY)
TauSpinner::set_vbfdistrModif
void set_vbfdistrModif(double(*function)(int, int, int, int, int, int, double[6][4], int, double))
Set vbfdistrModif function.
Definition: vbfdistr.cxx:46
TauSpinner::setPDFOpt
void setPDFOpt(int QCDdefault, int QCDvariant)
Definition: vbfdistr.cxx:51
TauSpinner::set_alphasModif
void set_alphasModif(void(*function)(double, int, int))
Set vbfdistrModif function.
Definition: vbfdistr.cxx:33
TauSpinner::calculateHH
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
Definition: tau_reweight_lib.cxx:1082
TauSpinner::alphasModif
void(* alphasModif)(double, int, int)
Pointer to a function that modify rezult of alpha_ calc. in vbfdist.
Definition: vbfdistr.cxx:30
TauSpinner::vbf_reinit_
void vbf_reinit_(int *key)
TauSpinner::vbfdistr
double vbfdistr(int I1, int I2, int I3, int I4, int H1, int H2, SimpleParticle &p1, SimpleParticle &p2, SimpleParticle &tau1, SimpleParticle &tau2, int KEY)
TauSpinner::prepareKinematicForHH
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
Definition: tau_reweight_lib.cxx:985
TauSpinner::getME2VBF
void getME2VBF(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X, SimpleParticle &tau1, SimpleParticle &tau2, double(&W)[2][2], int KEY)
Definition: vbfdistr.cxx:118
TauSpinner::vbfdistrModif
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
Generated by doxygen 1.9.1