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