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

#include <G4MuonDecayChannelWithSpin.hh>

+ Inheritance diagram for G4MuonDecayChannelWithSpin:

Public Member Functions

 G4MuonDecayChannelWithSpin (const G4String &theParentName, G4double theBR)
 
 ~G4MuonDecayChannelWithSpin () override=default
 
G4DecayProducts * DecayIt (G4double) override
 
- Public Member Functions inherited from G4MuonDecayChannel
 G4MuonDecayChannel (const G4String &parentName, G4double BR)
 
 ~G4MuonDecayChannel () override=default
 
- 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
 
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 Member Functions

 G4MuonDecayChannelWithSpin (const G4MuonDecayChannelWithSpin &)=default
 
G4MuonDecayChannelWithSpin & operator= (const G4MuonDecayChannelWithSpin &)
 
- Protected Member Functions inherited from G4MuonDecayChannel
 G4MuonDecayChannel ()=default
 
 G4MuonDecayChannel (const G4MuonDecayChannel &)=default
 
G4MuonDecayChannel & operator= (const G4MuonDecayChannel &)
 
- Protected Member Functions inherited from G4VDecayChannel
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChannel & operator= (const G4VDecayChannel &)
 
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
G4double DynamicalMass (G4double massPDG, G4double width, G4double maxDev=1.0) const
 

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 51 of file G4MuonDecayChannelWithSpin.hh.

Constructor & Destructor Documentation

◆ G4MuonDecayChannelWithSpin() [1/2]

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin ( const G4String & theParentName,
G4double theBR )

Definition at line 45 of file G4MuonDecayChannelWithSpin.cc.

47 : G4MuonDecayChannel(theParentName, theBR)
48{}
G4MuonDecayChannel::G4MuonDecayChannel
G4MuonDecayChannel()=default

Referenced by G4MuonDecayChannelWithSpin().

◆ ~G4MuonDecayChannelWithSpin()

G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin ( )
overridedefault

◆ G4MuonDecayChannelWithSpin() [2/2]

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin ( const G4MuonDecayChannelWithSpin & )
protecteddefault

Member Function Documentation

◆ DecayIt()

G4DecayProducts * G4MuonDecayChannelWithSpin::DecayIt ( G4double )
overridevirtual

Reimplemented from G4MuonDecayChannel.

Definition at line 78 of file G4MuonDecayChannelWithSpin.cc.

79{
80 // This version assumes V-A coupling with 1st order radiative correctons,
81 // the standard model Michel parameter values, but
82 // gives incorrect energy spectrum for neutrinos
83
84#ifdef G4VERBOSE
85 if (GetVerboseLevel() > 1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
86#endif
87
88 CheckAndFillParent();
89 CheckAndFillDaughters();
90
91 // parent mass
92 G4double parentmass = G4MT_parent->GetPDGMass();
93
94 G4double EMMU = parentmass;
95
96 // daughters'mass
97 G4double daughtermass[3];
98 // G4double sumofdaughtermass = 0.0;
99 for (G4int index = 0; index < 3; ++index) {
100 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
101 // sumofdaughtermass += daughtermass[index];
102 }
103
104 G4double EMASS = daughtermass[0];
105
106 // create parent G4DynamicParticle at rest
107 G4ThreeVector dummy;
108 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
109 // create G4Decayproducts
110 auto products = new G4DecayProducts(*parentparticle);
111 delete parentparticle;
112
113 // calculate electron energy
114
115 G4double michel_rho = 0.75; // Standard Model Michel rho
116 G4double michel_delta = 0.75; // Standard Model Michel delta
117 G4double michel_xsi = 1.00; // Standard Model Michel xsi
118 G4double michel_eta = 0.00; // Standard Model eta
119
120 G4double rndm, x, ctheta;
121
122 G4double FG;
123 G4double FG_max = 2.00;
124
125 G4double W_mue = (EMMU * EMMU + EMASS * EMASS) / (2. * EMMU);
126 G4double x0 = EMASS / W_mue;
127
128 G4double x0_squared = x0 * x0;
129
130 // ***************************************************
131 // x0 <= x <= 1. and -1 <= y <= 1
132 //
133 // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y
134 // ***************************************************
135
136 // ***** sampling F(x,y) directly (brute force) *****
137
138 const std::size_t MAX_LOOP = 10000;
139 for (std::size_t loop_count = 0; loop_count < MAX_LOOP; ++loop_count) {
140 // Sample the positron energy by sampling from F
141
142 rndm = G4UniformRand();
143
144 x = x0 + rndm * (1. - x0);
145
146 G4double x_squared = x * x;
147
148 G4double F_IS, F_AS, G_IS, G_AS;
149
150 F_IS = 1. / 6. * (-2. * x_squared + 3. * x - x0_squared);
151 F_AS = 1. / 6. * std::sqrt(x_squared - x0_squared) * (2. * x - 2. + std::sqrt(1. - x0_squared));
152
153 G_IS = 2. / 9. * (michel_rho - 0.75) * (4. * x_squared - 3. * x - x0_squared);
154 G_IS = G_IS + michel_eta * (1. - x) * x0;
155
156 G_AS = 3. * (michel_xsi - 1.) * (1. - x);
157 G_AS =
158 G_AS + 2. * (michel_xsi * michel_delta - 0.75) * (4. * x - 4. + std::sqrt(1. - x0_squared));
159 G_AS = 1. / 9. * std::sqrt(x_squared - x0_squared) * G_AS;
160
161 F_IS = F_IS + G_IS;
162 F_AS = F_AS + G_AS;
163
164 // *** Radiative Corrections ***
165
166 const G4double omega = std::log(EMMU / EMASS);
167 G4double R_IS = F_c(x, x0, omega);
168
169 G4double F = 6. * F_IS + R_IS / std::sqrt(x_squared - x0_squared);
170
171 // *** Radiative Corrections ***
172
173 G4double R_AS = F_theta(x, x0, omega);
174
175 rndm = G4UniformRand();
176
177 ctheta = 2. * rndm - 1.;
178
179 G4double G = 6. * F_AS - R_AS / std::sqrt(x_squared - x0_squared);
180
181 FG = std::sqrt(x_squared - x0_squared) * F * (1. + (G / F) * ctheta);
182
183 if (FG > FG_max) {
184 G4Exception("G4MuonDecayChannelWithSpin::DecayIt()", "PART113", JustWarning,
185 "Problem in Muon Decay: FG > FG_max");
186 FG_max = FG;
187 }
188
189 rndm = G4UniformRand();
190
191 if (FG >= rndm * FG_max) break;
192 }
193
194 G4double energy = x * W_mue;
195
196 rndm = G4UniformRand();
197
198 G4double phi = twopi * rndm;
199
200 if (energy < EMASS) energy = EMASS;
201
202 // Calculate daughter momentum
203 G4double daughtermomentum[3];
204
205 daughtermomentum[0] = std::sqrt(energy * energy - EMASS * EMASS);
206
207 G4double stheta = std::sqrt(1. - ctheta * ctheta);
208 G4double cphi = std::cos(phi);
209 G4double sphi = std::sin(phi);
210
211 // Coordinates of the decay positron with respect to the muon spin
212 G4double px = stheta * cphi;
213 G4double py = stheta * sphi;
214 G4double pz = ctheta;
215
216 G4ThreeVector direction0(px, py, pz);
217
218 direction0.rotateUz(parent_polarization);
219
220 auto daughterparticle0 =
221 new G4DynamicParticle(G4MT_daughters[0], daughtermomentum[0] * direction0);
222
223 products->PushProducts(daughterparticle0);
224
225 // daughter 1 ,2 (neutrinos)
226 // create neutrinos in the C.M frame of two neutrinos
227 G4double energy2 = parentmass - energy;
228 G4double vmass = std::sqrt((energy2 - daughtermomentum[0]) * (energy2 + daughtermomentum[0]));
229 G4double beta = -1.0 * daughtermomentum[0] / energy2;
230 G4double costhetan = 2. * G4UniformRand() - 1.0;
231 G4double sinthetan = std::sqrt((1.0 - costhetan) * (1.0 + costhetan));
232 G4double phin = twopi * G4UniformRand() * rad;
233 G4double sinphin = std::sin(phin);
234 G4double cosphin = std::cos(phin);
235
236 G4ThreeVector direction1(sinthetan * cosphin, sinthetan * sinphin, costhetan);
237 auto daughterparticle1 = new G4DynamicParticle(G4MT_daughters[1], direction1 * (vmass / 2.));
238 auto daughterparticle2 =
239 new G4DynamicParticle(G4MT_daughters[2], direction1 * (-1.0 * vmass / 2.));
240
241 // boost to the muon rest frame
242 G4LorentzVector p4;
243 p4 = daughterparticle1->Get4Momentum();
244 p4.boost(direction0.x() * beta, direction0.y() * beta, direction0.z() * beta);
245 daughterparticle1->Set4Momentum(p4);
246 p4 = daughterparticle2->Get4Momentum();
247 p4.boost(direction0.x() * beta, direction0.y() * beta, direction0.z() * beta);
248 daughterparticle2->Set4Momentum(p4);
249 products->PushProducts(daughterparticle1);
250 products->PushProducts(daughterparticle2);
251 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
252 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
253
254 // output message
255#ifdef G4VERBOSE
256 if (GetVerboseLevel() > 1) {
257 G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
258 G4cout << " create decay products in rest frame " << G4endl;
259 G4double TT = daughterparticle0->GetTotalEnergy() + daughterparticle1->GetTotalEnergy()
260 + daughterparticle2->GetTotalEnergy();
261 G4cout << "e " << daughterparticle0->GetTotalEnergy() / MeV << G4endl;
262 G4cout << "nu1" << daughterparticle1->GetTotalEnergy() / MeV << G4endl;
263 G4cout << "nu2" << daughterparticle2->GetTotalEnergy() / MeV << G4endl;
264 G4cout << "total" << (TT - parentmass) / keV << G4endl;
265 if (GetVerboseLevel() > 2) {
266 products->DumpInfo();
267 }
268 }
269#endif
270
271 return products;
272}
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition LorentzVector.cc:55
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition G4VDecayChannel.hh:164
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition G4VDecayChannel.hh:285
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition G4VDecayChannel.hh:163
G4VDecayChannel::parent_polarization
G4ThreeVector parent_polarization
Definition G4VDecayChannel.hh:156
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185

◆ operator=()

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

Definition at line 51 of file G4MuonDecayChannelWithSpin.cc.

52{
53 if (this != &right) {
54 kinematics_name = right.kinematics_name;
55 verboseLevel = right.verboseLevel;
56 rbranch = right.rbranch;
57
58 // copy parent name
59 delete parent_name;
60 parent_name = new G4String(*right.parent_name);
61
62 // clear daughters_name array
63 ClearDaughtersName();
64
65 // recreate array
66 numberOfDaughters = right.numberOfDaughters;
67 if (numberOfDaughters > 0) {
68 daughters_name = new G4String*[numberOfDaughters];
69 // copy daughters name
70 for (G4int index = 0; index < numberOfDaughters; ++index) {
71 daughters_name[index] = new G4String(*right.daughters_name[index]);
72 }
73 }
74 }
75 return *this;
76}
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition G4VDecayChannel.hh:150
The documentation for this class was generated from the following files: