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.2a"<<endl;
441  cout<<" -----------------"<<endl;
442  cout<<" 25.Apr.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 &Rzz, double &Rtx, double &Rty, double &Rtz)
521 {
522  Rzx = R31;
523  Rzy = R32;
524  Rzz = R33;
525  Rtx = R41;
526  Rty = R42;
527  Rtz = R43;
528 }
529 
530 /*******************************************************************************
531  Set coefficients for transverse components of Z/gamma spin density matrix multipliers
532 *******************************************************************************/
533 void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx)
534 {
535  RzXX = Rxx;
536  RzYY = Ryy;
537  RzXY = Rxy;
538  RzYX = Ryx;
539 }
540 
541 /*******************************************************************************
542  Get transverse components of Z/gamma spin density matrix
543 *******************************************************************************/
544 void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
545 {
546  Rxx = R11;
547  Ryy = R22;
548  Rxy = R12;
549  Ryx = R21;
550 }
551 /*******************************************************************************
552  Get transverse components of Z/gamma spin density matrix multipliers
553 *******************************************************************************/
554 void getZgamMultipliersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
555 {
556  Rxx = RzXX;
557  Ryy = RzYY;
558  Rxy = RzXY;
559  Ryx = RzYX;
560 }
561 
562 
563 
564 /*******************************************************************************
565  Get Higgs mass, width and normalization of Higgs born function
566 *******************************************************************************/
567 void getHiggsParameters(double *mass, double *width, double *normalization)
568 {
569  *mass = XMH;
570  *width = XGH;
571  *normalization = Xnorm;
572 }
573 
574 /*******************************************************************************
575  Set type of spin treatment used in the sample
576  Ipol = 0 sample was not polarised
577  Ipol = 1 sample was having complete longitudinal spin effects
578  Ipol = 2 sample was featuring longitudinal spin correlations only,
579  but not dependence on polarisation due to couplings of the Z
580  Ipol = 3 as in previous case, but only angular dependence of spin polarisation
581  was missing in the sample
582 *******************************************************************************/
583 void setSpinOfSample(int _Ipol)
584 {
585  Ipol = _Ipol;
586 }
587 
588 /*******************************************************************************
589  Turn nonSM calculation of Born cross-section on/off
590 *******************************************************************************/
591 void setNonSMkey(int _key)
592 {
593  nonSM2 = _key;
594 }
595 
596 /*******************************************************************************
597  Get nonSM weight
598 *******************************************************************************/
599 double getWtNonSM()
600 {
601  return WTnonSM;
602 }
603 /*******************************************************************************
604  Get weights for tau+ tau- decay matrix elements
605 *******************************************************************************/
606 double getWtamplitP(){return WTamplitP;}
607 double getWtamplitM(){return WTamplitM;}
608 
609 /*******************************************************************************
610  Get tau spin (helicities of tau+ tau- are 100% correlated)
611  Use after sample is reweighted to obtain information on attributed tau
612  longitudinal spin projection.
613 *******************************************************************************/
614 double getTauSpin()
615 {
616  return Polari;
617 }
618 
619 /*******************************************************************************
620  Calculate weights, case of event record vertex like W -> tau nu_tau decay.
621  Function for W+/- and H+/-
622 
623  Determines decay channel, calculates all necessary components for
624  calculation of all weights, calculates weights.
625  Input: X four momentum may be larger than sum of tau nu_tau, missing component
626  is assumed to be QED brem
627  Hidden input: none
628  Hidden output: WTamplitM or WTamplitP of tau decay (depending on tau charge)
629  NOTE: weight for sp_X production matrix elements is not calculated
630  for decays of charged intermediate W+-/H+-!
631  Explicit output: WT spin correlation weight
632 *******************************************************************************/
633 double calculateWeightFromParticlesWorHpn(SimpleParticle &sp_X, SimpleParticle &sp_tau, SimpleParticle &sp_nu_tau, vector<SimpleParticle> &sp_tau_daughters)
634 {
635  // Create Particles from SimpleParticles
636 
637  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
638  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
639  Particle nu_tau(sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
640 
641  vector<Particle> tau_daughters;
642 
643  // tau pdgid
644  int tau_pdgid = sp_tau.pdgid();
645 
646  // Create vector of tau daughters
647  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
648  {
649  Particle pp(sp_tau_daughters[i].px(),
650  sp_tau_daughters[i].py(),
651  sp_tau_daughters[i].pz(),
652  sp_tau_daughters[i].e(),
653  sp_tau_daughters[i].pdgid() );
654 
655  tau_daughters.push_back(pp);
656  }
657 
658  double phi2 = 0.0, theta2 = 0.0;
659 
660  // To calcluate matrix elements TAUOLA need tau decay products in tau rest-frame: we first boost all products
661  // 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),
662  // then rotate to have neutrino from tau decay along z axis;
663  // calculated for that purpose angles phi2, theta2 are stored for rotation back of HH
664  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
665 
666 
667  // Identify decay channel and then calculate polarimetric vector HH; calculates also WTamplit
668  double *HH = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
669 
670  double sign = 1.0; // tau from W is 100 % polarized, also from charged Higgs (but with opposite sign)
671  if ( abs(sp_X.pdgid()) == 24 ) { sign= 1.0; }
672  else if( abs(sp_X.pdgid()) == 37 ) { sign=-1.0; }
673  else
674  {
675  cout<<"wrong sp_W/H.pdgid()="<<sp_X.pdgid()<<endl;
676  exit(-1);
677  }
678  if (sp_X.pdgid() > 0 )
679  {WTamplitM = WTamplit;} // tau- decay matrix element^2, spin averaged.
680  else
681  {WTamplitP = WTamplit;} // tau+ decay matrix element^2, spin averaged.
682 
683  // spin correlation weight. Tau production matrix element is represented by `sign'
684  double WT = 1.0+sign*HH[2]; // [2] means 'pz' component
685 
686  // Print out some info about the channel
687  DEBUG
688  (
689  cout<<tau_pdgid<<" -> ";
690  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
691  cout<<" (HH: "<<HH[0]<<" "<<HH[1]<<" "<<HH[2]<<" "<<HH[3]<<") WT: "<<WT<<endl;
692  )
693 
694  // TP:31Nov2013 checks of possible problems: weight outside theoretically allowed range
695  if (WT<0.0) {
696  printf("TauSpinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 0.0\n",WT);
697  WT = 0.0;
698  }
699 
700  if (WT>2.0) {
701  printf("Tauspinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 2.0\n",WT);
702  WT = 2.0;
703  }
704 
705  delete HH;
706 
707  return WT;
708 }
709 
710 /*******************************************************************************
711  Calculate weights, case of event record vertex like Z/gamma/H ... -> tau tau decay.
712 
713  Determine decay channel, calculates all necessary components for
714  calculation of all weights, calculates weights.
715 
716  Input: X four momentum may be larger than sum of tau1 tau2, missing component
717  is assumed to be QED brem
718 
719  Hidden input: relWTnonS, nonSM2 (used only in getLongitudinalPolarization(...) )
720  Hidden output: WTamplitM or WTamplitP of tau1 tau2 decays
721  weight for sp_X production matrix element^2 is calculated inside
722  plzap2
723  Polari - helicity attributed to taus, 100% correlations between tau+ and tau-
724  WTnonSM
725  Explicit output: WT spin correlation weight
726 *******************************************************************************/
727 double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector<SimpleParticle> &sp_tau1_daughters, vector<SimpleParticle> &sp_tau2_daughters)
728 {
729  // cout << "sp_tau1_daughters = " << sp_tau1_daughters.size() << endl;
730  // cout << "sp_tau2_daughters = " << sp_tau2_daughters.size() << endl;
731  SimpleParticle sp_tau;
732  SimpleParticle sp_nu_tau;
733  vector<SimpleParticle> sp_tau_daughters;
734 
735  // First we calculate HH for tau+
736  // We enforce that sp_tau is tau+ so the 'nu_tau' is tau-
737  if (sp_tau1.pdgid() == -15 )
738  {
739  sp_tau = sp_tau1;
740  sp_nu_tau = sp_tau2;
741  sp_tau_daughters = sp_tau1_daughters;
742  }
743  else
744  {
745  sp_tau = sp_tau2;
746  sp_nu_tau = sp_tau1;
747  sp_tau_daughters = sp_tau2_daughters;
748  }
749 
750  double *HHp, *HHm;
751 
752  // We use artificial if(true){... } construction to separate namespace for tau+ and tau-
753  if(true)
754  {
755  // Create Particles from SimpleParticles
756  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
757  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
758  Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
759 
760  vector<Particle> tau_daughters;
761 
762  // tau pdgid
763  int tau_pdgid = sp_tau.pdgid();
764 
765  // Create list of tau daughters
766  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
767  {
768  Particle pp(sp_tau_daughters[i].px(),
769  sp_tau_daughters[i].py(),
770  sp_tau_daughters[i].pz(),
771  sp_tau_daughters[i].e(),
772  sp_tau_daughters[i].pdgid() );
773 
774  tau_daughters.push_back(pp);
775  }
776 
777  double phi2 = 0.0, theta2 = 0.0;
778 
779  // To calculate matrix elements TAUOLA need tau decay products in tau rest-frame: we first boost all products
780  // to tau rest frame (tau other tau of Z/H decay along z-axis, intermediate step for boost is tau-tau pair rest-frame),
781  // then rotate to have neutrino from tau decay along z axis;
782  // calculated for that purpose angles phi2, theta2 are stored for rotation back of HHp
783  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
784 
785 
786 
787 
788  // Identify decay channel and then calculate polarimetric vector HH; calculates also WTamplit
789  HHp = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
790 
791  DEBUG
792  (
793  cout<<tau_pdgid<<" -> ";
794  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
795  cout<<" (HHp: "<<HHp[0]<<" "<<HHp[1]<<" "<<HHp[2]<<" "<<HHp[3]<<") ";
796  cout<<endl;
797  )
798 
799  WTamplitP = WTamplit;
800  } // end of if(true); for tau+
801 
802 
803 
804  // Second we calculate HH for tau-
805  // We enforce that sp_tau is tau- so the 'nu_tau' is tau+
806  if(sp_tau1.pdgid() == 15 )
807  {
808  sp_tau = sp_tau1;
809  sp_nu_tau = sp_tau2;
810  sp_tau_daughters = sp_tau1_daughters;
811  }
812  else
813  {
814  sp_tau = sp_tau2;
815  sp_nu_tau = sp_tau1;
816  sp_tau_daughters = sp_tau2_daughters;
817  }
818 
819  // We use artificial if(true){... } construction to separate namespace for tau+ and tau-
820  if(true)
821  {
822  // Create Particles from SimpleParticles
823  Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
824  Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
825  Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
826 
827  vector<Particle> tau_daughters;
828 
829  // tau pdgid
830  int tau_pdgid = sp_tau.pdgid();
831 
832  // Create list of tau daughters
833  for(unsigned int i=0; i<sp_tau_daughters.size(); i++)
834  {
835  Particle pp(sp_tau_daughters[i].px(),
836  sp_tau_daughters[i].py(),
837  sp_tau_daughters[i].pz(),
838  sp_tau_daughters[i].e(),
839  sp_tau_daughters[i].pdgid() );
840 
841  tau_daughters.push_back(pp);
842  }
843 
844  double phi2 = 0.0, theta2 = 0.0;
845 
846 
847  // To calculate matrix elements TAUOLA need tau decay products in tau rest-frame: we first boost all products
848  // to tau rest frame (tau other tau of Z/H decay along z-axis, intermediate step for boost is tau-tau pair rest-frame),
849  // then rotate to have neutrino from tau decay along z axis;
850  // calculated for that purpose angles phi2, theta2 are stored for rotation back of HHm
851  prepareKinematicForHH (tau, nu_tau, tau_daughters, &phi2, &theta2);
852 
853 
854  // Identify decay channel and then calculate polarimetric vector HHm; calculates also WTamplit
855  HHm = calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
856 
857  DEBUG
858  (
859  cout<<tau_pdgid<<" -> ";
860  for(unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<" ";
861  cout<<" (HHm: "<<HHm[0]<<" "<<HHm[1]<<" "<<HHm[2]<<" "<<HHm[3]<<") ";
862  cout<<endl;
863  )
864 
865  WTamplitM = WTamplit;
866  } // end of if(true); for tau-
867 
868  // CALCULATION OF PRODUCTION MATRIX ELEMENTS, THEN SPIN WEIGHTS AND FURTHER HIDDEN OUTPUTS
869 
870  // sign, in ultrarelativistic limit it is component of spin density matrix:
871  // longitudinal spin correlation for intermediate vector state gamma*,
872  // for Z etc. sign= +1; for Higgs, other scalars, pseudoscalars sign= -1
873  double sign = 1.0;
874  if(sp_X.pdgid() == 25) { sign=-1.0; }
875  if(sp_X.pdgid() == 36) { sign=-1.0; }
876  if(sp_X.pdgid() ==553) { sign=-1.0; } // upsilon(1s) can be treated as scalar
877 
878  double WT = 0.0;
879 
880  Polari = 0.0;
881  if(sign == -1.0) // Case of scalar
882  {
883  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();
884  IfHiggs=true; //global variable
885  double pol = getLongitudinalPolarization(S, sp_tau, sp_nu_tau); // makes sense only for nonSM otherwise NaN
886 
887  if(nonSM2==1) // WARNING it is for spin=2 resonance!!
888  {
889 
890  double corrX2;
891  double polX2;
892 
893  // NOTE: in this case, sp_nu_tau is the 2nd tau
894  // nonSMHcorrPol(S, sp_tau, sp_nu_tau, &corrX2, &polX2); // for future use
895  // WARNING: may be for polX2*HHm[2] we need to fix sign!
896  polX2=pol;
897  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
898 
899  // 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
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  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];
902 
903  // we separate cross section into helicity parts. From this, we attribute helicity states to taus: ++ or --
904  double RRR = Tauola::randomDouble();
905  Polari=1.0;
906  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;
907  }
908  else // case of Higgs
909  {
910  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];
911 
912  // we separate cross section into helicity parts. From this, we attribute helicity states to taus: +- or -+
913  double RRR = Tauola::randomDouble();
914  Polari=1.0;
915  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;
916  }
917  }
918  else // Case of Drell Yan
919  {
920 
921  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();
922 
923  // Get Z polarization
924  // ( Variable names are misleading! sp_tau is tau+ and sp_nu_tau is tau- )
925 
926  IfHiggs=false;
927 
928  double pol = getLongitudinalPolarization(S, sp_tau, sp_nu_tau);
929  // 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
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  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];
932  /*
933  // HHm[3] = 1.0; // initialized elsewhere
934  // HHp[3] = 1.0;
935  printf("--WT = %13.10f \n", WT);
936  printf("---------polarimetric vectors---------------------------\n");
937  printf(" HHm(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
938  HHm[0], HHm[1], HHm[2], HHm[3] );
939  printf(" HHp(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
940  HHp[0], HHp[1], HHp[2], HHp[3] );
941  printf("---------RcorOld (adjusted)--------------------\n");
942  printf("SM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
943  R11,R12,0.0,0.0);
944  printf("SM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
945  - R21,- R22,0.0,0.0);
946  printf("SM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
947  0.0,0.0,1.0,pol);
948  printf("SM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
949  0.0,0.0,pol,1.0);
950  printf("-------------------------------------------\n");
951  */
952  double zero=0.0;
953  double jeden=1.0;
954  // printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",R11,R12,zero,zero);
955  //printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",R21,R22,zero,zero);
956  //printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",zero,zero,sign,pol);
957  //printf("old %13.10f. %13.10f. %13.10f. %13.10f. \n",zero,zero,pol,jeden);
958  // in future we may need extra factor for wt which is
959  // F=PLWEIGHT(IDE,IDF,SVAR,COSTHE,1)
960  // but it is then non standard version of the code.
961 
962  // to correct when in the sample only spin corr. are in, but no polarization
963  if(Ipol==2) WT = WT/(1.0+sign*HHp[2]*HHm[2]);
964 
965  // to correct sample when corr. are in, pol. is in, but angular
966  // dependence of pol is missing.
967  if(Ipol==3)
968  {
969  // valid for configurations close to Z peak, otherwise approximation is at use.
970  double polp1 = plzap2(11,sp_tau.pdgid(),S,0.0);
971  double pol1 =(2*(1-polp1)-1) ;
972  WT = WT/(1.0+sign*HHp[2]*HHm[2]+pol1*HHp[2]+pol1*HHm[2]);
973  }
974  double WTtest;
975  if(ifkorch==1){
976  WTtest=0;
977  HHp[3]=1.0; // WARNING: temporary fix, til now HHp/m[3] were not used.
978  HHm[3]=1.0; // They were not always initialized by calculateHH()!
979 
980  for (int k0=0;k0<4;k0++){
981  for (int k1=0;k1<4;k1++){
982  double F=1.0;
983  // coefficient to switch off some elements of correlation matrix
984  if(k0==0&&k1==0) F=RzXX;
985  if(k0==1&&k1==1) F=RzYY;
986  if(k0==0&&k1==1) F=RzXY;
987  if(k0==1&&k1==0) F=RzYX;
988  WTtest=WTtest+RcorExt[k0][k1]*HHp[k0]*HHm[k1]*F;
989  }
990  }
991  /*
992  printf("--WT = %13.10f \n", WTtest);
993  printf("---------polarimetric vectors---------------------------\n");
994  printf(" HHm(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
995  HHm[0], HHm[1], HHm[2], HHm[3] );
996  printf(" HHp(0, 1, 2, 3) = %13.10f %13.10f %13.10f %13.10f \n",
997  HHp[0], HHp[1], HHp[2], HHp[3] );
998  printf("---------RcorExt (adjusted)---------------------------\n");
999  printf("SM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1000  RcorExt[0][0], RcorExt[0][1], RcorExt[0][2],RcorExt[0][3]);
1001  printf("SM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1002  RcorExt[1][0], RcorExt[1][1], RcorExt[1][2], RcorExt[1][3]);
1003  printf("SM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1004  RcorExt[2][0],RcorExt[2][1],RcorExt[2][2],RcorExt[2][3]);
1005  printf("SM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
1006  RcorExt[3][0],RcorExt[3][1],RcorExt[3][2],RcorExt[3][3]);
1007  printf("-------------------------------------------\n");
1008  */
1009  R11 = RcorExt[0][0];
1010  R12 = RcorExt[0][1];
1011  R21 = RcorExt[1][0];
1012  R22 = RcorExt[1][1];
1013  R31 = RcorExt[2][0];
1014  R32 = RcorExt[2][1];
1015  R33 = RcorExt[2][2];
1016  R41 = RcorExt[3][0];
1017  R42 = RcorExt[3][1];
1018  R43 = RcorExt[3][2];
1019 
1020  // WTtest=WTtest; this line is to be activated/extened/modified at some tests.
1021  // cout << " why was minus sign in front of RzYY*R22*HHp[1]*HHm[1], WT (D-Y) formula?" << endl;
1022  // cout << " weight (non-korch) WT(nonSM="<<nonSM2<<")= " << WT << " spin-WTkorch= " << WTtest << " "<< endl;
1023  if(nonSM2==1) WT=WTtest; // to activate case of anomalous moment
1024  if(nonSM2==0&& Ipol==0) WT=WTtest; // to activate SM variant of full spin wt onto non polarized sample.
1025 
1026  }
1027  // we separate cross section into helicity parts. From this, we attribute helicity states to taus: ++ or --
1028  double RRR = Tauola::randomDouble();
1029  Polari=1.0;
1030  // ERW 6.04.25: added line for ifkorch
1031  // printf(" pol, RcorExt[3][2] = %13.10f %13.10f \n", pol, RcorExt[3][2]);
1032  if( ifkorch==1) pol = RcorExt[3][2];
1033  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;
1034  }
1035 
1036  // Print out some info about the channel
1037  DEBUG( cout<<" WT: "<<WT<<endl; )
1038 
1039  if (WT<0.0) {
1040  printf("Tauspinner::calculateWeightFromParticlesH WT is: %13.10f. Setting WT = 0.0\n",WT);
1041  WT = 0.0; // SwB:23Feb2013
1042  }
1043 
1044  if (WT>4.0) {
1045  printf("Tauspinner::calculateWeightFromParticlesH WT is: %13.10f. Setting WT = 4.0\n",WT);
1046  WT = 4.0; // SwB:23Feb2013
1047  }
1048 
1049  if( WT>4.0 || WT<0.0)
1050  {
1051  cout<<"Tauspinner::calculateWeightFromParticlesH ERROR: Z/gamma* or H, and WT not in range [0,4]."<<endl;
1052  exit(-1);
1053  }
1054 
1055  if (sign==-1.0 && nonSM2!=1) {
1056  if (WT>2.0) {
1057  WT = 2.0; // SwB:26Feb2013
1058  cout << "Tauspinner::calculateWeightFromParticlesH Setting WT to be 2.0" << endl;
1059  }
1060  }
1061 
1062  if( sign==-1.0 && (WT>2.0 || WT<0.0) && nonSM2!=1)
1063  {
1064  cout<<"Tauspinner::calculateWeightFromParticlesH ERROR: H and WT not in range [0,2]."<<endl;
1065  exit(-1);
1066  }
1067 
1068  delete[] HHp;
1069  delete[] HHm;
1070 
1071  return WT;
1072 }
1073 
1074 /*******************************************************************************
1075  Prepare kinematics for HH calculation
1076 
1077  Boost particles to effective bozon rest frame, and rotate them so that tau is on Z axis.
1078  Then rotate again with theta2 phi2 so neutrino from tau decay is along Z.
1079 *******************************************************************************/
1080 void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector<Particle> &tau_daughters, double *phi2, double *theta2)
1081 {
1082  Particle P_QQ( tau.px()+nu_tau.px(), tau.py()+nu_tau.py(), tau.pz()+nu_tau.pz(), tau.e()+nu_tau.e(), 0 );
1083 
1084  //cout<<endl<<"START: "<<endl;
1085  //print(P_QQ,nu_tau,tau,tau_daughters);
1086 
1087  // 1) boost tau, nu_tau and tau daughters to rest frame of P_QQ
1088 
1089  tau.boostToRestFrame(P_QQ);
1090  nu_tau.boostToRestFrame(P_QQ);
1091 
1092  for(unsigned int i=0; i<tau_daughters.size();i++)
1093  tau_daughters[i].boostToRestFrame(P_QQ);
1094 
1095  //cout<<endl<<"AFTER 1: "<<endl;
1096  //print(P_QQ,nu_tau,tau,tau_daughters);
1097 
1098  // 2) Rotate tau, nu_tau~, tau daughters to frame where tau is along Z
1099  // We set accompanying neutino in direction of Z+
1100 
1101  double phi = tau.getAnglePhi();
1102 
1103  tau.rotateXY(-phi);
1104 
1105  double theta = tau.getAngleTheta();
1106 
1107  tau.rotateXZ(M_PI-theta);
1108 
1109  nu_tau.rotateXY(-phi );
1110  nu_tau.rotateXZ(M_PI-theta);
1111 
1112  for(unsigned int i=0; i<tau_daughters.size();i++)
1113  {
1114  tau_daughters[i].rotateXY(-phi );
1115  tau_daughters[i].rotateXZ(M_PI-theta);
1116  }
1117 
1118  //cout<<endl<<"AFTER 2: "<<endl;
1119  //print(P_QQ,nu_tau,tau,tau_daughters);
1120 
1121  // 3) boost tau_daughters along Z to rest frame of tau
1122 
1123  for(unsigned int i=0; i<tau_daughters.size();i++)
1124  tau_daughters[i].boostAlongZ(-tau.pz(),tau.e());
1125 
1126  //cout<<endl<<"AFTER 3: "<<endl;
1127  //print(P_QQ,nu_tau,tau,tau_daughters);
1128 
1129  // 4) Now rotate tau daughters second time
1130  // so that nu_tau (from tau daughters list) is on Z axis
1131 
1132  // We can not be sure if tau_daughters[0] is neutrino !!!
1133  // That is the case, for tauola generated samples, but not in general!
1134  unsigned int i_stored = 0;
1135 
1136  *phi2=0;
1137  *theta2=0;
1138 
1139  for(unsigned int i=0; i<tau_daughters.size();i++)
1140  {
1141  if(abs(tau_daughters[i].pdgid())==16){
1142  *phi2 = tau_daughters[i].getAnglePhi();
1143 
1144  tau_daughters[i].rotateXY( -(*phi2) );
1145 
1146  *theta2 = tau_daughters[i].getAngleTheta();
1147 
1148  tau_daughters[i].rotateXZ( -(*theta2) );
1149 
1150  i_stored = i;
1151  break;
1152  }
1153  }
1154 
1155  for(unsigned int i=0; i<tau_daughters.size();i++)
1156  {
1157  if(i != i_stored) {
1158  tau_daughters[i].rotateXY( -(*phi2) );
1159  tau_daughters[i].rotateXZ( -(*theta2) );
1160  }
1161  }
1162 
1163  //cout<<endl<<"AFTER 4: "<<endl;
1164  //print(P_QQ,nu_tau,tau,tau_daughters);
1165 }
1166 
1167 /*******************************************************************************
1168  Calculates polarimetric vector HH of the tau decay.
1169 
1170  HH[3] is timelike because we use FORTRAN methods to calculate HH.
1171  First decide what is the channel. After that, 4-vectors
1172  are moved to tau rest frame of tau.
1173  Polarimetric vector HH is then rotated using angles phi and theta.
1174 
1175  Order of the tau decay products does not matter (it is adjusted)
1176 *******************************************************************************/
1177 double* calculateHH(int tau_pdgid, vector<Particle> &tau_daughters, double phi, double theta)
1178 {
1179  int channel = 0;
1180  double *HH = new double[4];
1181 
1182  HH[0]=HH[1]=HH[2]=HH[3]=0.0;
1183 
1184  vector<int> pdgid;
1185 
1186  // Create list of tau daughters pdgid
1187  for(unsigned int i=0; i<tau_daughters.size(); i++)
1188  pdgid.push_back( tau_daughters[i].pdgid() );
1189 
1190  // 17.04.2014: If Tauola++ is used for generation,
1191  // jaki_.ktom may be changed to 11 at the time of storing decay to event record
1192  // (case of full spin effects).
1193  // For Tauola++ jaki_.ktom is later of no use so 11 does not create problems.
1194  // For TauSpinner processing, jaki_.ktom should always be 1
1195  // This was the problem if Tauola++ generation and TauSpinner were used simultaneously.
1196  jaki_.ktom = 1;
1197 
1198  // tau^- --> pi^- nu_tau
1199  // tau^+ --> pi^+ anti_nu_tau
1200  // tau^- --> K^- nu_tau
1201  // tau^+ --> K^+ anti_nu_tau
1202  if( pdgid.size()==2 &&
1203  (
1204  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211) ) ||
1205  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211) ) ||
1206  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321) ) ||
1207  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321) )
1208  )
1209  ) {
1210  channel = 3; // channel: numbering convention of TAUOLA
1211  if(abs(pdgid[1])==321) channel = 6;
1212  DEBUG( cout<<"Channel "<<channel<<" : "; )
1213  // PXQ=AMTAU*EPI
1214  // PXN=AMTAU*ENU
1215  // QXN=PPI(4)*PNU(4)-PPI(1)*PNU(1)-PPI(2)*PNU(2)-PPI(3)*PNU(3)
1216 
1217  // BRAK=(GV**2+GA**2)*(2*PXQ*QXN-AMPI**2*PXN)
1218  // HV(I)=-ISGN*2*GA*GV*AMTAU*(2*PPI(I)*QXN-PNU(I)*AMPI**2)/BRAK
1219 
1220  const double AMTAU = 1.777;
1221  // is mass of the Pi+- OR K+-
1222  double AMPI = sqrt(tau_daughters[1].e() *tau_daughters[1].e()
1223  -tau_daughters[1].px()*tau_daughters[1].px()
1224  -tau_daughters[1].py()*tau_daughters[1].py()
1225  -tau_daughters[1].pz()*tau_daughters[1].pz());
1226 
1227  // two-body decay is so simple, that matrix element is calculated here
1228  double PXQ=AMTAU*tau_daughters[1].e();
1229  double PXN=AMTAU*tau_daughters[0].e();
1230  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();
1231  double BRAK=(2*PXQ*QXN-AMPI*AMPI*PXN);
1232 
1233  WTamplit = (1.16637E-5)*(1.16637E-5)*BRAK/2.;//AMPLIT=(GFERMI)**2*BRAK/2. //WARNING: Note for normalisation Cabbibo angle is missing!
1234  HH[0] = AMTAU*(2*tau_daughters[1].px()*QXN-tau_daughters[0].px()*AMPI*AMPI)/BRAK;
1235  HH[1] = AMTAU*(2*tau_daughters[1].py()*QXN-tau_daughters[0].py()*AMPI*AMPI)/BRAK;
1236  HH[2] = AMTAU*(2*tau_daughters[1].pz()*QXN-tau_daughters[0].pz()*AMPI*AMPI)/BRAK;
1237  HH[3] = 1.0;
1238  }
1239 
1240  // tau^- --> pi^- pi^0 nu_tau
1241  // tau^+ --> pi^+ pi^0 anti_nu_tau
1242  // tau^- --> K^- K^0 nu_tau; K^0 may be K_L K_S too.
1243  // tau^+ --> K^+ K^0 anti_nu_tau
1244  else if( pdgid.size()==3 &&
1245  (
1246  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 111) ) ||
1247  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 111) ) ||
1248  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 311) ) ||
1249  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 311) ) ||
1250  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 310) ) ||
1251  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 310) ) ||
1252  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 130) ) ||
1253  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 130) )
1254 
1255  )
1256  ) {
1257 
1258  channel = 4;
1259  DEBUG( cout<<"Channel "<<channel<<" : "; )
1260  // PRODPQ=PT(4)*QQ(4)
1261  // PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
1262  // PRODPN=PT(4)*PN(4)
1263  // BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
1264  // HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
1265 
1266  const double AMTAU = 1.777;
1267 
1268  int MNUM = 0;
1269  if(tau_daughters[2].pdgid() != 111) { MNUM=3; channel = 22;} // sub case of decay to K-K0
1270  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1271  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1272  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1273  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1274  float AMPLIT = 0.0;
1275  float HV[4] = { 0.0 };
1276 
1277  dam2pi_( &MNUM, PT, PN, PIM1, PIM2, &AMPLIT, HV );
1278 
1279  WTamplit = AMPLIT;
1280  HH[0] = -HV[0];
1281  HH[1] = -HV[1];
1282  HH[2] = -HV[2];
1283  HH[3] = HV[3];
1284  }
1285 
1286  // tau^- --> K^- pi^0 nu_tau
1287  // tau^+ --> K^+ pi^0 anti_nu_tau
1288  // tau^- --> pi^- K_S0 nu_tau
1289  // tau^+ --> pi^+ K_S0 anti_nu_tau
1290  // tau^- --> pi^- K_L0 nu_tau
1291  // tau^+ --> pi^+ K_L0 anti_nu_tau
1292  else if( pdgid.size()==3 &&
1293  (
1294  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 130) ) ||
1295  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 130) ) ||
1296  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 310) ) ||
1297  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 310) ) ||
1298  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 311) ) ||
1299  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 311) ) ||
1300  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321, 111) ) ||
1301  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 111) )
1302  )
1303  ) {
1304 
1305  channel = 7;
1306  DEBUG( cout<<"Channel "<<channel<<" : "; )
1307  // PRODPQ=PT(4)*QQ(4)
1308  // PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
1309  // PRODPN=PT(4)*PN(4)
1310  // BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
1311  // HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
1312 
1313  const double AMTAU = 1.777;
1314 
1315  double QQ[4];
1316  QQ[0]=tau_daughters[1].e() -tau_daughters[2].e() ;
1317  QQ[1]=tau_daughters[1].px()-tau_daughters[2].px();
1318  QQ[2]=tau_daughters[1].py()-tau_daughters[2].py();
1319  QQ[3]=tau_daughters[1].pz()-tau_daughters[2].pz();
1320 
1321  double PKS[4];
1322  PKS[0]=tau_daughters[1].e() +tau_daughters[2].e() ;
1323  PKS[1]=tau_daughters[1].px()+tau_daughters[2].px();
1324  PKS[2]=tau_daughters[1].py()+tau_daughters[2].py();
1325  PKS[3]=tau_daughters[1].pz()+tau_daughters[2].pz();
1326 
1327  // orthogonalization of QQ wr. to PKS
1328  double PKSD=PKS[0]*PKS[0]-PKS[1]*PKS[1]-PKS[2]*PKS[2]-PKS[3]*PKS[3];
1329  double QQPKS=QQ[0]*PKS[0]-QQ[1]*PKS[1]-QQ[2]*PKS[2]-QQ[3]*PKS[3];
1330 
1331  QQ[0]=QQ[0]-PKS[0]*QQPKS/PKSD;
1332  QQ[1]=QQ[1]-PKS[1]*QQPKS/PKSD;
1333  QQ[2]=QQ[2]-PKS[2]*QQPKS/PKSD;
1334  QQ[3]=QQ[3]-PKS[3]*QQPKS/PKSD;
1335 
1336  double PRODPQ=AMTAU*QQ[0];
1337  double PRODNQ=tau_daughters[0].e() *QQ[0]
1338  -tau_daughters[0].px()*QQ[1]
1339  -tau_daughters[0].py()*QQ[2]
1340  -tau_daughters[0].pz()*QQ[3];
1341  double PRODPN=AMTAU*tau_daughters[0].e();
1342  double QQ2 =QQ[0]*QQ[0]-QQ[1]*QQ[1]-QQ[2]*QQ[2]-QQ[3]*QQ[3];
1343 
1344  // in this case matrix element is calculated here
1345  double BRAK=(2*PRODPQ*PRODNQ-PRODPN*QQ2);
1346 
1347  WTamplit = (1.16637E-5)*(1.16637E-5)*BRAK/2.;//AMPLIT=(GFERMI)**2*BRAK/2. WARNING: Note for normalisation Cabibbo angle is missing!
1348  HH[0]=AMTAU*(2*PRODNQ*QQ[1]-QQ2*tau_daughters[0].px())/BRAK;
1349  HH[1]=AMTAU*(2*PRODNQ*QQ[2]-QQ2*tau_daughters[0].py())/BRAK;
1350  HH[2]=AMTAU*(2*PRODNQ*QQ[3]-QQ2*tau_daughters[0].pz())/BRAK;
1351  HH[3]=1.0;
1352  }
1353 
1354  // tau^- --> e^- anti_nu_e nu_tau
1355  // tau^+ --> e^+ nu_e anti_nu_tau
1356  else if( pdgid.size()==3 &&
1357  (
1358  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 11,-12) ) ||
1359  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-11, 12) )
1360  )
1361  ) {
1362  DEBUG( cout<<"Channel 1 : "; )
1363  channel = 1;
1364  // ITDKRC=0,XK0DEC=0.01 XK[4]={0},XA[4] nu_e, QP[4] e, XN[4] neutrino tauowe, AMPLIT, HH[4]
1365  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1366 
1367  int ITDKRC = 0;
1368  double XK0DEC = 0.01;
1369  double XK[4] = { 0.0 };
1370  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1371  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1372  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1373  double AMPLIT = 0.0;
1374  double HV[4] = { 0.0 };
1375 
1376  // We fix 4-momenta of electron and electron neutrino
1377  // Since electrons have small mass, they are prone to rounding errors
1378  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.511e-3*0.511e-3);
1379  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1380 
1381  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1382 
1383  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1384  HH[0] = -HV[0];
1385  HH[1] = -HV[1];
1386  HH[2] = -HV[2];
1387  HH[3] = HV[3];
1388  }
1389 
1390  // tau^- --> e^- anti_nu_e nu_tau + gamma
1391  // tau^+ --> e^+ nu_e anti_nu_tau + gamma
1392  else if( pdgid.size()==4 &&
1393  (
1394  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 11,-12, 22) ) ||
1395  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-11, 12, 22) )
1396  )
1397  ) {
1398  DEBUG( cout<<"Channel 1b : "; )
1399  channel = 1;
1400  // ITDKRC=0,XK0DEC=0.01 XK[4] gamma, XA[4] nu_e, QP[4] e, XN[4] neutrino tau , AMPLIT, HH[4]
1401  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1402 
1403  int ITDKRC = 1;
1404  double XK0DEC = 0.01;
1405  double XK[4] = { tau_daughters[3].px(), tau_daughters[3].py(), tau_daughters[3].pz(), tau_daughters[3].e() };
1406  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1407  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1408  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1409  double AMPLIT = 0.0;
1410  double HV[4] = { 0.0 };
1411 
1412  // We fix 4-momenta of electron and electron neutrino and photon
1413  // Since electrons have small mass, they are prone to rounding errors
1414  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.511e-3*0.511e-3);
1415  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1416  XK[3] = sqrt( XK[0]*XK[0] + XK[1]*XK[1] + XK[2]*XK[2] );
1417  // XK0DEC must be smaller in TauSpinner than what was used in generation. We do not use virt. corr anyway.
1418  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]);
1419 
1420  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1421 
1422  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1423  HH[0] = -HV[0];
1424  HH[1] = -HV[1];
1425  HH[2] = -HV[2];
1426  HH[3] = HV[3];
1427  }
1428 
1429  // tau^- --> mu^- antui_nu_mu nu_tau
1430  // tau^+ --> mu^+ nu_mu anti_nu_tau
1431  else if( pdgid.size()==3 &&
1432  (
1433  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 13,-14) ) ||
1434  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-13, 14) )
1435  )
1436  ) {
1437 
1438  DEBUG( cout<<"Channel 2 : "; )
1439  channel = 2;
1440  // ITDKRC=0,XK0DEC=0.01 XK[4]={0},XA[4] nu_mu, QP[4] mu, XN[4] neutrino tauowe, AMPLIT, HH[4]
1441  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1442 
1443  int ITDKRC = 0;
1444  double XK0DEC = 0.01;
1445  double XK[4] = { 0.0 };
1446  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1447  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1448  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1449  double AMPLIT = 0.0;
1450  double HV[4] = { 0.0 };
1451 
1452  // We fix 4-momenta of muon and muon neutrino
1453  // Since muon have small mass, they are prone to rounding errors
1454  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.105659*0.105659);
1455  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1456 
1457  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1458 
1459  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1460  HH[0] = -HV[0];
1461  HH[1] = -HV[1];
1462  HH[2] = -HV[2];
1463  HH[3] = HV[3];
1464  }
1465 
1466  // tau^- --> mu^- antui_nu_mu nu_tau + gamma
1467  // tau^+ --> mu^+ nu_mu anti_nu_tau + gamma
1468  else if( pdgid.size()==4 &&
1469  (
1470  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 13,-14, 22) ) ||
1471  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16,-13, 14, 22) )
1472  )
1473  ) {
1474 
1475  DEBUG( cout<<"Channel 2b : "; )
1476  channel = 2;
1477  // ITDKRC=0,XK0DEC=0.01 XK[4] gamma, XA[4] nu_mu, QP[4] mu, XN[4] neutrino tau, AMPLIT, HH[4]
1478  // SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1479 
1480  int ITDKRC = 1;
1481  double XK0DEC = 0.01;
1482  double XK[4] = { tau_daughters[3].px(), tau_daughters[3].py(), tau_daughters[3].pz(), tau_daughters[3].e() };
1483  double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1484  double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1485  double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1486  double AMPLIT = 0.0;
1487  double HV[4] = { 0.0 };
1488 
1489  // We fix 4-momenta of muon and muon neutrino and photon
1490  // Since muons have small mass, they are prone to rounding errors
1491  QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.105659*0.105659);
1492  XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1493  XK[3] = sqrt( XK[0]*XK[0] + XK[1]*XK[1] + XK[2]*XK[2] );
1494  // XK0DEC must be smaller in TauSpinner than what was used in generation. We do not use virt. corr anyway.
1495  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]);
1496 
1497 
1498  dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &AMPLIT, HV );
1499 
1500  WTamplit = AMPLIT; // WARNING: note XK0DEC dependence is not included in normalisation
1501  HH[0] = -HV[0];
1502  HH[1] = -HV[1];
1503  HH[2] = -HV[2];
1504  HH[3] = HV[3];
1505  }
1506 
1507  // tau^- --> pi^- pi^0 pi^0 nu_tau
1508  // tau^+ --> pi^+ pi^0 pi^0 anti_nu_tau
1509  else if( pdgid.size()==4 &&
1510  (
1511  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 111, 111,-211) ) ||
1512  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 111, 111, 211) )
1513  )
1514  ) {
1515  DEBUG( cout<<"Channel 5 : "; )
1516  channel = 5;
1517  // MNUM=0, PT[4] tau, PN[4] neutrino, pi0[4], pi0[4], pi[4], AMPLIT, HH[4]
1518  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1519 
1520  const double AMTAU = 1.777;
1521  int MNUM = 0;
1522  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1523  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1524  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1525  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1526  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1527  float AMPLIT = 0.0;
1528  float HV[4] = { 0.0 };
1529 
1530  // For RChL currents one needs to define 3-pi sub-channel used
1531  chanopt_.JJ=2;
1532 
1533  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1534 
1535  WTamplit = AMPLIT;
1536  HH[0] = -HV[0];
1537  HH[1] = -HV[1];
1538  HH[2] = -HV[2];
1539  HH[3] = HV[3];
1540  }
1541 
1542  // tau^- --> pi^+ pi^- pi^- nu_tau
1543  // tau^+ --> pi^- pi^+ pi^+ anti_nu_tau
1544  else if( pdgid.size()==4 &&
1545  (
1546  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211,-211, 211) ) ||
1547  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 211,-211) )
1548  )
1549  ) {
1550  DEBUG( cout<<"Channel 5 : "; )
1551  channel = 5;
1552  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1553  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1554 
1555  const double AMTAU = 1.777;
1556  int MNUM = 0;
1557  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1558  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1559  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1560  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1561  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1562  float AMPLIT = 0.0;
1563  float HV[4] = { 0.0 };
1564 
1565  // For RChL currents one needs to define 3-pi sub-channel used
1566  chanopt_.JJ=1;
1567 
1568  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1569 
1570  WTamplit = AMPLIT;
1571  HH[0] = -HV[0];
1572  HH[1] = -HV[1];
1573  HH[2] = -HV[2];
1574  HH[3] = HV[3];
1575  }
1576 
1577  // tau^- --> K^+ pi^- pi^+ nu_tau // prepared for modes with kaons
1578  // tau^+ --> K^- pi^- pi^+ anti_nu_tau
1579 
1580  // tau^- --> pi^+ K^- K^- nu_tau // prepared for modes with kaons
1581  // tau^+ --> pi^- K^+ K^+ anti_nu_tau
1582  // tau^- --> K^+ K^- pi^- nu_tau // prepared for modes with kaons
1583  // tau^+ --> K^- K^+ pi^+ anti_nu_tau
1584 
1585  // tau^- --> pi^- K^0 pi^0 nu_tau // prepared for modes with kaons
1586  // tau^+ --> pi^+ K^0 pi^0 anti_nu_tau
1587 
1588 
1589  // tau^- --> pi^- K^0 K^0 nu_tau // prepared for modes with kaons
1590  // tau^+ --> pi^+ K^0 K^0 anti_nu_tau
1591 
1592 
1593  // tau^- --> K^- K^0 pi^0 nu_tau // prepared for modes with kaons
1594  // tau^+ --> K^+ K^0 pi^0 anti_nu_tau
1595 
1596  // tau^- --> K^- pi^0 pi^0 nu_tau // prepared for modes with kaons
1597  // tau^+ --> K^+ pi^0 pi^0 anti_nu_tau
1598 
1599  // 3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
1600  // 4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
1601  // 5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
1602 
1603  else if( pdgid.size()==4 &&
1604  (
1605  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, -211, 321) ) ||
1606  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 211,-321) )
1607  )
1608  ) {
1609  DEBUG( cout<<"Channel 5 : "; )
1610  channel = 14;
1611  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1612  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1613 
1614  const double AMTAU = 1.777;
1615 
1616  // IF(I.EQ.14) NAMES(I-7)=' TAU- --> K-, PI-, K+ '
1617  int MNUM = 1;
1618  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1619  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1620  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1621  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1622  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1623  float AMPLIT = 0.0;
1624  float HV[4] = { 0.0 };
1625 
1626  // For RChL currents one needs to define 3-pi sub-channel used
1627  chanopt_.JJ=1;
1628 
1629  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1630 
1631  WTamplit = AMPLIT;
1632  HH[0] = -HV[0];
1633  HH[1] = -HV[1];
1634  HH[2] = -HV[2];
1635  HH[3] = HV[3];
1636  }
1637 
1638 
1639 
1640  else if( pdgid.size()==4 &&
1641  (
1642  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 311, -211, 311 ) ) ||
1643  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 311, 211, 311 ) ) ||
1644  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 311, -211, 310 ) ) ||
1645  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 311, 211, 310 ) ) ||
1646  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 311, -211, 130 ) ) ||
1647  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 311, 211, 130 ) ) ||
1648  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 310, -211, 311 ) ) ||
1649  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 310, 211, 311 ) ) ||
1650  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 310, -211, 310 ) ) ||
1651  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 310, 211, 310 ) ) ||
1652  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 310, -211, 130 ) ) ||
1653  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 310, 211, 130 ) ) ||
1654  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 130, -211, 311 ) ) ||
1655  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 130, 211, 311 ) ) ||
1656  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 130, -211, 310 ) ) ||
1657  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 130, 211, 310 ) ) ||
1658  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 130, -211, 130 ) ) ||
1659  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 130, 211, 130 ) )
1660  )
1661  ) {
1662  DEBUG( cout<<"Channel 5 : "; )
1663  channel = 15;
1664  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1665  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1666 
1667  const double AMTAU = 1.777;
1668  // IF(I.EQ.15) NAMES(I-7)=' TAU- --> K0, PI-, K0B '
1669  int MNUM = 2;
1670  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1671  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1672  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1673  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1674  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1675  float AMPLIT = 0.0;
1676  float HV[4] = { 0.0 };
1677 
1678  // For RChL currents one needs to define 3-pi sub-channel used
1679  chanopt_.JJ=1;
1680 
1681  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1682 
1683  WTamplit = AMPLIT;
1684  HH[0] = -HV[0];
1685  HH[1] = -HV[1];
1686  HH[2] = -HV[2];
1687  HH[3] = HV[3];
1688  }
1689 
1690 
1691 
1692  else if( pdgid.size()==4 &&
1693  (
1694  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, 311, 111) ) ||
1695  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 311, 111) ) ||
1696  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, 310, 111) ) ||
1697  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 310, 111) ) ||
1698  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, -321, 130, 111) ) ||
1699  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 130, 111) )
1700  )
1701  ) {
1702  DEBUG( cout<<"Channel 5 : "; )
1703  channel = 16;
1704  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1705  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1706 
1707  const double AMTAU = 1.777;
1708  // IF(I.EQ.16) NAMES(I-7)=' TAU- --> K-, K0, PI0 '
1709  int MNUM = 3;
1710  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1711  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1712  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1713  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1714  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1715  float AMPLIT = 0.0;
1716  float HV[4] = { 0.0 };
1717 
1718  // For RChL currents one needs to define 3-pi sub-channel used
1719  chanopt_.JJ=1;
1720 
1721  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1722 
1723  WTamplit = AMPLIT;
1724  HH[0] = -HV[0];
1725  HH[1] = -HV[1];
1726  HH[2] = -HV[2];
1727  HH[3] = HV[3];
1728  }
1729 
1730 
1731 
1732  else if( pdgid.size()==4 &&
1733  (
1734  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 111, 111,-321) ) ||
1735  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 111, 111, 321) )
1736  )
1737  ) {
1738  DEBUG( cout<<"Channel 5 : "; )
1739  channel = 17;
1740  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1741  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1742 
1743  const double AMTAU = 1.777;
1744  // IF(I.EQ.17) NAMES(I-7)=' TAU- --> PI0 PI0 K- ''
1745  int MNUM = 4;
1746  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1747  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1748  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1749  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1750  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1751  float AMPLIT = 0.0;
1752  float HV[4] = { 0.0 };
1753 
1754  // For RChL currents one needs to define 3-pi sub-channel used
1755  chanopt_.JJ=1;
1756 
1757  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1758 
1759  WTamplit = AMPLIT;
1760  HH[0] = -HV[0];
1761  HH[1] = -HV[1];
1762  HH[2] = -HV[2];
1763  HH[3] = HV[3];
1764  }
1765 
1766 
1767  else if( pdgid.size()==4 &&
1768  (
1769  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-321,-211, 211) ) ||
1770  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 321, 211,-211) )
1771  )
1772  ) {
1773  DEBUG( cout<<"Channel 5 : "; )
1774  channel = 18;
1775  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1776  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1777 
1778  const double AMTAU = 1.777;
1779  // IF(I.EQ.18) NAMES(I-7)=' TAU- --> K- PI- PI+ '
1780  int MNUM = 5;
1781  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1782  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1783  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1784  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1785  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1786  float AMPLIT = 0.0;
1787  float HV[4] = { 0.0 };
1788 
1789  // For RChL currents one needs to define 3-pi sub-channel used
1790  chanopt_.JJ=1;
1791 
1792  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1793 
1794  WTamplit = AMPLIT;
1795  HH[0] = -HV[0];
1796  HH[1] = -HV[1];
1797  HH[2] = -HV[2];
1798  HH[3] = HV[3];
1799  }
1800 
1801  else if( pdgid.size()==4 &&
1802  (
1803  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 311, 111) ) ||
1804  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 311, 111) ) ||
1805  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 310, 111) ) ||
1806  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 310, 111) ) ||
1807  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211, 130, 111) ) ||
1808  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 130, 111) )
1809  )
1810  ) {
1811  DEBUG( cout<<"Channel 5 : "; )
1812  channel = 19;
1813  // MNUM=0, PT[4] tau, PN[4] neutrino, pi[4], pi[4], pi[4], AMPLIT, HH[4]
1814  // CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HH)
1815 
1816  const double AMTAU = 1.777;
1817  // IF(I.EQ.19) NAMES(I-7)=' TAU- --> PI- K0B PI0 '
1818  int MNUM = 6;
1819  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1820  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1821  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1822  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1823  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1824  float AMPLIT = 0.0;
1825  float HV[4] = { 0.0 };
1826 
1827  // For RChL currents one needs to define 3-pi sub-channel used
1828  chanopt_.JJ=1;
1829 
1830  damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &AMPLIT, HV );
1831 
1832  WTamplit = AMPLIT;
1833  HH[0] = -HV[0];
1834  HH[1] = -HV[1];
1835  HH[2] = -HV[2];
1836  HH[3] = HV[3];
1837  }
1838  // tau^- --> pi^+ pi^+ pi^0 pi^- nu_tau
1839  // tau^+ --> pi^- pi^- pi^0 pi^+ anti_nu_tau
1840  else if( pdgid.size()==5 &&
1841  (
1842  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16,-211,-211, 211, 111) ) ||
1843  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 211, 211,-211, 111) )
1844  )
1845  ) {
1846  DEBUG( cout<<"Channel 8 : "; )
1847  channel = 8;
1848 
1849  const double AMTAU = 1.777;
1850  int MNUM = 1;
1851  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1852  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1853  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1854  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1855  float PIZ [4] = { (float)tau_daughters[4].px(), (float)tau_daughters[4].py(), (float)tau_daughters[4].pz(), (float)tau_daughters[4].e() };
1856  float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1857  float AMPLIT = 0.0;
1858  float HV[4] = { 0.0 };
1859 
1860  dam4pi_( &MNUM, PT, PN, PIM1, PIM2, PIZ, PIPL, &AMPLIT, HV );
1861 
1862  WTamplit = AMPLIT;
1863  HH[0] = -HV[0];
1864  HH[1] = -HV[1];
1865  HH[2] = -HV[2];
1866  HH[3] = HV[3];
1867  }
1868  // tau^- --> pi^0 pi^0 pi^0 pi^- nu_tau
1869  // tau^+ --> pi^0 pi^0 pi^0 pi^+ anti_nu_tau
1870  else if( pdgid.size()==5 &&
1871  (
1872  ( tau_pdgid== 15 && channelMatch(tau_daughters, 16, 111, 111, 111,-211) ) ||
1873  ( tau_pdgid==-15 && channelMatch(tau_daughters,-16, 111, 111, 111, 211) )
1874  )
1875  ) {
1876  DEBUG( cout<<"Channel 9 : "; )
1877  channel = 9;
1878 
1879  const double AMTAU = 1.777;
1880  int MNUM = 2;
1881  float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1882  float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1883  float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1884  float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1885  float PIZ [4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1886  float PIPL[4] = { (float)tau_daughters[4].px(), (float)tau_daughters[4].py(), (float)tau_daughters[4].pz(), (float)tau_daughters[4].e() };
1887  float AMPLIT = 0.0;
1888  float HV[4] = { 0.0 };
1889 
1890  dam4pi_( &MNUM, PT, PN, PIM1, PIM2, PIZ, PIPL, &AMPLIT, HV );
1891 
1892  WTamplit = AMPLIT;
1893  HH[0] = -HV[0];
1894  HH[1] = -HV[1];
1895  HH[2] = -HV[2];
1896  HH[3] = HV[3];
1897  }
1898  else {
1899 
1900  DEBUG( cout<<tau_daughters.size()<<"-part ???: "; )
1901 
1902  }
1903 
1904  // Now rotate vector HH using angles phi and theta
1905  Particle HHbuf(HH[0], HH[1], HH[2], HH[3], 0);
1906 
1907  HHbuf.rotateXZ(theta);
1908  HHbuf.rotateXY(phi);
1909 
1910  HH[0] = HHbuf.px();
1911  HH[1] = HHbuf.py();
1912  HH[2] = HHbuf.pz();
1913  HH[3] = HHbuf.e ();
1914 
1915  return HH;
1916 }
1917 
1918 /*******************************************************************************
1919  Returns longitudinal polarization of the single tau (in Z/gamma* -> tau+ tau- case) averaged over
1920  incoming configurations
1921  S: invariant mass^2 of the bozon
1922  &sp_tau: first tau
1923  &sp_nu_tau: second tau (in this case it is misleading name)
1924  Hidden output: WTnonSM
1925  Hidden input: relWTnonS, nonSM2
1926 *******************************************************************************/
1927 double getLongitudinalPolarization(double S, SimpleParticle &sp_tau, SimpleParticle &sp_nu_tau)
1928 {
1929 
1930  // std::cout << "nonSM2 = " << nonSM2 << std::endl;
1931 
1932  double ReA0=0; // 0 is added to A-coupling name, because SM part may be added internally
1933  double ImA0=0;
1934  double ReB=0;
1935  double ImB=0;
1936  double ReX=0;
1937  double ImX=0;
1938  double ReY=0;
1939  double ImY=0;
1940  double ReA00=0; // 0 is added to A-coupling name, because SM part may be added internally
1941  double ImA00=0;
1942  double ReB0=0;
1943  double ImB0=0;
1944  double ReX0=0;
1945  double ImX0=0;
1946  double ReY0=0;
1947  double ImY0=0;
1948 
1949  double GSWr[7],GSWi[7];
1950  double GSWr0[7],GSWi0[7];
1951  double Amz0;
1952  double Gamz0;
1953  double sin2W0;
1954  double alphaQED;
1955  double sumak=0.0;
1956  double xsec,xsec0,xsecsum;
1957  double GAMfrac=0.0;
1958  double GAMfrac2=0.0; // ratio of gamma gamma process to the rest, in generated sample and in rewieght
1959  double A0=0.0; // anomalous couplings for gamma gamma in sample
1960  double B0=0.0;
1961  double A=0.0; // in reweight
1962  double B=0.0;
1963 
1964  if(ifGSW==0){
1965  // ExtraEWparamsGet(&Amz0,&Gamz0,&sin2W0,&alfinv,&DeltSQ,&DeltV,&Gmu,&keyGSW);
1966  // these parameters should be as of effective Born, or whatever is used in sample preparation
1967  Amz0=91.1876;
1968  Gamz0=2.4952;// 2.5;
1969  sin2W0= 0.2315200; //0.231499;
1970  alphaQED= 1.0/128.86674175;//1.0/128.9503022;
1971 
1972  for (int k=0;k<7;k++){
1973  GSWr[k]=1.0; GSWr0[k]=1.0;
1974  GSWi[k]=0.0; GSWi0[k]=0.0;
1975  }
1976 
1977  }
1978 
1979  if(ifkorch==1 && nonSM2==0) {
1980 
1981  ReA0=0; // 0 is added to A-coupling name, because SM part may be added internally
1982  ImA0=0;
1983  ReB=0;
1984  ImB=0;
1985  ReX=0;
1986  ImX=0;
1987  ReY=0;
1988  ImY=0;
1989  GAMfrac=GAMfraci; // ratio of gamma gamma process to the rest, in generated sample
1990  GAMfrac2=GAMfrac2i; // and in rewieght
1991  A0=0.; //A0i; // anomalous couplings for gamma gamma in sample
1992  B0=0.; //B0i;
1993  A=0.; //Ai; // in reweight
1994  B=0.; //Bi;
1995 
1996 
1997  }
1998 
1999  if(ifkorch==1 && nonSM2==1) {
2000 
2001  ReA0=ReAini;
2002  ImA0=ImAini;
2003  ReB=ReBini;
2004  ImB=ImBini;
2005  ReX=ReXini;
2006  ImX=ImXini;
2007  ReY=ReYini;
2008  ImY=ImYini;
2009  GAMfrac=GAMfraci; // ratio of gamma gamma process to the rest, in generated sample
2010  GAMfrac2=GAMfrac2i; // and in rewieght
2011  A0=A0i; // anomalous couplings for gamma gamma in sample
2012  B0=B0i;
2013  A=Ai; // in reweight
2014  B=Bi;
2015 
2016  }
2017 
2018  // tau+ and tau- in lab frame
2019  Particle tau_plus ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
2020  Particle tau_minus( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
2021  // P_QQ = sum of tau+ and tau- in lab frame
2022  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 );
2023 
2024  Particle P_B1(0, 0, 1, 1, 0);
2025  Particle P_B2(0, 0,-1, 1, 0);
2026 
2027  tau_plus. boostToRestFrame(P_QQ);
2028  tau_minus.boostToRestFrame(P_QQ);
2029  P_B1. boostToRestFrame(P_QQ);
2030  P_B2. boostToRestFrame(P_QQ);
2031 
2032  double costheta1 = (tau_plus.px()*P_B1.px() +tau_plus.py()*P_B1.py() +tau_plus.pz()*P_B1.pz() ) /
2033  sqrt(tau_plus.px()*tau_plus.px()+tau_plus.py()*tau_plus.py()+tau_plus.pz()*tau_plus.pz()) /
2034  sqrt(P_B1.px() *P_B1.px() +P_B1.py() *P_B1.py() +P_B1.pz() *P_B1.pz() );
2035 
2036  double costheta2 = (tau_minus.px()*P_B2.px() +tau_minus.py()*P_B2.py() +tau_minus.pz()*P_B2.pz() ) /
2037  sqrt(tau_minus.px()*tau_minus.px()+tau_minus.py()*tau_minus.py()+tau_minus.pz()*tau_minus.pz()) /
2038  sqrt(P_B2.px() *P_B2.px() +P_B2.py() *P_B2.py() +P_B2.pz() *P_B2.pz() );
2039 
2040  double sintheta1 = sqrt(1-costheta1*costheta1);
2041  double sintheta2 = sqrt(1-costheta2*costheta2);
2042 
2043  // Cosine of hard scattering
2044  double costhe = (costheta1*sintheta2 + costheta2*sintheta1) / (sintheta1 + sintheta2);
2045  double costhem= (costheta1*sintheta2 + costheta2*sintheta1) / (sintheta1 + sintheta2);
2046 
2047  // Cosine of hard scattering Mustraal frame style draft of the code.
2048  if (FrameType==1){
2049  double RRR = Tauola::randomDouble();
2050  double wtMustraal;
2051  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));
2052  if( RRR< wtMustraal) {costhe=costheta1; costhem=costheta2;} // we will need second choice for interchanged beams
2053  else {costhe=costheta2; costhem=costheta1;}
2054 
2055  // wherever below "-costhe" or "-cost" is used "-costhem" should be used instead.
2056  // calls with costhe and -costhem should be grouped.
2057 }
2058  // Invariant mass^2 of, tau+tau- pair plus photons, system!
2059  double SS = S; // other option is for tests: P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
2060 
2061  /*
2062  We need to fix sign of costhe and attribute ID.
2063  calculate x1,x2; // x1*x2 = SS/CMSENE/CMSENE; // (x1-x2)/(x1+x2)=P_QQ/CMSENE*2 in lab;
2064  calculate weight WID[]=sig(ID,SS,+/-costhe)* f(x1,+/-ID,SS) * f(x2,-/+ID,SS) ; respectively for u d c s b
2065  f(x,ID,SS,CMSENE)=x*(1-x) // for the start it will invoke library
2066  on the basis of this generate ID and set sign for costhe.
2067  then we calculate polarization averaging over incoming states.
2068  */
2069 
2070  double x1x2 = SS/CMSENE/CMSENE;
2071  double x1Mx2 = P_QQ.pz()/CMSENE*2;
2072 
2073  double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
2074  double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
2075 
2076  double WID[13];
2077  WID[0] = f(x1, 0,SS,CMSENE)*f(x2, 0,SS,CMSENE) * sigborn(0,SS, costhe);
2078  WID[1] = f(x1, 1,SS,CMSENE)*f(x2,-1,SS,CMSENE) * sigborn(1,SS, costhe);
2079  WID[2] = f(x1,-1,SS,CMSENE)*f(x2, 1,SS,CMSENE) * sigborn(1,SS,-costhem);
2080  WID[3] = f(x1, 2,SS,CMSENE)*f(x2,-2,SS,CMSENE) * sigborn(2,SS, costhe);
2081  WID[4] = f(x1,-2,SS,CMSENE)*f(x2, 2,SS,CMSENE) * sigborn(2,SS,-costhem);
2082  WID[5] = f(x1, 3,SS,CMSENE)*f(x2,-3,SS,CMSENE) * sigborn(3,SS, costhe);
2083  WID[6] = f(x1,-3,SS,CMSENE)*f(x2, 3,SS,CMSENE) * sigborn(3,SS,-costhem);
2084  WID[7] = f(x1, 4,SS,CMSENE)*f(x2,-4,SS,CMSENE) * sigborn(4,SS, costhe);
2085  WID[8] = f(x1,-4,SS,CMSENE)*f(x2, 4,SS,CMSENE) * sigborn(4,SS,-costhem);
2086  WID[9] = f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(5,SS, costhe);
2087  WID[10]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(5,SS,-costhem);
2088  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]);
2089  WID[12]= 0.0;
2090 
2091  //ERW: Genis
2092  // WID[11]= f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(22,SS, costhe);
2093  // WID[12]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(22,SS,-costhe);
2094 
2095  // REMOVE !!!
2096  //for(int i=0;i<=10;i++) {if(!(i==3)) WID[i]=0.0;}
2097  // printf(" usunac widy nadpisane [3]= %13.10f \n",WID[3] );
2098  // printf(" usunac widy nadpisane [1]= %13.10f \n",WID[1] );
2099 
2100  // end REMOVE !!!
2101 
2102  double sum = 0.0; // normalize
2103  for(int i=0;i<=12;i++) sum+=WID[i];
2104 
2105  //ERW:: Genis
2106  // for(int i=0;i<=12;i++) WID[i] = 0;
2107  //WID[1] = sigborn(11, SS, costhe);
2108  //sum = WID[1];
2109 
2110 
2111 
2112  if( sum == 0.0 )
2113  {
2114  cout << "Tauspinner::calculateWeightFromParticlesH WARNING: sum of WID[0]-WID[10] is 0. Check LHAPDF configuration" << endl;
2115  }
2116 
2117  WTnonSM=1.0;
2118  if(relWTnonSM==0) WTnonSM=sum;
2119  if(nonSM2==1)
2120  {
2121  double WID2[13];
2122  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
2123  WID2[1] = f(x1, 1,SS,CMSENE)*f(x2,-1,SS,CMSENE) * sigborn(1,SS, costhe) * plweight(1,SS, costhe);
2124  WID2[2] = f(x1,-1,SS,CMSENE)*f(x2, 1,SS,CMSENE) * sigborn(1,SS,-costhem) * plweight(1,SS,-costhem);
2125  WID2[3] = f(x1, 2,SS,CMSENE)*f(x2,-2,SS,CMSENE) * sigborn(2,SS, costhe) * plweight(2,SS, costhe);
2126  WID2[4] = f(x1,-2,SS,CMSENE)*f(x2, 2,SS,CMSENE) * sigborn(2,SS,-costhem) * plweight(2,SS,-costhem);
2127  WID2[5] = f(x1, 3,SS,CMSENE)*f(x2,-3,SS,CMSENE) * sigborn(3,SS, costhe) * plweight(3,SS, costhe);
2128  WID2[6] = f(x1,-3,SS,CMSENE)*f(x2, 3,SS,CMSENE) * sigborn(3,SS,-costhem) * plweight(3,SS,-costhem);
2129  WID2[7] = f(x1, 4,SS,CMSENE)*f(x2,-4,SS,CMSENE) * sigborn(4,SS, costhe) * plweight(4,SS, costhe);
2130  WID2[8] = f(x1,-4,SS,CMSENE)*f(x2, 4,SS,CMSENE) * sigborn(4,SS,-costhem) * plweight(4,SS,-costhem);
2131  WID2[9] = f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(5,SS, costhe) * plweight(5,SS, costhe);
2132  WID2[10]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(5,SS,-costhem) * plweight(5,SS,-costhem);
2133  // WID2[11]= f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(22,SS, costhe)* plweight(22,SS, costhe);
2134  // WID2[12]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(22,SS,-costhem)* plweight(22,SS,-costhem);
2135  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]);
2136  WID2[12]= 0.0;
2137 
2138  // REMOVE !!!
2139  // for(int i=0;i<=10;i++){ if(!(i==3)) WID2[i]=0.0;}
2140  // printf(" usunac widy2 nadpisane [2]= %13.10f \n",WID2[3] );
2141  // printf(" usunac widy2 nadpisane [1]= %13.10f \n",WID2[1] );
2142  // end REMOVE !!!
2143 
2144  double sum2 = 0.0; // normalize
2145  for(int i=0;i<=12;i++) sum2+=WID2[i];
2146  WTnonSM=sum2/sum ;
2147  // WTnonSM= WTnonSM*1.2;
2148  // printf(" WTnonSM= %13.10f \n", WTnonSM);
2149  if(relWTnonSM==0) WTnonSM=sum2;
2150 
2151  }
2152 
2153  double pol = 0.0;
2154  // double Rxx = 0.0;
2155  // double Ryy = 0.0;
2156 
2157  if(IfHiggs && nonSM2==1) { // we assume that only glue glue process contributes for Higgs
2158  double polp = plzap2(0,15,S,costhe); // 0 means incoming gluon, 15 means outgoing tau
2159  pol += (2*(1-polp)-1);
2160  return pol;
2161  }
2162  if(IfHiggs) return NAN;
2163 
2164  // case of Z/gamma
2165  for(int i=0;i<=12;i++) WID[i]/=sum;
2166 
2167 
2168  R11 = 0.0;
2169  R22 = 0.0;
2170  R12 = 0.0;
2171  R21 = 0.0;
2172 
2173  int ICC = -1;
2174 
2175  sumak=0.0; // not used at present
2176  xsecsum=0.0;
2177  if(ifkorch==1){
2178  for(int i=0;i<=3;i++){
2179  for(int j=0;j<=3;j++){
2180  RcorExt[i][j]=0.0;
2181  }
2182  }
2183  }
2184 
2185  for(int i=1;i<=12;i++)
2186  {
2187  //printf("##### petla po i = %3i \n",i );
2188  ICC = i;
2189  double cost = costhe;
2190  // first beam quark or antiquark
2191  int IROT=0;
2192  if( ICC==2 || ICC==4 || ICC==6 || ICC==8 || ICC==10|| ICC==12 ) { cost = -costhem; IROT=1;}
2193 
2194  // ID of incoming quark (up or down type)
2195  int ID = 2;
2196  if( ICC==7 || ICC==8 ) ID = 4;
2197  if( ICC==1 || ICC==2 ) ID = 1;
2198  if( ICC==5 || ICC==6 ) ID = 3;
2199  if( ICC==9 || ICC==10 ) ID = 5;
2200  if( ICC==11 || ICC==12 ) ID = 22;
2201  int tau_pdgid = 15;
2202  if(i<11){
2203  // printf("##### kic \n");
2204  double polp = plzap2(ID,tau_pdgid,S,cost);
2205  // if(ifkorch==0) pol += (2*(1-polp)-1)*WID[i];
2206  pol += (2*(1-polp)-1)*WID[i];
2207  // printf("##### klic \n");
2208  // we obtain transverse spin components of density matrix from Tauolapp pre-tabulated O(alpha) EW results
2209  //
2210  Tauolapp::TauolaParticlePair pp;
2211  // printf("##### akikur \n");
2212  // Set them to 0 in case no tables are loaded by Tauola++
2213  pp.m_R[1][1] = pp.m_R[2][2] = 0.0;
2214  pp.m_R[1][2] = pp.m_R[2][1] = 0.0;
2215  pp.m_R[0][1] = pp.m_R[3][2] = 0.0;
2216  pp.m_R[0][2] = pp.m_R[3][1] = 0.0;
2217  pp.m_R[1][0] = pp.m_R[2][3] = 0.0;
2218  pp.m_R[2][0] = pp.m_R[1][3] = 0.0;
2219  pp.recalculateRij(ID,tau_pdgid,S,cost);
2220 
2221  // S=(91.17+10.)*(91.17+10.); //100.0;
2222  // costhe=0.5;
2223  //ERW:: Genis
2224  //pp.recalculateRij(11,tau_pdgid,S,costhe);
2225  //pol = pp.m_R[0][3]/pp.m_R[0][0];
2226 
2227  // ERW:: Genis for electrons
2228  double Rcor[4][4];
2229  // double RcorExt[4][4];
2230  double Energy;
2231  double theta;
2232 
2233 
2234  Energy=sqrt(S)/2.0;
2235  theta=acos(cost);
2236 
2237 
2238  if(ifkorch==1){
2239  // ========================= BEGIN ============
2240 
2241  int IFLAV; // 0 for LEPTONS, 2 -for UP QUARK, 1 - for DOWN QUARK
2242  int channel ; // 1 for LEPTONS, 2 -for UP QUARK, 3 - for DOWN QUARK
2243  channel = 1;
2244  IFLAV=0;
2245  if (ID==2||ID==4) {channel=2 ; IFLAV=2;}
2246  if (ID==1||ID==3||ID==5){channel=3 ; IFLAV=1;}
2247 
2248 
2249 
2250  //std::cout <<
2251  // std::cout << "Amz = " << Amz0 << std::endl;
2252  // std::cout << "Gamz = " << Gamz0 << std::endl;
2253  // std::cout << "sin2W = " << sin2W0 << std::endl;
2254 
2255  if(ifGSW==1){
2256  Amz0=Amz(IFLAV);
2257  Gamz0=Gamz(IFLAV);
2258  sin2W0=sin2W(IFLAV);
2259  alphaQED=7.2973525693e-3; // check the value, should be as for alpha(s=0)= 1/137.035...
2260  for (int k=0;k<7;k++){
2261  complex<double> rezu = EWFACT(IFLAV, k, S, cost);
2262 
2263 
2264  GSWr[k]=real(rezu);GSWr0[k]=1.0;//real(rezu);
2265  GSWi[k]=imag(rezu);GSWi0[k]=0.0;//imag(rezu);
2266  // std::cout << "EWFACT ( "<< k <<")= " << GSWr[k]<<" "<< GSWi[k] << std::endl;
2267  }
2268  }
2269  GSWadapt(IFLAV,S,cost,GSWr,GSWi);
2270  dipolqq_(&iqed,&Energy,&theta,&channel,&Amz0,&Gamz0,&sin2W0,&alphaQED,&ReA00,&ImA00, &ReB0,&ImB0, &ReX0,&ImX0,&ReY0,&ImY0, GSWr0,GSWi0,Rcor);
2271  xsec0=Rcor[0][0];
2272 
2273  dipolqq_(&iqed,&Energy,&theta,&channel,&Amz0,&Gamz0,&sin2W0,&alphaQED,&ReA0,&ImA0, &ReB,&ImB, &ReX,&ImX,&ReY,&ImY, GSWr,GSWi,Rcor);
2274  xsec=Rcor[0][0];
2275  Rcor[0][0]=1.0;
2276  // TauSpinner and Korchin convention change is performed in routine dipolqq_ which provides matrix R_ij in C++ conventions
2277  // 0123 <-> txyz
2278 
2279  if(IROT==1){// needed for cost=-costhem: rotate by pi around z-axis.
2280  /*
2281  //Rotate pp.m_R around z-axis as well? May be in future.
2282  //Not used at present not necessary adjustment of frame orientation for elements which equal zero.
2283  pp.m_R[1][3]=-pp.m_R[1][3];
2284  pp.m_R[3][1]=-pp.m_R[3][1];
2285  pp.m_R[0][1]=-pp.m_R[0][1];
2286  pp.m_R[1][0]=-pp.m_R[1][0];
2287  pp.m_R[2][3]=-pp.m_R[2][3];
2288  pp.m_R[3][2]=-pp.m_R[3][2];
2289  pp.m_R[2][0]=-pp.m_R[2][0];
2290  pp.m_R[0][2]=-pp.m_R[0][2];
2291  */
2292  Rcor[1][3]=-Rcor[1][3];
2293  Rcor[3][1]=-Rcor[3][1];
2294  Rcor[0][1]=-Rcor[0][1];
2295  Rcor[1][0]=-Rcor[1][0];
2296  Rcor[2][3]=-Rcor[2][3];
2297  Rcor[3][2]=-Rcor[3][2];
2298  Rcor[2][0]=-Rcor[2][0];
2299  Rcor[0][2]=-Rcor[0][2];
2300 
2301 // note that till Jan 3 2024 these components of R were not used in any applications.
2302  }
2303  // pol += WID[i]*pp.m_R[0][3]/pp.m_R[0][0];
2304  // std::cout << "przed daje pol("<<i<<") = " << pol << std::endl;
2305  // INFO:
2306  // pol += WID[i]*Rcor[2][3]; not needed there is:
2307  // if( ifkorch==1) pol = RcorExt[3][2];
2308  // in calculateWeightFromParticlesH(){...
2309  //
2310  // pp.m_R[0][3]/pp.m_R[0][0];
2311  // 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);
2312  // sumak +=WID[i]*pp.m_R[0][0]; --> replaced by xsecsum
2313  }
2314 
2315  // std::cout << "pol("<<i<<") = " << pol << std::endl;
2316  //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;
2317  //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;
2318  //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;
2319  //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;
2320 
2321  // ========================= END (korch==1) ============
2322  // printf("############# pp.m_R[1][1]= %13.10f \n", pp.m_R[1][1]);
2323  // These calculations are valid for Standard Model only
2324 
2325  //ERW 6.04.2025: pp.m_R jest w konwencji pochodzacej z fortranu.
2326  // to wypelnianie tylko dla (ifkorch==0)??
2327  R11 += WID[i]*pp.m_R[1][1];
2328  R22 += WID[i]*pp.m_R[2][2]; // note signs in line WT =
2329  R12 += WID[i]*pp.m_R[1][2]; // 1.0+sign*HHp[2]*HHm[2] ....- RzYY*R22 ....- RzYX*R21*HHp[1]*HHm[0];
2330  R21 += WID[i]*pp.m_R[2][1]; // sign operation with F=1/-1 like for RcorExt is not necessary.
2331 
2332  if(ifkorch==1){
2333  for(int i3=0;i3<=3;i3++){
2334  for(int j=0;j<=3;j++){
2335  int i4=i3-1;
2336  int j4=j-1;
2337  if(i4==-1) i4=3;
2338  if(j4==-1) j4=3;
2339  RcorExt[i4][j4]+=WID[i]*Rcor[i3][j]; //we move back to fortran convention xyzt because of f77-tauola polarimetric vectors
2340  }
2341  }
2342 
2343  xsecsum+=xsec/xsec0*WID[i];
2344 
2345  }
2346  }
2347 
2348  else{
2349  //ERW 6.04.2025: petla else tylko dla (ifkorch == 1)
2350  double Rcor[4][4];
2351  double polp = 0.5; // gamagama bring no polarization
2352  pol += (2*(1-polp)-1)*WID[i];
2353  // iqed QED part of anomalous switched on/off.
2354 
2355  double E;
2356  double theta;
2357  E=sqrt(S)/2.0;
2358  theta=acos(cost);
2359  dipolgamma_(&iqed, &E, &theta, &A0, &B0, Rcor);
2360  /*
2361  printf(" -------------------------------------------------\n");
2362  //ERW: rotation should be moved to dipolgamma OR NOT??
2363  printf("SM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2364  Rcor[0][0],Rcor[0][1],Rcor[0][2],Rcor[0][3]);
2365  printf("SM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2366  Rcor[1][0],Rcor[1][1],Rcor[1][2],Rcor[1][3]);
2367  printf("SM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2368  Rcor[2][0],Rcor[2][1],Rcor[2][2],Rcor[2][3]);
2369  printf("SM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2370  Rcor[3][0],Rcor[3][1],Rcor[3][2],Rcor[3][3]);
2371  */
2372  xsec0=Rcor[0][0];
2373 
2374 
2375  dipolgamma_(&iqed, &E, &theta, &A, &B, Rcor);
2376  /*
2377  printf(" -----------------------------------------------\n");
2378  //ERW: rotation should be moved to dipolgamma OR NOT??
2379  printf("BSM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2380  Rcor[0][0],Rcor[0][1],Rcor[0][2],Rcor[0][3]);
2381  printf("BSM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2382  Rcor[1][0],Rcor[1][1],Rcor[1][2],Rcor[1][3]);
2383  printf("BSM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2384  Rcor[2][0],Rcor[2][1],Rcor[2][2],Rcor[2][3]);
2385  printf("BSM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2386  Rcor[3][0],Rcor[3][1],Rcor[3][2],Rcor[3][3]);
2387  */
2388  xsec=Rcor[0][0];
2389  Rcor[0][0]=1.0;
2390  // TauSpinner and Korchin convention change is performed in routine dipolgamma_ which provides matrix R_ij in C++ conventions
2391  // 0123 <-> txyz
2392 
2393 
2394  /*
2395  // remnant print from the code infancy
2396  DEBUG
2397  (
2398 
2399  printf("dipol gamma przed r11 = %13.10f \n",R11 );
2400  printf("dipol gamma przed r22 = %13.10f \n",R22 );
2401  printf("dipol gamma przed r12 = %13.10f \n",R12 );
2402  printf("dipol gamma przed r21 = %13.10f \n",R21 );
2403  )
2404  */
2405  // ERW 6.04.2025: TO JEST ZLE??? i niepotrzebne, gamma gamma uzywane tylko jezeli (ifkorch==1)
2406  // i wowczas powinno byc wypelniane macierza wg konwecji RcorExt i bez normalizacji bo juz
2407  // zrobiona na poziomie fortranowskim
2408  R11 += WID[i]*Rcor[1][1]/Rcor[0][0]; // Historical --> division not necessary? Rcor[0][0]=1.0
2409  R22 += WID[i]*Rcor[2][2]/Rcor[0][0]; // note signs in line WT =
2410  R12 += WID[i]*Rcor[1][2]/Rcor[0][0]; // 1.0+sign*HHp[2]*HHm[2] ....- RzYY*R22 ....- RzYX*R21*HHp[1]*HHm[0];
2411  R21 += WID[i]*Rcor[2][1]/Rcor[0][0]; // sign operation with F=1/-1 like for RcorExt is not necessary.
2412 
2413  if(ifkorch==1){
2414  for(int i3=0;i3<=3;i3++){
2415  for(int j=0;j<=3;j++){
2416  int i4=i3-1;
2417  int j4=j-1;
2418  if(i4==-1) i4=3;
2419  if(j4==-1) j4=3;
2420  RcorExt[i4][j4]+=WID[i]*Rcor[i3][j]; // we move back to fortran convention xyzt because of f77-tauola polarimetric vectors
2421  }
2422  }
2423  xsecsum+=xsec/xsec0*WID[i];
2424  /*
2425  printf(" -----------------------------------------------\n");
2426 
2427  //ERW: rotation should be moved to dipolgamma OR NOT??
2428  printf("BSM R (0, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2429  RcorExt[0][0],RcorExt[0][1],RcorExt[0][2],RcorExt[0][3]);
2430  printf("BSM R (1, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2431  RcorExt[1][0],RcorExt[1][1],RcorExt[1][2],RcorExt[1][3]);
2432  printf("BSM R (2, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2433  RcorExt[2][0],RcorExt[2][1],RcorExt[2][2],RcorExt[2][3]);
2434  printf("BSM R (3, 0 1 2 3)= %13.10f %13.10f %13.10f %13.10f \n",
2435  RcorExt[3][0],RcorExt[3][1],RcorExt[3][2],RcorExt[3][3]);
2436 
2437  printf(" -----------------------------------------------\n");
2438  */
2439  }
2440 
2441  /*
2442  // remnant print from the code infancy
2443  DEBUG
2444  (
2445  printf("dipol gamma r11 = %13.10f \n",R11 );
2446  printf("dipol gamma r22 = %13.10f \n",R22 );
2447  printf("dipol gamma r12 = %13.10f \n",R12 );
2448  printf("dipol gamma r21 = %13.10f \n",R21 );
2449  )
2450  */
2451 
2452  }
2453  }
2454  if(ifkorch==1){ // Act on RcorExt after sum of all partonic contributions.
2455  for (int k0=0;k0<4;k0++){
2456  for (int k1=0;k1<4;k1++){
2457  double F=1.0;
2458  if(k0==1) F=-1.0; // possibly combination of two operation:
2459  // 1. sign issue equivalent to overall sign in front of HHP
2460  // 2. HHP pi-angle rotation around Y-axis
2461  // difference between KORLB-KKMC conventions of z axis and this of TauSpinner.
2462  // Checks for that are pending.
2463  // WARNING; ugly inconsisten fix, it should be through HHP rotation and sign adjustment not on R_ij.
2464  // ERW 6.04.25: absorb change in sign here to have RcorExt[k0][k1] ready ZBW but this may be about sign for HHp
2465 
2466  RcorExt[k0][k1]=RcorExt[k0][k1]*F;
2467  }
2468  }
2469  }
2470 
2471  //printf("BSM R (0,1,2,3)= %13.10f %13.10f %13.10f %13.10f %13.10f %13.10f %13.10f \n",
2472  // RcorExt[0][0],RcorExt[1][1],RcorExt[2][2],RcorExt[3][3], RcorExt[1][2],RcorExt[1][3],RcorExt[2][3]);
2473 
2474 
2475  if(ifkorch==1){
2476 
2477  if(relWTnonSM!=0) WTnonSM=xsecsum;
2478  if(relWTnonSM==0) {WTnonSM=xsecsum;
2479  std::cout << "input not supported: relWTnonSM=" <<relWTnonSM<< " ifkorch= "<<ifkorch << std::endl;
2480  std::cout << "X-sects nonSM2="<<nonSM2<<" sum= " << sum << " WTnonSM-korch= "<< xsecsum <<std::endl;
2481  }
2482  }
2483  // std::cout << "X-sects-korch (nonSM2="<<nonSM2<<")="<< WTnonSM <<std::endl;
2484  // Calculations are prepared only for pp collision.
2485  // Otherwise pol = 0.0
2486  if(!Ipp) pol=0.0;
2487 
2488  return pol;
2489 }
2490 
2491 /*******************************************************************************
2492  Check if pdg's of particles in the vector<Particle>&particles match the
2493  of (p1,p2,p3,p4,p5,p6)
2494 
2495  Returns true if 'particles' contain all of the listed pdgid-s.
2496  If it does - vector<Particle>&particles will be sorted in the order
2497  as listed pdgid-s.
2498 
2499  It is done so the order of particles is the same as the order used by
2500  TAUOLA Fortran routines.
2501 *******************************************************************************/
2502 bool channelMatch(vector<Particle> &particles, int p1, int p2, int p3, int p4, int p5, int p6)
2503 {
2504  // Copy pdgid-s of all particles
2505  vector<int> list;
2506 
2507  for(unsigned int i=0;i<particles.size();i++) list.push_back(particles[i].pdgid());
2508 
2509  // Create array out of pdgid-s
2510  int p[6] = { p1, p2, p3, p4, p5, p6 };
2511 
2512  // 1) Check if 'particles' contain all pdgid-s on the list 'p'
2513 
2514  for(int i=0;i<6;i++)
2515  {
2516  // if the pdgid is zero - finish
2517  if(p[i]==0) break;
2518 
2519  bool found = false;
2520 
2521  for(unsigned int j=0;j<list.size(); j++)
2522  {
2523  // if pdgid is found - erese it from the list and search for the next one
2524  if(list[j]==p[i])
2525  {
2526  found = true;
2527  list.erase(list.begin()+j);
2528  break;
2529  }
2530  }
2531 
2532  if(!found) return false;
2533  }
2534 
2535  // if there are more particles on the list - there is no match
2536  if(list.size()!=0) return false;
2537 
2538 
2539  // 2) Rearrange particles to match the order of pdgid-s listed in array 'p'
2540 
2541  vector<Particle> newList;
2542 
2543  for(int i=0;i<6;i++)
2544  {
2545  // if the pdgid is zero - finish
2546  if(p[i]==0) break;
2547 
2548  for(unsigned int j=0;j<particles.size(); j++)
2549  {
2550  // if pdgid is found - copy it to new list and erese from the old one
2551  if(particles[j].pdgid()==p[i])
2552  {
2553  newList.push_back(particles[j]);
2554  particles.erase(particles.begin()+j);
2555  break;
2556  }
2557  }
2558  }
2559 
2560  particles = newList;
2561 
2562  return true;
2563 }
2564 
2565 /*******************************************************************************
2566  Prints out two vertices:
2567  W -> tau, nu_tau
2568  tau -> tau_daughters
2569 *******************************************************************************/
2570 void print(Particle &W, Particle &nu_tau, Particle &tau, vector<Particle> &tau_daughters) {
2571 
2572  nu_tau.print();
2573  tau .print();
2574 
2575  double px=nu_tau.px()+tau.px();
2576  double py=nu_tau.py()+tau.py();
2577  double pz=nu_tau.pz()+tau.pz();
2578  double e =nu_tau.e ()+tau.e ();
2579 
2580  // Print out sum of tau and nu_tau and also W momentum for comparison
2581  cout<<"--------------------------------------------------------------------------------------------------------"<<endl;
2582  Particle sum1(px,py,pz,e,0);
2583  sum1.print();
2584  W .print();
2585 
2586  cout<<endl;
2587 
2588  // Print out tau daughters
2589  for(unsigned int i=0; i<tau_daughters.size();i++) tau_daughters[i].print();
2590 
2591  // Print out sum of tau decay products, and also tau momentum for comparison
2592  cout<<"--------------------------------------------------------------------------------------------------------"<<endl;
2593  Particle *sum2 = vector_sum(tau_daughters);
2594  sum2->print();
2595  tau.print();
2596  cout<<endl;
2597 
2598  delete sum2;
2599 }
2600 
2601 /*******************************************************************************
2602  Sums all 4-vectors of the particles on the list
2603 *******************************************************************************/
2604 Particle *vector_sum(vector<Particle> &x) {
2605 
2606  double px=0.0,py=0.0,pz=0.0,e=0.0;
2607 
2608  for(unsigned int i=0; i<x.size();i++)
2609  {
2610  px+=x[i].px();
2611  py+=x[i].py();
2612  pz+=x[i].pz();
2613  e +=x[i].e();
2614  }
2615 
2616  Particle *sum = new Particle(px,py,pz,e,0);
2617  return sum;
2618 }
2619 
2620 } // namespace TauSpinner
2621 
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:727
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:2502
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:633
TauSpinner::getZgamMultipliersTR
void getZgamMultipliersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
Definition: tau_reweight_lib.cxx:554
TauSpinner::getWtamplitP
double getWtamplitP()
Definition: tau_reweight_lib.cxx:606
TauSpinner::getZgamParametersTR
void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
Definition: tau_reweight_lib.cxx:544
TauSpinner::getLongitudinalPolarization
double getLongitudinalPolarization(double, SimpleParticle &, SimpleParticle &)
Definition: tau_reweight_lib.cxx:1927
TauSpinner::getHiggsParameters
void getHiggsParameters(double *mass, double *width, double *normalization)
Definition: tau_reweight_lib.cxx:567
TauSpinner::setZgamMultipliersTR
void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx)
Definition: tau_reweight_lib.cxx:533
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::getZgamParametersL
void getZgamParametersL(double &Rzx, double &Rzy, double &Rzz, double &Rtx, double &Rty, double &Rtz)
Definition: tau_reweight_lib.cxx:520
TauSpinner::vector_sum
Particle * vector_sum(vector< Particle > &x)
Definition: tau_reweight_lib.cxx:2604
TauSpinner::setNonSMkey
void setNonSMkey(int key)
Definition: tau_reweight_lib.cxx:591
TauSpinner::getWtNonSM
double getWtNonSM()
Definition: tau_reweight_lib.cxx:599
TauSpinner::calculateHH
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
Definition: tau_reweight_lib.cxx:1177
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:2570
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:1080
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:614
TauSpinner::getWtamplitM
double getWtamplitM()
Definition: tau_reweight_lib.cxx:607
Generated by doxygen 1.9.1