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