Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
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 
54 G4bool 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 
121 void G4NeutronHPPhotonDist::InitAngular(std::ifstream & aDataFile)
122 {
123 
124  G4int i, ii;
125  //angular distributions
126  aDataFile >> isoFlag;
127  if (isoFlag != 1)
128  {
129 if ( 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 
227 void 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 
252 void 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 
281 G4ReactionProductVector * 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 /*
317 G4double totalCascadeEnergy = 0.;
318 G4double lastCascadeEnergy = 0.;
319 G4double eGamm = 0;
320 G4int 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 (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) 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 (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) 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 (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) 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 (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) 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 (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) 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 (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) 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 
Generated on Sat May 25 2013 14:34:05 for Geant4 by   doxygen 1.8.4