| Geant4 Cross Reference |
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 // particle_hp -- source file 26 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // T.Koi 02-08-06 Modified Harmonise to resolv << 30 // 02-08-06 Modified Harmonise to reslove cross section trouble at high-end. T. KOI
31 // 31 //
32 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 // V. Ivanchenko July-2023 converted back capt <<
34 // 33 //
35 #include "G4ParticleHPElementData.hh" 34 #include "G4ParticleHPElementData.hh"
36 35
37 G4ParticleHPElementData::G4ParticleHPElementDa << 36 G4ParticleHPElementData::G4ParticleHPElementData()
38 { << 37 {
39 precision = 0.02; << 38 precision = 0.02;
40 theFissionData = new G4ParticleHPVector; << 39 theFissionData = new G4ParticleHPVector;
41 theCaptureData = new G4ParticleHPVector; << 40 theCaptureData = new G4ParticleHPVector;
42 theElasticData = new G4ParticleHPVector; << 41 theElasticData = new G4ParticleHPVector;
43 theInelasticData = new G4ParticleHPVector; << 42 theInelasticData = new G4ParticleHPVector;
44 theIsotopeWiseData = nullptr; << 43 theIsotopeWiseData = 0;
45 theBuffer = nullptr; << 44 theBuffer = nullptr;
46 } << 45 }
47 << 46
48 G4ParticleHPElementData::~G4ParticleHPElementD << 47 G4ParticleHPElementData::~G4ParticleHPElementData()
49 { << 48 {
50 delete theFissionData; << 49 delete theFissionData;
51 delete theCaptureData; << 50 delete theCaptureData;
52 delete theElasticData; << 51 delete theElasticData;
53 delete theInelasticData; << 52 delete theInelasticData;
54 delete[] theIsotopeWiseData; << 53 delete [] theIsotopeWiseData;
55 } <<
56 <<
57 void G4ParticleHPElementData::Init(G4Element* <<
58 G4ParticleD <<
59 const char* <<
60 { <<
61 auto nIso = (G4int)theElement->GetNumberOfIs <<
62 auto Z = theElement->GetZasInt(); <<
63 const std::size_t dsize = nIso > 0 ? nIso : <<
64 theIsotopeWiseData = new G4ParticleHPIsoData <<
65 <<
66 for (G4int i1 = 0; i1 < nIso; ++i1) { <<
67 G4int A = theElement->GetIsotope(i1)->GetN <<
68 G4int M = theElement->GetIsotope(i1)->Getm <<
69 G4double frac = theElement->GetRelativeAbu <<
70 UpdateData(A, Z, M, i1, frac, projectile, <<
71 } 54 }
72 //theElasticData->ThinOut(precision); << 55
73 //if (projectile == G4Neutron::Neutron()) th << 56 void G4ParticleHPElementData::Init(G4Element * theElement, G4ParticleDefinition* projectile, const char* dataDirVariable )
74 << 57 {
75 //theCaptureData->ThinOut(precision); << 58 G4int count = (G4int)theElement->GetNumberOfIsotopes();
76 //theFissionData->ThinOut(precision); << 59 if(count == 0)
77 } << 60 count += theStableOnes.GetNumberOfIsotopes((G4int)theElement->GetZ());
78 << 61 theIsotopeWiseData = new G4ParticleHPIsoData[count];
79 void G4ParticleHPElementData::UpdateData(G4int << 62 // filename = ein data-set je isotope.
80 G4dou << 63 count = 0;
81 G4Par << 64 G4int nIso = (G4int)theElement->GetNumberOfIsotopes();
82 const << 65 G4int Z = (G4int)theElement->GetZ();
83 { << 66
84 // Reads in the Data, using G4ParticleHPIsoD << 67 if(nIso!=0)
85 theIsotopeWiseData[index].Init(A, Z, M, abun << 68 {
86 << 69 for (G4int i1=0; i1<nIso; ++i1)
87 theBuffer = theIsotopeWiseData[index].MakeEl << 70 {
88 Harmonise(theElasticData, theBuffer); << 71 G4int A = theElement->GetIsotope(i1)->GetN();
89 delete theBuffer; << 72 G4int M = theElement->GetIsotope(i1)->Getm();
90 << 73 G4double frac = theElement->GetRelativeAbundanceVector()[i1]/CLHEP::perCent;
91 theBuffer = theIsotopeWiseData[index].MakeIn << 74 //UpdateData(A, Z, count++, frac);
92 Harmonise(theInelasticData, theBuffer); << 75 UpdateData(A, Z, M, count++, frac, projectile, dataDirVariable);
93 delete theBuffer; << 76 }
94 << 77 }
95 theBuffer = theIsotopeWiseData[index].MakeCa << 78 else
96 Harmonise(theCaptureData, theBuffer); << 79 {
97 delete theBuffer; << 80 G4int first = theStableOnes.GetFirstIsotope(Z);
98 << 81 for(G4int i1=0;
99 theBuffer = theIsotopeWiseData[index].MakeFi << 82 i1<theStableOnes.GetNumberOfIsotopes((G4int)theElement->GetZ()); ++i1)
100 Harmonise(theFissionData, theBuffer); << 83 {
101 delete theBuffer; << 84 G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
102 } << 85 G4double frac = theStableOnes.GetAbundance(first+i1);
>> 86 UpdateData(A, Z, count++, frac, projectile, dataDirVariable);
>> 87 }
>> 88 }
>> 89 theElasticData->ThinOut(precision);
>> 90 if( projectile == G4Neutron::Neutron() ) theInelasticData->ThinOut(precision);
103 91
104 void G4ParticleHPElementData::Harmonise(G4Part << 92 theCaptureData->ThinOut(precision);
105 { << 93 theFissionData->ThinOut(precision);
106 if (theNew == nullptr) { << 94 }
107 return; << 95
108 } << 96 void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance, G4ParticleDefinition* projectile, const char* dataDirVariable )
109 G4int s_tmp = 0, n = 0, m_tmp = 0; << 97 {
110 auto theMerge = new G4ParticleHPVector(theSt << 98 //Reads in the Data, using G4ParticleHPIsoData[], and its Init
111 while (theStore->GetEnergy(s_tmp) < theNew-> << 99 //theIsotopeWiseData[index].Init(A, Z, abundance);
112 && s_tmp < theStore->GetVectorLength( << 100 theIsotopeWiseData[index].Init(A, Z, M, abundance,projectile, dataDirVariable);
113 { << 101
114 theMerge->SetData(m_tmp++, theStore->GetEn << 102 theBuffer = theIsotopeWiseData[index].MakeElasticData();
115 ++s_tmp; << 103 Harmonise(theElasticData, theBuffer);
116 } << 104 delete theBuffer;
117 G4ParticleHPVector* active = theStore; << 105
118 G4ParticleHPVector* passive = theNew; << 106 theBuffer = theIsotopeWiseData[index].MakeInelasticData();
119 G4ParticleHPVector* tmp; << 107 Harmonise(theInelasticData, theBuffer);
120 G4int a = s_tmp, p = n, t; << 108 delete theBuffer;
121 while (a < active->GetVectorLength() << 109
122 && p < passive->GetVectorLength()) / << 110 theBuffer = theIsotopeWiseData[index].MakeCaptureData();
123 { << 111 Harmonise(theCaptureData, theBuffer);
124 if (active->GetEnergy(a) <= passive->GetEn << 112 delete theBuffer;
125 theMerge->SetData(m_tmp, active->GetEner << 113
126 G4double x = theMerge->GetEnergy(m_tmp); << 114 theBuffer = theIsotopeWiseData[index].MakeFissionData();
127 G4double y = std::max(0., passive->GetXs << 115 Harmonise(theFissionData, theBuffer);
128 theMerge->SetData(m_tmp, x, theMerge->Ge << 116 delete theBuffer;
129 ++m_tmp; << 117 }
>> 118
>> 119 void G4ParticleHPElementData::Harmonise(G4ParticleHPVector *& theStore, G4ParticleHPVector * theNew)
>> 120 {
>> 121 if(theNew == 0) { return; }
>> 122 G4int s_tmp = 0, n=0, m_tmp=0;
>> 123 G4ParticleHPVector * theMerge = new G4ParticleHPVector(theStore->GetVectorLength());
>> 124 while ( theStore->GetEnergy(s_tmp)<theNew->GetEnergy(0)&&s_tmp<theStore->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi
>> 125 {
>> 126 theMerge->SetData(m_tmp++, theStore->GetEnergy(s_tmp), theStore->GetXsec(s_tmp));
>> 127 ++s_tmp;
>> 128 }
>> 129 G4ParticleHPVector *active = theStore;
>> 130 G4ParticleHPVector * passive = theNew;
>> 131 G4ParticleHPVector * tmp;
>> 132 G4int a = s_tmp, p = n, t;
>> 133 while (a<active->GetVectorLength()&&p<passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 134 {
>> 135 if(active->GetEnergy(a) <= passive->GetEnergy(p))
>> 136 {
>> 137 theMerge->SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
>> 138 G4double x = theMerge->GetEnergy(m_tmp);
>> 139 G4double y = std::max(0., passive->GetXsec(x));
>> 140 theMerge->SetData(m_tmp, x, theMerge->GetXsec(m_tmp)+y);
>> 141 ++m_tmp;
>> 142 ++a;
>> 143 }
>> 144 else
>> 145 {
>> 146 tmp = active; t=a;
>> 147 active = passive; a=p;
>> 148 passive = tmp; p=t;
>> 149 }
>> 150 }
>> 151 while (a!=active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 152 {
>> 153 theMerge->SetData(m_tmp++, active->GetEnergy(a), active->GetXsec(a));
130 ++a; 154 ++a;
131 } 155 }
132 else { << 156 while (p!=passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
133 tmp = active; << 157 {
134 t = a; << 158 G4double x = passive->GetEnergy(p);
135 active = passive; << 159 G4double y = std::max(0., active->GetXsec(x));
136 a = p; << 160 theMerge->SetData(m_tmp++, x, passive->GetXsec(p)+y);
137 passive = tmp; << 161 ++p;
138 p = t; << 162 }
139 } << 163 delete theStore;
140 } << 164 theStore = theMerge;
141 while (a != active->GetVectorLength()) // L << 165 }
142 { <<
143 theMerge->SetData(m_tmp++, active->GetEner <<
144 ++a; <<
145 } <<
146 while (p != passive->GetVectorLength()) // <<
147 { <<
148 G4double x = passive->GetEnergy(p); <<
149 G4double y = std::max(0., active->GetXsec( <<
150 theMerge->SetData(m_tmp++, x, passive->Get <<
151 ++p; <<
152 } <<
153 delete theStore; <<
154 theStore = theMerge; <<
155 } <<
156 <<
157 G4ParticleHPVector* <<
158 G4ParticleHPElementData::MakePhysicsVector(G4E <<
159 G4P <<
160 G4P <<
161 cha <<
162 { <<
163 Init(theElement, projectile, dataDirVariable <<
164 return GetData(theSet); <<
165 } <<
166 <<
167 G4ParticleHPVector* <<
168 G4ParticleHPElementData::MakePhysicsVector(G4E <<
169 G4P <<
170 G4N <<
171 cha <<
172 { <<
173 Init(theElement, projectile, dataDirVariable <<
174 return GetData(theSet); <<
175 } <<
176 166
177 G4ParticleHPVector* << 167 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
178 G4ParticleHPElementData::MakePhysicsVector(G4E << 168 G4ParticleDefinition * projectile,
179 G4P << 169 G4ParticleHPFissionData* theSet,
180 G4P << 170 char* dataDirVariable)
181 cha << 171 {
182 { << 172 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
183 Init(theElement, projectile, dataDirVariable << 173 Init ( theElement, projectile, dataDirVariable );
184 return GetData(theSet); << 174 return GetData(theSet);
185 } << 175 }
186 G4ParticleHPVector* << 176 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
187 G4ParticleHPElementData::MakePhysicsVector(G4E << 177 G4ParticleDefinition * projectile,
188 G4P << 178 G4ParticleHPCaptureData * theSet,
189 G4P << 179 char* dataDirVariable)
190 cha << 180 {
191 { << 181 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
192 Init(theElement, projectile, dataDirVariable << 182 Init ( theElement, projectile, dataDirVariable );
193 return GetData(theSet); << 183 return GetData(theSet);
194 } << 184 }
>> 185 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
>> 186 G4ParticleDefinition * projectile,
>> 187 G4ParticleHPElasticData * theSet,
>> 188 char* dataDirVariable)
>> 189 {
>> 190 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
>> 191 Init ( theElement, projectile, dataDirVariable );
>> 192 return GetData(theSet);
>> 193 }
>> 194 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
>> 195 G4ParticleDefinition * projectile,
>> 196 G4ParticleHPInelasticData * theSet,
>> 197 char* dataDirVariable)
>> 198 {
>> 199 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
>> 200 Init ( theElement, projectile, dataDirVariable );
>> 201 return GetData(theSet);
>> 202 }
195 203