| 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 // neutron_hp -- source file 26 // neutron_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 // 070523 bug fix for G4FPE_DEBUG on by A. How 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 071031 bug fix T. Koi on behalf of A. Howar << 31 // 071031 bug fix T. Koi on behalf of A. Howard
32 // 081203 bug fix in Register method by T. Koi 32 // 081203 bug fix in Register method by T. Koi
33 // 33 //
34 // P. Arce, June-2014 Conversion neutron_hp to 34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 // 35 //
36 // June-2019 - E. Mendoza --> Modification to << 36 // June-2019 - E. Mendoza --> Modification to allow using an incomplete data library if the G4NEUTRONHP_SKIP_MISSING_ISOTOPES environmental flag is defined. The missing XS are set to 0.
37 // data library if the G4NEUTRONHP_SKIP_MISS <<
38 // flag is defined. The missing XS are set t <<
39 37
40 #include "G4ParticleHPChannel.hh" <<
41 38
42 #include "G4HadTmpUtil.hh" << 39 #include <stdlib.h>
43 #include "G4ParticleHPElasticFS.hh" << 40
44 #include "G4ParticleHPFinalState.hh" << 41 #include "G4ParticleHPChannel.hh"
45 #include "G4ParticleHPReactionWhiteBoard.hh" <<
46 #include "G4ParticleHPThermalBoost.hh" <<
47 #include "G4SystemOfUnits.hh" <<
48 #include "globals.hh" 42 #include "globals.hh"
>> 43 #include "G4SystemOfUnits.hh"
>> 44 #include "G4ParticleHPFinalState.hh"
>> 45 #include "G4HadTmpUtil.hh"
49 46
50 #include <cstdlib> << 47 #include "G4ParticleHPManager.hh"
>> 48 #include "G4ParticleHPReactionWhiteBoard.hh"
51 49
52 G4ParticleHPChannel::G4ParticleHPChannel(G4Par << 50 G4double G4ParticleHPChannel::GetXsec(G4double energy)
53 { << 51 {
54 fManager = G4ParticleHPManager::GetInstance( << 52 return std::max(0., theChannelData->GetXsec(energy));
55 if (fManager->GetUseWendtFissionModel()) { <<
56 wendtFissionGenerator = G4WendtFissionFrag <<
57 // Make sure both fission fragment models <<
58 fManager->SetProduceFissionFragments(false <<
59 } 53 }
60 theProjectile = (nullptr == p) ? G4Neutron:: << 54
61 theChannelData = new G4ParticleHPVector; << 55 G4double G4ParticleHPChannel::GetWeightedXsec(G4double energy, G4int isoNumber)
62 } << 56 {
63 << 57 return theIsotopeWiseData[isoNumber].GetXsec(energy);
64 G4ParticleHPChannel::~G4ParticleHPChannel() <<
65 { <<
66 delete theChannelData; <<
67 // Following statement disabled to avoid SEG <<
68 // theBuffer is also deleted as "theChannelD <<
69 delete[] theIsotopeWiseData; <<
70 if (theFinalStates != nullptr) { <<
71 for (G4int i = 0; i < niso; i++) { <<
72 delete theFinalStates[i]; <<
73 } <<
74 delete[] theFinalStates; <<
75 } 58 }
76 delete[] active; << 59
77 } << 60 G4double G4ParticleHPChannel::GetFSCrossSection(G4double energy, G4int isoNumber)
78 << 61 {
79 G4double G4ParticleHPChannel::GetXsec(G4double << 62 return theFinalStates[isoNumber]->GetXsec(energy);
80 { <<
81 return std::max(0., theChannelData->GetXsec( <<
82 } <<
83 <<
84 G4double G4ParticleHPChannel::GetWeightedXsec( <<
85 G4int isoNumber) const <<
86 { <<
87 return theIsotopeWiseData[isoNumber].GetXsec <<
88 } <<
89 <<
90 G4double G4ParticleHPChannel::GetFSCrossSectio <<
91 G4int isoNumber) const <<
92 { <<
93 return theFinalStates[isoNumber]->GetXsec(en <<
94 } <<
95 <<
96 void G4ParticleHPChannel::Init(G4Element* anEl <<
97 const G4String& dirName, const G4 <<
98 { <<
99 theFSType = aFSType; <<
100 Init(anElement, dirName); <<
101 } <<
102 <<
103 void G4ParticleHPChannel::Init(G4Element* anEl <<
104 { <<
105 theDir = dirName; <<
106 theElement = anElement; <<
107 } <<
108 <<
109 G4bool G4ParticleHPChannel::Register(G4Particl <<
110 { <<
111 ++registerCount; <<
112 G4int Z = theElement->GetZasInt(); <<
113 <<
114 niso = (G4int)theElement->GetNumberOfIsotope <<
115 const std::size_t nsize = niso > 0 ? niso : <<
116 <<
117 delete[] theIsotopeWiseData; <<
118 theIsotopeWiseData = new G4ParticleHPIsoData <<
119 delete[] active; <<
120 active = new G4bool[nsize]; <<
121 <<
122 delete[] theFinalStates; <<
123 theFinalStates = new G4ParticleHPFinalState* <<
124 delete theChannelData; <<
125 theChannelData = new G4ParticleHPVector; <<
126 for (G4int i = 0; i < niso; ++i) { <<
127 theFinalStates[i] = theFS->New(); <<
128 theFinalStates[i]->SetProjectile(theProjec <<
129 } <<
130 if (niso != 0 && registerCount == 0) { <<
131 for (G4int i1 = 0; i1 < niso; ++i1) { <<
132 G4int A = theElement->GetIsotope(i1)->Ge <<
133 G4int M = theElement->GetIsotope(i1)->Ge <<
134 //G4cout <<" Init: normal case i=" << i1 <<
135 // << " Z=" << Z << " A=" << A << G4 <<
136 G4double frac = theElement->GetRelativeA <<
137 theFinalStates[i1]->SetA_Z(A, Z, M); <<
138 UpdateData(A, Z, M, i1, frac, theProject <<
139 } <<
140 } 63 }
141 G4bool result = HasDataInAnyFinalState(); << 64
142 << 65 void G4ParticleHPChannel::
143 // To avoid issuing hash by worker threads << 66 Init(G4Element * anElement, const G4String dirName, const G4String aFSType)
144 if (result) theChannelData->Hash(); << 67 {
145 << 68 theFSType = aFSType;
146 return result; << 69 Init(anElement, dirName);
147 } <<
148 <<
149 void G4ParticleHPChannel::UpdateData(G4int A, <<
150 G4double <<
151 G4Particl <<
152 { <<
153 // Initialze the G4FissionFragment generator <<
154 if (wendtFissionGenerator != nullptr) { <<
155 wendtFissionGenerator->InitializeANucleus( <<
156 } <<
157 <<
158 theFinalStates[index]->Init(A, Z, M, theDir, <<
159 if (!theFinalStates[index]->HasAnyData()) re <<
160 // nothing there for exactly this isotope. <<
161 <<
162 // the above has put the X-sec into the FS <<
163 theBuffer = nullptr; <<
164 if (theFinalStates[index]->HasXsec()) { <<
165 theBuffer = theFinalStates[index]->GetXsec <<
166 theBuffer->Times(abundance / 100.); <<
167 theIsotopeWiseData[index].FillChannelData( <<
168 } 70 }
169 else // get data from CrossSection director << 71
>> 72 void G4ParticleHPChannel::Init(G4Element * anElement, const G4String dirName)
170 { 73 {
171 const G4String& tString = "/CrossSection"; << 74 theDir = dirName;
172 active[index] = theIsotopeWiseData[index]. << 75 theElement = anElement;
173 <<
174 if (active[index]) theBuffer = theIsotopeW <<
175 } 76 }
176 if (theBuffer != nullptr) Harmonise(theChann << 77
177 } << 78 G4bool G4ParticleHPChannel::Register(G4ParticleHPFinalState *theFS)
178 <<
179 void G4ParticleHPChannel::Harmonise(G4Particle <<
180 G4Particle <<
181 { <<
182 G4int s_tmp = 0, n = 0, m_tmp = 0; <<
183 auto theMerge = new G4ParticleHPVector; <<
184 G4ParticleHPVector* anActive = theStore; <<
185 G4ParticleHPVector* aPassive = theNew; <<
186 G4ParticleHPVector* tmp; <<
187 G4int a = s_tmp, p = n, t; <<
188 while (a < anActive->GetVectorLength() && p <<
189 // Loop checking, 11.05.2015, T. Koi <<
190 { 79 {
191 if (anActive->GetEnergy(a) <= aPassive->Ge << 80 registerCount++;
192 G4double xa = anActive->GetEnergy(a); << 81 G4int Z = G4lrint(theElement->GetZ());
193 theMerge->SetData(m_tmp, xa, anActive->G << 82
194 m_tmp++; << 83 Z = Z-registerCount;
195 a++; << 84 if ( registerCount > 5 ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material"); // for Elastic, Capture, Fission case
196 G4double xp = aPassive->GetEnergy(p); << 85 if ( Z < 1 ) return false;
197 if (std::abs(std::abs(xp - xa) / xa) < 0 << 86 /*
198 ++p; << 87 if(registerCount<5)
199 } << 88 {
>> 89 Z = Z-registerCount;
>> 90 }
>> 91 */
>> 92 //if(Z=theElement->GetZ()-5) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material");
>> 93 // Bug fix by TK on behalf of AH
>> 94 //if ( Z <=theElement->GetZ()-5 ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material");
>> 95 G4int count = 0;
>> 96 if(registerCount==0) count = theElement->GetNumberOfIsotopes();
>> 97 if(count == 0||registerCount!=0) count +=
>> 98 theStableOnes.GetNumberOfIsotopes(Z);
>> 99 niso = count;
>> 100 delete [] theIsotopeWiseData;
>> 101 theIsotopeWiseData = new G4ParticleHPIsoData [niso];
>> 102 delete [] active;
>> 103 active = new G4bool[niso];
>> 104
>> 105 delete [] theFinalStates;
>> 106 theFinalStates = new G4ParticleHPFinalState * [niso];
>> 107 delete theChannelData;
>> 108 theChannelData = new G4ParticleHPVector;
>> 109 for(G4int i=0; i<niso; i++)
>> 110 {
>> 111 theFinalStates[i] = theFS->New();
>> 112 theFinalStates[i]->SetProjectile(theProjectile);
200 } 113 }
201 else { << 114 count = 0;
202 tmp = anActive; << 115 G4int nIsos = niso;
203 t = a; << 116 if(theElement->GetNumberOfIsotopes()!=0&®isterCount==0)
204 anActive = aPassive; << 117 {
205 a = p; << 118 for (G4int i1=0; i1<nIsos; i1++)
206 aPassive = tmp; << 119 {
207 p = t; << 120 // G4cout <<" Init: normal case"<<G4endl;
>> 121 G4int A = theElement->GetIsotope(i1)->GetN();
>> 122 G4int M = theElement->GetIsotope(i1)->Getm();
>> 123 G4double frac = theElement->GetRelativeAbundanceVector()[i1]/perCent;
>> 124 //theFinalStates[i1]->SetA_Z(A, Z);
>> 125 //UpdateData(A, Z, count++, frac);
>> 126 theFinalStates[i1]->SetA_Z(A, Z, M);
>> 127 UpdateData(A, Z, M, count++, frac, theProjectile);
>> 128 }
>> 129 } else {
>> 130 //G4cout <<" Init: mean case: "
>> 131 // <<theStableOnes.GetNumberOfIsotopes(Z)<<" "
>> 132 // <<Z<<" "<<theElement
>> 133 // << G4endl;
>> 134 G4int first = theStableOnes.GetFirstIsotope(Z);
>> 135 for(G4int i1=0;
>> 136 i1<theStableOnes.GetNumberOfIsotopes(Z);
>> 137 i1++)
>> 138 {
>> 139 G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
>> 140 G4double frac = theStableOnes.GetAbundance(first+i1);
>> 141 theFinalStates[i1]->SetA_Z(A, Z);
>> 142 UpdateData(A, Z, count++, frac, theProjectile);
>> 143 }
208 } 144 }
>> 145 G4bool result = HasDataInAnyFinalState();
>> 146
>> 147 //To avoid issuing hash by worker threads
>> 148 if ( result ) theChannelData->Hash();
>> 149
>> 150 return result;
209 } 151 }
210 while (a != anActive->GetVectorLength()) // << 152
>> 153 //void G4ParticleHPChannel::UpdateData(G4int A, G4int Z, G4int index, G4double abundance)
>> 154 void G4ParticleHPChannel::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance, G4ParticleDefinition* projectile)
211 { 155 {
212 theMerge->SetData(m_tmp++, anActive->GetEn << 156 // Initialze the G4FissionFragment generator for this isomer if needed
213 ++a; << 157 if(wendtFissionGenerator)
>> 158 {
>> 159 wendtFissionGenerator->InitializeANucleus(A, Z, M, theDir);
>> 160 }
>> 161
>> 162 theFinalStates[index]->Init(A, Z, M, theDir, theFSType, projectile);
>> 163 if(!theFinalStates[index]->HasAnyData()) return; // nothing there for exactly this isotope.
>> 164
>> 165 // the above has put the X-sec into the FS
>> 166 theBuffer = 0;
>> 167 if(theFinalStates[index]->HasXsec())
>> 168 {
>> 169 theBuffer = theFinalStates[index]->GetXsec();
>> 170 theBuffer->Times(abundance/100.);
>> 171 theIsotopeWiseData[index].FillChannelData(theBuffer);
>> 172 }
>> 173 else // get data from CrossSection directory
>> 174 {
>> 175 G4String tString = "/CrossSection";
>> 176 //active[index] = theIsotopeWiseData[index].Init(A, Z, abundance, theDir, tString);
>> 177 active[index] = theIsotopeWiseData[index].Init(A, Z, M, abundance, theDir, tString);
>> 178 if(active[index]) theBuffer = theIsotopeWiseData[index].MakeChannelData();
>> 179 }
>> 180 if(theBuffer != 0) Harmonise(theChannelData, theBuffer);
214 } 181 }
215 while (p != aPassive->GetVectorLength()) // << 182
>> 183 void G4ParticleHPChannel::Harmonise(G4ParticleHPVector *& theStore, G4ParticleHPVector * theNew)
216 { 184 {
217 if (std::abs(theMerge->GetEnergy(std::max( << 185 G4int s_tmp = 0, n=0, m_tmp=0;
218 aPassive->GetEnergy(p)) / aPassive->GetEn << 186 G4ParticleHPVector * theMerge = new G4ParticleHPVector;
219 theMerge->SetData(m_tmp++, aPassive->Get << 187 G4ParticleHPVector * anActive = theStore;
220 ++p; << 188 G4ParticleHPVector * aPassive = theNew;
>> 189 G4ParticleHPVector * tmp;
>> 190 G4int a = s_tmp, p = n, t;
>> 191 while (a<anActive->GetVectorLength()&&p<aPassive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 192 {
>> 193 if(anActive->GetEnergy(a) <= aPassive->GetEnergy(p))
>> 194 {
>> 195 G4double xa = anActive->GetEnergy(a);
>> 196 theMerge->SetData(m_tmp, xa, anActive->GetXsec(a)+std::max(0., aPassive->GetXsec(xa)) );
>> 197 m_tmp++;
>> 198 a++;
>> 199 G4double xp = aPassive->GetEnergy(p);
>> 200 if( std::abs(std::abs(xp-xa)/xa)<0.001 )
>> 201 {
>> 202 p++;
>> 203 }
>> 204 } else {
>> 205 tmp = anActive; t=a;
>> 206 anActive = aPassive; a=p;
>> 207 aPassive = tmp; p=t;
>> 208 }
>> 209 }
>> 210 while (a!=anActive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 211 {
>> 212 theMerge->SetData(m_tmp++, anActive->GetEnergy(a), anActive->GetXsec(a));
>> 213 a++;
>> 214 }
>> 215 while (p!=aPassive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 216 {
>> 217 if(std::abs(theMerge->GetEnergy(std::max(0,m_tmp-1))-aPassive->GetEnergy(p))/aPassive->GetEnergy(p)>0.001)
>> 218 theMerge->SetData(m_tmp++, aPassive->GetEnergy(p), aPassive->GetXsec(p));
>> 219 p++;
>> 220 }
>> 221 delete theStore;
>> 222 theStore = theMerge;
221 } 223 }
222 delete theStore; <<
223 theStore = theMerge; <<
224 } <<
225 224
226 G4WendtFissionFragmentGenerator* G4ParticleHPC << 225 #include "G4ParticleHPThermalBoost.hh"
227 if ( wendtFissionGenerator ) return wendtFis <<
228 else return nullptr; <<
229 } <<
230 226
231 G4HadFinalState* << 227 G4HadFinalState * G4ParticleHPChannel::
232 G4ParticleHPChannel::ApplyYourself(const G4Had << 228 ApplyYourself(const G4HadProjectile & theTrack, G4int anIsotope)
233 G4int anIsotope, G4bool isElastic) << 229 {
234 { << 230 // G4cout << "G4ParticleHPChannel::ApplyYourself+"<<niso<<G4endl;
235 //G4cout << "G4ParticleHPChannel::ApplyYours << 231 if ( anIsotope != -1 && anIsotope != -2 )
236 // << " ni=" << anIsotope << " isElastic=" << 232 {
237 if (anIsotope != -1 && anIsotope != -2) { << 233 //Inelastic Case
238 // Inelastic Case << 234 //G4cout << "G4ParticleHPChannel Inelastic Case"
239 //G4cout << "G4ParticleHPChannel Inelastic << 235 //<< " Z= " << this->GetZ(it) << " A = " << this->GetN(it) << G4endl;
240 //<< " Z= " << GetZ(anIsotope) << " A = " << 236 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA( (G4int)this->GetN(anIsotope) );
241 fManager->GetReactionWhiteBoard()->SetTarg << 237 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ( (G4int)this->GetZ(anIsotope) );
242 fManager->GetReactionWhiteBoard()->SetTarg << 238 return theFinalStates[anIsotope]->ApplyYourself(theTrack);
243 return theFinalStates[anIsotope]->ApplyYou <<
244 } <<
245 G4double sum = 0; <<
246 G4int it = 0; <<
247 auto xsec = new G4double[niso]; <<
248 G4ParticleHPThermalBoost aThermalE; <<
249 for (G4int i = 0; i < niso; i++) { <<
250 if (theFinalStates[i]->HasAnyData()) { <<
251 /* <<
252 G4cout << "FS: " << i << theTrack.GetDef <<
253 << " Z=" << theFinalStates[i]->GetZ() <<
254 << " A=" << theFinalStates[i]->GetN() <<
255 << G4endl; <<
256 */ <<
257 xsec[i] = theIsotopeWiseData[i].GetXsec( <<
258 aThermalE.GetThermalEnergy(theTrack, t <<
259 theFinalSta <<
260 theTrack.Ge <<
261 sum += xsec[i]; <<
262 } <<
263 else { <<
264 xsec[i] = 0; <<
265 } <<
266 } <<
267 if (sum == 0) { <<
268 it = G4lrint(niso * G4UniformRand()); <<
269 } <<
270 else { <<
271 G4double random = G4UniformRand(); <<
272 G4double running = 0; <<
273 for (G4int ix = 0; ix < niso; ix++) { <<
274 running += xsec[ix]; <<
275 if (sum == 0 || random <= running / sum) <<
276 it = ix; <<
277 break; <<
278 } <<
279 } 239 }
280 if (it == niso) it--; << 240 G4double sum=0;
281 } << 241 G4int it=0;
282 delete[] xsec; << 242 G4double * xsec = new G4double[niso];
283 G4HadFinalState* theFinalState = nullptr; << 243 G4ParticleHPThermalBoost aThermalE;
284 const auto A = (G4int)this->GetN(it); << 244 for (G4int i=0; i<niso; i++)
285 const auto Z = (G4int)this->GetZ(it); << 245 {
286 const auto M = (G4int)this->GetM(it); << 246 if(theFinalStates[i]->HasAnyData())
287 << 247 {
288 //-2:Marker for Fission << 248 xsec[i] = theIsotopeWiseData[i].GetXsec(aThermalE.GetThermalEnergy(theTrack,
289 if ((wendtFissionGenerator != nullptr) && an << 249 theFinalStates[i]->GetN(),
290 theFinalState = wendtFissionGenerator->App << 250 theFinalStates[i]->GetZ(),
291 } << 251 theTrack.GetMaterial()->GetTemperature()));
292 << 252 sum += xsec[i];
293 // Use the standard procedure if the G4Fissi <<
294 if (theFinalState == nullptr) { <<
295 G4int icounter = 0; <<
296 G4int icounter_max = 1024; <<
297 while (theFinalState == nullptr) // Loop <<
298 { <<
299 icounter++; <<
300 if (icounter > icounter_max) { <<
301 G4cout << "Loop-counter exceeded the t <<
302 << __LINE__ << "th line of " << <<
303 break; <<
304 } 253 }
305 if (isElastic) { << 254 else
306 // Register 0 K cross-section for DBRC << 255 {
307 G4ParticleHPVector* xsforFS = theIsoto << 256 xsec[i]=0;
308 // Only G4ParticleHPElasticFS has the <<
309 static_cast<G4ParticleHPElasticFS*>(th <<
310 } 257 }
311 theFinalState = theFinalStates[it]->Appl << 258 }
>> 259 if(sum == 0)
>> 260 {
>> 261 // G4cout << "G4ParticleHPChannel::ApplyYourself theFinalState->Initialize+"<<G4endl;
>> 262 // G4cout << "G4ParticleHPChannel::ApplyYourself theFinalState->Initialize-"<<G4endl;
>> 263 it = static_cast<G4int>(niso*G4UniformRand());
312 } 264 }
313 } << 265 else
>> 266 {
>> 267 // G4cout << "Are we still here? "<<sum<<G4endl;
>> 268 // G4cout << "TESTHP 23 NISO="<<niso<<G4endl;
>> 269 G4double random = G4UniformRand();
>> 270 G4double running=0;
>> 271 // G4cout << "G4ParticleHPChannel::ApplyYourself Done the sum"<<niso<<G4endl;
>> 272 // G4cout << "TESTHP 24 NISO="<<niso<<G4endl;
>> 273 for (G4int ix=0; ix<niso; ix++)
>> 274 {
>> 275 running += xsec[ix];
>> 276 //if(random<=running/sum)
>> 277 if( sum == 0 || random <= running/sum )
>> 278 {
>> 279 it = ix;
>> 280 break;
>> 281 }
>> 282 }
>> 283 if(it==niso) it--;
>> 284 }
>> 285 delete [] xsec;
>> 286 G4HadFinalState * theFinalState=0;
>> 287 const G4int A = (G4int)this->GetN(it);
>> 288 const G4int Z = (G4int)this->GetZ(it);
>> 289 const G4int M = (G4int)this->GetM(it);
314 290
315 // G4cout <<"THE IMPORTANT RETURN"<<G4endl; << 291 //-2:Marker for Fission
316 // G4cout << "TK G4ParticleHPChannel Elastic << 292 if(wendtFissionGenerator&&anIsotope==-2)
317 //<< " Z= " << this->GetZ(it) << " A = " << << 293 {
318 fManager->GetReactionWhiteBoard()->SetTargA( << 294 theFinalState = wendtFissionGenerator->ApplyYourself(theTrack, Z, A);
319 fManager->GetReactionWhiteBoard()->SetTargZ( << 295 }
320 fManager->GetReactionWhiteBoard()->SetTargM( << 296
>> 297 // Use the standard procedure if the G4FissionFragmentGenerator model fails
>> 298 if (!theFinalState)
>> 299 {
>> 300
>> 301 G4int icounter=0;
>> 302 G4int icounter_max=1024;
>> 303 while(theFinalState==0) // Loop checking, 11.05.2015, T. Koi
>> 304 {
>> 305 icounter++;
>> 306 if ( icounter > icounter_max ) {
>> 307 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
>> 308 break;
>> 309 }
>> 310 // G4cout << "TESTHP 24 it="<<it<<G4endl;
>> 311 theFinalState = theFinalStates[it]->ApplyYourself(theTrack);
>> 312 }
>> 313 }
>> 314
>> 315 //G4cout <<"THE IMPORTANT RETURN"<<G4endl;
>> 316 //G4cout << "TK G4ParticleHPChannel Elastic, Capture and Fission Cases "
>> 317 //<< " Z= " << this->GetZ(it) << " A = " << this->GetN(it) << G4endl;
>> 318 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA( A );
>> 319 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ( Z );
>> 320 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargM( M );
>> 321
>> 322 return theFinalState;
>> 323 }
321 324
322 return theFinalState; <<
323 } <<
324 325
325 void G4ParticleHPChannel::DumpInfo() const << 326 void G4ParticleHPChannel::DumpInfo(){
326 { << 327
327 G4cout << " Element: " << theElement->GetNam << 328 G4cout<<" Element: "<<theElement->GetName()<<G4endl;
328 G4cout << " Directory name: " << theDir << G << 329 G4cout<<" Directory name: "<<theDir<<G4endl;
329 G4cout << " FS name: " << theFSType << G4end << 330 G4cout<<" FS name: "<<theFSType<<G4endl;
330 G4cout << " Number of Isotopes: " << niso << << 331 G4cout<<" Number of Isotopes: "<<niso<<G4endl;
331 G4cout << " Have cross sections: " << G4endl << 332 G4cout<<" Have cross sections: "<<G4endl;
332 for (int i = 0; i < niso; i++) { << 333 for(int i=0;i<niso;i++){
333 G4cout << theFinalStates[i]->HasXsec() << << 334 G4cout<<theFinalStates[i]->HasXsec()<<" ";
334 } << 335 }
335 G4cout << G4endl; << 336 G4cout<<G4endl;
336 if (theChannelData != nullptr) { << 337 if(theChannelData){
337 G4cout << " Cross Section (total for this << 338 G4cout<<" Cross Section (total for this channel):"<<G4endl;
338 int np = theChannelData->GetVectorLength() << 339 int np=theChannelData->GetVectorLength();
339 G4cout << np << G4endl; << 340 G4cout<<np<<G4endl;
340 for (int i = 0; i < np; i++) { << 341 for(int i=0;i<np;i++){
341 G4cout << theChannelData->GetEnergy(i) / << 342 G4cout<<theChannelData->GetEnergy(i)/eV<<" "<<theChannelData->GetXsec(i)<<G4endl;
342 } 343 }
343 } 344 }
>> 345
344 } 346 }
>> 347
>> 348
>> 349
345 350