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.6.p1)


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