The function to calculate a quark decay amplitude.
142{
143 static EvtId B0 = EvtPDL::getId("B0");
144 static EvtId B0B = EvtPDL::getId("anti-B0");
145
146
147
148 double rndm = EvtRandom::random();
149
150 if (getNDaug() == 4) {
151 EvtId tempDaug[2];
152
153 if (rndm < 0.5) {
154 tempDaug[0] = getDaug(0);
155 tempDaug[1] = getDaug(3);
156 } else {
157 tempDaug[0] = getDaug(2);
158 tempDaug[1] = getDaug(1);
159 }
160
161 p->initializePhaseSpace(2, tempDaug);
162 } else {
163 p->initializePhaseSpace(2, getDaugs());
164 }
165
166 EvtParticle* s1, *s2;
167
168 s1 = p->getDaug(0);
169 s2 = p->getDaug(1);
170
171
172 if (s1->getNDaug() > 0) {
173 s1->deleteDaughters();
174 }
175 if (s2->getNDaug() > 0) {
176 s2->deleteDaughters();
177 }
178
179 EvtVector4R p1 = s1->getP4();
180 EvtVector4R p2 = s2->getP4();
181
182
183 rndm = EvtRandom::random();
184 int mixed = (rndm < 0.5) ? 1 : 0;
185
186
187
188 if (mixed == 1) {
189 EvtId mixedId = (rndm < 0.25) ? getDaug(0) : getDaug(1);
190 EvtId mixedId2 = mixedId;
191 if (getNDaug() == 4 && rndm < 0.25)
192 mixedId2 = getDaug(2);
193 if (getNDaug() == 4 && rndm > 0.25)
194 mixedId2 = getDaug(3);
195 s1->init(mixedId, p1);
196 s2->init(mixedId2, p2);
197 }
198
199
200
201 if (mixed == 0) {
202 EvtId unmixedId = (rndm < 0.75) ? getDaug(0) : getDaug(1);
203 EvtId unmixedId2 = (rndm < 0.75) ? getDaug(1) : getDaug(0);
204 if (getNDaug() == 4 && rndm < 0.75)
205 unmixedId2 = getDaug(3);
206 if (getNDaug() == 4 && rndm > 0.75)
207 unmixedId2 = getDaug(2);
208 s1->init(unmixedId, p1);
209 s2->init(unmixedId2, p2);
210 }
211
212
213
214
215 s1->setLifetime();
216 s2->setLifetime();
217
218
219 EvtComplex osc_amp;
220
221 if (EvtRandom::random() <
_sdprob) {
222 double t1 = s1->getLifetime();
223 double t2 = s2->getLifetime();
224
225 EvtComplex exp1(0, 0.5 *
_freq * t1);
226 EvtComplex exp2(0, 0.5 *
_freq * t2);
227
228 EvtId finalBsig = (EvtRandom::random() > 0.5) ? B0B : B0;
229
230 EvtComplex g1p = 0.5 * (exp(-exp1) + exp(exp1));
231 EvtComplex g1m = 0.5 * (exp(-exp1) - exp(exp1));
232 EvtComplex g2p = 0.5 * (exp(-exp2) + exp(exp2));
233 EvtComplex g2m = 0.5 * (exp(-exp2) - exp(exp2));
234
235 EvtComplex BB = g1p * g2p;
236 EvtComplex barBB = g1m * g2p;
237 EvtComplex BbarB = g1p * g2m;
238 EvtComplex barBbarB = g1m * g2m;
239
240 if (!mixed && finalBsig == B0) {
241 osc_amp = BB;
242 }
243 if (!mixed && finalBsig == B0B) {
244 osc_amp = barBbarB;
245 }
246 if (mixed && finalBsig == B0) {
247 osc_amp = barBB;
248 }
249 if (mixed && finalBsig == B0B) {
250 osc_amp = BbarB;
251 }
252 } else {
253 double dct = s1->getLifetime() - s2->getLifetime();
254
255 EvtComplex exp1(0, 0.5 *
_freq * dct);
256
257 EvtId stateAtDeltaTeq0 = (s2->getId() == B0) ? B0B : B0;
258
259
260 EvtComplex gp = 0.5 * (exp(-1.0 * exp1) + exp(exp1));
261 EvtComplex gm = 0.5 * (exp(-1.0 * exp1) - exp(exp1));
262 EvtComplex sqz = exp(EvtComplex(0, 0.5));
263
264 EvtComplex BB = gp;
265 EvtComplex barBB = -sqz * gm;
266 EvtComplex BbarB = -sqz * gm;
267 EvtComplex barBbarB = gp;
268
269 if (!mixed && stateAtDeltaTeq0 == B0) {
270 osc_amp = BB;
271 }
272 if (!mixed && stateAtDeltaTeq0 == B0B) {
273 osc_amp = barBbarB;
274 }
275
276 if (mixed && stateAtDeltaTeq0 == B0) {
277 osc_amp = barBB;
278 }
279 if (mixed && stateAtDeltaTeq0 == B0B) {
280 osc_amp = BbarB;
281 }
282
283 }
284
285 double norm = 1.0 / p1.d3mag();
286 vertex(0, norm * osc_amp * p1 * (p->eps(0)));
287 vertex(1, norm * osc_amp * p1 * (p->eps(1)));
288 vertex(2, norm * osc_amp * p1 * (p->eps(2)));
289
290 return;
291}
double _sdprob
probability of the spontaneous disentanglement [0,1]
double _freq
mixing frequency in hbar/mm