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
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "G4QNGamma.hh"
00043 #include "G4HadronicDeprecate.hh"
00044
00045
00046
00047 std::vector<G4int> G4QNGamma::ElementZ;
00048 std::vector<G4double> G4QNGamma::ElProbInMat;
00049 std::vector<std::vector<G4int>*> G4QNGamma::ElIsoN;
00050 std::vector<std::vector<G4double>*>G4QNGamma::IsoProbInEl;
00051
00052
00053 G4QNGamma::G4QNGamma(const G4String& processName)
00054 : G4VDiscreteProcess(processName, fHadronic)
00055 {
00056 G4HadronicDeprecate("G4QNGamma");
00057
00058 EnMomConservation = G4LorentzVector(0.,0.,0.,0.);
00059 nOfNeutrons = 0;
00060 #ifdef debug
00061 G4cout<<"G4QNGamma::Constructor is called"<<G4endl;
00062 #endif
00063 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
00064 }
00065
00066
00067 G4QNGamma::~G4QNGamma()
00068 {
00069
00070
00071 G4int IPIE=IsoProbInEl.size();
00072 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)
00073 {
00074 std::vector<G4double>* SPI=IsoProbInEl[ip];
00075 SPI->clear();
00076 delete SPI;
00077 std::vector<G4int>* IsN=ElIsoN[ip];
00078 IsN->clear();
00079 delete IsN;
00080 }
00081 ElProbInMat.clear();
00082 ElementZ.clear();
00083 IsoProbInEl.clear();
00084 ElIsoN.clear();
00085 }
00086
00087
00088 G4LorentzVector G4QNGamma::GetEnegryMomentumConservation()
00089 {
00090 return EnMomConservation;
00091 }
00092
00093 G4int G4QNGamma::GetNumberOfNeutronsInTarget()
00094 {
00095 return nOfNeutrons;
00096 }
00097
00098
00099
00100
00101
00102 G4double G4QNGamma::GetMeanFreePath(const G4Track& aTrack,G4double,G4ForceCondition* Fc)
00103 {
00104 #ifdef debug
00105 G4cout<<"G4QNGamma::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
00106 #endif
00107 *Fc = NotForced;
00108 #ifdef debug
00109 G4cout<<"G4QNGamma::GetMeanFreePath: Before GetDynPart"<<G4endl;
00110 #endif
00111 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
00112 #ifdef debug
00113 G4cout<<"G4QNGamma::GetMeanFreePath: Before GetDef"<<G4endl;
00114 #endif
00115 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
00116 G4double Momentum = incidentParticle->GetTotalMomentum();
00117 if( !IsApplicable(*incidentParticleDefinition))
00118 {
00119 G4cout<<"-W-G4QNGamma::GetMeanFreePath called for not implemented particle"<<G4endl;
00120 return DBL_MAX;
00121 }
00122 #ifdef debug
00123 G4cout<<"G4QNGamma::GetMeanFreePath: BeforeGetMaterial P="<<Momentum<<G4endl;
00124 #endif
00125
00126 if(Momentum > 500.) return DBL_MAX;
00127
00128 const G4Material* material = aTrack.GetMaterial();
00129 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
00130 const G4ElementVector* theElementVector = material->GetElementVector();
00131 G4int nE=material->GetNumberOfElements();
00132 #ifdef debug
00133 G4cout<<"G4QNGamma::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
00134 #endif
00135
00136 G4VQCrossSection* CSmanager = 0;
00137 G4QNeutronCaptureRatio* capMan = 0;
00138 G4int pPDG =0;
00139 G4double sigma =0.;
00140 if(incidentParticleDefinition == G4Neutron::Neutron())
00141 {
00142 CSmanager=G4QNeutronNuclearCrossSection::GetPointer();
00143 capMan=G4QNeutronCaptureRatio::GetPointer();
00144 #ifdef debug
00145 G4cout<<"G4QNGamma::GetMeanFreePath: CSmanager is defined for neutrons"<<G4endl;
00146 #endif
00147 pPDG=2112;
00148 }
00149 else
00150 {
00151 G4cout<<"-Warning-G4QNGamma::GetMeanFreePath:Particle "
00152 <<incidentParticleDefinition->GetPDGEncoding()<<" isn't a neutron"<<G4endl;
00153 return DBL_MAX;
00154 }
00155
00156 G4QIsotope* Isotopes = G4QIsotope::Get();
00157 G4int IPIE=IsoProbInEl.size();
00158 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)
00159 {
00160 std::vector<G4double>* SPI=IsoProbInEl[ip];
00161 SPI->clear();
00162 delete SPI;
00163 std::vector<G4int>* IsN=ElIsoN[ip];
00164 IsN->clear();
00165 delete IsN;
00166 }
00167 ElProbInMat.clear();
00168 ElementZ.clear();
00169 IsoProbInEl.clear();
00170 ElIsoN.clear();
00171 for(G4int i=0; i<nE; ++i)
00172 {
00173 G4Element* pElement=(*theElementVector)[i];
00174 G4int Z = static_cast<G4int>(pElement->GetZ());
00175 ElementZ.push_back(Z);
00176 G4int isoSize=0;
00177 G4int indEl=0;
00178 G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
00179 if(isoVector) isoSize=isoVector->size();
00180 #ifdef debug
00181 G4cout<<"G4QNGamma::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl;
00182 #endif
00183 if(isoSize)
00184 {
00185 indEl=pElement->GetIndex()+1;
00186 if(!Isotopes->IsDefined(Z,indEl))
00187 {
00188 std::vector<std::pair<G4int,G4double>*>* newAbund =
00189 new std::vector<std::pair<G4int,G4double>*>;
00190 G4double* abuVector=pElement->GetRelativeAbundanceVector();
00191 for(G4int j=0; j<isoSize; j++)
00192 {
00193 G4int N=pElement->GetIsotope(j)->GetN()-Z;
00194 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QNGamma::GetMeanFreePath"
00195 <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00196 G4double abund=abuVector[j];
00197 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00198 #ifdef debug
00199 G4cout<<"G4QNGamma::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
00200 #endif
00201 newAbund->push_back(pr);
00202 }
00203 #ifdef debug
00204 G4cout<<"G4QNGamma::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
00205 #endif
00206 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund);
00207 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
00208 delete newAbund;
00209 }
00210 }
00211 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);
00212 std::vector<G4double>* SPI = new std::vector<G4double>;
00213 IsoProbInEl.push_back(SPI);
00214 std::vector<G4int>* IsN = new std::vector<G4int>;
00215 ElIsoN.push_back(IsN);
00216 G4int nIs=cs->size();
00217 G4double susi=0.;
00218 #ifdef debug
00219 G4cout<<"G4QNGamma::GetMeanFreePath: Before Loop nIs="<<nIs<<G4endl;
00220 #endif
00221 if(nIs) for(G4int j=0; j<nIs; j++)
00222 {
00223 std::pair<G4int,G4double>* curIs=(*cs)[j];
00224 G4int N=curIs->first;
00225 IsN->push_back(N);
00226 #ifdef debug
00227 G4cout<<"G4QNGam::GetMeanFrP: Before CS, P="<<Momentum<<",Z="<<Z<<",N="<<N<<G4endl;
00228 #endif
00229
00230 G4double CSI=CSmanager->GetCrossSection(true, Momentum, Z, N, pPDG) *
00231 capMan->GetRatio(Momentum, Z, N);
00232 #ifdef debug
00233 G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
00234 #endif
00235 curIs->second = CSI;
00236 susi+=CSI;
00237 SPI->push_back(susi);
00238 }
00239 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];
00240 ElProbInMat.push_back(sigma);
00241 }
00242 #ifdef debug
00243 G4cout<<"G4QNGam::GetMeanFrPa: Sigma="<<sigma<<G4endl;
00244 #endif
00245 if(sigma > 0.) return 1./sigma;
00246 return DBL_MAX;
00247 }
00248
00249
00250 G4bool G4QNGamma::IsApplicable(const G4ParticleDefinition& particle)
00251 {
00252 if ( particle == *( G4Neutron::Neutron() ) ) return true;
00253 #ifdef debug
00254 G4cout<<"***G4QNGamma::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
00255 #endif
00256 return false;
00257 }
00258
00259 G4VParticleChange* G4QNGamma::PostStepDoIt(const G4Track& track, const G4Step& step)
00260 {
00261 #ifdef debug
00262 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00263 #endif
00264 static const G4LorentzVector vacuum4M(0.,0.,0.,0.);
00265
00266 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
00267 const G4ParticleDefinition* particle=projHadron->GetDefinition();
00268 #ifdef debug
00269 G4cout<<"G4QNGamma::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
00270 #endif
00271 G4ForceCondition cond=NotForced;
00272 GetMeanFreePath(track, 1., &cond);
00273 #ifdef debug
00274 G4cout<<"G4QNGamma::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
00275 #endif
00276 G4LorentzVector proj4M=projHadron->Get4Momentum();
00277 G4double momentum = projHadron->GetTotalMomentum();
00278 G4double Momentum=proj4M.rho();
00279 if(std::fabs(Momentum-momentum)>.001)
00280 G4cerr<<"*G4QNGamma::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl;
00281 #ifdef debug
00282 G4double mp=proj4M.m();
00283 if(std::fabs(mp-mNeut)>.001)G4cerr<<"*G4QNGamma::PostStDoIt: M="<<mp<<"#"<<mNeut<<G4endl;
00284 G4cout<<"->G4QNGam::PostStDoIt:*called*,4M="<<proj4M<<",P="<<Momentum<<",m="<<mp<<G4endl;
00285 #endif
00286
00287 if (!IsApplicable(*particle) || Momentum > 500.)
00288 {
00289 G4cerr<<"G4QNGamma::PostStepDoIt: Only neutrons with P="<<Momentum<<" < 500"<<G4endl;
00290 return 0;
00291 }
00292 const G4Material* material = track.GetMaterial();
00293 G4int Z=0;
00294 const G4ElementVector* theElementVector = material->GetElementVector();
00295 G4int nE=material->GetNumberOfElements();
00296 #ifdef debug
00297 G4cout<<"G4QNGamma::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
00298 #endif
00299 G4int EPIM=ElProbInMat.size();
00300 #ifdef debug
00301 G4cout<<"G4QNGam::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
00302 #endif
00303 G4int i=0;
00304 if(EPIM>1)
00305 {
00306 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
00307 for(i=0; i<nE; ++i)
00308 {
00309 #ifdef debug
00310 G4cout<<"G4QNGamma::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
00311 #endif
00312 if (rnd<ElProbInMat[i]) break;
00313 }
00314 if(i>=nE) i=nE-1;
00315 }
00316 G4Element* pElement=(*theElementVector)[i];
00317 Z=static_cast<G4int>(pElement->GetZ());
00318 #ifdef debug
00319 G4cout<<"G4QNGamma::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
00320 #endif
00321 if(Z <= 0)
00322 {
00323 G4cerr<<"---Warning---G4QNGamma::PostStepDoIt: Element with Z="<<Z<<G4endl;
00324 if(Z<0) return 0;
00325 }
00326 std::vector<G4double>* SPI = IsoProbInEl[i];
00327 std::vector<G4int>* IsN = ElIsoN[i];
00328 G4int nofIsot=SPI->size();
00329 #ifdef debug
00330 G4cout<<"G4QNGam::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
00331 #endif
00332 G4int j=0;
00333 if(nofIsot>1)
00334 {
00335 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand();
00336 for(j=0; j<nofIsot; ++j)
00337 {
00338 #ifdef debug
00339 G4cout<<"G4QNGamma::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
00340 #endif
00341 if(rndI < (*SPI)[j]) break;
00342 }
00343 if(j>=nofIsot) j=nofIsot-1;
00344 }
00345 G4int N =(*IsN)[j];
00346 #ifdef debug
00347 G4cout<<"G4QNGamma::PostStepDoIt: Z="<<Z<<", j="<<i<<", N(isotope)="<<N<<G4endl;
00348 #endif
00349 G4double kinEnergy= projHadron->GetKineticEnergy();
00350 G4ParticleMomentum dir = projHadron->GetMomentumDirection();
00351
00352
00353
00354
00355
00356
00357
00358
00359 if(N<0)
00360 {
00361 G4cerr<<"-Warning-G4QNGamma::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
00362 return 0;
00363 }
00364 nOfNeutrons=N;
00365 #ifdef debug
00366 G4cout<<"G4QNGamma::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
00367 #endif
00368 aParticleChange.Initialize(track);
00369 G4double weight = track.GetWeight();
00370 #ifdef debug
00371 G4cout<<"G4QNGamma::PostStepDoIt: weight="<<weight<<G4endl;
00372 #endif
00373 G4double localtime = track.GetGlobalTime();
00374 #ifdef debug
00375 G4cout<<"G4QNGamma::PostStepDoIt: localtime="<<localtime<<G4endl;
00376 #endif
00377 G4ThreeVector position = track.GetPosition();
00378 G4TouchableHandle trTouchable = track.GetTouchableHandle();
00379 #ifdef debug
00380 G4cout<<"G4QNGamma::PostStepDoIt: position="<<position<<G4endl;
00381 #endif
00382 G4int targPDG = 90000000 + Z*1000 + N;
00383 G4QPDGCode targQPDG(targPDG);
00384 G4double tM = targQPDG.GetMass();
00385 #ifdef debug
00386 G4cout<<"G4QNGamma::PostStepDoIt: n + targPDG="<<targPDG<<G4endl;
00387 #endif
00388
00389 G4LorentzVector tot4M=G4LorentzVector(0.,0.,0.,tM)+proj4M;
00390 G4double totM2=tot4M.m2();
00391 G4int tZ=Z;
00392 G4int tN=N+1;
00393 G4int resPDG = targPDG + 1;
00394 G4double rM=G4QPDGCode(resPDG).GetMass();
00395 G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,rM);
00396 G4LorentzVector g4M=G4LorentzVector(0.,0.,0.,0.);
00397 #ifdef debug
00398 G4cout<<"G4QNGamma::PostStepDoIt: tM="<<tM << ", rM="<<rM << ", Q="<<tM+mNeut-rM<<G4endl;
00399 #endif
00400 if(!G4QHadron(tot4M).DecayIn2(r4M, g4M))
00401 {
00402
00403
00404
00405 G4cerr<<"-Warning-G4QNGamma::PostStDoIt: tM="<<std::sqrt(totM2)<<" < rM="<<rM<<G4endl;
00406 aParticleChange.ProposeEnergy(kinEnergy);
00407 aParticleChange.ProposeLocalEnergyDeposit(0.);
00408 aParticleChange.ProposeMomentumDirection(dir);
00409 aParticleChange.ProposeTrackStatus(fAlive);
00410 return G4VDiscreteProcess::PostStepDoIt(track,step);
00411 }
00412 #ifdef debug
00413 G4cout<<"G4QNGam::PStDoIt: RA="<<r4M.rho()<<r4M<<", Gamma="<<g4M.rho()<<g4M<<G4endl;
00414 #endif
00415 EnMomConservation = tot4M - r4M - g4M;
00416 aParticleChange.ProposeEnergy(0.);
00417 aParticleChange.ProposeTrackStatus(fStopAndKill);
00418 aParticleChange.SetNumberOfSecondaries(2);
00419
00420 G4ParticleDefinition* theDefinition = G4Gamma::Gamma();
00421 G4DynamicParticle* theGam = new G4DynamicParticle(theDefinition, g4M);
00422 G4Track* capGamma = new G4Track(theGam, localtime, position );
00423 capGamma->SetWeight(weight);
00424 capGamma->SetTouchableHandle(trTouchable);
00425 aParticleChange.AddSecondary(capGamma);
00426
00427
00428 G4int tA=tZ+tN;
00429 if (resPDG==90000001) theDefinition = G4Neutron::Neutron();
00430 else if(resPDG==90001000) theDefinition = G4Proton::Proton();
00431 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(tZ, tA, 0, tZ);
00432 G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition, r4M);
00433 G4Track* scatReN = new G4Track(theReN, localtime, position );
00434 scatReN->SetWeight(weight);
00435 scatReN->SetTouchableHandle(trTouchable);
00436 aParticleChange.AddSecondary(scatReN);
00437
00438 return G4VDiscreteProcess::PostStepDoIt(track, step);
00439 }