00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "G4NeutronHPChannelList.hh"
00033 #include "G4Element.hh"
00034 #include "G4HadFinalState.hh"
00035 #include "G4HadProjectile.hh"
00036 #include "G4NeutronHPFinalState.hh"
00037
00038 G4int G4NeutronHPChannelList::trycounter = 0;
00039
00040 G4NeutronHPChannelList::G4NeutronHPChannelList(G4int n)
00041 {
00042 nChannels = n;
00043 theChannels = new G4NeutronHPChannel * [n];
00044 allChannelsCreated = false;
00045 theInitCount = 0;
00046 }
00047
00048 G4NeutronHPChannelList::G4NeutronHPChannelList()
00049 {
00050 nChannels = 0;
00051 theChannels = 0;
00052 allChannelsCreated = false;
00053 theInitCount = 0;
00054 }
00055
00056 G4NeutronHPChannelList::~G4NeutronHPChannelList()
00057 {
00058 if(theChannels!=0)
00059 {
00060 for(G4int i=0;i<nChannels; i++)
00061 {
00062 delete theChannels[i];
00063 }
00064 delete [] theChannels;
00065 }
00066 }
00067
00068 #include "G4NeutronHPThermalBoost.hh"
00069 #include "G4NeutronHPManager.hh"
00070 G4HadFinalState * G4NeutronHPChannelList::ApplyYourself(const G4Element * , const G4HadProjectile & aTrack)
00071 {
00072 G4NeutronHPThermalBoost aThermalE;
00073 G4int i, ii;
00074
00075 G4int numberOfIsos(0);
00076 for(ii=0; ii<nChannels; ii++)
00077 {
00078 numberOfIsos = theChannels[ii]->GetNiso();
00079 if(numberOfIsos!=0) break;
00080 }
00081 G4double * running= new G4double [numberOfIsos];
00082 running[0] = 0;
00083 for(i=0;i<numberOfIsos; i++)
00084 {
00085 if(i!=0) running[i] = running[i-1];
00086 for(ii=0; ii<nChannels; ii++)
00087 {
00088 if(theChannels[ii]->HasAnyData(i))
00089 {
00090 running[i] +=theChannels[ii]->GetWeightedXsec(aThermalE.GetThermalEnergy(aTrack,
00091 theChannels[ii]->GetN(i),
00092 theChannels[ii]->GetZ(i),
00093 aTrack.GetMaterial()->GetTemperature()),
00094 i);
00095 }
00096 }
00097 }
00098 G4int isotope=nChannels-1;
00099 G4double random=G4UniformRand();
00100 for(i=0;i<numberOfIsos; i++)
00101 {
00102 isotope = i;
00103
00104 if(running[numberOfIsos-1] != 0) if(random<running[i]/running[numberOfIsos-1]) break;
00105 }
00106 delete [] running;
00107
00108
00109 running = new G4double[nChannels];
00110 running[0]=0;
00111 G4int targA=-1;
00112 G4int targZ=-1;
00113 for(i=0; i<nChannels; i++)
00114 {
00115 if(i!=0) running[i] = running[i-1];
00116 if(theChannels[i]->HasAnyData(isotope))
00117 {
00118 running[i] += theChannels[i]->GetFSCrossSection(aThermalE.GetThermalEnergy(aTrack,
00119 theChannels[i]->GetN(isotope),
00120 theChannels[i]->GetZ(isotope),
00121 aTrack.GetMaterial()->GetTemperature()),
00122 isotope);
00123 targA=(G4int)theChannels[i]->GetN(isotope);
00124 targZ=(G4int)theChannels[i]->GetZ(isotope);
00125 }
00126 }
00127
00128
00129 if ( running[nChannels-1] == 0 )
00130 {
00131
00132 if ( targA == -1 && targZ == -1 ) {
00133 throw G4HadronicException(__FILE__, __LINE__, "NeutronHP model encounter lethal discrepancy with cross section data");
00134 }
00135
00136
00137 G4cout << "Warning from NeutronHP: could not find proper reaction channel. This may cause by inconsistency between cross section and model. Unchanged final states are returned." << G4endl;
00138 unChanged.Clear();
00139
00140
00141 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( targZ , targA , 0.0 );
00142 G4DynamicParticle* targ_dp = new G4DynamicParticle( targ_pd , G4ThreeVector(1,0,0), 0.0 );
00143 unChanged.SetEnergyChange(aTrack.GetKineticEnergy());
00144 unChanged.SetMomentumChange(aTrack.Get4Momentum().vect() );
00145 unChanged.AddSecondary(targ_dp);
00146
00147 G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA( targA );
00148 G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ( targZ );
00149 return &unChanged;
00150 }
00151
00152
00153
00154 G4int lChan=0;
00155 random=G4UniformRand();
00156 for(i=0; i<nChannels; i++)
00157 {
00158 lChan = i;
00159 if(running[nChannels-1] != 0) if(random<running[i]/running[nChannels-1]) break;
00160 }
00161 delete [] running;
00162 return theChannels[lChan]->ApplyYourself(aTrack, isotope);
00163 }
00164
00165 void G4NeutronHPChannelList::Init(G4Element * anElement, const G4String & dirName)
00166 {
00167 theDir = dirName;
00168
00169 theElement = anElement;
00170
00171 ;
00172 }
00173
00174 void G4NeutronHPChannelList::Register(G4NeutronHPFinalState * theFS,
00175 const G4String & aName)
00176 {
00177 if(!allChannelsCreated)
00178 {
00179 if(nChannels!=0)
00180 {
00181 G4NeutronHPChannel ** theBuffer = new G4NeutronHPChannel * [nChannels+1];
00182 G4int i;
00183 for(i=0; i<nChannels; i++)
00184 {
00185 theBuffer[i] = theChannels[i];
00186 }
00187 delete [] theChannels;
00188 theChannels = theBuffer;
00189 }
00190 else
00191 {
00192 theChannels = new G4NeutronHPChannel * [nChannels+1];
00193 }
00194 G4String name;
00195 name = aName+"/";
00196 theChannels[nChannels] = new G4NeutronHPChannel;
00197 theChannels[nChannels]->Init(theElement, theDir, name);
00198 nChannels++;
00199 }
00200
00201
00202
00203
00204 theChannels[theInitCount]->Register(theFS);
00205 theInitCount++;
00206 }