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