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

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

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