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