Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/particles/management/src/G4KL3DecayChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /particles/management/src/G4KL3DecayChannel.cc (Version 11.3.0) and /particles/management/src/G4KL3DecayChannel.cc (Version 9.4.p4)


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