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