#include <G4QDiffraction.hh>
Inheritance diagram for G4QDiffraction:
Public Member Functions | |
G4QDiffraction (const G4String &processName="CHIPS_DiffractiveInteraction") | |
~G4QDiffraction () | |
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 73 of file G4QDiffraction.hh.
G4QDiffraction::G4QDiffraction | ( | const G4String & | processName = "CHIPS_DiffractiveInteraction" |
) |
Definition at line 64 of file G4QDiffraction.cc.
References G4cout, G4endl, G4HadronicDeprecate, G4QCHIPSWorld::Get(), G4QCHIPSWorld::GetParticles(), G4VProcess::GetProcessName(), and G4VProcess::verboseLevel.
00064 : 00065 G4VDiscreteProcess(processName, fHadronic) 00066 { 00067 G4HadronicDeprecate("G4QDiffraction"); 00068 00069 #ifdef debug 00070 G4cout<<"G4QDiffraction::Constructor is called processName="<<processName<<G4endl; 00071 #endif 00072 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl; 00073 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max) 00074 }
G4QDiffraction::~G4QDiffraction | ( | ) |
G4LorentzVector G4QDiffraction::GetEnegryMomentumConservation | ( | ) |
G4double G4QDiffraction::GetMeanFreePath | ( | const G4Track & | aTrack, | |
G4double | previousStepSize, | |||
G4ForceCondition * | condition | |||
) | [virtual] |
Implements G4VDiscreteProcess.
Definition at line 87 of file G4QDiffraction.cc.
References DBL_MAX, G4cerr, G4cout, G4endl, G4QIsotope::Get(), G4QIsotope::GetCSVector(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterial(), G4QIsotope::GetMeanCrossSection(), G4Material::GetNumberOfElements(), G4DynamicParticle::GetTotalMomentum(), G4Material::GetVecNbOfAtomsPerVolume(), G4QIsotope::InitElement(), IsApplicable(), G4QIsotope::IsDefined(), G4Neutron::Neutron(), NotForced, and G4Proton::Proton().
Referenced by PostStepDoIt().
00088 { 00089 *F = NotForced; 00090 const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle(); 00091 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition(); 00092 if( !IsApplicable(*incidentParticleDefinition)) 00093 G4cout<<"-Warning-G4QDiffraction::GetMeanFreePath for notImplemented Particle"<<G4endl; 00094 // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material 00095 G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle 00096 #ifdef debug 00097 G4double KinEn = incidentParticle->GetKineticEnergy(); 00098 G4cout<<"G4QDiffraction::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl; 00099 #endif 00100 const G4Material* material = Track.GetMaterial(); // Get the current material 00101 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume(); 00102 const G4ElementVector* theElementVector = material->GetElementVector(); 00103 G4int nE=material->GetNumberOfElements(); 00104 #ifdef debug 00105 G4cout<<"G4QDiffraction::GetMeanFreePath:"<<nE<<" Elems in Material="<<*material<<G4endl; 00106 #endif 00107 G4int pPDG=0; 00108 // @@ At present it is made only for n & p, but can be extended if inXS are available 00109 if (incidentParticleDefinition == G4Proton::Proton() ) pPDG=2212; 00110 else if(incidentParticleDefinition == G4Neutron::Neutron()) pPDG=2112; 00111 else G4cout<<"G4QDiffraction::GetMeanFreePath: only nA & pA are implemented"<<G4endl; 00112 00113 G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton 00114 G4double sigma=0.; // Sums over elements for the material 00115 G4int IPIE=IsoProbInEl.size(); // How many old elements? 00116 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) 00117 { 00118 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector 00119 SPI->clear(); 00120 delete SPI; 00121 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector 00122 IsN->clear(); 00123 delete IsN; 00124 } 00125 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) 00126 ElementZ.clear(); // Clear the body vector for Z of Elements 00127 IsoProbInEl.clear(); // Clear the body vector for SPI 00128 ElIsoN.clear(); // Clear the body vector for N of Isotopes 00129 for(G4int i=0; i<nE; ++i) 00130 { 00131 G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element 00132 G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element 00133 ElementZ.push_back(Z); // Remember Z of the Element 00134 G4int isoSize=0; // The default for the isoVectorLength is 0 00135 G4int indEl=0; // Index of non-natural element or 0(default) 00136 G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect 00137 if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector 00138 #ifdef debug 00139 G4cout<<"G4QDiffraction::GetMeanFreePath: isovector Length="<<isoSize<<G4endl; 00140 #endif 00141 if(isoSize) // The Element has non-trivial abundance set 00142 { 00143 indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order 00144 #ifdef debug 00145 G4cout<<"G4QDiffr::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl; 00146 #endif 00147 if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define 00148 { 00149 std::vector<std::pair<G4int,G4double>*>* newAbund = 00150 new std::vector<std::pair<G4int,G4double>*>; 00151 G4double* abuVector=pElement->GetRelativeAbundanceVector(); 00152 for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes 00153 { 00154 G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z ! 00155 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z=" 00156 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl; 00157 G4double abund=abuVector[j]; 00158 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund); 00159 #ifdef debug 00160 G4cout<<"G4QDiffract::GetMeanFreePath:pair#"<<j<<",N="<<N<<",ab="<<abund<<G4endl; 00161 #endif 00162 newAbund->push_back(pr); 00163 } 00164 #ifdef debug 00165 G4cout<<"G4QDiffract::GetMeanFreePath:pairVectorLength="<<newAbund->size()<<G4endl; 00166 #endif 00167 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd 00168 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary 00169 delete newAbund; // Was "new" in the beginning of the name space 00170 } 00171 } 00172 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer 00173 std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector 00174 IsoProbInEl.push_back(SPI); 00175 std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector 00176 ElIsoN.push_back(IsN); 00177 G4int nIs=cs->size(); // A#Of Isotopes in the Element 00178 #ifdef debug 00179 G4cout<<"G4QDiffract::GetMFP:=***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl; 00180 #endif 00181 G4double susi=0.; // sum of CS over isotopes 00182 if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El 00183 { 00184 std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice 00185 G4int N=curIs->first; // #of Neuterons in the isotope j of El i 00186 IsN->push_back(N); // Remember Min N for the Element 00187 #ifdef debug 00188 G4cout<<"G4Q::GMFP:j="<<j<<",P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PD="<<pPDG<<G4endl; 00189 #endif 00190 G4double CSI=CalculateXS(Momentum, Z, N, pPDG); // XS(j,i) for theIsotope 00191 00192 #ifdef debug 00193 G4cout<<"G4QDiffraction::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom=" 00194 <<Momentu<<", XSec="<<CSI/millibarn<<G4endl; 00195 #endif 00196 curIs->second = CSI; 00197 susi+=CSI; // Make a sum per isotopes 00198 SPI->push_back(susi); // Remember summed cross-section 00199 } // End of temporary initialization of the cross sections in the G4QIsotope singeltone 00200 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV) 00201 #ifdef debug 00202 G4cout<<"G4QDiffraction::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl) 00203 <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl; 00204 #endif 00205 ElProbInMat.push_back(sigma); 00206 } // End of LOOP over Elements 00207 // Check that cross section is not zero and return the mean free path 00208 #ifdef debug 00209 G4cout<<"G4QDiffraction::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl; 00210 #endif 00211 if(sigma > 0.) return 1./sigma; // Mean path [distance] 00212 return DBL_MAX; 00213 }
G4int G4QDiffraction::GetNumberOfNeutronsInTarget | ( | ) |
G4bool G4QDiffraction::IsApplicable | ( | const G4ParticleDefinition & | particle | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 215 of file G4QDiffraction.cc.
References G4cout, G4endl, G4ParticleDefinition::GetPDGEncoding(), G4Neutron::Neutron(), and G4Proton::Proton().
Referenced by GetMeanFreePath(), and PostStepDoIt().
00216 { 00217 if (particle == *( G4Proton::Proton() )) return true; 00218 else if (particle == *( G4Neutron::Neutron() )) return true; 00219 //else if (particle == *( G4MuonMinus::MuonMinus() )) return true; 00220 //else if (particle == *( G4TauPlus::TauPlus() )) return true; 00221 //else if (particle == *( G4TauMinus::TauMinus() )) return true; 00222 //else if (particle == *( G4Electron::Electron() )) return true; 00223 //else if (particle == *( G4Positron::Positron() )) return true; 00224 //else if (particle == *( G4Gamma::Gamma() )) return true; 00225 //else if (particle == *( G4MuonPlus::MuonPlus() )) return true; 00226 //else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true; 00227 //else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true; 00228 //else if (particle == *( G4PionMinus::PionMinus() )) return true; 00229 //else if (particle == *( G4PionPlus::PionPlus() )) return true; 00230 //else if (particle == *( G4KaonPlus::KaonPlus() )) return true; 00231 //else if (particle == *( G4KaonMinus::KaonMinus() )) return true; 00232 //else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true; 00233 //else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true; 00234 //else if (particle == *( G4Lambda::Lambda() )) return true; 00235 //else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true; 00236 //else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true; 00237 //else if (particle == *( G4SigmaZero::SigmaZero() )) return true; 00238 //else if (particle == *( G4XiMinus::XiMinus() )) return true; 00239 //else if (particle == *( G4XiZero::XiZero() )) return true; 00240 //else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true; 00241 //else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true; 00242 //else if (particle == *( G4AntiProton::AntiProton() )) return true; 00243 #ifdef debug 00244 G4cout<<"***>>G4QDiffraction::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<G4endl; 00245 #endif 00246 return false; 00247 }
G4VParticleChange * G4QDiffraction::PostStepDoIt | ( | const G4Track & | aTrack, | |
const G4Step & | aStep | |||
) | [virtual] |
Reimplemented from G4VDiscreteProcess.
Definition at line 249 of file G4QDiffraction.cc.
References G4ParticleChange::AddSecondary(), G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4VProcess::aParticleChange, G4Electron::Electron(), G4ParticleTable::FindIon(), fStopAndKill, G4cerr, G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4QCHIPSWorld::Get(), G4QHadron::Get4Momentum(), G4DynamicParticle::Get4Momentum(), G4QHadron::GetBaryonNumber(), G4QHadron::GetCharge(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4Track::GetGlobalTime(), G4DynamicParticle::GetKineticEnergy(), G4QPDGCode::GetMass(), G4Track::GetMaterial(), GetMeanFreePath(), G4DynamicParticle::GetMomentumDirection(), G4Material::GetNumberOfElements(), G4QCHIPSWorld::GetParticles(), G4ParticleDefinition::GetParticleSubType(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetParticleType(), G4QHadron::GetPDGCode(), G4ParticleDefinition::GetPDGEncoding(), G4QDiffractionRatio::GetPointer(), G4Track::GetPosition(), G4QHadron::GetStrangeness(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetTouchableHandle(), G4Track::GetWeight(), G4Element::GetZ(), G4ParticleChange::Initialize(), IsApplicable(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), G4Neutron::Neutron(), NotForced, G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), position, G4Positron::Positron(), G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4Proton::Proton(), G4DynamicParticle::Set4Momentum(), G4DynamicParticle::SetDefinition(), G4Track::SetTouchableHandle(), G4Track::SetWeight(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4QDiffractionRatio::TargFragment(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().
00250 { 00251 static const G4double mProt= G4QPDGCode(2212).GetMass(); // CHIPS proton Mass in MeV 00252 static const G4double mNeut= G4QPDGCode(2112).GetMass(); // CHIPS neutron Mass in MeV 00253 static const G4double mPion= G4QPDGCode(111).GetMass(); // CHIPS Pi0 Mass in MeV 00254 static G4QDiffractionRatio* diffRatio; 00255 // 00256 //------------------------------------------------------------------------------------- 00257 static G4bool CWinit = true; // CHIPS Warld needs to be initted 00258 if(CWinit) 00259 { 00260 CWinit=false; 00261 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) 00262 diffRatio=G4QDiffractionRatio::GetPointer(); 00263 } 00264 //------------------------------------------------------------------------------------- 00265 const G4DynamicParticle* projHadron = track.GetDynamicParticle(); 00266 const G4ParticleDefinition* particle=projHadron->GetDefinition(); 00267 #ifdef debug 00268 G4cout<<"G4QDiffraction::PostStepDoIt: Before the GetMeanFreePath is called In4M=" 00269 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type=" 00270 <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl; 00271 #endif 00272 G4ForceCondition cond=NotForced; 00273 GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters? 00274 #ifdef debug 00275 G4cout<<"G4QDiffraction::PostStepDoIt: After GetMeanFreePath is called"<<G4endl; 00276 #endif 00277 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV! 00278 G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV 00279 G4double Momentum = proj4M.rho(); // @@ Just for the test purposes 00280 if(std::fabs(Momentum-momentum)>.000001) 00281 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl; 00282 #ifdef pdebug 00283 G4cout<<"G4QDiffraction::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum 00284 <<",proj4M="<<proj4M<<", projM="<<proj4M.m()<<G4endl; 00285 #endif 00286 if (!IsApplicable(*particle)) // Check applicability 00287 { 00288 G4cerr<<"G4QDiffraction::PostStepDoIt: Only NA is implemented."<<G4endl; 00289 return 0; 00290 } 00291 const G4Material* material = track.GetMaterial(); // Get the current material 00292 G4int Z=0; 00293 const G4ElementVector* theElementVector = material->GetElementVector(); 00294 G4int nE=material->GetNumberOfElements(); 00295 #ifdef debug 00296 G4cout<<"G4QDiffraction::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl; 00297 #endif 00298 G4int projPDG=0; // PDG Code prototype for the captured hadron 00299 // Not all these particles are implemented yet (see Is Applicable) 00300 if (particle == G4Proton::Proton() ) projPDG= 2212; 00301 else if (particle == G4Neutron::Neutron() ) projPDG= 2112; 00302 //else if (particle == G4PionMinus::PionMinus() ) projPDG= -211; 00303 //else if (particle == G4PionPlus::PionPlus() ) projPDG= 211; 00304 //else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321; 00305 //else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321; 00306 //else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130; 00307 //else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310; 00308 //else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13; 00309 //else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13; 00310 //else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14; 00311 //else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14; 00312 //else if (particle == G4Electron::Electron() ) projPDG= 11; 00313 //else if (particle == G4Positron::Positron() ) projPDG= -11; 00314 //else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12; 00315 //else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12; 00316 //else if (particle == G4Gamma::Gamma() ) projPDG= 22; 00317 //else if (particle == G4TauPlus::TauPlus() ) projPDG= -15; 00318 //else if (particle == G4TauMinus::TauMinus() ) projPDG= 15; 00319 //else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16; 00320 //else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16; 00321 //else if (particle == G4Lambda::Lambda() ) projPDG= 3122; 00322 //else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222; 00323 //else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112; 00324 //else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212; 00325 //else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312; 00326 //else if (particle == G4XiZero::XiZero() ) projPDG= 3322; 00327 //else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334; 00328 //else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112; 00329 //else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212; 00330 #ifdef debug 00331 G4int prPDG=particle->GetPDGEncoding(); 00332 G4cout<<"G4QDiffraction::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl; 00333 #endif 00334 if(!projPDG) 00335 { 00336 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl; 00337 return 0; 00338 } 00339 //G4double pM2=proj4M.m2(); // in MeV^2 00340 //G4double pM=std::sqrt(pM2); // in MeV 00341 G4double pM=mNeut; 00342 if(projPDG==2112) pM=mProt; 00343 // Element treatment 00344 G4int EPIM=ElProbInMat.size(); 00345 #ifdef debug 00346 G4cout<<"G4QDiffra::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl; 00347 #endif 00348 G4int i=0; 00349 if(EPIM>1) 00350 { 00351 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand(); 00352 for(i=0; i<nE; ++i) 00353 { 00354 #ifdef debug 00355 G4cout<<"G4QDiffra::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl; 00356 #endif 00357 if (rnd<ElProbInMat[i]) break; 00358 } 00359 if(i>=nE) i=nE-1; // Top limit for the Element 00360 } 00361 G4Element* pElement=(*theElementVector)[i]; 00362 Z=static_cast<G4int>(pElement->GetZ()); 00363 #ifdef debug 00364 G4cout<<"G4QDiffraction::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl; 00365 #endif 00366 if(Z<=0) 00367 { 00368 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Element with Z="<<Z<<G4endl; 00369 if(Z<0) return 0; 00370 } 00371 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes 00372 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i] 00373 G4int nofIsot=SPI->size(); // #of isotopes in the element i 00374 #ifdef debug 00375 G4cout<<"G4QDiffraction::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl; 00376 #endif 00377 G4int j=0; 00378 if(nofIsot>1) 00379 { 00380 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element 00381 for(j=0; j<nofIsot; ++j) 00382 { 00383 #ifdef debug 00384 G4cout<<"G4QDiffraction::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl; 00385 #endif 00386 if(rndI < (*SPI)[j]) break; 00387 } 00388 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope 00389 } 00390 G4int N =(*IsN)[j]; ; // Randomized number of neutrons 00391 #ifdef debug 00392 G4cout<<"G4QDiffraction::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<", MeV="<<MeV<<G4endl; 00393 #endif 00394 if(N<0) 00395 { 00396 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Isotope Z="<<Z<<" has 0>N="<<N<<G4endl; 00397 return 0; 00398 } 00399 nOfNeutrons=N; // Remember it for the energy-momentum check 00400 #ifdef debug 00401 G4cout<<"G4QDiffraction::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl; 00402 #endif 00403 if(N<0) 00404 { 00405 G4cerr<<"*Warning*G4QDiffraction::PostStepDoIt:Element with N="<<N<< G4endl; 00406 return 0; 00407 } 00408 aParticleChange.Initialize(track); 00409 #ifdef debug 00410 G4cout<<"G4QDiffraction::PostStepDoIt: track is initialized"<<G4endl; 00411 #endif 00412 G4double weight = track.GetWeight(); 00413 G4double localtime = track.GetGlobalTime(); 00414 G4ThreeVector position = track.GetPosition(); 00415 #ifdef debug 00416 G4cout<<"G4QDiffraction::PostStepDoIt: before Touchable extraction"<<G4endl; 00417 #endif 00418 G4TouchableHandle trTouchable = track.GetTouchableHandle(); 00419 #ifdef debug 00420 G4cout<<"G4QDiffraction::PostStepDoIt: Touchable is extracted"<<G4endl; 00421 #endif 00422 G4int targPDG=90000000+Z*1000+N; // CHIPS PDG Code of the target nucleus 00423 G4QPDGCode targQPDG(targPDG); // @@ use G4Ion and get rid of CHIPS World 00424 G4double tM=targQPDG.GetMass(); // CHIPS final nucleus mass in MeV 00425 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?) 00426 G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector 00427 G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction 00428 #ifdef debug 00429 G4cout<<"G4QDiffraction::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl; 00430 #endif 00431 EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation 00432 // @@ Probably this is not necessary any more 00433 #ifdef debug 00434 G4cout<<"G4QDiff::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl; 00435 #endif 00436 G4double xSec=CalculateXS(Momentum, Z, N, projPDG); // Recalculate CrossSection 00437 #ifdef debug 00438 G4cout<<"G4QDiffra::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",XS="<<xSec/millibarn<<G4endl; 00439 #endif 00440 #ifdef nandebug 00441 if(xSec>0. || xSec<0. || xSec==0); 00442 else G4cout<<"-Warning-G4QDiffraction::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl; 00443 #endif 00444 // @@ check a possibility to separate p, n, or alpha (!) 00445 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing 00446 { 00447 #ifdef pdebug 00448 G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Zero cross-section* PDG="<<projPDG 00449 <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl; 00450 #endif 00451 //Do Nothing Action insead of the reaction 00452 aParticleChange.ProposeEnergy(kinEnergy); 00453 aParticleChange.ProposeLocalEnergyDeposit(0.); 00454 aParticleChange.ProposeMomentumDirection(dir) ; 00455 return G4VDiscreteProcess::PostStepDoIt(track,step); 00456 } 00457 G4double totCMMass=tot4M.m(); // Total CM mass, pM=projectileMass, tM=targetMass 00458 if(totCMMass < mPion+pM+tM) // The diffraction reaction is impossible -> Do Nothing 00459 { 00460 #ifdef pdebug 00461 G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Below Diffraction Threshold* cmM="<<totCMMass 00462 <<">pM="<<pM<<"+tM="<<tM<<"+pi0="<<mPion<<"=="<<pM+tM+mPion<<G4endl; 00463 #endif 00464 //Do Nothing Action insead of the reaction 00465 aParticleChange.ProposeEnergy(kinEnergy); 00466 aParticleChange.ProposeLocalEnergyDeposit(0.); 00467 aParticleChange.ProposeMomentumDirection(dir) ; 00468 return G4VDiscreteProcess::PostStepDoIt(track,step); 00469 } 00470 // Kill interacting hadron 00471 aParticleChange.ProposeTrackStatus(fStopAndKill); 00472 G4QHadronVector* out=diffRatio->TargFragment(projPDG, proj4M, Z, N); 00473 G4int nSec=out->size(); // #of secondaries in the diffraction reaction 00474 G4DynamicParticle* theSec=0; // A prototype for secondary for the secondary 00475 G4LorentzVector dif4M(0.,0.,0.,0.); // Prototype for the secondary 4-momentum 00476 G4int difPDG=0; // PDG code of the secondary 00477 G4QHadron* difQH=0; // Prototype for a Q-secondary 00478 #ifdef pdebug 00479 G4cout<<"G4QDiffraction::PostStepDoIt: =---=found=---= nSecondaries="<<nSec<<G4endl; 00480 #endif 00481 for(i=0; i<nSec; i++) 00482 { 00483 difQH = (*out)[i]; 00484 difPDG= difQH->GetPDGCode(); 00485 G4ParticleDefinition* theDefinition=0; 00486 if (difPDG==2212 || difPDG==90001000) theDefinition=G4Proton::Proton(); 00487 else if(difPDG==2112 || difPDG==90000001) theDefinition=G4Neutron::Neutron(); 00488 else if(difPDG== 22) theDefinition=G4Gamma::Gamma(); 00489 else if(difPDG== 111) theDefinition=G4PionZero::PionZero(); 00490 else if(difPDG==-211 || difPDG==89999001) theDefinition=G4PionMinus::PionMinus(); 00491 else if(difPDG== 211 || difPDG==90000999) theDefinition=G4PionPlus::PionPlus(); 00492 else if(difPDG== 321 || difPDG==89001000) theDefinition=G4KaonPlus::KaonPlus(); 00493 else if(difPDG==-321 || difPDG==90999000) theDefinition=G4KaonMinus::KaonMinus(); 00494 else if(difPDG== 130 || difPDG==-311 || difPDG==89000001) 00495 theDefinition=G4KaonZeroLong::KaonZeroLong(); 00496 else if(difPDG== 310 || difPDG== 311 || difPDG==90999999) 00497 theDefinition=G4KaonZeroShort::KaonZeroShort(); 00498 else if(difPDG==3122 || difPDG==91000000) theDefinition=G4Lambda::Lambda(); 00499 else if(difPDG== 3222) theDefinition=G4SigmaPlus::SigmaPlus(); 00500 else if(difPDG== 3112) theDefinition=G4SigmaMinus::SigmaMinus(); 00501 else if(difPDG== 3212) theDefinition=G4SigmaZero::SigmaZero(); 00502 else if(difPDG== 3312) theDefinition=G4XiMinus::XiMinus(); 00503 else if(difPDG== 3322) theDefinition=G4XiZero::XiZero(); 00504 else if(difPDG== 3334) theDefinition=G4OmegaMinus::OmegaMinus(); 00505 else if(difPDG==-2112) theDefinition=G4AntiNeutron::AntiNeutron(); 00506 else if(difPDG==-2212) theDefinition=G4AntiProton::AntiProton(); 00507 else if(difPDG==-3122) theDefinition=G4AntiLambda::AntiLambda(); 00508 else if(difPDG==-3222) theDefinition=G4AntiSigmaPlus::AntiSigmaPlus(); 00509 else if(difPDG==-3112) theDefinition=G4AntiSigmaMinus::AntiSigmaMinus(); 00510 else if(difPDG==-3212) theDefinition=G4AntiSigmaZero::AntiSigmaZero(); 00511 else if(difPDG==-3312) theDefinition=G4AntiXiMinus::AntiXiMinus(); 00512 else if(difPDG==-3322) theDefinition=G4AntiXiZero::AntiXiZero(); 00513 else if(difPDG==-3334) theDefinition=G4AntiOmegaMinus::AntiOmegaMinus(); 00514 else if(difPDG== -11) theDefinition=G4Electron::Electron(); 00515 else if(difPDG== -13) theDefinition=G4MuonMinus::MuonMinus(); 00516 else if(difPDG== 11) theDefinition=G4Positron::Positron(); 00517 else if(difPDG== 13) theDefinition=G4MuonPlus::MuonPlus(); 00518 else 00519 { 00520 Z = difQH->GetCharge(); 00521 G4int B = difQH->GetBaryonNumber(); 00522 G4int S = difQH->GetStrangeness(); 00523 if(S||Z>B||Z<0)G4cout<<"-Warning-G4QDif::PoStDoIt:Z="<<Z<<",A="<<B<<",S="<<S<<G4endl; 00524 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(Z,B,0,0); 00525 #ifdef pdebug 00526 G4cout<<"G4QDiffraction::PoStDoIt:Ion,Z="<<Z<<",A="<<B<<",D="<<theDefinition<<G4endl; 00527 #endif 00528 } 00529 if(theDefinition) 00530 { 00531 theSec = new G4DynamicParticle; // A secondary for the recoil hadron 00532 theSec->SetDefinition(theDefinition); 00533 dif4M = difQH->Get4Momentum(); 00534 EnMomConservation-=dif4M; 00535 theSec->Set4Momentum(dif4M); 00536 G4Track* aNewTrack = new G4Track(theSec, localtime, position ); 00537 aNewTrack->SetWeight(weight); // weighted 00538 aNewTrack->SetTouchableHandle(trTouchable); 00539 aParticleChange.AddSecondary( aNewTrack ); 00540 #ifdef pdebug 00541 G4cout<<"G4QDiffraction::PostStepDoIt: Filled 4M="<<dif4M<<", PDG="<<difPDG<<G4endl; 00542 #endif 00543 } 00544 else G4cout<<"-Warning-G4QDif::PSDI: Lost PDG="<<difPDG<<", Z="<<difQH->GetCharge() 00545 <<", A="<<difQH->GetBaryonNumber()<<",S ="<<difQH->GetStrangeness()<<G4endl; 00546 delete difQH; // Clean up the output QHadrons 00547 } 00548 delete out; // Delete the output QHadron-vector 00549 #ifdef debug 00550 G4cout<<"G4QDiffraction::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl; 00551 #endif 00552 return G4VDiscreteProcess::PostStepDoIt(track, step); 00553 }