Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4KL3DecayChannel class implementation << 27 // 23 // 28 // Author: H.Kurashige, 30 May 1997 << 24 // $Id: G4KL3DecayChannel.cc,v 1.5 2001/07/11 10:02:00 gunter Exp $ 29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-04-00 $ 30 << 26 // 31 #include "G4KL3DecayChannel.hh" << 27 // >> 28 // ------------------------------------------------------------ >> 29 // GEANT 4 class header file >> 30 // >> 31 // History: first implementation, based on object model of >> 32 // 30 May 1997 H.Kurashige >> 33 // ------------------------------------------------------------ 32 34 33 #include "G4DecayProducts.hh" << 34 #include "G4LorentzRotation.hh" << 35 #include "G4LorentzVector.hh" << 36 #include "G4ParticleDefinition.hh" 35 #include "G4ParticleDefinition.hh" 37 #include "G4PhysicalConstants.hh" << 36 #include "G4DecayProducts.hh" 38 #include "G4SystemOfUnits.hh" << 39 #include "G4VDecayChannel.hh" 37 #include "G4VDecayChannel.hh" >> 38 #include "G4KL3DecayChannel.hh" 40 #include "Randomize.hh" 39 #include "Randomize.hh" >> 40 #include "G4LorentzVector.hh" >> 41 #include "G4LorentzRotation.hh" 41 42 42 G4KL3DecayChannel::G4KL3DecayChannel(const G4S << 43 const G4S << 44 const G4S << 45 : G4VDecayChannel("KL3 Decay", theParentName << 46 theNutrinoName) << 47 { << 48 static const G4String K_plus("kaon+"); << 49 static const G4String K_minus("kaon-"); << 50 static const G4String K_L("kaon0L"); << 51 static const G4String Mu_plus("mu+"); << 52 static const G4String Mu_minus("mu-"); << 53 static const G4String E_plus("e+"); << 54 static const G4String E_minus("e-"); << 55 43 >> 44 G4KL3DecayChannel::G4KL3DecayChannel( >> 45 const G4String& theParentName, >> 46 G4double theBR, >> 47 const G4String& thePionName, >> 48 const G4String& theLeptonName, >> 49 const G4String& theNutrinoName) >> 50 :G4VDecayChannel("KL3 Decay",theParentName, >> 51 theBR, 3, >> 52 thePionName,theLeptonName,theNutrinoName) >> 53 { >> 54 //#ifdef G4VERBOSE >> 55 //if (GetVerboseLevel()>1) { >> 56 // G4cout << "G4KL3DecayChannel:: constructor "; >> 57 // G4cout << "addr[" << this << "]" << G4endl; >> 58 //} >> 59 //#endif 56 // check modes 60 // check modes 57 if (((theParentName == K_plus) && (theLepton << 61 if ( ((theParentName == "kaon+")&&(theLeptonName == "e+")) || 58 || ((theParentName == K_minus) && (theLe << 62 ((theParentName == "kaon-")&&(theLeptonName == "e-")) ) { 59 { << 60 // K+- (Ke3) 63 // K+- (Ke3) 61 pLambda = 0.0286; 64 pLambda = 0.0286; 62 pXi0 = -0.35; << 65 pXi0 = -0.35; 63 } << 66 } else if ( ((theParentName == "kaon+")&&(theLeptonName == "mu+")) || 64 else if (((theParentName == K_plus) && (theL << 67 ((theParentName == "kaon-")&&(theLeptonName == "mu-")) ) { 65 || ((theParentName == K_minus) && ( << 66 { << 67 // K+- (Kmu3) 68 // K+- (Kmu3) 68 pLambda = 0.033; 69 pLambda = 0.033; 69 pXi0 = -0.35; << 70 pXi0 = -0.35; 70 } << 71 } else if ( (theParentName == "kaon0L") && 71 else if ((theParentName == K_L) && ((theLept << 72 ((theLeptonName == "e+") ||(theLeptonName == "e-")) ){ 72 // K0L (Ke3) 73 // K0L (Ke3) 73 pLambda = 0.0300; 74 pLambda = 0.0300; 74 pXi0 = -0.11; << 75 pXi0 = -0.11; 75 } << 76 } else if ( (theParentName == "kaon0L") && 76 else if ((theParentName == K_L) && ((theLept << 77 ((theLeptonName == "mu+") ||(theLeptonName == "mu-")) ){ 77 // K0L (Kmu3) 78 // K0L (Kmu3) 78 pLambda = 0.034; 79 pLambda = 0.034; 79 pXi0 = -0.11; << 80 pXi0 = -0.11; 80 } << 81 } else { 81 else { << 82 //#ifdef G4VERBOSE 82 #ifdef G4VERBOSE << 83 //if (GetVerboseLevel()>0) { 83 if (GetVerboseLevel() > 2) { << 84 // G4cout << "G4KL3DecayChannel:: constructor :"; 84 G4cout << "G4KL3DecayChannel:: construct << 85 // G4cout << "illegal arguments " << G4endl;; 85 G4cout << "illegal arguments " << G4endl << 86 // DumpInfo(); 86 ; << 87 // } 87 DumpInfo(); << 88 //#endif 88 } << 89 #endif << 90 // set values for K0L (Ke3) temporarily 89 // set values for K0L (Ke3) temporarily 91 pLambda = 0.0300; 90 pLambda = 0.0300; 92 pXi0 = -0.11; << 91 pXi0 = -0.11; 93 } 92 } 94 } 93 } 95 94 96 G4KL3DecayChannel& G4KL3DecayChannel::operator << 95 G4KL3DecayChannel::~G4KL3DecayChannel() 97 { 96 { 98 if (this != &right) { << 99 kinematics_name = right.kinematics_name; << 100 verboseLevel = right.verboseLevel; << 101 rbranch = right.rbranch; << 102 << 103 // copy parent name << 104 parent_name = new G4String(*right.parent_n << 105 << 106 // clear daughters_name array << 107 ClearDaughtersName(); << 108 << 109 // recreate array << 110 numberOfDaughters = right.numberOfDaughter << 111 if (numberOfDaughters > 0) { << 112 if (daughters_name != nullptr) ClearDaug << 113 daughters_name = new G4String*[numberOfD << 114 // copy daughters name << 115 for (G4int index = 0; index < numberOfDa << 116 daughters_name[index] = new G4String(* << 117 } << 118 } << 119 pLambda = right.pLambda; << 120 pXi0 = right.pXi0; << 121 } << 122 return *this; << 123 } 97 } 124 98 125 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4 << 99 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4double) 126 { 100 { 127 // this version neglects muon polarization << 101 // this version neglects muon polarization 128 // assumes the pure V-A couplin 102 // assumes the pure V-A coupling 129 // gives incorrect energy spect << 103 // gives incorrect energy spectrum for Nutrinos 130 #ifdef G4VERBOSE 104 #ifdef G4VERBOSE 131 if (GetVerboseLevel() > 1) G4cout << "G4KL3D << 105 if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl; 132 #endif 106 #endif 133 << 134 // fill parent particle and its mass 107 // fill parent particle and its mass 135 CheckAndFillParent(); << 108 if (parent == 0) { 136 G4double massK = G4MT_parent->GetPDGMass(); << 109 FillParent(); >> 110 } >> 111 massK = parent->GetPDGMass(); 137 112 138 // fill daughter particles and their mass 113 // fill daughter particles and their mass 139 CheckAndFillDaughters(); << 114 if (daughters == 0) { 140 G4double daughterM[3]; << 115 FillDaughters(); 141 daughterM[idPi] = G4MT_daughters[idPi]->GetP << 116 } 142 daughterM[idLepton] = G4MT_daughters[idLepto << 117 daughterM[idPi] = daughters[idPi]->GetPDGMass(); 143 daughterM[idNutrino] = G4MT_daughters[idNutr << 118 daughterM[idLepton] = daughters[idLepton]->GetPDGMass(); >> 119 daughterM[idNutrino] = daughters[idNutrino]->GetPDGMass(); 144 120 145 // determine momentum/energy of daughters ac << 121 // determine momentum/energy of daughters >> 122 // according to DalitzDensity 146 G4double daughterP[3], daughterE[3]; 123 G4double daughterP[3], daughterE[3]; 147 G4double w; 124 G4double w; 148 G4double r; 125 G4double r; 149 const size_t MAX_LOOP = 10000; << 126 do { 150 for (std::size_t loop_counter = 0; loop_coun << 151 r = G4UniformRand(); 127 r = G4UniformRand(); 152 PhaseSpace(massK, &daughterM[0], &daughter 128 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]); 153 w = DalitzDensity(massK, daughterE[idPi], << 129 w = DalitzDensity(daughterE[idPi],daughterE[idLepton],daughterE[idNutrino]); 154 daughterM[idPi], daughte << 130 } while ( r > w); 155 if (r <= w) break; << 156 } << 157 131 158 // output message 132 // output message 159 #ifdef G4VERBOSE 133 #ifdef G4VERBOSE 160 if (GetVerboseLevel() > 1) { << 134 if (GetVerboseLevel()>1) { 161 G4cout << *daughters_name[0] << ":" << dau << 135 G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]" <<G4endl; 162 G4cout << *daughters_name[1] << ":" << dau << 136 G4cout << *daughters_name[1] << ":" << daughterP[1]/GeV << "[GeV/c]" <<G4endl; 163 G4cout << *daughters_name[2] << ":" << dau << 137 G4cout << *daughters_name[2] << ":" << daughterP[2]/GeV << "[GeV/c]" <<G4endl; 164 } 138 } 165 #endif 139 #endif 166 << 140 //create parent G4DynamicParticle at rest 167 // create parent G4DynamicParticle at rest << 141 G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0); 168 auto direction = new G4ThreeVector(1.0, 0.0, << 142 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, *direction, 0.0); 169 auto parentparticle = new G4DynamicParticle( << 170 delete direction; 143 delete direction; 171 144 172 // create G4Decayproducts << 145 //create G4Decayproducts 173 auto products = new G4DecayProducts(*parentp << 146 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 174 delete parentparticle; 147 delete parentparticle; 175 148 176 // create daughter G4DynamicParticle << 149 //create daughter G4DynamicParticle 177 G4double costheta, sintheta, phi, sinphi, co << 150 G4double costheta, sintheta, phi, sinphi, cosphi; 178 G4double costhetan, sinthetan, phin, sinphin 151 G4double costhetan, sinthetan, phin, sinphin, cosphin; 179 << 152 180 // pion 153 // pion 181 costheta = 2. * G4UniformRand() - 1.0; << 154 costheta = 2.*G4UniformRand()-1.0; 182 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 155 sintheta = sqrt((1.0-costheta)*(1.0+costheta)); 183 phi = twopi * G4UniformRand() * rad; << 156 phi = 2.0*M_PI*G4UniformRand()*rad; 184 sinphi = std::sin(phi); << 157 sinphi = sin(phi); 185 cosphi = std::cos(phi); << 158 cosphi = cos(phi); 186 direction = new G4ThreeVector(sintheta * cos << 159 direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta); 187 G4ThreeVector momentum0 = (*direction) * dau << 160 G4ThreeVector momentum0 = (*direction)*daughterP[0]; 188 auto daughterparticle = new G4DynamicParticl << 161 G4DynamicParticle * daughterparticle >> 162 = new G4DynamicParticle( daughters[0], momentum0); 189 products->PushProducts(daughterparticle); 163 products->PushProducts(daughterparticle); 190 164 191 // neutrino 165 // neutrino 192 costhetan = << 166 costhetan = (daughterP[1]*daughterP[1]-daughterP[2]*daughterP[2]-daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]); 193 (daughterP[1] * daughterP[1] - daughterP[2 << 167 sinthetan = sqrt((1.0-costhetan)*(1.0+costhetan)); 194 / (2.0 * daughterP[2] * daughterP[0]); << 168 phin = 2.0*M_PI*G4UniformRand()*rad; 195 sinthetan = std::sqrt((1.0 - costhetan) * (1 << 169 sinphin = sin(phin); 196 phin = twopi * G4UniformRand() * rad; << 170 cosphin = cos(phin); 197 sinphin = std::sin(phin); << 171 direction->setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 198 cosphin = std::cos(phin); << 172 direction->setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 199 direction->setX(sinthetan * cosphin * costhe << 173 direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta); 200 + costhetan * sintheta * cos << 201 direction->setY(sinthetan * cosphin * costhe << 202 + costhetan * sintheta * sin << 203 direction->setZ(-sinthetan * cosphin * sinth << 204 174 205 G4ThreeVector momentum2 = (*direction) * dau << 175 G4ThreeVector momentum2 = (*direction)*daughterP[2]; 206 daughterparticle = new G4DynamicParticle(G4M << 176 daughterparticle = new G4DynamicParticle( daughters[2], momentum2); 207 products->PushProducts(daughterparticle); 177 products->PushProducts(daughterparticle); 208 178 209 // lepton << 179 //lepton 210 G4ThreeVector momentum1 = (momentum0 + momen 180 G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0); 211 daughterparticle = new G4DynamicParticle(G4M << 181 daughterparticle = >> 182 new G4DynamicParticle( daughters[1], momentum1); 212 products->PushProducts(daughterparticle); 183 products->PushProducts(daughterparticle); 213 184 214 #ifdef G4VERBOSE 185 #ifdef G4VERBOSE 215 if (GetVerboseLevel() > 1) { << 186 if (GetVerboseLevel()>1) { 216 G4cout << "G4KL3DecayChannel::DecayIt "; << 187 G4cout << "G4KL3DecayChannel::DecayIt "; 217 G4cout << " create decay products in rest << 188 G4cout << " create decay products in rest frame " <<G4endl; 218 G4cout << " decay products address=" << p << 189 G4cout << " decay products address=" << products << G4endl; 219 products->DumpInfo(); << 190 products->DumpInfo(); 220 } 191 } 221 #endif 192 #endif 222 delete direction; 193 delete direction; 223 return products; 194 return products; 224 } 195 } 225 196 226 void G4KL3DecayChannel::PhaseSpace(G4double pa << 197 void G4KL3DecayChannel::PhaseSpace(G4double parentM, >> 198 const G4double* M, >> 199 G4double* E, >> 200 G4double* P ) >> 201 // algorism of this code is originally written in GDECA3 of GEANT3 227 { 202 { 228 // Algorithm in this code was originally wri << 203 229 << 204 //sum of daughters'mass 230 // sum of daughters'mass << 231 G4double sumofdaughtermass = 0.0; 205 G4double sumofdaughtermass = 0.0; 232 G4int index; 206 G4int index; 233 const G4int N_DAUGHTER = 3; << 207 for (index=0; index<3; index++){ 234 << 235 for (index = 0; index < N_DAUGHTER; ++index) << 236 sumofdaughtermass += M[index]; 208 sumofdaughtermass += M[index]; 237 } 209 } 238 210 239 // calculate daughter momentum. Generate two << 211 //calculate daughter momentum >> 212 // Generate two 240 G4double rd1, rd2, rd; 213 G4double rd1, rd2, rd; 241 G4double momentummax = 0.0, momentumsum = 0. << 214 G4double momentummax=0.0, momentumsum = 0.0; 242 G4double energy; 215 G4double energy; 243 const size_t MAX_LOOP = 10000; << 216 244 for (std::size_t loop_counter = 0; loop_coun << 217 do { 245 rd1 = G4UniformRand(); 218 rd1 = G4UniformRand(); 246 rd2 = G4UniformRand(); 219 rd2 = G4UniformRand(); 247 if (rd2 > rd1) { 220 if (rd2 > rd1) { 248 rd = rd1; << 221 rd = rd1; 249 rd1 = rd2; 222 rd1 = rd2; 250 rd2 = rd; 223 rd2 = rd; 251 } << 224 } 252 momentummax = 0.0; 225 momentummax = 0.0; 253 momentumsum = 0.0; 226 momentumsum = 0.0; 254 // daughter 0 227 // daughter 0 255 energy = rd2 * (parentM - sumofdaughtermas << 228 energy = rd2*(parentM - sumofdaughtermass); 256 P[0] = std::sqrt(energy * energy + 2.0 * e << 229 P[0] = sqrt(energy*energy + 2.0*energy*M[0]); 257 E[0] = energy; 230 E[0] = energy; 258 if (P[0] > momentummax) momentummax = P[0] << 231 if ( P[0] >momentummax )momentummax = P[0]; 259 momentumsum += P[0]; << 232 momentumsum += P[0]; 260 // daughter 1 233 // daughter 1 261 energy = (1. - rd1) * (parentM - sumofdaug << 234 energy = (1.-rd1)*(parentM - sumofdaughtermass); 262 P[1] = std::sqrt(energy * energy + 2.0 * e << 235 P[1] = sqrt(energy*energy + 2.0*energy*M[1]); 263 E[1] = energy; 236 E[1] = energy; 264 if (P[1] > momentummax) momentummax = P[1] << 237 if ( P[1] >momentummax )momentummax = P[1]; 265 momentumsum += P[1]; << 238 momentumsum += P[1]; 266 // daughter 2 239 // daughter 2 267 energy = (rd1 - rd2) * (parentM - sumofdau << 240 energy = (rd1-rd2)*(parentM - sumofdaughtermass); 268 P[2] = std::sqrt(energy * energy + 2.0 * e << 241 P[2] = sqrt(energy*energy + 2.0*energy*M[2]); 269 E[2] = energy; 242 E[2] = energy; 270 if (P[2] > momentummax) momentummax = P[2] << 243 if ( P[2] >momentummax )momentummax = P[2]; 271 momentumsum += P[2]; << 244 momentumsum += P[2]; 272 if (momentummax <= momentumsum - momentumm << 245 } while (momentummax > momentumsum - momentummax ); 273 } << 246 274 #ifdef G4VERBOSE 247 #ifdef G4VERBOSE 275 if (GetVerboseLevel() > 2) { << 248 if (GetVerboseLevel()>2) { 276 G4cout << "G4KL3DecayChannel::PhaseSpace << 249 G4cout << "G4KL3DecayChannel::PhaseSpace "; 277 G4cout << "Kon mass:" << parentM / GeV << << 250 G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl; 278 for (index = 0; index < 3; ++index) { << 251 for (index=0; index<3; index++){ 279 G4cout << index << " : " << M[index] / G << 252 G4cout << index << " : " << M[index]/GeV << "GeV/c/c "; 280 G4cout << " : " << E[index] / GeV << "Ge << 253 G4cout << " : " << E[index]/GeV << "GeV "; 281 G4cout << " : " << P[index] / GeV << "Ge << 254 G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl; 282 } << 255 } 283 } 256 } 284 #endif 257 #endif 285 } 258 } 286 259 287 G4double G4KL3DecayChannel::DalitzDensity(G4do << 260 288 G4do << 261 G4double G4KL3DecayChannel::DalitzDensity(G4double Epi, G4double El, G4double Enu) 289 { 262 { 290 // KL3 decay - Dalitz Plot Density, see Chou << 263 // KL3 decay Dalitz Plot Density 291 // Arguments << 264 // see Chounet et al Phys. Rep. 4, 201 >> 265 // arguments 292 // Epi: kinetic enregy of pion 266 // Epi: kinetic enregy of pion 293 // El: kinetic enregy of lepton (e or mu 267 // El: kinetic enregy of lepton (e or mu) 294 // Enu: kinetic energy of nutrino 268 // Enu: kinetic energy of nutrino 295 // Constants << 269 // constants 296 // pLambda : linear energy dependence of 270 // pLambda : linear energy dependence of f+ 297 // pXi0 : = f+(0)/f- 271 // pXi0 : = f+(0)/f- 298 // pNorm : normalization factor 272 // pNorm : normalization factor 299 // Variables << 273 // variables 300 // Epi: total energy of pion 274 // Epi: total energy of pion 301 // El: total energy of lepton (e or mu) 275 // El: total energy of lepton (e or mu) 302 // Enu: total energy of nutrino 276 // Enu: total energy of nutrino 303 277 304 // calculate total energy << 278 // mass of daughters >> 279 G4double massPi = daughterM[idPi]; >> 280 G4double massL = daughterM[idLepton]; >> 281 G4double massNu = daughterM[idNutrino]; >> 282 >> 283 // calcurate total energy 305 Epi = Epi + massPi; 284 Epi = Epi + massPi; 306 El = El + massL; << 285 El = El + massL; 307 Enu = Enu + massNu; 286 Enu = Enu + massNu; >> 287 >> 288 G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK; >> 289 G4double E = Epi_max - Epi; >> 290 G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi; 308 291 309 G4double Epi_max = (massK * massK + massPi * << 292 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; 293 G4double Fmax = 1.0; 315 if (pLambda > 0.0) Fmax = (1.0 + pLambda * ( << 294 if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0)); 316 << 317 G4double Xi = pXi0 * (1.0 + pLambda * q2 / m << 318 295 319 G4double coeffA = massK * (2.0 * El * Enu - << 296 G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi); 320 G4double coeffB = massL * massL * (Enu - E / << 321 G4double coeffC = massL * massL * E / 4.0; << 322 297 323 G4double RhoMax = (Fmax * Fmax) * (massK * m << 298 G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu); >> 299 G4double coeffB = massL*massL*(Enu-E/2.0); >> 300 G4double coeffC = massL*massL*E/4.0; 324 301 325 G4double Rho = (F * F) * (coeffA + coeffB * << 302 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0); 326 303 >> 304 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi); >> 305 327 #ifdef G4VERBOSE 306 #ifdef G4VERBOSE 328 if (GetVerboseLevel() > 2) { << 307 if (GetVerboseLevel()>2) { 329 G4cout << "G4KL3DecayChannel::DalitzDensit << 308 G4cout << "G4KL3DecayChannel::DalitzDensity " <<G4endl; 330 G4cout << " Pi[" << massPi / GeV << "GeV/c << 309 G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl; 331 G4cout << " L[" << massL / GeV << "GeV/c/c << 310 G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl; 332 G4cout << " Nu[" << massNu / GeV << "GeV/c << 311 G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl; 333 G4cout << " F :" << F << " Fmax :" << Fmax << 312 G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl; 334 G4cout << " A :" << coeffA << " B :" << c << 313 G4cout << " A :" << coeffA << " B :" << coeffB << " C :"<< coeffC <<G4endl; 335 G4cout << " Rho :" << Rho << " RhoMax :" << 314 G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl; 336 } 315 } 337 #endif 316 #endif 338 return (Rho / RhoMax); << 317 return (Rho/RhoMax); 339 } 318 } >> 319 >> 320 340 321