#include <G4QIonIonElastic.hh>
Inheritance diagram for G4QIonIonElastic:
Public Member Functions | |
G4QIonIonElastic (const G4String &processName="CHIPS_IonIonElasticScattering") | |
~G4QIonIonElastic () | |
G4bool | IsApplicable (const G4ParticleDefinition &particle) |
G4double | GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) |
G4VParticleChange * | PostStepDoIt (const G4Track &aTrack, const G4Step &aStep) |
G4LorentzVector | GetEnegryMomentumConservation () |
G4int | GetNumberOfNeutronsInTarget () |
Definition at line 66 of file G4QIonIonElastic.hh.
G4QIonIonElastic::G4QIonIonElastic | ( | const G4String & | processName = "CHIPS_IonIonElasticScattering" |
) |
Definition at line 60 of file G4QIonIonElastic.cc.
References G4cout, G4endl, G4HadronicDeprecate, G4VProcess::GetProcessName(), and G4VProcess::verboseLevel.
00060 : 00061 G4VDiscreteProcess(processName, fHadronic) 00062 { 00063 G4HadronicDeprecate("G4QIonIonElastic"); 00064 00065 #ifdef debug 00066 G4cout<<"G4QIonIonElastic::Constructor is called processName="<<processName<<G4endl; 00067 #endif 00068 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl; 00069 //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max) 00070 }
G4QIonIonElastic::~G4QIonIonElastic | ( | ) |
G4LorentzVector G4QIonIonElastic::GetEnegryMomentumConservation | ( | ) | [inline] |
G4double G4QIonIonElastic::GetMeanFreePath | ( | const G4Track & | aTrack, | |
G4double | previousStepSize, | |||
G4ForceCondition * | condition | |||
) | [virtual] |
Implements G4VDiscreteProcess.
Definition at line 78 of file G4QIonIonElastic.cc.
References G4Alpha::Alpha(), DBL_MAX, G4Deuteron::Deuteron(), G4cerr, G4cout, G4endl, G4GenericIon::GenericIon(), G4QIsotope::Get(), G4ParticleDefinition::GetBaryonNumber(), G4VQCrossSection::GetCrossSection(), G4QIsotope::GetCSVector(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterial(), G4QIsotope::GetMeanCrossSection(), G4Material::GetNumberOfElements(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGEncoding(), G4QProtonElasticCrossSection::GetPointer(), G4QIonIonCrossSection::GetPointer(), G4DynamicParticle::GetTotalMomentum(), G4Material::GetVecNbOfAtomsPerVolume(), G4He3::He3(), G4QIsotope::InitElement(), IsApplicable(), G4QIsotope::IsDefined(), NotForced, and G4Triton::Triton().
Referenced by PostStepDoIt().
00080 { 00081 static const G4double mProt = G4QPDGCode(2212).GetMass()/MeV;// CHIPS proton Mass in MeV 00082 *Fc = NotForced; 00083 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle(); 00084 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition(); 00085 if( !IsApplicable(*incidentParticleDefinition)) 00086 G4cout<<"-Warning-G4QIonIonElastic::GetMeanFreePath: notImplementedParticle"<<G4endl; 00087 // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material 00088 G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle 00089 #ifdef debug 00090 G4double KinEn = incidentParticle->GetKineticEnergy(); 00091 G4cout<<"G4QIonIonElastic::GetMeanFreePath: kinE="<<KinEn<<",Mom="<<Momentum<<G4endl; 00092 #endif 00093 const G4Material* material = aTrack.GetMaterial(); // Get the current material 00094 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume(); 00095 const G4ElementVector* theElementVector = material->GetElementVector(); 00096 G4int nE=material->GetNumberOfElements(); 00097 #ifdef debug 00098 G4cout<<"G4QIonIonElastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl; 00099 #endif 00100 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer(); 00101 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); 00102 G4int pPDG=0; 00103 // Probably enough: pPDG=incidentParticleDefinition->GetPDGEncoding(); 00104 G4int pZ=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge()); 00105 G4int pA=incidentParticleDefinition->GetBaryonNumber(); 00106 if (incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020; 00107 else if (incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 1000020040; 00108 else if (incidentParticleDefinition == G4Triton::Triton() ) pPDG = 1000010030; 00109 else if (incidentParticleDefinition == G4He3::He3() ) pPDG = 1000020030; 00110 else if (incidentParticleDefinition == G4GenericIon::GenericIon() || (pZ > 0 && pA > 0)) 00111 { 00112 pPDG=incidentParticleDefinition->GetPDGEncoding(); 00113 #ifdef debug 00114 G4int prPDG=1000000000+10000*pZ+10*pA; 00115 G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl; 00116 #endif 00117 } 00118 else G4cout<<"-Warning-G4QIonIonElastic::GetMeanFreePath:Unknown projectile Ion"<<G4endl; 00119 Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A 00120 G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton 00121 G4double sigma=0.; // Sums over elements for the material 00122 G4int IPIE=IsoProbInEl.size(); // How many old elements? 00123 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) 00124 { 00125 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector 00126 SPI->clear(); 00127 delete SPI; 00128 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector 00129 IsN->clear(); 00130 delete IsN; 00131 } 00132 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) 00133 ElementZ.clear(); // Clear the body vector for Z of Elements 00134 IsoProbInEl.clear(); // Clear the body vector for SPI 00135 ElIsoN.clear(); // Clear the body vector for N of Isotopes 00136 for(G4int i=0; i<nE; ++i) 00137 { 00138 G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element 00139 G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element 00140 ElementZ.push_back(Z); // Remember Z of the Element 00141 G4int isoSize=0; // The default for the isoVectorLength is 0 00142 G4int indEl=0; // Index of non-natural element or 0(default) 00143 G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect 00144 if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector 00145 #ifdef debug 00146 G4cout<<"G4QIonIonElastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; 00147 #endif 00148 if(isoSize) // The Element has non-trivial abundance set 00149 { 00150 indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order 00151 #ifdef debug 00152 G4cout<<"G4QIIEl::GetMFP:iE="<<indEl<<", def="<<Isotopes->IsDefined(Z,indEl)<<G4endl; 00153 #endif 00154 if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define 00155 { 00156 std::vector<std::pair<G4int,G4double>*>* newAbund = 00157 new std::vector<std::pair<G4int,G4double>*>; 00158 G4double* abuVector=pElement->GetRelativeAbundanceVector(); 00159 for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes 00160 { 00161 G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z ! 00162 if(pElement->GetIsotope(j)->GetZ()!=Z) G4cerr<<"G4QIonIonEl::GetMeanFreePath Z=" 00163 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl; 00164 G4double abund=abuVector[j]; 00165 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund); 00166 #ifdef debug 00167 G4cout<<"G4QIonIonElastic::GetMeanFP:pair#="<<j<<",N="<<N<<",ab="<<abund<<G4endl; 00168 #endif 00169 newAbund->push_back(pr); 00170 } 00171 #ifdef debug 00172 G4cout<<"G4QIonIonElastic::GetMeanFP: pairVectorLength="<<newAbund->size()<<G4endl; 00173 #endif 00174 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd 00175 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary 00176 delete newAbund; // Was "new" in the beginning of the name space 00177 } 00178 } 00179 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer 00180 std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector 00181 IsoProbInEl.push_back(SPI); 00182 std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector 00183 ElIsoN.push_back(IsN); 00184 G4int nIs=cs->size(); // A#Of Isotopes in the Element 00185 #ifdef debug 00186 G4cout<<"G4QIonIonEl::GetMFP:=***=> #isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl; 00187 #endif 00188 G4double susi=0.; // sum of CS over isotopes 00189 if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El 00190 { 00191 std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice 00192 G4int N=curIs->first; // #of Neuterons in the isotope j of El i 00193 IsN->push_back(N); // Remember Min N for the Element 00194 #ifdef debug 00195 G4cout<<"G4QIIE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl; 00196 #endif 00197 G4bool ccsf=false; // Extract elastic Ion-Ion cross-section 00198 #ifdef debug 00199 G4cout<<"G4QIonIonElastic::GMFP: GetCS #1 j="<<j<<G4endl; 00200 #endif 00201 G4double CSI=0.; 00202 if(Z==1 && !N) // Proton target. Reversed kinematics. 00203 { 00204 G4double projM=G4QPDGCode(2212).GetNuclMass(pZ,pA-pZ,0); // Mass of the projNucleus 00205 CSI=PCSmanager->GetCrossSection(true, Momentum*mProt/projM, Z, N, 2212); // Ap CS 00206 } 00207 else CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG); // CS(j,i) for isotope 00208 #ifdef debug 00209 G4cout<<"G4QIonIonElastic::GMFP: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="<<Momentum 00210 <<", XSec="<<CSI/millibarn<<G4endl; 00211 #endif 00212 curIs->second = CSI; 00213 susi+=CSI; // Make a sum per isotopes 00214 SPI->push_back(susi); // Remember summed cross-section 00215 } // End of temporary initialization of the cross sections in the G4QIsotope singeltone 00216 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV) 00217 #ifdef debug 00218 G4cout<<"G4QIonIonEl::GMFP:<S>="<<Isotopes->GetMeanCrossSection(Z,indEl)<<",AddToSig=" 00219 <<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl; 00220 #endif 00221 ElProbInMat.push_back(sigma); 00222 } // End of LOOP over Elements 00223 // Check that cross section is not zero and return the mean free path 00224 #ifdef debug 00225 G4cout<<"G4QIonIonElastic::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl; 00226 #endif 00227 if(sigma > 0.00000000001) return 1./sigma; // Mean path [distance] 00228 return DBL_MAX; 00229 }
G4int G4QIonIonElastic::GetNumberOfNeutronsInTarget | ( | ) | [inline] |
G4bool G4QIonIonElastic::IsApplicable | ( | const G4ParticleDefinition & | particle | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 232 of file G4QIonIonElastic.cc.
References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4cout, G4endl, G4GenericIon::GenericIon(), G4ParticleDefinition::GetBaryonNumber(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGEncoding(), G4He3::He3(), and G4Triton::Triton().
Referenced by GetMeanFreePath(), and PostStepDoIt().
00233 { 00234 G4int Z=static_cast<G4int>(particle.GetPDGCharge()); 00235 G4int A=particle.GetBaryonNumber(); 00236 if (particle == *( G4Deuteron::Deuteron() )) return true; 00237 else if (particle == *( G4Alpha::Alpha() )) return true; 00238 else if (particle == *( G4Triton::Triton() )) return true; 00239 else if (particle == *( G4He3::He3() )) return true; 00240 else if (particle == *( G4GenericIon::GenericIon() )) return true; 00241 else if (Z > 0 && A > 0) return true; 00242 #ifdef debug 00243 G4cout<<"***>>G4QIonIonElastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<", A=" 00244 <<A<<", Z="<<Z<<G4endl; 00245 #endif 00246 return false; 00247 }
G4VParticleChange * G4QIonIonElastic::PostStepDoIt | ( | const G4Track & | aTrack, | |
const G4Step & | aStep | |||
) | [virtual] |
Reimplemented from G4VDiscreteProcess.
Definition at line 249 of file G4QIonIonElastic.cc.
References G4ParticleChange::AddSecondary(), G4Alpha::Alpha(), G4VProcess::aParticleChange, G4Deuteron::Deuteron(), G4ParticleTable::FindIon(), fStopButAlive, G4cerr, G4cout, G4endl, G4UniformRand, G4GenericIon::GenericIon(), G4QCHIPSWorld::Get(), G4DynamicParticle::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4VQCrossSection::GetCrossSection(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4VQCrossSection::GetExchangeT(), G4Track::GetGlobalTime(), G4VQCrossSection::GetHMaxT(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4QPDGCode::GetMass(), G4Track::GetMaterial(), GetMeanFreePath(), G4DynamicParticle::GetMomentumDirection(), G4Material::GetNumberOfElements(), G4ParticleDefinition::GetParticleName(), G4QCHIPSWorld::GetParticles(), G4ParticleDefinition::GetParticleSubType(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetParticleType(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGEncoding(), G4QIonIonCrossSection::GetPointer(), G4QNeutronElasticCrossSection::GetPointer(), G4QProtonElasticCrossSection::GetPointer(), G4Track::GetPosition(), G4VQCrossSection::GetSlope(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetTouchableHandle(), G4Track::GetWeight(), G4Element::GetZ(), G4He3::He3(), G4ParticleChange::Initialize(), IsApplicable(), NotForced, position, G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4Proton::Proton(), G4DynamicParticle::Set4Momentum(), G4DynamicParticle::SetDefinition(), G4Track::SetTouchableHandle(), G4Track::SetWeight(), and G4Triton::Triton().
00250 { 00251 static const G4double mProt= G4QPDGCode(2212).GetMass(); // CHIPS proton mass in MeV 00252 static const G4double fm2MeV2 = 3*38938./1.09; // (3/1.09)*(hc)^2 in fm^2*MeV^2 00253 static G4bool CWinit = true; // CHIPS Warld needs to be initted 00254 if(CWinit) 00255 { 00256 CWinit=false; 00257 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) 00258 } 00259 //------------------------------------------------------------------------------------- 00260 const G4DynamicParticle* projHadron = track.GetDynamicParticle(); 00261 const G4ParticleDefinition* particle=projHadron->GetDefinition(); 00262 #ifdef debug 00263 G4cout<<"G4QIonIonElastic::PostStepDoIt: Before the GetMeanFreePath is called In4M=" 00264 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type=" 00265 <<particle->GetParticleType()<<", Subtp="<<particle->GetParticleSubType()<<G4endl; 00266 #endif 00267 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV! 00268 G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV 00269 G4double Momentum = proj4M.rho(); // @@ Just for the test purposes 00270 if(std::fabs(Momentum-momentum)>.000001) 00271 G4cerr<<"-Warn-G4QIonIonElastic::PostStepDoIt:P(IU)="<<Momentum<<"#"<<momentum<<G4endl; 00272 G4double pM2=proj4M.m2(); // in MeV^2 00273 G4double pM=std::sqrt(pM2); // in MeV 00274 #ifdef pdebug 00275 G4cout<<"G4QIonIonEl::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",pM="<<pM<<G4endl; 00276 #endif 00277 if (!IsApplicable(*particle)) // Check applicability 00278 { 00279 G4cerr<<"G4QIonIonElastic::PostStepDoIt: Only NA elastic is implemented."<<G4endl; 00280 return 0; 00281 } 00282 const G4Material* material = track.GetMaterial(); // Get the current material 00283 const G4ElementVector* theElementVector = material->GetElementVector(); 00284 G4int nE=material->GetNumberOfElements(); 00285 #ifdef debug 00286 G4cout<<"G4QIonIonElastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl; 00287 #endif 00288 // Probably enough: projPDG=particle->GetPDGEncoding(); 00289 G4int projPDG=0; // CHIPS PDG Code for the captured hadron 00290 G4int pZ=static_cast<G4int>(particle->GetPDGCharge()); 00291 G4int pA=particle->GetBaryonNumber(); 00292 if (particle == G4Deuteron::Deuteron() ) projPDG= 1000010020; 00293 else if (particle == G4Alpha::Alpha() ) projPDG= 1000020040; 00294 else if (particle == G4Triton::Triton() ) projPDG= 1000010030; 00295 else if (particle == G4He3::He3() ) projPDG= 1000020030; 00296 else if (particle == G4GenericIon::GenericIon() || (pZ > 0 && pA > 0)) 00297 { 00298 projPDG=particle->GetPDGEncoding(); 00299 #ifdef debug 00300 G4int prPDG=1000000000+10000*pZ+10*pA; 00301 G4cout<<"G4QIonIonElastic::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl; 00302 #endif 00303 } 00304 else G4cout<<"-Warning-G4QIonIonElastic::PostStepDoIt:Unknown projectile Ion"<<G4endl; 00305 #ifdef debug 00306 G4int prPDG=particle->GetPDGEncoding(); 00307 G4cout<<"G4QIonIonElastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl; 00308 #endif 00309 if(!projPDG) 00310 { 00311 G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:Undefined interactingNucleus"<<G4endl; 00312 return 0; 00313 } 00314 G4int pN=pA-pZ; // Projectile N 00315 G4int EPIM=ElProbInMat.size(); 00316 #ifdef debug 00317 G4cout<<"G4QIonIonElastic::PSDI:m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl; 00318 #endif 00319 G4int i=0; 00320 if(EPIM>1) 00321 { 00322 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand(); 00323 for(i=0; i<nE; ++i) 00324 { 00325 #ifdef debug 00326 G4cout<<"G4QIonIonElastic::PSDI: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl; 00327 #endif 00328 if (rnd<ElProbInMat[i]) break; 00329 } 00330 if(i>=nE) i=nE-1; // Top limit for the Element 00331 } 00332 G4Element* pElement=(*theElementVector)[i]; 00333 G4int tZ=static_cast<G4int>(pElement->GetZ()); 00334 #ifdef debug 00335 G4cout<<"G4QIonIonElastic::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl; 00336 #endif 00337 if(tZ<=0) 00338 { 00339 G4cerr<<"---Warning---G4QIonIonElastic::PostStepDoIt: Element with Z="<<tZ<<G4endl; 00340 if(tZ<0) return 0; 00341 } 00342 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes 00343 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i] 00344 G4int nofIsot=SPI->size(); // #of isotopes in the element i 00345 #ifdef debug 00346 G4cout<<"G4QIonIonElastic::PosStDoIt: nI="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl; 00347 #endif 00348 G4int j=0; 00349 if(nofIsot>1) 00350 { 00351 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element 00352 for(j=0; j<nofIsot; ++j) 00353 { 00354 #ifdef debug 00355 G4cout<<"G4QIonIonElastic::PostStDI: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl; 00356 #endif 00357 if(rndI < (*SPI)[j]) break; 00358 } 00359 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope 00360 } 00361 G4int tN =(*IsN)[j]; ; // Randomized number of neutrons 00362 #ifdef debug 00363 G4cout<<"G4QIonIonElastic::PostStepDoIt:j="<<i<<",N(isotope)="<<tN<<",MeV="<<MeV<<G4endl; 00364 #endif 00365 if(tN<0) 00366 { 00367 G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:IsotopeZ="<<tZ<<" & 0>N="<<tN<<G4endl; 00368 return 0; 00369 } 00370 nOfNeutrons=tN; // Remember it for the energy-momentum check 00371 #ifdef debug 00372 G4cout<<"G4QIonIonElastic::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl; 00373 #endif 00374 aParticleChange.Initialize(track); 00375 #ifdef debug 00376 G4cout<<"G4QIonIonElastic::PostStepDoIt: track is initialized"<<G4endl; 00377 #endif 00378 G4double weight = track.GetWeight(); 00379 G4double localtime = track.GetGlobalTime(); 00380 G4ThreeVector position = track.GetPosition(); 00381 #ifdef debug 00382 G4cout<<"G4QIonIonElastic::PostStepDoIt: before Touchable extraction"<<G4endl; 00383 #endif 00384 G4TouchableHandle trTouchable = track.GetTouchableHandle(); 00385 #ifdef debug 00386 G4cout<<"G4QIonIonElastic::PostStepDoIt: Touchable is extracted"<<G4endl; 00387 #endif 00388 // Definitions for do nothing 00389 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?) 00390 G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector 00391 // !! Exception for the proton target !! 00392 G4bool revkin=false; 00393 G4ThreeVector bvel(0.,0.,0.); 00394 G4int tA=tZ+tN; 00395 G4int targPDG=90000000+tZ*1000+tN; // CHIPS PDG Code of the target nucleus 00396 G4QPDGCode targQPDG(targPDG); // @@ one can use G4Ion & get rid of CHIPS World 00397 G4double tM=targQPDG.GetMass(); // CHIPS target mass in MeV 00398 G4LorentzVector targ4M(0.,0.,0.,tM); 00399 if(tZ==1 && !tN) // Update parameters in DB and trans to Antilab 00400 { 00401 G4ForceCondition cond=NotForced; 00402 GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters? 00403 #ifdef debug 00404 G4cout<<"G4QIonIonElastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl; 00405 #endif 00406 revkin=true; 00407 tZ=pZ; 00408 tN=pN; 00409 tA=tZ+tN; 00410 tM=pM; 00411 pZ=1; // @@ Is that necessary ?? 00412 pN=0; // @@ Is that necessary ?? 00413 pA=1; // @@ Is that necessary ?? 00414 pM=mProt; 00415 bvel=proj4M.vect()/proj4M.e(); // Lab->Antilab transition boost velocity 00416 proj4M=targ4M.boost(-bvel); // Proton 4-mom in Antilab 00417 targ4M=G4LorentzVector(0.,0.,0.,tM); // Projectile nucleus 4-mom in Antilab 00418 Momentum = proj4M.rho(); // Recalculate Momentum in Antilab 00419 } 00420 G4LorentzVector tot4M=proj4M+targ4M; // Total 4-mom of the reaction 00421 #ifdef debug 00422 G4cout<<"G4QIonIonElastic::PostStDI: tM="<<tM<<", p4M="<<proj4M<<", t4M="<<tot4M<<G4endl; 00423 #endif 00424 EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation 00425 G4VQCrossSection* PELmanager=G4QProtonElasticCrossSection::GetPointer(); 00426 G4VQCrossSection* NELmanager=G4QNeutronElasticCrossSection::GetPointer(); 00427 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer(); 00428 // @@ Probably this is not necessary any more 00429 #ifdef debug 00430 G4cout<<"G4QIIEl::PSDI:f, P="<<Momentum<<",Z="<<tZ<<",N="<<tN<<",tPDG="<<projPDG<<G4endl; 00431 #endif 00432 // false means elastic cross-section 00433 G4double xSec=0.; // Proto of Recalculated Cross Section 00434 if(revkin) xSec=PELmanager->GetCrossSection(false, Momentum, tZ, tN, 2212); 00435 else xSec=CSmanager->GetCrossSection(false, Momentum, tZ, tN, projPDG); 00436 #ifdef debug 00437 G4cout<<"G4QIIEl::PSDI: pPDG="<<projPDG<<",P="<<Momentum<<",CS="<<xSec/millibarn<<G4endl; 00438 #endif 00439 #ifdef nandebug 00440 if(xSec>0. || xSec<0. || xSec==0); 00441 else G4cout<<"-NaN-Warning-G4QIonIonElastic::PostStDoIt: xSec="<<xSec/millibarn<<G4endl; 00442 #endif 00443 // @@ check a possibility to separate p, n, or alpha (!) 00444 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing 00445 { 00446 #ifdef pdebug 00447 G4cerr<<"-Warning-G4QIonIonElastic::PostStDoIt: *Zero cross-section* PDG="<<projPDG 00448 <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl; 00449 #endif 00450 //Do Nothing Action insead of the reaction 00451 aParticleChange.ProposeEnergy(kinEnergy); 00452 aParticleChange.ProposeLocalEnergyDeposit(0.); 00453 aParticleChange.ProposeMomentumDirection(dir) ; 00454 return G4VDiscreteProcess::PostStepDoIt(track,step); 00455 } 00456 G4double mint=0; 00457 G4double maxt=0; 00458 G4double dtM=tM+tM; 00459 if(revkin) 00460 { 00461 mint=PELmanager->GetExchangeT(tZ,tN,2212); // functional randomized -t in MeV^2 00462 maxt=PELmanager->GetHMaxT(); 00463 } 00464 else 00465 { 00466 G4double PA=Momentum*pA; 00467 G4double PA2=PA*PA; 00468 maxt=dtM*PA2/(std::sqrt(PA2+pM2)+tM/2+pM2/dtM); 00469 #ifdef pdebug 00470 G4cout<<"G4QIonIonElastic::PostStDoIt:pPDG="<<projPDG<<",tPDG="<<targPDG<<",P=" 00471 <<Momentum<<",CS="<<xSec<<",maxt="<<maxt<<G4endl; 00472 #endif 00473 xSec=PELmanager->GetCrossSection(false, Momentum, 1, 0, 2212);// pp=nn 00474 G4double B1=PELmanager->GetSlope(1,0,2212); // slope for pp=nn 00475 xSec=NELmanager->GetCrossSection(false, Momentum, 1, 0, 2112);// np=pn 00476 G4double B2 =NELmanager->GetSlope(1,0,2112); // slope for np=pn 00477 G4double mB =((pZ*tZ+pN*tN)*B1+(pZ*tN+pN*tZ)*B2)/(pA+tA); 00478 G4double pR2=std::pow(pA+4.,.305)/fm2MeV2; 00479 G4double tR2=std::pow(tA+4.,.305)/fm2MeV2; 00480 G4double eB =mB+pR2+tR2; 00481 mint=-std::log(1.-G4UniformRand()*(1.-std::exp(-eB*maxt)))/eB; 00482 mint+=mint; 00483 #ifdef pdebug 00484 G4cout<<"G4QIonIonElastic::PostStDoIt:B1="<<B1<<",B2="<<B2<<",mB="<<mB 00485 <<",pR2="<<pR2<<",tR2="<<tR2<<",eB="<<eB<<",mint="<<mint<<G4endl; 00486 #endif 00487 } 00488 #ifdef nandebug 00489 if(mint>-.0000001); 00490 else G4cout<<"-Warning-G4QIonIonElastic::PostStDoIt:-t="<<mint<<G4endl; 00491 #endif 00492 G4double cost=1.-mint/maxt; // cos(theta) in CMS 00493 // 00494 #ifdef ppdebug 00495 G4cout<<"G4QIonIonElastic::PoStDoI:t="<<mint<<",dpcm2="<<CSmanager->GetHMaxT()<<",Ek=" 00496 <<kinEnergy<<",tM="<<tM<<",pM="<<pM<<",cost="<<cost<<G4endl; 00497 #endif 00498 if(cost>1. || cost<-1. || !(cost>-1. || cost<1.)) 00499 { 00500 if(cost>1.000001 || cost<-1.000001 || !(cost>-1. || cost<1.)) 00501 { 00502 G4double tM2=tM*tM; // Squared target mass 00503 G4double pEn=pM+kinEnergy; // tot projectile Energy in MeV 00504 G4double sM=dtM*pEn+tM2+pM2; // Mondelstam s 00505 G4double twop2cm=(tM2+tM2)*(pEn*pEn-pM2)/sM;// Max_t/2 (2*p^2_cm) 00506 G4cout<<"-Warning-G4QIonIonElastic::PoStDI:cos="<<cost<<",t="<<mint<<",T="<<kinEnergy 00507 <<",tM="<<tM<<",tmax="<<2*kinEnergy*tM<<",p="<<projPDG<<",t="<<targPDG<<G4endl; 00508 G4cout<<"G4QIonIonElastic::PSDI:dpcm2="<<twop2cm<<"="<<CSmanager->GetHMaxT()<<G4endl; 00509 } 00510 if (cost>1.) cost=1.; 00511 else if(cost<-1.) cost=-1.; 00512 } 00513 G4LorentzVector scat4M(0.,0.,0.,pM); // 4-mom of the scattered projectile 00514 G4LorentzVector reco4M(0.,0.,0.,tM); // 4-mom of the recoil target 00515 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-tM-pM)*.01); 00516 if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) 00517 { 00518 G4cout<<"-Warning-G4QIonIonE::PSDI:t4M="<<tot4M<<",pM="<<pM<<",tM="<<tM<<",cost=" 00519 <<cost<<G4endl; 00520 } 00521 if(revkin) 00522 { 00523 G4LorentzVector tmpLV=scat4M.boost(bvel); // Recoil target Back to LS 00524 scat4M=reco4M.boost(bvel); // Scattered Progectile Back to LS 00525 reco4M=tmpLV; 00526 pM=tM; 00527 tM=mProt; 00528 } 00529 #ifdef debug 00530 G4cout<<"G4QIonIonElast::PSDI:s4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl; 00531 G4cout<<"G4QIonIonElastic::PSDI:scatE="<<scat4M.e()-pM<<",recoE="<<reco4M.e()-tM<<",d4M=" 00532 <<tot4M-scat4M-reco4M<<G4endl; 00533 #endif 00534 // Update G4VParticleChange for the scattered projectile 00535 G4double finE=scat4M.e()-pM; // Final kinetic energy of the scattered proton 00536 if(finE>0.0) aParticleChange.ProposeEnergy(finE); 00537 else 00538 { 00539 if(finE<-1.e-8 || !(finE>-1.||finE<1.)) // NAN or negative 00540 G4cerr<<"*Warning*G4QIonIonElastic::PostStDoIt: Zero or negative scattered E="<<finE 00541 <<", s4M="<<scat4M<<", r4M="<<reco4M<<", d4M="<<tot4M-scat4M-reco4M<<G4endl; 00542 aParticleChange.ProposeEnergy(0.) ; 00543 aParticleChange.ProposeTrackStatus(fStopButAlive); 00544 } 00545 G4ThreeVector findir=scat4M.vect()/scat4M.rho(); // Unit vector in new direction 00546 aParticleChange.ProposeMomentumDirection(findir); // new direction for the scattered part 00547 EnMomConservation-=scat4M; // It must be initialized by (pE+tM,pP) 00548 // This is how in general the secondary should be identified 00549 G4DynamicParticle* theSec = new G4DynamicParticle; // A secondary for the recoil hadron 00550 #ifdef pdebug 00551 G4cout<<"G4QIonIonElastic::PostStepDoIt: Ion tZ="<<tZ<<", tA="<<tA<<G4endl; 00552 #endif 00553 G4ParticleDefinition* theDefinition=0; 00554 if(revkin) theDefinition = G4Proton::Proton(); 00555 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(tZ,tA,0,tZ); 00556 if(!theDefinition)G4cout<<"-Warning-G4QIonIonElastic::PoStDI:drop PDG="<<targPDG<<G4endl; 00557 #ifdef pdebug 00558 G4cout<<"G4QIonIonElastic::PoStDI:RecoilName="<<theDefinition->GetParticleName()<<G4endl; 00559 #endif 00560 theSec->SetDefinition(theDefinition); 00561 EnMomConservation-=reco4M; 00562 #ifdef tdebug 00563 G4cout<<"G4QIonIonElastic::PSD:"<<targPDG<<reco4M<<reco4M.m()<<EnMomConservation<<G4endl; 00564 #endif 00565 theSec->Set4Momentum(reco4M); 00566 #ifdef debug 00567 G4ThreeVector curD=theSec->GetMomentumDirection(); 00568 G4double curM=theSec->GetMass(); 00569 G4double curE=theSec->GetKineticEnergy()+curM; 00570 G4cout<<"G4QIonIonElastic::PSDI: p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl; 00571 #endif 00572 // Make a recoil nucleus 00573 G4Track* aNewTrack = new G4Track(theSec, localtime, position ); 00574 aNewTrack->SetWeight(weight); // weighted 00575 aNewTrack->SetTouchableHandle(trTouchable); 00576 aParticleChange.AddSecondary( aNewTrack ); 00577 #ifdef debug 00578 G4cout<<"G4QIonIonElastic::PostStepDoIt: **** PostStepDoIt is done ****"<<G4endl; 00579 #endif 00580 return G4VDiscreteProcess::PostStepDoIt(track, step); 00581 }