]>
Commit | Line | Data |
---|---|---|
0ca57c2f | 1 | //-------------------------------------------------------------------------- |
2 | // | |
3 | // Module: EvtD0gammaDalitz.cc | |
4 | // | |
5 | // Modification history: | |
6 | // | |
7 | // JGT February 13, 2012 Module created | |
8 | // | |
9 | //------------------------------------------------------------------------ | |
10 | ||
11 | #include <cstdlib> | |
12 | #include <string> | |
13 | #include <cmath> | |
14 | ||
15 | #include "EvtGenBase/EvtPatches.hh" | |
16 | ||
17 | #include "EvtGenBase/EvtParticle.hh" | |
18 | #include "EvtGenBase/EvtGenKine.hh" | |
19 | #include "EvtGenBase/EvtPDL.hh" | |
20 | #include "EvtGenBase/EvtReport.hh" | |
21 | #include "EvtGenBase/EvtResonance.hh" | |
22 | #include "EvtGenBase/EvtResonance2.hh" | |
23 | #include "EvtGenModels/EvtD0gammaDalitz.hh" | |
24 | #include "EvtGenBase/EvtConst.hh" | |
25 | #include "EvtGenBase/EvtFlatte.hh" | |
26 | #include "EvtGenBase/EvtDecayTable.hh" | |
27 | ||
28 | ||
29 | // Initialize the static variables. | |
30 | const EvtSpinType::spintype& EvtD0gammaDalitz::_SCALAR = EvtSpinType::SCALAR; | |
31 | const EvtSpinType::spintype& EvtD0gammaDalitz::_VECTOR = EvtSpinType::VECTOR; | |
32 | const EvtSpinType::spintype& EvtD0gammaDalitz::_TENSOR = EvtSpinType::TENSOR; | |
33 | ||
34 | const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_EtaPic = EvtDalitzReso::EtaPic; | |
35 | const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_PicPicKK = EvtDalitzReso::PicPicKK; | |
36 | ||
37 | const EvtDalitzReso::NumType& EvtD0gammaDalitz::_RBW = EvtDalitzReso::RBW_CLEO_ZEMACH; | |
38 | const EvtDalitzReso::NumType& EvtD0gammaDalitz::_GS = EvtDalitzReso::GS_CLEO_ZEMACH; | |
39 | const EvtDalitzReso::NumType& EvtD0gammaDalitz::_KMAT = EvtDalitzReso::K_MATRIX; | |
40 | ||
41 | const EvtCyclic3::Pair& EvtD0gammaDalitz::_AB = EvtCyclic3::AB; | |
42 | const EvtCyclic3::Pair& EvtD0gammaDalitz::_AC = EvtCyclic3::AC; | |
43 | const EvtCyclic3::Pair& EvtD0gammaDalitz::_BC = EvtCyclic3::BC; | |
44 | ||
45 | ||
46 | EvtD0gammaDalitz::EvtD0gammaDalitz() | |
47 | { | |
48 | /* Empty constructor. */ | |
49 | } | |
50 | ||
51 | ||
52 | EvtD0gammaDalitz::~EvtD0gammaDalitz() | |
53 | { | |
54 | /* Empty destructor. */ | |
55 | } | |
56 | ||
57 | ||
58 | std::string EvtD0gammaDalitz::getName() | |
59 | { | |
60 | return "D0GAMMADALITZ"; | |
61 | } | |
62 | ||
63 | ||
64 | EvtDecayBase* EvtD0gammaDalitz::clone() | |
65 | { | |
66 | return new EvtD0gammaDalitz; | |
67 | } | |
68 | ||
69 | ||
70 | void EvtD0gammaDalitz::init() | |
71 | { | |
72 | // check that there are 0 arguments | |
73 | checkNArg(0); | |
74 | ||
75 | // Check that this model is valid for the specified decay. | |
76 | checkNDaug( 3 ); | |
77 | checkSpinParent ( _SCALAR ); | |
78 | checkSpinDaughter( 0, _SCALAR ); | |
79 | checkSpinDaughter( 1, _SCALAR ); | |
80 | checkSpinDaughter( 2, _SCALAR ); | |
81 | ||
82 | // Get the values of the EvtId objects from the data files. | |
83 | readPDGValues(); | |
84 | ||
85 | // Get the EvtId of the D0 and its 3 daughters. | |
86 | getParentId(); | |
87 | ||
88 | EvtId dau[ 3 ]; | |
89 | for ( int index = 0; index < 3; index++ ) | |
90 | { | |
91 | dau[ index ] = getDaug( index ); | |
92 | } | |
93 | ||
94 | // Look for K0bar h+ h-. The order will be K[0SL] h+ h- | |
95 | for ( int index = 0; index < 3; index++ ) | |
96 | { | |
97 | if ( ( dau[ index ] == _K0B ) || ( dau[ index ] == _KS ) || ( dau[ index ] == _KL ) ) | |
98 | { | |
99 | _d1 = index; | |
100 | } | |
101 | else if ( ( dau[ index ] == _PIP ) || ( dau[ index ] == _KP ) ) | |
102 | { | |
103 | _d2 = index; | |
104 | } | |
105 | else if ( ( dau[ index ] == _PIM ) || ( dau[ index ] == _KM ) ) | |
106 | { | |
107 | _d3 = index; | |
108 | } | |
109 | else | |
110 | { | |
111 | reportInvalidAndExit(); | |
112 | } | |
113 | } | |
114 | ||
115 | // Check if we're dealing with Ks pi pi or with Ks K K. | |
116 | _isKsPiPi = false; | |
117 | if ( dau[ _d2 ] == _PIP || dau[ _d2 ] == _PIM ) | |
118 | { | |
119 | _isKsPiPi = true; | |
120 | } | |
121 | ||
122 | } | |
123 | ||
124 | ||
125 | void EvtD0gammaDalitz::initProbMax() | |
126 | { | |
127 | setProbMax( 5200. ); | |
128 | } | |
129 | ||
130 | ||
131 | void EvtD0gammaDalitz::decay( EvtParticle* part ) | |
132 | { | |
133 | // Check if the D is from a B+- -> D0 K+- decay with the appropriate model. | |
134 | EvtParticle* parent = part->getParent(); // If there are no mistakes, should be B+ or B-. | |
135 | if (parent != 0 && EvtDecayTable::getInstance()->getDecayFunc( parent )->getName() == "BTODDALITZCPK" ) | |
136 | { | |
137 | EvtId parId = parent->getId(); | |
138 | if ( ( parId == _BP ) || ( parId == _BM ) || | |
139 | ( parId == _B0 ) || ( parId == _B0B) ) | |
140 | { | |
141 | _bFlavor = parId; | |
142 | } | |
143 | else | |
144 | { | |
145 | reportInvalidAndExit(); | |
146 | } | |
147 | } | |
148 | else | |
149 | { | |
150 | reportInvalidAndExit(); | |
151 | } | |
152 | ||
153 | // Read the D decay parameters from the B decay model. | |
154 | // Gamma angle in rad. | |
155 | double gamma = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 0 ); | |
156 | // Strong phase in rad. | |
157 | double delta = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 1 ); | |
158 | // Ratio between B->D0K and B->D0barK | |
159 | double rB = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 2 ); | |
160 | ||
161 | // Same structure for all of these decays. | |
162 | part->initializePhaseSpace( getNDaug(), getDaugs() ); | |
163 | EvtVector4R pA = part->getDaug( _d1 )->getP4(); | |
164 | EvtVector4R pB = part->getDaug( _d2 )->getP4(); | |
165 | EvtVector4R pC = part->getDaug( _d3 )->getP4(); | |
166 | ||
167 | // Squared invariant masses. | |
168 | double mSqAB = ( pA + pB ).mass2(); | |
169 | double mSqAC = ( pA + pC ).mass2(); | |
170 | double mSqBC = ( pB + pC ).mass2(); | |
171 | ||
172 | EvtComplex amp( 1.0, 0.0 ); | |
173 | ||
174 | // Direct and conjugated amplitudes. | |
175 | EvtComplex ampDir; | |
176 | EvtComplex ampCnj; | |
177 | ||
178 | if ( _isKsPiPi ) | |
179 | { | |
180 | // Direct and conjugated Dalitz points. | |
181 | EvtDalitzPoint pointDir( _mKs, _mPi, _mPi, mSqAB, mSqBC, mSqAC ); | |
182 | EvtDalitzPoint pointCnj( _mKs, _mPi, _mPi, mSqAC, mSqBC, mSqAB ); | |
183 | ||
184 | // Direct and conjugated amplitudes. | |
185 | ampDir = dalitzKsPiPi( pointDir ); | |
186 | ampCnj = dalitzKsPiPi( pointCnj ); | |
187 | } | |
188 | else | |
189 | { | |
190 | // Direct and conjugated Dalitz points. | |
191 | EvtDalitzPoint pointDir( _mKs, _mK, _mK, mSqAB, mSqBC, mSqAC ); | |
192 | EvtDalitzPoint pointCnj( _mKs, _mK, _mK, mSqAC, mSqBC, mSqAB ); | |
193 | ||
194 | // Direct and conjugated amplitudes. | |
195 | ampDir = dalitzKsKK( pointDir ); | |
196 | ampCnj = dalitzKsKK( pointCnj ); | |
197 | } | |
198 | ||
199 | if ( _bFlavor == _BP || _bFlavor == _B0 ) | |
200 | { | |
201 | amp = ampCnj + rB * exp( EvtComplex( 0., delta + gamma ) ) * ampDir; | |
202 | } | |
203 | else | |
204 | { | |
205 | amp = ampDir + rB * exp( EvtComplex( 0., delta - gamma ) ) * ampCnj; | |
206 | } | |
207 | ||
208 | vertex( amp ); | |
209 | ||
210 | return; | |
211 | ||
212 | } | |
213 | ||
214 | ||
215 | EvtComplex EvtD0gammaDalitz::dalitzKsPiPi( const EvtDalitzPoint& point ) const | |
216 | { | |
217 | static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 ); | |
218 | ||
219 | EvtComplex amp = 0.; | |
220 | ||
221 | // This corresponds to relativistic Breit-Wigner distributions. Not K-matrix. | |
222 | // Defining resonances. | |
223 | static EvtDalitzReso KStarm ( plot, _BC, _AC, _VECTOR, 0.893606, 0.0463407, _RBW ); | |
224 | static EvtDalitzReso KStarp ( plot, _BC, _AB, _VECTOR, 0.893606, 0.0463407, _RBW ); | |
225 | static EvtDalitzReso rho0 ( plot, _AC, _BC, _VECTOR, 0.7758 , 0.1464 , _GS ); | |
226 | static EvtDalitzReso omega ( plot, _AC, _BC, _VECTOR, 0.78259 , 0.00849 , _RBW ); | |
227 | static EvtDalitzReso f0_980 ( plot, _AC, _BC, _SCALAR, 0.975 , 0.044 , _RBW ); | |
228 | static EvtDalitzReso f0_1370 ( plot, _AC, _BC, _SCALAR, 1.434 , 0.173 , _RBW ); | |
229 | static EvtDalitzReso f2_1270 ( plot, _AC, _BC, _TENSOR, 1.2754 , 0.1851 , _RBW ); | |
230 | static EvtDalitzReso K0Starm_1430( plot, _BC, _AC, _SCALAR, 1.459 , 0.175 , _RBW ); | |
231 | static EvtDalitzReso K0Starp_1430( plot, _BC, _AB, _SCALAR, 1.459 , 0.175 , _RBW ); | |
232 | static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256 , 0.0985 , _RBW ); | |
233 | static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256 , 0.0985 , _RBW ); | |
234 | static EvtDalitzReso sigma ( plot, _AC, _BC, _SCALAR, 0.527699, 0.511861 , _RBW ); | |
235 | static EvtDalitzReso sigma2 ( plot, _AC, _BC, _SCALAR, 1.03327 , 0.0987890, _RBW ); | |
236 | static EvtDalitzReso KStarm_1680 ( plot, _BC, _AC, _VECTOR, 1.677 , 0.205 , _RBW ); | |
237 | ||
238 | // Adding terms to the amplitude with their corresponding amplitude and phase terms. | |
239 | amp += EvtComplex( .848984 , .893618 ); | |
240 | amp += EvtComplex( -1.16356 , 1.19933 ) * KStarm .evaluate( point ); | |
241 | amp += EvtComplex( .106051 , - .118513 ) * KStarp .evaluate( point ); | |
242 | amp += EvtComplex( 1.0 , 0.0 ) * rho0 .evaluate( point ); | |
243 | amp += EvtComplex( - .0249569, .0388072 ) * omega .evaluate( point ); | |
244 | amp += EvtComplex( - .423586 , - .236099 ) * f0_980 .evaluate( point ); | |
245 | amp += EvtComplex( -2.16486 , 3.62385 ) * f0_1370 .evaluate( point ); | |
246 | amp += EvtComplex( .217748 , - .133327 ) * f2_1270 .evaluate( point ); | |
247 | amp += EvtComplex( 1.62128 , 1.06816 ) * K0Starm_1430.evaluate( point ); | |
248 | amp += EvtComplex( .148802 , .0897144 ) * K0Starp_1430.evaluate( point ); | |
249 | amp += EvtComplex( 1.15489 , - .773363 ) * K2Starm_1430.evaluate( point ); | |
250 | amp += EvtComplex( .140865 , - .165378 ) * K2Starp_1430.evaluate( point ); | |
251 | amp += EvtComplex( -1.55556 , - .931685 ) * sigma .evaluate( point ); | |
252 | amp += EvtComplex( - .273791 , - .0535596 ) * sigma2 .evaluate( point ); | |
253 | amp += EvtComplex( -1.69720 , .128038 ) * KStarm_1680 .evaluate( point ); | |
254 | ||
255 | return amp; | |
256 | } | |
257 | ||
258 | ||
259 | EvtComplex EvtD0gammaDalitz::dalitzKsKK( const EvtDalitzPoint& point ) const | |
260 | { | |
261 | static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 ); | |
262 | ||
263 | // Defining resonances. | |
264 | static EvtDalitzReso a00_980 ( plot, _AC, _BC, _SCALAR, 0.999 , _RBW, .550173, .324, _EtaPic ); | |
265 | static EvtDalitzReso phi ( plot, _AC, _BC, _VECTOR, 1.01943, .00459319 , _RBW ); | |
266 | static EvtDalitzReso a0p_980 ( plot, _AC, _AB, _SCALAR, 0.999 , _RBW, .550173, .324, _EtaPic ); | |
267 | static EvtDalitzReso f0_1370 ( plot, _AC, _BC, _SCALAR, 1.350 , .265 , _RBW ); | |
268 | static EvtDalitzReso a0m_980 ( plot, _AB, _AC, _SCALAR, 0.999 , _RBW, .550173, .324, _EtaPic ); | |
269 | static EvtDalitzReso f0_980 ( plot, _AC, _BC, _SCALAR, 0.965 , _RBW, .695 , .165, _PicPicKK ); | |
270 | static EvtDalitzReso f2_1270 ( plot, _AC, _BC, _TENSOR, 1.2754 , .1851 , _RBW ); | |
271 | static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474 , .265 , _RBW ); | |
272 | static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474 , .265 , _RBW ); | |
273 | static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474 , .265 , _RBW ); | |
274 | ||
275 | // Adding terms to the amplitude with their corresponding amplitude and phase terms. | |
276 | EvtComplex amp( 0., 0. ); // Phase space amplitude. | |
277 | amp += EvtComplex( 1.0 , 0.0 ) * a00_980 .evaluate( point ); | |
278 | amp += EvtComplex( -.126314 , .188701 ) * phi .evaluate( point ); | |
279 | amp += EvtComplex( -.561428 , .0135338 ) * a0p_980 .evaluate( point ); | |
280 | amp += EvtComplex( .035 , -.00110488 ) * f0_1370 .evaluate( point ); | |
281 | amp += EvtComplex( -.0872735 , .0791190 ) * a0m_980 .evaluate( point ); | |
282 | amp += EvtComplex( 0. , 0. ) * f0_980 .evaluate( point ); | |
283 | amp += EvtComplex( .257341 , -.0408343 ) * f2_1270 .evaluate( point ); | |
284 | amp += EvtComplex( -.0614342 , -.649930 ) * a00_1450.evaluate( point ); | |
285 | amp += EvtComplex( -.104629 , .830120 ) * a0p_1450.evaluate( point ); | |
286 | amp += EvtComplex( 0. , 0. ) * a0m_1450.evaluate( point ); | |
287 | ||
288 | return 2.8 * amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi. | |
289 | } | |
290 | ||
291 | ||
292 | void EvtD0gammaDalitz::readPDGValues() | |
293 | { | |
294 | // Define the EvtIds. | |
295 | _BP = EvtPDL::getId( "B+" ); | |
296 | _BM = EvtPDL::getId( "B-" ); | |
297 | _B0 = EvtPDL::getId( "B0" ); | |
298 | _B0B = EvtPDL::getId( "anti-B0" ); | |
299 | _D0 = EvtPDL::getId( "D0" ); | |
300 | _D0B = EvtPDL::getId( "anti-D0" ); | |
301 | _KM = EvtPDL::getId( "K-" ); | |
302 | _KP = EvtPDL::getId( "K+" ); | |
303 | _K0 = EvtPDL::getId( "K0" ); | |
304 | _K0B = EvtPDL::getId( "anti-K0" ); | |
305 | _KL = EvtPDL::getId( "K_L0" ); | |
306 | _KS = EvtPDL::getId( "K_S0" ); | |
307 | _PIM = EvtPDL::getId( "pi-" ); | |
308 | _PIP = EvtPDL::getId( "pi+" ); | |
309 | ||
310 | // Read the relevant masses. | |
311 | _mD0 = EvtPDL::getMass( _D0 ); | |
312 | _mKs = EvtPDL::getMass( _KS ); | |
313 | _mPi = EvtPDL::getMass( _PIP ); | |
314 | _mK = EvtPDL::getMass( _KP ); | |
315 | } | |
316 | ||
317 | ||
318 | void EvtD0gammaDalitz::reportInvalidAndExit() const | |
319 | { | |
320 | report( ERROR, "EvtD0gammaDalitz" ) << "EvtD0gammaDalitz: Invalid mode." << std::endl; | |
321 | exit( 1 ); | |
322 | } |