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 10.3)


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