Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Protected Types | Protected Member Functions | List of all members
G4KL3DecayChannel Class Reference

#include <G4KL3DecayChannel.hh>

+ Inheritance diagram for G4KL3DecayChannel:

Public Member Functions

 G4KL3DecayChannel (const G4String &theParentName, G4double theBR, const G4String &thePionName, const G4String &theLeptonName, const G4String &theNutrinoName)
 
virtual ~G4KL3DecayChannel ()
 
virtual G4DecayProducts * DecayIt (G4double)
 
void SetDalitzParameter (G4double aLambda, G4double aXi)
 
G4double GetDalitzParameterLambda () const
 
G4double GetDalitzParameterXi () const
 
- Public Member Functions inherited from G4VDecayChannel
 G4VDecayChannel (const G4String &aName, G4int Verbose=1)
 
 G4VDecayChannel (const G4String &aName, const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="", const G4String &theDaughterName4="", const G4String &theDaughterName5="")
 
virtual ~G4VDecayChannel ()
 
G4bool operator== (const G4VDecayChannel &r) const
 
G4bool operator!= (const G4VDecayChannel &r) const
 
G4bool operator< (const G4VDecayChannel &right) const
 
virtual G4DecayProducts * DecayIt (G4double parentMass=-1.0)=0
 
const G4String & GetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinition * GetParent ()
 
G4ParticleDefinition * GetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4String & GetParentName () const
 
const G4String & GetDaughterName (G4int anIndex) const
 
G4double GetParentMass () const
 
G4double GetDaughterMass (G4int anIndex) const
 
void SetParent (const G4ParticleDefinition *particle_type)
 
void SetParent (const G4String &particle_name)
 
void SetBR (G4double value)
 
void SetNumberOfDaughters (G4int value)
 
void SetDaughter (G4int anIndex, const G4ParticleDefinition *particle_type)
 
void SetDaughter (G4int anIndex, const G4String &particle_name)
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
void DumpInfo ()
 
G4double GetRangeMass () const
 
void SetRangeMass (G4double val)
 
virtual G4bool IsOKWithParentMass (G4double parentMass)
 
void SetPolarization (const G4ThreeVector &)
 
const G4ThreeVector & GetPolarization () const
 

Protected Types

enum  { idPi =0 , idLepton =1 , idNutrino =2 }
 

Protected Member Functions

 G4KL3DecayChannel (const G4KL3DecayChannel &)
 
G4KL3DecayChannel & operator= (const G4KL3DecayChannel &)
 
void PhaseSpace (G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
 
G4double DalitzDensity (G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
 
- Protected Member Functions inherited from G4VDecayChannel
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
G4double DynamicalMass (G4double massPDG, G4double width, G4double maxDev=1.0) const
 
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChannel & operator= (const G4VDecayChannel &)
 

Additional Inherited Members

- Protected Attributes inherited from G4VDecayChannel
G4String kinematics_name = ""
 
G4double rbranch = 0.0
 
G4String * parent_name = nullptr
 
G4String ** daughters_name = nullptr
 
G4double rangeMass = 2.5
 
G4ThreeVector parent_polarization
 
G4ParticleTable * particletable = nullptr
 
G4ParticleDefinition * G4MT_parent = nullptr
 
G4ParticleDefinition ** G4MT_daughters = nullptr
 
G4double G4MT_parent_mass = 0.0
 
G4double * G4MT_daughters_mass = nullptr
 
G4double * G4MT_daughters_width = nullptr
 
G4Mutex daughtersMutex
 
G4Mutex parentMutex
 
G4int numberOfDaughters = 0
 
G4int verboseLevel = 1
 
- Static Protected Attributes inherited from G4VDecayChannel
static const G4String noName = " "
 

Detailed Description

Definition at line 37 of file G4KL3DecayChannel.hh.

Member Enumeration Documentation

◆ anonymous enum

anonymous enum
protected
Enumerator
idPi 
idLepton 
idNutrino 

Definition at line 61 of file G4KL3DecayChannel.hh.

61{ idPi=0, idLepton=1, idNutrino=2 };

Constructor & Destructor Documentation

◆ G4KL3DecayChannel() [1/2]

G4KL3DecayChannel::G4KL3DecayChannel ( const G4String &  theParentName,
G4double  theBR,
const G4String &  thePionName,
const G4String &  theLeptonName,
const G4String &  theNutrinoName 
)

Definition at line 47 of file G4KL3DecayChannel.cc.

53 : G4VDecayChannel("KL3 Decay", theParentName, theBR, 3,
54 thePionName, theLeptonName, theNutrinoName)
55{
56 static const G4String K_plus("kaon+");
57 static const G4String K_minus("kaon-");
58 static const G4String K_L("kaon0L");
59 static const G4String Mu_plus("mu+");
60 static const G4String Mu_minus("mu-");
61 static const G4String E_plus("e+");
62 static const G4String E_minus("e-");
63
64 // check modes
65 if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
66 ((theParentName == K_minus)&&(theLeptonName == E_minus)) )
67 {
68 // K+- (Ke3)
69 pLambda = 0.0286;
70 pXi0 = -0.35;
71 }
72 else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
73 ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) )
74 {
75 // K+- (Kmu3)
76 pLambda = 0.033;
77 pXi0 = -0.35;
78 }
79 else if ( (theParentName == K_L) &&
80 ((theLeptonName == E_plus) || (theLeptonName == E_minus)) )
81 {
82 // K0L (Ke3)
83 pLambda = 0.0300;
84 pXi0 = -0.11;
85 }
86 else if ( (theParentName == K_L) &&
87 ((theLeptonName == Mu_plus) || (theLeptonName == Mu_minus)) )
88 {
89 // K0L (Kmu3)
90 pLambda = 0.034;
91 pXi0 = -0.11;
92 }
93 else
94 {
95#ifdef G4VERBOSE
96 if (GetVerboseLevel()>2)
97 {
98 G4cout << "G4KL3DecayChannel:: constructor :";
99 G4cout << "illegal arguments " << G4endl;;
100 DumpInfo();
101 }
102#endif
103 // set values for K0L (Ke3) temporarily
104 pLambda = 0.0300;
105 pXi0 = -0.11;
106 }
107}
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:310

Referenced by G4KL3DecayChannel().

◆ ~G4KL3DecayChannel()

G4KL3DecayChannel::~G4KL3DecayChannel ( )
virtual

Definition at line 109 of file G4KL3DecayChannel.cc.

110{
111}

◆ G4KL3DecayChannel() [2/2]

G4KL3DecayChannel::G4KL3DecayChannel ( const G4KL3DecayChannel &  right)
protected

Definition at line 113 of file G4KL3DecayChannel.cc.

114 : G4VDecayChannel(right),
115 pLambda(right.pLambda),
116 pXi0(right.pXi0)
117{
118}

Member Function Documentation

◆ DalitzDensity()

G4double G4KL3DecayChannel::DalitzDensity ( G4double  parentmass,
G4double  Epi,
G4double  El,
G4double  Enu,
G4double  massPi,
G4double  massL,
G4double  massNu 
)
protected

Definition at line 333 of file G4KL3DecayChannel.cc.

336{
337 // KL3 decay - Dalitz Plot Density, see Chounet et al Phys. Rep. 4, 201
338 // Arguments
339 // Epi: kinetic enregy of pion
340 // El: kinetic enregy of lepton (e or mu)
341 // Enu: kinetic energy of nutrino
342 // Constants
343 // pLambda : linear energy dependence of f+
344 // pXi0 : = f+(0)/f-
345 // pNorm : normalization factor
346 // Variables
347 // Epi: total energy of pion
348 // El: total energy of lepton (e or mu)
349 // Enu: total energy of nutrino
350
351 // calculate total energy
352 Epi = Epi + massPi;
353 El = El + massL;
354 Enu = Enu + massNu;
355
356 G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
357 G4double E = Epi_max - Epi;
358 G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
359
360 G4double F = 1.0 + pLambda*q2/massPi/massPi;
361 G4double Fmax = 1.0;
362 if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
363
364 G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
365
366 G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
367 G4double coeffB = massL*massL*(Enu-E/2.0);
368 G4double coeffC = massL*massL*E/4.0;
369
370 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
371
372 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
373
374#ifdef G4VERBOSE
375 if (GetVerboseLevel()>2)
376 {
377 G4cout << "G4KL3DecayChannel::DalitzDensity " <<G4endl;
378 G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl;
379 G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl;
380 G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl;
381 G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl;
382 G4cout << " A :" << coeffA << " B :" << coeffB << " C :"<< coeffC
383 << G4endl;
384 G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl;
385 }
386#endif
387 return (Rho/RhoMax);
388}
G4double
double G4double
Definition: G4Types.hh:83

Referenced by DecayIt().

◆ DecayIt()

G4DecayProducts * G4KL3DecayChannel::DecayIt ( G4double  )
virtual

Implements G4VDecayChannel.

Definition at line 152 of file G4KL3DecayChannel.cc.

153{
154 // this version neglects muon polarization
155 // assumes the pure V-A coupling
156 // gives incorrect energy spectrum for Neutrinos
157#ifdef G4VERBOSE
158 if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
159#endif
160
161 // fill parent particle and its mass
162 CheckAndFillParent();
163 G4double massK = G4MT_parent->GetPDGMass();
164
165 // fill daughter particles and their mass
166 CheckAndFillDaughters();
167 G4double daughterM[3];
168 daughterM[idPi] = G4MT_daughters[idPi]->GetPDGMass();
169 daughterM[idLepton] = G4MT_daughters[idLepton]->GetPDGMass();
170 daughterM[idNutrino] = G4MT_daughters[idNutrino]->GetPDGMass();
171
172 // determine momentum/energy of daughters according to DalitzDensity
173 G4double daughterP[3], daughterE[3];
174 G4double w;
175 G4double r;
176 const size_t MAX_LOOP = 10000;
177 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
178 {
179 r = G4UniformRand();
180 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
181 w = DalitzDensity(massK, daughterE[idPi], daughterE[idLepton],
182 daughterE[idNutrino], daughterM[idPi],
183 daughterM[idLepton], daughterM[idNutrino]);
184 if ( r <= w) break;
185 }
186
187 // output message
188#ifdef G4VERBOSE
189 if (GetVerboseLevel()>1)
190 {
191 G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]"
192 << G4endl;
193 G4cout << *daughters_name[1] << ":" << daughterP[1]/GeV << "[GeV/c]"
194 << G4endl;
195 G4cout << *daughters_name[2] << ":" << daughterP[2]/GeV << "[GeV/c]"
196 << G4endl;
197 }
198#endif
199
200 // create parent G4DynamicParticle at rest
201 G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0);
202 G4DynamicParticle* parentparticle
203 = new G4DynamicParticle( G4MT_parent, *direction, 0.0);
204 delete direction;
205
206 // create G4Decayproducts
207 G4DecayProducts *products
208 = new G4DecayProducts(*parentparticle);
209 delete parentparticle;
210
211 // create daughter G4DynamicParticle
212 G4double costheta, sintheta, phi, sinphi, cosphi;
213 G4double costhetan, sinthetan, phin, sinphin, cosphin;
214
215 // pion
216 costheta = 2.*G4UniformRand()-1.0;
217 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
218 phi = twopi*G4UniformRand()*rad;
219 sinphi = std::sin(phi);
220 cosphi = std::cos(phi);
221 direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
222 G4ThreeVector momentum0 = (*direction)*daughterP[0];
223 G4DynamicParticle * daughterparticle
224 = new G4DynamicParticle(G4MT_daughters[0], momentum0);
225 products->PushProducts(daughterparticle);
226
227 // neutrino
228 costhetan = (daughterP[1]*daughterP[1]
229 - daughterP[2]*daughterP[2]
230 - daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
231 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
232 phin = twopi*G4UniformRand()*rad;
233 sinphin = std::sin(phin);
234 cosphin = std::cos(phin);
235 direction->setX( sinthetan*cosphin*costheta*cosphi
236 - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
237 direction->setY( sinthetan*cosphin*costheta*sinphi
238 + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
239 direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
240
241 G4ThreeVector momentum2 = (*direction)*daughterP[2];
242 daughterparticle = new G4DynamicParticle( G4MT_daughters[2], momentum2);
243 products->PushProducts(daughterparticle);
244
245 // lepton
246 G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
247 daughterparticle = new G4DynamicParticle( G4MT_daughters[1], momentum1);
248 products->PushProducts(daughterparticle);
249
250#ifdef G4VERBOSE
251 if (GetVerboseLevel()>1)
252 {
253 G4cout << "G4KL3DecayChannel::DecayIt ";
254 G4cout << " create decay products in rest frame " <<G4endl;
255 G4cout << " decay products address=" << products << G4endl;
256 products->DumpInfo();
257 }
258#endif
259 delete direction;
260 return products;
261}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:334
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:151
G4KL3DecayChannel::PhaseSpace
void PhaseSpace(G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
Definition: G4KL3DecayChannel.cc:263
G4KL3DecayChannel::DalitzDensity
G4double DalitzDensity(G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
Definition: G4KL3DecayChannel.cc:333
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:173
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:158
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition: G4VDecayChannel.hh:172
G4VDecayChannel::CheckAndFillDaughters
void CheckAndFillDaughters()
Definition: G4VDecayChannel.hh:337

◆ GetDalitzParameterLambda()

G4double G4KL3DecayChannel::GetDalitzParameterLambda ( ) const
inline

Definition at line 107 of file G4KL3DecayChannel.hh.

108{
109 return pLambda;
110}

◆ GetDalitzParameterXi()

G4double G4KL3DecayChannel::GetDalitzParameterXi ( ) const
inline

Definition at line 113 of file G4KL3DecayChannel.hh.

114{
115 return pXi0;
116}

◆ operator=()

G4KL3DecayChannel & G4KL3DecayChannel::operator= ( const G4KL3DecayChannel &  right)
protected

Definition at line 120 of file G4KL3DecayChannel.cc.

121{
122 if (this != &right)
123 {
124 kinematics_name = right.kinematics_name;
125 verboseLevel = right.verboseLevel;
126 rbranch = right.rbranch;
127
128 // copy parent name
129 parent_name = new G4String(*right.parent_name);
130
131 // clear daughters_name array
132 ClearDaughtersName();
133
134 // recreate array
135 numberOfDaughters = right.numberOfDaughters;
136 if ( numberOfDaughters >0 )
137 {
138 if (daughters_name !=0) ClearDaughtersName();
139 daughters_name = new G4String*[numberOfDaughters];
140 //copy daughters name
141 for (G4int index=0; index<numberOfDaughters; ++index)
142 {
143 daughters_name[index] = new G4String(*right.daughters_name[index]);
144 }
145 }
146 pLambda = right.pLambda;
147 pXi0 = right.pXi0;
148 }
149 return *this;
150}
G4int
int G4int
Definition: G4Types.hh:85
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:156
G4VDecayChannel::kinematics_name
G4String kinematics_name
Definition: G4VDecayChannel.hh:152

◆ PhaseSpace()

void G4KL3DecayChannel::PhaseSpace ( G4double  Mparent,
const G4double *  Mdaughter,
G4double *  Edaughter,
G4double *  Pdaughter 
)
protected

Definition at line 263 of file G4KL3DecayChannel.cc.

267{
268 // Algorithm in this code was originally written in GDECA3 in GEANT3
269
270 // sum of daughters'mass
271 G4double sumofdaughtermass = 0.0;
272 G4int index;
273 const G4int N_DAUGHTER=3;
274
275 for (index=0; index<N_DAUGHTER; ++index)
276 {
277 sumofdaughtermass += M[index];
278 }
279
280 // calculate daughter momentum. Generate two
281 G4double rd1, rd2, rd;
282 G4double momentummax=0.0, momentumsum = 0.0;
283 G4double energy;
284 const size_t MAX_LOOP=10000;
285 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
286 {
287 rd1 = G4UniformRand();
288 rd2 = G4UniformRand();
289 if (rd2 > rd1)
290 {
291 rd = rd1;
292 rd1 = rd2;
293 rd2 = rd;
294 }
295 momentummax = 0.0;
296 momentumsum = 0.0;
297 // daughter 0
298 energy = rd2*(parentM - sumofdaughtermass);
299 P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
300 E[0] = energy;
301 if ( P[0] >momentummax )momentummax = P[0];
302 momentumsum += P[0];
303 // daughter 1
304 energy = (1.-rd1)*(parentM - sumofdaughtermass);
305 P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
306 E[1] = energy;
307 if ( P[1] >momentummax )momentummax = P[1];
308 momentumsum += P[1];
309 // daughter 2
310 energy = (rd1-rd2)*(parentM - sumofdaughtermass);
311 P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
312 E[2] = energy;
313 if ( P[2] >momentummax )momentummax = P[2];
314 momentumsum += P[2];
315 if (momentummax <= momentumsum - momentummax ) break;
316 }
317#ifdef G4VERBOSE
318 if (GetVerboseLevel()>2)
319 {
320 G4cout << "G4KL3DecayChannel::PhaseSpace ";
321 G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl;
322 for (index=0; index<3; ++index)
323 {
324 G4cout << index << " : " << M[index]/GeV << "GeV/c/c ";
325 G4cout << " : " << E[index]/GeV << "GeV ";
326 G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl;
327 }
328 }
329#endif
330}
M
#define M(row, col)
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:158

Referenced by DecayIt().

◆ SetDalitzParameter()

void G4KL3DecayChannel::SetDalitzParameter ( G4double  aLambda,
G4double  aXi 
)
inline

Definition at line 100 of file G4KL3DecayChannel.hh.

101{
102 pLambda = aLambda;
103 pXi0 = aXi;
104}
The documentation for this class was generated from the following files: