| 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 = NULL;
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 = theElement->GetNumberOfIsotopes();
76 //theFissionData->ThinOut(precision); << 59 if(count == 0) count +=
77 } << 60 theStableOnes.GetNumberOfIsotopes(static_cast<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 = theElement->GetNumberOfIsotopes();
82 const << 65 G4int Z = static_cast<G4int> (theElement->GetZ());
83 { << 66 // G4int i1;
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 // G4cout <<" Init: normal case"<<G4endl;
89 delete theBuffer; << 72 G4int A = theElement->GetIsotope(i1)->GetN();
90 << 73 G4int M = theElement->GetIsotope(i1)->Getm();
91 theBuffer = theIsotopeWiseData[index].MakeIn << 74 G4double frac = theElement->GetRelativeAbundanceVector()[i1]/CLHEP::perCent;
92 Harmonise(theInelasticData, theBuffer); << 75 //UpdateData(A, Z, count++, frac);
93 delete theBuffer; << 76 UpdateData(A, Z, M, count++, frac, projectile, dataDirVariable);
94 << 77 }
95 theBuffer = theIsotopeWiseData[index].MakeCa << 78 }else{
96 Harmonise(theCaptureData, theBuffer); << 79 // G4cout <<" Init: theStableOnes case: Z="<<Z<<G4endl;
97 delete theBuffer; << 80 G4int first = theStableOnes.GetFirstIsotope(Z);
98 << 81 // G4cout <<"first="<<first<<" "<<theStableOnes.GetNumberOfIsotopes(theElement->GetZ())<<G4endl;
99 theBuffer = theIsotopeWiseData[index].MakeFi << 82 for(G4int i1=0;
100 Harmonise(theFissionData, theBuffer); << 83 i1<theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()) );
101 delete theBuffer; << 84 i1++)
102 } << 85 {
103 << 86 // G4cout <<" Init: theStableOnes in the loop"<<G4endl;
104 void G4ParticleHPElementData::Harmonise(G4Part << 87 G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
105 { << 88 G4double frac = theStableOnes.GetAbundance(first+i1);
106 if (theNew == nullptr) { << 89 // G4cout <<" Init: theStableOnes in the loop: "<<A<<G4endl;
107 return; << 90 UpdateData(A, Z, count++, frac, projectile, dataDirVariable);
108 } << 91 }
109 G4int s_tmp = 0, n = 0, m_tmp = 0; <<
110 auto theMerge = new G4ParticleHPVector(theSt <<
111 while (theStore->GetEnergy(s_tmp) < theNew-> <<
112 && s_tmp < theStore->GetVectorLength( <<
113 { <<
114 theMerge->SetData(m_tmp++, theStore->GetEn <<
115 ++s_tmp; <<
116 } <<
117 G4ParticleHPVector* active = theStore; <<
118 G4ParticleHPVector* passive = theNew; <<
119 G4ParticleHPVector* tmp; <<
120 G4int a = s_tmp, p = n, t; <<
121 while (a < active->GetVectorLength() <<
122 && p < passive->GetVectorLength()) / <<
123 { <<
124 if (active->GetEnergy(a) <= passive->GetEn <<
125 theMerge->SetData(m_tmp, active->GetEner <<
126 G4double x = theMerge->GetEnergy(m_tmp); <<
127 G4double y = std::max(0., passive->GetXs <<
128 theMerge->SetData(m_tmp, x, theMerge->Ge <<
129 ++m_tmp; <<
130 ++a; <<
131 } <<
132 else { <<
133 tmp = active; <<
134 t = a; <<
135 active = passive; <<
136 a = p; <<
137 passive = tmp; <<
138 p = t; <<
139 } 92 }
>> 93 theElasticData->ThinOut(precision);
>> 94 if( projectile == G4Neutron::Neutron() ) theInelasticData->ThinOut(precision);
>> 95
>> 96 theCaptureData->ThinOut(precision);
>> 97 theFissionData->ThinOut(precision);
140 } 98 }
141 while (a != active->GetVectorLength()) // L << 99
>> 100 //void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int index, G4double abundance)
>> 101 void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance, G4ParticleDefinition* projectile, const char* dataDirVariable )
142 { 102 {
143 theMerge->SetData(m_tmp++, active->GetEner << 103 //Reads in the Data, using G4ParticleHPIsoData[], and its Init
144 ++a; << 104 // G4cout << "entered: ElementWiseData::UpdateData"<<G4endl;
>> 105 //theIsotopeWiseData[index].Init(A, Z, abundance);
>> 106 theIsotopeWiseData[index].Init(A, Z, M, abundance,projectile, dataDirVariable);
>> 107 // G4cout << "ElementWiseData::UpdateData Init finished"<<G4endl;
>> 108
>> 109 theBuffer = theIsotopeWiseData[index].MakeElasticData();
>> 110 // G4cout << "ElementWiseData::UpdateData MakeElasticData finished: "
>> 111 // <<theBuffer->GetVectorLength()<<G4endl;
>> 112 Harmonise(theElasticData, theBuffer);
>> 113 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
>> 114 // <<theElasticData->GetVectorLength()<<G4endl;
>> 115 delete theBuffer;
>> 116
>> 117 theBuffer = theIsotopeWiseData[index].MakeInelasticData();
>> 118 // G4cout << "ElementWiseData::UpdateData MakeInelasticData finished: "
>> 119 // <<theBuffer->GetVectorLength()<<G4endl;
>> 120 Harmonise(theInelasticData, theBuffer);
>> 121 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
>> 122 // <<theInelasticData->GetVectorLength()<<G4endl;
>> 123 delete theBuffer;
>> 124
>> 125 theBuffer = theIsotopeWiseData[index].MakeCaptureData();
>> 126 // G4cout << "ElementWiseData::UpdateData MakeCaptureData finished: "
>> 127 // <<theBuffer->GetVectorLength()<<G4endl;
>> 128 Harmonise(theCaptureData, theBuffer);
>> 129 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
>> 130 // <<theCaptureData->GetVectorLength()<<G4endl;
>> 131 delete theBuffer;
>> 132
>> 133 theBuffer = theIsotopeWiseData[index].MakeFissionData();
>> 134 // G4cout << "ElementWiseData::UpdateData MakeFissionData finished: "
>> 135 // <<theBuffer->GetVectorLength()<<G4endl;
>> 136 Harmonise(theFissionData, theBuffer);
>> 137 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
>> 138 // <<theFissionData->GetVectorLength()<<G4endl;
>> 139 delete theBuffer;
>> 140
>> 141 // G4cout << "ElementWiseData::UpdateData finished"<endl;
145 } 142 }
146 while (p != passive->GetVectorLength()) // << 143
>> 144 void G4ParticleHPElementData::Harmonise(G4ParticleHPVector *& theStore, G4ParticleHPVector * theNew)
147 { 145 {
148 G4double x = passive->GetEnergy(p); << 146 if(theNew == 0) { return; }
149 G4double y = std::max(0., active->GetXsec( << 147 G4int s_tmp = 0, n=0, m_tmp=0;
150 theMerge->SetData(m_tmp++, x, passive->Get << 148 G4ParticleHPVector * theMerge = new G4ParticleHPVector(theStore->GetVectorLength());
151 ++p; << 149 // G4cout << "Harmonise 1: "<<theStore->GetEnergy(s_tmp)<<" "<<theNew->GetEnergy(0)<<G4endl;
>> 150 while ( theStore->GetEnergy(s_tmp)<theNew->GetEnergy(0)&&s_tmp<theStore->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi
>> 151 {
>> 152 theMerge->SetData(m_tmp++, theStore->GetEnergy(s_tmp), theStore->GetXsec(s_tmp));
>> 153 s_tmp++;
>> 154 }
>> 155 G4ParticleHPVector *active = theStore;
>> 156 G4ParticleHPVector * passive = theNew;
>> 157 G4ParticleHPVector * tmp;
>> 158 G4int a = s_tmp, p = n, t;
>> 159 // G4cout << "Harmonise 2: "<<active->GetVectorLength()<<" "<<passive->GetVectorLength()<<G4endl;
>> 160 while (a<active->GetVectorLength()&&p<passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 161 {
>> 162 if(active->GetEnergy(a) <= passive->GetEnergy(p))
>> 163 {
>> 164 theMerge->SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
>> 165 G4double x = theMerge->GetEnergy(m_tmp);
>> 166 G4double y = std::max(0., passive->GetXsec(x));
>> 167 theMerge->SetData(m_tmp, x, theMerge->GetXsec(m_tmp)+y);
>> 168 m_tmp++;
>> 169 a++;
>> 170 } else {
>> 171 // G4cout << "swapping in Harmonise"<<G4endl;
>> 172 tmp = active; t=a;
>> 173 active = passive; a=p;
>> 174 passive = tmp; p=t;
>> 175 }
>> 176 }
>> 177 // G4cout << "Harmonise 3: "<< a <<" "<<active->GetVectorLength()<<" "<<m<<G4endl;
>> 178 while (a!=active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 179 {
>> 180 theMerge->SetData(m_tmp++, active->GetEnergy(a), active->GetXsec(a));
>> 181 a++;
>> 182 }
>> 183 // G4cout << "Harmonise 4: "<< p <<" "<<passive->GetVectorLength()<<" "<<m<<G4endl;
>> 184 while (p!=passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 185 {
>> 186 // Modified by T. KOI
>> 187 //theMerge->SetData(m++, passive->GetEnergy(p), passive->GetXsec(p));
>> 188 G4double x = passive->GetEnergy(p);
>> 189 G4double y = std::max(0., active->GetXsec(x));
>> 190 theMerge->SetData(m_tmp++, x, passive->GetXsec(p)+y);
>> 191 p++;
>> 192 }
>> 193 // G4cout <<"Harmonise 5: "<< theMerge->GetVectorLength() << " " << m << G4endl;
>> 194 delete theStore;
>> 195 theStore = theMerge;
>> 196 // G4cout <<"Harmonise 6: "<< theStore->GetVectorLength() << " " << m << G4endl;
152 } 197 }
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 198
167 G4ParticleHPVector* << 199 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
168 G4ParticleHPElementData::MakePhysicsVector(G4E << 200 G4ParticleDefinition * projectile,
169 G4P << 201 G4ParticleHPFissionData* theSet,
170 G4N << 202 char* dataDirVariable)
171 cha << 203 {
172 { << 204 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
173 Init(theElement, projectile, dataDirVariable << 205 Init ( theElement, projectile, dataDirVariable );
174 return GetData(theSet); << 206 return GetData(theSet);
175 } << 207 }
176 << 208 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
177 G4ParticleHPVector* << 209 G4ParticleDefinition * projectile,
178 G4ParticleHPElementData::MakePhysicsVector(G4E << 210 G4ParticleHPCaptureData * theSet,
179 G4P << 211 char* dataDirVariable)
180 G4P << 212 {
181 cha << 213 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
182 { << 214 Init ( theElement, projectile, dataDirVariable );
183 Init(theElement, projectile, dataDirVariable << 215 return GetData(theSet);
184 return GetData(theSet); << 216 }
185 } << 217 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
186 G4ParticleHPVector* << 218 G4ParticleDefinition * projectile,
187 G4ParticleHPElementData::MakePhysicsVector(G4E << 219 G4ParticleHPElasticData * theSet,
188 G4P << 220 char* dataDirVariable)
189 G4P << 221 {
190 cha << 222 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
191 { << 223 Init ( theElement, projectile, dataDirVariable );
192 Init(theElement, projectile, dataDirVariable << 224 return GetData(theSet);
193 return GetData(theSet); << 225 }
194 } << 226 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
>> 227 G4ParticleDefinition * projectile,
>> 228 G4ParticleHPInelasticData * theSet,
>> 229 char* dataDirVariable)
>> 230 {
>> 231 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
>> 232 Init ( theElement, projectile, dataDirVariable );
>> 233 return GetData(theSet);
>> 234 }
195 235