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

Last change on this file since 1228 was 1055, checked in by garnier, 15 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.