| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4DalitzDecayChannel class implementation
27 //
28 // Author: H.Kurashige, 30 May 1997
29 // --------------------------------------------------------------------
30
31 #include "G4DalitzDecayChannel.hh"
32
33 #include "G4DecayProducts.hh"
34 #include "G4LorentzRotation.hh"
35 #include "G4LorentzVector.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4PhaseSpaceDecayChannel.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4VDecayChannel.hh"
41 #include "Randomize.hh"
42
43 G4DalitzDecayChannel::G4DalitzDecayChannel(const G4String& theParentName, G4double theBR,
44 const G4String& theLeptonName,
45 const G4String& theAntiLeptonName)
46 : G4VDecayChannel("Dalitz Decay", 1)
47 {
48 // set names for daughter particles
49 SetParent(theParentName);
50 SetBR(theBR);
51 SetNumberOfDaughters(3);
52 G4String gammaName = "gamma";
53 SetDaughter(idGamma, gammaName);
54 SetDaughter(idLepton, theLeptonName);
55 SetDaughter(idAntiLepton, theAntiLeptonName);
56 }
57
58 G4DalitzDecayChannel& G4DalitzDecayChannel::operator=(const G4DalitzDecayChannel& right)
59 {
60 if (this != &right) {
61 kinematics_name = right.kinematics_name;
62 verboseLevel = right.verboseLevel;
63 rbranch = right.rbranch;
64
65 // copy parent name
66 parent_name = new G4String(*right.parent_name);
67
68 // clear daughters_name array
69 ClearDaughtersName();
70
71 // recreate array
72 numberOfDaughters = right.numberOfDaughters;
73 if (numberOfDaughters > 0) {
74 if (daughters_name != nullptr) ClearDaughtersName();
75 daughters_name = new G4String*[numberOfDaughters];
76 // copy daughters name
77 for (G4int index = 0; index < numberOfDaughters; ++index) {
78 daughters_name[index] = new G4String(*right.daughters_name[index]);
79 }
80 }
81 }
82 return *this;
83 }
84
85 G4DecayProducts* G4DalitzDecayChannel::DecayIt(G4double)
86 {
87 #ifdef G4VERBOSE
88 if (GetVerboseLevel() > 1) G4cout << "G4DalitzDecayChannel::DecayIt ";
89 #endif
90 CheckAndFillParent();
91 CheckAndFillDaughters();
92
93 // parent mass
94 G4double parentmass = G4MT_parent->GetPDGMass();
95
96 // create parent G4DynamicParticle at rest
97 G4ThreeVector dummy;
98 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
99
100 // daughters' mass
101 G4double leptonmass = G4MT_daughters[idLepton]->GetPDGMass();
102
103 // Generate t ( = std::exp(x):mass Square of (l+ l-) system)
104 G4double xmin = 2.0 * std::log(2.0 * leptonmass);
105 G4double xmax = 2.0 * std::log(parentmass);
106 G4double wmax = 1.5;
107 G4double x, w, ww, w1, w2, w3, t;
108 const std::size_t MAX_LOOP = 10000;
109 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
110 x = G4UniformRand() * (xmax - xmin) + xmin;
111 w = G4UniformRand() * wmax;
112 t = std::exp(x);
113 w1 = (1.0 - 4.0 * leptonmass * leptonmass / t);
114 if (w1 > 0.0) {
115 w2 = (1.0 + 2.0 * leptonmass * leptonmass / t);
116 w3 = (1.0 - t / parentmass / parentmass);
117 w3 = w3 * w3 * w3;
118 ww = w3 * w2 * std::sqrt(w1);
119 }
120 else {
121 ww = 0.0;
122 }
123 if (w <= ww) break;
124 }
125
126 // calculate gamma momentum
127 G4double Pgamma = G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, std::sqrt(t));
128 G4double costheta = 2. * G4UniformRand() - 1.0;
129 G4double sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
130 G4double phi = twopi * G4UniformRand() * rad;
131 G4ThreeVector gdirection(sintheta * std::cos(phi), sintheta * std::sin(phi), costheta);
132
133 // create G4DynamicParticle for gamma
134 auto gammaparticle = new G4DynamicParticle(G4MT_daughters[idGamma], gdirection, Pgamma);
135
136 // calculate beta of (l+ l-)system
137 G4double beta = Pgamma / (parentmass - Pgamma);
138
139 // calculate lepton momentum in the rest frame of (l+ l-)system
140 G4double Plepton = G4PhaseSpaceDecayChannel::Pmx(std::sqrt(t), leptonmass, leptonmass);
141 G4double Elepton = std::sqrt(Plepton * Plepton + leptonmass * leptonmass);
142 costheta = 2. * G4UniformRand() - 1.0;
143 sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
144 phi = twopi * G4UniformRand() * rad;
145 G4ThreeVector ldirection(sintheta * std::cos(phi), sintheta * std::sin(phi), costheta);
146 // create G4DynamicParticle for leptons in the rest frame of (l+ l-)system
147 auto leptonparticle =
148 new G4DynamicParticle(G4MT_daughters[idLepton], ldirection, Elepton - leptonmass);
149 auto antileptonparticle =
150 new G4DynamicParticle(G4MT_daughters[idAntiLepton], -1.0 * ldirection, Elepton - leptonmass);
151 // boost leptons in the rest frame of the parent
152 G4LorentzVector p4 = leptonparticle->Get4Momentum();
153 p4.boost(-1.0 * gdirection.x() * beta, -1.0 * gdirection.y() * beta,
154 -1.0 * gdirection.z() * beta);
155 leptonparticle->Set4Momentum(p4);
156 p4 = antileptonparticle->Get4Momentum();
157 p4.boost(-1.0 * gdirection.x() * beta, -1.0 * gdirection.y() * beta,
158 -1.0 * gdirection.z() * beta);
159 antileptonparticle->Set4Momentum(p4);
160
161 // create G4Decayproducts
162 auto products = new G4DecayProducts(*parentparticle);
163 delete parentparticle;
164 products->PushProducts(gammaparticle);
165 products->PushProducts(leptonparticle);
166 products->PushProducts(antileptonparticle);
167
168 #ifdef G4VERBOSE
169 if (GetVerboseLevel() > 1) {
170 G4cout << "G4DalitzDecayChannel::DecayIt ";
171 G4cout << " create decay products in rest frame " << G4endl;
172 products->DumpInfo();
173 }
174 #endif
175 return products;
176 }
177