Geant4  10.02.p01
G4FermiPhaseSpaceDecay.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4ExcitationHandler.hh,v 1.13 2010-11-17 16:20:31 vnivanch Exp $
27 //
28 // Hadronic Process: Phase space decay for the Fermi BreakUp model
29 // by V. Lara
30 //
31 // Modifications:
32 // 01.04.2011 General cleanup by V.Ivanchenko:
33 // - IsotropicVector is inlined
34 // - Momentum computation return zero or positive value
35 // - DumpProblem method is added providing more information
36 // - Reduced usage of exotic std functions
37 
39 #include "G4SystemOfUnits.hh"
40 #include "G4HadronicException.hh"
41 
43 {
45 }
46 
48 {}
49 
50 std::vector<G4LorentzVector*>*
52  const std::vector<G4double>& mr) const
53  // Calculates momentum for N fragments (Kopylov's method of sampling is used)
54 {
55  size_t N = mr.size();
56 
57  std::vector<G4LorentzVector*>* P =
58  new std::vector<G4LorentzVector*>(N, 0);
59 
60  G4double mtot = 0.0;
61  for(size_t k=0; k<N; ++k) { mtot += mr[k]; }
62  G4double mu = mtot;
63  G4double PFragMagCM = 0.0;
64  G4double Mass = M;
65  G4double T = Mass-mtot;
66  G4LorentzVector PFragCM(0.0,0.0,0.0,0.0);
67  G4LorentzVector PRestCM(0.0,0.0,0.0,0.0);
68  G4LorentzVector PRestLab(0.0,0.0,0.0,Mass);
69 
70  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
71 
72  for (size_t k = N-1; k>0; --k)
73  {
74  mu -= mr[k];
75  if (k>1) { T *= BetaKopylov(k, rndmEngine); }
76  else { T = 0.0; }
77 
78  G4double RestMass = mu + T;
79 
80  PFragMagCM = PtwoBody(Mass,mr[k],RestMass);
81 
82  // Create a unit vector with a random direction isotropically distributed
83  G4ThreeVector RandVector(IsotropicVector(PFragMagCM, rndmEngine));
84 
85  PFragCM.setVect(RandVector);
86  PFragCM.setE(std::sqrt(PFragMagCM*PFragMagCM + mr[k]*mr[k]));
87 
88  PRestCM.setVect(-RandVector);
89  PRestCM.setE(std::sqrt(PFragMagCM*PFragMagCM + RestMass*RestMass));
90 
91 
92  G4ThreeVector BoostV = PRestLab.boostVector();
93 
94  PFragCM.boost(BoostV);
95  PRestCM.boost(BoostV);
96  PRestLab = PRestCM;
97 
98  (*P)[k] = new G4LorentzVector(PFragCM);
99 
100  Mass = RestMass;
101  }
102 
103  (*P)[0] = new G4LorentzVector(PRestLab);
104 
105  return P;
106 }
107 
108 void
110  G4double P) const
111 {
112  G4cout << "G4FermiPhaseSpaceDecay: problem of decay of M(GeV)= " << E/GeV
113  << " on M1(GeV)= " << P1/GeV << " and M2(GeV)= " << P2/GeV
114  << " P(MeV)= " << P/MeV << " < 0" << G4endl;
115  // exception only if the problem is numerically significant
116  if(P < -CLHEP::eV) {
117  throw G4HadronicException(__FILE__, __LINE__,"Error in decay kinematics");
118  }
119 }
120 
121 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static const double MeV
Definition: G4SIunits.hh:211
static const G4double * P1[nN]
CLHEP::Hep3Vector G4ThreeVector
static double P[]
void DumpProblem(G4double E, G4double P1, G4double P2, G4double P) const
G4GLOB_DLL std::ostream G4cout
static const double GeV
Definition: G4SIunits.hh:214
static const double eV
Definition: G4SIunits.hh:212
static const G4double * P2[nN]
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
std::vector< G4LorentzVector * > * KopylovNBodyDecay(G4double, const std::vector< G4double > &) const
CLHEP::HepLorentzVector G4LorentzVector