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