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

Last change on this file since 893 was 819, checked in by garnier, 17 years ago

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