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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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