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