source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPPhotonDist.cc@ 1142

Last change on this file since 1142 was 1055, checked in by garnier, 17 years ago

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