| 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 // G4KL3DecayChannel class implementation 26 // G4KL3DecayChannel class implementation
27 // 27 //
28 // Author: H.Kurashige, 30 May 1997 << 28 // Author: H.Kurashige, 30 May 1997
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4KL3DecayChannel.hh" <<
32 <<
33 #include "G4DecayProducts.hh" <<
34 #include "G4LorentzRotation.hh" <<
35 #include "G4LorentzVector.hh" <<
36 #include "G4ParticleDefinition.hh" 31 #include "G4ParticleDefinition.hh"
37 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
>> 34 #include "G4DecayProducts.hh"
39 #include "G4VDecayChannel.hh" 35 #include "G4VDecayChannel.hh"
>> 36 #include "G4KL3DecayChannel.hh"
40 #include "Randomize.hh" 37 #include "Randomize.hh"
>> 38 #include "G4LorentzVector.hh"
>> 39 #include "G4LorentzRotation.hh"
41 40
42 G4KL3DecayChannel::G4KL3DecayChannel(const G4S << 41 G4KL3DecayChannel::G4KL3DecayChannel()
43 const G4S << 42 : G4VDecayChannel()
44 const G4S << 43 {
45 : G4VDecayChannel("KL3 Decay", theParentName << 44 }
46 theNutrinoName) << 45
>> 46
>> 47 G4KL3DecayChannel::
>> 48 G4KL3DecayChannel(const G4String& theParentName,
>> 49 G4double theBR,
>> 50 const G4String& thePionName,
>> 51 const G4String& theLeptonName,
>> 52 const G4String& theNutrinoName)
>> 53 : G4VDecayChannel("KL3 Decay", theParentName, theBR, 3,
>> 54 thePionName, theLeptonName, theNutrinoName)
47 { 55 {
48 static const G4String K_plus("kaon+"); 56 static const G4String K_plus("kaon+");
49 static const G4String K_minus("kaon-"); 57 static const G4String K_minus("kaon-");
50 static const G4String K_L("kaon0L"); 58 static const G4String K_L("kaon0L");
51 static const G4String Mu_plus("mu+"); 59 static const G4String Mu_plus("mu+");
52 static const G4String Mu_minus("mu-"); 60 static const G4String Mu_minus("mu-");
53 static const G4String E_plus("e+"); 61 static const G4String E_plus("e+");
54 static const G4String E_minus("e-"); 62 static const G4String E_minus("e-");
55 << 63
56 // check modes 64 // check modes
57 if (((theParentName == K_plus) && (theLepton << 65 if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
58 || ((theParentName == K_minus) && (theLe << 66 ((theParentName == K_minus)&&(theLeptonName == E_minus)) )
59 { 67 {
60 // K+- (Ke3) 68 // K+- (Ke3)
61 pLambda = 0.0286; 69 pLambda = 0.0286;
62 pXi0 = -0.35; << 70 pXi0 = -0.35;
63 } 71 }
64 else if (((theParentName == K_plus) && (theL << 72 else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
65 || ((theParentName == K_minus) && ( << 73 ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) )
66 { 74 {
67 // K+- (Kmu3) 75 // K+- (Kmu3)
68 pLambda = 0.033; 76 pLambda = 0.033;
69 pXi0 = -0.35; << 77 pXi0 = -0.35;
70 } 78 }
71 else if ((theParentName == K_L) && ((theLept << 79 else if ( (theParentName == K_L) &&
>> 80 ((theLeptonName == E_plus) || (theLeptonName == E_minus)) )
>> 81 {
72 // K0L (Ke3) 82 // K0L (Ke3)
73 pLambda = 0.0300; 83 pLambda = 0.0300;
74 pXi0 = -0.11; << 84 pXi0 = -0.11;
75 } 85 }
76 else if ((theParentName == K_L) && ((theLept << 86 else if ( (theParentName == K_L) &&
>> 87 ((theLeptonName == Mu_plus) || (theLeptonName == Mu_minus)) )
>> 88 {
77 // K0L (Kmu3) 89 // K0L (Kmu3)
78 pLambda = 0.034; 90 pLambda = 0.034;
79 pXi0 = -0.11; << 91 pXi0 = -0.11;
80 } 92 }
81 else { << 93 else
>> 94 {
82 #ifdef G4VERBOSE 95 #ifdef G4VERBOSE
83 if (GetVerboseLevel() > 2) { << 96 if (GetVerboseLevel()>2)
>> 97 {
84 G4cout << "G4KL3DecayChannel:: construct 98 G4cout << "G4KL3DecayChannel:: constructor :";
85 G4cout << "illegal arguments " << G4endl << 99 G4cout << "illegal arguments " << G4endl;;
86 ; <<
87 DumpInfo(); 100 DumpInfo();
88 } 101 }
89 #endif 102 #endif
90 // set values for K0L (Ke3) temporarily 103 // set values for K0L (Ke3) temporarily
91 pLambda = 0.0300; 104 pLambda = 0.0300;
92 pXi0 = -0.11; << 105 pXi0 = -0.11;
93 } 106 }
94 } 107 }
95 108
>> 109 G4KL3DecayChannel::~G4KL3DecayChannel()
>> 110 {
>> 111 }
>> 112
>> 113 G4KL3DecayChannel::G4KL3DecayChannel(const G4KL3DecayChannel& right)
>> 114 : G4VDecayChannel(right),
>> 115 pLambda(right.pLambda),
>> 116 pXi0(right.pXi0)
>> 117 {
>> 118 }
>> 119
96 G4KL3DecayChannel& G4KL3DecayChannel::operator 120 G4KL3DecayChannel& G4KL3DecayChannel::operator=(const G4KL3DecayChannel& right)
97 { 121 {
98 if (this != &right) { << 122 if (this != &right)
>> 123 {
99 kinematics_name = right.kinematics_name; 124 kinematics_name = right.kinematics_name;
100 verboseLevel = right.verboseLevel; 125 verboseLevel = right.verboseLevel;
101 rbranch = right.rbranch; 126 rbranch = right.rbranch;
102 127
103 // copy parent name 128 // copy parent name
104 parent_name = new G4String(*right.parent_n 129 parent_name = new G4String(*right.parent_name);
105 130
106 // clear daughters_name array 131 // clear daughters_name array
107 ClearDaughtersName(); 132 ClearDaughtersName();
108 133
109 // recreate array 134 // recreate array
110 numberOfDaughters = right.numberOfDaughter 135 numberOfDaughters = right.numberOfDaughters;
111 if (numberOfDaughters > 0) { << 136 if ( numberOfDaughters >0 )
112 if (daughters_name != nullptr) ClearDaug << 137 {
>> 138 if (daughters_name !=0) ClearDaughtersName();
113 daughters_name = new G4String*[numberOfD 139 daughters_name = new G4String*[numberOfDaughters];
114 // copy daughters name << 140 //copy daughters name
115 for (G4int index = 0; index < numberOfDa << 141 for (G4int index=0; index<numberOfDaughters; ++index)
>> 142 {
116 daughters_name[index] = new G4String(* 143 daughters_name[index] = new G4String(*right.daughters_name[index]);
117 } 144 }
118 } 145 }
119 pLambda = right.pLambda; << 146 pLambda = right.pLambda;
120 pXi0 = right.pXi0; 147 pXi0 = right.pXi0;
121 } 148 }
122 return *this; 149 return *this;
123 } 150 }
124 151
125 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4 << 152 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4double)
126 { 153 {
127 // this version neglects muon polarization << 154 // this version neglects muon polarization
128 // assumes the pure V-A couplin 155 // assumes the pure V-A coupling
129 // gives incorrect energy spect 156 // gives incorrect energy spectrum for Neutrinos
130 #ifdef G4VERBOSE 157 #ifdef G4VERBOSE
131 if (GetVerboseLevel() > 1) G4cout << "G4KL3D << 158 if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
132 #endif 159 #endif
133 160
134 // fill parent particle and its mass 161 // fill parent particle and its mass
135 CheckAndFillParent(); 162 CheckAndFillParent();
136 G4double massK = G4MT_parent->GetPDGMass(); 163 G4double massK = G4MT_parent->GetPDGMass();
137 164
138 // fill daughter particles and their mass 165 // fill daughter particles and their mass
139 CheckAndFillDaughters(); 166 CheckAndFillDaughters();
140 G4double daughterM[3]; 167 G4double daughterM[3];
141 daughterM[idPi] = G4MT_daughters[idPi]->GetP 168 daughterM[idPi] = G4MT_daughters[idPi]->GetPDGMass();
142 daughterM[idLepton] = G4MT_daughters[idLepto 169 daughterM[idLepton] = G4MT_daughters[idLepton]->GetPDGMass();
143 daughterM[idNutrino] = G4MT_daughters[idNutr 170 daughterM[idNutrino] = G4MT_daughters[idNutrino]->GetPDGMass();
144 171
145 // determine momentum/energy of daughters ac << 172 // determine momentum/energy of daughters according to DalitzDensity
146 G4double daughterP[3], daughterE[3]; 173 G4double daughterP[3], daughterE[3];
147 G4double w; 174 G4double w;
148 G4double r; 175 G4double r;
149 const size_t MAX_LOOP = 10000; 176 const size_t MAX_LOOP = 10000;
150 for (std::size_t loop_counter = 0; loop_coun << 177 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
>> 178 {
151 r = G4UniformRand(); 179 r = G4UniformRand();
152 PhaseSpace(massK, &daughterM[0], &daughter 180 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
153 w = DalitzDensity(massK, daughterE[idPi], << 181 w = DalitzDensity(massK, daughterE[idPi], daughterE[idLepton],
154 daughterM[idPi], daughte << 182 daughterE[idNutrino], daughterM[idPi],
155 if (r <= w) break; << 183 daughterM[idLepton], daughterM[idNutrino]);
>> 184 if ( r <= w) break;
156 } 185 }
157 186
158 // output message 187 // output message
159 #ifdef G4VERBOSE 188 #ifdef G4VERBOSE
160 if (GetVerboseLevel() > 1) { << 189 if (GetVerboseLevel()>1)
161 G4cout << *daughters_name[0] << ":" << dau << 190 {
162 G4cout << *daughters_name[1] << ":" << dau << 191 G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]"
163 G4cout << *daughters_name[2] << ":" << dau << 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;
164 } 197 }
165 #endif 198 #endif
166 199
167 // create parent G4DynamicParticle at rest 200 // create parent G4DynamicParticle at rest
168 auto direction = new G4ThreeVector(1.0, 0.0, << 201 G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0);
169 auto parentparticle = new G4DynamicParticle( << 202 G4DynamicParticle* parentparticle
>> 203 = new G4DynamicParticle( G4MT_parent, *direction, 0.0);
170 delete direction; 204 delete direction;
171 205
172 // create G4Decayproducts 206 // create G4Decayproducts
173 auto products = new G4DecayProducts(*parentp << 207 G4DecayProducts *products
>> 208 = new G4DecayProducts(*parentparticle);
174 delete parentparticle; 209 delete parentparticle;
175 210
176 // create daughter G4DynamicParticle << 211 // create daughter G4DynamicParticle
177 G4double costheta, sintheta, phi, sinphi, co << 212 G4double costheta, sintheta, phi, sinphi, cosphi;
178 G4double costhetan, sinthetan, phin, sinphin 213 G4double costhetan, sinthetan, phin, sinphin, cosphin;
179 << 214
180 // pion 215 // pion
181 costheta = 2. * G4UniformRand() - 1.0; << 216 costheta = 2.*G4UniformRand()-1.0;
182 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 217 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
183 phi = twopi * G4UniformRand() * rad; << 218 phi = twopi*G4UniformRand()*rad;
184 sinphi = std::sin(phi); 219 sinphi = std::sin(phi);
185 cosphi = std::cos(phi); 220 cosphi = std::cos(phi);
186 direction = new G4ThreeVector(sintheta * cos << 221 direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
187 G4ThreeVector momentum0 = (*direction) * dau << 222 G4ThreeVector momentum0 = (*direction)*daughterP[0];
188 auto daughterparticle = new G4DynamicParticl << 223 G4DynamicParticle * daughterparticle
>> 224 = new G4DynamicParticle(G4MT_daughters[0], momentum0);
189 products->PushProducts(daughterparticle); 225 products->PushProducts(daughterparticle);
190 226
191 // neutrino 227 // neutrino
192 costhetan = << 228 costhetan = (daughterP[1]*daughterP[1]
193 (daughterP[1] * daughterP[1] - daughterP[2 << 229 - daughterP[2]*daughterP[2]
194 / (2.0 * daughterP[2] * daughterP[0]); << 230 - daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
195 sinthetan = std::sqrt((1.0 - costhetan) * (1 << 231 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
196 phin = twopi * G4UniformRand() * rad; << 232 phin = twopi*G4UniformRand()*rad;
197 sinphin = std::sin(phin); 233 sinphin = std::sin(phin);
198 cosphin = std::cos(phin); 234 cosphin = std::cos(phin);
199 direction->setX(sinthetan * cosphin * costhe << 235 direction->setX( sinthetan*cosphin*costheta*cosphi
200 + costhetan * sintheta * cos << 236 - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
201 direction->setY(sinthetan * cosphin * costhe << 237 direction->setY( sinthetan*cosphin*costheta*sinphi
202 + costhetan * sintheta * sin << 238 + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
203 direction->setZ(-sinthetan * cosphin * sinth << 239 direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
204 240
205 G4ThreeVector momentum2 = (*direction) * dau << 241 G4ThreeVector momentum2 = (*direction)*daughterP[2];
206 daughterparticle = new G4DynamicParticle(G4M << 242 daughterparticle = new G4DynamicParticle( G4MT_daughters[2], momentum2);
207 products->PushProducts(daughterparticle); 243 products->PushProducts(daughterparticle);
208 244
209 // lepton 245 // lepton
210 G4ThreeVector momentum1 = (momentum0 + momen 246 G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
211 daughterparticle = new G4DynamicParticle(G4M << 247 daughterparticle = new G4DynamicParticle( G4MT_daughters[1], momentum1);
212 products->PushProducts(daughterparticle); 248 products->PushProducts(daughterparticle);
213 249
214 #ifdef G4VERBOSE 250 #ifdef G4VERBOSE
215 if (GetVerboseLevel() > 1) { << 251 if (GetVerboseLevel()>1)
216 G4cout << "G4KL3DecayChannel::DecayIt "; << 252 {
217 G4cout << " create decay products in rest << 253 G4cout << "G4KL3DecayChannel::DecayIt ";
218 G4cout << " decay products address=" << p << 254 G4cout << " create decay products in rest frame " <<G4endl;
219 products->DumpInfo(); << 255 G4cout << " decay products address=" << products << G4endl;
>> 256 products->DumpInfo();
220 } 257 }
221 #endif 258 #endif
222 delete direction; 259 delete direction;
223 return products; 260 return products;
224 } 261 }
225 262
226 void G4KL3DecayChannel::PhaseSpace(G4double pa << 263 void G4KL3DecayChannel::PhaseSpace(G4double parentM,
>> 264 const G4double* M,
>> 265 G4double* E,
>> 266 G4double* P )
227 { 267 {
228 // Algorithm in this code was originally wri 268 // Algorithm in this code was originally written in GDECA3 in GEANT3
229 << 269
230 // sum of daughters'mass 270 // sum of daughters'mass
231 G4double sumofdaughtermass = 0.0; 271 G4double sumofdaughtermass = 0.0;
232 G4int index; 272 G4int index;
233 const G4int N_DAUGHTER = 3; << 273 const G4int N_DAUGHTER=3;
234 << 274
235 for (index = 0; index < N_DAUGHTER; ++index) << 275 for (index=0; index<N_DAUGHTER; ++index)
>> 276 {
236 sumofdaughtermass += M[index]; 277 sumofdaughtermass += M[index];
237 } 278 }
238 279
239 // calculate daughter momentum. Generate two << 280 // calculate daughter momentum. Generate two
240 G4double rd1, rd2, rd; 281 G4double rd1, rd2, rd;
241 G4double momentummax = 0.0, momentumsum = 0. << 282 G4double momentummax=0.0, momentumsum = 0.0;
242 G4double energy; 283 G4double energy;
243 const size_t MAX_LOOP = 10000; << 284 const size_t MAX_LOOP=10000;
244 for (std::size_t loop_counter = 0; loop_coun << 285 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
>> 286 {
245 rd1 = G4UniformRand(); 287 rd1 = G4UniformRand();
246 rd2 = G4UniformRand(); 288 rd2 = G4UniformRand();
247 if (rd2 > rd1) { << 289 if (rd2 > rd1)
248 rd = rd1; << 290 {
>> 291 rd = rd1;
249 rd1 = rd2; 292 rd1 = rd2;
250 rd2 = rd; 293 rd2 = rd;
251 } << 294 }
252 momentummax = 0.0; 295 momentummax = 0.0;
253 momentumsum = 0.0; 296 momentumsum = 0.0;
254 // daughter 0 297 // daughter 0
255 energy = rd2 * (parentM - sumofdaughtermas << 298 energy = rd2*(parentM - sumofdaughtermass);
256 P[0] = std::sqrt(energy * energy + 2.0 * e << 299 P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
257 E[0] = energy; 300 E[0] = energy;
258 if (P[0] > momentummax) momentummax = P[0] << 301 if ( P[0] >momentummax )momentummax = P[0];
259 momentumsum += P[0]; << 302 momentumsum += P[0];
260 // daughter 1 303 // daughter 1
261 energy = (1. - rd1) * (parentM - sumofdaug << 304 energy = (1.-rd1)*(parentM - sumofdaughtermass);
262 P[1] = std::sqrt(energy * energy + 2.0 * e << 305 P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
263 E[1] = energy; 306 E[1] = energy;
264 if (P[1] > momentummax) momentummax = P[1] << 307 if ( P[1] >momentummax )momentummax = P[1];
265 momentumsum += P[1]; << 308 momentumsum += P[1];
266 // daughter 2 309 // daughter 2
267 energy = (rd1 - rd2) * (parentM - sumofdau << 310 energy = (rd1-rd2)*(parentM - sumofdaughtermass);
268 P[2] = std::sqrt(energy * energy + 2.0 * e << 311 P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
269 E[2] = energy; 312 E[2] = energy;
270 if (P[2] > momentummax) momentummax = P[2] << 313 if ( P[2] >momentummax )momentummax = P[2];
271 momentumsum += P[2]; << 314 momentumsum += P[2];
272 if (momentummax <= momentumsum - momentumm << 315 if (momentummax <= momentumsum - momentummax ) break;
273 } 316 }
274 #ifdef G4VERBOSE 317 #ifdef G4VERBOSE
275 if (GetVerboseLevel() > 2) { << 318 if (GetVerboseLevel()>2)
>> 319 {
276 G4cout << "G4KL3DecayChannel::PhaseSpace 320 G4cout << "G4KL3DecayChannel::PhaseSpace ";
277 G4cout << "Kon mass:" << parentM / GeV << << 321 G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl;
278 for (index = 0; index < 3; ++index) { << 322 for (index=0; index<3; ++index)
279 G4cout << index << " : " << M[index] / G << 323 {
280 G4cout << " : " << E[index] / GeV << "Ge << 324 G4cout << index << " : " << M[index]/GeV << "GeV/c/c ";
281 G4cout << " : " << P[index] / GeV << "Ge << 325 G4cout << " : " << E[index]/GeV << "GeV ";
>> 326 G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl;
282 } 327 }
283 } 328 }
284 #endif 329 #endif
285 } 330 }
286 331
287 G4double G4KL3DecayChannel::DalitzDensity(G4do << 332 G4double
288 G4do << 333 G4KL3DecayChannel::DalitzDensity(G4double massK, G4double Epi, G4double El,
>> 334 G4double Enu, G4double massPi, G4double massL,
>> 335 G4double massNu)
289 { 336 {
290 // KL3 decay - Dalitz Plot Density, see Chou 337 // KL3 decay - Dalitz Plot Density, see Chounet et al Phys. Rep. 4, 201
291 // Arguments 338 // Arguments
292 // Epi: kinetic enregy of pion 339 // Epi: kinetic enregy of pion
293 // El: kinetic enregy of lepton (e or mu 340 // El: kinetic enregy of lepton (e or mu)
294 // Enu: kinetic energy of nutrino 341 // Enu: kinetic energy of nutrino
295 // Constants 342 // Constants
296 // pLambda : linear energy dependence of 343 // pLambda : linear energy dependence of f+
297 // pXi0 : = f+(0)/f- 344 // pXi0 : = f+(0)/f-
298 // pNorm : normalization factor 345 // pNorm : normalization factor
299 // Variables 346 // Variables
300 // Epi: total energy of pion 347 // Epi: total energy of pion
301 // El: total energy of lepton (e or mu) 348 // El: total energy of lepton (e or mu)
302 // Enu: total energy of nutrino 349 // Enu: total energy of nutrino
303 350
304 // calculate total energy 351 // calculate total energy
305 Epi = Epi + massPi; 352 Epi = Epi + massPi;
306 El = El + massL; << 353 El = El + massL;
307 Enu = Enu + massNu; 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;
308 359
309 G4double Epi_max = (massK * massK + massPi * << 360 G4double F = 1.0 + pLambda*q2/massPi/massPi;
310 G4double E = Epi_max - Epi; <<
311 G4double q2 = massK * massK + massPi * massP <<
312 <<
313 G4double F = 1.0 + pLambda * q2 / massPi / m <<
314 G4double Fmax = 1.0; 361 G4double Fmax = 1.0;
315 if (pLambda > 0.0) Fmax = (1.0 + pLambda * ( << 362 if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
316 363
317 G4double Xi = pXi0 * (1.0 + pLambda * q2 / m << 364 G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
318 365
319 G4double coeffA = massK * (2.0 * El * Enu - << 366 G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
320 G4double coeffB = massL * massL * (Enu - E / << 367 G4double coeffB = massL*massL*(Enu-E/2.0);
321 G4double coeffC = massL * massL * E / 4.0; << 368 G4double coeffC = massL*massL*E/4.0;
322 369
323 G4double RhoMax = (Fmax * Fmax) * (massK * m << 370 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
324 <<
325 G4double Rho = (F * F) * (coeffA + coeffB * <<
326 371
>> 372 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
>> 373
327 #ifdef G4VERBOSE 374 #ifdef G4VERBOSE
328 if (GetVerboseLevel() > 2) { << 375 if (GetVerboseLevel()>2)
329 G4cout << "G4KL3DecayChannel::DalitzDensit << 376 {
330 G4cout << " Pi[" << massPi / GeV << "GeV/c << 377 G4cout << "G4KL3DecayChannel::DalitzDensity " <<G4endl;
331 G4cout << " L[" << massL / GeV << "GeV/c/c << 378 G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl;
332 G4cout << " Nu[" << massNu / GeV << "GeV/c << 379 G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl;
333 G4cout << " F :" << F << " Fmax :" << Fmax << 380 G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl;
334 G4cout << " A :" << coeffA << " B :" << c << 381 G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl;
335 G4cout << " Rho :" << Rho << " RhoMax :" << 382 G4cout << " A :" << coeffA << " B :" << coeffB << " C :"<< coeffC
>> 383 << G4endl;
>> 384 G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl;
336 } 385 }
337 #endif 386 #endif
338 return (Rho / RhoMax); << 387 return (Rho/RhoMax);
339 } 388 }
340 389