G4NeutronHPPhotonDist Class Reference

#include <G4NeutronHPPhotonDist.hh>


Public Member Functions

 G4NeutronHPPhotonDist ()
 ~G4NeutronHPPhotonDist ()
G4bool InitMean (std::ifstream &aDataFile)
void InitAngular (std::ifstream &aDataFile)
void InitEnergies (std::ifstream &aDataFile)
void InitPartials (std::ifstream &aDataFile)
G4ReactionProductVectorGetPhotons (G4double anEnergy)
G4double GetTargetMass ()
G4bool NeedsCascade ()
G4double GetLevelEnergy ()


Detailed Description

Definition at line 55 of file G4NeutronHPPhotonDist.hh.


Constructor & Destructor Documentation

G4NeutronHPPhotonDist::G4NeutronHPPhotonDist (  )  [inline]

Definition at line 59 of file G4NeutronHPPhotonDist.hh.

00060       : repFlag( 0 ) 
00061       , targetMass( 0.0 ) 
00062       , nDiscrete( 0 ) 
00063       , isoFlag( 0 )
00064       , tabulationType( 0 )
00065       , nDiscrete2( 0 ) 
00066       , nIso( 0 ) 
00067       , nPartials( 0 )
00068       , theInternalConversionFlag( 0 )
00069       , nGammaEnergies( 0 )
00070       , theBaseEnergy( 0.0 )
00071   {
00072 
00073      disType = 0;
00074      energy = 0;
00075      theYield = 0;
00076      thePartialXsec = 0;
00077      isPrimary = 0;
00078      theShells = 0;
00079      theGammas = 0;
00080      nNeu = 0;
00081      theLegendre = 0;
00082      theAngular = 0;
00083      distribution = 0;
00084      probs = 0;
00085      partials = 0;
00086      actualMult = 0;
00087 
00088      theLevelEnergies = 0;
00089      theTransitionProbabilities = 0;
00090      thePhotonTransitionFraction = 0;
00091 
00092   }

G4NeutronHPPhotonDist::~G4NeutronHPPhotonDist (  )  [inline]

Definition at line 94 of file G4NeutronHPPhotonDist.hh.

00095   {
00096      delete [] disType;
00097      delete [] energy;
00098      delete [] theYield;
00099      delete [] thePartialXsec;
00100      delete [] isPrimary;
00101      delete [] theShells;
00102      delete [] theGammas;
00103      delete [] nNeu;
00104      delete [] theAngular;
00105      delete [] distribution;
00106      delete [] probs;
00107 
00108      if ( theLegendre != NULL )
00109      {
00110         for ( G4int i = 0 ; i < (nDiscrete2-nIso) ; i++ )
00111            if ( theLegendre[i] != NULL ) delete[] theLegendre[i]; 
00112 
00113         delete [] theLegendre;
00114      }
00115 
00116      if ( partials != 0 ) 
00117      {
00118         for ( G4int i = 0 ; i < nPartials ; i++ )
00119            { delete partials[i]; }
00120 
00121         delete [] partials;
00122      }
00123 
00124      delete [] actualMult;
00125 
00126      // delete theLevelEnergies;
00127      // delete theTransitionProbabilities;
00128      // delete thePhotonTransitionFraction;
00129 // TKDB
00130      delete [] theLevelEnergies;
00131      delete [] theTransitionProbabilities;
00132      delete [] thePhotonTransitionFraction;
00133   }


Member Function Documentation

G4double G4NeutronHPPhotonDist::GetLevelEnergy (  )  [inline]

Definition at line 149 of file G4NeutronHPPhotonDist.hh.

Referenced by G4NeutronHPInelasticCompFS::CompositeApply().

00149 {return theBaseEnergy;}

G4ReactionProductVector * G4NeutronHPPhotonDist::GetPhotons ( G4double  anEnergy  ) 

Definition at line 281 of file G4NeutronHPPhotonDist.cc.

References DBL_MAX, G4Electron::Electron(), G4cout, G4endl, G4Poisson(), G4UniformRand, G4Gamma::Gamma(), G4NeutronHPAngularP::GetCosTh(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4NeutronHPVector::GetVectorLength(), G4NeutronHPVector::GetX(), G4NeutronHPPartial::GetY(), G4INCL::Math::pi, G4NeutronHPVector::Sample(), G4NeutronHPLegendreStore::SampleMax(), G4NeutronHPLegendreStore::SetCoeff(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetMomentum(), and G4ReactionProduct::SetTotalEnergy().

Referenced by G4NeutronHPCaptureFS::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::CompositeApply(), and G4NeutronHPFSFissionFS::GetPhotons().

00282 {
00283 
00284   //G4cout << "G4NeutronHPPhotonDist::GetPhotons repFlag " << repFlag << G4endl;
00285   // the partial cross-section case is not in this yet. @@@@  << 070601 TK add partial 
00286   G4int i, ii, iii;
00287   G4int nSecondaries = 0;
00288   G4ReactionProductVector * thePhotons = new G4ReactionProductVector;
00289   if(repFlag==1)
00290   {
00291     G4double current=0;
00292     for(i=0; i<nDiscrete; i++)
00293     {
00294       current = theYield[i].GetY(anEnergy);
00295       actualMult[i] = G4Poisson(current); // max cut-off still missing @@@
00296       if(nDiscrete==1&&current<1.0001) 
00297       {
00298         actualMult[i] = static_cast<G4int>(current);
00299         if(current<1) 
00300         {
00301           actualMult[i] = 0;
00302           if(G4UniformRand()<current) actualMult[i] = 1;
00303         }
00304       }
00305       nSecondaries += actualMult[i];
00306     }
00307     //G4cout << "nSecondaries " << nSecondaries  << " anEnergy " << anEnergy/eV << G4endl;
00308     for(i=0;i<nSecondaries;i++)
00309     {
00310       G4ReactionProduct * theOne = new G4ReactionProduct;
00311       theOne->SetDefinition(G4Gamma::Gamma());
00312       thePhotons->push_back(theOne);
00313     }
00314     G4int count=0;
00315 
00316 /*
00317 G4double totalCascadeEnergy = 0.;
00318 G4double lastCascadeEnergy = 0.;
00319 G4double eGamm = 0;
00320 G4int maxEnergyIndex = 0;
00321 */
00322     //Gcout << "nDiscrete " << nDiscrete << " nPartials " << nPartials << G4endl;
00323 //3456
00324       if ( nDiscrete == 1 && nPartials == 1 )  
00325       {
00326          if ( actualMult[ 0 ] > 0 ) 
00327          {
00328             if ( disType[0] == 1 ) // continuum
00329             {
00330 
00331 /*
00332       for(ii=0; ii< actualMult[0]; ii++)
00333       {   
00334 
00335           G4double  sum=0, run=0;
00336           for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
00337           G4double random = G4UniformRand();
00338           G4int theP = 0;
00339           for(iii=0; iii<nPartials; iii++)
00340           {
00341             run+=probs[iii].GetY(anEnergy);
00342             theP = iii;
00343             if(random<run/sum) break;
00344           }
00345           if(theP==nPartials) theP=nPartials-1; // das sortiert J aus.
00346           sum=0; 
00347           G4NeutronHPVector * temp;
00348           temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
00349           // Looking for TotalCascdeEnergy or LastMaxEnergy
00350           if (ii == 0)
00351           {
00352             maxEnergyIndex = temp->GetVectorLength()-1;
00353             totalCascadeEnergy = temp->GetX(maxEnergyIndex);
00354             lastCascadeEnergy = totalCascadeEnergy;
00355           }
00356           lastCascadeEnergy -= eGamm;
00357           if (ii != actualMult[i]-1) eGamm = temp->SampleWithMax(lastCascadeEnergy);
00358           else eGamm = lastCascadeEnergy;
00359           thePhotons->operator[](count)->SetKineticEnergy(eGamm);
00360           delete temp;
00361 
00362      }
00363 */
00364                G4NeutronHPVector * temp;
00365                temp = partials[ 0 ]->GetY(anEnergy); //@@@ look at, seems fishy
00366                G4double maximumE = temp->GetX( temp->GetVectorLength()-1 ); // This is an assumption.
00367 
00368                //G4cout << "start " << actualMult[ 0 ] << " maximumE " << maximumE/eV << G4endl;
00369 
00370                std::vector< G4double > photons_e_best( actualMult[ 0 ] , 0.0 );
00371                G4double best = DBL_MAX;
00372                G4int maxTry = 1000; 
00373                for ( G4int j = 0 ; j < maxTry ; j++ )
00374                {
00375                   std::vector< G4double > photons_e( actualMult[ 0 ] , 0.0 );
00376                   for ( std::vector< G4double >::iterator 
00377                       it = photons_e.begin() ; it < photons_e.end() ; it++ ) 
00378                  {
00379                      *it = temp->Sample();   
00380                  }
00381                  if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) > maximumE ) 
00382                  {
00383                     if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) < best )
00384                        photons_e_best = photons_e;
00385                     continue;
00386                  }
00387                  else
00388                  {
00389                     for ( std::vector< G4double >::iterator 
00390                         it = photons_e.begin() ; it < photons_e.end() ; it++ ) 
00391                     {
00392                        thePhotons->operator[](count)->SetKineticEnergy( *it );
00393                     }  
00394                     //G4cout << "OK " << actualMult[0] << " j " << j << " total photons E  " 
00395                     //          << std::accumulate( photons_e.begin() , photons_e.end() , 0.0 )/eV << " ratio " << std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) / maximumE 
00396                     //          << G4endl;
00397                     
00398                     break;
00399                  }
00400                  G4cout << "NeutronHPPhotonDist could not find fitted energy set for multiplicity of " <<  actualMult[0] << "." << G4endl; 
00401                  G4cout << "NeutronHPPhotonDist will use the best set." << G4endl; 
00402                  for ( std::vector< G4double >::iterator 
00403                      it = photons_e_best.begin() ; it < photons_e_best.end() ; it++ ) 
00404                  {
00405                      thePhotons->operator[](count)->SetKineticEnergy( *it );
00406                  }
00407                  //G4cout << "Not Good " << actualMult[0] << " j " << j << " total photons E  " 
00408                  //       << best/eV << " ratio " << best / maximumE 
00409                  //       << G4endl;
00410                }
00411                // TKDB
00412                delete temp;
00413             }
00414             else    // discrete
00415             {
00416                thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
00417             }
00418             count++;
00419             if(count    > nSecondaries)  throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
00420         }
00421          
00422       }
00423       else
00424       {
00425     for(i=0; i<nDiscrete; i++)
00426     { 
00427       for(ii=0; ii< actualMult[i]; ii++)
00428       {   
00429         if(disType[i]==1) // continuum
00430         {
00431           G4double  sum=0, run=0;
00432           for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
00433           G4double random = G4UniformRand();
00434           G4int theP = 0;
00435           for(iii=0; iii<nPartials; iii++)
00436           {
00437             run+=probs[iii].GetY(anEnergy);
00438             theP = iii;
00439             if(random<run/sum) break;
00440           }
00441           if(theP==nPartials) theP=nPartials-1; // das sortiert J aus.
00442           sum=0; 
00443           G4NeutronHPVector * temp;
00444           temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
00445           G4double eGamm = temp->Sample();
00446           thePhotons->operator[](count)->SetKineticEnergy(eGamm);
00447           delete temp;
00448         }
00449         else // discrete
00450         {
00451           thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
00452         }
00453         count++;
00454         if(count > nSecondaries)  throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
00455       }
00456     }
00457       }
00458     // now do the angular distributions...
00459     if( isoFlag == 1)
00460     {
00461       for (i=0; i< nSecondaries; i++)
00462       {
00463         G4double costheta = 2.*G4UniformRand()-1;
00464         G4double theta = std::acos(costheta);
00465         G4double phi = twopi*G4UniformRand();
00466         G4double sinth = std::sin(theta);
00467         G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00468         G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00469         thePhotons->operator[](i)->SetMomentum( temp ) ;
00470   //      G4cout << "Isotropic distribution in PhotonDist"<<temp<<G4endl;
00471       }
00472     }
00473     else
00474     {
00475       for(i=0; i<nSecondaries; i++)
00476       { 
00477         G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
00478         for(ii=0; ii<nDiscrete2; ii++) 
00479         {
00480           if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
00481         }
00482         if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
00483         if(ii<nIso)
00484         {
00485           // isotropic distribution
00486           G4double theta = pi*G4UniformRand();
00487           G4double phi = twopi*G4UniformRand();
00488           G4double sinth = std::sin(theta);
00489           G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00490           G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00491           thePhotons->operator[](i)->SetMomentum( tempVector ) ;
00492         }
00493         else if(tabulationType==1)
00494         {
00495           // legendre polynomials
00496           G4int it(0);
00497           for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
00498           {
00499             it = iii;
00500             if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
00501               break;
00502           }
00503           G4NeutronHPLegendreStore aStore(2);
00504           aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));  
00505           //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 
00506           //TKDB 110512
00507           if ( it > 0 ) 
00508           {
00509              aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 
00510           }
00511           else
00512           {
00513              aStore.SetCoeff(0, &(theLegendre[ii-nIso][it])); 
00514           }
00515           G4double cosTh = aStore.SampleMax(anEnergy);
00516           G4double theta = std::acos(cosTh);
00517           G4double phi = twopi*G4UniformRand();
00518           G4double sinth = std::sin(theta);
00519           G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00520           G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00521           thePhotons->operator[](i)->SetMomentum( tempVector ) ;
00522         }
00523         else
00524         {
00525           // tabulation of probabilities.
00526           G4int it(0);
00527           for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
00528           {
00529             it = iii;
00530             if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
00531               break;
00532           }
00533           G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@
00534           G4double theta = std::acos(costh);
00535           G4double phi = twopi*G4UniformRand();
00536           G4double sinth = std::sin(theta);
00537           G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00538           G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
00539           thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
00540         }
00541       }  
00542     } 
00543   }
00544   else if(repFlag == 2)
00545   {
00546     G4double * running = new G4double[nGammaEnergies];
00547     running[0]=theTransitionProbabilities[0];
00548     //G4int i; //declaration at 284th
00549     for(i=1; i<nGammaEnergies; i++)
00550     {
00551       running[i]=running[i-1]+theTransitionProbabilities[i];
00552     }
00553     G4double random = G4UniformRand();
00554     G4int it=0;
00555     for(i=0; i<nGammaEnergies; i++)
00556     {
00557       it = i;
00558       if(random < running[i]/running[nGammaEnergies-1]) break;
00559     }
00560     delete [] running;
00561     G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
00562     G4ReactionProduct * theOne = new G4ReactionProduct;
00563     theOne->SetDefinition(G4Gamma::Gamma());
00564     random = G4UniformRand();
00565     if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
00566     {
00567       theOne->SetDefinition(G4Electron::Electron());
00568       //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00569       //But never enter at least with G4NDL3.13
00570       totalEnergy += G4Electron::Electron()->GetPDGMass(); //proposed correction: add this line for electron
00571     }
00572     theOne->SetTotalEnergy(totalEnergy);
00573     if( isoFlag == 1 )
00574     {
00575       G4double costheta = 2.*G4UniformRand()-1;
00576       G4double theta = std::acos(costheta);
00577       G4double phi = twopi*G4UniformRand();
00578       G4double sinth = std::sin(theta);
00579       //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00580       //G4double en = theOne->GetTotalEnergy();
00581       G4double en = theOne->GetTotalMomentum();
00582       //But never cause real effect at least with G4NDL3.13 TK
00583       G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00584       theOne->SetMomentum( temp ) ;
00585     }
00586     else
00587     {
00588       G4double currentEnergy = theOne->GetTotalEnergy();
00589       for(ii=0; ii<nDiscrete2; ii++) 
00590       {
00591         if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
00592       }
00593       if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
00594       if(ii<nIso)
00595       {
00596         //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00597         // isotropic distribution
00598         //G4double theta = pi*G4UniformRand();
00599         G4double theta = std::acos(2.*G4UniformRand()-1.);
00600         //But this is alos never cause real effect at least with G4NDL3.13 TK  not repFlag == 2 AND isoFlag != 1
00601         G4double phi = twopi*G4UniformRand();
00602         G4double sinth = std::sin(theta);
00603         //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00604         //G4double en = theOne->GetTotalEnergy();
00605         G4double en = theOne->GetTotalMomentum();
00606         //But never cause real effect at least with G4NDL3.13 TK
00607         G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00608         theOne->SetMomentum( tempVector ) ;
00609       }
00610       else if(tabulationType==1)
00611       {
00612         // legendre polynomials
00613         G4int itt(0);
00614         for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
00615         {
00616           itt = iii;
00617           if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
00618             break;
00619         }
00620         G4NeutronHPLegendreStore aStore(2);
00621         aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt]));  
00622         //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 
00623         //TKDB 110512
00624         if ( itt > 0 ) 
00625         {
00626            aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1])); 
00627         }
00628         else
00629         {
00630            aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt])); 
00631         }
00632         G4double cosTh = aStore.SampleMax(anEnergy);
00633         G4double theta = std::acos(cosTh);
00634         G4double phi = twopi*G4UniformRand();
00635         G4double sinth = std::sin(theta);
00636         //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00637         //G4double en = theOne->GetTotalEnergy();
00638         G4double en = theOne->GetTotalMomentum();
00639         //But never cause real effect at least with G4NDL3.13 TK
00640         G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00641         theOne->SetMomentum( tempVector ) ;
00642       }
00643       else
00644       {
00645         // tabulation of probabilities.
00646         G4int itt(0);
00647         for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
00648         {
00649           itt = iii;
00650           if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
00651             break;
00652         }
00653         G4double costh = theAngular[ii-nIso][itt].GetCosTh(); // no interpolation yet @@
00654         G4double theta = std::acos(costh);
00655         G4double phi = twopi*G4UniformRand();
00656         G4double sinth = std::sin(theta);
00657         //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
00658         //G4double en = theOne->GetTotalEnergy();
00659         G4double en = theOne->GetTotalMomentum();
00660         //But never cause real effect at least with G4NDL3.13 TK
00661         G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
00662         theOne->SetMomentum( tmpVector ) ;
00663       }
00664     }
00665     thePhotons->push_back(theOne);
00666   }
00667   else if( repFlag==0 )
00668   {
00669 
00670 // TK add 
00671       if ( thePartialXsec == 0 ) 
00672       {
00673          //G4cout << "repFlag is 0, but no PartialXsec data" << G4endl;
00674          //G4cout << "This is not support yet." << G4endl;
00675          return thePhotons;
00676       }
00677 
00678 // Partial Case 
00679 
00680       G4ReactionProduct * theOne = new G4ReactionProduct;
00681       theOne->SetDefinition( G4Gamma::Gamma() );
00682       thePhotons->push_back( theOne );
00683 
00684 // Energy 
00685 
00686       //G4cout << "Partial Case nDiscrete " << nDiscrete << G4endl;
00687       G4double sum = 0.0; 
00688       std::vector < G4double > dif( nDiscrete , 0.0 ); 
00689       for ( G4int j = 0 ; j < nDiscrete ; j++ ) 
00690       {
00691          G4double x = thePartialXsec[ j ].GetXsec( anEnergy );  // x in barn 
00692          if ( x > 0 ) 
00693          {
00694             sum += x;   
00695          } 
00696          dif [ j ] = sum; 
00697          //G4cout << "j " << j << ", x " << x << ", dif " << dif [ j ] << G4endl;
00698       }
00699       
00700       G4double rand = G4UniformRand();
00701 
00702       G4int iphoton = 0; 
00703       for ( G4int j = 0 ; j < nDiscrete ; j++ ) 
00704       {
00705          G4double y = rand*sum; 
00706          if ( dif [ j ] > y ) 
00707          {
00708             iphoton = j; 
00709             break;  
00710          } 
00711       }
00712       //G4cout << "iphoton " << iphoton << G4endl;
00713       //G4cout << "photon energy " << theGammas[ iphoton ] /eV  << G4endl;
00714 
00715 // Angle 
00716       G4double cosTheta = 0.0; // mu
00717 
00718       if ( isoFlag == 1 )
00719       {
00720 
00721 //       Isotropic Case
00722 
00723          cosTheta = 2.*G4UniformRand()-1;
00724 
00725       }
00726       else
00727       {
00728 
00729          if ( iphoton < nIso )
00730          {
00731 
00732 //          still Isotropic 
00733 
00734             cosTheta = 2.*G4UniformRand()-1;
00735 
00736          }
00737          else
00738          {
00739 
00740             //G4cout << "Not Isotropic and isoFlag " << isoFlag  << G4endl;
00741             //G4cout << "tabulationType " << tabulationType << G4endl;
00742             //G4cout << "nDiscrete2  " << nDiscrete2 << G4endl;
00743             //G4cout << "nIso " << nIso << G4endl;
00744             //G4cout << "size of nNeu " << nDiscrete2-nIso << G4endl;
00745             //G4cout << "nNeu[iphoton-nIso] " << nNeu[iphoton-nIso] << G4endl;
00746 
00747             if ( tabulationType == 1 )
00748             {
00749 //             legendre polynomials
00750 
00751                G4int iangle = 0; 
00752                for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
00753                {
00754                   iangle = j;
00755                   if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
00756                }
00757  
00758                G4NeutronHPLegendreStore aStore( 2 );
00759                aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );  
00760                aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) ); 
00761 
00762                cosTheta = aStore.SampleMax( anEnergy );
00763 
00764             }
00765             else if ( tabulationType == 2 )
00766             {
00767         
00768 //             tabulation of probabilities.
00769 
00770                G4int iangle = 0; 
00771                for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
00772                {
00773                   iangle = j;
00774                   if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
00775                }
00776 
00777                cosTheta = theAngular[iphoton-nIso][ iangle ].GetCosTh(); // no interpolation yet @@
00778 
00779             }
00780          }
00781       }
00782       
00783 // Set 
00784       G4double phi = twopi*G4UniformRand();
00785       G4double theta = std::acos( cosTheta );
00786       G4double sinTheta = std::sin( theta );
00787 
00788       G4double photonE = theGammas[ iphoton ];
00789       G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
00790       G4ThreeVector photonP = photonE * direction;
00791       thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
00792     
00793   }
00794   else
00795   {
00796     delete thePhotons;
00797     thePhotons = 0; // no gamma data available; some work needed @@@@@@@
00798   }    
00799   return thePhotons;
00800 }

G4double G4NeutronHPPhotonDist::GetTargetMass (  )  [inline]

Definition at line 145 of file G4NeutronHPPhotonDist.hh.

Referenced by G4NeutronHPCaptureFS::Init().

00145 {return targetMass;}

void G4NeutronHPPhotonDist::InitAngular ( std::ifstream &  aDataFile  ) 

Definition at line 121 of file G4NeutronHPPhotonDist.cc.

References G4cout, G4endl, and G4InterpolationManager::Init().

Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPCaptureFS::Init(), and G4FissionLibrary::Init().

00122 {
00123 
00124   G4int i, ii;
00125   //angular distributions
00126   aDataFile >> isoFlag;
00127   if (isoFlag != 1)
00128   {
00129 if ( repFlag == 2 ) G4cout << "G4NeutronHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 Hyper News. Thanks." << G4endl;
00130     aDataFile >> tabulationType >> nDiscrete2 >> nIso;
00131 //080731
00132       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;
00133 
00134       // The order of cross section (InitPartials) and distribution (InitAngular here) data are different, we have to re-coordinate consistent data order.
00135       std::vector < G4double > vct_gammas_par; 
00136       std::vector < G4double > vct_shells_par; 
00137       std::vector < G4int > vct_primary_par; 
00138       std::vector < G4int > vct_distype_par; 
00139       std::vector < G4NeutronHPVector* > vct_pXS_par;
00140       if ( theGammas != NULL ) 
00141       {
00142          //copy the cross section data 
00143          for ( i = 0 ; i < nDiscrete ; i++ )
00144          {
00145             vct_gammas_par.push_back( theGammas[ i ] );
00146             vct_shells_par.push_back( theShells[ i ] );
00147             vct_primary_par.push_back( isPrimary[ i ] );
00148             vct_distype_par.push_back( disType[ i ] );
00149             G4NeutronHPVector* hpv = new G4NeutronHPVector;
00150             *hpv = thePartialXsec[ i ];
00151             vct_pXS_par.push_back( hpv );
00152          }
00153       }
00154      if ( theGammas == NULL ) theGammas = new G4double[nDiscrete2];
00155      if ( theShells == NULL ) theShells = new G4double[nDiscrete2];
00156 //080731
00157 
00158     for (i=0; i< nIso; i++) // isotropic photons
00159     {
00160        aDataFile >> theGammas[i] >> theShells[i];
00161        theGammas[i]*=eV;
00162        theShells[i]*=eV;
00163     }
00164     nNeu = new G4int [nDiscrete2-nIso];
00165     if(tabulationType==1)theLegendre=new G4NeutronHPLegendreTable *[nDiscrete2-nIso];
00166     if(tabulationType==2)theAngular =new G4NeutronHPAngularP *[nDiscrete2-nIso];
00167     for(i=nIso; i< nDiscrete2; i++)
00168     {
00169       if(tabulationType==1) 
00170       {
00171         aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
00172         theGammas[i]*=eV;
00173         theShells[i]*=eV;
00174         theLegendre[i-nIso]=new G4NeutronHPLegendreTable[nNeu[i-nIso]];
00175         theLegendreManager.Init(aDataFile); 
00176         for (ii=0; ii<nNeu[i-nIso]; ii++)
00177         {
00178           theLegendre[i-nIso][ii].Init(aDataFile);
00179         }
00180       }
00181       else if(tabulationType==2)
00182       {
00183         aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
00184         theGammas[i]*=eV;
00185         theShells[i]*=eV;
00186         theAngular[i-nIso]=new G4NeutronHPAngularP[nNeu[i-nIso]];
00187         for (ii=0; ii<nNeu[i-nIso]; ii++)
00188         {
00189           theAngular[i-nIso][ii].Init(aDataFile);
00190         }
00191       }
00192       else
00193       {
00194         G4cout << "tabulation type: tabulationType"<<G4endl;
00195         throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions.");
00196       }
00197     }
00198 //080731
00199      if ( vct_gammas_par.size() > 0 ) 
00200      {
00201         //Reordering cross section data to corrsponding distribution data 
00202         for ( i = 0 ; i < nDiscrete ; i++ ) 
00203         {
00204            for ( G4int j = 0 ; j < nDiscrete ; j++ )  
00205            {
00206               // Checking gamma and shell to identification 
00207               if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] )
00208               {
00209                  isPrimary [ i ] = vct_primary_par [ j ];
00210                  disType [ i ] = vct_distype_par [ j ];
00211                  thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) );
00212               }
00213            }
00214         }
00215         //Garbage collection 
00216         for ( std::vector < G4NeutronHPVector* >::iterator 
00217            it = vct_pXS_par.begin() ; it != vct_pXS_par.end() ; it++ )
00218         {
00219            delete *it;
00220         }
00221      }
00222 //080731
00223   }
00224 }

void G4NeutronHPPhotonDist::InitEnergies ( std::ifstream &  aDataFile  ) 

Definition at line 227 of file G4NeutronHPPhotonDist.cc.

References G4NeutronHPPartial::Init(), G4NeutronHPVector::Init(), and G4NeutronHPPartial::InitInterpolation().

Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPCaptureFS::Init(), and G4FissionLibrary::Init().

00228 {
00229   G4int i, energyDistributionsNeeded = 0;
00230   for (i=0; i<nDiscrete; i++)
00231   {
00232     if( disType[i]==1) energyDistributionsNeeded =1;
00233   }
00234   if(!energyDistributionsNeeded) return;
00235   aDataFile >>  nPartials;
00236   distribution = new G4int[nPartials];
00237   probs = new G4NeutronHPVector[nPartials];
00238   partials = new G4NeutronHPPartial * [nPartials];
00239   G4int nen;
00240   G4int dummy;
00241   for (i=0; i<nPartials; i++)
00242   {
00243     aDataFile >> dummy;
00244     probs[i].Init(aDataFile, eV);
00245     aDataFile >> nen;
00246     partials[i] = new G4NeutronHPPartial(nen);
00247     partials[i]->InitInterpolation(aDataFile);
00248     partials[i]->Init(aDataFile);
00249   }
00250 }

G4bool G4NeutronHPPhotonDist::InitMean ( std::ifstream &  aDataFile  ) 

Definition at line 54 of file G4NeutronHPPhotonDist.cc.

References G4cout, G4endl, and G4NeutronHPVector::Init().

Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPCaptureFS::Init(), and G4FissionLibrary::Init().

00055 {
00056 
00057   G4bool result = true;
00058   if(aDataFile >> repFlag)
00059   {
00060 
00061     aDataFile >> targetMass;
00062     if(repFlag==1)
00063     {
00064     // multiplicities
00065       aDataFile >> nDiscrete;
00066       disType = new G4int[nDiscrete];
00067       energy = new G4double[nDiscrete];
00068       actualMult = new G4int[nDiscrete];
00069       theYield = new G4NeutronHPVector[nDiscrete];
00070       for (G4int i=0; i<nDiscrete; i++)
00071       {
00072         aDataFile >> disType[i]>>energy[i];
00073         energy[i]*=eV;
00074         theYield[i].Init(aDataFile, eV);
00075       }
00076     }
00077     else if(repFlag == 2)
00078     {
00079        aDataFile >> theInternalConversionFlag;
00080        aDataFile >> theBaseEnergy;
00081        theBaseEnergy*=eV;
00082        aDataFile >> theInternalConversionFlag;
00083        // theInternalConversionFlag == 1 No IC, theInternalConversionFlag == 2 with IC 
00084        aDataFile >> nGammaEnergies;
00085        theLevelEnergies = new G4double[nGammaEnergies];
00086        theTransitionProbabilities = new G4double[nGammaEnergies];
00087        if(theInternalConversionFlag == 2) thePhotonTransitionFraction = new G4double[nGammaEnergies];
00088        for(G4int  ii=0; ii<nGammaEnergies; ii++)
00089        {
00090          if(theInternalConversionFlag == 1)
00091          {
00092            aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii];
00093            theLevelEnergies[ii]*=eV;
00094          }
00095          else if(theInternalConversionFlag == 2)
00096          {
00097            aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii];
00098            theLevelEnergies[ii]*=eV;
00099          }
00100          else
00101          {
00102            throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Unknown conversion flag");
00103          }
00104       }
00105        // Note, that this is equivalent to using the 'Gamma' classes.
00106       // throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Transition probability array not sampled for the moment.");
00107     }
00108     else
00109     {
00110       G4cout << "Data representation in G4NeutronHPPhotonDist: "<<repFlag<<G4endl;
00111       throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: This data representation is not implemented.");
00112     }
00113   }
00114   else
00115   {
00116     result = false;
00117   }
00118   return result;
00119 }

void G4NeutronHPPhotonDist::InitPartials ( std::ifstream &  aDataFile  ) 

Definition at line 252 of file G4NeutronHPPhotonDist.cc.

References G4NeutronHPVector::Init().

Referenced by G4NeutronHPInelasticCompFS::Init(), and G4NeutronHPInelasticBaseFS::Init().

00253 {
00254 
00255   //G4cout << "G4NeutronHPPhotonDist::InitPartials " << G4endl;
00256   aDataFile >> nDiscrete >> targetMass;
00257   if(nDiscrete != 1)
00258   {
00259     theTotalXsec.Init(aDataFile, eV);
00260   }
00261   G4int i;
00262   theGammas = new G4double[nDiscrete];
00263   theShells = new G4double[nDiscrete];
00264   isPrimary = new G4int[nDiscrete];
00265   disType = new G4int[nDiscrete];
00266   thePartialXsec = new G4NeutronHPVector[nDiscrete];
00267   for(i=0; i<nDiscrete; i++)
00268   {
00269     aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
00270     theGammas[i]*=eV;
00271     theShells[i]*=eV;
00272     thePartialXsec[i].Init(aDataFile, eV);
00273   }  
00274 
00275   //G4cout << "G4NeutronHPPhotonDist::InitPartials Test " << G4endl;
00276   //G4cout << "G4NeutronHPPhotonDist::InitPartials nDiscrete " << nDiscrete << G4endl;
00277   //G4NeutronHPVector* aHP = new G4NeutronHPVector;
00278   //aHP->Check(1);
00279 }

G4bool G4NeutronHPPhotonDist::NeedsCascade (  )  [inline]

Definition at line 147 of file G4NeutronHPPhotonDist.hh.

00147 {return repFlag==2;}


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:42 2013 for Geant4 by  doxygen 1.4.7