11 #include "EvtGenBase/EvtParticle.hh"
12 #include "EvtGenBase/EvtPDL.hh"
13 #include "EvtGenBase/EvtScalarParticle.hh"
14 #include "EvtGenBase/EvtDiracSpinor.hh"
15 #include "EvtGenBase/EvtId.hh"
16 #include "EvtGenBase/EvtAmp.hh"
18 #include "generators/evtgen/models/EvtBSemiTauonicAmplitude.h"
19 #include "generators/evtgen/models/EvtBSemiTauonicHelicityAmplitudeCalculator.h"
33 EvtVector4R p4Boosted = boostTo(p->getP4(), p4boost,
true);
34 const double theta = acos(p4Boosted.get(3) / p4Boosted.d3mag());
35 const double phi = atan2(p4Boosted.get(2), p4Boosted.get(1));
38 EvtDiracSpinor spRest[2] = {p->sp(0), p->sp(1)};
40 sp[0] = boostTo(spRest[0], p4boost);
41 sp[1] = boostTo(spRest[1], p4boost);
43 EvtDiracSpinor spplus;
44 EvtDiracSpinor spminus;
48 if (EvtPDL::getStdHep(p->getId()) > 0) {
49 spplus.set(1.0, 0.0, 0.0, 0.0);
50 spminus.set(0.0, 1.0, 0.0, 0.0);
51 norm = sqrt(real(sp[0].get_spinor(0) * sp[0].get_spinor(0) + sp[0].get_spinor(1) * sp[0].get_spinor(1)));
53 spplus.set(0.0, 0.0, 0.0, 1.0);
54 spminus.set(0.0, 0.0, 1.0, 0.0);
55 norm = sqrt(real(sp[0].get_spinor(2) * sp[0].get_spinor(2) + sp[0].get_spinor(3) * sp[0].get_spinor(3)));
58 spplus.applyRotateEuler(phi, theta, -phi);
59 spminus.applyRotateEuler(phi, theta, -phi);
64 for (
int i = 0; i < 2; i++) {
65 if (EvtPDL::getStdHep(p->getId()) > 0) {
66 R.set(0, i, (spplus * sp[i]) / norm);
67 R.set(1, i, (spminus * sp[i]) / norm);
69 R.set(0, i, (sp[i]*spplus) / norm);
70 R.set(1, i, (sp[i]*spminus) / norm);
80 EvtId lepton, EvtId nudaug,
81 EvtBSemiTauonicHelicityAmplitudeCalculator* CalcHelAmp)
94 EvtScalarParticle* scalar_part;
95 EvtParticle* root_part;
97 scalar_part =
new EvtScalarParticle;
100 scalar_part->noLifeTime();
104 p_init.set(EvtPDL::getMass(parent), 0.0, 0.0, 0.0);
105 scalar_part->init(parent, p_init);
107 root_part =
static_cast<EvtParticle*
>(scalar_part);
109 root_part->setDiagonalSpinDensity();
111 EvtParticle* daughter, *lep, *trino;
117 listdaug[1] = lepton;
118 listdaug[2] = nudaug;
120 amp.init(parent, 3, listdaug);
122 root_part->makeDaughters(3, listdaug);
123 daughter = root_part->getDaug(0);
124 lep = root_part->getDaug(1);
125 trino = root_part->getDaug(2);
128 daughter->noLifeTime();
136 rho.setDiag(root_part->getSpinStates());
140 double m = root_part->mass();
142 EvtVector4R p4meson, p4lepton, p4nu, p4w;
144 double q2, elepton, plepton;
146 double erho, prho, costl;
148 double maxfoundprob = 0.0;
152 for (massiter = 0; massiter < 3; massiter++) {
154 mass[0] = EvtPDL::getMeanMass(meson);
155 mass[1] = EvtPDL::getMeanMass(lepton);
156 mass[2] = EvtPDL::getMeanMass(nudaug);
158 mass[0] = EvtPDL::getMinMass(meson);
161 mass[0] = EvtPDL::getMaxMass(meson);
162 if ((mass[0] + mass[1] + mass[2]) > m) mass[0] = m - mass[1] - mass[2] - 0.00001;
165 double q2min = mass[1] * mass[1];
166 double q2max = (m - mass[0]) * (m - mass[0]);
167 double dq2 = (q2max - q2min) / 25;
172 for (i = 0; i < 25; i++) {
173 q2 = q2min + (i + 0.5) * dq2;
175 erho = (m * m + mass[0] * mass[0] - q2) / (2.0 * m);
177 prho = sqrt(erho * erho - mass[0] * mass[0]);
179 p4meson.set(erho, 0.0, 0.0, -1.0 * prho);
180 p4w.set(m - erho, 0.0, 0.0, prho);
186 elepton = (q2 + mass[1] * mass[1]) / (2.0 * sqrt(q2));
187 plepton = sqrt(elepton * elepton - mass[1] * mass[1]);
192 for (j = 0; j < 3; j++) {
194 costl = 0.99 * (j - 1.0);
198 p4lepton.set(elepton, 0.0,
199 plepton * sqrt(1.0 - costl * costl), plepton * costl);
200 p4nu.set(plepton, 0.0,
201 -1.0 * plepton * sqrt(1.0 - costl * costl), -1.0 * plepton * costl);
203 EvtVector4R boost((m - erho), 0.0, 0.0, 1.0 * prho);
204 p4lepton = boostTo(p4lepton, boost);
205 p4nu = boostTo(p4nu, boost);
209 daughter->init(meson, p4meson);
210 lep->init(lepton, p4lepton);
211 trino->init(nudaug, p4nu);
213 CalcAmp(root_part, amp, CalcHelAmp);
219 prob = rho.normalizedProb(amp.getSpinDensity());
227 double a = probctl[1];
228 double b = 0.5 * (probctl[2] - probctl[0]);
229 double c = 0.5 * (probctl[2] + probctl[0]) - probctl[1];
232 if (probctl[1] > prob) prob = probctl[1];
233 if (probctl[2] > prob) prob = probctl[2];
235 if (fabs(c) > 1e-20) {
236 double ctlx = -0.5 * b / c;
237 if (fabs(ctlx) < 1.0) {
238 double probtmp = a + b * ctlx + c * ctlx * ctlx;
239 if (probtmp > prob) prob = probtmp;
244 if (prob > maxfoundprob) {
249 if (EvtPDL::getWidth(meson) <= 0.0) {
255 root_part->deleteTree();