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