evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtDalitzPoint.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
22 
23 #include "EvtGenBase/EvtPatches.hh"
24 
25 #include <assert.h>
26 #include <math.h>
27 #include <stdio.h>
28 using namespace EvtCyclic3;
29 
31  _mA( -1. ), _mB( -1. ), _mC( -1. ), _qAB( -1. ), _qBC( -1. ), _qCA( -1. )
32 {
33 }
34 
35 EvtDalitzPoint::EvtDalitzPoint( double mA, double mB, double mC, double qAB,
36  double qBC, double qCA ) :
37  _mA( mA ), _mB( mB ), _mC( mC ), _qAB( qAB ), _qBC( qBC ), _qCA( qCA )
38 {
39 }
40 
41 // Constructor from Zemach coordinates
42 
43 EvtDalitzPoint::EvtDalitzPoint( double mA, double mB, double mC,
44  EvtCyclic3::Pair i, double qres, double qhel,
45  double qsum ) :
46  _mA( mA ), _mB( mB ), _mC( mC )
47 {
48  double qi = qres + qsum / 3.;
49  double qj = -qres / 2. + qhel + qsum / 3.;
50  double qk = -qres / 2. - qhel + qsum / 3.;
51 
52  if ( i == AB ) {
53  _qAB = qi;
54  _qBC = qj;
55  _qCA = qk;
56  } else if ( i == BC ) {
57  _qAB = qk;
58  _qBC = qi;
59  _qCA = qj;
60  } else if ( i == CA ) {
61  _qAB = qj;
62  _qBC = qk;
63  _qCA = qi;
64  }
65 }
66 
68  const EvtDalitzCoord& x ) :
69  _mA( dp.m( A ) ), _mB( dp.m( B ) ), _mC( dp.m( C ) )
70 {
71  if ( x.pair1() == AB )
72  _qAB = x.q1();
73  else if ( x.pair2() == AB )
74  _qAB = x.q2();
75  else
76  _qAB = dp.sum() - x.q1() - x.q2();
77 
78  if ( x.pair1() == BC )
79  _qBC = x.q1();
80  else if ( x.pair2() == BC )
81  _qBC = x.q2();
82  else
83  _qBC = dp.sum() - x.q1() - x.q2();
84 
85  if ( x.pair1() == CA )
86  _qCA = x.q1();
87  else if ( x.pair2() == CA )
88  _qCA = x.q2();
89  else
90  _qCA = dp.sum() - x.q1() - x.q2();
91 }
92 
94 {
95  double ret = _qAB;
96  if ( BC == i )
97  ret = _qBC;
98  else if ( CA == i )
99  ret = _qCA;
100 
101  return ret;
102 }
103 
105 {
106  double ret = _mA;
107  if ( B == i )
108  ret = _mB;
109  else if ( C == i )
110  ret = _mC;
111 
112  return ret;
113 }
114 
115 // Zemach variables
116 
118 {
119  return ( 2. * q( i ) - q( EvtCyclic3::prev( i ) ) -
120  q( EvtCyclic3::next( i ) ) ) /
121  3.;
122 }
124 {
125  Pair j = next( i );
126  Pair k = prev( i );
127  return ( q( j ) - q( k ) ) / 2.;
128 }
129 double EvtDalitzPoint::qsum() const
130 {
131  return _qAB + _qBC + _qCA;
132 }
133 
135 {
137  return dp.qMin( i, j, q( j ) );
138 }
139 
141 {
143  return dp.qMax( i, j, q( j ) );
144 }
145 
147 {
148  if ( i == j )
149  return m( i ) * m( i );
150  else
151  return ( q( combine( i, j ) ) - m( i ) * m( i ) - m( j ) * m( j ) ) / 2.;
152 }
153 
155 {
157  return dp.e( i, j, q( j ) );
158 }
159 
161 {
163  return dp.p( i, j, q( j ) );
164 }
165 
167  EvtCyclic3::Pair pairRes ) const
168 {
170  return dp.cosTh( pairAng, q( pairAng ), pairRes, q( pairRes ) );
171 }
172 
174  EvtCyclic3::Pair j ) const
175 {
176  return EvtDalitzCoord( i, q( i ), j, q( j ) );
177 }
178 
180 {
181  return EvtDalitzPlot( _mA, _mB, _mC, bigM() );
182 }
183 
185 {
186  // Check masses
187 
188  double M = bigM();
189  if ( _mA < 0 || _mB < 0 || _mC < 0 || M <= 0 )
190  return false;
191  if ( M < _mA + _mB + _mC )
192  return false;
193 
194  // Check that first coordinate is within absolute limits
195 
196  bool inside = false;
198 
199  if ( dp.qAbsMin( AB ) <= _qAB && _qAB <= dp.qAbsMax( AB ) )
200  if ( qMin( BC, AB ) <= _qBC && _qBC <= qMax( BC, AB ) )
201  inside = true;
202 
203  return inside;
204 }
205 
206 double EvtDalitzPoint::bigM() const
207 {
208  return sqrt( _qAB + _qBC + _qCA - _mA * _mA - _mB * _mB - _mC * _mC );
209 }
210 
212 {
213  getDalitzPlot().print();
214  printf( "%f %f %f\n", _qAB, _qBC, _qCA );
215 }
bool isValid() const
Index prev(Index i)
Definition: EvtCyclic3.cpp:128
double bigM() const
double qhel(EvtCyclic3::Pair i) const
EvtDalitzCoord getDalitzPoint(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
double qres(EvtCyclic3::Pair i) const
void print() const
double qAbsMin(EvtCyclic3::Pair i) const
double m(EvtCyclic3::Index) const
Pair combine(Index i, Index j)
Definition: EvtCyclic3.cpp:208
double qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j, double q) const
EvtDalitzPlot getDalitzPlot() const
double qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
double sum() const
double pp(EvtCyclic3::Index i, EvtCyclic3::Index j) const
double q1() const
Index next(Index i)
Definition: EvtCyclic3.cpp:142
double q2() const
double e(EvtCyclic3::Index i, EvtCyclic3::Pair j, double q) const
double p(EvtCyclic3::Index i, EvtCyclic3::Pair j, double q) const
double q(EvtCyclic3::Pair) const
double qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j, double q) const
double qsum() const
EvtCyclic3::Pair pair1() const
double qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
double p(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
double qAbsMax(EvtCyclic3::Pair i) const
double cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const
EvtCyclic3::Pair pair2() const
double e(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
double cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
void print() const