Belle II Software  release-05-01-25
BSTD_calc_dgam.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 "generators/evtgen/EvtBSemiTauonicDecayRateCalculator.h"
12 #include "generators/evtgen/models/EvtBSemiTauonicHelicityAmplitudeCalculator.h"
13 
14 #include <cstdlib>
15 #include <iostream>
16 #include <fstream>
17 
18 #include <unistd.h>
19 
20 using namespace Belle2;
21 
22 void usage()
23 {
24  std::cout << "Usage: BSTD_calc_dgam [-hr] m_B m_D* m_D m_tau m_lep rho1^2 rhoA1^2 R1(1) R2(1) aS1 aR3 m_b m_c" << std::endl;
25  std::cout << "Calculate differential decay rate for each helicity state, " << std::endl;
26  std::cout << "the ratios R(D) and R(D*), " << std::endl;
27  std::cout << "and polarization P_tau and P_D* for B->D(*)taunu decays " << std::endl;
28  std::cout << "using the BSTD model." << std::endl;
29  std::cout << std::endl;
30  std::cout << "The results are written to " << std::endl;
31  std::cout << "dstar.{w,costau,2D}.dat, d.{w,costau,2D}.dat " << std::endl;
32  std::cout << "and r_pol.dat in the ASCII text format." << std::endl;
33  std::cout << std::endl;
34  std::cout << "Options: " << std::endl;
35  std::cout << " -h\tShow this usage." << std::endl;
36  std::cout << " -r\tCalculate only R and polarizations, without writing output files." << std::endl;
37 }
38 
40 int main(int argc, char* argv[])
41 {
42 
43  bool calcROnly = false;
44 
45  while (1) {
46  // get options
47  int c = getopt(argc, argv, "rh");
48  if (c == -1)break; //end of parameter.
49  switch (c) {
50  case 'h':
51  usage();
52  exit(0);
53  break;
54  case 'r':
55  calcROnly = true;
56  break;
57  default:
58  case '?':
59  std::cerr << "Unknown option" << std::endl;
60  usage();
61  exit(1);
62  }
63  }
64 
65  const int narg = argc - optind;
66  if (narg != 13) {
67  usage();
68  exit(1);
69  }
70 
71  const double m_B(atof(argv[optind++]));
72  const double m_Dst(atof(argv[optind++]));
73  const double m_D(atof(argv[optind++]));
74  const double m_tau(atof(argv[optind++]));
75  const double m_lep(atof(argv[optind++]));
76 
77  const double rho12(atof(argv[optind++]));
78  const double rhoA12(atof(argv[optind++]));
79  const double ffR11(atof(argv[optind++]));
80  const double ffR21(atof(argv[optind++]));
81  const double aS1(atof(argv[optind++]));
82  const double aR3(atof(argv[optind++]));
83  const double m_b(atof(argv[optind++]));
84  const double m_c(atof(argv[optind++]));
85 
86  // double m_tau(1.777),mB(5.2795),mDstar(2.01027),mD(1.86962); // B0
87  // double m_tau(1.777),mB(5.2791),mDstar(2.00697),mD(1.86484); // B+
88  double AbsCV1, ArgCV1, AbsCV2, ArgCV2, AbsCS1, ArgCS1, AbsCS2, ArgCS2, AbsCT, ArgCT;
89 
90  std::cout << "Input Abs(CV1) Arg(CV1) Abs(CV2) Arg(CV2) Abs(CS1) Arg(CS1) Abs(CS2) Arg(CS2) Abs(CT) Arg(CT) >" << std::endl;
91  std::cin >> AbsCV1 >> ArgCV1 >> AbsCV2 >> ArgCV2 >> AbsCS1 >> ArgCS1 >> AbsCS2 >> ArgCS2 >> AbsCT >> ArgCT;
92  EvtComplex Coeffs[5];
93  Coeffs[0] = EvtComplex(AbsCV1 * cos(ArgCV1), AbsCV1 * sin(ArgCV1));
94  Coeffs[1] = EvtComplex(AbsCV2 * cos(ArgCV2), AbsCV2 * sin(ArgCV2));
95  Coeffs[2] = EvtComplex(AbsCS1 * cos(ArgCS1), AbsCS1 * sin(ArgCS1));
96  Coeffs[3] = EvtComplex(AbsCS2 * cos(ArgCS2), AbsCS2 * sin(ArgCS2));
97  Coeffs[4] = EvtComplex(AbsCT * cos(ArgCT), AbsCT * sin(ArgCT));
98 
99  EvtBSemiTauonicHelicityAmplitudeCalculator BSTD(rho12, rhoA12, ffR11, ffR21, aS1, aR3, m_b, m_c,
100  Coeffs[0], Coeffs[1], Coeffs[2], Coeffs[3], Coeffs[4],
101  m_B, m_D, m_Dst);
102  EvtBSemiTauonicDecayRateCalculator CalcDecayRate;
103 
104  std::cout << "tau mass: " << m_tau << std::endl;
105  std::cout << "light lepton mass: " << m_lep << std::endl;
106 
107  if (!calcROnly) {
108  // calculate and write differential decay rate to files
109  const double totgamd = CalcDecayRate.GammaD(BSTD, m_tau);
110 
111  const double totgamds = CalcDecayRate.GammaDstar(BSTD, m_tau);
112 
113  int i_tauhel[2] = { -1, 1}; // left, right
114  int i_dhel[4] = { -1, 0, 1, 2};
115  double dgam[2][4];
116  double wmin = 1.0; //BSTD.wmin();
117  double wmax = 1.35; //BSTD.wmax(m_tau,1);
118  double dw = (wmax - wmin) / 70.;
119 
120  // dGam/dw
121  std::cout << "Calculating dGam/dw ..." << std::endl;
122  std::ofstream ofile1("dstar.w.dat");
123  ofile1 << "# w dgamlm dgamrm dgaml0 dgamr0 dgamlp dgamrp total" << std::endl;
124  for (double w = wmin + dw / 2.;
125  w < wmax; w += dw) {
126  double tot(0);
127  for (int i = 0; i < 2; i++) {
128  for (int j = 0; j < 3; j++) {
129  dgam[i][j] = CalcDecayRate.dGammadw(BSTD, m_tau, i_tauhel[i], i_dhel[j], w) / totgamds * dw;
130  tot += dgam[i][j];
131  }
132  }
133  ofile1 << w
134  << " " << dgam[0][0] << " " << dgam[1][0]
135  << " " << dgam[0][1] << " " << dgam[1][1]
136  << " " << dgam[0][2] << " " << dgam[1][2]
137  << " " << tot << std::endl;
138  }
139 
140  std::ofstream ofile2("d.w.dat");
141  ofile2 << "# w dgaml dgamr total" << std::endl;
142  wmin = 1.0; //BSTD.wmin();
143  wmax = 1.42; //BSTD.wmax(m_tau,2);
144  dw = (wmax - wmin) / 84;
145  for (double w = wmin + dw / 2.;
146  w < wmax; w += dw) {
147  double tot(0);
148  for (int i = 0; i < 2; i++) {
149  int j = 3;
150  dgam[i][j] = CalcDecayRate.dGammadw(BSTD, m_tau, i_tauhel[i], i_dhel[j], w) / totgamd * dw;
151  tot += dgam[i][j];
152  }
153  ofile2 << w
154  << " " << dgam[0][3] << " " << dgam[1][3]
155  << " " << tot << std::endl;
156  }
157 
158  // dGam/dcostau
159  std::cout << "Calculating dGam/dcostau ..." << std::endl;
160  double dcostau = (1 - (-1)) / 100.;
161  std::ofstream ofilecostau1("dstar.costau.dat");
162  ofilecostau1 << "# costau dgamlm dgamrm dgaml0 dgamr0 dgamlp dgamrp total" << std::endl;
163  for (double costau = -1 + dcostau / 2.;
164  costau < 1; costau += dcostau) {
165  double tot(0);
166  for (int i = 0; i < 2; i++) {
167  for (int j = 0; j < 3; j++) {
168  dgam[i][j] = CalcDecayRate.dGammadcostau(BSTD, m_tau, i_tauhel[i], i_dhel[j], costau) / totgamds * dcostau;
169  tot += dgam[i][j];
170  }
171  }
172  ofilecostau1 << costau
173  << " " << dgam[0][0] << " " << dgam[1][0]
174  << " " << dgam[0][1] << " " << dgam[1][1]
175  << " " << dgam[0][2] << " " << dgam[1][2]
176  << " " << tot << std::endl;
177  }
178 
179  std::ofstream ofilecostau2("d.costau.dat");
180  ofilecostau2 << "# costau dgaml dgamr total" << std::endl;
181  for (double costau = -1 + dcostau / 2.;
182  costau < 1; costau += dcostau) {
183  double tot(0);
184  for (int i = 0; i < 2; i++) {
185  int j = 3;
186  dgam[i][j] = CalcDecayRate.dGammadcostau(BSTD, m_tau, i_tauhel[i], i_dhel[j], costau) / totgamd * dcostau;
187  tot += dgam[i][j];
188  }
189  ofilecostau2 << costau
190  << " " << dgam[0][3] << " " << dgam[1][3]
191  << " " << tot << std::endl;
192  }
193 
194  // 2D
195  std::cout << "Calculating dGam/dw/dcostau ..." << std::endl;
196  wmin = 1.0;
197  wmax = 1.35;
198  dw = (wmax - wmin) / 20.;
199  dcostau = (1 - (-1)) / 20.;
200  std::ofstream ofile2d1("dstar.2d.dat");
201  ofile2d1 << "# w costau dgamlm dgamrm dgaml0 dgamr0 dgamlp dgamrp total" << std::endl;
202  for (double w = wmin + dw / 2.; w < wmax; w += dw) {
203  for (double costau = -1 + dcostau / 2.; costau < 1; costau += dcostau) {
204  double tot(0);
205  for (int i = 0; i < 2; i++) {
206  for (int j = 0; j < 3; j++) {
207  dgam[i][j] = CalcDecayRate.dGammadwdcostau(BSTD, m_tau, i_tauhel[i], i_dhel[j], w, costau) / totgamds * dcostau * dw;
208  tot += dgam[i][j];
209  }
210  }
211  ofile2d1 << w << " " << costau
212  << " " << dgam[0][0] << " " << dgam[1][0]
213  << " " << dgam[0][1] << " " << dgam[1][1]
214  << " " << dgam[0][2] << " " << dgam[1][2]
215  << " " << tot << std::endl;
216  }
217  }
218 
219  wmin = 1.0;
220  wmax = 1.42;
221  dw = (wmax - wmin) / 20.;
222  dcostau = (1 - (-1)) / 20.;
223  std::ofstream ofile2d2("d.2d.dat");
224  ofile2d2 << "# w costau dgaml dgamr total" << std::endl;
225  for (double w = wmin + dw / 2.; w < wmax; w += dw) {
226  for (double costau = -1 + dcostau / 2.; costau < 1; costau += dcostau) {
227  double tot(0);
228  for (int i = 0; i < 2; i++) {
229  int j = 3;
230  dgam[i][j] = CalcDecayRate.dGammadwdcostau(BSTD, m_tau, i_tauhel[i], i_dhel[j], w, costau) / totgamd * dcostau * dw;
231  tot += dgam[i][j];
232  }
233  ofile2d2 << w << " " << costau
234  << " " << dgam[0][3] << " " << dgam[1][3]
235  << " " << tot << std::endl;
236  }
237  }
238  }
239 
240  // R and polarization data
241  /* cppcheck-suppress duplicateCondition */
242  if (!calcROnly) {
243  std::ofstream ofile3("r_pol.dat");
244  ofile3 << "# rgamd rgamds ptaud ptauds pds" << std::endl;
245  ofile3 << CalcDecayRate.RGammaD(BSTD, m_tau, m_lep) << " ";
246  ofile3 << CalcDecayRate.RGammaDstar(BSTD, m_tau, m_lep) << " ";
247  ofile3 << CalcDecayRate.PtauD(BSTD, m_tau) << " ";
248  ofile3 << CalcDecayRate.PtauDstar(BSTD, m_tau) << " ";
249  ofile3 << CalcDecayRate.PDstar(BSTD, m_tau);
250  ofile3 << std::endl;
251  }
252 
253  std::cout << "# rgamd rgamds ptaud ptauds pds" << std::endl;
254  std::cout << CalcDecayRate.RGammaD(BSTD, m_tau, m_lep) << " ";
255  std::cout << CalcDecayRate.RGammaDstar(BSTD, m_tau, m_lep) << " ";
256  std::cout << CalcDecayRate.PtauD(BSTD, m_tau) << " ";
257  std::cout << CalcDecayRate.PtauDstar(BSTD, m_tau) << " ";
258  std::cout << CalcDecayRate.PDstar(BSTD, m_tau);
259  std::cout << std::endl;
260 
261  return 0;
262 }
263 
Belle2::EvtBSemiTauonicDecayRateCalculator::PtauD
double PtauD(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau)
Function calculates the polarization of tau, (RH - LH)/(LH + RH), in B->Dtaunu decay.
Definition: EvtBSemiTauonicDecayRateCalculator.cc:180
Belle2::EvtBSemiTauonicDecayRateCalculator::RGammaD
double RGammaD(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau, const double mlep=0.0005110)
Function calculates the ratio of Br(B->Dtaunu)/Br(B->Dlnu), R(D).
Definition: EvtBSemiTauonicDecayRateCalculator.cc:159
Belle2::EvtBSemiTauonicDecayRateCalculator::GammaDstar
double GammaDstar(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau)
Function calculates the differential decay rate Gamma for D*taunu decay, integrated for w and costau ...
Definition: EvtBSemiTauonicDecayRateCalculator.cc:116
Belle2::EvtBSemiTauonicDecayRateCalculator::dGammadw
double dGammadw(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double w)
Function calculates the differential decay rate dGamma/dw, integrated for costau.
Definition: EvtBSemiTauonicDecayRateCalculator.cc:51
Belle2::EvtBSemiTauonicDecayRateCalculator::PtauDstar
double PtauDstar(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau)
Function calculates the polarization of tau, (RH - LH)/(LH + RH), in B->D*taunu decay.
Definition: EvtBSemiTauonicDecayRateCalculator.cc:187
Belle2::EvtBSemiTauonicDecayRateCalculator::dGammadwdcostau
double dGammadwdcostau(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double w, double costau)
Function calculates the differential decay rate dGamma/dw/dcostau.
Definition: EvtBSemiTauonicDecayRateCalculator.cc:41
Belle2::EvtBSemiTauonicDecayRateCalculator::RGammaDstar
double RGammaDstar(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau, const double mlep=0.0005110)
Function calculates the ratio of Br(B->Dtaunu)/Br(B->Dlnu), R(D*).
Definition: EvtBSemiTauonicDecayRateCalculator.cc:168
Belle2::EvtBSemiTauonicDecayRateCalculator::GammaD
double GammaD(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau)
Function calculates the decay rate Gamma for Dtaunu decay, integrated for w and costau and summed for...
Definition: EvtBSemiTauonicDecayRateCalculator.cc:108
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::EvtBSemiTauonicDecayRateCalculator::PDstar
double PDstar(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau)
Function calculates the polarization of D*, longitudinal/(longitudinal + transverse),...
Definition: EvtBSemiTauonicDecayRateCalculator.cc:197
Belle2::EvtBSemiTauonicDecayRateCalculator::dGammadcostau
double dGammadcostau(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double costau)
Function calculates the differential decay rate dGamma/dcostau, integrated for w.
Definition: EvtBSemiTauonicDecayRateCalculator.cc:69
Belle2::EvtBSemiTauonicDecayRateCalculator
Class to calculate the differential decay rate, R(D), R(D*), polarizations of tau and D* using BSTD m...
Definition: EvtBSemiTauonicDecayRateCalculator.h:32
Belle2::EvtBSemiTauonicHelicityAmplitudeCalculator
The class calculates the helicity amplitude of semi-tauonic B decays including new physics effects ba...
Definition: EvtBSemiTauonicHelicityAmplitudeCalculator.h:38