ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPPhotonDist.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPPhotonDist.cc
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 
48 #include "G4ParticleHPPhotonDist.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4ParticleHPLegendreStore.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  G4int i, ii;
125  //angular distributions
126  aDataFile >> isoFlag;
127  if (isoFlag != 1)
128  {
129  if (repFlag == 2) G4cout << "G4ParticleHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 HyperNews. " << G4endl;
130  aDataFile >> tabulationType >> nDiscrete2 >> nIso;
131 //080731
132  if (theGammas != NULL && nDiscrete2 != nDiscrete)
133  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
136  // (InitAngular here) data are different, we have to re-coordinate
137  // consistent data order.
138  std::vector < G4double > vct_gammas_par;
139  std::vector < G4double > vct_shells_par;
140  std::vector < G4int > vct_primary_par;
141  std::vector < G4int > vct_distype_par;
142  std::vector < G4ParticleHPVector* > vct_pXS_par;
143  if ( theGammas != NULL && theShells != NULL )
144  {
145  //copy the cross section data
146  for ( i = 0 ; i < nDiscrete ; i++ )
147  {
148  vct_gammas_par.push_back( theGammas[ i ] );
149  vct_shells_par.push_back( theShells[ i ] );
150  vct_primary_par.push_back( isPrimary[ i ] );
151  vct_distype_par.push_back( disType[ i ] );
152  G4ParticleHPVector* hpv = new G4ParticleHPVector;
153  *hpv = thePartialXsec[ i ];
154  vct_pXS_par.push_back( hpv );
155  }
156  }
157  if ( theGammas == NULL ) theGammas = new G4double[nDiscrete2];
158  if ( theShells == NULL ) theShells = new G4double[nDiscrete2];
159 //080731
160 
161  for (i=0; i< nIso; i++) // isotropic photons
162  {
163  aDataFile >> theGammas[i] >> theShells[i];
164  theGammas[i]*=eV;
165  theShells[i]*=eV;
166  }
167  nNeu = new G4int [nDiscrete2-nIso];
168  if(tabulationType==1)theLegendre=new G4ParticleHPLegendreTable *[nDiscrete2-nIso];
169  if(tabulationType==2)theAngular =new G4ParticleHPAngularP *[nDiscrete2-nIso];
170  for(i=nIso; i< nDiscrete2; i++)
171  {
172  if(tabulationType==1)
173  {
174  aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
175  theGammas[i]*=eV;
176  theShells[i]*=eV;
177  theLegendre[i-nIso]=new G4ParticleHPLegendreTable[nNeu[i-nIso]];
178  theLegendreManager.Init(aDataFile);
179  for (ii=0; ii<nNeu[i-nIso]; ii++)
180  {
181  theLegendre[i-nIso][ii].Init(aDataFile);
182  }
183  }
184  else if(tabulationType==2)
185  {
186  aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
187  theGammas[i]*=eV;
188  theShells[i]*=eV;
189  theAngular[i-nIso]=new G4ParticleHPAngularP[nNeu[i-nIso]];
190  for (ii=0; ii<nNeu[i-nIso]; ii++)
191  {
192  theAngular[i-nIso][ii].Init(aDataFile);
193  }
194  }
195  else
196  {
197  G4cout << "tabulation type: tabulationType"<<G4endl;
198  throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions.");
199  }
200  }
201 //080731
202  if ( vct_gammas_par.size() > 0 )
203  {
204  //Reordering cross section data to corrsponding distribution data
205  for ( i = 0 ; i < nDiscrete ; i++ )
206  {
207  for ( G4int j = 0 ; j < nDiscrete ; j++ )
208  {
209  // Checking gamma and shell to identification
210  if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] )
211  {
212  isPrimary [ i ] = vct_primary_par [ j ];
213  disType [ i ] = vct_distype_par [ j ];
214  thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) );
215  }
216  }
217  }
218  //Garbage collection
219  for ( std::vector < G4ParticleHPVector* >::iterator
220  it = vct_pXS_par.begin() ; it != vct_pXS_par.end() ; it++ )
221  {
222  delete *it;
223  }
224  }
225 //080731
226  }
227 }
228 
229 
230 void G4ParticleHPPhotonDist::InitEnergies(std::istream & aDataFile)
231 {
232  G4int i, energyDistributionsNeeded = 0;
233  for (i=0; i<nDiscrete; i++)
234  {
235  if( disType[i]==1) energyDistributionsNeeded =1;
236  }
237  if(!energyDistributionsNeeded) return;
238  aDataFile >> nPartials;
239  distribution = new G4int[nPartials];
240  probs = new G4ParticleHPVector[nPartials];
241  partials = new G4ParticleHPPartial * [nPartials];
242  G4int nen;
243  G4int dummy;
244  for (i=0; i<nPartials; i++)
245  {
246  aDataFile >> dummy;
247  probs[i].Init(aDataFile, eV);
248  aDataFile >> nen;
249  partials[i] = new G4ParticleHPPartial(nen);
250  partials[i]->InitInterpolation(aDataFile);
251  partials[i]->Init(aDataFile);
252  }
253 }
254 
255 void G4ParticleHPPhotonDist::InitPartials(std::istream& aDataFile,
256  G4ParticleHPVector* theXsec)
257 {
258  if (theXsec) theReactionXsec = theXsec;
259 
260  aDataFile >> nDiscrete >> targetMass;
261  if(nDiscrete != 1)
262  {
263  theTotalXsec.Init(aDataFile, eV);
264  }
265  G4int i;
266  theGammas = new G4double[nDiscrete];
267  theShells = new G4double[nDiscrete];
268  isPrimary = new G4int[nDiscrete];
269  disType = new G4int[nDiscrete];
270  thePartialXsec = new G4ParticleHPVector[nDiscrete];
271  for(i=0; i<nDiscrete; i++)
272  {
273  aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
274  theGammas[i]*=eV;
275  theShells[i]*=eV;
276  thePartialXsec[i].Init(aDataFile, eV);
277  }
278 
279  //G4cout << "G4ParticleHPPhotonDist::InitPartials Test " << G4endl;
280  //G4cout << "G4ParticleHPPhotonDist::InitPartials nDiscrete " << nDiscrete << G4endl;
281  //G4ParticleHPVector* aHP = new G4ParticleHPVector;
282  //aHP->Check(1);
283 }
284 
285 G4ReactionProductVector * G4ParticleHPPhotonDist::GetPhotons(G4double anEnergy)
286 {
287 
288  //G4cout << "G4ParticleHPPhotonDist::GetPhotons repFlag " << repFlag << G4endl;
289  // the partial cross-section case is not in this yet. @@@@ << 070601 TK add partial
290  if ( actualMult.Get() == NULL ) {
291  actualMult.Get() = new std::vector<G4int>( nDiscrete );
292  }
293  G4int i, ii, iii;
294  G4int nSecondaries = 0;
295  G4ReactionProductVector* thePhotons = new G4ReactionProductVector;
296 
297  if (repFlag==1) {
298  G4double current=0;
299  for (i = 0; i < nDiscrete; i++) {
300  current = theYield[i].GetY(anEnergy);
301  actualMult.Get()->at(i) = G4Poisson(current); // max cut-off still missing @@@
302  if (nDiscrete == 1 && current < 1.0001) {
303  actualMult.Get()->at(i) = static_cast<G4int>(current);
304  if(current<1)
305  {
306  actualMult.Get()->at(i) = 0;
307  if(G4UniformRand()<current) actualMult.Get()->at(i) = 1;
308  }
309  }
310  nSecondaries += actualMult.Get()->at(i);
311  }
312  //G4cout << "nSecondaries " << nSecondaries << " anEnergy " << anEnergy/eV << G4endl;
313  for (i = 0; i < nSecondaries; i++) {
314  G4ReactionProduct * theOne = new G4ReactionProduct;
315  theOne->SetDefinition(G4Gamma::Gamma());
316  thePhotons->push_back(theOne);
317  }
318 
319  G4int count = 0;
320 
321  if (nDiscrete == 1 && nPartials == 1) {
322  if (actualMult.Get()->at(0) > 0) {
323  if (disType[0] == 1) {
324  // continuum
325  G4ParticleHPVector* temp;
326  temp = partials[ 0 ]->GetY(anEnergy); //@@@ look at, seems fishy
327  G4double maximumE = temp->GetX( temp->GetVectorLength()-1 ); // This is an assumption.
328 
329  //G4cout << "start " << actualMult[ 0 ] << " maximumE " << maximumE/eV << G4endl;
330 
331  std::vector< G4double > photons_e_best( actualMult.Get()->at(0) , 0.0 );
332  G4double best = DBL_MAX;
333  G4int maxTry = 1000;
334  for (G4int j = 0; j < maxTry; j++) {
335  std::vector<G4double> photons_e(actualMult.Get()->at(0), 0.0);
336  for (std::vector<G4double>::iterator it = photons_e.begin(); it < photons_e.end(); it++) {
337  *it = temp->Sample();
338  }
339 
340  if (std::accumulate(photons_e.begin(), photons_e.end(), 0.0) > maximumE) {
341  if (std::accumulate(photons_e.begin(), photons_e.end(), 0.0) < best)
342  photons_e_best = photons_e;
343  continue;
344 
345  } else {
346  G4int iphot = 0;
347  for (std::vector<G4double>::iterator it = photons_e.begin(); it < photons_e.end(); it++) {
348  thePhotons->operator[](iphot)->SetKineticEnergy(*it); // Replace index count, which was not incremented,
349  // with iphot, which is, as per Artem Zontikov,
350  // bug report 2167
351  iphot++;
352  }
353  // G4cout << "OK " << actualMult[0] << " j " << j << " total photons E "
354  // << std::accumulate(photons_e.begin(), photons_e.end(), 0.0)/eV << " ratio "
355  // << std::accumulate(photons_e.begin(), photons_e.end(), 0.0)/maximumE
356  // << G4endl;
357 
358  break;
359  }
360  //G4cout << "Not Good " << actualMult[0] << " j " << j << " total photons E "
361  // << best/eV << " ratio " << best / maximumE
362  // << G4endl;
363  }
364  // TKDB
365  delete temp;
366 
367  } else {
368  // discrete
369  thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
370  }
371  count++;
372  if (count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPPhotonDist::GetPhotons inconsistancy");
373  }
374 
375  } else { // nDiscrete != 1 or nPartials != 1
376  for (i=0; i<nDiscrete; i++) {
377  for (ii=0; ii< actualMult.Get()->at(i); ii++) {
378  if (disType[i] == 1) {
379  // continuum
380  G4double sum=0, run=0;
381  for (iii = 0; iii < nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
382  G4double random = G4UniformRand();
383  G4int theP = 0;
384  for (iii = 0; iii < nPartials; iii++) {
385  run+=probs[iii].GetY(anEnergy);
386  theP = iii;
387  if(random<run/sum) break;
388  }
389 
390  if (theP == nPartials) theP=nPartials-1; // das sortiert J aus.
391  sum = 0;
392  G4ParticleHPVector * temp;
393  temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
394  G4double eGamm = temp->Sample();
395  thePhotons->operator[](count)->SetKineticEnergy(eGamm);
396  delete temp;
397 
398  } else {
399  // discrete
400  thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
401  }
402  count++;
403  if (count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPPhotonDist::GetPhotons inconsistancy");
404  }
405  }
406  }
407 
408  // now do the angular distributions...
409  if (isoFlag == 1) {
410  for (i=0; i< nSecondaries; i++)
411  {
412  G4double costheta = 2.*G4UniformRand()-1;
413  G4double theta = std::acos(costheta);
414  G4double phi = twopi*G4UniformRand();
415  G4double sinth = std::sin(theta);
416  G4double en = thePhotons->operator[](i)->GetTotalEnergy();
417  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
418  thePhotons->operator[](i)->SetMomentum( temp ) ;
419  // G4cout << "Isotropic distribution in PhotonDist"<<temp<<G4endl;
420  }
421  }
422  else
423  {
424  for(i=0; i<nSecondaries; i++)
425  {
426  G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
427  for(ii=0; ii<nDiscrete2; ii++)
428  {
429  if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
430  }
431  if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
432  if(ii<nIso)
433  {
434  // isotropic distribution
435  //
436  //Fix Bugzilla report #1745
437  //G4double theta = pi*G4UniformRand();
438  G4double costheta = 2.*G4UniformRand()-1;
439  G4double theta = std::acos(costheta);
440  G4double phi = twopi*G4UniformRand();
441  G4double sinth = std::sin(theta);
442  G4double en = thePhotons->operator[](i)->GetTotalEnergy();
443  // DHW G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
444  G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costheta );
445  thePhotons->operator[](i)->SetMomentum( tempVector ) ;
446  }
447  else if(tabulationType==1)
448  {
449  // legendre polynomials
450  G4int it(0);
451  for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
452  {
453  it = iii;
454  if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
455  break;
456  }
457  G4ParticleHPLegendreStore aStore(2);
458  aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));
459  //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
460  //TKDB 110512
461  if ( it > 0 )
462  {
463  aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
464  }
465  else
466  {
467  aStore.SetCoeff(0, &(theLegendre[ii-nIso][it]));
468  }
469  G4double cosTh = aStore.SampleMax(anEnergy);
470  G4double theta = std::acos(cosTh);
471  G4double phi = twopi*G4UniformRand();
472  G4double sinth = std::sin(theta);
473  G4double en = thePhotons->operator[](i)->GetTotalEnergy();
474  G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
475  thePhotons->operator[](i)->SetMomentum( tempVector ) ;
476  }
477  else
478  {
479  // tabulation of probabilities.
480  G4int it(0);
481  for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
482  {
483  it = iii;
484  if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
485  break;
486  }
487  G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@
488  G4double theta = std::acos(costh);
489  G4double phi = twopi*G4UniformRand();
490  G4double sinth = std::sin(theta);
491  G4double en = thePhotons->operator[](i)->GetTotalEnergy();
492  G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
493  thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
494  }
495  }
496  }
497 
498  } else if (repFlag == 2) {
499  G4double * running = new G4double[nGammaEnergies];
500  running[0]=theTransitionProbabilities[0];
501  //G4int i; //declaration at 284th
502  for(i=1; i<nGammaEnergies; i++)
503  {
504  running[i]=running[i-1]+theTransitionProbabilities[i];
505  }
506  G4double random = G4UniformRand();
507  G4int it=0;
508  for(i=0; i<nGammaEnergies; i++)
509  {
510  it = i;
511  if(random < running[i]/running[nGammaEnergies-1]) break;
512  }
513  delete [] running;
514  G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
515  G4ReactionProduct * theOne = new G4ReactionProduct;
516  theOne->SetDefinition(G4Gamma::Gamma());
517  random = G4UniformRand();
518  if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
519  {
520  theOne->SetDefinition(G4Electron::Electron());
521  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
522  //But never enter at least with G4NDL3.13
523  totalEnergy += G4Electron::Electron()->GetPDGMass(); //proposed correction: add this line for electron
524  }
525  theOne->SetTotalEnergy(totalEnergy);
526  if( isoFlag == 1 )
527  {
528  G4double costheta = 2.*G4UniformRand()-1;
529  G4double theta = std::acos(costheta);
530  G4double phi = twopi*G4UniformRand();
531  G4double sinth = std::sin(theta);
532  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
533  //G4double en = theOne->GetTotalEnergy();
534  G4double en = theOne->GetTotalMomentum();
535  //But never cause real effect at least with G4NDL3.13 TK
536  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
537  theOne->SetMomentum( temp ) ;
538  }
539  else
540  {
541  G4double currentEnergy = theOne->GetTotalEnergy();
542  for(ii=0; ii<nDiscrete2; ii++)
543  {
544  if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
545  }
546  if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
547  if(ii<nIso)
548  {
549  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
550  // isotropic distribution
551  //G4double theta = pi*G4UniformRand();
552  G4double theta = std::acos(2.*G4UniformRand()-1.);
553  //But this is alos never cause real effect at least with G4NDL3.13 TK not repFlag == 2 AND isoFlag != 1
554  G4double phi = twopi*G4UniformRand();
555  G4double sinth = std::sin(theta);
556  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
557  //G4double en = theOne->GetTotalEnergy();
558  G4double en = theOne->GetTotalMomentum();
559  //But never cause real effect at least with G4NDL3.13 TK
560  G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
561  theOne->SetMomentum( tempVector ) ;
562  }
563  else if(tabulationType==1)
564  {
565  // legendre polynomials
566  G4int itt(0);
567  for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
568  {
569  itt = iii;
570  if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
571  break;
572  }
573  G4ParticleHPLegendreStore aStore(2);
574  aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt]));
575  //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
576  //TKDB 110512
577  if ( itt > 0 )
578  {
579  aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1]));
580  }
581  else
582  {
583  aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt]));
584  }
585  G4double cosTh = aStore.SampleMax(anEnergy);
586  G4double theta = std::acos(cosTh);
587  G4double phi = twopi*G4UniformRand();
588  G4double sinth = std::sin(theta);
589  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
590  //G4double en = theOne->GetTotalEnergy();
591  G4double en = theOne->GetTotalMomentum();
592  //But never cause real effect at least with G4NDL3.13 TK
593  G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
594  theOne->SetMomentum( tempVector ) ;
595  }
596  else
597  {
598  // tabulation of probabilities.
599  G4int itt(0);
600  for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
601  {
602  itt = iii;
603  if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
604  break;
605  }
606  G4double costh = theAngular[ii-nIso][itt].GetCosTh(); // no interpolation yet @@
607  G4double theta = std::acos(costh);
608  G4double phi = twopi*G4UniformRand();
609  G4double sinth = std::sin(theta);
610  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
611  //G4double en = theOne->GetTotalEnergy();
612  G4double en = theOne->GetTotalMomentum();
613  //But never cause real effect at least with G4NDL3.13 TK
614  G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
615  theOne->SetMomentum( tmpVector ) ;
616  }
617  }
618  thePhotons->push_back(theOne);
619  }
620  else if( repFlag==0 )
621  {
622 
623 // TK add
624  if ( thePartialXsec == 0 )
625  {
626  //G4cout << "repFlag is 0, but no PartialXsec data" << G4endl;
627  //G4cout << "This is not support yet." << G4endl;
628  return thePhotons;
629  }
630 
631 // Partial Case
632 
633  G4ReactionProduct * theOne = new G4ReactionProduct;
634  theOne->SetDefinition( G4Gamma::Gamma() );
635  thePhotons->push_back( theOne );
636 
637 // Energy
638 
639  //G4cout << "Partial Case nDiscrete " << nDiscrete << G4endl;
640  G4double sum = 0.0;
641  std::vector < G4double > dif( nDiscrete , 0.0 );
642  for ( G4int j = 0 ; j < nDiscrete ; j++ )
643  {
644  G4double x = thePartialXsec[ j ].GetXsec( anEnergy ); // x in barn
645  if ( x > 0 )
646  {
647  sum += x;
648  }
649  dif [ j ] = sum;
650  //G4cout << "j " << j << ", x " << x << ", dif " << dif [ j ] << G4endl;
651  }
652 
653  G4double rand = G4UniformRand();
654 
655  G4int iphoton = 0;
656  for ( G4int j = 0 ; j < nDiscrete ; j++ )
657  {
658  G4double y = rand*sum;
659  if ( dif [ j ] > y )
660  {
661  iphoton = j;
662  break;
663  }
664  }
665  //G4cout << "iphoton " << iphoton << G4endl;
666  //G4cout << "photon energy " << theGammas[ iphoton ] /eV << G4endl;
667 
668  // Statistically suppress the photon according to reaction cross section
669  // Fix proposed by Artem Zontikov, Bug report #1824
670  if (theReactionXsec) {
671  if (thePartialXsec[iphoton].GetXsec(anEnergy)/theReactionXsec->GetXsec(anEnergy) < G4UniformRand() ) {
672  delete thePhotons;
673  thePhotons = 0;
674  return thePhotons;
675  }
676  }
677 
678 // Angle
679  G4double cosTheta = 0.0; // mu
680 
681  if ( isoFlag == 1 )
682  {
683 
684 // Isotropic Case
685 
686  cosTheta = 2.*G4UniformRand()-1;
687 
688  }
689  else
690  {
691 
692  if ( iphoton < nIso )
693  {
694 
695 // still Isotropic
696 
697  cosTheta = 2.*G4UniformRand()-1;
698 
699  }
700  else
701  {
702 
703  //G4cout << "Not Isotropic and isoFlag " << isoFlag << G4endl;
704  //G4cout << "tabulationType " << tabulationType << G4endl;
705  //G4cout << "nDiscrete2 " << nDiscrete2 << G4endl;
706  //G4cout << "nIso " << nIso << G4endl;
707  //G4cout << "size of nNeu " << nDiscrete2-nIso << G4endl;
708  //G4cout << "nNeu[iphoton-nIso] " << nNeu[iphoton-nIso] << G4endl;
709 
710  if ( tabulationType == 1 )
711  {
712 // legendre polynomials
713 
714  G4int iangle = 0;
715  for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
716  {
717  iangle = j;
718  if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
719  }
720 
721  G4ParticleHPLegendreStore aStore( 2 );
722  aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );
723  aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) );
724 
725  cosTheta = aStore.SampleMax( anEnergy );
726 
727  }
728  else if ( tabulationType == 2 )
729  {
730 
731 // tabulation of probabilities.
732 
733  G4int iangle = 0;
734  for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
735  {
736  iangle = j;
737  if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
738  }
739 
740  cosTheta = theAngular[iphoton-nIso][ iangle ].GetCosTh(); // no interpolation yet @@
741 
742  }
743  }
744  }
745 
746 // Set
747  G4double phi = twopi*G4UniformRand();
748  G4double theta = std::acos( cosTheta );
749  G4double sinTheta = std::sin( theta );
750 
751  G4double photonE = theGammas[ iphoton ];
752  G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
753  G4ThreeVector photonP = photonE * direction;
754  thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
755 
756  }
757  else
758  {
759  delete thePhotons;
760  thePhotons = 0; // no gamma data available; some work needed @@@@@@@
761  }
762  return thePhotons;
763 }
764