| 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 // GEM de-excitation model
27 // by V. Ivanchenko (July 2019)
28 //
29 #include "G4GEMProbabilityVI.hh"
30 #include "G4NuclearLevelData.hh"
31 #include "G4LevelManager.hh"
32 #include "G4PairingCorrection.hh"
33 #include "G4NucleiProperties.hh"
34 #include "G4RandomDirection.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "Randomize.hh"
38 #include "G4Pow.hh"
39 #include "G4Exp.hh"
40
41
42 G4GEMProbabilityVI::G4GEMProbabilityVI(G4int anA, G4int aZ, const G4LevelManager* p)
43 : G4VEmissionProbability(aZ, anA), lManager(p)
44 {
45 fragA = fragZ = 0;
46 resA13 = U = delta0 = delta1 = a0 = a1 = probmax = alphaP = betaP = 0.0;
47 Umax = bCoulomb = 0.0;
48 Gamma = 1.0;
49 pcoeff = Gamma*pEvapMass*CLHEP::millibarn
50 /((CLHEP::pi*CLHEP::hbarc)*(CLHEP::pi*CLHEP::hbarc));
51 coeff = CLHEP::fermi*CLHEP::fermi/(CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
52
53 isExcited = (!lManager || 0.0 == lManager->MaxLevelEnergy()) ? false : true;
54 A13 = pG4pow->Z13(theA);
55
56 if(0 == aZ) {
57 ResetIntegrator(30, 0.25*CLHEP::MeV, 0.02);
58 } else {
59 ResetIntegrator(30, 0.5*CLHEP::MeV, 0.03);
60 }
61 }
62
63 G4double G4GEMProbabilityVI::TotalProbability(
64 const G4Fragment& fragment,
65 const G4double tmin, const G4double tmax,
66 const G4double CB, const G4double exEnergy,
67 const G4double exEvap)
68 {
69 fragA = fragment.GetA_asInt();
70 fragZ = fragment.GetZ_asInt();
71
72 bCoulomb = CB;
73 U = fragment.GetExcitationEnergy();
74 delta0 = pNuclearLevelData->GetPairingCorrection(fragZ,fragA);
75 delta1 = pNuclearLevelData->GetPairingCorrection(resZ,resA);
76 Umax = pMass - pEvapMass - pResMass - CB;
77 if(0.0 >= Umax) { return 0.0; }
78
79 resA13 = pG4pow->Z13(resA);
80 a0 = pNuclearLevelData->GetLevelDensity(fragZ,fragA,U);
81 const G4double twoMass = pMass + pMass;
82 const G4double evapMass2 = pEvapMass*pEvapMass;
83 G4double ekinmax =
84 ((pMass-pResMass)*(pMass+pResMass) + evapMass2)/twoMass - pEvapMass;
85 G4double ekinmin =
86 std::max((CB*(twoMass - CB) + evapMass2)/twoMass - pEvapMass,0.0);
87 if(ekinmax <= ekinmin) { return 0.0; }
88 pProbability = IntegrateProbability(ekinmin, ekinmax, CB);
89 pProbability += tmax - tmin + exEnergy -exEvap;
90 /*
91 G4cout << "G4GEMProbabilityVI: Z= " << theZ << " A= " << theA
92 << " resZ= " << resZ << " resA= " << resA
93 << " fragZ= " << fragZ << " fragA= " << fragA
94 << " prob= " << pProbability
95 << "\n U= " << U << " Umax= " << Umax << " d0= " << delta0
96 << " a0= " << a0 << G4endl;
97 */
98 return pProbability;
99 }
100
101 G4double G4GEMProbabilityVI::ComputeProbability(G4double ekin, G4double)
102 {
103 // abnormal case - should never happens
104 if(pMass < pEvapMass + pResMass) { return 0.0; }
105
106 const G4double m02 = pMass*pMass;
107 const G4double m12 = pEvapMass*pEvapMass;
108 const G4double mres = std::sqrt(m02 + m12 - 2.*pMass*(pEvapMass + ekin));
109
110 G4double excRes = std::max(mres - pResMass, 0.0);
111 a1 = pNuclearLevelData->GetLevelDensity(resZ,resA,excRes);
112 G4double prob = 0.5; //CrossSection(0.0, excRes);
113
114 //G4cout<<"### G4GEMProbabilityVI::ComputeProbability: Ekin(MeV)= "<<ekin
115 //<< " excRes(MeV)= " << excRes << " prob= " << prob << << G4endl;
116 return prob;
117 }
118
119 G4double G4GEMProbabilityVI::SampleEnergy(
120 const G4double tmin, const G4double tmax,
121 const G4double CB, const G4double exEnergy,
122 const G4double exEvap)
123 {
124 G4double ekin = tmax - tmin - CB -exEnergy + exEvap;
125 return ekin;
126 }
127