Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4NeutronHPPhotonDist.cc
Go to the documentation of this file.
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// neutron_hp -- source file
27// J.P. Wellisch, Nov-1996
28// A prototype of the low energy neutron transport model.
29//
30// 070523 Try to limit sum of secondary photon energy while keeping distribution shape
31// in the of nDiscrete = 1 an nPartial = 1. Most case are satisfied.
32// T. Koi
33// 070606 Add Partial case by T. Koi
34// 070618 fix memory leaking by T. Koi
35// 080801 fix memory leaking by T. Koi
36// 080801 Correcting data disorder which happened when both InitPartial
37// and InitAnglurar methods was called in a same instance by T. Koi
38// 090514 Fix bug in IC electron emission case
39// Contribution from Chao Zhang ([email protected]) and Dongming Mei([email protected])
40// But it looks like never cause real effect in G4NDL3.13 (at least Natural elements) TK
41// 101111 Change warning message for "repFlag == 2 && isoFlag != 1" case
42//
43// there is a lot of unused (and undebugged) code in this file. Kept for the moment just in case. @@
44
45#include <numeric>
46
47#include "G4NeutronHPPhotonDist.hh"
48#include "G4PhysicalConstants.hh"
49#include "G4SystemOfUnits.hh"
50#include "G4NeutronHPLegendreStore.hh"
51#include "G4Electron.hh"
52#include "G4Poisson.hh"
53
54G4bool G4NeutronHPPhotonDist::InitMean(std::ifstream & aDataFile)
55{
56
57 G4bool result = true;
58 if(aDataFile >> repFlag)
59 {
60
61 aDataFile >> targetMass;
62 if(repFlag==1)
63 {
64 // multiplicities
65 aDataFile >> nDiscrete;
66 disType = new G4int[nDiscrete];
67 energy = new G4double[nDiscrete];
68 actualMult = new G4int[nDiscrete];
69 theYield = new G4NeutronHPVector[nDiscrete];
70 for (G4int i=0; i<nDiscrete; i++)
71 {
72 aDataFile >> disType[i]>>energy[i];
73 energy[i]*=eV;
74 theYield[i].Init(aDataFile, eV);
75 }
76 }
77 else if(repFlag == 2)
78 {
79 aDataFile >> theInternalConversionFlag;
80 aDataFile >> theBaseEnergy;
81 theBaseEnergy*=eV;
82 aDataFile >> theInternalConversionFlag;
83 // theInternalConversionFlag == 1 No IC, theInternalConversionFlag == 2 with IC
84 aDataFile >> nGammaEnergies;
85 theLevelEnergies = new G4double[nGammaEnergies];
86 theTransitionProbabilities = new G4double[nGammaEnergies];
87 if(theInternalConversionFlag == 2) thePhotonTransitionFraction = new G4double[nGammaEnergies];
88 for(G4int ii=0; ii<nGammaEnergies; ii++)
89 {
90 if(theInternalConversionFlag == 1)
91 {
92 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii];
93 theLevelEnergies[ii]*=eV;
94 }
95 else if(theInternalConversionFlag == 2)
96 {
97 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii];
98 theLevelEnergies[ii]*=eV;
99 }
100 else
101 {
102 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Unknown conversion flag");
103 }
104 }
105 // Note, that this is equivalent to using the 'Gamma' classes.
106 // throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Transition probability array not sampled for the moment.");
107 }
108 else
109 {
110 G4cout << "Data representation in G4NeutronHPPhotonDist: "<<repFlag<<G4endl;
111 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: This data representation is not implemented.");
112 }
113 }
114 else
115 {
116 result = false;
117 }
118 return result;
119}
120
121void G4NeutronHPPhotonDist::InitAngular(std::ifstream & aDataFile)
122{
123
124 G4int i, ii;
125 //angular distributions
126 aDataFile >> isoFlag;
127 if (isoFlag != 1)
128 {
129if ( repFlag == 2 ) G4cout << "G4NeutronHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 Hyper News. Thanks." << G4endl;
130 aDataFile >> tabulationType >> nDiscrete2 >> nIso;
131//080731
132 if ( theGammas != NULL && nDiscrete2 != nDiscrete ) G4cout << "080731c G4NeutronHPPhotonDist nDiscrete2 != nDiscrete, It looks like something wrong in your NDL files. Please update the latest. If you still have this messages after the update, then please report to Geant4 Hyper News." << G4endl;
133
134 // The order of cross section (InitPartials) and distribution (InitAngular here) data are different, we have to re-coordinate consistent data order.
135 std::vector < G4double > vct_gammas_par;
136 std::vector < G4double > vct_shells_par;
137 std::vector < G4int > vct_primary_par;
138 std::vector < G4int > vct_distype_par;
139 std::vector < G4NeutronHPVector* > vct_pXS_par;
140 if ( theGammas != NULL )
141 {
142 //copy the cross section data
143 for ( i = 0 ; i < nDiscrete ; i++ )
144 {
145 vct_gammas_par.push_back( theGammas[ i ] );
146 vct_shells_par.push_back( theShells[ i ] );
147 vct_primary_par.push_back( isPrimary[ i ] );
148 vct_distype_par.push_back( disType[ i ] );
149 G4NeutronHPVector* hpv = new G4NeutronHPVector;
150 *hpv = thePartialXsec[ i ];
151 vct_pXS_par.push_back( hpv );
152 }
153 }
154 if ( theGammas == NULL ) theGammas = new G4double[nDiscrete2];
155 if ( theShells == NULL ) theShells = new G4double[nDiscrete2];
156//080731
157
158 for (i=0; i< nIso; i++) // isotropic photons
159 {
160 aDataFile >> theGammas[i] >> theShells[i];
161 theGammas[i]*=eV;
162 theShells[i]*=eV;
163 }
164 nNeu = new G4int [nDiscrete2-nIso];
165 if(tabulationType==1)theLegendre=new G4NeutronHPLegendreTable *[nDiscrete2-nIso];
166 if(tabulationType==2)theAngular =new G4NeutronHPAngularP *[nDiscrete2-nIso];
167 for(i=nIso; i< nDiscrete2; i++)
168 {
169 if(tabulationType==1)
170 {
171 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
172 theGammas[i]*=eV;
173 theShells[i]*=eV;
174 theLegendre[i-nIso]=new G4NeutronHPLegendreTable[nNeu[i-nIso]];
175 theLegendreManager.Init(aDataFile);
176 for (ii=0; ii<nNeu[i-nIso]; ii++)
177 {
178 theLegendre[i-nIso][ii].Init(aDataFile);
179 }
180 }
181 else if(tabulationType==2)
182 {
183 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
184 theGammas[i]*=eV;
185 theShells[i]*=eV;
186 theAngular[i-nIso]=new G4NeutronHPAngularP[nNeu[i-nIso]];
187 for (ii=0; ii<nNeu[i-nIso]; ii++)
188 {
189 theAngular[i-nIso][ii].Init(aDataFile);
190 }
191 }
192 else
193 {
194 G4cout << "tabulation type: tabulationType"<<G4endl;
195 throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions.");
196 }
197 }
198//080731
199 if ( vct_gammas_par.size() > 0 )
200 {
201 //Reordering cross section data to corrsponding distribution data
202 for ( i = 0 ; i < nDiscrete ; i++ )
203 {
204 for ( G4int j = 0 ; j < nDiscrete ; j++ )
205 {
206 // Checking gamma and shell to identification
207 if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] )
208 {
209 isPrimary [ i ] = vct_primary_par [ j ];
210 disType [ i ] = vct_distype_par [ j ];
211 thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) );
212 }
213 }
214 }
215 //Garbage collection
216 for ( std::vector < G4NeutronHPVector* >::iterator
217 it = vct_pXS_par.begin() ; it != vct_pXS_par.end() ; it++ )
218 {
219 delete *it;
220 }
221 }
222//080731
223 }
224}
225
226
227void G4NeutronHPPhotonDist::InitEnergies(std::ifstream & aDataFile)
228{
229 G4int i, energyDistributionsNeeded = 0;
230 for (i=0; i<nDiscrete; i++)
231 {
232 if( disType[i]==1) energyDistributionsNeeded =1;
233 }
234 if(!energyDistributionsNeeded) return;
235 aDataFile >> nPartials;
236 distribution = new G4int[nPartials];
237 probs = new G4NeutronHPVector[nPartials];
238 partials = new G4NeutronHPPartial * [nPartials];
239 G4int nen;
240 G4int dummy;
241 for (i=0; i<nPartials; i++)
242 {
243 aDataFile >> dummy;
244 probs[i].Init(aDataFile, eV);
245 aDataFile >> nen;
246 partials[i] = new G4NeutronHPPartial(nen);
247 partials[i]->InitInterpolation(aDataFile);
248 partials[i]->Init(aDataFile);
249 }
250}
251
252void G4NeutronHPPhotonDist::InitPartials(std::ifstream & aDataFile)
253{
254
255 //G4cout << "G4NeutronHPPhotonDist::InitPartials " << G4endl;
256 aDataFile >> nDiscrete >> targetMass;
257 if(nDiscrete != 1)
258 {
259 theTotalXsec.Init(aDataFile, eV);
260 }
261 G4int i;
262 theGammas = new G4double[nDiscrete];
263 theShells = new G4double[nDiscrete];
264 isPrimary = new G4int[nDiscrete];
265 disType = new G4int[nDiscrete];
266 thePartialXsec = new G4NeutronHPVector[nDiscrete];
267 for(i=0; i<nDiscrete; i++)
268 {
269 aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
270 theGammas[i]*=eV;
271 theShells[i]*=eV;
272 thePartialXsec[i].Init(aDataFile, eV);
273 }
274
275 //G4cout << "G4NeutronHPPhotonDist::InitPartials Test " << G4endl;
276 //G4cout << "G4NeutronHPPhotonDist::InitPartials nDiscrete " << nDiscrete << G4endl;
277 //G4NeutronHPVector* aHP = new G4NeutronHPVector;
278 //aHP->Check(1);
279}
280
281G4ReactionProductVector * G4NeutronHPPhotonDist::GetPhotons(G4double anEnergy)
282{
283
284 //G4cout << "G4NeutronHPPhotonDist::GetPhotons repFlag " << repFlag << G4endl;
285 // the partial cross-section case is not in this yet. @@@@ << 070601 TK add partial
286 G4int i, ii, iii;
287 G4int nSecondaries = 0;
288 G4ReactionProductVector * thePhotons = new G4ReactionProductVector;
289 if(repFlag==1)
290 {
291 G4double current=0;
292 for(i=0; i<nDiscrete; i++)
293 {
294 current = theYield[i].GetY(anEnergy);
295 actualMult[i] = G4Poisson(current); // max cut-off still missing @@@
296 if(nDiscrete==1&&current<1.0001)
297 {
298 actualMult[i] = static_cast<G4int>(current);
299 if(current<1)
300 {
301 actualMult[i] = 0;
302 if(G4UniformRand()<current) actualMult[i] = 1;
303 }
304 }
305 nSecondaries += actualMult[i];
306 }
307 //G4cout << "nSecondaries " << nSecondaries << " anEnergy " << anEnergy/eV << G4endl;
308 for(i=0;i<nSecondaries;i++)
309 {
310 G4ReactionProduct * theOne = new G4ReactionProduct;
311 theOne->SetDefinition(G4Gamma::Gamma());
312 thePhotons->push_back(theOne);
313 }
314 G4int count=0;
315
316/*
317G4double totalCascadeEnergy = 0.;
318G4double lastCascadeEnergy = 0.;
319G4double eGamm = 0;
320G4int maxEnergyIndex = 0;
321*/
322 //Gcout << "nDiscrete " << nDiscrete << " nPartials " << nPartials << G4endl;
323//3456
324 if ( nDiscrete == 1 && nPartials == 1 )
325 {
326 if ( actualMult[ 0 ] > 0 )
327 {
328 if ( disType[0] == 1 ) // continuum
329 {
330
331/*
332 for(ii=0; ii< actualMult[0]; ii++)
333 {
334
335 G4double sum=0, run=0;
336 for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
337 G4double random = G4UniformRand();
338 G4int theP = 0;
339 for(iii=0; iii<nPartials; iii++)
340 {
341 run+=probs[iii].GetY(anEnergy);
342 theP = iii;
343 if(random<run/sum) break;
344 }
345 if(theP==nPartials) theP=nPartials-1; // das sortiert J aus.
346 sum=0;
347 G4NeutronHPVector * temp;
348 temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
349 // Looking for TotalCascdeEnergy or LastMaxEnergy
350 if (ii == 0)
351 {
352 maxEnergyIndex = temp->GetVectorLength()-1;
353 totalCascadeEnergy = temp->GetX(maxEnergyIndex);
354 lastCascadeEnergy = totalCascadeEnergy;
355 }
356 lastCascadeEnergy -= eGamm;
357 if (ii != actualMult[i]-1) eGamm = temp->SampleWithMax(lastCascadeEnergy);
358 else eGamm = lastCascadeEnergy;
359 thePhotons->operator[](count)->SetKineticEnergy(eGamm);
360 delete temp;
361
362 }
363*/
364 G4NeutronHPVector * temp;
365 temp = partials[ 0 ]->GetY(anEnergy); //@@@ look at, seems fishy
366 G4double maximumE = temp->GetX( temp->GetVectorLength()-1 ); // This is an assumption.
367
368 //G4cout << "start " << actualMult[ 0 ] << " maximumE " << maximumE/eV << G4endl;
369
370 std::vector< G4double > photons_e_best( actualMult[ 0 ] , 0.0 );
371 G4double best = DBL_MAX;
372 G4int maxTry = 1000;
373 for ( G4int j = 0 ; j < maxTry ; j++ )
374 {
375 std::vector< G4double > photons_e( actualMult[ 0 ] , 0.0 );
376 for ( std::vector< G4double >::iterator
377 it = photons_e.begin() ; it < photons_e.end() ; it++ )
378 {
379 *it = temp->Sample();
380 }
381 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) > maximumE )
382 {
383 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) < best )
384 photons_e_best = photons_e;
385 continue;
386 }
387 else
388 {
389 for ( std::vector< G4double >::iterator
390 it = photons_e.begin() ; it < photons_e.end() ; it++ )
391 {
392 thePhotons->operator[](count)->SetKineticEnergy( *it );
393 }
394 //G4cout << "OK " << actualMult[0] << " j " << j << " total photons E "
395 // << std::accumulate( photons_e.begin() , photons_e.end() , 0.0 )/eV << " ratio " << std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) / maximumE
396 // << G4endl;
397
398 break;
399 }
400 G4cout << "NeutronHPPhotonDist could not find fitted energy set for multiplicity of " << actualMult[0] << "." << G4endl;
401 G4cout << "NeutronHPPhotonDist will use the best set." << G4endl;
402 for ( std::vector< G4double >::iterator
403 it = photons_e_best.begin() ; it < photons_e_best.end() ; it++ )
404 {
405 thePhotons->operator[](count)->SetKineticEnergy( *it );
406 }
407 //G4cout << "Not Good " << actualMult[0] << " j " << j << " total photons E "
408 // << best/eV << " ratio " << best / maximumE
409 // << G4endl;
410 }
411 // TKDB
412 delete temp;
413 }
414 else // discrete
415 {
416 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
417 }
418 count++;
419 if(count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
420 }
421
422 }
423 else
424 {
425 for(i=0; i<nDiscrete; i++)
426 {
427 for(ii=0; ii< actualMult[i]; ii++)
428 {
429 if(disType[i]==1) // continuum
430 {
431 G4double sum=0, run=0;
432 for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
433 G4double random = G4UniformRand();
434 G4int theP = 0;
435 for(iii=0; iii<nPartials; iii++)
436 {
437 run+=probs[iii].GetY(anEnergy);
438 theP = iii;
439 if(random<run/sum) break;
440 }
441 if(theP==nPartials) theP=nPartials-1; // das sortiert J aus.
442 sum=0;
443 G4NeutronHPVector * temp;
444 temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
445 G4double eGamm = temp->Sample();
446 thePhotons->operator[](count)->SetKineticEnergy(eGamm);
447 delete temp;
448 }
449 else // discrete
450 {
451 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
452 }
453 count++;
454 if(count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
455 }
456 }
457 }
458 // now do the angular distributions...
459 if( isoFlag == 1)
460 {
461 for (i=0; i< nSecondaries; i++)
462 {
463 G4double costheta = 2.*G4UniformRand()-1;
464 G4double theta = std::acos(costheta);
465 G4double phi = twopi*G4UniformRand();
466 G4double sinth = std::sin(theta);
467 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
468 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
469 thePhotons->operator[](i)->SetMomentum( temp ) ;
470 // G4cout << "Isotropic distribution in PhotonDist"<<temp<<G4endl;
471 }
472 }
473 else
474 {
475 for(i=0; i<nSecondaries; i++)
476 {
477 G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
478 for(ii=0; ii<nDiscrete2; ii++)
479 {
480 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
481 }
482 if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
483 if(ii<nIso)
484 {
485 // isotropic distribution
486 G4double theta = pi*G4UniformRand();
487 G4double phi = twopi*G4UniformRand();
488 G4double sinth = std::sin(theta);
489 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
490 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
491 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
492 }
493 else if(tabulationType==1)
494 {
495 // legendre polynomials
496 G4int it(0);
497 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
498 {
499 it = iii;
500 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
501 break;
502 }
503 G4NeutronHPLegendreStore aStore(2);
504 aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));
505 //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
506 //TKDB 110512
507 if ( it > 0 )
508 {
509 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
510 }
511 else
512 {
513 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it]));
514 }
515 G4double cosTh = aStore.SampleMax(anEnergy);
516 G4double theta = std::acos(cosTh);
517 G4double phi = twopi*G4UniformRand();
518 G4double sinth = std::sin(theta);
519 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
520 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
521 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
522 }
523 else
524 {
525 // tabulation of probabilities.
526 G4int it(0);
527 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
528 {
529 it = iii;
530 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
531 break;
532 }
533 G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@
534 G4double theta = std::acos(costh);
535 G4double phi = twopi*G4UniformRand();
536 G4double sinth = std::sin(theta);
537 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
538 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
539 thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
540 }
541 }
542 }
543 }
544 else if(repFlag == 2)
545 {
546 G4double * running = new G4double[nGammaEnergies];
547 running[0]=theTransitionProbabilities[0];
548 //G4int i; //declaration at 284th
549 for(i=1; i<nGammaEnergies; i++)
550 {
551 running[i]=running[i-1]+theTransitionProbabilities[i];
552 }
553 G4double random = G4UniformRand();
554 G4int it=0;
555 for(i=0; i<nGammaEnergies; i++)
556 {
557 it = i;
558 if(random < running[i]/running[nGammaEnergies-1]) break;
559 }
560 delete [] running;
561 G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
562 G4ReactionProduct * theOne = new G4ReactionProduct;
563 theOne->SetDefinition(G4Gamma::Gamma());
564 random = G4UniformRand();
565 if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
566 {
567 theOne->SetDefinition(G4Electron::Electron());
568 //Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25, 2009
569 //But never enter at least with G4NDL3.13
570 totalEnergy += G4Electron::Electron()->GetPDGMass(); //proposed correction: add this line for electron
571 }
572 theOne->SetTotalEnergy(totalEnergy);
573 if( isoFlag == 1 )
574 {
575 G4double costheta = 2.*G4UniformRand()-1;
576 G4double theta = std::acos(costheta);
577 G4double phi = twopi*G4UniformRand();
578 G4double sinth = std::sin(theta);
579 //Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25, 2009
580 //G4double en = theOne->GetTotalEnergy();
581 G4double en = theOne->GetTotalMomentum();
582 //But never cause real effect at least with G4NDL3.13 TK
583 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
584 theOne->SetMomentum( temp ) ;
585 }
586 else
587 {
588 G4double currentEnergy = theOne->GetTotalEnergy();
589 for(ii=0; ii<nDiscrete2; ii++)
590 {
591 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
592 }
593 if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
594 if(ii<nIso)
595 {
596 //Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25, 2009
597 // isotropic distribution
598 //G4double theta = pi*G4UniformRand();
599 G4double theta = std::acos(2.*G4UniformRand()-1.);
600 //But this is alos never cause real effect at least with G4NDL3.13 TK not repFlag == 2 AND isoFlag != 1
601 G4double phi = twopi*G4UniformRand();
602 G4double sinth = std::sin(theta);
603 //Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25, 2009
604 //G4double en = theOne->GetTotalEnergy();
605 G4double en = theOne->GetTotalMomentum();
606 //But never cause real effect at least with G4NDL3.13 TK
607 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
608 theOne->SetMomentum( tempVector ) ;
609 }
610 else if(tabulationType==1)
611 {
612 // legendre polynomials
613 G4int itt(0);
614 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
615 {
616 itt = iii;
617 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
618 break;
619 }
620 G4NeutronHPLegendreStore aStore(2);
621 aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt]));
622 //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
623 //TKDB 110512
624 if ( itt > 0 )
625 {
626 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1]));
627 }
628 else
629 {
630 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt]));
631 }
632 G4double cosTh = aStore.SampleMax(anEnergy);
633 G4double theta = std::acos(cosTh);
634 G4double phi = twopi*G4UniformRand();
635 G4double sinth = std::sin(theta);
636 //Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25, 2009
637 //G4double en = theOne->GetTotalEnergy();
638 G4double en = theOne->GetTotalMomentum();
639 //But never cause real effect at least with G4NDL3.13 TK
640 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
641 theOne->SetMomentum( tempVector ) ;
642 }
643 else
644 {
645 // tabulation of probabilities.
646 G4int itt(0);
647 for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
648 {
649 itt = iii;
650 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
651 break;
652 }
653 G4double costh = theAngular[ii-nIso][itt].GetCosTh(); // no interpolation yet @@
654 G4double theta = std::acos(costh);
655 G4double phi = twopi*G4UniformRand();
656 G4double sinth = std::sin(theta);
657 //Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25, 2009
658 //G4double en = theOne->GetTotalEnergy();
659 G4double en = theOne->GetTotalMomentum();
660 //But never cause real effect at least with G4NDL3.13 TK
661 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
662 theOne->SetMomentum( tmpVector ) ;
663 }
664 }
665 thePhotons->push_back(theOne);
666 }
667 else if( repFlag==0 )
668 {
669
670// TK add
671 if ( thePartialXsec == 0 )
672 {
673 //G4cout << "repFlag is 0, but no PartialXsec data" << G4endl;
674 //G4cout << "This is not support yet." << G4endl;
675 return thePhotons;
676 }
677
678// Partial Case
679
680 G4ReactionProduct * theOne = new G4ReactionProduct;
681 theOne->SetDefinition( G4Gamma::Gamma() );
682 thePhotons->push_back( theOne );
683
684// Energy
685
686 //G4cout << "Partial Case nDiscrete " << nDiscrete << G4endl;
687 G4double sum = 0.0;
688 std::vector < G4double > dif( nDiscrete , 0.0 );
689 for ( G4int j = 0 ; j < nDiscrete ; j++ )
690 {
691 G4double x = thePartialXsec[ j ].GetXsec( anEnergy ); // x in barn
692 if ( x > 0 )
693 {
694 sum += x;
695 }
696 dif [ j ] = sum;
697 //G4cout << "j " << j << ", x " << x << ", dif " << dif [ j ] << G4endl;
698 }
699
700 G4double rand = G4UniformRand();
701
702 G4int iphoton = 0;
703 for ( G4int j = 0 ; j < nDiscrete ; j++ )
704 {
705 G4double y = rand*sum;
706 if ( dif [ j ] > y )
707 {
708 iphoton = j;
709 break;
710 }
711 }
712 //G4cout << "iphoton " << iphoton << G4endl;
713 //G4cout << "photon energy " << theGammas[ iphoton ] /eV << G4endl;
714
715// Angle
716 G4double cosTheta = 0.0; // mu
717
718 if ( isoFlag == 1 )
719 {
720
721// Isotropic Case
722
723 cosTheta = 2.*G4UniformRand()-1;
724
725 }
726 else
727 {
728
729 if ( iphoton < nIso )
730 {
731
732// still Isotropic
733
734 cosTheta = 2.*G4UniformRand()-1;
735
736 }
737 else
738 {
739
740 //G4cout << "Not Isotropic and isoFlag " << isoFlag << G4endl;
741 //G4cout << "tabulationType " << tabulationType << G4endl;
742 //G4cout << "nDiscrete2 " << nDiscrete2 << G4endl;
743 //G4cout << "nIso " << nIso << G4endl;
744 //G4cout << "size of nNeu " << nDiscrete2-nIso << G4endl;
745 //G4cout << "nNeu[iphoton-nIso] " << nNeu[iphoton-nIso] << G4endl;
746
747 if ( tabulationType == 1 )
748 {
749// legendre polynomials
750
751 G4int iangle = 0;
752 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
753 {
754 iangle = j;
755 if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
756 }
757
758 G4NeutronHPLegendreStore aStore( 2 );
759 aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );
760 aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) );
761
762 cosTheta = aStore.SampleMax( anEnergy );
763
764 }
765 else if ( tabulationType == 2 )
766 {
767
768// tabulation of probabilities.
769
770 G4int iangle = 0;
771 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
772 {
773 iangle = j;
774 if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
775 }
776
777 cosTheta = theAngular[iphoton-nIso][ iangle ].GetCosTh(); // no interpolation yet @@
778
779 }
780 }
781 }
782
783// Set
784 G4double phi = twopi*G4UniformRand();
785 G4double theta = std::acos( cosTheta );
786 G4double sinTheta = std::sin( theta );
787
788 G4double photonE = theGammas[ iphoton ];
789 G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
790 G4ThreeVector photonP = photonE * direction;
791 thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
792
793 }
794 else
795 {
796 delete thePhotons;
797 thePhotons = 0; // no gamma data available; some work needed @@@@@@@
798 }
799 return thePhotons;
800}
801
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
G4ReactionProductVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
Definition: G4ReactionProductVector.hh:43
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:85
G4NeutronHPAngularP::Init
void Init(std::ifstream &aDataFile)
Definition: G4NeutronHPAngularP.hh:52
G4NeutronHPAngularP::GetCosTh
G4double GetCosTh(G4int l)
Definition: G4NeutronHPAngularP.hh:80
G4NeutronHPLegendreStore::SetCoeff
void SetCoeff(G4int i, G4int l, G4double coeff)
Definition: G4NeutronHPLegendreStore.hh:59
G4NeutronHPLegendreStore::SampleMax
G4double SampleMax(G4double energy)
Definition: G4NeutronHPLegendreStore.cc:117
G4NeutronHPLegendreTable::Init
void Init(std::ifstream &aDataFile)
Definition: G4NeutronHPLegendreTable.hh:67
G4NeutronHPPartial::Init
void Init(std::ifstream &aDataFile)
Definition: G4NeutronHPPartial.hh:71
G4NeutronHPPartial::InitInterpolation
void InitInterpolation(G4int i, std::ifstream &aDataFile)
Definition: G4NeutronHPPartial.hh:61
G4NeutronHPPartial::GetY
G4double GetY(G4int i, G4int j)
Definition: G4NeutronHPPartial.hh:116
G4NeutronHPPhotonDist::InitPartials
void InitPartials(std::ifstream &aDataFile)
Definition: G4NeutronHPPhotonDist.cc:252
G4NeutronHPPhotonDist::InitMean
G4bool InitMean(std::ifstream &aDataFile)
Definition: G4NeutronHPPhotonDist.cc:54
G4NeutronHPPhotonDist::InitEnergies
void InitEnergies(std::ifstream &aDataFile)
Definition: G4NeutronHPPhotonDist.cc:227
G4NeutronHPPhotonDist::GetPhotons
G4ReactionProductVector * GetPhotons(G4double anEnergy)
Definition: G4NeutronHPPhotonDist.cc:281
G4NeutronHPPhotonDist::InitAngular
void InitAngular(std::ifstream &aDataFile)
Definition: G4NeutronHPPhotonDist.cc:121
G4NeutronHPVector::Init
void Init(std::ifstream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
Definition: G4NeutronHPVector.hh:215
G4NeutronHPVector::GetVectorLength
G4int GetVectorLength() const
Definition: G4NeutronHPVector.hh:193
G4NeutronHPVector::GetX
G4double GetX(G4int i) const
Definition: G4NeutronHPVector.hh:121
G4NeutronHPVector::GetXsec
G4double GetXsec(G4int i)
Definition: G4NeutronHPVector.hh:120
G4NeutronHPVector::GetY
G4double GetY(G4double x)
Definition: G4NeutronHPVector.hh:192
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:166
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:155
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83