| Geant4 Cross Reference |
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 //
27 // Hadronic Process: Phase space decay for the Fermi BreakUp model
28 // by V. Lara
29 //
30 // Modifications:
31 // 01.04.2011 General cleanup by V.Ivanchenko:
32 // - IsotropicVector is inlined
33 // - Momentum computation return zero or positive value
34 // - DumpProblem method is added providing more information
35 // - Reduced usage of exotic std functions
36
37 #include "G4FermiPhaseSpaceDecay.hh"
38
39 #include "G4RandomDirection.hh"
40 #include "G4Pow.hh"
41
42 #include <CLHEP/Units/SystemOfUnits.h>
43 #include <CLHEP/Units/PhysicalConstants.h>
44 #include <CLHEP/Random/RandomEngine.h>
45
46 G4FermiPhaseSpaceDecay::G4FermiPhaseSpaceDecay()
47 {
48 g4calc = G4Pow::GetInstance();
49 }
50
51 G4FermiPhaseSpaceDecay::~G4FermiPhaseSpaceDecay()
52 {}
53
54 std::vector<G4LorentzVector*>* G4FermiPhaseSpaceDecay::Decay(G4double M,
55 const std::vector<G4double>& mr) const
56 // Calculates momentum for N fragments (Kopylov's method of sampling is used)
57 {
58 std::size_t N = mr.size();
59
60 std::vector<G4LorentzVector*>* P =
61 new std::vector<G4LorentzVector*>(N, nullptr);
62
63 G4double mtot = 0.0;
64 for(std::size_t k=0; k<N; ++k) { mtot += mr[k]; }
65
66 G4double mu = mtot;
67 G4double PFragMagCM = 0.0;
68
69 // Primary mass is above the sum of mass of components
70 G4double Mass = std::max(M, mtot + CLHEP::eV);
71 G4double T = Mass-mtot;
72
73 G4LorentzVector PFragCM(0.0,0.0,0.0,0.0);
74 G4LorentzVector PRestCM(0.0,0.0,0.0,0.0);
75 G4LorentzVector PRestLab(0.0,0.0,0.0,Mass);
76
77 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
78
79 for (G4int k = (G4int)N-1; k>0; --k)
80 {
81 mu -= mr[k];
82 if (k>1) { T *= BetaKopylov(k, rndmEngine); }
83 else { T = 0.0; }
84
85 G4double RestMass = mu + T;
86
87 PFragMagCM = PtwoBody(Mass,mr[k],RestMass);
88
89 // Create a unit vector with a random direction isotropically distributed
90 G4ThreeVector RandVector = PFragMagCM*G4RandomDirection();
91
92 PFragCM.setVect(RandVector);
93 PFragCM.setE(std::sqrt(PFragMagCM*PFragMagCM + mr[k]*mr[k]));
94
95 PRestCM.setVect(-RandVector);
96 PRestCM.setE(std::sqrt(PFragMagCM*PFragMagCM + RestMass*RestMass));
97
98 G4ThreeVector BoostV = PRestLab.boostVector();
99
100 PFragCM.boost(BoostV);
101 (*P)[k] = new G4LorentzVector(PFragCM);
102
103 PRestCM.boost(BoostV);
104 PRestLab = PRestCM;
105
106 Mass = RestMass;
107 }
108
109 (*P)[0] = new G4LorentzVector(PRestLab);
110
111 return P;
112 }
113
114 G4double G4FermiPhaseSpaceDecay::BetaKopylov(G4int K,
115 CLHEP::HepRandomEngine* rndmEngine) const
116 {
117 G4int N = 3*K - 5;
118 G4double xN = (G4double)N;
119 G4double xN1= (G4double)(N + 1);
120 G4double F;
121 // VI variant
122 G4double Fmax = std::sqrt(g4calc->powN(xN/xN1,N)/xN1);
123 G4double chi;
124 do {
125 chi = rndmEngine->flat();
126 F = std::sqrt(g4calc->powN(chi,N)*(1-chi));
127 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
128 } while ( Fmax*rndmEngine->flat() > F);
129 return chi;
130 }
131