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