| 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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
28 //
29 #ifndef G4ParticleHPNBodyPhaseSpace_h
30 #define G4ParticleHPNBodyPhaseSpace_h 1
31
32 #include "G4Neutron.hh"
33 #include "G4Pow.hh"
34 #include "G4ReactionProduct.hh"
35 #include "G4VParticleHPEnergyAngular.hh"
36 #include "G4ios.hh"
37 #include "globals.hh"
38
39 #include <CLHEP/Units/PhysicalConstants.h>
40
41 #include <fstream>
42
43 class G4ParticleHPNBodyPhaseSpace : public G4VParticleHPEnergyAngular
44 {
45 public:
46 G4ParticleHPNBodyPhaseSpace()
47 {
48 theTotalMass = 0.0;
49 theTotalCount = 0;
50 }
51 ~G4ParticleHPNBodyPhaseSpace() override = default;
52
53 public:
54 void Init(G4double aMass, G4int aCount)
55 {
56 theTotalMass = aMass;
57 theTotalCount = aCount;
58 }
59
60 void Init(std::istream& aDataFile) override
61 {
62 aDataFile >> theTotalMass >> theTotalCount;
63 theTotalMass *= G4Neutron::Neutron()->GetPDGMass();
64 }
65
66 G4ReactionProduct* Sample(G4double anEnergy, G4double massCode, G4double mass) override;
67
68 private:
69 inline G4double Prob(G4double anEnergy, G4double eMax, G4int n)
70 {
71 G4double result;
72 result = std::sqrt(anEnergy) * G4Pow::GetInstance()->powA(eMax - anEnergy, 3. * n / 2. - 4.);
73 return result;
74 }
75
76 inline G4double C(G4double anEnergy, G4double mass)
77 {
78 G4double result(0);
79 if (theTotalCount == 3)
80 result = 4. / CLHEP::pi / G4Pow::GetInstance()->powN(GetEmax(anEnergy, mass), 2);
81 if (theTotalCount == 4)
82 result = 105. / 32. / G4Pow::GetInstance()->powA(GetEmax(anEnergy, mass), 3.5);
83 // if(theTotalCount==5) result =
84 // 256./14./CLHEP::pi/G4Pow::GetInstance()->powA(GetEmax(anEnergy, mass), 5.);
85 if (theTotalCount == 5)
86 result = 256. / 14. / CLHEP::pi / G4Pow::GetInstance()->powN(GetEmax(anEnergy, mass), 5);
87 return result;
88 }
89
90 inline G4double GetEmax(G4double anEnergy, G4double mass)
91 {
92 G4double result;
93 G4double tMass = GetTarget()->GetMass();
94 G4double pMass = GetProjectileRP()->GetMass();
95 G4double availableEnergy = GetQValue() + anEnergy * tMass / (pMass + tMass);
96 result = availableEnergy * (theTotalMass - mass) / theTotalMass;
97 return result;
98 }
99
100 G4double MeanEnergyOfThisInteraction() override { return -1; }
101
102 private:
103 G4double theTotalMass;
104 G4int theTotalCount;
105 };
106 #endif
107