Belle II Software  release-05-01-25
EvtBSemiTauonicAmplitude.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Koji Hara *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
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"
17 
18 #include "generators/evtgen/models/EvtBSemiTauonicAmplitude.h"
19 #include "generators/evtgen/models/EvtBSemiTauonicHelicityAmplitudeCalculator.h"
20 
21 namespace Belle2 {
27  using std::endl;
28 
29  EvtSpinDensity EvtBSemiTauonicAmplitude::RotateToHelicityBasisInBoostedFrame(const EvtParticle* p,
30  EvtVector4R p4boost)
31  {
32  // theta and phi of p momentum in p4boost rest frame
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));
36 
37  // here p must be EvtDiracParticle (if not EvtGen will abort)
38  EvtDiracSpinor spRest[2] = {p->sp(0), p->sp(1)};
39  EvtDiracSpinor sp[2];
40  sp[0] = boostTo(spRest[0], p4boost);
41  sp[1] = boostTo(spRest[1], p4boost);
42 
43  EvtDiracSpinor spplus;
44  EvtDiracSpinor spminus;
45 
46  double norm;
47 
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)));
52  } else {
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)));
56  }
57 
58  spplus.applyRotateEuler(phi, theta, -phi);
59  spminus.applyRotateEuler(phi, theta, -phi);
60 
61  EvtSpinDensity R;
62  R.setDim(2);
63 
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);
68  } else {
69  R.set(0, i, (sp[i]*spplus) / norm);
70  R.set(1, i, (sp[i]*spminus) / norm);
71  }
72  }
73 
74  return R;
75 
76  }
77 
78 // copied from EvtSemileptonicAmp.cpp and modified for BSemiTauonic
79  double EvtBSemiTauonicAmplitude::CalcMaxProb(EvtId parent, EvtId meson,
80  EvtId lepton, EvtId nudaug,
81  EvtBSemiTauonicHelicityAmplitudeCalculator* CalcHelAmp)
82  {
83  //This routine takes the arguements parent, meson, and lepton
84  //number, and a form factor model, and returns a maximum
85  //probability for this semileptonic form factor model. A
86  //brute force method is used. The 2D cos theta lepton and
87  //q2 phase space is probed.
88 
89  //Start by declaring a particle at rest.
90 
91  //It only makes sense to have a scalar parent. For now.
92  //This should be generalized later.
93 
94  EvtScalarParticle* scalar_part;
95  EvtParticle* root_part;
96 
97  scalar_part = new EvtScalarParticle;
98 
99  //cludge to avoid generating random numbers!
100  scalar_part->noLifeTime();
101 
102  EvtVector4R p_init;
103 
104  p_init.set(EvtPDL::getMass(parent), 0.0, 0.0, 0.0);
105  scalar_part->init(parent, p_init);
106  //root_part = (EvtParticle*)scalar_part;
107  root_part = static_cast<EvtParticle*>(scalar_part);
108 
109  root_part->setDiagonalSpinDensity();
110 
111  EvtParticle* daughter, *lep, *trino;
112 
113  EvtAmp amp;
114 
115  EvtId listdaug[3];
116  listdaug[0] = meson;
117  listdaug[1] = lepton;
118  listdaug[2] = nudaug;
119 
120  amp.init(parent, 3, listdaug);
121 
122  root_part->makeDaughters(3, listdaug);
123  daughter = root_part->getDaug(0);
124  lep = root_part->getDaug(1);
125  trino = root_part->getDaug(2);
126 
127  //cludge to avoid generating random numbers!
128  daughter->noLifeTime();
129  lep->noLifeTime();
130  trino->noLifeTime();
131 
132 
133  //Initial particle is unpolarized, well it is a scalar so it is
134  //trivial
135  EvtSpinDensity rho;
136  rho.setDiag(root_part->getSpinStates());
137 
138  double mass[3];
139 
140  double m = root_part->mass();
141 
142  EvtVector4R p4meson, p4lepton, p4nu, p4w;
143 
144  double q2, elepton, plepton;
145  int i, j;
146  double erho, prho, costl;
147 
148  double maxfoundprob = 0.0;
149  double prob;
150  int massiter;
151 
152  for (massiter = 0; massiter < 3; massiter++) {
153 
154  mass[0] = EvtPDL::getMeanMass(meson);
155  mass[1] = EvtPDL::getMeanMass(lepton);
156  mass[2] = EvtPDL::getMeanMass(nudaug);
157  if (massiter == 1) {
158  mass[0] = EvtPDL::getMinMass(meson);
159  }
160  if (massiter == 2) {
161  mass[0] = EvtPDL::getMaxMass(meson);
162  if ((mass[0] + mass[1] + mass[2]) > m) mass[0] = m - mass[1] - mass[2] - 0.00001;
163  }
164 
165  double q2min = mass[1] * mass[1]; // limit to minimum=lepton mass square
166  double q2max = (m - mass[0]) * (m - mass[0]);
167  double dq2 = (q2max - q2min) / 25;
168 // std::cout<<"m: "<<m<<" mass[0]: "<<mass[0]<<" q2min: "<<q2min<<" q2max: "<<q2max<<std::endl;
169 
170  //loop over q2
171 
172  for (i = 0; i < 25; i++) {
173  q2 = q2min + (i + 0.5) * dq2; // <-- !! not start from unphysical q2=0 !!
174 
175  erho = (m * m + mass[0] * mass[0] - q2) / (2.0 * m);
176 
177  prho = sqrt(erho * erho - mass[0] * mass[0]);
178 
179  p4meson.set(erho, 0.0, 0.0, -1.0 * prho);
180  p4w.set(m - erho, 0.0, 0.0, prho);
181 
182 // std::cout<<"q2: "<<q2<<std::endl;
183 // std::cout<<"p4meson: "<<p4meson<<std::endl;
184 
185  //This is in the W rest frame
186  elepton = (q2 + mass[1] * mass[1]) / (2.0 * sqrt(q2));
187  plepton = sqrt(elepton * elepton - mass[1] * mass[1]);
188 // std::cout<<"elepton: "<<elepton<<" plepton: "<<plepton<<std::endl;
189 
190  double probctl[3];
191 
192  for (j = 0; j < 3; j++) {
193 
194  costl = 0.99 * (j - 1.0);
195 
196  //These are in the W rest frame. Need to boost out into
197  //the B frame.
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);
202 
203  EvtVector4R boost((m - erho), 0.0, 0.0, 1.0 * prho);
204  p4lepton = boostTo(p4lepton, boost);
205  p4nu = boostTo(p4nu, boost);
206 
207  //Now initialize the daughters...
208 
209  daughter->init(meson, p4meson);
210  lep->init(lepton, p4lepton);
211  trino->init(nudaug, p4nu);
212 
213  CalcAmp(root_part, amp, CalcHelAmp);
214 
215  //Now find the probability at this q2 and cos theta lepton point
216  //and compare to maxfoundprob.
217 
218  //Do a little magic to get the probability!!
219  prob = rho.normalizedProb(amp.getSpinDensity());
220 
221  probctl[j] = prob;
222  }
223 
224  //probclt contains prob at ctl=-1,0,1.
225  //prob=a+b*ctl+c*ctl^2
226 
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];
230 
231  prob = probctl[0];
232  if (probctl[1] > prob) prob = probctl[1];
233  if (probctl[2] > prob) prob = probctl[2];
234 
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;
240  }
241 
242  }
243 
244  if (prob > maxfoundprob) {
245  maxfoundprob = prob;
246  }
247 
248  }
249  if (EvtPDL::getWidth(meson) <= 0.0) {
250  //if the particle is narrow dont bother with changing the mass.
251  massiter = 4;
252  }
253 
254  }
255  root_part->deleteTree();
256 
257  maxfoundprob *= 1.1;
258  return maxfoundprob;
259 
260  }
261 
263 } // Belle 2 Namespace
Belle2::EvtBSemiTauonicAmplitude::CalcMaxProb
double CalcMaxProb(EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator)
The function calculates the maximum probability.
Definition: EvtBSemiTauonicAmplitude.cc:87
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::EvtBSemiTauonicAmplitude::CalcAmp
virtual void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator)=0
The function calculates the spin dependent amplitude.
Belle2::EvtBSemiTauonicAmplitude::RotateToHelicityBasisInBoostedFrame
EvtSpinDensity RotateToHelicityBasisInBoostedFrame(const EvtParticle *p, EvtVector4R p4boost)
The function calculates the rotation matrix to convert the spin basis to the helicity basis in the bo...
Definition: EvtBSemiTauonicAmplitude.cc:37