C++ Interface to Tauola
tau_reweight_lib.cxx
1 #include "TauSpinner/tau_reweight_lib.h"
2 #include "TauSpinner/nonSM.h"
3 #include "TauSpinner/Tauola_wrapper.h"
4 #include "TauSpinner/EWtables.h"
5 #include <complex>
6 #include <stdlib.h>
7 #include "Tauola/TauolaParticlePair.h"
8 using namespace Tauolapp;
9 namespace TauSpinner {
10 
11 /***** GLOBAL VARIABLES AND THEIR PRE-INITIALISATION *****/
12 
13 double CMSENE = 7000.0;// Center of mass system energy (used by PDF's)
14 bool Ipp = true; // pp collisions
15 int Ipol = 0; // Is the sample in use polarized? (relevant for Z/gamma case only)
16 int nonSM2 = 1; // Turn on/off nonSM calculations
17 int nonSMN = 0; // Turn on, if calculations of nonSM weight, is to be used to modify shapes only
18 int relWTnonSM = 1; // 1: relWTnonSM is relative to SM; 0: absolute
19 
20 double WTnonSM=1.0; // nonSM weight
21 int keyGSW = 1; // variant of electroweak formfactor treatment, complete or partial use.
22 int ifkorch=1; // key to switch GSW+anomalous dipoles code and calculation of R_ij
23  // with dipolqq_(). Then this routine is used both for nonSM2=0 or 1;
24  // for 0 ew-forfactors are set to 1 and Amz0,Gamz0,sin2W0
25 
26 int ifGSW=1; // key to switch electroweak-formfactors, works only with ifkorch=1
27 int FrameType=0; // key to switch (1) Mustraal (NLO improving) (0) Collins-Soper like, frames in evaluation of born kinematics.
28 
29 double ReAini=0.0; // anomalous dipole moments magnetic/electric (A/B) for gamma* and for Z (X/Y)
30 double ImAini=0.0; // also work for ifkorch=1 only
31 double ReBini=0.0;
32 double ImBini=0.0;
33 double ReXini=0.0;
34 double ImXini=0.0;
35 double ReYini=0.0;
36 double ImYini=0.;
37  double GSWN1r = 1.0;
38  double GSWN1i = 0.0;
39  double GSWN2r = 1.0;
40  double GSWN2i = 0.0;
41  double GSWN3r = 1.0;
42  double GSWN3i = 0.0;
43  double GSWN4r = 1.0;
44  double GSWN4i = 0.0;
45  double GSWN5r = 1.0;
46  double GSWN5i = 0.0;
47  double GSWN6r = 1.0;
48  double GSWN6i = 0.0;
49  double GAMfraci=0.0;
50  double GAMfrac2i=0.0; // ratio of gamma gamma process to the rest, in generated sample and in rewieght
51  double A0i=0.0; // anomalous couplings for gamma gamma in sample
52  double B0i=0.0;
53  double Ai=0.0; // in reweight
54  double Bi=0.0;
55  int iqed = 0; // 1/0 QED part of dipole mom. A on/off
56 
57 double Polari =0.0; // Helicity, attributed to the tau pair. If not attributed then 0.
58  // Polari is attributed for the case when spin effects are taken into account.
59  // Polari is available for the user program with getTauSpin()
60 bool IfHiggs=false; // flag used in sigborn()
61 double WTamplit=1; // AMPLIT weight for the decay last invoked
62 double WTamplitP=1; // AMPLIT weight for tau+ decay
63 double WTamplitM=1; // AMPLIT weight for tau- decay
64 
65 // Higgs parameters
66 int IfHsimple=0; // switch for simple Higgs (just Breit-Wigner)
67 double XMH = 125.0; // mass (GeV)
68 double XGH = 1.0; // width, should be detector width, analysis dependent.
69 double Xnorm = 0.15; // normalization of Higgs Born cross-section, at hoc
70 
71 // Transverse components of Higgs spin density matrix
72 double RXX = 0.0; //-1.0;
73 double RYY = 0.0; // 1.0;
74 double RXY = 0.0;
75 double RYX = 0.0;
76 
77 // Coefficients for transverse components of Z/gamma spin density matrix
78 double RzXX = 0.0; //-1.0;
79 double RzYY = 0.0; // 1.0;
80 double RzXY = 0.0;
81 double RzYX = 0.0;
82 
83 // Initial values of logitudinal/transverse components of spin density matrix to be calculated inside getLongitudinalPolarization
84 double R11 = 0.0;
85 double R12 = 0.0;
86 double R21 = 0.0;
87 double R22 = 0.0;
88 double R31 = 0.0;
89 double R32 = 0.0;
90 double R33 = 0.0;
91 double R41 = 0.0;
92 double R42 = 0.0;
93 double R43 = 0.0;
94 double R44 = 0.0;
95 // R34 R43 is not used
96 
97 //Tauolapp::TauolaParticlePair pp;
98  double RcorExt[4][4];
99 /***** END: GLOBAL VARIABLES AND THEIR PRE-INITIALISATION *****/
100 
101 
102 double f(double x, int ID, double SS, double cmsene)
103 // PDF's parton density function divided by x;
104 // x - fraction of parton momenta
105 // ID flavour of incoming quark
106 // SS scale of hard process
107 // cmsene center of mass for pp collision.
108 {
109  // LHAPDF manual: http://projects.hepforge.org/lhapdf/manual
110  // double xfx(const double &x;, const double &Q;, int fl);
111  // returns xf(x, Q) for flavour fl - this time the flavour encoding
112  // is as in the LHAPDF manual...
113  // -6..-1 = tbar,...,ubar, dbar
114  // 1..6 = duscbt
115  // 0 = g
116  // xfx is the C++ wrapper for fortran evolvePDF(x,Q,f)
117 
118  // note that SS=Q^2, make the proper choice of PDFs arguments.
119  return LHAPDF::xfx(x, sqrt(SS), ID)/x;
120 
121  //return x*(1-x);
122 
123 }
124 
125  // Calculates Born cross-section summed over final taus spins.
126  // Input parameters:
127  // incoming flavour ID
128  // invariant mass^2 SS
129  // scattering angle costhe
130  // Hidden input (type of the process): IfHiggs, IfHsimple
131 double sigborn(int ID, double SS, double costhe)
132 {
133  // cout << "ID : " << ID << " HgsWTnonSM = " << HgsWTnonSM << " IfHsimple = " << IfHsimple << endl;
134  // BORN x-section.
135  // WARNING: overall sign of costheta must be fixed
136  int tauID = 15;
137 
138  // case of the Higgs boson
139  if (IfHiggs) {
140  double SIGggHiggs=0.;
141  // for the time being only for gluon it is non-zero.
142  if(ID==0){
143  int IPOne = 1;
144  int IMOne = -1;
145  SIGggHiggs=disth_(&SS, &costhe, &IPOne, &IPOne)+disth_(&SS, &costhe, &IPOne, &IMOne)+
146  disth_(&SS, &costhe, &IMOne, &IPOne)+disth_(&SS, &costhe, &IMOne, &IMOne);
147 
148 
149  double PI=3.14159265358979324;
150  SIGggHiggs *= XMH * XMH * XMH *XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) + XGH*XGH*XMH*XMH);
151  // cout << "JK: SIGggHiggs = " << SS << " " << XMH << " " << XGH << " " << XMH * XMH * XMH * XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) + XGH*XGH*XMH*XMH) << " " << SIGggHiggs << endl;
152 
153  if(IfHsimple==1) SIGggHiggs = Xnorm * XMH * XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) +XGH*XGH*XMH*XMH);
154  // cout << "ZW: SIGggHiggs = " << SS << " " << costhe << " " << SIGggHiggs << endl;
155  }
156  return SIGggHiggs;
157  }
158 
159  // case of Drell-Yan
160 
161  if (ID==0) return 0.0 ; // for the time being for gluon it is zero.
162  if (ID*ID==22*22){
163  int iZero = 0;
164  double dOne = 1.0;
165  double dMOne = -1.0;
166  // ERW:: Genis
167  // sum DY Born over all tau helicity configurations:
168  // return ( t_gamm_(&iZero, &SS, &costhe, &dOne , &dOne) + t_gamm_(&iZero, &SS, &costhe, &dOne , &dMOne)
169  // + t_gamm_(&iZero, &SS, &costhe, &dMOne, &dOne) + t_gamm_(&iZero, &SS, &costhe, &dMOne, &dMOne))/SS/128.86674175/128.86674175/3.141559265358979324; //123231.;
170 // overall norm. factor .../SS/123231 most probably it is alpha_QED**2/pi/2/SS is from comparison between Born we use and Born used in Warsaw group.
171  }
172  if (ID>0) initwk_( &ID, &tauID, &SS);
173  else
174  {
175  ID = -ID;
176  initwk_( &ID, &tauID, &SS);
177  }
178 
179  int iZero = 0;
180  double dOne = 1.0;
181  double dMOne = -1.0;
182  // sum DY Born over all tau helicity configurations:
183  return ( t_born_(&iZero, &SS, &costhe, &dOne , &dOne) + t_born_(&iZero, &SS, &costhe, &dOne , &dMOne)
184  + t_born_(&iZero, &SS, &costhe, &dMOne, &dOne) + t_born_(&iZero, &SS, &costhe, &dMOne, &dMOne))/SS/128.86674175/128.86674175/3.141559265358979324; //123231.;
185 // overall norm. factor .../SS/123231 most probably it is alpha_QED**2/pi/2/SS is from comparison between Born we use and Born used in Warsaw group.
186 }
187 
188 /*******************************************************************************
189  Initialize anomalous dipole moments 8=2*2*2 cnstants:
190  magnetic/electric real/imaginary parts
191  for virtual photon/ for Z,
192  and its switches:
193  ifGSW electroweak form-factors should be on for Korchin calculations?
194  ifkorch if calculation by Korchin et al. should be used or not.
195 *******************************************************************************/
196 void initialize_GSW( int _ifGSW, int _ifkorch, int _iqed,
197  double _ReAini,
198  double _ImAini,
199  double _ReBini,
200  double _ImBini,
201  double _ReXini,
202  double _ImXini,
203  double _ReYini,
204  double _ImYini){
205  ifGSW =_ifGSW;
206  ifkorch=_ifkorch;
207  iqed=_iqed;
208  ReAini=_ReAini;
209  ImAini=_ImAini;
210  ReBini=_ReBini;
211  ImBini=_ImBini;
212  ReXini=_ReXini;
213  ImXini=_ImXini;
214  ReYini=_ReYini;
215  ImYini=_ImYini;
216 
217  }
218 
219 /*******************************************************************************
220  Initialize anomalous gamma gamma couplings target A and B, in sample A0i B0i
221  if sample is of no anonmalous couplings set A0i B0i to zero.
222 Also how much of gamma gamma process there should be with respect to quarks
223 in sample GAMfraci in target GAMfrac2i
224 *******************************************************************************/
225  void initialize_gamagama( double _GAMfraci,
226  double _GAMfrac2i,
227  double _A0i,
228  double _B0i,
229  double _Ai,
230  double _Bi
231  ){
232  GAMfraci=_GAMfraci;
233  GAMfrac2i=_GAMfrac2i;
234  A0i= _A0i;
235  B0i=_B0i;
236  Ai=_Ai;
237  Bi=_Bi;
238  }
239 /*************************************************************
240  Provided by user complex numbers (ReGSW_i,ImGSW_i) will
241  multiply form factors.
242  keyGSW will provide further changes, compatible to what
243  was used in work on effective-improved Born paper
244 *************************************************************/
245  void initialize_GSW_norm( int _keyGSW,
246  double _ReGSW1,
247  double _ImGSW1,
248  double _ReGSW2,
249  double _ImGSW2,
250  double _ReGSW3,
251  double _ImGSW3,
252  double _ReGSW4,
253  double _ImGSW4,
254  double _ReGSW5,
255  double _ImGSW5,
256  double _ReGSW6,
257  double _ImGSW6){
258 
259  keyGSW = _keyGSW;
260 
261  GSWN1r = _ReGSW1;
262  GSWN1i = _ImGSW1;
263  GSWN2r = _ReGSW2;
264  GSWN2i = _ImGSW2;
265  GSWN3r = _ReGSW3;
266  GSWN3i = _ImGSW3;
267  GSWN4r = _ReGSW4;
268  GSWN4i = _ImGSW4;
269  GSWN5r = _ReGSW5;
270  GSWN5i = _ImGSW5;
271  GSWN6r = _ReGSW6;
272  GSWN6i = _ImGSW6;
273  }
274 /******************************************************
275 First we multiply form-factors by complex numbers,
276 then we go after keyGSW variants.
277 ******************************************************/
278  void GSWadapt(int ID,double S, double cost, double GSWr[7],double GSWi[7]){
279 
280  double SS=S;
281  double SSS=S;
282  double costhe=cost;
283 
284  if(keyGSW!=1){
285  double AMZi= Amz(ID);
286  costhe=0.0;SSS=AMZi*AMZi;
287  GSWr[0]= 1.0;
288  GSWr[1]= 1.0;
289  GSWr[2]= 1.0;
290  GSWr[3]= 1.0;
291  GSWr[5]= 1.0;
292 
293  GSWi[0]= 0.0; //imGSW[1];
294  GSWi[1]= 0.0;
295  GSWi[2]= 0.0;
296  GSWi[3]= 0.0;
297  GSWi[5]= 0.0;
298 
299  }
300  if(keyGSW==1){ //options do nothing
301  // GSWr[5] = real(EWFACT(ID,5,SS,costhe)); // reGSW[6]; //part only
302  // GSWi[5] = imag(EWFACT(ID,5,SS,costhe)); //imGSW[6]; //part only
303  // GSWi[5] = imag(EWFACT(ID,6,SS,costhe)); //imGSW[6];
304  }
305  else if(keyGSW==3){
306  GSWr[0]= 1.0; // reGSW[1];
307  GSWr[1]= 1.0;
308  GSWr[2]= 1.0;
309  GSWr[3]= 1.0;
310  GSWr[5]= real(EWFACT(ID,6,SSS,costhe)); // reGSW[6];
311 
312  GSWi[0]= 0.0; //imGSW[1];
313  GSWi[1]= 0.0;
314  GSWi[2]= 0.0;
315  GSWi[3]= 0.0;
316  GSWi[5]= imag(EWFACT(ID,6,SSS,costhe)); //imGSW[6];
317 
318  }
319  else if(keyGSW==4){
320  GSWr[0] = real(EWFACT(ID,0,SSS,costhe)); // reGSW[1];
321  GSWr[1] = 1.0;
322  GSWr[2] = 1.0;
323  GSWr[3] = 1.0;
324  GSWr[5] = real(EWFACT(ID,6,SSS,costhe)); // reGSW[6];
325 
326  GSWi[0] = imag(EWFACT(ID,0,SSS,costhe)); //imGSW[1];
327  GSWi[1] = 0.0;
328  GSWi[2] = 0.0;
329  GSWi[3] = 0.0;
330  GSWi[5] = imag(EWFACT(ID,6,SSS,costhe)); //imGSW[6];
331 
332  }
333  else if(keyGSW==5){
334  GSWr[0]= real(EWFACT(ID,0,SSS,costhe)); // reGSW[1];
335  GSWr[1]= real(EWFACT(ID,1,SSS,costhe)); // reGSW[2];
336  GSWr[2]= real(EWFACT(ID,2,SSS,costhe)); // reGSW[3];
337  GSWr[3]= real(EWFACT(ID,3,SSS,costhe)); //1.0;
338  GSWr[5]= real(EWFACT(ID,6,SSS,costhe)); // reGSW[6];
339 
340  GSWi[0]= imag(EWFACT(ID,0,SSS,costhe)); //imGSW[1];
341  GSWi[1]= imag(EWFACT(ID,1,SSS,costhe)); //imGSW[2];
342  GSWi[2]= imag(EWFACT(ID,2,SSS,costhe)); //imGSW[3];
343  GSWi[3]= imag(EWFACT(ID,3,SSS,costhe)); //0.0;
344  GSWi[5]= imag(EWFACT(ID,6,SSS,costhe)); //imGSW[6];
345 
346  /*
347  // real v, vv
348  ImGSW2 = 0.0;
349  ImGSW3 = 0.0;
350  ImGSW4 = 0.0;
351  */
352  /*
353  // real v, vv=1
354  ReGSW4 = 1.0;
355  ImGSW2 = 0.0;
356  ImGSW3 = 0.0;
357  ImGSW4 = 0.0;
358  */
359  /*
360  // real pigamgam
361  ReGSW4 = 1.0;
362  ImGSW1 = 0.0;
363  ImGSW2 = 0.0;
364  ImGSW3 = 0.0;
365  ImGSW4 = 0.0;
366  ImGSW6 = 0.0; // make Pi_gammgamma real: formally it is higher order. But
367  // see bardin+grunewald+pasarino-9902452.pdf
368  // why it is inconsistent with the LEP choice
369  */
370  /*
371  // LEP2005 style
372  ReGSW4 = 1.0;
373  ImGSW1 = 0.0;
374  ImGSW2 = 0.0;
375  ImGSW3 = 0.0;
376  ImGSW4 = 0.0;
377  */
378  }
379  else if(keyGSW==12){
380  GSWr[0]= real(EWFACT(0,0,SSS,costhe)); // for (v1)
381  }
382  else if(keyGSW==13){
383  GSWr[0]=real(EWFACT(ID,0,SSS,costhe)); // for (v2)
384  GSWr[1]=real(EWFACT(ID,1,SSS,costhe))/real(EWFACT(0,1,SSS,costhe)); // for (v2)
385  GSWr[2]=real(EWFACT(ID,2,SSS,costhe))/real(EWFACT(0,2,SSS,costhe)); // for (v2)
386  }
387 
388  // rotations and size change:
389  GSWr[0] = GSWr[0]*GSWN1r- GSWi[0]*GSWN1i;
390  GSWi[0] = GSWr[0]*GSWN1i+ GSWi[0]*GSWN1r;
391  GSWr[1] = GSWr[1]*GSWN2r- GSWi[1]*GSWN2i;
392  GSWi[1] = GSWr[1]*GSWN2i+ GSWi[1]*GSWN2r;
393  GSWr[2] = GSWr[2]*GSWN3r- GSWi[2]*GSWN3i;
394  GSWi[2] = GSWr[2]*GSWN3i+ GSWi[2]*GSWN3r;
395  GSWr[3] = GSWr[3]*GSWN4r- GSWi[3]*GSWN4i;
396  GSWi[3] = GSWr[3]*GSWN4i+ GSWi[3]*GSWN4r;
397  GSWr[4] = GSWr[4]*GSWN5r- GSWi[4]*GSWN5i;
398  GSWi[4] = GSWr[4]*GSWN5i+ GSWi[4]*GSWN5r;
399  GSWr[5] = GSWr[5]*GSWN6r- GSWi[5]*GSWN6i;
400  GSWi[5] = GSWr[5]*GSWN6i+ GSWi[5]*GSWN6r;
401 
402  }
403 
404  complex<double> Frezu(int NO){
405  complex<double> F;
406  if( NO==0 )
407  { F=complex<double>(GSWN1r,GSWN1i);}
408  else if(NO==1 )
409  { F=complex<double>(GSWN2r,GSWN2i);}
410  else if(NO==2 )
411  { F=complex<double>(GSWN3r,GSWN3i);}
412  else if(NO==3 )
413  { F=complex<double>(GSWN4r,GSWN4i);}
414  else if(NO==4 )
415  { F=complex<double>(GSWN5r,GSWN5i);}
416  else if(NO==5 ) // leptonic alpha_qed
417  { F=complex<double>(GSWN6r,GSWN6i);}
418  else if(NO==6 ) // complete alpha_qed
419  { F=complex<double>(GSWN6r,GSWN6i);}
420  else
421  {cout<<"wrong NO="<<NO<<endl; exit(-1);}
422  return F;
423  }
424 
425 /*******************************************************************************
426  Initialize TauSpinner
427 
428  Print info and set global variables
429 *******************************************************************************/
430 void initialize_spinner(bool _Ipp, int _Ipol, int _nonSM2, int _nonSMN, double _CMSENE)
431 {
432  Ipp = _Ipp;
433  Ipol = _Ipol;
434  nonSM2 = _nonSM2;
435  nonSMN = _nonSMN;
436 
437  CMSENE = _CMSENE;
438 
439  cout<<" ------------------------------------------------------"<<endl;
440  cout<<" TauSpinner v2.1.2"<<endl;
441  cout<<" -----------------"<<endl;
442  cout<<" 17.Feb.2025 "<<endl;
443  cout<<" by Z. Czyczula (until 2015), T. Przedzinski, E. Richter-Was, Z. Was,"<<endl;
444  cout<<" matrix elements implementations "<<endl;
445  cout<<" also J. Kalinowski, W. Kotlarski and M. Bachmani"<<endl;
446  cout<<" ------------------------------------------------------"<<endl;
447  cout<<" Ipp - true for pp collision; otherwise polarization"<<endl;
448  cout<<" of individual taus from Z/gamma* is set to 0.0"<<endl;
449  cout<<" Ipp = "<<Ipp<<endl;
450  cout<<" CMSENE - used in PDF calculations; only if Ipp = true"<<endl;
451  cout<<" and only for Z/gamma*"<<endl;
452  cout<<" CMSENE = "<<CMSENE<<endl;
453  cout<<" Ipol - relevant for Z/gamma* decays "<<endl;
454  cout<<" 0 - events generated without spin effects "<<endl;
455  cout<<" 1 - events generated with all spin effects "<<endl;
456  cout<<" 2 - events generated with spin correlations and <pol>=0 "<<endl;
457  cout<<" 3 - events generated with spin correlations and"<<endl;
458  cout<<" polarization but missing angular dependence of <pol>"<<endl;
459  cout<<" Ipol = "<<Ipol<<endl;
460  cout<<" Ipol - relevant for Z/gamma* decays "<<endl;
461  cout<<" NOTE: For Ipol=0,1 algorithm is identical. "<<endl;
462  cout<<" However in user program role of wt need change. "<<endl;
463  cout<<" nonSM2 = "<<nonSM2<<endl;
464  cout<<" 1/0 extra term in cross section, density matrix on/off "<<endl;
465  cout<<" nonSMN = "<<nonSMN<<endl;
466  cout<<" 1/0 extra term in cross section, for shapes only? on/off "<<endl;
467  cout<<" note KEY - for options of matrix elements calculations "<<endl;
468  cout<<" in cases of final states with two jets "<<endl;
469  cout<<" ------------------------------------------------------ "<<endl;
470 }
471 /*******************************************************************************
472  Set flag for type of effective Born kinematic approximation.
473  1: Mustraal frame (NLO like)
474  0: old default, similar to Collins-Soper.
475 *******************************************************************************/
476  void setFrameType(int _FrameType)
477 {
478  FrameType = _FrameType;
479 }
480 
481 /*******************************************************************************
482  Set flag for calculating relative(NONSM-SM)/absolute weight for X-section
483  calculated as by product in longitudinal polarization method.
484  1: relWTnonSM is relative to SM (default)
485  0: absolute
486 *******************************************************************************/
487 void setRelWTnonSM(int _relWTnonSM)
488 {
489  relWTnonSM = _relWTnonSM;
490 }
491 
492 /*******************************************************************************
493  Set Higgs mass, width and normalization of Higgs born function
494  Default is mass = 125, width = 1.0, normalization = 0.15
495 *******************************************************************************/
496  void setHiggsParameters(int jak, double mass, double width, double normalization)
497 {
498  IfHsimple=jak;
499  XMH = mass;
500  XGH = width;
501  Xnorm = normalization;
502 }
503 
504 /*******************************************************************************
505  Set transverse components of Higgs spin density matrix
506 *******************************************************************************/
507 void setHiggsParametersTR(double Rxx, double Ryy, double Rxy, double Ryx)
508 {
509 
510  RXX = Rxx;
511  RYY = Ryy;
512  RXY = Rxy;
513  RYX = Ryx;
514 }
515 
516 
517 /************************************************************************************
518  Get longitudinal (or time) - transverse components of Z/gamma spin density matrix
519 *************************************************************************************/
520 void getZgamParametersL(double &Rzx, double &Rzy, double &Rtx, double &Rty, double &Rtz)
521 {
522  Rzx = R31;
523  Rzy = R32;
524  Rtx = R41;
525  Rty = R42;
526  Rtz = R43;
527 }
528 
529 /*******************************************************************************
530  Set coefficients for transverse components of Z/gamma spin density matrix multipliers
531 *******************************************************************************/
532 void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx)
533 {
534  RzXX = Rxx;
535  RzYY = Ryy;
536  RzXY = Rxy;
537  RzYX = Ryx;
538 }
539 
540 /*******************************************************************************
541  Get transverse components of Z/gamma spin density matrix
542 *******************************************************************************/
543 void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
544 {
545  Rxx = R11;
546  Ryy = R22;
547  Rxy = R12;
548  Ryx = R21;
549 }
550 /*******************************************************************************
551  Get transverse components of Z/gamma spin density matrix multipliers
552 *******************************************************************************/
553 void getZgamMultipliersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
554 {
555  Rxx = RzXX;
556  Ryy = RzYY;
557  Rxy = RzXY;
558  Ryx = RzYX;
559 }
560 
561 
562 
563 /*******************************************************************************
564  Get Higgs mass, width and normalization of Higgs born function
565 *******************************************************************************/
566 void getHiggsParameters(double *mass, double *width, double *normalization)
567 {
568  *mass = XMH;
569  *width = XGH;
570  *normalization = Xnorm;
571 }
572 
573 /*******************************************************************************
574  Set type of spin treatment used in the sample
575  Ipol = 0 sample was not polarised
576  Ipol = 1 sample was having complete longitudinal spin effects
577  Ipol = 2 sample was featuring longitudinal spin correlations only,
578  but not dependence on polarisation due to couplings of the Z
579  Ipol = 3 as in previous case, but only angular dependence of spin polarisation
580  was missing in the sample
581 *******************************************************************************/
582 void setSpinOfSample(int _Ipol)
583 {
584  Ipol = _Ipol;
585 }
586 
587 /*******************************************************************************
588  Turn nonSM calculation of Born cross-section on/off
589 *******************************************************************************/
590 void setNonSMkey(int _key)
591 {
592  nonSM2 = _key;
593 }
594 
595 /*******************************************************************************
596  Get nonSM weight
597 *******************************************************************************/
598 double getWtNonSM()
599 {
600  return WTnonSM;
601 }
602 /*******************************************************************************
603  Get weights for tau+ tau- decay matrix elements
604 *******************************************************************************/
605 double getWtamplitP(){return WTamplitP;}
606 double getWtamplitM(){return WTamplitM;}
607 
608 /*******************************************************************************
609  Get tau spin (helicities of tau+ tau- are 100% correlated)
610  Use after sample is reweighted to obtain information on attributed tau
611  longitudinal spin projection.
612 *******************************************************************************/
613 double getTauSpin()
614 {
615  return Polari;
616 }
617 
618 /*******************************************************************************
619  Calculate weights, case of event record vertex like W -> tau nu_tau decay.
620  Function for W+/- and H+/-
621 
622  Determines decay channel, calculates all necessary components for
623  calculation of all weights, calculates weights.
624  Input: X four momentum may be larger than sum of tau nu_tau, missing component
625  is assumed to be QED brem
626  Hidden input: none
627  Hidden output: WTamplitM or WTamplitP of tau decay (depending on tau charge)
628  NOTE: weight for sp_X production matrix elements is not calculated
629  for decays of charged intermediate W+-/H+-!
630  Explicit output: WT spin correlation weight
631 *******************************************************************************/
632 double calculateWeightFromParticlesWorHpn(SimpleParticle &sp_X, SimpleParticle &sp_tau, SimpleParticle &sp_nu_tau, vector<SimpleParticle> &sp_tau_daughters)
633 {
634  // Create Particles from SimpleParticles
635 
636  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
637  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
638  Particle nu_tau(sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
639 
640  vector<Particle> tau_daughters;
641 
642  // tau pdgid
643  int tau_pdgid = sp_tau.pdgid();
644 
645  // Create vector of tau daughters
646  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
647  {
648  Particle pp(sp_tau_daughters[i].px(),
649  sp_tau_daughters[i].py(),
650  sp_tau_daughters[i].pz(),
651  sp_tau_daughters[i].e(),
652  sp_tau_daughters[i].pdgid() );
653 
654  tau_daughters.push_back(pp);
655  }
656 
657  double phi2 = 0.0, theta2 = 0.0;
658 
659  // To calcluate matrix elements TAUOLA need tau decay products in tau rest-frame: we first boost all products
660  // to tau rest frame (tau nu_tau of W decay along z-axis, intermediate step for boost is tau nu_tau (of W decay) rest-frame),
661  // then rotate to have neutrino from tau decay along z axis;
662  // calculated for that purpose angles phi2, theta2 are stored for rotation back of HH
663  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
664 
665 
666  // Identify decay channel and then calculate polarimetric vector HH; calculates also WTamplit
667  double *HH = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
668 
669  double sign = 1.0; // tau from W is 100 % polarized, also from charged Higgs (but with opposite sign)
670  if ( abs(sp_X.pdgid()) == 24 ) { sign= 1.0; }
671  else if( abs(sp_X.pdgid()) == 37 ) { sign=-1.0; }
672  else
673  {
674  cout<<"wrong sp_W/H.pdgid()="<<sp_X.pdgid()<<endl;
675  exit(-1);
676  }
677  if (sp_X.pdgid() > 0 )
678  {WTamplitM = WTamplit;} // tau- decay matrix element^2, spin averaged.
679  else
680  {WTamplitP = WTamplit;} // tau+ decay matrix element^2, spin averaged.
681 
682  // spin correlation weight. Tau production matrix element is represented by `sign'
683  double WT = 1.0+sign*HH[2]; // [2] means 'pz' component
684 
685  // Print out some info about the channel
686  DEBUG
687  (
688  cout<<tau_pdgid<<" -> ";
689  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
690  cout<<" (HH: "<<HH[0]<<" "<<HH[1]<<" "<<HH[2]<<" "<<HH[3]<<") WT: "<<WT<<endl;
691  )
692 
693  // TP:31Nov2013 checks of possible problems: weight outside theoretically allowed range
694  if (WT<0.0) {
695  printf("TauSpinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 0.0\n",WT);
696  WT = 0.0;
697  }
698 
699  if (WT>2.0) {
700  printf("Tauspinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 2.0\n",WT);
701  WT = 2.0;
702  }
703 
704  delete HH;
705 
706  return WT;
707 }
708 
709 /*******************************************************************************
710  Calculate weights, case of event record vertex like Z/gamma/H ... -> tau tau decay.
711 
712  Determine decay channel, calculates all necessary components for
713  calculation of all weights, calculates weights.
714 
715  Input: X four momentum may be larger than sum of tau1 tau2, missing component
716  is assumed to be QED brem
717 
718  Hidden input: relWTnonS, nonSM2 (used only in getLongitudinalPolarization(...) )
719  Hidden output: WTamplitM or WTamplitP of tau1 tau2 decays
720  weight for sp_X production matrix element^2 is calculated inside
721  plzap2
722  Polari - helicity attributed to taus, 100% correlations between tau+ and tau-
723  WTnonSM
724  Explicit output: WT spin correlation weight
725 *******************************************************************************/
726 double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector<SimpleParticle> &sp_tau1_daughters, vector<SimpleParticle> &sp_tau2_daughters)
727 {
728  // cout << "sp_tau1_daughters = " << sp_tau1_daughters.size() << endl;
729  // cout << "sp_tau2_daughters = " << sp_tau2_daughters.size() << endl;
730  SimpleParticle sp_tau;
731  SimpleParticle sp_nu_tau;
732  vector<SimpleParticle> sp_tau_daughters;
733 
734  // First we calculate HH for tau+
735  // We enforce that sp_tau is tau+ so the 'nu_tau' is tau-
736  if (sp_tau1.pdgid() == -15 )
737  {
738  sp_tau = sp_tau1;
739  sp_nu_tau = sp_tau2;
740  sp_tau_daughters = sp_tau1_daughters;
741  }
742  else
743  {
744  sp_tau = sp_tau2;
745  sp_nu_tau = sp_tau1;
746  sp_tau_daughters = sp_tau2_daughters;
747  }
748 
749  double *HHp, *HHm;
750 
751  // We use artificial if(true){... } construction to separate namespace for tau+ and tau-
752  if(true)
753  {
754  // Create Particles from SimpleParticles
755  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
756  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
757  Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
758 
759  vector<Particle> tau_daughters;
760 
761  // tau pdgid
762  int tau_pdgid = sp_tau.pdgid();
763 
764  // Create list of tau daughters
765  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
766  {
767  Particle pp(sp_tau_daughters[i].px(),
768  sp_tau_daughters[i].py(),
769  sp_tau_daughters[i].pz(),
770  sp_tau_daughters[i].e(),
771  sp_tau_daughters[i].pdgid() );
772 
773  tau_daughters.push_back(pp);
774  }
775 
776  double phi2 = 0.0, theta2 = 0.0;
777 
778  // To calculate matrix elements TAUOLA need tau decay products in tau rest-frame: we first boost all products
779  // to tau rest frame (tau other tau of Z/H decay along z-axis, intermediate step for boost is tau-tau pair rest-frame),
780  // then rotate to have neutrino from tau decay along z axis;
781  // calculated for that purpose angles phi2, theta2 are stored for rotation back of HHp
782  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
783 
784 
785 
786 
787  // Identify decay channel and then calculate polarimetric vector HH; calculates also WTamplit
788  HHp = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
789 
790  DEBUG
791  (
792  cout<<tau_pdgid<<" -> ";
793  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
794  cout<<" (HHp: "<<HHp[0]<<" "<<HHp[1]<<" "<<HHp[2]<<" "<<HHp[3]<<") ";
795  cout<<endl;
796  )
797 
798  WTamplitP = WTamplit;
799  } // end of if(true); for tau+
800 
801 
802 
803  // Second we calculate HH for tau-
804  // We enforce that sp_tau is tau- so the 'nu_tau' is tau+
805  if(sp_tau1.pdgid() == 15 )
806  {
807  sp_tau = sp_tau1;
808  sp_nu_tau = sp_tau2;
809  sp_tau_daughters = sp_tau1_daughters;
810  }
811  else
812  {
813  sp_tau = sp_tau2;
814  sp_nu_tau = sp_tau1;
815  sp_tau_daughters = sp_tau2_daughters;
816  }
817 
818  // We use artificial if(true){... } construction to separate namespace for tau+ and tau-
819  if(true)
820  {
821  // Create Particles from SimpleParticles
822  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
823  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
824  Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
825 
826  vector<Particle> tau_daughters;
827 
828  // tau pdgid
829  int tau_pdgid = sp_tau.pdgid();
830 
831  // Create list of tau daughters
832  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
833  {
834  Particle pp(sp_tau_daughters[i].px(),
835  sp_tau_daughters[i].py(),
836  sp_tau_daughters[i].pz(),
837  sp_tau_daughters[i].e(),
838  sp_tau_daughters[i].pdgid() );
839 
840  tau_daughters.push_back(pp);
841  }
842 
843  double phi2 = 0.0, theta2 = 0.0;
844 
845 
846  // To calculate matrix elements TAUOLA need tau decay products in tau rest-frame: we first boost all products
847  // to tau rest frame (tau other tau of Z/H decay along z-axis, intermediate step for boost is tau-tau pair rest-frame),
848  // then rotate to have neutrino from tau decay along z axis;
849  // calculated for that purpose angles phi2, theta2 are stored for rotation back of HHm
850  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
851 
852 
853  // Identify decay channel and then calculate polarimetric vector HHm; calculates also WTamplit
854  HHm = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
855 
856  DEBUG
857  (
858  cout<<tau_pdgid<<" -> ";
859  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
860  cout<<" (HHm: "<<HHm[0]<<" "<<HHm[1]<<" "<<HHm[2]<<" "<<HHm[3]<<") ";
861  cout<<endl;
862  )
863 
864  WTamplitM = WTamplit;
865  } // end of if(true); for tau-
866 
867  // CALCULATION OF PRODUCTION MATRIX ELEMENTS, THEN SPIN WEIGHTS AND FURTHER HIDDEN OUTPUTS
868 
869  // sign, in ultrarelativistic limit it is component of spin density matrix:
870  // longitudinal spin correlation for intermediate vector state gamma*,
871  // for Z etc. sign= +1; for Higgs, other scalars, pseudoscalars sign= -1
872  double sign = 1.0;
873  if(sp_X.pdgid() == 25) { sign=-1.0; }
874  if(sp_X.pdgid() == 36) { sign=-1.0; }
875  if(sp_X.pdgid() ==553) { sign=-1.0; } // upsilon(1s) can be treated as scalar
876 
877  double WT = 0.0;
878 
879  Polari = 0.0;
880  if(sign == -1.0) // Case of scalar
881  {
882  double S = sp_X.e()*sp_X.e() - sp_X.px()*sp_X.px() - sp_X.py()*sp_X.py() - sp_X.pz()*sp_X.pz();
883  IfHiggs=true; //global variable
884  double pol = getLongitudinalPolarization(S, sp_tau, sp_nu_tau); // makes sense only for nonSM otherwise NaN
885 
886  if(nonSM2==1) // WARNING it is for spin=2 resonance!!
887  {
888 
889  double corrX2;
890  double polX2;
891 
892  // NOTE: in this case, sp_nu_tau is the 2nd tau
893  // nonSMHcorrPol(S, sp_tau, sp_nu_tau, &corrX2, &polX2); // for future use
894  // WARNING: may be for polX2*HHm[2] we need to fix sign!
895  polX2=pol;
896  corrX2=-sign; // if X2 is of spin=2, spin correlation like for Z, we use RzXX,RzYY,RzXY,RzYX as transverse components of density matrix
897 
898  // possibly sign issue equivalent to overall sign in front of HHP (not cancelled sign of R with respect to L coupling) and pi-angle rotation around Y-axis
899  // WT = 1.0+corrX2*HHp[2]*HHm[2]+polX2*HHp[2]+polX2*HHm[2] + RzXX*HHp[0]*HHm[0] + RzYY*HHp[1]*HHm[1] + RzXY*HHp[0]*HHm[1] + RzYX*HHp[1]*HHm[0];
900  WT = 1.0+corrX2*HHp[2]*HHm[2]+polX2*HHp[2]+polX2*HHm[2] + RzXX*HHp[0]*HHm[0] - RzYY*HHp[1]*HHm[1] + RzXY*HHp[0]*HHm[1] - RzYX*HHp[1]*HHm[0];
901 
902  // we separate cross section into helicity parts. From this, we attribute helicity states to taus: ++ or --
903  double RRR = Tauola::randomDouble();
904  Polari=1.0;
905  if (RRR<(1.0+polX2)*(1.0+corrX2*HHp[2]*HHm[2]+HHp[2]+HHm[2])/(2.0+2.0*corrX2*HHp[2]*HHm[2]+2.0*polX2*HHp[2]+2.0*polX2*HHm[2])) Polari=-1.0;
906  }
907  else // case of Higgs
908  {
909  WT = 1.0 + sign*HHp[2]*HHm[2] + RXX*HHp[0]*HHm[0] + RYY*HHp[1]*HHm[1] + RXY*HHp[0]*HHm[1] + RYX*HHp[1]*HHm[0];
910 
911  // we separate cross section into helicity parts. From this, we attribute helicity states to taus: +- or -+
912  double RRR = Tauola::randomDouble();
913  Polari=1.0;
914  if (RRR<(1.0+sign*HHp[2]*HHm[2]+HHp[2]-HHm[2])/(2.0+2.0*sign*HHp[2]*HHm[2])) Polari=-1.0;
915  }
916  }
917  else // Case of Drell Yan
918  {
919 
920  double S = sp_X.e()*sp_X.e() - sp_X.px()*sp_X.px() - sp_X.py()*sp_X.py() - sp_X.pz()*sp_X.pz();
921 
922  // Get Z polarization
923  // ( Variable names are misleading! sp_tau is tau+ and sp_nu_tau is tau- )
924 
925  IfHiggs=false;
926 
927  double pol = getLongitudinalPolarization(S, sp_tau, sp_nu_tau);
928  // possibly sign issue equivalent to overall sign in front of HHP (not cancelled sign of R with respect to L coupling) and pi-angle rotation around Y-axis
929  // WT = 1.0+sign*HHp[2]*HHm[2]+pol*HHp[2]+pol*HHm[2] + RzXX*R11*HHp[0]*HHm[0] + RzYY*R22*HHp[1]*HHm[1] + RzXY*R12*HHp[0]*HHm[1] + RzYX*R21*HHp[1]*HHm[0];
930  WT = 1.0+sign*HHp[2]*HHm[2]+pol*HHp[2]+pol*HHm[2] + RzXX*R11*HHp[0]*HHm[0] - RzYY*R22*HHp[1]*HHm[1] + RzXY*R12*HHp[0]*HHm[1] - RzYX*R21*HHp[1]*HHm[0];
931  /*
932  // HHm[3] = 1.0; // initialized elsewhere
933  // HHp[3] = 1.0;
934  printf("--WT = %13.10f \n", WT);
935  printf("---------polarimetric vectors---------------------------\n");
936  printf(" HHm(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
937  HHm[0], HHm[1], HHm[2], HHm[3] );
938  printf(" HHp(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
939  HHp[0], HHp[1], HHp[2], HHp[3] );
940  printf("---------RcorOld (adjusted)--------------------\n");
941  printf("SM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
942  R11,R12,0.0,0.0);
943  printf("SM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
944  - R21,- R22,0.0,0.0);
945  printf("SM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
946  0.0,0.0,1.0,pol);
947  printf("SM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
948  0.0,0.0,pol,1.0);
949  printf("-------------------------------------------\n");
950  */
951  double zero=0.0;
952  double jeden=1.0;
953  // printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",R11,R12,zero,zero);
954  //printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",R21,R22,zero,zero);
955  //printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",zero,zero,sign,pol);
956  //printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",zero,zero,pol,jeden);
957  // in future we may need extra factor for wt which is
958  // F=PLWEIGHT(IDE,IDF,SVAR,COSTHE,1)
959  // but it is then non standard version of the code.
960 
961  // to correct when in the sample only spin corr. are in, but no polarization
962  if(Ipol==2) WT = WT/(1.0+sign*HHp[2]*HHm[2]);
963 
964  // to correct sample when corr. are in, pol. is in, but angular
965  // dependence of pol is missing.
966  if(Ipol==3)
967  {
968  // valid for configurations close to Z peak, otherwise approximation is at use.
969  double polp1 = plzap2(11,sp_tau.pdgid(),S,0.0);
970  double pol1 =(2*(1-polp1)-1) ;
971  WT = WT/(1.0+sign*HHp[2]*HHm[2]+pol1*HHp[2]+pol1*HHm[2]);
972  }
973  double WTtest;
974  if(ifkorch==1){
975  WTtest=0;
976  HHp[3]=1.0; // WARNING: temporary fix, til now HHp/m[3] were not used.
977  HHm[3]=1.0; // They were not always initialized by calculateHH()!
978 
979  for (int k0=0;k0<4;k0++){
980  for (int k1=0;k1<4;k1++){
981  double F=1.0;
982  // coefficient to switch off some elements of correlation matrix
983  if(k0==0&&k1==0) F=RzXX;
984  if(k0==1&&k1==1) F=RzYY;
985  if(k0==0&&k1==1) F=RzXY;
986  if(k0==1&&k1==0) F=RzYX;
987  WTtest=WTtest+RcorExt[k0][k1]*HHp[k0]*HHm[k1]*F;
988  }
989  }
990  /*
991  printf("--WT = %13.10f \n", WTtest);
992  printf("---------polarimetric vectors---------------------------\n");
993  printf(" HHm(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
994  HHm[0], HHm[1], HHm[2], HHm[3] );
995  printf(" HHp(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
996  HHp[0], HHp[1], HHp[2], HHp[3] );
997  printf("---------RcorExt (adjusted)---------------------------\n");
998  printf("SM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
999  RcorExt[0][0], RcorExt[0][1], RcorExt[0][2],RcorExt[0][3]);
1000  printf("SM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1001  RcorExt[1][0], RcorExt[1][1], RcorExt[1][2], RcorExt[1][3]);
1002  printf("SM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1003  RcorExt[2][0],RcorExt[2][1],RcorExt[2][2],RcorExt[2][3]);
1004  printf("SM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1005  RcorExt[3][0],RcorExt[3][1],RcorExt[3][2],RcorExt[3][3]);
1006  printf("-------------------------------------------\n");
1007  */
1008  R11 = RcorExt[0][0];
1009  R12 = RcorExt[0][1];
1010  R21 = RcorExt[1][0];
1011  R22 = RcorExt[1][1];
1012  R31 = RcorExt[2][0];
1013  R32 = RcorExt[2][1];
1014  R33 = RcorExt[2][2];
1015  R41 = RcorExt[3][0];
1016  R42 = RcorExt[3][1];
1017  R43 = RcorExt[3][2];
1018 
1019  // WTtest=WTtest; this line is to be activated/extened/modified at some tests.
1020  // cout << " why was minus sign in front of RzYY*R22*HHp[1]*HHm[1], WT (D-Y) formula?" << endl;
1021  // cout << " weight (non-korch) WT(nonSM="<<nonSM2<<")= " << WT << " spin-WTkorch= " << WTtest << " "<< endl;
1022  if(nonSM2==1) WT=WTtest; // to activate case of anomalous moment
1023  if(nonSM2==0&& Ipol==0) WT=WTtest; // to activate SM variant of full spin wt onto non polarized sample.
1024 
1025  }
1026  // we separate cross section into helicity parts. From this, we attribute helicity states to taus: ++ or --
1027  double RRR = Tauola::randomDouble();
1028  Polari=1.0;
1029  // ERW 6.04.25: added line for ifkorch
1030  // printf(" pol, RcorExt[3][2] = %13.10f %13.10f \n", pol, RcorExt[3][2]);
1031  if( ifkorch==1) pol = RcorExt[3][2];
1032  if (RRR<(1.0+pol)*(1.0+sign*HHp[2]*HHm[2]+HHp[2]+HHm[2])/(2.0+2.0*sign*HHp[2]*HHm[2]+2.0*pol*HHp[2]+2.0*pol*HHm[2])) Polari=-1.0;
1033  }
1034 
1035  // Print out some info about the channel
1036  DEBUG( cout<<" WT: "<<WT<<endl; )
1037 
1038  if (WT<0.0) {
1039  printf("Tauspinner::calculateWeightFromParticlesH WT is: %13.10f. Setting WT = 0.0\n",WT);
1040  WT = 0.0; // SwB:23Feb2013
1041  }
1042 
1043  if (WT>4.0) {
1044  printf("Tauspinner::calculateWeightFromParticlesH WT is: %13.10f. Setting WT = 4.0\n",WT);
1045  WT = 4.0; // SwB:23Feb2013
1046  }
1047 
1048  if( WT>4.0 || WT<0.0)
1049  {
1050  cout<<"Tauspinner::calculateWeightFromParticlesH ERROR: Z/gamma* or H, and WT not in range [0,4]."<<endl;
1051  exit(-1);
1052  }
1053 
1054  if (sign==-1.0 && nonSM2!=1) {
1055  if (WT>2.0) {
1056  WT = 2.0; // SwB:26Feb2013
1057  cout << "Tauspinner::calculateWeightFromParticlesH Setting WT to be 2.0" << endl;
1058  }
1059  }
1060 
1061  if( sign==-1.0 && (WT>2.0 || WT<0.0) && nonSM2!=1)
1062  {
1063  cout<<"Tauspinner::calculateWeightFromParticlesH ERROR: H and WT not in range [0,2]."<<endl;
1064  exit(-1);
1065  }
1066 
1067  delete[] HHp;
1068  delete[] HHm;
1069 
1070  return WT;
1071 }
1072 
1073 /*******************************************************************************
1074  Prepare kinematics for HH calculation
1075 
1076  Boost particles to effective bozon rest frame, and rotate them so that tau is on Z axis.
1077  Then rotate again with theta2 phi2 so neutrino from tau decay is along Z.
1078 *******************************************************************************/
1079 void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector<Particle> &tau_daughters, double *phi2, double *theta2)
1080 {
1081  Particle P_QQ( tau.px()+nu_tau.px(), tau.py()+nu_tau.py(), tau.pz()+nu_tau.pz(), tau.e()+nu_tau.e(), 0 );
1082 
1083  //cout<<endl<<"START: "<<endl;
1084  //print(P_QQ,nu_tau,tau,tau_daughters);
1085 
1086  // 1) boost tau, nu_tau and tau daughters to rest frame of P_QQ
1087 
1088  tau.boostToRestFrame(P_QQ);
1089  nu_tau.boostToRestFrame(P_QQ);
1090 
1091  for(unsigned int i=0; i<tau_daughters.size();i++)
1092  tau_daughters[i].boostToRestFrame(P_QQ);
1093 
1094  //cout<<endl<<"AFTER 1: "<<endl;
1095  //print(P_QQ,nu_tau,tau,tau_daughters);
1096 
1097  // 2) Rotate tau, nu_tau~, tau daughters to frame where tau is along Z
1098  // We set accompanying neutino in direction of Z+
1099 
1100  double phi = tau.getAnglePhi();
1101 
1102  tau.rotateXY(-phi);
1103 
1104  double theta = tau.getAngleTheta();
1105 
1106  tau.rotateXZ(M_PI-theta);
1107 
1108  nu_tau.rotateXY(-phi );
1109  nu_tau.rotateXZ(M_PI-theta);
1110 
1111  for(unsigned int i=0; i<tau_daughters.size();i++)
1112  {
1113  tau_daughters[i].rotateXY(-phi );
1114  tau_daughters[i].rotateXZ(M_PI-theta);
1115  }
1116 
1117  //cout<<endl<<"AFTER 2: "<<endl;
1118  //print(P_QQ,nu_tau,tau,tau_daughters);
1119 
1120  // 3) boost tau_daughters along Z to rest frame of tau
1121 
1122  for(unsigned int i=0; i<tau_daughters.size();i++)
1123  tau_daughters[i].boostAlongZ(-tau.pz(),tau.e());
1124 
1125  //cout<<endl<<"AFTER 3: "<<endl;
1126  //print(P_QQ,nu_tau,tau,tau_daughters);
1127 
1128  // 4) Now rotate tau daughters second time
1129  // so that nu_tau (from tau daughters list) is on Z axis
1130 
1131  // We can not be sure if tau_daughters[0] is neutrino !!!
1132  // That is the case, for tauola generated samples, but not in general!
1133  unsigned int i_stored = 0;
1134 
1135  *phi2=0;
1136  *theta2=0;
1137 
1138  for(unsigned int i=0; i<tau_daughters.size();i++)
1139  {
1140  if(abs(tau_daughters[i].pdgid())==16){
1141  *phi2 = tau_daughters[i].getAnglePhi();
1142 
1143  tau_daughters[i].rotateXY( -(*phi2) );
1144 
1145  *theta2 = tau_daughters[i].getAngleTheta();
1146 
1147  tau_daughters[i].rotateXZ( -(*theta2) );
1148 
1149  i_stored = i;
1150  break;
1151  }
1152  }
1153 
1154  for(unsigned int i=0; i<tau_daughters.size();i++)
1155  {
1156  if(i != i_stored) {
1157  tau_daughters[i].rotateXY( -(*phi2) );
1158  tau_daughters[i].rotateXZ( -(*theta2) );
1159  }
1160  }
1161 
1162  //cout<<endl<<"AFTER 4: "<<endl;
1163  //print(P_QQ,nu_tau,tau,tau_daughters);
1164 }
1165 
1166 /*******************************************************************************
1167  Calculates polarimetric vector HH of the tau decay.
1168 
1169  HH[3] is timelike because we use FORTRAN methods to calculate HH.
1170  First decide what is the channel. After that, 4-vectors
1171  are moved to tau rest frame of tau.
1172  Polarimetric vector HH is then rotated using angles phi and theta.
1173 
1174  Order of the tau decay products does not matter (it is adjusted)
1175 *******************************************************************************/
1176 double* calculateHH(int tau_pdgid, vector<Particle> &tau_daughters, double phi, double theta)
1177 {
1178  int channel = 0;
1179  double *HH = new double[4];
1180 
1181  HH[0]=HH[1]=HH[2]=HH[3]=0.0;
1182 
1183  vector<int> pdgid;
1184 
1185  // Create list of tau daughters pdgid
1186  for(unsigned int i=0; i<tau_daughters.size(); i++)
1187  pdgid.push_back( tau_daughters[i].pdgid() );
1188 
1189  // 17.04.2014: If Tauola++ is used for generation,
1190  // jaki_.ktom may be changed to 11 at the time of storing decay to event record
1191  // (case of full spin effects).
1192  // For Tauola++ jaki_.ktom is later of no use so 11 does not create problems.
1193  // For TauSpinner processing, jaki_.ktom should always be 1
1194  // This was the problem if Tauola++ generation and TauSpinner were used simultaneously.
1195  jaki_.ktom = 1;
1196 
1197  // tau^- --> pi^- nu_tau
1198  // tau^+ --> pi^+ anti_nu_tau
1199  // tau^- --> K^- nu_tau
1200  // tau^+ --> K^+ anti_nu_tau
1201  if( pdgid.size()==2 &&
1202  (
1203  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211) ) ||
1204  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211) ) ||
1205  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321) ) ||
1206  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321) )
1207  )
1208  ) {
1209  channel = 3; // channel: numbering convention of TAUOLA
1210  if(abs(pdgid[1])==321) channel = 6;
1211  DEBUG( cout<<"Channel "<<channel<<" : "; )
1212  // PXQ=AMTAU*EPI
1213  // PXN=AMTAU*ENU
1214  // QXN=PPI(4)*PNU(4)-PPI(1)*PNU(1)-PPI(2)*PNU(2)-PPI(3)*PNU(3)
1215 
1216  // BRAK=(GV**2+GA**2)*(2*PXQ*QXN-AMPI**2*PXN)
1217  // HV(I)=-ISGN*2*GA*GV*AMTAU*(2*PPI(I)*QXN-PNU(I)*AMPI**2)/BRAK
1218 
1219  const double AMTAU = 1.777;
1220  // is mass of the Pi+- OR K+-
1221  double AMPI = sqrt(tau_daughters[1].e() *tau_daughters[1].e()
1222  -tau_daughters[1].px()*tau_daughters[1].px()
1223  -tau_daughters[1].py()*tau_daughters[1].py()
1224  -tau_daughters[1].pz()*tau_daughters[1].pz());
1225 
1226  // two-body decay is so simple, that matrix element is calculated here
1227  double PXQ=AMTAU*tau_daughters[1].e();
1228  double PXN=AMTAU*tau_daughters[0].e();
1229  double QXN=tau_daughters[1].e()*tau_daughters[0].e()-tau_daughters[1].px()*tau_daughters[0].px()-tau_daughters[1].py()*tau_daughters[0].py()-tau_daughters[1].pz()*tau_daughters[0].pz();
1230  double BRAK=(2*PXQ*QXN-AMPI*AMPI*PXN);
1231 
1232  WTamplit = (1.16637E-5)*(1.16637E-5)*BRAK/2.;//AMPLIT=(GFERMI)**2*BRAK/2. //WARNING: Note for normalisation Cabbibo angle is missing!
1233  HH[0] = AMTAU*(2*tau_daughters[1].px()*QXN-tau_daughters[0].px()*AMPI*AMPI)/BRAK;
1234  HH[1] = AMTAU*(2*tau_daughters[1].py()*QXN-tau_daughters[0].py()*AMPI*AMPI)/BRAK;
1235  HH[2] = AMTAU*(2*tau_daughters[1].pz()*QXN-tau_daughters[0].pz()*AMPI*AMPI)/BRAK;
1236  HH[3] = 1.0;
1237  }
1238 
1239  // tau^- --> pi^- pi^0 nu_tau
1240  // tau^+ --> pi^+ pi^0 anti_nu_tau
1241  // tau^- --> K^- K^0 nu_tau; K^0 may be K_L K_S too.
1242  // tau^+ --> K^+ K^0 anti_nu_tau
1243  else if( pdgid.size()==3 &&
1244  (
1245  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 111) ) ||
1246  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 111) ) ||
1247  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 311) ) ||
1248  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 311) ) ||
1249  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 310) ) ||
1250  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 310) ) ||
1251  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 130) ) ||
1252  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 130) )
1253 
1254  )
1255  ) {
1256 
1257  channel = 4;
1258  DEBUG( cout<<"Channel "<<channel<<" : "; )
1259  // PRODPQ=PT(4)*QQ(4)
1260  // PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
1261  // PRODPN=PT(4)*PN(4)
1262  // BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
1263  // HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
1264 
1265  const double AMTAU = 1.777;
1266 
1267  int MNUM = 0;
1268  if(tau_daughters[2].pdgid() != 111) { MNUM=3; channel = 22;} // sub case of decay to K-K0
1269  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1270  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1271  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1272  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1273  float AMPLIT = 0.0;
1274  float HV[4] = { 0.0 };
1275 
1276  dam2pi_( &MNUM, PT, PN, PIM1, PIM2, &AMPLIT, HV );
1277 
1278  WTamplit = AMPLIT;
1279  HH[0] = -HV[0];
1280  HH[1] = -HV[1];
1281  HH[2] = -HV[2];
1282  HH[3] = HV[3];
1283  }
1284 
1285  // tau^- --> K^- pi^0 nu_tau
1286  // tau^+ --> K^+ pi^0 anti_nu_tau
1287  // tau^- --> pi^- K_S0 nu_tau
1288  // tau^+ --> pi^+ K_S0 anti_nu_tau
1289  // tau^- --> pi^- K_L0 nu_tau
1290  // tau^+ --> pi^+ K_L0 anti_nu_tau
1291  else if( pdgid.size()==3 &&
1292  (
1293  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 130) ) ||
1294  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 130) ) ||
1295  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 310) ) ||
1296  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 310) ) ||
1297  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 311) ) ||
1298  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 311) ) ||
1299  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 111) ) ||
1300  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 111) )
1301  )
1302  ) {
1303 
1304  channel = 7;
1305  DEBUG( cout<<"Channel "<<channel<<" : "; )
1306  // PRODPQ=PT(4)*QQ(4)
1307  // PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
1308  // PRODPN=PT(4)*PN(4)
1309  // BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
1310  // HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
1311 
1312  const double AMTAU = 1.777;
1313 
1314  double QQ[4];
1315  QQ[0]=tau_daughters[1].e() -tau_daughters[2].e() ;
1316  QQ[1]=tau_daughters[1].px()-tau_daughters[2].px();
1317  QQ[2]=tau_daughters[1].py()-tau_daughters[2].py();
1318  QQ[3]=tau_daughters[1].pz()-tau_daughters[2].pz();
1319 
1320  double PKS[4];
1321  PKS[0]=tau_daughters[1].e() +tau_daughters[2].e() ;
1322  PKS[1]=tau_daughters[1].px()+tau_daughters[2].px();
1323  PKS[2]=tau_daughters[1].py()+tau_daughters[2].py();
1324  PKS[3]=tau_daughters[1].pz()+tau_daughters[2].pz();
1325 
1326  // orthogonalization of QQ wr. to PKS
1327  double PKSD=PKS[0]*PKS[0]-PKS[1]*PKS[1]-PKS[2]*PKS[2]-PKS[3]*PKS[3];
1328  double QQPKS=QQ[0]*PKS[0]-QQ[1]*PKS[1]-QQ[2]*PKS[2]-QQ[3]*PKS[3];
1329 
1330  QQ[0]=QQ[0]-PKS[0]*QQPKS/PKSD;
1331  QQ[1]=QQ[1]-PKS[1]*QQPKS/PKSD;
1332  QQ[2]=QQ[2]-PKS[2]*QQPKS/PKSD;
1333  QQ[3]=QQ[3]-PKS[3]*QQPKS/PKSD;
1334 
1335  double PRODPQ=AMTAU*QQ[0];
1336  double PRODNQ=tau_daughters[0].e() *QQ[0]
1337  -tau_daughters[0].px()*QQ[1]
1338  -tau_daughters[0].py()*QQ[2]
1339  -tau_daughters[0].pz()*QQ[3];
1340  double PRODPN=AMTAU*tau_daughters[0].e();
1341  double QQ2 =QQ[0]*QQ[0]-QQ[1]*QQ[1]-QQ[2]*QQ[2]-QQ[3]*QQ[3];
1342 
1343  // in this case matrix element is calculated here
1344  double BRAK=(2*PRODPQ*PRODNQ-PRODPN*QQ2);
1345 
1346  WTamplit = (1.16637E-5)*(1.16637E-5)*BRAK/2.;//AMPLIT=(GFERMI)**2*BRAK/2. WARNING: Note for normalisation Cabibbo angle is missing!
1347  HH[0]=AMTAU*(2*PRODNQ*QQ[1]-QQ2*tau_daughters[0].px())/BRAK;
1348  HH[1]=AMTAU*(2*PRODNQ*QQ[2]-QQ2*tau_daughters[0].py())/BRAK;
1349  HH[2]=AMTAU*(2*PRODNQ*QQ[3]-QQ2*tau_daughters[0].pz())/BRAK;
1350  HH[3]=1.0;
1351  }
1352 
1353  // tau^- --> e^- anti_nu_e nu_tau
1354  // tau^+ --> e^+ nu_e anti_nu_tau
1355  else if( pdgid.size()==3 &&
1356  (
1357  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 11,-12) ) ||
1358  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-11, 12) )
1359  )
1360  ) {
1361  DEBUG( cout<<"Channel 1 : "; )
1362  channel = 1;
1363  // ITDKRC=0,XK0DEC=0.01 XK[4]={0},XA[4] nu_e, QP[4] e, XN[4] neutrino tauowe, AMPLIT, HH[4]
1364  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1365 
1366  int ITDKRC = 0;
1367  double XK0DEC = 0.01;
1368  double XK[4] = { 0.0 };
1369  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1370  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1371  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1372  double AMPLIT = 0.0;
1373  double HV[4] = { 0.0 };
1374 
1375  // We fix 4-momenta of electron and electron neutrino
1376  // Since electrons have small mass, they are prone to rounding errors
1377  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.511e-3*0.511e-3);
1378  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1379 
1380  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1381 
1382  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1383  HH[0] = -HV[0];
1384  HH[1] = -HV[1];
1385  HH[2] = -HV[2];
1386  HH[3] = HV[3];
1387  }
1388 
1389  // tau^- --> e^- anti_nu_e nu_tau + gamma
1390  // tau^+ --> e^+ nu_e anti_nu_tau + gamma
1391  else if( pdgid.size()==4 &&
1392  (
1393  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 11,-12, 22) ) ||
1394  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-11, 12, 22) )
1395  )
1396  ) {
1397  DEBUG( cout<<"Channel 1b : "; )
1398  channel = 1;
1399  // ITDKRC=0,XK0DEC=0.01 XK[4] gamma, XA[4] nu_e, QP[4] e, XN[4] neutrino tau , AMPLIT, HH[4]
1400  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1401 
1402  int ITDKRC = 1;
1403  double XK0DEC = 0.01;
1404  double XK[4] = { tau_daughters[3].px(), tau_daughters[3].py(), tau_daughters[3].pz(), tau_daughters[3].e() };
1405  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1406  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1407  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1408  double AMPLIT = 0.0;
1409  double HV[4] = { 0.0 };
1410 
1411  // We fix 4-momenta of electron and electron neutrino and photon
1412  // Since electrons have small mass, they are prone to rounding errors
1413  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.511e-3*0.511e-3);
1414  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1415  XK[3] = sqrt( XK[0]*XK[0] + XK[1]*XK[1] + XK[2]*XK[2] );
1416  // XK0DEC must be smaller in TauSpinner than what was used in generation. We do not use virt. corr anyway.
1417  if(XK0DEC > XK[3]/(XK[3]+XA[3]+QP[3]+XN[3])) XK0DEC=0.5*XK[3]/(XK[3]+XA[3]+QP[3]+XN[3]);
1418 
1419  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1420 
1421  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1422  HH[0] = -HV[0];
1423  HH[1] = -HV[1];
1424  HH[2] = -HV[2];
1425  HH[3] = HV[3];
1426  }
1427 
1428  // tau^- --> mu^- antui_nu_mu nu_tau
1429  // tau^+ --> mu^+ nu_mu anti_nu_tau
1430  else if( pdgid.size()==3 &&
1431  (
1432  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 13,-14) ) ||
1433  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-13, 14) )
1434  )
1435  ) {
1436 
1437  DEBUG( cout<<"Channel 2 : "; )
1438  channel = 2;
1439  // ITDKRC=0,XK0DEC=0.01 XK[4]={0},XA[4] nu_mu, QP[4] mu, XN[4] neutrino tauowe, AMPLIT, HH[4]
1440  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1441 
1442  int ITDKRC = 0;
1443  double XK0DEC = 0.01;
1444  double XK[4] = { 0.0 };
1445  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1446  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1447  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1448  double AMPLIT = 0.0;
1449  double HV[4] = { 0.0 };
1450 
1451  // We fix 4-momenta of muon and muon neutrino
1452  // Since muon have small mass, they are prone to rounding errors
1453  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.105659*0.105659);
1454  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1455 
1456  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1457 
1458  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1459  HH[0] = -HV[0];
1460  HH[1] = -HV[1];
1461  HH[2] = -HV[2];
1462  HH[3] = HV[3];
1463  }
1464 
1465  // tau^- --> mu^- antui_nu_mu nu_tau + gamma
1466  // tau^+ --> mu^+ nu_mu anti_nu_tau + gamma
1467  else if( pdgid.size()==4 &&
1468  (
1469  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 13,-14, 22) ) ||
1470  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-13, 14, 22) )
1471  )
1472  ) {
1473 
1474  DEBUG( cout<<"Channel 2b : "; )
1475  channel = 2;
1476  // ITDKRC=0,XK0DEC=0.01 XK[4] gamma, XA[4] nu_mu, QP[4] mu, XN[4] neutrino tau, AMPLIT, HH[4]
1477  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1478 
1479  int ITDKRC = 1;
1480  double XK0DEC = 0.01;
1481  double XK[4] = { tau_daughters[3].px(), tau_daughters[3].py(), tau_daughters[3].pz(), tau_daughters[3].e() };
1482  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1483  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1484  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1485  double AMPLIT = 0.0;
1486  double HV[4] = { 0.0 };
1487 
1488  // We fix 4-momenta of muon and muon neutrino and photon
1489  // Since muons have small mass, they are prone to rounding errors
1490  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.105659*0.105659);
1491  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1492  XK[3] = sqrt( XK[0]*XK[0] + XK[1]*XK[1] + XK[2]*XK[2] );
1493  // XK0DEC must be smaller in TauSpinner than what was used in generation. We do not use virt. corr anyway.
1494  if(XK0DEC > XK[3]/(XK[3]+XA[3]+QP[3]+XN[3])) XK0DEC=0.5*XK[3]/(XK[3]+XA[3]+QP[3]+XN[3]);
1495 
1496 
1497  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1498 
1499  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1500  HH[0] = -HV[0];
1501  HH[1] = -HV[1];
1502  HH[2] = -HV[2];
1503  HH[3] = HV[3];
1504  }
1505 
1506  // tau^- --> pi^- pi^0 pi^0 nu_tau
1507  // tau^+ --> pi^+ pi^0 pi^0 anti_nu_tau
1508  else if( pdgid.size()==4 &&
1509  (
1510  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 111, 111,-211) ) ||
1511  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 111, 111, 211) )
1512  )
1513  ) {
1514  DEBUG( cout<<"Channel 5 : "; )
1515  channel = 5;
1516  // MNUM=0, PT[4] tau, PN[4] neutrino, pi0[4], pi0[4], pi[4], AMPLIT, HH[4]
1517  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1518 
1519  const double AMTAU = 1.777;
1520  int MNUM = 0;
1521  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1522  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1523  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1524  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1525  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1526  float AMPLIT = 0.0;
1527  float HV[4] = { 0.0 };
1528 
1529  // For RChL currents one needs to define 3-pi sub-channel used
1530  chanopt_.JJ=2;
1531 
1532  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1533 
1534  WTamplit = AMPLIT;
1535  HH[0] = -HV[0];
1536  HH[1] = -HV[1];
1537  HH[2] = -HV[2];
1538  HH[3] = HV[3];
1539  }
1540 
1541  // tau^- --> pi^+ pi^- pi^- nu_tau
1542  // tau^+ --> pi^- pi^+ pi^+ anti_nu_tau
1543  else if( pdgid.size()==4 &&
1544  (
1545  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211,-211, 211) ) ||
1546  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 211,-211) )
1547  )
1548  ) {
1549  DEBUG( cout<<"Channel 5 : "; )
1550  channel = 5;
1551  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1552  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1553 
1554  const double AMTAU = 1.777;
1555  int MNUM = 0;
1556  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1557  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1558  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1559  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1560  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1561  float AMPLIT = 0.0;
1562  float HV[4] = { 0.0 };
1563 
1564  // For RChL currents one needs to define 3-pi sub-channel used
1565  chanopt_.JJ=1;
1566 
1567  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1568 
1569  WTamplit = AMPLIT;
1570  HH[0] = -HV[0];
1571  HH[1] = -HV[1];
1572  HH[2] = -HV[2];
1573  HH[3] = HV[3];
1574  }
1575 
1576  // tau^- --> K^+ pi^- pi^+ nu_tau // prepared for modes with kaons
1577  // tau^+ --> K^- pi^- pi^+ anti_nu_tau
1578 
1579  // tau^- --> pi^+ K^- K^- nu_tau // prepared for modes with kaons
1580  // tau^+ --> pi^- K^+ K^+ anti_nu_tau
1581  // tau^- --> K^+ K^- pi^- nu_tau // prepared for modes with kaons
1582  // tau^+ --> K^- K^+ pi^+ anti_nu_tau
1583 
1584  // tau^- --> pi^- K^0 pi^0 nu_tau // prepared for modes with kaons
1585  // tau^+ --> pi^+ K^0 pi^0 anti_nu_tau
1586 
1587 
1588  // tau^- --> pi^- K^0 K^0 nu_tau // prepared for modes with kaons
1589  // tau^+ --> pi^+ K^0 K^0 anti_nu_tau
1590 
1591 
1592  // tau^- --> K^- K^0 pi^0 nu_tau // prepared for modes with kaons
1593  // tau^+ --> K^+ K^0 pi^0 anti_nu_tau
1594 
1595  // tau^- --> K^- pi^0 pi^0 nu_tau // prepared for modes with kaons
1596  // tau^+ --> K^+ pi^0 pi^0 anti_nu_tau
1597 
1598  // 3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
1599  // 4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
1600  // 5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
1601 
1602  else if( pdgid.size()==4 &&
1603  (
1604  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, -211, 321) ) ||
1605  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 211,-321) )
1606  )
1607  ) {
1608  DEBUG( cout<<"Channel 5 : "; )
1609  channel = 14;
1610  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1611  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1612 
1613  const double AMTAU = 1.777;
1614 
1615  // IF(I.EQ.14) NAMES(I-7)=' TAU- --> K-, PI-, K+ '
1616  int MNUM = 1;
1617  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1618  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1619  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1620  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1621  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1622  float AMPLIT = 0.0;
1623  float HV[4] = { 0.0 };
1624 
1625  // For RChL currents one needs to define 3-pi sub-channel used
1626  chanopt_.JJ=1;
1627 
1628  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1629 
1630  WTamplit = AMPLIT;
1631  HH[0] = -HV[0];
1632  HH[1] = -HV[1];
1633  HH[2] = -HV[2];
1634  HH[3] = HV[3];
1635  }
1636 
1637 
1638 
1639  else if( pdgid.size()==4 &&
1640  (
1641  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 311, -211, 311 ) ) ||
1642  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 311, 211, 311 ) ) ||
1643  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 311, -211, 310 ) ) ||
1644  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 311, 211, 310 ) ) ||
1645  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 311, -211, 130 ) ) ||
1646  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 311, 211, 130 ) ) ||
1647  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 310, -211, 311 ) ) ||
1648  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 310, 211, 311 ) ) ||
1649  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 310, -211, 310 ) ) ||
1650  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 310, 211, 310 ) ) ||
1651  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 310, -211, 130 ) ) ||
1652  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 310, 211, 130 ) ) ||
1653  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 130, -211, 311 ) ) ||
1654  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 130, 211, 311 ) ) ||
1655  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 130, -211, 310 ) ) ||
1656  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 130, 211, 310 ) ) ||
1657  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 130, -211, 130 ) ) ||
1658  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 130, 211, 130 ) )
1659  )
1660  ) {
1661  DEBUG( cout<<"Channel 5 : "; )
1662  channel = 15;
1663  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1664  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1665 
1666  const double AMTAU = 1.777;
1667  // IF(I.EQ.15) NAMES(I-7)=' TAU- --> K0, PI-, K0B '
1668  int MNUM = 2;
1669  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1670  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1671  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1672  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1673  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1674  float AMPLIT = 0.0;
1675  float HV[4] = { 0.0 };
1676 
1677  // For RChL currents one needs to define 3-pi sub-channel used
1678  chanopt_.JJ=1;
1679 
1680  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1681 
1682  WTamplit = AMPLIT;
1683  HH[0] = -HV[0];
1684  HH[1] = -HV[1];
1685  HH[2] = -HV[2];
1686  HH[3] = HV[3];
1687  }
1688 
1689 
1690 
1691  else if( pdgid.size()==4 &&
1692  (
1693  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, 311, 111) ) ||
1694  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 311, 111) ) ||
1695  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, 310, 111) ) ||
1696  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 310, 111) ) ||
1697  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, 130, 111) ) ||
1698  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 130, 111) )
1699  )
1700  ) {
1701  DEBUG( cout<<"Channel 5 : "; )
1702  channel = 16;
1703  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1704  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1705 
1706  const double AMTAU = 1.777;
1707  // IF(I.EQ.16) NAMES(I-7)=' TAU- --> K-, K0, PI0 '
1708  int MNUM = 3;
1709  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1710  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1711  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1712  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1713  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1714  float AMPLIT = 0.0;
1715  float HV[4] = { 0.0 };
1716 
1717  // For RChL currents one needs to define 3-pi sub-channel used
1718  chanopt_.JJ=1;
1719 
1720  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1721 
1722  WTamplit = AMPLIT;
1723  HH[0] = -HV[0];
1724  HH[1] = -HV[1];
1725  HH[2] = -HV[2];
1726  HH[3] = HV[3];
1727  }
1728 
1729 
1730 
1731  else if( pdgid.size()==4 &&
1732  (
1733  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 111, 111,-321) ) ||
1734  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 111, 111, 321) )
1735  )
1736  ) {
1737  DEBUG( cout<<"Channel 5 : "; )
1738  channel = 17;
1739  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1740  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1741 
1742  const double AMTAU = 1.777;
1743  // IF(I.EQ.17) NAMES(I-7)=' TAU- --> PI0 PI0 K- ''
1744  int MNUM = 4;
1745  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1746  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1747  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1748  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1749  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1750  float AMPLIT = 0.0;
1751  float HV[4] = { 0.0 };
1752 
1753  // For RChL currents one needs to define 3-pi sub-channel used
1754  chanopt_.JJ=1;
1755 
1756  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1757 
1758  WTamplit = AMPLIT;
1759  HH[0] = -HV[0];
1760  HH[1] = -HV[1];
1761  HH[2] = -HV[2];
1762  HH[3] = HV[3];
1763  }
1764 
1765 
1766  else if( pdgid.size()==4 &&
1767  (
1768  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321,-211, 211) ) ||
1769  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 211,-211) )
1770  )
1771  ) {
1772  DEBUG( cout<<"Channel 5 : "; )
1773  channel = 18;
1774  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1775  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1776 
1777  const double AMTAU = 1.777;
1778  // IF(I.EQ.18) NAMES(I-7)=' TAU- --> K- PI- PI+ '
1779  int MNUM = 5;
1780  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1781  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1782  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1783  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1784  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1785  float AMPLIT = 0.0;
1786  float HV[4] = { 0.0 };
1787 
1788  // For RChL currents one needs to define 3-pi sub-channel used
1789  chanopt_.JJ=1;
1790 
1791  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1792 
1793  WTamplit = AMPLIT;
1794  HH[0] = -HV[0];
1795  HH[1] = -HV[1];
1796  HH[2] = -HV[2];
1797  HH[3] = HV[3];
1798  }
1799 
1800  else if( pdgid.size()==4 &&
1801  (
1802  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 311, 111) ) ||
1803  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 311, 111) ) ||
1804  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 310, 111) ) ||
1805  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 310, 111) ) ||
1806  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 130, 111) ) ||
1807  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 130, 111) )
1808  )
1809  ) {
1810  DEBUG( cout<<"Channel 5 : "; )
1811  channel = 19;
1812  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1813  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1814 
1815  const double AMTAU = 1.777;
1816  // IF(I.EQ.19) NAMES(I-7)=' TAU- --> PI- K0B PI0 '
1817  int MNUM = 6;
1818  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1819  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1820  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1821  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1822  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1823  float AMPLIT = 0.0;
1824  float HV[4] = { 0.0 };
1825 
1826  // For RChL currents one needs to define 3-pi sub-channel used
1827  chanopt_.JJ=1;
1828 
1829  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1830 
1831  WTamplit = AMPLIT;
1832  HH[0] = -HV[0];
1833  HH[1] = -HV[1];
1834  HH[2] = -HV[2];
1835  HH[3] = HV[3];
1836  }
1837  // tau^- --> pi^+ pi^+ pi^0 pi^- nu_tau
1838  // tau^+ --> pi^- pi^- pi^0 pi^+ anti_nu_tau
1839  else if( pdgid.size()==5 &&
1840  (
1841  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211,-211, 211, 111) ) ||
1842  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 211,-211, 111) )
1843  )
1844  ) {
1845  DEBUG( cout<<"Channel 8 : "; )
1846  channel = 8;
1847 
1848  const double AMTAU = 1.777;
1849  int MNUM = 1;
1850  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1851  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1852  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1853  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1854  float PIZ [4] = { (float)tau_daughters[4].px(), (float)tau_daughters[4].py(), (float)tau_daughters[4].pz(), (float)tau_daughters[4].e() };
1855  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1856  float AMPLIT = 0.0;
1857  float HV[4] = { 0.0 };
1858 
1859  dam4pi_( &MNUM, PT, PN, PIM1, PIM2, PIZ, PIPL, &AMPLIT, HV );
1860 
1861  WTamplit = AMPLIT;
1862  HH[0] = -HV[0];
1863  HH[1] = -HV[1];
1864  HH[2] = -HV[2];
1865  HH[3] = HV[3];
1866  }
1867  // tau^- --> pi^0 pi^0 pi^0 pi^- nu_tau
1868  // tau^+ --> pi^0 pi^0 pi^0 pi^+ anti_nu_tau
1869  else if( pdgid.size()==5 &&
1870  (
1871  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 111, 111, 111,-211) ) ||
1872  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 111, 111, 111, 211) )
1873  )
1874  ) {
1875  DEBUG( cout<<"Channel 9 : "; )
1876  channel = 9;
1877 
1878  const double AMTAU = 1.777;
1879  int MNUM = 2;
1880  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1881  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1882  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1883  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1884  float PIZ [4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1885  float PIPL[4] = { (float)tau_daughters[4].px(), (float)tau_daughters[4].py(), (float)tau_daughters[4].pz(), (float)tau_daughters[4].e() };
1886  float AMPLIT = 0.0;
1887  float HV[4] = { 0.0 };
1888 
1889  dam4pi_( &MNUM, PT, PN, PIM1, PIM2, PIZ, PIPL, &AMPLIT, HV );
1890 
1891  WTamplit = AMPLIT;
1892  HH[0] = -HV[0];
1893  HH[1] = -HV[1];
1894  HH[2] = -HV[2];
1895  HH[3] = HV[3];
1896  }
1897  else {
1898 
1899  DEBUG( cout<<tau_daughters.size()<<"-part ???: "; )
1900 
1901  }
1902 
1903  // Now rotate vector HH using angles phi and theta
1904  Particle HHbuf(HH[0], HH[1], HH[2], HH[3], 0);
1905 
1906  HHbuf.rotateXZ(theta);
1907  HHbuf.rotateXY(phi);
1908 
1909  HH[0] = HHbuf.px();
1910  HH[1] = HHbuf.py();
1911  HH[2] = HHbuf.pz();
1912  HH[3] = HHbuf.e ();
1913 
1914  return HH;
1915 }
1916 
1917 /*******************************************************************************
1918  Returns longitudinal polarization of the single tau (in Z/gamma* -> tau+ tau- case) averaged over
1919  incoming configurations
1920  S: invariant mass^2 of the bozon
1921  &sp_tau: first tau
1922  &sp_nu_tau: second tau (in this case it is misleading name)
1923  Hidden output: WTnonSM
1924  Hidden input: relWTnonS, nonSM2
1925 *******************************************************************************/
1926 double getLongitudinalPolarization(double S, SimpleParticle &sp_tau, SimpleParticle &sp_nu_tau)
1927 {
1928 
1929  // std::cout << "nonSM2 = " << nonSM2 << std::endl;
1930 
1931  double ReA0=0; // 0 is added to A-coupling name, because SM part may be added internally
1932  double ImA0=0;
1933  double ReB=0;
1934  double ImB=0;
1935  double ReX=0;
1936  double ImX=0;
1937  double ReY=0;
1938  double ImY=0;
1939  double ReA00=0; // 0 is added to A-coupling name, because SM part may be added internally
1940  double ImA00=0;
1941  double ReB0=0;
1942  double ImB0=0;
1943  double ReX0=0;
1944  double ImX0=0;
1945  double ReY0=0;
1946  double ImY0=0;
1947 
1948  double GSWr[7],GSWi[7];
1949  double GSWr0[7],GSWi0[7];
1950  double Amz0;
1951  double Gamz0;
1952  double sin2W0;
1953  double alphaQED;
1954  double sumak=0.0;
1955  double xsec,xsec0,xsecsum;
1956  double GAMfrac=0.0;
1957  double GAMfrac2=0.0; // ratio of gamma gamma process to the rest, in generated sample and in rewieght
1958  double A0=0.0; // anomalous couplings for gamma gamma in sample
1959  double B0=0.0;
1960  double A=0.0; // in reweight
1961  double B=0.0;
1962 
1963  if(ifGSW==0){
1964  // ExtraEWparamsGet(&Amz0,&Gamz0,&sin2W0,&alfinv,&DeltSQ,&DeltV,&Gmu,&keyGSW);
1965  // these parameters should be as of effective Born, or whatever is used in sample preparation
1966  Amz0=91.1876;
1967  Gamz0=2.4952;// 2.5;
1968  sin2W0= 0.2315200; //0.231499;
1969  alphaQED= 1.0/128.86674175;//1.0/128.9503022;
1970 
1971  for (int k=0;k<7;k++){
1972  GSWr[k]=1.0; GSWr0[k]=1.0;
1973  GSWi[k]=0.0; GSWi0[k]=0.0;
1974  }
1975 
1976  }
1977 
1978  if(ifkorch==1 && nonSM2==0) {
1979 
1980  ReA0=0; // 0 is added to A-coupling name, because SM part may be added internally
1981  ImA0=0;
1982  ReB=0;
1983  ImB=0;
1984  ReX=0;
1985  ImX=0;
1986  ReY=0;
1987  ImY=0;
1988  GAMfrac=GAMfraci; // ratio of gamma gamma process to the rest, in generated sample
1989  GAMfrac2=GAMfrac2i; // and in rewieght
1990  A0=0.; //A0i; // anomalous couplings for gamma gamma in sample
1991  B0=0.; //B0i;
1992  A=0.; //Ai; // in reweight
1993  B=0.; //Bi;
1994 
1995 
1996  }
1997 
1998  if(ifkorch==1 && nonSM2==1) {
1999 
2000  ReA0=ReAini;
2001  ImA0=ImAini;
2002  ReB=ReBini;
2003  ImB=ImBini;
2004  ReX=ReXini;
2005  ImX=ImXini;
2006  ReY=ReYini;
2007  ImY=ImYini;
2008  GAMfrac=GAMfraci; // ratio of gamma gamma process to the rest, in generated sample
2009  GAMfrac2=GAMfrac2i; // and in rewieght
2010  A0=A0i; // anomalous couplings for gamma gamma in sample
2011  B0=B0i;
2012  A=Ai; // in reweight
2013  B=Bi;
2014 
2015  }
2016 
2017  // tau+ and tau- in lab frame
2018  Particle tau_plus ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
2019  Particle tau_minus( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
2020  // P_QQ = sum of tau+ and tau- in lab frame
2021  Particle P_QQ( tau_plus.px()+tau_minus.px(), tau_plus.py()+tau_minus.py(), tau_plus.pz()+tau_minus.pz(), tau_plus.e()+tau_minus.e(), 0 );
2022 
2023  Particle P_B1(0, 0, 1, 1, 0);
2024  Particle P_B2(0, 0,-1, 1, 0);
2025 
2026  tau_plus. boostToRestFrame(P_QQ);
2027  tau_minus.boostToRestFrame(P_QQ);
2028  P_B1. boostToRestFrame(P_QQ);
2029  P_B2. boostToRestFrame(P_QQ);
2030 
2031  double costheta1 = (tau_plus.px()*P_B1.px() +tau_plus.py()*P_B1.py() +tau_plus.pz()*P_B1.pz() ) /
2032  sqrt(tau_plus.px()*tau_plus.px()+tau_plus.py()*tau_plus.py()+tau_plus.pz()*tau_plus.pz()) /
2033  sqrt(P_B1.px() *P_B1.px() +P_B1.py() *P_B1.py() +P_B1.pz() *P_B1.pz() );
2034 
2035  double costheta2 = (tau_minus.px()*P_B2.px() +tau_minus.py()*P_B2.py() +tau_minus.pz()*P_B2.pz() ) /
2036  sqrt(tau_minus.px()*tau_minus.px()+tau_minus.py()*tau_minus.py()+tau_minus.pz()*tau_minus.pz()) /
2037  sqrt(P_B2.px() *P_B2.px() +P_B2.py() *P_B2.py() +P_B2.pz() *P_B2.pz() );
2038 
2039  double sintheta1 = sqrt(1-costheta1*costheta1);
2040  double sintheta2 = sqrt(1-costheta2*costheta2);
2041 
2042  // Cosine of hard scattering
2043  double costhe = (costheta1*sintheta2 + costheta2*sintheta1) / (sintheta1 + sintheta2);
2044  double costhem= (costheta1*sintheta2 + costheta2*sintheta1) / (sintheta1 + sintheta2);
2045 
2046  // Cosine of hard scattering Mustraal frame style draft of the code.
2047  if (FrameType==1){
2048  double RRR = Tauola::randomDouble();
2049  double wtMustraal;
2050  wtMustraal= P_B1.e()*P_B1.e()*(1+costheta1*costheta1)/(P_B1.e()*P_B1.e()*(1+costheta1*costheta1)+P_B2.e()*P_B2.e()*(1+costheta2*costheta2));
2051  if( RRR< wtMustraal) {costhe=costheta1; costhem=costheta2;} // we will need second choice for interchanged beams
2052  else {costhe=costheta2; costhem=costheta1;}
2053 
2054  // wherever below "-costhe" or "-cost" is used "-costhem" should be used instead.
2055  // calls with costhe and -costhem should be grouped.
2056 }
2057  // Invariant mass^2 of, tau+tau- pair plus photons, system!
2058  double SS = S; // other option is for tests: P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
2059 
2060  /*
2061  We need to fix sign of costhe and attribute ID.
2062  calculate x1,x2; // x1*x2 = SS/CMSENE/CMSENE; // (x1-x2)/(x1+x2)=P_QQ/CMSENE*2 in lab;
2063  calculate weight WID[]=sig(ID,SS,+/-costhe)* f(x1,+/-ID,SS) * f(x2,-/+ID,SS) ; respectively for u d c s b
2064  f(x,ID,SS,CMSENE)=x*(1-x) // for the start it will invoke library
2065  on the basis of this generate ID and set sign for costhe.
2066  then we calculate polarization averaging over incoming states.
2067  */
2068 
2069  double x1x2 = SS/CMSENE/CMSENE;
2070  double x1Mx2 = P_QQ.pz()/CMSENE*2;
2071 
2072  double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
2073  double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
2074 
2075  double WID[13];
2076  WID[0] = f(x1, 0,SS,CMSENE)*f(x2, 0,SS,CMSENE) * sigborn(0,SS, costhe);
2077  WID[1] = f(x1, 1,SS,CMSENE)*f(x2,-1,SS,CMSENE) * sigborn(1,SS, costhe);
2078  WID[2] = f(x1,-1,SS,CMSENE)*f(x2, 1,SS,CMSENE) * sigborn(1,SS,-costhem);
2079  WID[3] = f(x1, 2,SS,CMSENE)*f(x2,-2,SS,CMSENE) * sigborn(2,SS, costhe);
2080  WID[4] = f(x1,-2,SS,CMSENE)*f(x2, 2,SS,CMSENE) * sigborn(2,SS,-costhem);
2081  WID[5] = f(x1, 3,SS,CMSENE)*f(x2,-3,SS,CMSENE) * sigborn(3,SS, costhe);
2082  WID[6] = f(x1,-3,SS,CMSENE)*f(x2, 3,SS,CMSENE) * sigborn(3,SS,-costhem);
2083  WID[7] = f(x1, 4,SS,CMSENE)*f(x2,-4,SS,CMSENE) * sigborn(4,SS, costhe);
2084  WID[8] = f(x1,-4,SS,CMSENE)*f(x2, 4,SS,CMSENE) * sigborn(4,SS,-costhem);
2085  WID[9] = f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(5,SS, costhe);
2086  WID[10]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(5,SS,-costhem);
2087  WID[11]= GAMfrac*(WID[0]+WID[1]+WID[2]+WID[3]+WID[4]+WID[5]+WID[6]+WID[7]+WID[8]+WID[9]+WID[10]);
2088  WID[12]= 0.0;
2089 
2090  //ERW: Genis
2091  // WID[11]= f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(22,SS, costhe);
2092  // WID[12]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(22,SS,-costhe);
2093 
2094  // REMOVE !!!
2095  //for(int i=0;i<=10;i++) {if(!(i==3)) WID[i]=0.0;}
2096  // printf(" usunac widy nadpisane [3]= %13.10f \n",WID[3] );
2097  // printf(" usunac widy nadpisane [1]= %13.10f \n",WID[1] );
2098 
2099  // end REMOVE !!!
2100 
2101  double sum = 0.0; // normalize
2102  for(int i=0;i<=12;i++) sum+=WID[i];
2103 
2104  //ERW:: Genis
2105  // for(int i=0;i<=12;i++) WID[i] = 0;
2106  //WID[1] = sigborn(11, SS, costhe);
2107  //sum = WID[1];
2108 
2109 
2110 
2111  if( sum == 0.0 )
2112  {
2113  cout << "Tauspinner::calculateWeightFromParticlesH WARNING: sum of WID[0]-WID[10] is 0. Check LHAPDF configuration" << endl;
2114  }
2115 
2116  WTnonSM=1.0;
2117  if(relWTnonSM==0) WTnonSM=sum;
2118  if(nonSM2==1)
2119  {
2120  double WID2[13];
2121  WID2[0] = f(x1, 0,SS,CMSENE)*f(x2, 0,SS,CMSENE) * sigborn(0,SS, costhe) * plweight(0,SS, costhe); // plweight = ratio of cross-section nonSM/SM
2122  WID2[1] = f(x1, 1,SS,CMSENE)*f(x2,-1,SS,CMSENE) * sigborn(1,SS, costhe) * plweight(1,SS, costhe);
2123  WID2[2] = f(x1,-1,SS,CMSENE)*f(x2, 1,SS,CMSENE) * sigborn(1,SS,-costhem) * plweight(1,SS,-costhem);
2124  WID2[3] = f(x1, 2,SS,CMSENE)*f(x2,-2,SS,CMSENE) * sigborn(2,SS, costhe) * plweight(2,SS, costhe);
2125  WID2[4] = f(x1,-2,SS,CMSENE)*f(x2, 2,SS,CMSENE) * sigborn(2,SS,-costhem) * plweight(2,SS,-costhem);
2126  WID2[5] = f(x1, 3,SS,CMSENE)*f(x2,-3,SS,CMSENE) * sigborn(3,SS, costhe) * plweight(3,SS, costhe);
2127  WID2[6] = f(x1,-3,SS,CMSENE)*f(x2, 3,SS,CMSENE) * sigborn(3,SS,-costhem) * plweight(3,SS,-costhem);
2128  WID2[7] = f(x1, 4,SS,CMSENE)*f(x2,-4,SS,CMSENE) * sigborn(4,SS, costhe) * plweight(4,SS, costhe);
2129  WID2[8] = f(x1,-4,SS,CMSENE)*f(x2, 4,SS,CMSENE) * sigborn(4,SS,-costhem) * plweight(4,SS,-costhem);
2130  WID2[9] = f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(5,SS, costhe) * plweight(5,SS, costhe);
2131  WID2[10]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(5,SS,-costhem) * plweight(5,SS,-costhem);
2132  // WID2[11]= f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(22,SS, costhe)* plweight(22,SS, costhe);
2133  // WID2[12]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(22,SS,-costhem)* plweight(22,SS,-costhem);
2134  WID2[11]= GAMfrac2*(WID2[0]+WID2[1]+WID2[2]+WID2[3]+WID2[4]+WID2[5]+WID2[6]+WID2[7]+WID2[8]+WID2[9]+WID2[10]);
2135  WID2[12]= 0.0;
2136 
2137  // REMOVE !!!
2138  // for(int i=0;i<=10;i++){ if(!(i==3)) WID2[i]=0.0;}
2139  // printf(" usunac widy2 nadpisane [2]= %13.10f \n",WID2[3] );
2140  // printf(" usunac widy2 nadpisane [1]= %13.10f \n",WID2[1] );
2141  // end REMOVE !!!
2142 
2143  double sum2 = 0.0; // normalize
2144  for(int i=0;i<=12;i++) sum2+=WID2[i];
2145  WTnonSM=sum2/sum ;
2146  // WTnonSM= WTnonSM*1.2;
2147  // printf(" WTnonSM= %13.10f \n", WTnonSM);
2148  if(relWTnonSM==0) WTnonSM=sum2;
2149 
2150  }
2151 
2152  double pol = 0.0;
2153  // double Rxx = 0.0;
2154  // double Ryy = 0.0;
2155 
2156  if(IfHiggs && nonSM2==1) { // we assume that only glue glue process contributes for Higgs
2157  double polp = plzap2(0,15,S,costhe); // 0 means incoming gluon, 15 means outgoing tau
2158  pol += (2*(1-polp)-1);
2159  return pol;
2160  }
2161  if(IfHiggs) return NAN;
2162 
2163  // case of Z/gamma
2164  for(int i=0;i<=12;i++) WID[i]/=sum;
2165 
2166 
2167  R11 = 0.0;
2168  R22 = 0.0;
2169  R12 = 0.0;
2170  R21 = 0.0;
2171 
2172  int ICC = -1;
2173 
2174  sumak=0.0; // not used at present
2175  xsecsum=0.0;
2176  if(ifkorch==1){
2177  for(int i=0;i<=3;i++){
2178  for(int j=0;j<=3;j++){
2179  RcorExt[i][j]=0.0;
2180  }
2181  }
2182  }
2183 
2184  for(int i=1;i<=12;i++)
2185  {
2186  //printf("##### petla po i = %3i \n",i );
2187  ICC = i;
2188  double cost = costhe;
2189  // first beam quark or antiquark
2190  int IROT=0;
2191  if( ICC==2 || ICC==4 || ICC==6 || ICC==8 || ICC==10|| ICC==12 ) { cost = -costhem; IROT=1;}
2192 
2193  // ID of incoming quark (up or down type)
2194  int ID = 2;
2195  if( ICC==7 || ICC==8 ) ID = 4;
2196  if( ICC==1 || ICC==2 ) ID = 1;
2197  if( ICC==5 || ICC==6 ) ID = 3;
2198  if( ICC==9 || ICC==10 ) ID = 5;
2199  if( ICC==11 || ICC==12 ) ID = 22;
2200  int tau_pdgid = 15;
2201  if(i<11){
2202  // printf("##### kic \n");
2203  double polp = plzap2(ID,tau_pdgid,S,cost);
2204  // if(ifkorch==0) pol += (2*(1-polp)-1)*WID[i];
2205  pol += (2*(1-polp)-1)*WID[i];
2206  // printf("##### klic \n");
2207  // we obtain transverse spin components of density matrix from Tauolapp pre-tabulated O(alpha) EW results
2208  //
2209  Tauolapp::TauolaParticlePair pp;
2210  // printf("##### akikur \n");
2211  // Set them to 0 in case no tables are loaded by Tauola++
2212  pp.m_R[1][1] = pp.m_R[2][2] = 0.0;
2213  pp.m_R[1][2] = pp.m_R[2][1] = 0.0;
2214  pp.m_R[0][1] = pp.m_R[3][2] = 0.0;
2215  pp.m_R[0][2] = pp.m_R[3][1] = 0.0;
2216  pp.m_R[1][0] = pp.m_R[2][3] = 0.0;
2217  pp.m_R[2][0] = pp.m_R[1][3] = 0.0;
2218  pp.recalculateRij(ID,tau_pdgid,S,cost);
2219 
2220  // S=(91.17+10.)*(91.17+10.); //100.0;
2221  // costhe=0.5;
2222  //ERW:: Genis
2223  //pp.recalculateRij(11,tau_pdgid,S,costhe);
2224  //pol = pp.m_R[0][3]/pp.m_R[0][0];
2225 
2226  // ERW:: Genis for electrons
2227  double Rcor[4][4];
2228  // double RcorExt[4][4];
2229  double Energy;
2230  double theta;
2231 
2232 
2233  Energy=sqrt(S)/2.0;
2234  theta=acos(cost);
2235 
2236 
2237  if(ifkorch==1){
2238  // ========================= BEGIN ============
2239 
2240  int IFLAV; // 0 for LEPTONS, 2 -for UP QUARK, 1 - for DOWN QUARK
2241  int channel ; // 1 for LEPTONS, 2 -for UP QUARK, 3 - for DOWN QUARK
2242  channel = 1;
2243  IFLAV=0;
2244  if (ID==2||ID==4) {channel=2 ; IFLAV=2;}
2245  if (ID==1||ID==3||ID==5){channel=3 ; IFLAV=1;}
2246 
2247 
2248 
2249  //std::cout <<
2250  // std::cout << "Amz = " << Amz0 << std::endl;
2251  // std::cout << "Gamz = " << Gamz0 << std::endl;
2252  // std::cout << "sin2W = " << sin2W0 << std::endl;
2253 
2254  if(ifGSW==1){
2255  Amz0=Amz(IFLAV);
2256  Gamz0=Gamz(IFLAV);
2257  sin2W0=sin2W(IFLAV);
2258  alphaQED=7.2973525693e-3; // check the value, should be as for alpha(s=0)= 1/137.035...
2259  for (int k=0;k<7;k++){
2260  complex<double> rezu = EWFACT(IFLAV, k, S, cost);
2261 
2262 
2263  GSWr[k]=real(rezu);GSWr0[k]=1.0;//real(rezu);
2264  GSWi[k]=imag(rezu);GSWi0[k]=0.0;//imag(rezu);
2265  // std::cout << "EWFACT ( "<< k <<")= " << GSWr[k]<<" "<< GSWi[k] << std::endl;
2266  }
2267  }
2268  GSWadapt(IFLAV,S,cost,GSWr,GSWi);
2269  dipolqq_(&iqed,&Energy,&theta,&channel,&Amz0,&Gamz0,&sin2W0,&alphaQED,&ReA00,&ImA00, &ReB0,&ImB0, &ReX0,&ImX0,&ReY0,&ImY0, GSWr0,GSWi0,Rcor);
2270  xsec0=Rcor[0][0];
2271 
2272  dipolqq_(&iqed,&Energy,&theta,&channel,&Amz0,&Gamz0,&sin2W0,&alphaQED,&ReA0,&ImA0, &ReB,&ImB, &ReX,&ImX,&ReY,&ImY, GSWr,GSWi,Rcor);
2273  xsec=Rcor[0][0];
2274  Rcor[0][0]=1.0;
2275  // TauSpinner and Korchin convention change is performed in routine dipolqq_ which provides matrix R_ij in C++ conventions
2276  // 0123 <-> txyz
2277 
2278  if(IROT==1){// needed for cost=-costhem: rotate by pi around z-axis.
2279  /*
2280  //Rotate pp.m_R around z-axis as well? May be in future.
2281  //Not used at present not necessary adjustment of frame orientation for elements which equal zero.
2282  pp.m_R[1][3]=-pp.m_R[1][3];
2283  pp.m_R[3][1]=-pp.m_R[3][1];
2284  pp.m_R[0][1]=-pp.m_R[0][1];
2285  pp.m_R[1][0]=-pp.m_R[1][0];
2286  pp.m_R[2][3]=-pp.m_R[2][3];
2287  pp.m_R[3][2]=-pp.m_R[3][2];
2288  pp.m_R[2][0]=-pp.m_R[2][0];
2289  pp.m_R[0][2]=-pp.m_R[0][2];
2290  */
2291  Rcor[1][3]=-Rcor[1][3];
2292  Rcor[3][1]=-Rcor[3][1];
2293  Rcor[0][1]=-Rcor[0][1];
2294  Rcor[1][0]=-Rcor[1][0];
2295  Rcor[2][3]=-Rcor[2][3];
2296  Rcor[3][2]=-Rcor[3][2];
2297  Rcor[2][0]=-Rcor[2][0];
2298  Rcor[0][2]=-Rcor[0][2];
2299 
2300 // note that till Jan 3 2024 these components of R were not used in any applications.
2301  }
2302  // pol += WID[i]*pp.m_R[0][3]/pp.m_R[0][0];
2303  // std::cout << "przed daje pol("<<i<<") = " << pol << std::endl;
2304  // INFO:
2305  // pol += WID[i]*Rcor[2][3]; not needed there is:
2306  // if( ifkorch==1) pol = RcorExt[3][2];
2307  // in calculateWeightFromParticlesH(){...
2308  //
2309  // pp.m_R[0][3]/pp.m_R[0][0];
2310  // printf(" i= %3i Wid= %13.10f pol(i)= %13.10f pol= %13.10f\n",i, WID[i],pp.m_R[0][3]/pp.m_R[0][0],pol);
2311  // sumak +=WID[i]*pp.m_R[0][0]; --> replaced by xsecsum
2312  }
2313 
2314  // std::cout << "pol("<<i<<") = " << pol << std::endl;
2315  //std::cout << " m_R = " << pp.m_R[0][0] << " " << pp.m_R[0][1] << " " << pp.m_R[0][2] << " " << pp.m_R[0][3] << std::endl;
2316  //std::cout << " m_R = " << pp.m_R[1][0] << " " << pp.m_R[1][1] << " " << pp.m_R[1][2] << " " << pp.m_R[1][3] << std::endl;
2317  //std::cout << " m_R = " << pp.m_R[2][0] << " " << pp.m_R[2][1] << " " << pp.m_R[2][2] << " " << pp.m_R[2][3] << std::endl;
2318  //std::cout << " m_R = " << pp.m_R[3][0] << " " << pp.m_R[3][1] << " " << pp.m_R[3][2] << " " << pp.m_R[3][3] << std::endl;
2319 
2320  // ========================= END (korch==1) ============
2321  // printf("############# pp.m_R[1][1]= %13.10f \n", pp.m_R[1][1]);
2322  // These calculations are valid for Standard Model only
2323 
2324  //ERW 6.04.2025: pp.m_R jest w konwencji pochodzacej z fortranu.
2325  // to wypelnianie tylko dla (ifkorch==0)??
2326  R11 += WID[i]*pp.m_R[1][1];
2327  R22 += WID[i]*pp.m_R[2][2]; // note signs in line WT =
2328  R12 += WID[i]*pp.m_R[1][2]; // 1.0+sign*HHp[2]*HHm[2] ....- RzYY*R22 ....- RzYX*R21*HHp[1]*HHm[0];
2329  R21 += WID[i]*pp.m_R[2][1]; // sign operation with F=1/-1 like for RcorExt is not necessary.
2330 
2331  if(ifkorch==1){
2332  for(int i3=0;i3<=3;i3++){
2333  for(int j=0;j<=3;j++){
2334  int i4=i3-1;
2335  int j4=j-1;
2336  if(i4==-1) i4=3;
2337  if(j4==-1) j4=3;
2338  RcorExt[i4][j4]+=WID[i]*Rcor[i3][j]; //we move back to fortran convention xyzt because of f77-tauola polarimetric vectors
2339  }
2340  }
2341 
2342  xsecsum+=xsec/xsec0*WID[i];
2343 
2344  }
2345  }
2346 
2347  else{
2348  //ERW 6.04.2025: petla else tylko dla (ifkorch == 1)
2349  double Rcor[4][4];
2350  double polp = 0.5; // gamagama bring no polarization
2351  pol += (2*(1-polp)-1)*WID[i];
2352  // iqed QED part of anomalous switched on/off.
2353 
2354  double E;
2355  double theta;
2356  E=sqrt(S)/2.0;
2357  theta=acos(cost);
2358  dipolgamma_(&iqed, &E, &theta, &A0, &B0, Rcor);
2359  /*
2360  printf(" -------------------------------------------------\n");
2361  //ERW: rotation should be moved to dipolgamma OR NOT??
2362  printf("SM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2363  Rcor[0][0],Rcor[0][1],Rcor[0][2],Rcor[0][3]);
2364  printf("SM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2365  Rcor[1][0],Rcor[1][1],Rcor[1][2],Rcor[1][3]);
2366  printf("SM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2367  Rcor[2][0],Rcor[2][1],Rcor[2][2],Rcor[2][3]);
2368  printf("SM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2369  Rcor[3][0],Rcor[3][1],Rcor[3][2],Rcor[3][3]);
2370  */
2371  xsec0=Rcor[0][0];
2372 
2373 
2374  dipolgamma_(&iqed, &E, &theta, &A, &B, Rcor);
2375  /*
2376  printf(" -----------------------------------------------\n");
2377  //ERW: rotation should be moved to dipolgamma OR NOT??
2378  printf("BSM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2379  Rcor[0][0],Rcor[0][1],Rcor[0][2],Rcor[0][3]);
2380  printf("BSM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2381  Rcor[1][0],Rcor[1][1],Rcor[1][2],Rcor[1][3]);
2382  printf("BSM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2383  Rcor[2][0],Rcor[2][1],Rcor[2][2],Rcor[2][3]);
2384  printf("BSM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2385  Rcor[3][0],Rcor[3][1],Rcor[3][2],Rcor[3][3]);
2386  */
2387  xsec=Rcor[0][0];
2388  Rcor[0][0]=1.0;
2389  // TauSpinner and Korchin convention change is performed in routine dipolgamma_ which provides matrix R_ij in C++ conventions
2390  // 0123 <-> txyz
2391 
2392 
2393  /*
2394  // remnant print from the code infancy
2395  DEBUG
2396  (
2397 
2398  printf("dipol gamma przed r11 = %13.10f \n",R11 );
2399  printf("dipol gamma przed r22 = %13.10f \n",R22 );
2400  printf("dipol gamma przed r12 = %13.10f \n",R12 );
2401  printf("dipol gamma przed r21 = %13.10f \n",R21 );
2402  )
2403  */
2404  // ERW 6.04.2025: TO JEST ZLE??? i niepotrzebne, gamma gamma uzywane tylko jezeli (ifkorch==1)
2405  // i wowczas powinno byc wypelniane macierza wg konwecji RcorExt i bez normalizacji bo juz
2406  // zrobiona na poziomie fortranowskim
2407  R11 += WID[i]*Rcor[1][1]/Rcor[0][0]; // Historical --> division not necessary? Rcor[0][0]=1.0
2408  R22 += WID[i]*Rcor[2][2]/Rcor[0][0]; // note signs in line WT =
2409  R12 += WID[i]*Rcor[1][2]/Rcor[0][0]; // 1.0+sign*HHp[2]*HHm[2] ....- RzYY*R22 ....- RzYX*R21*HHp[1]*HHm[0];
2410  R21 += WID[i]*Rcor[2][1]/Rcor[0][0]; // sign operation with F=1/-1 like for RcorExt is not necessary.
2411 
2412  if(ifkorch==1){
2413  for(int i3=0;i3<=3;i3++){
2414  for(int j=0;j<=3;j++){
2415  int i4=i3-1;
2416  int j4=j-1;
2417  if(i4==-1) i4=3;
2418  if(j4==-1) j4=3;
2419  RcorExt[i4][j4]+=WID[i]*Rcor[i3][j]; // we move back to fortran convention xyzt because of f77-tauola polarimetric vectors
2420  }
2421  }
2422  xsecsum+=xsec/xsec0*WID[i];
2423  /*
2424  printf(" -----------------------------------------------\n");
2425 
2426  //ERW: rotation should be moved to dipolgamma OR NOT??
2427  printf("BSM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2428  RcorExt[0][0],RcorExt[0][1],RcorExt[0][2],RcorExt[0][3]);
2429  printf("BSM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2430  RcorExt[1][0],RcorExt[1][1],RcorExt[1][2],RcorExt[1][3]);
2431  printf("BSM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2432  RcorExt[2][0],RcorExt[2][1],RcorExt[2][2],RcorExt[2][3]);
2433  printf("BSM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2434  RcorExt[3][0],RcorExt[3][1],RcorExt[3][2],RcorExt[3][3]);
2435 
2436  printf(" -----------------------------------------------\n");
2437  */
2438  }
2439 
2440  /*
2441  // remnant print from the code infancy
2442  DEBUG
2443  (
2444  printf("dipol gamma r11 = %13.10f \n",R11 );
2445  printf("dipol gamma r22 = %13.10f \n",R22 );
2446  printf("dipol gamma r12 = %13.10f \n",R12 );
2447  printf("dipol gamma r21 = %13.10f \n",R21 );
2448  )
2449  */
2450 
2451  }
2452  }
2453  if(ifkorch==1){ // Act on RcorExt after sum of all partonic contributions.
2454  for (int k0=0;k0<4;k0++){
2455  for (int k1=0;k1<4;k1++){
2456  double F=1.0;
2457  if(k0==1) F=-1.0; // possibly combination of two operation:
2458  // 1. sign issue equivalent to overall sign in front of HHP
2459  // 2. HHP pi-angle rotation around Y-axis
2460  // difference between KORLB-KKMC conventions of z axis and this of TauSpinner.
2461  // Checks for that are pending.
2462  // WARNING; ugly inconsisten fix, it should be through HHP rotation and sign adjustment not on R_ij.
2463  // ERW 6.04.25: absorb change in sign here to have RcorExt[k0][k1] ready ZBW but this may be about sign for HHp
2464 
2465  RcorExt[k0][k1]=RcorExt[k0][k1]*F;
2466  }
2467  }
2468  }
2469 
2470  //printf("BSM R (0,1,2,3)= %13.10f %13.10f %13.10f %13.10f %13.10f %13.10f %13.10f \n",
2471  // RcorExt[0][0],RcorExt[1][1],RcorExt[2][2],RcorExt[3][3], RcorExt[1][2],RcorExt[1][3],RcorExt[2][3]);
2472 
2473 
2474  if(ifkorch==1){
2475 
2476  if(relWTnonSM!=0) WTnonSM=xsecsum;
2477  if(relWTnonSM==0) {WTnonSM=xsecsum;
2478  std::cout << "input not supported: relWTnonSM=" <<relWTnonSM<< " ifkorch= "<<ifkorch << std::endl;
2479  std::cout << "X-sects nonSM2="<<nonSM2<<" sum= " << sum << " WTnonSM-korch= "<< xsecsum <<std::endl;
2480  }
2481  }
2482  // std::cout << "X-sects-korch (nonSM2="<<nonSM2<<")="<< WTnonSM <<std::endl;
2483  // Calculations are prepared only for pp collision.
2484  // Otherwise pol = 0.0
2485  if(!Ipp) pol=0.0;
2486 
2487  return pol;
2488 }
2489 
2490 /*******************************************************************************
2491  Check if pdg's of particles in the vector<Particle>&particles match the
2492  of (p1,p2,p3,p4,p5,p6)
2493 
2494  Returns true if 'particles' contain all of the listed pdgid-s.
2495  If it does - vector<Particle>&particles will be sorted in the order
2496  as listed pdgid-s.
2497 
2498  It is done so the order of particles is the same as the order used by
2499  TAUOLA Fortran routines.
2500 *******************************************************************************/
2501 bool channelMatch(vector<Particle> &particles, int p1, int p2, int p3, int p4, int p5, int p6)
2502 {
2503  // Copy pdgid-s of all particles
2504  vector<int> list;
2505 
2506  for(unsigned int i=0;i<particles.size();i++) list.push_back(particles[i].pdgid());
2507 
2508  // Create array out of pdgid-s
2509  int p[6] = { p1, p2, p3, p4, p5, p6 };
2510 
2511  // 1) Check if 'particles' contain all pdgid-s on the list 'p'
2512 
2513  for(int i=0;i<6;i++)
2514  {
2515  // if the pdgid is zero - finish
2516  if(p[i]==0) break;
2517 
2518  bool found = false;
2519 
2520  for(unsigned int j=0;j<list.size(); j++)
2521  {
2522  // if pdgid is found - erese it from the list and search for the next one
2523  if(list[j]==p[i])
2524  {
2525  found = true;
2526  list.erase(list.begin()+j);
2527  break;
2528  }
2529  }
2530 
2531  if(!found) return false;
2532  }
2533 
2534  // if there are more particles on the list - there is no match
2535  if(list.size()!=0) return false;
2536 
2537 
2538  // 2) Rearrange particles to match the order of pdgid-s listed in array 'p'
2539 
2540  vector<Particle> newList;
2541 
2542  for(int i=0;i<6;i++)
2543  {
2544  // if the pdgid is zero - finish
2545  if(p[i]==0) break;
2546 
2547  for(unsigned int j=0;j<particles.size(); j++)
2548  {
2549  // if pdgid is found - copy it to new list and erese from the old one
2550  if(particles[j].pdgid()==p[i])
2551  {
2552  newList.push_back(particles[j]);
2553  particles.erase(particles.begin()+j);
2554  break;
2555  }
2556  }
2557  }
2558 
2559  particles = newList;
2560 
2561  return true;
2562 }
2563 
2564 /*******************************************************************************
2565  Prints out two vertices:
2566  W -> tau, nu_tau
2567  tau -> tau_daughters
2568 *******************************************************************************/
2569 void print(Particle &W, Particle &nu_tau, Particle &tau, vector<Particle> &tau_daughters) {
2570 
2571  nu_tau.print();
2572  tau .print();
2573 
2574  double px=nu_tau.px()+tau.px();
2575  double py=nu_tau.py()+tau.py();
2576  double pz=nu_tau.pz()+tau.pz();
2577  double e =nu_tau.e ()+tau.e ();
2578 
2579  // Print out sum of tau and nu_tau and also W momentum for comparison
2580  cout<<"--------------------------------------------------------------------------------------------------------"<<endl;
2581  Particle sum1(px,py,pz,e,0);
2582  sum1.print();
2583  W .print();
2584 
2585  cout<<endl;
2586 
2587  // Print out tau daughters
2588  for(unsigned int i=0; i<tau_daughters.size();i++) tau_daughters[i].print();
2589 
2590  // Print out sum of tau decay products, and also tau momentum for comparison
2591  cout<<"--------------------------------------------------------------------------------------------------------"<<endl;
2592  Particle *sum2 = vector_sum(tau_daughters);
2593  sum2->print();
2594  tau.print();
2595  cout<<endl;
2596 
2597  delete sum2;
2598 }
2599 
2600 /*******************************************************************************
2601  Sums all 4-vectors of the particles on the list
2602 *******************************************************************************/
2603 Particle *vector_sum(vector<Particle> &x) {
2604 
2605  double px=0.0,py=0.0,pz=0.0,e=0.0;
2606 
2607  for(unsigned int i=0; i<x.size();i++)
2608  {
2609  px+=x[i].px();
2610  py+=x[i].py();
2611  pz+=x[i].pz();
2612  e +=x[i].e();
2613  }
2614 
2615  Particle *sum = new Particle(px,py,pz,e,0);
2616  return sum;
2617 }
2618 
2619 } // namespace TauSpinner
2620 
Tauolapp::TauolaParticlePair::recalculateRij
void recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta)
Definition: TauolaParticlePair.cxx:985
TauSpinner::plweight
double plweight(int ide, double svar, double costhe)
Definition: nonSM.cxx:185
TauSpinner::sin2W
double sin2W(int FLAV)
Definition: EWtables.cxx:453
TauSpinner::plzap2
double plzap2(int ide, int idf, double svar, double costhe)
Definition: nonSM.cxx:125
TauSpinner::initialize_GSW_norm
void initialize_GSW_norm(int _keyGSW, double _ReGSW1, double _ImGSW1, double _ReGSW2, double _ImGSW2, double _ReGSW3, double _ImGSW3, double _ReGSW4, double _ImGSW4, double _ReGSW5, double _ImGSW5, double _ReGSW6, double _ImGSW6)
Definition: tau_reweight_lib.cxx:245
TauSpinner::calculateWeightFromParticlesH
double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector< SimpleParticle > &sp_tau1_daughters, vector< SimpleParticle > &sp_tau2_daughters)
Definition: tau_reweight_lib.cxx:726
TauSpinner::Amz
double Amz(int FLAV)
Definition: EWtables.cxx:425
TauSpinner::channelMatch
bool channelMatch(vector< Particle > &particles, int p1, int p2=0, int p3=0, int p4=0, int p5=0, int p6=0)
Definition: tau_reweight_lib.cxx:2501
TauSpinner::Gamz
double Gamz(int FLAV)
Definition: EWtables.cxx:439
TauSpinner::disth_
double disth_(double *SVAR, double *COSTHE, int *TA, int *TB)
TauSpinner::calculateWeightFromParticlesWorHpn
double calculateWeightFromParticlesWorHpn(SimpleParticle &W, SimpleParticle &tau, SimpleParticle &nu_tau, vector< SimpleParticle > &tau_daughters)
Definition: tau_reweight_lib.cxx:632
TauSpinner::getZgamMultipliersTR
void getZgamMultipliersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
Definition: tau_reweight_lib.cxx:553
TauSpinner::getWtamplitP
double getWtamplitP()
Definition: tau_reweight_lib.cxx:605
TauSpinner::getZgamParametersTR
void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
Definition: tau_reweight_lib.cxx:543
TauSpinner::getLongitudinalPolarization
double getLongitudinalPolarization(double, SimpleParticle &, SimpleParticle &)
Definition: tau_reweight_lib.cxx:1926
TauSpinner::getHiggsParameters
void getHiggsParameters(double *mass, double *width, double *normalization)
Definition: tau_reweight_lib.cxx:566
TauSpinner::setZgamMultipliersTR
void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx)
Definition: tau_reweight_lib.cxx:532
TauSpinner::setFrameType
void setFrameType(int _FrameType)
Definition: tau_reweight_lib.cxx:476
TauSpinner::EWFACT
complex< double > EWFACT(int FLAV, int NO, double s, double costhe)
Definition: EWtables.cxx:356
TauSpinner::vector_sum
Particle * vector_sum(vector< Particle > &x)
Definition: tau_reweight_lib.cxx:2603
TauSpinner::setNonSMkey
void setNonSMkey(int key)
Definition: tau_reweight_lib.cxx:590
TauSpinner::getWtNonSM
double getWtNonSM()
Definition: tau_reweight_lib.cxx:598
TauSpinner::calculateHH
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
Definition: tau_reweight_lib.cxx:1176
TauSpinner::setHiggsParametersTR
void setHiggsParametersTR(double Rxx, double Ryy, double Rxy, double Ryx)
Definition: tau_reweight_lib.cxx:507
TauSpinner::GSWadapt
void GSWadapt(int ID, double S, double cost, double GSWr[7], double GSWi[7])
Definition: tau_reweight_lib.cxx:278
TauSpinner::setRelWTnonSM
void setRelWTnonSM(int _relWTnonSM)
Definition: tau_reweight_lib.cxx:487
TauSpinner::print
void print(Particle &W, Particle &nu_tau, Particle &tau, vector< Particle > &tau_daughters)
Definition: tau_reweight_lib.cxx:2569
TauSpinner::initialize_spinner
void initialize_spinner(bool _Ipp, int _Ipol, int _nonSM2, int _nonSMN, double _CMSENE)
Definition: tau_reweight_lib.cxx:430
TauSpinner::prepareKinematicForHH
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
Definition: tau_reweight_lib.cxx:1079
TauSpinner::getZgamParametersL
void getZgamParametersL(double &Rzx, double &Rzy, double &Rtx, double &Rty, double &Rtz)
Definition: tau_reweight_lib.cxx:520
TauSpinner::setHiggsParameters
void setHiggsParameters(int jak, double mass, double width, double normalization)
Definition: tau_reweight_lib.cxx:496
TauSpinner::getTauSpin
double getTauSpin()
Definition: tau_reweight_lib.cxx:613
TauSpinner::getWtamplitM
double getWtamplitM()
Definition: tau_reweight_lib.cxx:606
Generated by doxygen 1.9.1