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