Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/particle_hp/src/G4ParticleHPChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/particle_hp/src/G4ParticleHPChannel.cc (Version 11.3.0) and /processes/hadronic/models/particle_hp/src/G4ParticleHPChannel.cc (Version 10.7)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // neutron_hp -- source file                       26 // neutron_hp -- source file
 27 // J.P. Wellisch, Nov-1996                         27 // J.P. Wellisch, Nov-1996
 28 // A prototype of the low energy neutron trans     28 // A prototype of the low energy neutron transport model.
 29 //                                                 29 //
 30 // 070523 bug fix for G4FPE_DEBUG on by A. How     30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
 31 // 071031 bug fix T. Koi on behalf of A. Howar <<  31 // 071031 bug fix T. Koi on behalf of A. Howard 
 32 // 081203 bug fix in Register method by T. Koi     32 // 081203 bug fix in Register method by T. Koi
 33 //                                                 33 //
 34 // P. Arce, June-2014 Conversion neutron_hp to     34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
 35 //                                                 35 //
 36 // June-2019 - E. Mendoza --> Modification to      36 // June-2019 - E. Mendoza --> Modification to allow using an incomplete
 37 //   data library if the G4NEUTRONHP_SKIP_MISS     37 //   data library if the G4NEUTRONHP_SKIP_MISSING_ISOTOPES environmental
 38 //   flag is defined. The missing XS are set t     38 //   flag is defined. The missing XS are set to 0.
 39                                                    39 
 40 #include "G4ParticleHPChannel.hh"              << 
 41                                                << 
 42 #include "G4HadTmpUtil.hh"                     << 
 43 #include "G4ParticleHPElasticFS.hh"            << 
 44 #include "G4ParticleHPFinalState.hh"           << 
 45 #include "G4ParticleHPReactionWhiteBoard.hh"   << 
 46 #include "G4ParticleHPThermalBoost.hh"         << 
 47 #include "G4SystemOfUnits.hh"                  << 
 48 #include "globals.hh"                          << 
 49                                                    40 
 50 #include <cstdlib>                             <<  41 #include <stdlib.h>
 51                                                    42 
 52 G4ParticleHPChannel::G4ParticleHPChannel(G4Par <<  43 #include "G4ParticleHPChannel.hh"
 53 {                                              <<  44 #include "globals.hh"
 54   fManager = G4ParticleHPManager::GetInstance( <<  45 #include "G4SystemOfUnits.hh"
 55   if (fManager->GetUseWendtFissionModel()) {   <<  46 #include "G4ParticleHPFinalState.hh"
 56     wendtFissionGenerator = G4WendtFissionFrag <<  47 #include "G4HadTmpUtil.hh"
 57     // Make sure both fission fragment models  << 
 58     fManager->SetProduceFissionFragments(false << 
 59   }                                            << 
 60   theProjectile = (nullptr == p) ? G4Neutron:: << 
 61   theChannelData = new G4ParticleHPVector;     << 
 62 }                                              << 
 63                                                    48 
 64 G4ParticleHPChannel::~G4ParticleHPChannel()    <<  49 #include "G4ParticleHPReactionWhiteBoard.hh"
 65 {                                              << 
 66   delete theChannelData;                       << 
 67   // Following statement disabled to avoid SEG << 
 68   // theBuffer is also deleted as "theChannelD << 
 69   delete[] theIsotopeWiseData;                 << 
 70   if (theFinalStates != nullptr) {             << 
 71     for (G4int i = 0; i < niso; i++) {         << 
 72       delete theFinalStates[i];                << 
 73     }                                          << 
 74     delete[] theFinalStates;                   << 
 75   }                                            << 
 76   delete[] active;                             << 
 77 }                                              << 
 78                                                    50 
 79 G4double G4ParticleHPChannel::GetXsec(G4double <<  51 G4double G4ParticleHPChannel::GetXsec(G4double energy)
 80 {                                                  52 {
 81   return std::max(0., theChannelData->GetXsec(     53   return std::max(0., theChannelData->GetXsec(energy));
 82 }                                                  54 }
 83                                                <<  55   
 84 G4double G4ParticleHPChannel::GetWeightedXsec( <<  56 G4double G4ParticleHPChannel::GetWeightedXsec(G4double energy, G4int isoNumber)
 85                 G4int isoNumber) const         << 
 86 {                                                  57 {
 87   return theIsotopeWiseData[isoNumber].GetXsec     58   return theIsotopeWiseData[isoNumber].GetXsec(energy);
 88 }                                                  59 }
 89                                                <<  60   
 90 G4double G4ParticleHPChannel::GetFSCrossSectio <<  61 G4double G4ParticleHPChannel::GetFSCrossSection(G4double energy, G4int isoNumber)
 91             G4int isoNumber) const             << 
 92 {                                                  62 {
 93   return theFinalStates[isoNumber]->GetXsec(en     63   return theFinalStates[isoNumber]->GetXsec(energy);
 94 }                                                  64 }
 95                                                <<  65   
 96 void G4ParticleHPChannel::Init(G4Element* anEl <<  66   void G4ParticleHPChannel::
 97              const G4String& dirName, const G4 <<  67   Init(G4Element * anElement, const G4String dirName, const G4String aFSType) 
 98 {                                              <<  68   {
 99   theFSType = aFSType;                         <<  69     theFSType = aFSType;
100   Init(anElement, dirName);                    <<  70     Init(anElement, dirName);
101 }                                              << 
102                                                << 
103 void G4ParticleHPChannel::Init(G4Element* anEl << 
104 {                                              << 
105   theDir = dirName;                            << 
106   theElement = anElement;                      << 
107 }                                              << 
108                                                << 
109 G4bool G4ParticleHPChannel::Register(G4Particl << 
110 {                                              << 
111   ++registerCount;                             << 
112   G4int Z = theElement->GetZasInt();           << 
113                                                << 
114   niso = (G4int)theElement->GetNumberOfIsotope << 
115   const std::size_t nsize = niso > 0 ? niso :  << 
116                                                << 
117   delete[] theIsotopeWiseData;                 << 
118   theIsotopeWiseData = new G4ParticleHPIsoData << 
119   delete[] active;                             << 
120   active = new G4bool[nsize];                  << 
121                                                << 
122   delete[] theFinalStates;                     << 
123   theFinalStates = new G4ParticleHPFinalState* << 
124   delete theChannelData;                       << 
125   theChannelData = new G4ParticleHPVector;     << 
126   for (G4int i = 0; i < niso; ++i) {           << 
127     theFinalStates[i] = theFS->New();          << 
128     theFinalStates[i]->SetProjectile(theProjec << 
129   }                                                71   }
130   if (niso != 0 && registerCount == 0) {       <<  72    
131     for (G4int i1 = 0; i1 < niso; ++i1) {      <<  73   void G4ParticleHPChannel::Init(G4Element * anElement, const G4String dirName)  
132       G4int A = theElement->GetIsotope(i1)->Ge <<  74   {
133       G4int M = theElement->GetIsotope(i1)->Ge <<  75     theDir = dirName;
134       //G4cout <<" Init: normal case i=" << i1 <<  76     theElement = anElement;
135       //     << " Z=" << Z << " A=" << A << G4 << 
136       G4double frac = theElement->GetRelativeA << 
137       theFinalStates[i1]->SetA_Z(A, Z, M);     << 
138       UpdateData(A, Z, M, i1, frac, theProject << 
139     }                                          << 
140   }                                                77   }
141   G4bool result = HasDataInAnyFinalState();    <<  78   
                                                   >>  79   G4bool G4ParticleHPChannel::Register(G4ParticleHPFinalState *theFS)
                                                   >>  80   {
                                                   >>  81   registerCount++;
                                                   >>  82     G4int Z = G4lrint(theElement->GetZ());
142                                                    83 
143   // To avoid issuing hash by worker threads   <<  84     Z = Z-registerCount;
144   if (result) theChannelData->Hash();          <<  85     if ( registerCount > 5 ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material"); // for Elastic, Capture, Fission case 
                                                   >>  86     if ( Z < 1 ) return false; 
                                                   >>  87 /*
                                                   >>  88     if(registerCount<5)
                                                   >>  89     {
                                                   >>  90       Z = Z-registerCount;
                                                   >>  91     }
                                                   >>  92 */
                                                   >>  93     //if(Z=theElement->GetZ()-5) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material");
                                                   >>  94     // Bug fix by TK on behalf of AH
                                                   >>  95     //if ( Z <=theElement->GetZ()-5 ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material");
                                                   >>  96     G4int count = 0;
                                                   >>  97     if(registerCount==0) count = theElement->GetNumberOfIsotopes();
                                                   >>  98     if(count == 0||registerCount!=0) count +=
                                                   >>  99          theStableOnes.GetNumberOfIsotopes(Z);
                                                   >> 100     niso = count;
                                                   >> 101     delete [] theIsotopeWiseData;
                                                   >> 102     theIsotopeWiseData = new G4ParticleHPIsoData [niso];
                                                   >> 103     delete [] active;
                                                   >> 104     active = new G4bool[niso];
                                                   >> 105 
                                                   >> 106     delete [] theFinalStates;
                                                   >> 107     theFinalStates = new G4ParticleHPFinalState * [niso];
                                                   >> 108     delete theChannelData;
                                                   >> 109     theChannelData = new G4ParticleHPVector; 
                                                   >> 110     for(G4int i=0; i<niso; i++)
                                                   >> 111     {
                                                   >> 112       theFinalStates[i] = theFS->New();
                                                   >> 113       theFinalStates[i]->SetProjectile(theProjectile);
                                                   >> 114     }
                                                   >> 115     count = 0;
                                                   >> 116     G4int nIsos = niso;
                                                   >> 117     if(theElement->GetNumberOfIsotopes()!=0&&registerCount==0)
                                                   >> 118     {
                                                   >> 119       for (G4int i1=0; i1<nIsos; i1++)
                                                   >> 120       {
                                                   >> 121         // G4cout <<" Init: normal case"<<G4endl;
                                                   >> 122         G4int A = theElement->GetIsotope(i1)->GetN();
                                                   >> 123         G4int M = theElement->GetIsotope(i1)->Getm();
                                                   >> 124         G4double frac = theElement->GetRelativeAbundanceVector()[i1]/perCent;
                                                   >> 125         //theFinalStates[i1]->SetA_Z(A, Z);
                                                   >> 126   //UpdateData(A, Z, count++, frac);
                                                   >> 127         theFinalStates[i1]->SetA_Z(A, Z, M);
                                                   >> 128   UpdateData(A, Z, M, count++, frac, theProjectile);
                                                   >> 129       }
                                                   >> 130     } else {
                                                   >> 131       //G4cout <<" Init: mean case: "
                                                   >> 132       //       <<theStableOnes.GetNumberOfIsotopes(Z)<<" "
                                                   >> 133   //     <<Z<<" "<<theElement
                                                   >> 134   //     << G4endl;
                                                   >> 135       G4int first = theStableOnes.GetFirstIsotope(Z);
                                                   >> 136       for(G4int i1=0; 
                                                   >> 137         i1<theStableOnes.GetNumberOfIsotopes(Z);
                                                   >> 138         i1++)
                                                   >> 139       {
                                                   >> 140         G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
                                                   >> 141         G4double frac = theStableOnes.GetAbundance(first+i1);
                                                   >> 142         theFinalStates[i1]->SetA_Z(A, Z);
                                                   >> 143         UpdateData(A, Z, count++, frac, theProjectile);
                                                   >> 144       }
                                                   >> 145     }
                                                   >> 146     G4bool result = HasDataInAnyFinalState();
145                                                   147 
146   return result;                               << 148     //To avoid issuing hash by worker threads
147 }                                              << 149     if ( result ) theChannelData->Hash();
148                                                   150 
149 void G4ParticleHPChannel::UpdateData(G4int A,  << 151     return result;
150                                      G4double  << 
151                                      G4Particl << 
152 {                                              << 
153   // Initialze the G4FissionFragment generator << 
154   if (wendtFissionGenerator != nullptr) {      << 
155     wendtFissionGenerator->InitializeANucleus( << 
156   }                                               152   }
157                                                << 153   
158   theFinalStates[index]->Init(A, Z, M, theDir, << 154   //void G4ParticleHPChannel::UpdateData(G4int A, G4int Z, G4int index, G4double abundance)
159   if (!theFinalStates[index]->HasAnyData()) re << 155   void G4ParticleHPChannel::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance, G4ParticleDefinition* projectile)
160   // nothing there for exactly this isotope.   << 
161                                                << 
162   // the above has put the X-sec into the FS   << 
163   theBuffer = nullptr;                         << 
164   if (theFinalStates[index]->HasXsec()) {      << 
165     theBuffer = theFinalStates[index]->GetXsec << 
166     theBuffer->Times(abundance / 100.);        << 
167     theIsotopeWiseData[index].FillChannelData( << 
168   }                                            << 
169   else  // get data from CrossSection director << 
170   {                                               156   {
171     const G4String& tString = "/CrossSection"; << 157     // Initialze the G4FissionFragment generator for this isomer if needed
172     active[index] = theIsotopeWiseData[index]. << 158   if(wendtFissionGenerator)
173                                                << 159   {
174     if (active[index]) theBuffer = theIsotopeW << 160     wendtFissionGenerator->InitializeANucleus(A, Z, M, theDir);
                                                   >> 161   }
                                                   >> 162 
                                                   >> 163   theFinalStates[index]->Init(A, Z, M, theDir, theFSType, projectile);
                                                   >> 164     if(!theFinalStates[index]->HasAnyData()) return; // nothing there for exactly this isotope.
                                                   >> 165 
                                                   >> 166     // the above has put the X-sec into the FS
                                                   >> 167     theBuffer = 0;
                                                   >> 168     if(theFinalStates[index]->HasXsec())
                                                   >> 169     {
                                                   >> 170       theBuffer = theFinalStates[index]->GetXsec();
                                                   >> 171       theBuffer->Times(abundance/100.);
                                                   >> 172       theIsotopeWiseData[index].FillChannelData(theBuffer);
                                                   >> 173     }
                                                   >> 174     else // get data from CrossSection directory
                                                   >> 175     {
                                                   >> 176       G4String tString = "/CrossSection";
                                                   >> 177       //active[index] = theIsotopeWiseData[index].Init(A, Z, abundance, theDir, tString);
                                                   >> 178       active[index] = theIsotopeWiseData[index].Init(A, Z, M, abundance, theDir, tString);
                                                   >> 179       if(active[index]) theBuffer = theIsotopeWiseData[index].MakeChannelData();
                                                   >> 180     }
                                                   >> 181     if(theBuffer != 0) Harmonise(theChannelData, theBuffer);
175   }                                               182   }
176   if (theBuffer != nullptr) Harmonise(theChann << 183   
177 }                                              << 184   void G4ParticleHPChannel::Harmonise(G4ParticleHPVector *& theStore, G4ParticleHPVector * theNew)
178                                                << 
179 void G4ParticleHPChannel::Harmonise(G4Particle << 
180                                     G4Particle << 
181 {                                              << 
182   G4int s_tmp = 0, n = 0, m_tmp = 0;           << 
183   auto theMerge = new G4ParticleHPVector;      << 
184   G4ParticleHPVector* anActive = theStore;     << 
185   G4ParticleHPVector* aPassive = theNew;       << 
186   G4ParticleHPVector* tmp;                     << 
187   G4int a = s_tmp, p = n, t;                   << 
188   while (a < anActive->GetVectorLength() && p  << 
189     // Loop checking, 11.05.2015, T. Koi       << 
190   {                                               185   {
191     if (anActive->GetEnergy(a) <= aPassive->Ge << 186     G4int s_tmp = 0, n=0, m_tmp=0;
192       G4double xa = anActive->GetEnergy(a);    << 187     G4ParticleHPVector * theMerge = new G4ParticleHPVector;
193       theMerge->SetData(m_tmp, xa, anActive->G << 188     G4ParticleHPVector * anActive = theStore;
194       m_tmp++;                                 << 189     G4ParticleHPVector * aPassive = theNew;
195       a++;                                     << 190     G4ParticleHPVector * tmp;
196       G4double xp = aPassive->GetEnergy(p);    << 191     G4int a = s_tmp, p = n, t;
197       if (std::abs(std::abs(xp - xa) / xa) < 0 << 192     while (a<anActive->GetVectorLength()&&p<aPassive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
198         ++p;                                   << 193     {
                                                   >> 194       if(anActive->GetEnergy(a) <= aPassive->GetEnergy(p))
                                                   >> 195       {
                                                   >> 196         G4double xa  = anActive->GetEnergy(a);
                                                   >> 197         theMerge->SetData(m_tmp, xa, anActive->GetXsec(a)+std::max(0., aPassive->GetXsec(xa)) );
                                                   >> 198         m_tmp++;
                                                   >> 199         a++;
                                                   >> 200         G4double xp = aPassive->GetEnergy(p);
                                                   >> 201         if( std::abs(std::abs(xp-xa)/xa)<0.001 )
                                                   >> 202         {
                                                   >> 203           p++;
                                                   >> 204         }
                                                   >> 205       } else {
                                                   >> 206         tmp = anActive; t=a;
                                                   >> 207         anActive = aPassive; a=p;
                                                   >> 208         aPassive = tmp; p=t;
199       }                                           209       }
200     }                                             210     }
201     else {                                     << 211     while (a!=anActive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
202       tmp = anActive;                          << 212     {
203       t = a;                                   << 213       theMerge->SetData(m_tmp++, anActive->GetEnergy(a), anActive->GetXsec(a));
204       anActive = aPassive;                     << 214       a++;
205       a = p;                                   << 215     }
206       aPassive = tmp;                          << 216     while (p!=aPassive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
207       p = t;                                   << 217     {
                                                   >> 218       if(std::abs(theMerge->GetEnergy(std::max(0,m_tmp-1))-aPassive->GetEnergy(p))/aPassive->GetEnergy(p)>0.001)
                                                   >> 219         theMerge->SetData(m_tmp++, aPassive->GetEnergy(p), aPassive->GetXsec(p));
                                                   >> 220       p++;
208     }                                             221     }
                                                   >> 222     delete theStore;
                                                   >> 223     theStore = theMerge;
209   }                                               224   }
210   while (a != anActive->GetVectorLength())  // << 
211   {                                            << 
212     theMerge->SetData(m_tmp++, anActive->GetEn << 
213     ++a;                                       << 
214   }                                            << 
215   while (p != aPassive->GetVectorLength())  // << 
216   {                                            << 
217     if (std::abs(theMerge->GetEnergy(std::max( << 
218      aPassive->GetEnergy(p)) / aPassive->GetEn << 
219       theMerge->SetData(m_tmp++, aPassive->Get << 
220     ++p;                                       << 
221   }                                            << 
222   delete theStore;                             << 
223   theStore = theMerge;                         << 
224 }                                              << 
225                                                   225 
226 G4WendtFissionFragmentGenerator* G4ParticleHPC << 226 #include "G4ParticleHPThermalBoost.hh"
227   if ( wendtFissionGenerator ) return wendtFis << 
228   else                         return nullptr; << 
229 }                                              << 
230                                                   227 
231 G4HadFinalState*                               << 228   G4HadFinalState * G4ParticleHPChannel::
232 G4ParticleHPChannel::ApplyYourself(const G4Had << 229   ApplyYourself(const G4HadProjectile & theTrack, G4int anIsotope)
233            G4int anIsotope, G4bool isElastic)  << 230   {
234 {                                              << 231 //    G4cout << "G4ParticleHPChannel::ApplyYourself+"<<niso<<G4endl;
235   //G4cout << "G4ParticleHPChannel::ApplyYours << 232     if ( anIsotope != -1 && anIsotope != -2 ) 
236   //   << " ni=" << anIsotope << " isElastic=" << 233     {
237   if (anIsotope != -1 && anIsotope != -2) {    << 234        //Inelastic Case
238     // Inelastic Case                          << 235        //G4cout << "G4ParticleHPChannel Inelastic Case" 
239     //G4cout << "G4ParticleHPChannel Inelastic << 236        //<< " Z= " << this->GetZ(it) << " A = " << this->GetN(it) << G4endl;
240     //<< " Z= " << GetZ(anIsotope) << " A = "  << 237        G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA( (G4int)this->GetN(anIsotope) ); 
241     fManager->GetReactionWhiteBoard()->SetTarg << 238        G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ( (G4int)this->GetZ(anIsotope) ); 
242     fManager->GetReactionWhiteBoard()->SetTarg << 239        return theFinalStates[anIsotope]->ApplyYourself(theTrack);
243     return theFinalStates[anIsotope]->ApplyYou << 
244   }                                            << 
245   G4double sum = 0;                            << 
246   G4int it = 0;                                << 
247   auto xsec = new G4double[niso];              << 
248   G4ParticleHPThermalBoost aThermalE;          << 
249   for (G4int i = 0; i < niso; i++) {           << 
250     if (theFinalStates[i]->HasAnyData()) {     << 
251       /*                                       << 
252       G4cout << "FS: " << i << theTrack.GetDef << 
253        << " Z=" << theFinalStates[i]->GetZ()   << 
254        << " A=" << theFinalStates[i]->GetN()   << 
255        << G4endl;                              << 
256       */                                       << 
257       xsec[i] = theIsotopeWiseData[i].GetXsec( << 
258         aThermalE.GetThermalEnergy(theTrack, t << 
259                                    theFinalSta << 
260                                    theTrack.Ge << 
261       sum += xsec[i];                          << 
262     }                                             240     }
263     else {                                     << 241     G4double sum=0;
264       xsec[i] = 0;                             << 242     G4int it=0;
                                                   >> 243     G4double * xsec = new G4double[niso];
                                                   >> 244     G4ParticleHPThermalBoost aThermalE;
                                                   >> 245     for (G4int i=0; i<niso; i++)
                                                   >> 246     {
                                                   >> 247       if(theFinalStates[i]->HasAnyData())
                                                   >> 248       {
                                                   >> 249         xsec[i] = theIsotopeWiseData[i].GetXsec(aThermalE.GetThermalEnergy(theTrack,
                                                   >> 250                                                                theFinalStates[i]->GetN(),
                                                   >> 251                      theFinalStates[i]->GetZ(),
                                                   >> 252                              theTrack.GetMaterial()->GetTemperature()));
                                                   >> 253         sum += xsec[i];
                                                   >> 254       }
                                                   >> 255       else
                                                   >> 256       {
                                                   >> 257         xsec[i]=0;
                                                   >> 258       }
                                                   >> 259     } 
                                                   >> 260     if(sum == 0) 
                                                   >> 261     {
                                                   >> 262 //      G4cout << "G4ParticleHPChannel::ApplyYourself theFinalState->Initialize+"<<G4endl;
                                                   >> 263 //      G4cout << "G4ParticleHPChannel::ApplyYourself theFinalState->Initialize-"<<G4endl;
                                                   >> 264       it = static_cast<G4int>(niso*G4UniformRand());
265     }                                             265     }
266   }                                            << 266     else
267   if (sum == 0) {                              << 267     {
268     it = G4lrint(niso * G4UniformRand());      << 268 //      G4cout << "Are we still here? "<<sum<<G4endl;
269   }                                            << 269 //      G4cout << "TESTHP 23 NISO="<<niso<<G4endl;
270   else {                                       << 270       G4double random = G4UniformRand();
271     G4double random = G4UniformRand();         << 271       G4double running=0;
272     G4double running = 0;                      << 272 //      G4cout << "G4ParticleHPChannel::ApplyYourself Done the sum"<<niso<<G4endl;
273     for (G4int ix = 0; ix < niso; ix++) {      << 273 //      G4cout << "TESTHP 24 NISO="<<niso<<G4endl;
274       running += xsec[ix];                     << 274       for (G4int ix=0; ix<niso; ix++)
275       if (sum == 0 || random <= running / sum) << 275       {
276         it = ix;                               << 276         running += xsec[ix];
277         break;                                 << 277         //if(random<=running/sum) 
                                                   >> 278         if( sum == 0 || random <= running/sum ) 
                                                   >> 279         { 
                                                   >> 280           it = ix;
                                                   >> 281     break;
                                                   >> 282         }
278       }                                           283       }
                                                   >> 284       if(it==niso) it--;
279     }                                             285     }
280     if (it == niso) it--;                      << 286     delete [] xsec;
281   }                                            << 287     G4HadFinalState * theFinalState=0;
282   delete[] xsec;                               << 288     const G4int A = (G4int)this->GetN(it);
283   G4HadFinalState* theFinalState = nullptr;    << 289     const G4int Z = (G4int)this->GetZ(it);
284   const auto A = (G4int)this->GetN(it);        << 290     const G4int M = (G4int)this->GetM(it);
285   const auto Z = (G4int)this->GetZ(it);        << 
286   const auto M = (G4int)this->GetM(it);        << 
287                                                << 
288   //-2:Marker for Fission                      << 
289   if ((wendtFissionGenerator != nullptr) && an << 
290     theFinalState = wendtFissionGenerator->App << 
291   }                                            << 
292                                                   291 
293   // Use the standard procedure if the G4Fissi << 292                                        //-2:Marker for Fission
294   if (theFinalState == nullptr) {              << 293     if(wendtFissionGenerator&&anIsotope==-2)
295     G4int icounter = 0;                        << 294     {
296     G4int icounter_max = 1024;                 << 295       theFinalState = wendtFissionGenerator->ApplyYourself(theTrack, Z, A);
297     while (theFinalState == nullptr)  // Loop  << 
298     {                                          << 
299       icounter++;                              << 
300       if (icounter > icounter_max) {           << 
301         G4cout << "Loop-counter exceeded the t << 
302                << __LINE__ << "th line of " << << 
303         break;                                 << 
304       }                                        << 
305       if (isElastic) {                         << 
306         // Register 0 K cross-section for DBRC << 
307         G4ParticleHPVector* xsforFS = theIsoto << 
308         // Only G4ParticleHPElasticFS has the  << 
309         static_cast<G4ParticleHPElasticFS*>(th << 
310       }                                        << 
311       theFinalState = theFinalStates[it]->Appl << 
312     }                                             296     }
313   }                                            << 
314                                                   297 
315   // G4cout <<"THE IMPORTANT RETURN"<<G4endl;  << 298     // Use the standard procedure if the G4FissionFragmentGenerator model fails
316   // G4cout << "TK G4ParticleHPChannel Elastic << 299     if (!theFinalState)
317   //<< " Z= " << this->GetZ(it) << " A = " <<  << 300     {
318   fManager->GetReactionWhiteBoard()->SetTargA( << 
319   fManager->GetReactionWhiteBoard()->SetTargZ( << 
320   fManager->GetReactionWhiteBoard()->SetTargM( << 
321                                                   301 
322   return theFinalState;                        << 302        G4int icounter=0;
323 }                                              << 303        G4int icounter_max=1024;
                                                   >> 304        while(theFinalState==0) // Loop checking, 11.05.2015, T. Koi
                                                   >> 305        {
                                                   >> 306           icounter++;
                                                   >> 307           if ( icounter > icounter_max ) {
                                                   >> 308        G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
                                                   >> 309              break;
                                                   >> 310           }
                                                   >> 311 //        G4cout << "TESTHP 24 it="<<it<<G4endl;
                                                   >> 312           theFinalState = theFinalStates[it]->ApplyYourself(theTrack);
                                                   >> 313        }
                                                   >> 314     }
324                                                   315 
325 void G4ParticleHPChannel::DumpInfo() const     << 316     //G4cout <<"THE IMPORTANT RETURN"<<G4endl;
326 {                                              << 317     //G4cout << "TK G4ParticleHPChannel Elastic, Capture and Fission Cases "
327   G4cout << " Element: " << theElement->GetNam << 318     //<< " Z= " << this->GetZ(it) << " A = " << this->GetN(it) << G4endl;
328   G4cout << " Directory name: " << theDir << G << 319     G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA( A );
329   G4cout << " FS name: " << theFSType << G4end << 320     G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ( Z );
330   G4cout << " Number of Isotopes: " << niso << << 321     G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargM( M );
331   G4cout << " Have cross sections: " << G4endl << 322 
332   for (int i = 0; i < niso; i++) {             << 323     return theFinalState;
333     G4cout << theFinalStates[i]->HasXsec() <<  << 
334   }                                               324   }
335   G4cout << G4endl;                            << 325 
336   if (theChannelData != nullptr) {             << 326 
337     G4cout << " Cross Section (total for this  << 327 void G4ParticleHPChannel::DumpInfo(){
338     int np = theChannelData->GetVectorLength() << 328 
339     G4cout << np << G4endl;                    << 329   G4cout<<" Element: "<<theElement->GetName()<<G4endl;
340     for (int i = 0; i < np; i++) {             << 330   G4cout<<" Directory name: "<<theDir<<G4endl;
341       G4cout << theChannelData->GetEnergy(i) / << 331   G4cout<<" FS name: "<<theFSType<<G4endl;
                                                   >> 332   G4cout<<" Number of Isotopes: "<<niso<<G4endl;
                                                   >> 333   G4cout<<" Have cross sections: "<<G4endl;
                                                   >> 334   for(int i=0;i<niso;i++){
                                                   >> 335     G4cout<<theFinalStates[i]->HasXsec()<<"  ";
                                                   >> 336   }
                                                   >> 337   G4cout<<G4endl;
                                                   >> 338   if(theChannelData){
                                                   >> 339     G4cout<<" Cross Section (total for this channel):"<<G4endl;
                                                   >> 340     int np=theChannelData->GetVectorLength();
                                                   >> 341     G4cout<<np<<G4endl;
                                                   >> 342     for(int i=0;i<np;i++){
                                                   >> 343       G4cout<<theChannelData->GetEnergy(i)/eV<<"  "<<theChannelData->GetXsec(i)<<G4endl;
342     }                                             344     }
343   }                                               345   }
                                                   >> 346 
344 }                                                 347 }
                                                   >> 348  
                                                   >> 349 
                                                   >> 350 
345                                                   351