Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/particles/management/src/G4PhaseSpaceDecayChannel.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/G4PhaseSpaceDecayChannel.cc (Version 11.3.0) and /particles/management/src/G4PhaseSpaceDecayChannel.cc (Version 8.2.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 // G4PhaseSpaceDecayChannel class implementati << 
 27 //                                                 26 //
 28 // Author: H.Kurashige, 27 July 1996           <<  27 // $Id: G4PhaseSpaceDecayChannel.cc,v 1.12 2006/09/28 14:30:33 kurasige Exp $
 29 // ------------------------------------------- <<  28 // GEANT4 tag $Name: geant4-08-02-patch-01-ref $
 30                                                <<  29 //
 31 #include "G4PhaseSpaceDecayChannel.hh"         <<  30 // 
                                                   >>  31 // ------------------------------------------------------------
                                                   >>  32 //      GEANT 4 class header file
                                                   >>  33 //
                                                   >>  34 //      History: first implementation, based on object model of
                                                   >>  35 //      27 July 1996 H.Kurashige
                                                   >>  36 //      20 Oct 1996 H.Kurashige
                                                   >>  37 //      30 May 1997 H.Kurashige
                                                   >>  38 // ------------------------------------------------------------
 32                                                    39 
 33 #include "G4DecayProducts.hh"                  << 
 34 #include "G4LorentzRotation.hh"                << 
 35 #include "G4LorentzVector.hh"                  << 
 36 #include "G4ParticleDefinition.hh"                 40 #include "G4ParticleDefinition.hh"
 37 #include "G4PhysicalConstants.hh"              <<  41 #include "G4DecayProducts.hh"
 38 #include "G4SystemOfUnits.hh"                  << 
 39 #include "G4VDecayChannel.hh"                      42 #include "G4VDecayChannel.hh"
                                                   >>  43 #include "G4PhaseSpaceDecayChannel.hh"
 40 #include "Randomize.hh"                            44 #include "Randomize.hh"
                                                   >>  45 #include "G4LorentzVector.hh"
                                                   >>  46 #include "G4LorentzRotation.hh"
                                                   >>  47 
 41                                                    48 
 42 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayCha     49 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
 43   : G4VDecayChannel("Phase Space", Verbose)    <<  50  :G4VDecayChannel("Phase Space", Verbose)
 44 {}                                             <<  51 {
 45                                                    52 
 46 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayCha <<  53 }
 47                                                << 
 48                                                << 
 49                                                << 
 50                                                << 
 51                                                << 
 52                                                << 
 53   : G4VDecayChannel("Phase Space", theParentNa << 
 54                     theDaughterName2, theDaugh << 
 55 {}                                             << 
 56                                                    54 
 57 G4DecayProducts* G4PhaseSpaceDecayChannel::Dec <<  55 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
                                                   >>  56              const G4String& theParentName,
                                                   >>  57              G4double        theBR,
                                                   >>  58              G4int           theNumberOfDaughters,
                                                   >>  59              const G4String& theDaughterName1,
                                                   >>  60              const G4String& theDaughterName2,
                                                   >>  61              const G4String& theDaughterName3,
                                                   >>  62              const G4String& theDaughterName4    )
                                                   >>  63                          :G4VDecayChannel("Phase Space",
                                                   >>  64             theParentName,theBR,
                                                   >>  65             theNumberOfDaughters,
                                                   >>  66             theDaughterName1,
                                                   >>  67             theDaughterName2,
                                                   >>  68             theDaughterName3,
                                                   >>  69             theDaughterName4)
 58 {                                                  70 {
 59 #ifdef G4VERBOSE                               << 
 60   if (GetVerboseLevel() > 1) G4cout << "G4Phas << 
 61 #endif                                         << 
 62                                                    71 
 63   G4DecayProducts* products = nullptr;         <<  72 }
 64                                                << 
 65   CheckAndFillParent();                        << 
 66   CheckAndFillDaughters();                     << 
 67                                                    73 
 68   if (parentMass > 0.0)                        <<  74 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
 69     current_parent_mass.Put(parentMass);       <<  75 {
 70   else                                         <<  76 }
 71     current_parent_mass.Put(G4MT_parent_mass); << 
 72                                                    77 
 73   switch (numberOfDaughters) {                 <<  78 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
 74     case 0:                                    <<  79 {
 75 #ifdef G4VERBOSE                                   80 #ifdef G4VERBOSE
 76       if (GetVerboseLevel() > 0) {             <<  81   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
 77         G4cout << "G4PhaseSpaceDecayChannel::D << 
 78         G4cout << " daughters not defined " << << 
 79       }                                        << 
 80 #endif                                             82 #endif
 81       break;                                   <<  83   
 82     case 1:                                    <<  84   G4DecayProducts * products = 0;
 83       products = OneBodyDecayIt();             <<  85  
 84       break;                                   <<  86   if (parent == 0) FillParent();  
 85     case 2:                                    <<  87   if (daughters == 0) FillDaughters();
 86       products = TwoBodyDecayIt();             <<  88 
 87       break;                                   <<  89   if (parentMass >0.0) parent_mass = parentMass;  
 88     case 3:                                    <<  90   switch (numberOfDaughters){
 89       products = ThreeBodyDecayIt();           <<  91   case 0:
 90       break;                                   <<  92 #ifdef G4VERBOSE
 91     default:                                   <<  93     if (GetVerboseLevel()>0) {
 92       products = ManyBodyDecayIt();            <<  94       G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
 93       break;                                   <<  95       G4cout << " daughters not defined " <<G4endl;
 94   }                                            <<  96     }
 95 #ifdef G4VERBOSE                               <<  97 #endif
 96   if ((products == nullptr) && (GetVerboseLeve <<  98     break;
 97     G4cout << "G4PhaseSpaceDecayChannel::Decay <<  99   case 1:
 98     G4cout << *parent_name << " cannot decay " << 100     products =  OneBodyDecayIt();
                                                   >> 101     break;
                                                   >> 102   case 2:
                                                   >> 103     products =  TwoBodyDecayIt();
                                                   >> 104     break;
                                                   >> 105   case 3:
                                                   >> 106     products =  ThreeBodyDecayIt();
                                                   >> 107     break;
                                                   >> 108   default:
                                                   >> 109     products =  ManyBodyDecayIt();
                                                   >> 110     break;
                                                   >> 111   }
                                                   >> 112 #ifdef G4VERBOSE
                                                   >> 113   if ((products == 0) && (GetVerboseLevel()>0)) {
                                                   >> 114     G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
                                                   >> 115     G4cout << *parent_name << " can not decay " << G4endl;
 99     DumpInfo();                                   116     DumpInfo();
100   }                                               117   }
101 #endif                                            118 #endif
102   return products;                                119   return products;
103 }                                                 120 }
104                                                   121 
105 G4DecayProducts* G4PhaseSpaceDecayChannel::One << 122 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
106 {                                                 123 {
107 #ifdef G4VERBOSE                                  124 #ifdef G4VERBOSE
108   if (GetVerboseLevel() > 1) G4cout << "G4Phas << 125   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
109 #endif                                            126 #endif
110   // parent mass                               << 
111   G4double parentmass = current_parent_mass.Ge << 
112                                                   127 
113   // create parent G4DynamicParticle at rest   << 128   //create parent G4DynamicParticle at rest
114   G4ThreeVector dummy;                            129   G4ThreeVector dummy;
115   auto parentparticle = new G4DynamicParticle( << 130   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
116   // create G4Decayproducts                    << 131   //create G4Decayproducts
117   auto products = new G4DecayProducts(*parentp << 132   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
118   delete parentparticle;                          133   delete parentparticle;
119                                                   134 
120   // create daughter G4DynamicParticle at rest << 135   //create daughter G4DynamicParticle at rest
121   auto daughterparticle = new G4DynamicParticl << 136   G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
122   if (useGivenDaughterMass) daughterparticle-> << 
123   products->PushProducts(daughterparticle);       137   products->PushProducts(daughterparticle);
124                                                   138 
125 #ifdef G4VERBOSE                                  139 #ifdef G4VERBOSE
126   if (GetVerboseLevel() > 1) {                 << 140   if (GetVerboseLevel()>1) {
127     G4cout << "G4PhaseSpaceDecayChannel::OneBo << 141      G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
128     G4cout << " create decay products in rest  << 142      G4cout << "  create decay products in rest frame " <<G4endl;
129     products->DumpInfo();                      << 143      products->DumpInfo();
130   }                                               144   }
131 #endif                                            145 #endif
132   return products;                                146   return products;
133 }                                                 147 }
134                                                   148 
135 G4DecayProducts* G4PhaseSpaceDecayChannel::Two << 149 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
136 {                                                 150 {
137 #ifdef G4VERBOSE                                  151 #ifdef G4VERBOSE
138   if (GetVerboseLevel() > 1) G4cout << "G4Phas << 152   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
139 #endif                                            153 #endif
140   // parent mass                                  154   // parent mass
141   G4double parentmass = current_parent_mass.Ge << 155   G4double parentmass = parent_mass;
                                                   >> 156   
                                                   >> 157   //daughters'mass
                                                   >> 158   G4double daughtermass[2]; 
                                                   >> 159   G4double daughtermomentum;
                                                   >> 160   daughtermass[0] = daughters_mass[0];
                                                   >> 161   daughtermass[1] = daughters_mass[1];
142                                                   162 
143   // daughters'mass, width                     << 163   //create parent G4DynamicParticle at rest
144   G4double daughtermass[2], daughterwidth[2];  << 
145   daughtermass[0] = G4MT_daughters_mass[0];    << 
146   daughtermass[1] = G4MT_daughters_mass[1];    << 
147   daughterwidth[0] = G4MT_daughters_width[0];  << 
148   daughterwidth[1] = G4MT_daughters_width[1];  << 
149                                                << 
150   // create parent G4DynamicParticle at rest   << 
151   G4ThreeVector dummy;                            164   G4ThreeVector dummy;
152   auto parentparticle = new G4DynamicParticle( << 165   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
153   // create G4Decayproducts                    << 166   //create G4Decayproducts
154   auto products = new G4DecayProducts(*parentp << 167   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
155   delete parentparticle;                          168   delete parentparticle;
156                                                   169 
157   if (!useGivenDaughterMass) {                 << 170   //calculate daughter momentum
158     G4bool withWidth = (daughterwidth[0] > 1.0 << 171   daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
159                        || (daughterwidth[1] >  << 172   if (daughtermomentum <0.0) {
160     if (withWidth) {                           << 173 #ifdef G4VERBOSE
161       G4double sumofdaughterwidthsq =          << 174     if (GetVerboseLevel()>0) {
162         daughterwidth[0] * daughterwidth[0] +  << 175       G4cerr << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt " 
163       // check parent mass and daughter mass   << 176              << "sum of daughter mass is larger than parent mass" << G4endl;
164       G4double maxDev =                        << 177       G4cerr << "parent :" << parent->GetParticleName() << "  " << parent_mass/GeV << G4endl;
165         (parentmass - daughtermass[0] - daught << 178       G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << "  " << daughtermass[0]/GeV << G4endl;
166       if (maxDev <= -1.0 * rangeMass) {        << 179       G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << "  " << daughtermass[1]/GeV << G4endl;
167 #ifdef G4VERBOSE                               << 
168         if (GetVerboseLevel() > 0) {           << 
169           G4cout << "G4PhaseSpaceDecayChannel: << 
170                  << "Sum of daughter mass is l << 
171           G4cout << "Parent :" << G4MT_parent- << 
172                  << current_parent_mass.Get()  << 
173           G4cout << "Daughter 1 :" << G4MT_dau << 
174                  << daughtermass[0] / GeV << G << 
175           G4cout << "Daughter 2:" << G4MT_daug << 
176                  << daughtermass[1] / GeV << G << 
177         }                                      << 
178 #endif                                         << 
179         G4Exception("G4PhaseSpaceDecayChannel: << 
180                     "Cannot create decay produ << 
181                      larger than parent mass!" << 
182         return products;                       << 
183       }                                        << 
184       G4double dm1 = daughtermass[0];          << 
185       if (daughterwidth[0] > 0.) dm1 = Dynamic << 
186       G4double dm2 = daughtermass[1];          << 
187       if (daughterwidth[1] > 0.) dm2 = Dynamic << 
188       while (dm1 + dm2 > parentmass)  // Loop  << 
189       {                                        << 
190         dm1 = DynamicalMass(daughtermass[0], d << 
191         dm2 = DynamicalMass(daughtermass[1], d << 
192       }                                        << 
193       daughtermass[0] = dm1;                   << 
194       daughtermass[1] = dm2;                   << 
195     }                                          << 
196   }                                            << 
197   else {                                       << 
198     // use given daughter mass;                << 
199     daughtermass[0] = givenDaughterMasses[0];  << 
200     daughtermass[1] = givenDaughterMasses[1];  << 
201   }                                            << 
202   if (parentmass < daughtermass[0] + daughterm << 
203 #ifdef G4VERBOSE                               << 
204     if (GetVerboseLevel() > 0) {               << 
205       G4cout << "G4PhaseSpaceDecayChannel::Two << 
206              << "Sum of daughter mass is large << 
207       G4cout << "Parent :" << G4MT_parent->Get << 
208              << current_parent_mass.Get() / Ge << 
209       G4cout << "Daughter 1 :" << G4MT_daughte << 
210              << daughtermass[0] / GeV << G4end << 
211       G4cout << "Daughter 2:" << G4MT_daughter << 
212              << daughtermass[1] / GeV << G4end << 
213       if (useGivenDaughterMass) {              << 
214         G4cout << "Daughter Mass is given." << << 
215       }                                        << 
216     }                                             180     }
217 #endif                                            181 #endif
218     G4Exception("G4PhaseSpaceDecayChannel::Two << 182     G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
219                 "Cannot create decay products: << 183                 "can not create decay products", JustWarning,
220                  larger than parent mass!");   << 184                 "sum of daughter mass is larger than parent mass");
221     return products;                              185     return products;
222   }                                               186   }
223                                                   187 
224   // calculate daughter momentum               << 188   G4double costheta = 2.*G4UniformRand()-1.0;
225   G4double daughtermomentum = Pmx(parentmass,  << 189   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
                                                   >> 190   G4double phi  = twopi*G4UniformRand()*rad;
                                                   >> 191   G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
226                                                   192 
227   G4double costheta = 2. * G4UniformRand() - 1 << 193   //create daughter G4DynamicParticle 
228   G4double sintheta = std::sqrt((1.0 - costhet << 194   G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
229   G4double phi = twopi * G4UniformRand() * rad << 
230   G4ThreeVector direction(sintheta * std::cos( << 
231                                                << 
232   // create daughter G4DynamicParticle         << 
233   G4double Ekin = std::sqrt(daughtermomentum * << 
234                   - daughtermass[0];           << 
235   auto daughterparticle =                      << 
236     new G4DynamicParticle(G4MT_daughters[0], d << 
237   products->PushProducts(daughterparticle);       195   products->PushProducts(daughterparticle);
238   Ekin = std::sqrt(daughtermomentum * daughter << 196   daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
239          - daughtermass[1];                    << 
240   daughterparticle =                           << 
241     new G4DynamicParticle(G4MT_daughters[1], - << 
242   products->PushProducts(daughterparticle);       197   products->PushProducts(daughterparticle);
243                                                   198 
244 #ifdef G4VERBOSE                                  199 #ifdef G4VERBOSE
245   if (GetVerboseLevel() > 1) {                 << 200   if (GetVerboseLevel()>1) {
246     G4cout << "G4PhaseSpaceDecayChannel::TwoBo << 201      G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
247     G4cout << " Create decay products in rest  << 202      G4cout << "  create decay products in rest frame " <<G4endl;
248     products->DumpInfo();                      << 203      products->DumpInfo();
249   }                                               204   }
250 #endif                                            205 #endif
251   return products;                                206   return products;
252 }                                                 207 }
253                                                   208 
254 G4DecayProducts* G4PhaseSpaceDecayChannel::Thr << 209 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
                                                   >> 210 // algorism of this code is originally written in GDECA3 of GEANT3
255 {                                                 211 {
256   // Algorithm of this code originally written << 
257                                                << 
258 #ifdef G4VERBOSE                                  212 #ifdef G4VERBOSE
259   if (GetVerboseLevel() > 1) G4cout << "G4Phas << 213   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
260 #endif                                            214 #endif
261   // parent mass                                  215   // parent mass
262   G4double parentmass = current_parent_mass.Ge << 216   G4double parentmass = parent_mass;
263   // daughters'mass                            << 217   //daughters'mass
264   G4double daughtermass[3], daughterwidth[3];  << 218   G4double daughtermass[3]; 
265   G4double sumofdaughtermass = 0.0;               219   G4double sumofdaughtermass = 0.0;
266   G4double sumofdaughterwidthsq = 0.0;         << 220   for (G4int index=0; index<3; index++){
267   G4bool withWidth = false;                    << 221     daughtermass[index] = daughters_mass[index];
268   for (G4int index = 0; index < 3; ++index) {  << 
269     daughtermass[index] = G4MT_daughters_mass[ << 
270     sumofdaughtermass += daughtermass[index];     222     sumofdaughtermass += daughtermass[index];
271     daughterwidth[index] = G4MT_daughters_widt << 
272     sumofdaughterwidthsq += daughterwidth[inde << 
273     withWidth = withWidth || (daughterwidth[in << 
274   }                                               223   }
275                                                << 224   
276   // create parent G4DynamicParticle at rest   << 225    //create parent G4DynamicParticle at rest
277   G4ThreeVector dummy;                            226   G4ThreeVector dummy;
278   auto parentparticle = new G4DynamicParticle( << 227   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
                                                   >> 228 
279                                                   229 
280   // create G4Decayproducts                    << 230   //create G4Decayproducts
281   auto products = new G4DecayProducts(*parentp << 231   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
282   delete parentparticle;                          232   delete parentparticle;
283                                                   233 
284   if (!useGivenDaughterMass) {                 << 234   if (sumofdaughtermass >parentmass) {
285     if (withWidth) {                           << 235 #ifdef G4VERBOSE
286       G4double maxDev = (parentmass - sumofdau << 236     if (GetVerboseLevel()>0) {
287       if (maxDev <= -1.0 * rangeMass) {        << 237       G4cerr << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt " 
288 #ifdef G4VERBOSE                               << 238              << "sum of daughter mass is larger than parent mass" << G4endl;
289         if (GetVerboseLevel() > 0) {           << 239       G4cerr << "parent :" << parent->GetParticleName() << "  " << parent_mass/GeV << G4endl;
290           G4cout << "G4PhaseSpaceDecayChannel: << 240       G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << "  " << daughtermass[0]/GeV << G4endl;
291                  << "Sum of daughter mass is l << 241       G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << "  " << daughtermass[1]/GeV << G4endl;
292           G4cout << "Parent :" << G4MT_parent- << 242       G4cerr << "daughter 3:" << daughters[2]->GetParticleName() << "  " << daughtermass[2]/GeV << G4endl;
293                  << current_parent_mass.Get()  << 
294           G4cout << "Daughter 1 :" << G4MT_dau << 
295                  << daughtermass[0] / GeV << G << 
296           G4cout << "Daughter 2:" << G4MT_daug << 
297                  << daughtermass[1] / GeV << G << 
298           G4cout << "Daughter 3:" << G4MT_daug << 
299                  << daughtermass[2] / GeV << G << 
300         }                                      << 
301 #endif                                         << 
302         G4Exception("G4PhaseSpaceDecayChannel: << 
303                     "Cannot create decay produ << 
304                      is larger than parent mas << 
305         return products;                       << 
306       }                                        << 
307       G4double dm1 = daughtermass[0];          << 
308       if (daughterwidth[0] > 0.) dm1 = Dynamic << 
309       G4double dm2 = daughtermass[1];          << 
310       if (daughterwidth[1] > 0.) dm2 = Dynamic << 
311       G4double dm3 = daughtermass[2];          << 
312       if (daughterwidth[2] > 0.) dm3 = Dynamic << 
313       while (dm1 + dm2 + dm3 > parentmass)  // << 
314       {                                        << 
315         dm1 = DynamicalMass(daughtermass[0], d << 
316         dm2 = DynamicalMass(daughtermass[1], d << 
317         dm3 = DynamicalMass(daughtermass[2], d << 
318       }                                        << 
319       daughtermass[0] = dm1;                   << 
320       daughtermass[1] = dm2;                   << 
321       daughtermass[2] = dm3;                   << 
322       sumofdaughtermass = dm1 + dm2 + dm3;     << 
323     }                                          << 
324   }                                            << 
325   else {                                       << 
326     // use given daughter mass;                << 
327     daughtermass[0] = givenDaughterMasses[0];  << 
328     daughtermass[1] = givenDaughterMasses[1];  << 
329     daughtermass[2] = givenDaughterMasses[2];  << 
330     sumofdaughtermass = daughtermass[0] + daug << 
331   }                                            << 
332                                                << 
333   if (sumofdaughtermass > parentmass) {        << 
334 #ifdef G4VERBOSE                               << 
335     if (GetVerboseLevel() > 0) {               << 
336       G4cout << "G4PhaseSpaceDecayChannel::Thr << 
337              << "Sum of daughter mass is large << 
338       G4cout << "Parent :" << G4MT_parent->Get << 
339              << current_parent_mass.Get() / Ge << 
340       G4cout << "Daughter 1 :" << G4MT_daughte << 
341              << daughtermass[0] / GeV << G4end << 
342       G4cout << "Daughter 2:" << G4MT_daughter << 
343              << daughtermass[1] / GeV << G4end << 
344       G4cout << "Daughter 3:" << G4MT_daughter << 
345              << daughtermass[2] / GeV << G4end << 
346     }                                          << 
347     if (useGivenDaughterMass) {                << 
348       G4cout << "Daughter Mass is given." << G << 
349     }                                             243     }
350 #endif                                            244 #endif
351     G4Exception("G4PhaseSpaceDecayChannel::Thr << 245     G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
352                 "Can not create decay products << 246                 "can not create decay products", JustWarning,
353                  is larger than parent mass!") << 247                 "sum of daughter mass is larger than parent mass");
354     return products;                              248     return products;
355   }                                               249   }
356                                                   250 
357   // calculate daughter momentum               << 251 
358   // Generate two                              << 252 
                                                   >> 253   //calculate daughter momentum
                                                   >> 254   //  Generate two 
359   G4double rd1, rd2, rd;                          255   G4double rd1, rd2, rd;
360   G4double daughtermomentum[3];                   256   G4double daughtermomentum[3];
361   G4double momentummax = 0.0, momentumsum = 0. << 257   G4double momentummax=0.0, momentumsum = 0.0;
362   G4double energy;                                258   G4double energy;
363   const std::size_t MAX_LOOP = 10000;          << 
364                                                   259 
365   for (std::size_t loop_counter = 0; loop_coun << 260   do {
366     rd1 = G4UniformRand();                        261     rd1 = G4UniformRand();
367     rd2 = G4UniformRand();                        262     rd2 = G4UniformRand();
368     if (rd2 > rd1) {                              263     if (rd2 > rd1) {
369       rd = rd1;                                << 264       rd  = rd1;
370       rd1 = rd2;                                  265       rd1 = rd2;
371       rd2 = rd;                                   266       rd2 = rd;
372     }                                          << 267     } 
373     momentummax = 0.0;                            268     momentummax = 0.0;
374     momentumsum = 0.0;                            269     momentumsum = 0.0;
375     // daughter 0                                 270     // daughter 0
376     energy = rd2 * (parentmass - sumofdaughter << 271     energy = rd2*(parentmass - sumofdaughtermass);
377     daughtermomentum[0] = std::sqrt(energy * e << 272     daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
378     if (daughtermomentum[0] > momentummax) mom << 273     if ( daughtermomentum[0] >momentummax )momentummax =  daughtermomentum[0];
379     momentumsum += daughtermomentum[0];        << 274     momentumsum  +=  daughtermomentum[0];
380     // daughter 1                                 275     // daughter 1
381     energy = (1. - rd1) * (parentmass - sumofd << 276     energy = (1.-rd1)*(parentmass - sumofdaughtermass);
382     daughtermomentum[1] = std::sqrt(energy * e << 277     daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
383     if (daughtermomentum[1] > momentummax) mom << 278     if ( daughtermomentum[1] >momentummax )momentummax =  daughtermomentum[1];
384     momentumsum += daughtermomentum[1];        << 279     momentumsum  +=  daughtermomentum[1];
385     // daughter 2                                 280     // daughter 2
386     energy = (rd1 - rd2) * (parentmass - sumof << 281     energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
387     daughtermomentum[2] = std::sqrt(energy * e << 282     daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
388     if (daughtermomentum[2] > momentummax) mom << 283     if ( daughtermomentum[2] >momentummax )momentummax =  daughtermomentum[2];
389     momentumsum += daughtermomentum[2];        << 284     momentumsum  +=  daughtermomentum[2];
390     if (momentummax <= momentumsum - momentumm << 285   } while (momentummax >  momentumsum - momentummax );
391   }                                            << 
392                                                << 
393   // output message                               286   // output message
394 #ifdef G4VERBOSE                                  287 #ifdef G4VERBOSE
395   if (GetVerboseLevel() > 1) {                 << 288   if (GetVerboseLevel()>1) {
396     G4cout << "     daughter 0:" << daughtermo << 289     G4cout << "     daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
397     G4cout << "     daughter 1:" << daughtermo << 290     G4cout << "     daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
398     G4cout << "     daughter 2:" << daughtermo << 291     G4cout << "     daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
399     G4cout << "   momentum sum:" << momentumsu << 292     G4cout << "   momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
400   }                                               293   }
401 #endif                                            294 #endif
402                                                   295 
403   // create daughter G4DynamicParticle         << 296   //create daughter G4DynamicParticle 
404   G4double costheta, sintheta, phi, sinphi, co << 297   G4double costheta, sintheta, phi, sinphi, cosphi; 
405   G4double costhetan, sinthetan, phin, sinphin << 298   G4double costhetan, sinthetan, phin, sinphin, cosphin; 
406   costheta = 2. * G4UniformRand() - 1.0;       << 299   costheta = 2.*G4UniformRand()-1.0;
407   sintheta = std::sqrt((1.0 - costheta) * (1.0 << 300   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
408   phi = twopi * G4UniformRand() * rad;         << 301   phi  = twopi*G4UniformRand()*rad;
409   sinphi = std::sin(phi);                         302   sinphi = std::sin(phi);
410   cosphi = std::cos(phi);                         303   cosphi = std::cos(phi);
411                                                << 304   G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
412   G4ThreeVector direction0(sintheta * cosphi,  << 305   G4DynamicParticle * daughterparticle 
413   G4double Ekin =                              << 306          = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
414     std::sqrt(daughtermomentum[0] * daughtermo << 
415     - daughtermass[0];                         << 
416   auto daughterparticle =                      << 
417     new G4DynamicParticle(G4MT_daughters[0], d << 
418   products->PushProducts(daughterparticle);       307   products->PushProducts(daughterparticle);
419                                                   308 
420   costhetan = (daughtermomentum[1] * daughterm << 309   costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
421                - daughtermomentum[0] * daughte << 310   sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
422               / (2.0 * daughtermomentum[2] * d << 311   phin  = twopi*G4UniformRand()*rad;
423   sinthetan = std::sqrt((1.0 - costhetan) * (1 << 
424   phin = twopi * G4UniformRand() * rad;        << 
425   sinphin = std::sin(phin);                       312   sinphin = std::sin(phin);
426   cosphin = std::cos(phin);                       313   cosphin = std::cos(phin);
427   G4ThreeVector direction2;                       314   G4ThreeVector direction2;
428   direction2.setX(sinthetan * cosphin * costhe << 315   direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 
429                   + costhetan * sintheta * cos << 316   direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 
430   direction2.setY(sinthetan * cosphin * costhe << 317   direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
431                   + costhetan * sintheta * sin << 318   daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
432   direction2.setZ(-sinthetan * cosphin * sinth << 
433   G4ThreeVector pmom = daughtermomentum[2] * d << 
434   Ekin = std::sqrt(pmom.mag2() + daughtermass[ << 
435   daughterparticle =                           << 
436     new G4DynamicParticle(G4MT_daughters[2], p << 
437   products->PushProducts(daughterparticle);       319   products->PushProducts(daughterparticle);
438                                                   320 
439   pmom = (direction0 * daughtermomentum[0] + d << 321   daughterparticle = 
440          * (-1.0);                             << 322        new G4DynamicParticle( 
441   Ekin = std::sqrt(pmom.mag2() + daughtermass[ << 323           daughters[1],
442   daughterparticle =                           << 324          (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)   
443     new G4DynamicParticle(G4MT_daughters[1], p << 325     );
444   products->PushProducts(daughterparticle);       326   products->PushProducts(daughterparticle);
445                                                << 327   
446 #ifdef G4VERBOSE                                  328 #ifdef G4VERBOSE
447   if (GetVerboseLevel() > 1) {                 << 329   if (GetVerboseLevel()>1) {
448     G4cout << "G4PhaseSpaceDecayChannel::Three << 330      G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
449     G4cout << " create decay products in rest  << 331      G4cout << "  create decay products in rest frame " <<G4endl;
450     products->DumpInfo();                      << 332      products->DumpInfo();
451   }                                               333   }
452 #endif                                            334 #endif
453   return products;                                335   return products;
454 }                                                 336 }
455                                                   337 
456 G4DecayProducts* G4PhaseSpaceDecayChannel::Man << 338 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
                                                   >> 339 // algorithm of this code is originally written in FORTRAN by M.Asai
                                                   >> 340 //******************************************************************
                                                   >> 341 //  NBODY
                                                   >> 342 //   N-body phase space Monte-Carlo generator
                                                   >> 343 //  Makoto Asai 
                                                   >> 344 //   Hiroshima Institute of Technology
                                                   >> 345 //    (asai@kekvax.kek.jp)
                                                   >> 346 //  Revised release : 19/Apr/1995
                                                   >> 347 //
457 {                                                 348 {
458   // Algorithm of this code originally written << 
459   // ***************************************** << 
460   // NBODY - N-body phase space Monte-Carlo ge << 
461   // Makoto Asai - Hiroshima Institute of Tech << 
462   // Revised release : 19/Apr/1995             << 
463                                                << 
464   G4int index, index2;                            349   G4int index, index2;
465                                                   350 
                                                   >> 351   //return value
                                                   >> 352   G4DecayProducts *products;
                                                   >> 353 
466 #ifdef G4VERBOSE                                  354 #ifdef G4VERBOSE
467   if (GetVerboseLevel() > 1) G4cout << "G4Phas << 355   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
468 #endif                                            356 #endif
469                                                << 
470   // parent mass                                  357   // parent mass
471   G4double parentmass = current_parent_mass.Ge << 358   G4double parentmass = parent_mass;
472                                                << 359   //daughters'mass
473   // parent particle                           << 360   G4double *daughtermass = new G4double[numberOfDaughters]; 
474   G4ThreeVector dummy;                         << 
475   auto parentparticle = new G4DynamicParticle( << 
476                                                << 
477   // create G4Decayproducts                    << 
478   auto products = new G4DecayProducts(*parentp << 
479   delete parentparticle;                       << 
480                                                << 
481   // daughters'mass                            << 
482   auto daughtermass = new G4double[numberOfDau << 
483                                                << 
484   G4double sumofdaughtermass = 0.0;               361   G4double sumofdaughtermass = 0.0;
485   for (index = 0; index < numberOfDaughters; + << 362   for (index=0; index<numberOfDaughters; index++){
486     if (!useGivenDaughterMass) {               << 363     daughtermass[index] = daughters_mass[index];
487       daughtermass[index] = G4MT_daughters_mas << 
488     }                                          << 
489     else {                                     << 
490       daughtermass[index] = givenDaughterMasse << 
491     }                                          << 
492     sumofdaughtermass += daughtermass[index];     364     sumofdaughtermass += daughtermass[index];
493   }                                               365   }
494   if (sumofdaughtermass > parentmass) {        << 366   
495 #ifdef G4VERBOSE                               << 367   //Calculate daughter momentum
496     if (GetVerboseLevel() > 0) {               << 368   G4double *daughtermomentum = new G4double[numberOfDaughters];
497       G4cout << "G4PhaseSpaceDecayChannel::Man << 369   G4ThreeVector direction;  
498              << "Sum of daughter mass is large << 370   G4DynamicParticle **daughterparticle;
499       G4cout << "Parent :" << G4MT_parent->Get << 371   G4double *sm = new G4double[numberOfDaughters];
500              << current_parent_mass.Get() / Ge << 
501       G4cout << "Daughter 1 :" << G4MT_daughte << 
502              << daughtermass[0] / GeV << G4end << 
503       G4cout << "Daughter 2:" << G4MT_daughter << 
504              << daughtermass[1] / GeV << G4end << 
505       G4cout << "Daughter 3:" << G4MT_daughter << 
506              << daughtermass[2] / GeV << G4end << 
507       G4cout << "Daughter 4:" << G4MT_daughter << 
508              << daughtermass[3] / GeV << G4end << 
509       G4cout << "Daughter 5:" << G4MT_daughter << 
510              << daughtermass[4] / GeV << G4end << 
511     }                                          << 
512 #endif                                         << 
513     G4Exception("G4PhaseSpaceDecayChannel::Man << 
514                 "Can not create decay products << 
515                  is larger than parent mass!") << 
516     delete[] daughtermass;                     << 
517     return products;                           << 
518   }                                            << 
519                                                << 
520   // Calculate daughter momentum               << 
521   auto daughtermomentum = new G4double[numberO << 
522   G4ThreeVector direction;                     << 
523   G4DynamicParticle** daughterparticle;        << 
524   auto sm = new G4double[numberOfDaughters];   << 
525   G4double tmas;                                  372   G4double tmas;
526   G4double weight = 1.0;                          373   G4double weight = 1.0;
527   G4int numberOfTry = 0;                       << 374   G4int    numberOfTry = 0;
528   const std::size_t MAX_LOOP = 10000;          << 375  
529                                                << 376   do {
530   for (std::size_t loop_counter = 0; loop_coun << 377     //Generate rundom number in descending order 
531     // Generate rundom number in descending or << 
532     G4double temp;                                378     G4double temp;
533     auto rd = new G4double[numberOfDaughters]; << 379     G4double *rd = new G4double[numberOfDaughters];
534     rd[0] = 1.0;                                  380     rd[0] = 1.0;
535     for (index = 1; index < numberOfDaughters  << 381     for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand(); 
536       rd[index] = G4UniformRand();             << 382     rd[ numberOfDaughters -1] = 0.0;
537     }                                          << 383     for(index =1; index < numberOfDaughters -1; index++) {
538     rd[numberOfDaughters - 1] = 0.0;           << 384       for(index2 = index+1; index2 < numberOfDaughters; index2++) {
539     for (index = 1; index < numberOfDaughters  << 385         if (rd[index] < rd[index2]){
540       for (index2 = index + 1; index2 < number << 386           temp        = rd[index];
541         if (rd[index] < rd[index2]) {          << 387           rd[index]  = rd[index2];
542           temp = rd[index];                    << 
543           rd[index] = rd[index2];              << 
544           rd[index2] = temp;                      388           rd[index2] = temp;
545         }                                         389         }
546       }                                           390       }
547     }                                             391     }
548                                                << 392     
549     // calculate virtual mass                  << 393     //calcurate virtual mass 
550     tmas = parentmass - sumofdaughtermass;     << 394     tmas = parentmass -  sumofdaughtermass;
551     temp = sumofdaughtermass;                  << 395     temp =  sumofdaughtermass; 
552     for (index = 0; index < numberOfDaughters; << 396     for(index =0; index < numberOfDaughters; index++) {
553       sm[index] = rd[index] * tmas + temp;     << 397       sm[index] = rd[index]*tmas + temp;
554       temp -= daughtermass[index];                398       temp -= daughtermass[index];
555       if (GetVerboseLevel() > 1) {             << 399       if (GetVerboseLevel()>1) {
556         G4cout << index << "  rundom number:"  << 400         G4cout << index << "  rundom number:" << rd[index]; 
557         G4cout << "   virtual mass:" << sm[ind << 401         G4cout << "   virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl; 
558       }                                           402       }
559     }                                             403     }
560     delete[] rd;                               << 404     delete [] rd;
561                                                   405 
562     // Calculate daughter momentum             << 406     //Calculate daughter momentum
563     weight = 1.0;                                 407     weight = 1.0;
564     G4bool smOK = true;                        << 408     index =numberOfDaughters-1;
565     for (index = 0; index < numberOfDaughters  << 409     daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
566       smOK = (sm[index] - daughtermass[index]  << 
567     }                                          << 
568     if (!smOK) continue;                       << 
569                                                << 
570     index = numberOfDaughters - 1;             << 
571     daughtermomentum[index] = Pmx(sm[index - 1 << 
572 #ifdef G4VERBOSE                                  410 #ifdef G4VERBOSE
573     if (GetVerboseLevel() > 1) {               << 411     if (GetVerboseLevel()>1) {
574       G4cout << "     daughter " << index << "    412       G4cout << "     daughter " << index << ":" << *daughters_name[index];
575       G4cout << " momentum:" << daughtermoment << 413       G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
576     }                                             414     }
577 #endif                                            415 #endif
578     for (index = numberOfDaughters - 2; index  << 416     for(index =numberOfDaughters-2; index>=0; index--) {
579       // calculate                             << 417       // calculate 
580       daughtermomentum[index] = Pmx(sm[index], << 418       daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
581       if (daughtermomentum[index] < 0.0) {     << 419       if(daughtermomentum[index] < 0.0) {
582         // !!! illegal momentum !!!               420         // !!! illegal momentum !!!
583 #ifdef G4VERBOSE                                  421 #ifdef G4VERBOSE
584         if (GetVerboseLevel() > 0) {           << 422         if (GetVerboseLevel()>0) {
585           G4cout << "G4PhaseSpaceDecayChannel: << 423           G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
586           G4cout << "     Cannot calculate dau << 424           G4cout << "     can not calculate daughter momentum " <<G4endl;
587           G4cout << "     parent:" << *parent_    425           G4cout << "     parent:" << *parent_name;
588           G4cout << " mass:" << parentmass / G << 426           G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
589           G4cout << "     daughter " << index     427           G4cout << "     daughter " << index << ":" << *daughters_name[index];
590           G4cout << " mass:" << daughtermass[i << 428           G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
591           G4cout << " mass:" << daughtermoment << 429           G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
592           if (useGivenDaughterMass) {          << 
593             G4cout << "Daughter Mass is given. << 
594           }                                    << 
595         }                                         430         }
596 #endif                                            431 #endif
597         delete[] sm;                           << 432   delete [] sm;
598         delete[] daughtermass;                 << 433   delete [] daughtermass;
599         delete[] daughtermomentum;             << 434   delete [] daughtermomentum;
600         delete products;                       << 435 
601         G4Exception("G4PhaseSpaceDecayChannel: << 436   G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
602                     "Can not create decay prod << 437                 "can not create decay products", JustWarning,
603                      is larger than parent mas << 438                 "sum of daughter mass is larger than parent mass");
604         return nullptr;  // Error detection    << 439 
605       }                                        << 440   return 0;   // Error detection
606                                                << 441 
607       // calculate weight of this events       << 442       } else {
608       weight *= daughtermomentum[index] / sm[i << 443   // calculate weight of this events
609 #ifdef G4VERBOSE                               << 444         weight *=  daughtermomentum[index]/sm[index];
610       if (GetVerboseLevel() > 1) {             << 445         if (GetVerboseLevel()>1) {
611         G4cout << "     daughter " << index << << 446           G4cout << "     daughter " << index << ":" << *daughters_name[index];
612         G4cout << " momentum:" << daughtermome << 447           G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
                                                   >> 448         }
613       }                                           449       }
614 #endif                                         << 
615     }                                             450     }
616                                                << 451     if (GetVerboseLevel()>1) {
617 #ifdef G4VERBOSE                               << 452       G4cout << "    weight: " << weight <<G4endl;
618     if (GetVerboseLevel() > 1) {               << 
619       G4cout << "    weight: " << weight << G4 << 
620     }                                             453     }
621 #endif                                         << 454     
622                                                << 
623     // exit if number of Try exceeds 100          455     // exit if number of Try exceeds 100
624     if (++numberOfTry > 100) {                 << 456     if (numberOfTry++ >100) {
625 #ifdef G4VERBOSE                               << 457       if (GetVerboseLevel()>0) {
626       if (GetVerboseLevel() > 0) {             << 458         G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
627         G4cout << "G4PhaseSpaceDecayChannel::M << 459   G4cout << " can not determine Decay Kinematics " << G4endl;
628         G4cout << "Cannot determine Decay Kine << 
629         G4cout << "parent : " << *parent_name  << 
630         G4cout << "daughters : ";              << 
631         for (index = 0; index < numberOfDaught << 
632           G4cout << *daughters_name[index] <<  << 
633         }                                      << 
634         G4cout << G4endl;                      << 
635       }                                           460       }
636 #endif                                         << 461       delete [] sm;
637       G4Exception("G4PhaseSpaceDecayChannel::M << 462       delete [] daughtermass;
638                   "Cannot decay: Decay Kinemat << 463       delete [] daughtermomentum;
639                                                << 464       return 0;  // Error detection
640       delete[] sm;                             << 
641       delete[] daughtermass;                   << 
642       delete[] daughtermomentum;               << 
643       delete products;                         << 
644       return nullptr;  // Error detection      << 
645     }                                             465     }
646     if (weight < G4UniformRand()) break;       << 466   } while ( weight > G4UniformRand());
647   }                                            << 
648                                                   467 
649 #ifdef G4VERBOSE                                  468 #ifdef G4VERBOSE
650   if (GetVerboseLevel() > 1) {                 << 469   if (GetVerboseLevel()>1) {
651     G4cout << "Start calculation of daughters  << 470       G4cout << "Start calulation of daughters momentum vector "<<G4endl;
652   }                                               471   }
653 #endif                                         << 472 #endif  
654                                                << 
655   G4double costheta, sintheta, phi;               473   G4double costheta, sintheta, phi;
656   G4double beta;                                  474   G4double beta;
657   daughterparticle = new G4DynamicParticle*[nu    475   daughterparticle = new G4DynamicParticle*[numberOfDaughters];
658                                                   476 
659   index = numberOfDaughters - 2;               << 477   index = numberOfDaughters -2;
660   costheta = 2. * G4UniformRand() - 1.0;       << 478   costheta = 2.*G4UniformRand()-1.0;
661   sintheta = std::sqrt((1.0 - costheta) * (1.0 << 479   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
662   phi = twopi * G4UniformRand() * rad;         << 480   phi  = twopi*G4UniformRand()*rad;
663   direction.setZ(costheta);                       481   direction.setZ(costheta);
664   direction.setY(sintheta * std::sin(phi));    << 482   direction.setY(sintheta*std::sin(phi));
665   direction.setX(sintheta * std::cos(phi));    << 483   direction.setX(sintheta*std::cos(phi));
666   daughterparticle[index] =                    << 484   daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
667     new G4DynamicParticle(G4MT_daughters[index << 485   daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
668   daughterparticle[index + 1] =                << 486 
669     new G4DynamicParticle(G4MT_daughters[index << 487   for (index = numberOfDaughters -3;  index >= 0; index--) {
670                                                << 488     //calculate momentum direction
671   for (index = numberOfDaughters - 3; index >= << 489     costheta = 2.*G4UniformRand()-1.0;
672     // calculate momentum direction            << 490     sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
673     costheta = 2. * G4UniformRand() - 1.0;     << 491     phi  = twopi*G4UniformRand()*rad;
674     sintheta = std::sqrt((1.0 - costheta) * (1 << 
675     phi = twopi * G4UniformRand() * rad;       << 
676     direction.setZ(costheta);                     492     direction.setZ(costheta);
677     direction.setY(sintheta * std::sin(phi));  << 493     direction.setY(sintheta*std::sin(phi));
678     direction.setX(sintheta * std::cos(phi));  << 494     direction.setX(sintheta*std::cos(phi));
679                                                   495 
680     // boost already created particles         << 496     // boost already created particles 
681     beta = daughtermomentum[index];               497     beta = daughtermomentum[index];
682     beta /=                                    << 498     beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
683       std::sqrt(daughtermomentum[index] * daug << 499     for (G4int index2 = index+1; index2<numberOfDaughters; index2++) {
684     for (index2 = index + 1; index2 < numberOf << 
685       G4LorentzVector p4;                         500       G4LorentzVector p4;
686       // make G4LorentzVector for secondaries     501       // make G4LorentzVector for secondaries
687       p4 = daughterparticle[index2]->Get4Momen    502       p4 = daughterparticle[index2]->Get4Momentum();
688                                                   503 
689       // boost secondaries to  new frame       << 504       // boost secondaries to  new frame 
690       p4.boost(direction.x() * beta, direction << 505       p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
691                                                   506 
692       // change energy/momentum                   507       // change energy/momentum
693       daughterparticle[index2]->Set4Momentum(p    508       daughterparticle[index2]->Set4Momentum(p4);
694     }                                             509     }
695     // create daughter G4DynamicParticle       << 510     //create daughter G4DynamicParticle 
696     daughterparticle[index] =                  << 511     daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
697       new G4DynamicParticle(G4MT_daughters[ind << 
698   }                                               512   }
699                                                   513 
700   // add daughters to G4Decayproducts          << 514   //create G4Decayproducts
701   for (index = 0; index < numberOfDaughters; + << 515   G4DynamicParticle *parentparticle; 
                                                   >> 516   direction.setX(1.0);  direction.setY(0.0); direction.setZ(0.0);
                                                   >> 517   parentparticle = new G4DynamicParticle( parent, direction, 0.0);
                                                   >> 518   products = new G4DecayProducts(*parentparticle);
                                                   >> 519   delete parentparticle;
                                                   >> 520   for (index = 0; index<numberOfDaughters; index++) {
702     products->PushProducts(daughterparticle[in    521     products->PushProducts(daughterparticle[index]);
703   }                                               522   }
704                                                << 
705 #ifdef G4VERBOSE                                  523 #ifdef G4VERBOSE
706   if (GetVerboseLevel() > 1) {                 << 524   if (GetVerboseLevel()>1) { 
707     G4cout << "G4PhaseSpaceDecayChannel::ManyB << 525     G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
708     G4cout << " create decay products in rest  << 526     G4cout << "  create decay products in rest frame " << G4endl;
709     products->DumpInfo();                         527     products->DumpInfo();
710   }                                               528   }
711 #endif                                            529 #endif
712                                                << 530   delete [] daughterparticle;
713   delete[] daughterparticle;                   << 531   delete [] daughtermomentum;
714   delete[] daughtermomentum;                   << 532   delete [] daughtermass;
715   delete[] daughtermass;                       << 533   delete [] sm;
716   delete[] sm;                                 << 534   
717                                                << 
718   return products;                                535   return products;
719 }                                                 536 }
720                                                   537 
721 G4bool G4PhaseSpaceDecayChannel::SetDaughterMa << 
722 {                                              << 
723   for (G4int idx = 0; idx < numberOfDaughters; << 
724     givenDaughterMasses[idx] = masses[idx];    << 
725   }                                            << 
726   useGivenDaughterMass = true;                 << 
727   return useGivenDaughterMass;                 << 
728 }                                              << 
729                                                   538 
730 G4bool G4PhaseSpaceDecayChannel::SampleDaughte << 539 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
731 {                                                 540 {
732   useGivenDaughterMass = false;                << 541    // calcurate momentum of daughter particles in two-body decay
733   return useGivenDaughterMass;                 << 542    G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
                                                   >> 543    if (ppp>0) return std::sqrt(ppp);
                                                   >> 544    else       return -1.;
734 }                                                 545 }
735                                                   546 
736 G4bool G4PhaseSpaceDecayChannel::IsOKWithParen << 
737 {                                              << 
738   if (!useGivenDaughterMass) return G4VDecayCh << 
739                                                   547 
740   CheckAndFillParent();                        << 
741   CheckAndFillDaughters();                     << 
742                                                   548 
743   G4double sumOfDaughterMassMin = 0.0;         << 
744   for (G4int index = 0; index < numberOfDaught << 
745     sumOfDaughterMassMin += givenDaughterMasse << 
746   }                                            << 
747   return (parentMass >= sumOfDaughterMassMin); << 
748 }                                              << 
749                                                << 
750 G4double G4PhaseSpaceDecayChannel::Pmx(G4doubl << 
751 {                                              << 
752   // calcurate momentum of daughter particles  << 
753   G4double ppp = (e + p1 + p2) * (e + p1 - p2) << 
754   if (ppp > 0) return std::sqrt(ppp);          << 
755   return -1.;                                  << 
756 }                                              << 
757                                                   549