EvtGen 2.2.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
Loading...
Searching...
No Matches
EvtIncoherentMixing.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/EvtId.hh"
24#include "EvtGenBase/EvtPDL.hh"
26
27#include <stdlib.h>
28
29//-----------------------------------------------------------------------------
30// Implementation file for class : EvtIncoherentMixing
31//
32// 2003-10-09 : Patrick Robbe
33//-----------------------------------------------------------------------------
34
39double EvtIncoherentMixing::m_deltamd = 0.502e12;
40// dGamma_s corresponds to DeltaGamma / Gamma = 10 %
41double EvtIncoherentMixing::m_dGammas = 6.852e10;
42double EvtIncoherentMixing::m_deltams = 20.e12;
43
44//=============================================================================
45// Standard constructor, initializes variables
46//=============================================================================
48{
49 m_doB0Mixing = false;
50 m_doBsMixing = false;
51 m_dGammad = 0.;
52 // dGammas corresponds to DeltaGamma / Gamma = 10 %
53 m_dGammas = 6.852e10;
54 m_deltamd = 0.502e12;
55 m_deltams = 20.e12;
56 m_enableFlip = false;
57}
58//=============================================================================
59void EvtIncoherentMixing::incoherentB0Mix( const EvtId id, double& t, int& mix )
60{
61 static const EvtId B0 = EvtPDL::getId( "B0" );
62 static const EvtId B0B = EvtPDL::getId( "anti-B0" );
63
64 if ( ( B0 != id ) && ( B0B != id ) ) {
65 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
66 << "Bad configuration in incoherentB0Mix" << std::endl;
67 ::abort();
68 }
69
70 double x = getdeltamd() * EvtPDL::getctau( B0 ) / EvtConst::c;
71
72 double y = getdGammad() * ( EvtPDL::getctau( B0 ) / EvtConst::c ) / 2.;
73
74 double fac = 1.; // No CP violation
75
76 double mixprob = ( x * x + y * y ) /
77 ( x * x + y * y + ( 1. / fac ) * ( 2. + x * x - y * y ) );
78
79 int mixsign;
80
81 // decide if state is mixed
82 mixsign = ( mixprob > EvtRandom::Flat( 0., 1. ) ) ? -1 : 1;
83
84 double prob;
85
86 do {
87 t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( B0 ) / ( 1. - y );
88 prob = ( 1. + exp( -2. * y * t / EvtPDL::getctau( B0 ) ) +
89 mixsign * 2. * exp( -y * t / EvtPDL::getctau( B0 ) ) *
90 cos( getdeltamd() * t / EvtConst::c ) ) /
91 2.;
92 } while ( prob < 2. * EvtRandom::Flat() );
93
94 mix = 0;
95 if ( mixsign == -1 )
96 mix = 1;
97
98 return;
99}
100// ============================================================================
101void EvtIncoherentMixing::incoherentBsMix( const EvtId id, double& t, int& mix )
102{
103 static const EvtId BS = EvtPDL::getId( "B_s0" );
104 static const EvtId BSB = EvtPDL::getId( "anti-B_s0" );
105
106 if ( ( BS != id ) && ( BSB != id ) ) {
107 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
108 << "Bad configuration in incoherentBsMix" << std::endl;
109 ::abort();
110 }
111
112 double x = getdeltams() * EvtPDL::getctau( BS ) / EvtConst::c;
113
114 double y = getdGammas() * ( EvtPDL::getctau( BS ) / EvtConst::c ) / 2.;
115
116 double fac = 1.; // No CP violation
117
118 double mixprob = ( x * x + y * y ) /
119 ( x * x + y * y + ( 1. / fac ) * ( 2. + x * x - y * y ) );
120
121 int mixsign;
122
123 // decide if state is mixed
124 mixsign = ( mixprob > EvtRandom::Flat( 0., 1. ) ) ? -1 : 1;
125
126 double prob;
127
128 do {
129 t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( BS ) / ( 1. - y );
130 prob = ( 1. + exp( -2. * y * t / EvtPDL::getctau( BS ) ) +
131 mixsign * 2. * exp( -y * t / EvtPDL::getctau( BS ) ) *
132 cos( getdeltams() * t / EvtConst::c ) ) /
133 2.;
134 } while ( prob < 2. * EvtRandom::Flat() );
135
136 mix = 0;
137 if ( mixsign == -1 )
138 mix = 1;
139
140 return;
141}
142
143// ========================================================================
145{
146 if ( !( p->getParent() ) )
147 return false;
148
149 static const EvtId BS0 = EvtPDL::getId( "B_s0" );
150 static const EvtId BSB = EvtPDL::getId( "anti-B_s0" );
151
152 if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) )
153 return false;
154
155 if ( ( p->getParent()->getId() == BS0 ) || ( p->getParent()->getId() == BSB ) )
156 return true;
157
158 return false;
159}
160
161// ========================================================================
163{
164 if ( !( p->getParent() ) )
165 return false;
166
167 static const EvtId B0 = EvtPDL::getId( "B0" );
168 static const EvtId B0B = EvtPDL::getId( "anti-B0" );
169
170 if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) )
171 return false;
172
173 if ( ( p->getParent()->getId() == B0 ) || ( p->getParent()->getId() == B0B ) )
174 return true;
175
176 return false;
177}
178//============================================================================
179// Return the tag of the event (ie the anti-flavour of the produced
180// B meson). Flip the flavour of the event with probB probability
181//============================================================================
182void EvtIncoherentMixing::OtherB( EvtParticle* p, double& t, EvtId& otherb,
183 double probB )
184{
185 //if(p->getId() == B0 || p->getId() == B0B)
186 //added by liming Zhang
187 enableFlip();
188 if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
189 p->getParent()->setLifetime();
190 t = p->getParent()->getLifetime();
191 } else {
192 p->setLifetime();
193 t = p->getLifetime();
194 }
195
196 if ( flipIsEnabled() ) {
197 //std::cout << " liming << flipIsEnabled " << std::endl;
198 // Flip the flavour of the particle with probability probB
199 bool isFlipped = ( EvtRandom::Flat( 0., 1. ) < probB );
200
201 if ( isFlipped ) {
202 if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
203 p->getParent()->setId(
205 p->setId( EvtPDL::chargeConj( p->getId() ) );
206 } else {
207 p->setId( EvtPDL::chargeConj( p->getId() ) );
208 }
209 }
210 }
211
212 if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
213 // if B has mixed, tag flavour is charge conjugate of parent of B-meson
214 otherb = EvtPDL::chargeConj( p->getParent()->getId() );
215 } else {
216 // else it is opposite flavour than this B hadron
217 otherb = EvtPDL::chargeConj( p->getId() );
218 }
219
220 return;
221}
222//============================================================================
223// Return the tag of the event (ie the anti-flavour of the produced
224// B meson). No flip
225//============================================================================
226void EvtIncoherentMixing::OtherB( EvtParticle* p, double& t, EvtId& otherb )
227{
228 if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
229 p->getParent()->setLifetime();
230 t = p->getParent()->getLifetime();
231 } else {
232 p->setLifetime();
233 t = p->getLifetime();
234 }
235
236 if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
237 // if B has mixed, tag flavour is charge conjugate of parent of B-meson
238 otherb = EvtPDL::chargeConj( p->getParent()->getId() );
239 } else {
240 // else it is opposite flavour than this B hadron
241 otherb = EvtPDL::chargeConj( p->getId() );
242 }
243
244 return;
245}
246
247// activate or desactivate the Bs mixing
256
257// activate or desactivate the B0 mixing
266
267// is mixing activated ?
276
277// set values for the mixing
279{
280 m_dGammad = value;
281}
283{
284 m_deltamd = value;
285}
287{
288 m_dGammas = value;
289}
291{
292 m_deltams = value;
293}
294
295// get parameters for mixing
297{
298 return m_dGammad;
299}
301{
302 return m_deltamd;
303}
305{
306 return m_dGammas;
307}
309{
310 return m_deltams;
311}
312
322{
323 m_enableFlip = false;
324}
EvtComplex exp(const EvtComplex &c)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
Definition EvtReport.cpp:32
@ EVTGEN_ERROR
Definition EvtReport.hh:49
static const double c
Definition EvtConst.hh:30
Definition EvtId.hh:27
static void setdeltams(double value)
static void setdGammad(double value)
static bool isB0Mixed(EvtParticle *)
static void setdGammas(double value)
static void setdeltamd(double value)
static bool isBsMixed(EvtParticle *)
static void OtherB(EvtParticle *p, double &t, EvtId &otherb, double probB)
static void incoherentBsMix(const EvtId id, double &t, int &mix)
EvtIncoherentMixing()
Standard constructor.
static void incoherentB0Mix(const EvtId id, double &t, int &mix)
static EvtId chargeConj(EvtId id)
Definition EvtPDL.cpp:201
static double getctau(EvtId i)
Definition EvtPDL.cpp:351
static EvtId getId(const std::string &name)
Definition EvtPDL.cpp:283
EvtId getId() const
void setLifetime(double tau)
double getLifetime() const
void setId(EvtId id)
EvtParticle * getParent() const
static double Flat()
Definition EvtRandom.cpp:95