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