C++ Interface to Tauola
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"
5using namespace Tauolapp;
6
7namespace TauSpinner {
8
9/***** GLOBAL VARIABLES AND THEIR PRE-INITIALISATION *****/
10
11double CMSENE = 7000.0;// Center of mass system energy (used by PDF's)
12bool Ipp = true; // pp collisions
13int Ipol = 1; // Is the sample in use polarized? (relevant for Z/gamma case only)
14int nonSM2 = 0; // Turn on/off nonSM calculations
15int nonSMN = 0; // Turn on, if calculations of nonSM weight, is to be used to modify shapes only
16int relWTnonSM = 1; // 1: relWTnonSM is relative to SM; 0: absolute
17double WTnonSM=1.0; // nonSM weight
18double 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()
21bool IfHiggs=false; // flag used in sigborn()
22double WTamplit=1; // AMPLIT weight for the decay last invoked
23double WTamplitP=1; // AMPLIT weight for tau+ decay
24double WTamplitM=1; // AMPLIT weight for tau- decay
25
26// Higgs parameters
27int IfHsimple=0; // switch for simple Higgs (just Breit-Wigner)
28double XMH = 125.0; // mass (GeV)
29double XGH = 1.0; // width, should be detector width, analysis dependent.
30double Xnorm = 0.15; // normalization of Higgs Born cross-section, at hoc
31
32// Transverse components of Higgs spin density matrix
33double RXX = 0.0; //-1.0;
34double RYY = 0.0; // 1.0;
35double RXY = 0.0;
36double RYX = 0.0;
37
38// Coefficients for transverse components of Z/gamma spin density matrix
39double RzXX = 0.0; //-1.0;
40double RzYY = 0.0; // 1.0;
41double RzXY = 0.0;
42double RzYX = 0.0;
43
44// Values of transverse components of Z/gamma spin density matrix calculated inside getLongitudinalPolarization
45double R11 = 0.0;
46double R22 = 0.0;
47double R12 = 0.0; // for future use
48double R21 = 0.0; // for future use
49
50/***** END: GLOBAL VARIABLES AND THEIR PRE-INITIALISATION *****/
51
52double 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
81double 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*******************************************************************************/
133void 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*******************************************************************************/
181void 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*******************************************************************************/
201void 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*******************************************************************************/
214void 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*******************************************************************************/
225void 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*******************************************************************************/
236void 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*******************************************************************************/
252void setSpinOfSample(int _Ipol)
253{
254 Ipol = _Ipol;
255}
256
257/*******************************************************************************
258 Turn nonSM calculation of Born cross-section on/off
259*******************************************************************************/
260void setNonSMkey(int _key)
261{
262 nonSM2 = _key;
263}
264
265/*******************************************************************************
266 Get nonSM weight
267*******************************************************************************/
269{
270 return WTnonSM;
271}
272/*******************************************************************************
273 Get weights for tau+ tau- decay matrix elements
274*******************************************************************************/
275double getWtamplitP(){return WTamplitP;}
276double 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*******************************************************************************/
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*******************************************************************************/
302double 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*******************************************************************************/
396double 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*******************************************************************************/
667void 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*******************************************************************************/
764double* 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*******************************************************************************/
1514double 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*******************************************************************************/
1681bool 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*******************************************************************************/
1749void 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*******************************************************************************/
1783Particle *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
void recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta)
double plweight(int ide, double svar, double costhe)
Definition: nonSM.cxx:175
double plzap2(int ide, int idf, double svar, double costhe)
Definition: nonSM.cxx:115
double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector< SimpleParticle > &sp_tau1_daughters, vector< SimpleParticle > &sp_tau2_daughters)
bool channelMatch(vector< Particle > &particles, int p1, int p2=0, int p3=0, int p4=0, int p5=0, int p6=0)
double disth_(double *SVAR, double *COSTHE, int *TA, int *TB)
double calculateWeightFromParticlesWorHpn(SimpleParticle &W, SimpleParticle &tau, SimpleParticle &nu_tau, vector< SimpleParticle > &tau_daughters)
void setSpinOfSample(bool _Ipol)
double getWtamplitP()
void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
double getLongitudinalPolarization(double, SimpleParticle &, SimpleParticle &)
void getHiggsParameters(double *mass, double *width, double *normalization)
void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx)
Particle * vector_sum(vector< Particle > &x)
void setNonSMkey(int key)
double getWtNonSM()
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
void setHiggsParametersTR(double Rxx, double Ryy, double Rxy, double Ryx)
void setRelWTnonSM(int _relWTnonSM)
void print(Particle &W, Particle &nu_tau, Particle &tau, vector< Particle > &tau_daughters)
void initialize_spinner(bool _Ipp, int _Ipol, int _nonSM2, int _nonSMN, double _CMSENE)
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
void setHiggsParameters(int jak, double mass, double width, double normalization)
double getTauSpin()
double getWtamplitM()