#include <G4QLowEnergy.hh>
Inheritance diagram for G4QLowEnergy:
Public Member Functions | |
G4QLowEnergy (const G4String &processName="CHIPS_LowEnergyIonIonInelastic") | |
~G4QLowEnergy () | |
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 () |
void | SwitchOnEvaporation () |
void | SwitchOffEvaporation () |
Definition at line 72 of file G4QLowEnergy.hh.
G4QLowEnergy::G4QLowEnergy | ( | const G4String & | processName = "CHIPS_LowEnergyIonIonInelastic" |
) |
Definition at line 60 of file G4QLowEnergy.cc.
References G4cout, G4endl, G4HadronicDeprecate, G4QCHIPSWorld::Get(), G4QCHIPSWorld::GetParticles(), G4VProcess::GetProcessName(), and G4VProcess::verboseLevel.
00060 : 00061 G4VDiscreteProcess(processName, fHadronic), evaporate(true) 00062 { 00063 G4HadronicDeprecate("G4QLowEnergy"); 00064 00065 #ifdef debug 00066 G4cout<<"G4QLowEnergy::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 }
G4QLowEnergy::~G4QLowEnergy | ( | ) |
G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation | ( | ) |
G4double G4QLowEnergy::GetMeanFreePath | ( | const G4Track & | aTrack, | |
G4double | previousStepSize, | |||
G4ForceCondition * | condition | |||
) | [virtual] |
Implements G4VDiscreteProcess.
Definition at line 83 of file G4QLowEnergy.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(), G4QProtonNuclearCrossSection::GetPointer(), G4QIonIonCrossSection::GetPointer(), G4DynamicParticle::GetTotalMomentum(), G4Material::GetVecNbOfAtomsPerVolume(), G4He3::He3(), G4QIsotope::InitElement(), IsApplicable(), G4QIsotope::IsDefined(), NotForced, G4Proton::Proton(), and G4Triton::Triton().
00084 { 00085 *F = NotForced; 00086 const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle(); 00087 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition(); 00088 if( !IsApplicable(*incidentParticleDefinition)) 00089 G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath for notImplemented Particle"<<G4endl; 00090 // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material 00091 G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle 00092 #ifdef debug 00093 G4double KinEn = incidentParticle->GetKineticEnergy(); 00094 G4cout<<"G4QLowEnergy::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl; 00095 #endif 00096 const G4Material* material = Track.GetMaterial(); // Get the current material 00097 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume(); 00098 const G4ElementVector* theElementVector = material->GetElementVector(); 00099 G4int nE=material->GetNumberOfElements(); 00100 #ifdef debug 00101 G4cout<<"G4QLowEnergy::GetMeanFreePath:"<<nE<<" Elems"<<G4endl; 00102 #endif 00103 G4int pPDG=0; 00104 G4int Z=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge()); 00105 G4int A=incidentParticleDefinition->GetBaryonNumber(); 00106 if ( incidentParticleDefinition == G4Proton::Proton() ) pPDG = 2212; 00107 else if ( incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020; 00108 else if ( incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 1000020040; 00109 else if ( incidentParticleDefinition == G4Triton::Triton() ) pPDG = 1000010030; 00110 else if ( incidentParticleDefinition == G4He3::He3() ) pPDG = 1000020030; 00111 else if ( incidentParticleDefinition == G4GenericIon::GenericIon() || (Z > 0 && A > 0)) 00112 { 00113 pPDG=incidentParticleDefinition->GetPDGEncoding(); 00114 #ifdef debug 00115 G4int prPDG=1000000000+10000*A+10*Z; 00116 G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl; 00117 #endif 00118 } 00119 else G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath: only AA & pA implemented"<<G4endl; 00120 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer(); 00121 if(pPDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer(); // just to test 00122 Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A 00123 G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton 00124 G4double sigma=0.; // Sums over elements for the material 00125 G4int IPIE=IsoProbInEl.size(); // How many old elements? 00126 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) 00127 { 00128 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector 00129 SPI->clear(); 00130 delete SPI; 00131 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector 00132 IsN->clear(); 00133 delete IsN; 00134 } 00135 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) 00136 ElementZ.clear(); // Clear the body vector for Z of Elements 00137 IsoProbInEl.clear(); // Clear the body vector for SPI 00138 ElIsoN.clear(); // Clear the body vector for N of Isotopes 00139 for(G4int i=0; i<nE; ++i) 00140 { 00141 G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element 00142 Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element 00143 ElementZ.push_back(Z); // Remember Z of the Element 00144 G4int isoSize=0; // The default for the isoVectorLength is 0 00145 G4int indEl=0; // Index of non-natural element or 0(default) 00146 G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect 00147 if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector 00148 #ifdef debug 00149 G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl; 00150 #endif 00151 if(isoSize) // The Element has non-trivial abundance set 00152 { 00153 indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order 00154 #ifdef debug 00155 G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl; 00156 #endif 00157 if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define 00158 { 00159 std::vector<std::pair<G4int,G4double>*>* newAbund = 00160 new std::vector<std::pair<G4int,G4double>*>; 00161 G4double* abuVector=pElement->GetRelativeAbundanceVector(); 00162 for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes 00163 { 00164 G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z ! 00165 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z=" 00166 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl; 00167 G4double abund=abuVector[j]; 00168 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund); 00169 #ifdef debug 00170 G4cout<<"G4QLowEnergy::GetMeanFreePath:pair#"<<j<<",N="<<N<<",a="<<abund<<G4endl; 00171 #endif 00172 newAbund->push_back(pr); 00173 } 00174 #ifdef debug 00175 G4cout<<"G4QLowEnergy::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl; 00176 #endif 00177 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd 00178 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary 00179 delete newAbund; // Was "new" in the beginning of the name space 00180 } 00181 } 00182 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer 00183 std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector 00184 IsoProbInEl.push_back(SPI); 00185 std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector 00186 ElIsoN.push_back(IsN); 00187 G4int nIs=cs->size(); // A#Of Isotopes in the Element 00188 #ifdef debug 00189 G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl; 00190 #endif 00191 G4double susi=0.; // sum of CS over isotopes 00192 if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El 00193 { 00194 std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice 00195 G4int N=curIs->first; // #of Neuterons in the isotope j of El i 00196 IsN->push_back(N); // Remember Min N for the Element 00197 #ifdef debug 00198 G4cout<<"G4QLoE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl; 00199 #endif 00200 G4bool ccsf=true; // Extract inelastic Ion-Ion cross-section 00201 #ifdef debug 00202 G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl; 00203 #endif 00204 G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope 00205 #ifdef debug 00206 G4cout<<"G4QLowEnergy::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom=" 00207 <<Momentum<<", XSec="<<CSI/millibarn<<G4endl; 00208 #endif 00209 curIs->second = CSI; 00210 susi+=CSI; // Make a sum per isotopes 00211 SPI->push_back(susi); // Remember summed cross-section 00212 } // End of temporary initialization of the cross sections in the G4QIsotope singeltone 00213 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV) 00214 #ifdef debug 00215 G4cout<<"G4QLowEnergy::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl) 00216 <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl; 00217 #endif 00218 ElProbInMat.push_back(sigma); 00219 } // End of LOOP over Elements 00220 // Check that cross section is not zero and return the mean free path 00221 #ifdef debug 00222 G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl; 00223 #endif 00224 if(sigma > 0.000000001) return 1./sigma; // Mean path [distance] 00225 return DBL_MAX; 00226 }
G4int G4QLowEnergy::GetNumberOfNeutronsInTarget | ( | ) |
G4bool G4QLowEnergy::IsApplicable | ( | const G4ParticleDefinition & | particle | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 228 of file G4QLowEnergy.cc.
References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4cout, G4endl, G4GenericIon::GenericIon(), G4ParticleDefinition::GetBaryonNumber(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGEncoding(), G4He3::He3(), G4Neutron::Neutron(), G4Proton::Proton(), and G4Triton::Triton().
Referenced by GetMeanFreePath(), and PostStepDoIt().
00229 { 00230 G4int Z=static_cast<G4int>(particle.GetPDGCharge()); 00231 G4int A=particle.GetBaryonNumber(); 00232 if (particle == *( G4Proton::Proton() )) return true; 00233 else if (particle == *( G4Neutron::Neutron() )) return true; 00234 else if (particle == *( G4Deuteron::Deuteron() )) return true; 00235 else if (particle == *( G4Alpha::Alpha() )) return true; 00236 else if (particle == *( G4Triton::Triton() )) return true; 00237 else if (particle == *( G4He3::He3() )) return true; 00238 else if (particle == *( G4GenericIon::GenericIon() )) return true; 00239 else if (Z > 0 && A > 0) return true; 00240 #ifdef debug 00241 G4cout<<"***>>G4QLowEnergy::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<", A=" 00242 <<A<<", Z="<<Z<<G4endl; 00243 #endif 00244 return false; 00245 }
G4VParticleChange * G4QLowEnergy::PostStepDoIt | ( | const G4Track & | aTrack, | |
const G4Step & | aStep | |||
) | [virtual] |
Reimplemented from G4VDiscreteProcess.
Definition at line 247 of file G4QLowEnergy.cc.
References G4ParticleChange::AddSecondary(), G4Alpha::Alpha(), G4VProcess::aParticleChange, G4QNucleus::CoulombBarGen(), G4Deuteron::Deuteron(), G4QNucleus::EvaporateNucleus(), FatalException, G4ParticleTable::FindIon(), fStopAndKill, G4cerr, G4cout, G4endl, G4Exception(), G4RandomDirection(), G4UniformRand, G4Gamma::Gamma(), G4GenericIon::GenericIon(), G4QCHIPSWorld::Get(), G4DynamicParticle::Get4Momentum(), G4ParticleDefinition::GetAntiQuarkContent(), G4ParticleDefinition::GetBaryonNumber(), G4VQCrossSection::GetCrossSection(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4VQCrossSection::GetExchangeT(), G4Track::GetGlobalTime(), G4VQCrossSection::GetHMaxT(), G4DynamicParticle::GetKineticEnergy(), G4QPDGCode::GetMass(), G4Track::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4QPDGCode::GetNuclMass(), G4Material::GetNumberOfElements(), G4QCHIPSWorld::GetParticles(), G4ParticleDefinition::GetParticleSubType(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetParticleType(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGEncoding(), G4QNeutronElasticCrossSection::GetPointer(), G4QProtonElasticCrossSection::GetPointer(), G4Track::GetPosition(), G4ParticleDefinition::GetQuarkContent(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetTouchableHandle(), G4Track::GetWeight(), G4Element::GetZ(), G4He3::He3(), G4ParticleChange::Initialize(), IsApplicable(), JustWarning, G4Lambda::Lambda(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), position, G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4Proton::Proton(), G4DynamicParticle::Set4Momentum(), G4DynamicParticle::SetDefinition(), G4VParticleChange::SetNumberOfSecondaries(), G4Track::SetTouchableHandle(), G4Track::SetWeight(), and G4Triton::Triton().
00248 { 00249 static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV; // CHIPS proton Mass in MeV 00250 static const G4double mPro2= mProt*mProt; // CHIPS sq proton Mass 00251 static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV; // CHIPS neutron Mass in MeV 00252 static const G4double mNeu2= mNeut*mNeut; // CHIPS sq neutron Mass 00253 static const G4double mLamb= G4QPDGCode(3122).GetMass()/MeV; // CHIPS Lambda Mass in MeV 00254 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0)/MeV; 00255 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0)/MeV; 00256 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0)/MeV; 00257 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0)/MeV; 00258 static const G4double mFm= 250*MeV; 00259 static const G4double third= 1./3.; 00260 static const G4ThreeVector zeroMom(0.,0.,0.); 00261 static G4ParticleDefinition* aGamma = G4Gamma::Gamma(); 00262 static G4ParticleDefinition* aPiZero = G4PionZero::PionZero(); 00263 static G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus(); 00264 static G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus(); 00265 static G4ParticleDefinition* aProton = G4Proton::Proton(); 00266 static G4ParticleDefinition* aNeutron = G4Neutron::Neutron(); 00267 static G4ParticleDefinition* aLambda = G4Lambda::Lambda(); 00268 static G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron(); 00269 static G4ParticleDefinition* aTriton = G4Triton::Triton(); 00270 static G4ParticleDefinition* aHe3 = G4He3::He3(); 00271 static G4ParticleDefinition* anAlpha = G4Alpha::Alpha(); 00272 static const G4int nCh=26; // #of combinations 00273 static G4QNucleus Nuc; // A fake nucleus to call Evaporation 00274 // 00275 //------------------------------------------------------------------------------------- 00276 static G4bool CWinit = true; // CHIPS Warld needs to be initted 00277 if(CWinit) 00278 { 00279 CWinit=false; 00280 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) 00281 } 00282 //------------------------------------------------------------------------------------- 00283 const G4DynamicParticle* projHadron = track.GetDynamicParticle(); 00284 const G4ParticleDefinition* particle=projHadron->GetDefinition(); 00285 #ifdef pdebug 00286 G4cout<<"G4QLowEnergy::PostStepDoIt: *** Called *** In4M=" 00287 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type=" 00288 <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl; 00289 #endif 00290 //G4ForceCondition cond=NotForced; 00291 //GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters? 00292 #ifdef debug 00293 G4cout<<"G4QLowEnergy::PostStepDoIt: After GetMeanFreePath is called"<<G4endl; 00294 #endif 00295 std::vector<G4Track*> result; 00296 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV! 00297 G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV 00298 G4double Momentum = proj4M.rho(); // @@ Just for the test purposes 00299 if(std::fabs(Momentum-momentum)>.000001) 00300 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl; 00301 #ifdef debug 00302 G4cout<<"G4QLowEnergy::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",proj4M,m=" 00303 <<proj4M<<proj4M.m()<<G4endl; 00304 #endif 00305 if (!IsApplicable(*particle)) // Check applicability 00306 { 00307 G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl; 00308 return 0; 00309 } 00310 const G4Material* material = track.GetMaterial(); // Get the current material 00311 const G4ElementVector* theElementVector = material->GetElementVector(); 00312 G4int nE=material->GetNumberOfElements(); 00313 #ifdef debug 00314 G4cout<<"G4QLowEnergy::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl; 00315 #endif 00316 G4int projPDG=0; // PDG Code prototype for the captured hadron 00317 // Not all these particles are implemented yet (see Is Applicable) 00318 G4int Z=static_cast<G4int>(particle->GetPDGCharge()); 00319 G4int A=particle->GetBaryonNumber(); 00320 if (particle == G4Proton::Proton() ) projPDG= 2212; 00321 else if (particle == G4Neutron::Neutron() ) projPDG= 2112; 00322 else if (particle == G4Deuteron::Deuteron() ) projPDG= 1000010020; 00323 else if (particle == G4Alpha::Alpha() ) projPDG= 1000020040; 00324 else if (particle == G4Triton::Triton() ) projPDG= 1000010030; 00325 else if (particle == G4He3::He3() ) projPDG= 1000020030; 00326 else if (particle == G4GenericIon::GenericIon() || (Z > 0 && A > 0)) 00327 { 00328 projPDG=particle->GetPDGEncoding(); 00329 #ifdef debug 00330 G4int prPDG=1000000000+10000*Z+10*A; // just for the testing purposes 00331 G4cout<<"G4QLowEnergy::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl; 00332 #endif 00333 } 00334 else G4cout<<"-Warning-G4QLowEnergy::PostStepDoIt:Unknown projectile Ion"<<G4endl; 00335 #ifdef pdebug 00336 G4int prPDG=particle->GetPDGEncoding(); 00337 G4cout<<"G4QLowEnergy::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl; 00338 #endif 00339 if(!projPDG) 00340 { 00341 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl; 00342 return 0; 00343 } 00344 // Element treatment 00345 G4int EPIM=ElProbInMat.size(); 00346 #ifdef debug 00347 G4cout<<"G4QLowEn::PostStDoIt: m="<<EPIM<<", n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl; 00348 #endif 00349 G4int i=0; 00350 if(EPIM>1) 00351 { 00352 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand(); 00353 for(i=0; i<nE; ++i) 00354 { 00355 #ifdef debug 00356 G4cout<<"G4QLowEn::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl; 00357 #endif 00358 if (rnd<ElProbInMat[i]) break; 00359 } 00360 if(i>=nE) i=nE-1; // Top limit for the Element 00361 } 00362 G4Element* pElement=(*theElementVector)[i]; 00363 G4int tZ=static_cast<G4int>(pElement->GetZ()); 00364 #ifdef debug 00365 G4cout<<"G4QLowEnergy::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl; 00366 #endif 00367 if(tZ<=0) 00368 { 00369 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Element with Z="<<tZ<<G4endl; 00370 if(tZ<0) return 0; 00371 } 00372 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes 00373 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i] 00374 G4int nofIsot=SPI->size(); // #of isotopes in the element i 00375 #ifdef debug 00376 G4cout<<"G4QLowEnergy::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl; 00377 #endif 00378 G4int j=0; 00379 if(nofIsot>1) 00380 { 00381 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element 00382 for(j=0; j<nofIsot; ++j) 00383 { 00384 #ifdef debug 00385 G4cout<<"G4QLowEnergy::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl; 00386 #endif 00387 if(rndI < (*SPI)[j]) break; 00388 } 00389 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope 00390 } 00391 G4int tN =(*IsN)[j]; ; // Randomized number of neutrons 00392 #ifdef debug 00393 G4cout<<"G4QLowEnergy::PostStepDoIt: j="<<i<<", N(isotope)="<<tN<<", MeV="<<MeV<<G4endl; 00394 #endif 00395 if(tN<0) 00396 { 00397 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Isotope Z="<<tZ<<" has 0>N="<<tN<<G4endl; 00398 return 0; 00399 } 00400 nOfNeutrons=tN; // Remember it for the energy-momentum check 00401 #ifdef debug 00402 G4cout<<"G4QLowEnergy::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl; 00403 #endif 00404 if(tN<0) 00405 { 00406 G4cerr<<"***Warning***G4QLowEnergy::PostStepDoIt: Element with N="<<tN<< G4endl; 00407 return 0; 00408 } 00409 aParticleChange.Initialize(track); 00410 #ifdef debug 00411 G4cout<<"G4QLowEnergy::PostStepDoIt: track is initialized"<<G4endl; 00412 #endif 00413 G4double weight = track.GetWeight(); 00414 G4double localtime = track.GetGlobalTime(); 00415 G4ThreeVector position = track.GetPosition(); 00416 #ifdef debug 00417 G4cout<<"G4QLowEnergy::PostStepDoIt: before Touchable extraction"<<G4endl; 00418 #endif 00419 G4TouchableHandle trTouchable = track.GetTouchableHandle(); 00420 #ifdef debug 00421 G4cout<<"G4QLowEnergy::PostStepDoIt: Touchable is extracted"<<G4endl; 00422 #endif 00423 G4int targPDG=90000000+tZ*1000+tN; // Target PDG code 00424 G4QPDGCode targQPDG(targPDG); // To get CHIPS properties 00425 G4double tM=targQPDG.GetMass(); // CHIPS target nucleus mass in MeV 00426 G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3); // Strangeness 00427 G4int pZ=static_cast<G4int>(particle->GetPDGCharge()); // Charge of the projectile 00428 G4int pN=particle->GetBaryonNumber()-pZ-pL;// #of neutrons in projectile 00429 G4double pM=targQPDG.GetNuclMass(pZ,pN,0); // CHIPS projectile nucleus mass in MeV 00430 G4double cosp=-14*Momentum*(tM-pM)/tM/pM; // Asymmetry power for angular distribution 00431 #ifdef debug 00432 G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ 00433 <<","<<tN<<"), cosp="<<cosp<<G4endl; 00434 #endif 00435 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?) 00436 G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector 00437 G4LorentzVector targ4M=G4LorentzVector(0.,0.,0.,tM); // Target's 4-mom 00438 G4LorentzVector tot4M =proj4M+targ4M; // Total 4-mom of the reaction 00439 #ifdef pdebug 00440 G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl; 00441 #endif 00442 EnMomConservation=tot4M; // Total 4mom of reaction for E/M conservation 00443 // --- Start of Coulomb barrier check --- 00444 G4double dER = tot4M.e() - tM - pM; // Energy of the reaction 00445 G4double pA = pZ+pN; // Atomic weight of the projectile (chged blw) 00446 G4double tA = tZ+tN; // Atomic weight of the target (changed below) 00447 G4double CBE = Nuc.CoulombBarGen(tZ, tA, pZ, pA); // Coulomb Barrier of the reaction 00448 #ifdef debug 00449 G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ 00450 <<","<<tN<<"), dE="<<dER<<", CB="<<CBE<<G4endl; 00451 #endif 00452 if(dER<CBE) // The cross-section iz 0 -> Do Nothing 00453 { 00454 #ifdef debug 00455 G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Below Coulomb barrier* PDG="<<projPDG 00456 <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl; 00457 #endif 00458 //Do Nothing Action insead of the reaction 00459 aParticleChange.ProposeEnergy(kinEnergy); 00460 aParticleChange.ProposeLocalEnergyDeposit(0.); 00461 aParticleChange.ProposeMomentumDirection(dir) ; 00462 return G4VDiscreteProcess::PostStepDoIt(track,step); 00463 } 00464 // --- End of Coulomb barrier check --- 00465 // @@ Probably this is not necessary any more 00466 #ifdef debug 00467 G4cout<<"G4QLE::PSDI:false,P="<<Momentum<<",Z="<<pZ<<",N="<<pN<<",PDG="<<projPDG<<G4endl; 00468 #endif 00469 G4double xSec=CalculateXS(Momentum, tZ, tN, projPDG); // Recalculate CrossSection 00470 #ifdef pdebug 00471 G4cout<<"G4QLowEn::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",tZ="<<tZ<<",N="<<tN<<",XS=" 00472 <<xSec/millibarn<<G4endl; 00473 #endif 00474 #ifdef nandebug 00475 if(xSec>0. || xSec<0. || xSec==0); 00476 else G4cout<<"-Warning-G4QLowEnergy::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl; 00477 #endif 00478 // @@ check a possibility to separate p, n, or alpha (!) 00479 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing 00480 { 00481 #ifdef debug 00482 G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Zero cross-section* PDG="<<projPDG 00483 <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl; 00484 #endif 00485 //Do Nothing Action insead of the reaction 00486 aParticleChange.ProposeEnergy(kinEnergy); 00487 aParticleChange.ProposeLocalEnergyDeposit(0.); 00488 aParticleChange.ProposeMomentumDirection(dir) ; 00489 return G4VDiscreteProcess::PostStepDoIt(track,step); 00490 } 00491 // Kill interacting hadron 00492 aParticleChange.ProposeEnergy(0.); 00493 aParticleChange.ProposeTrackStatus(fStopAndKill); 00494 G4int tB=tZ+tN; 00495 G4int pB=pZ+pN; 00496 #ifdef pdebug 00497 G4cout<<"G4QLowEn::PSDI: Projectile track is killed"<<", tA="<<tB<<", pA="<<pB<<G4endl; 00498 #endif 00499 // algorithm implementation STARTS HERE --- All calculations are in IU -------- 00500 tA=tB; 00501 pA=pB; 00502 G4double tR=1.1; // target nucleus R in fm 00503 if(tB > 1) tR*=std::pow(tA,third); // in fm 00504 G4double pR=1.1; // projectile nucleus R in fm 00505 if(pB > 1) pR*=std::pow(pA,third); // in fm 00506 G4double R=tR+pR; // total radius 00507 G4double R2=R*R; 00508 G4int tD=0; 00509 G4int pD=0; 00510 G4int nAt=0; 00511 G4int nAtM=27; 00512 G4int nSec = 0; 00513 G4double tcM=0.; 00514 G4double tnM=1.; 00515 #ifdef edebug 00516 G4int totChg=0; 00517 G4int totBaN=0; 00518 G4LorentzVector tch4M =tot4M; // Total 4-mom of the reaction 00519 #endif 00520 while((!tD && !pD && ++nAt<nAtM) || tcM<tnM) 00521 { 00522 #ifdef edebug 00523 totChg=tZ+pZ; 00524 totBaN=tB+pB; 00525 tch4M =tot4M; 00526 G4cout<<">G4QLEn::PSDI:#"<<nAt<<",tC="<<totChg<<",tA="<<totBaN<<",t4M="<<tch4M<<G4endl; 00527 #endif 00528 G4LorentzVector tt4M=tot4M; 00529 G4int resultSize=result.size(); 00530 if(resultSize) 00531 { 00532 for(i=0; i<resultSize; ++i) delete result[i]; 00533 result.clear(); 00534 } 00535 G4double D=std::sqrt(R2*G4UniformRand()); 00536 #ifdef pdebug 00537 G4cout<<"G4QLowEn::PSDI: D="<<D<<", tR="<<tR<<", pR="<<pR<<G4endl; 00538 #endif 00539 if(D > std::fabs(tR-pR)) // leading parts can be separated 00540 { 00541 nSec = 0; 00542 G4double DtR=D-tR; 00543 G4double DpR=D-pR; 00544 G4double DtR2=DtR*DtR; 00545 G4double DpR2=DpR*DpR; 00546 //G4double DtR3=DtR2*DtR; 00547 G4double DpR3=DpR2*DpR; 00548 //G4double DtR4=DtR3*DtR; 00549 G4double DpR4=DpR3*DpR; 00550 G4double tR2=tR*tR; 00551 G4double pR2=pR*pR; 00552 G4double tR3=tR2*tR; 00553 //G4double pR3=pR2*pR; 00554 G4double tR4=tR3*tR; 00555 //G4double pR4=pR3*pR; 00556 G4double HD=16.*D; 00557 G4double tF=tA*(6*tR2*(pR2-DtR2)-4*D*(tR3-DpR3)+3*(tR4-DpR4))/HD/tR3; 00558 G4double pF=tF; 00559 tD=static_cast<G4int>(tF); 00560 pD=static_cast<G4int>(pF); 00561 //if(G4UniformRand() < tF-tD) ++tD; // Simple solution 00562 //if(G4UniformRand() < pF-pD) ++pD; 00563 // Enhance alphas solution 00564 if(std::fabs(tF-4.) < 1.) tD=4; // @@ 1. is the enhancement parameter 00565 else if(G4UniformRand() < tF-tD) ++tD; 00566 if(std::fabs(pF-4.) < 1.) pD=4; 00567 else if(G4UniformRand() < pF-pD) ++pD; 00568 if(tD > tB) tD=tB; 00569 if(pD > pB) pD=tB; 00570 // @@ Quasi Free is not debugged @@ The following close it 00571 if(tD < 1) tD=1; 00572 if(pD < 1) pD=1; 00573 #ifdef pdebug 00574 G4cout<<"G4QLE::PSDI:#"<<nAt<<",pF="<<pF<<",tF="<<tF<<",pD="<<pD<<",tD="<<tD<<G4endl; 00575 #endif 00576 G4int pC=0; // charge of the projectile fraction 00577 G4int tC=0; // charge of the target fraction 00578 if((tD || pD) && tD<tB && pD<pB)// Periferal interaction 00579 { 00580 if(!tD || !pD) // Quasi-Elastic: B+A->(B-1)+N+A || ->B+N+(A-1) 00581 { 00582 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); 00583 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); 00584 G4int pPDG=2112; // Proto of the nucleon PDG (proton) 00585 G4double prM =mNeut; // Proto of the nucleon mass 00586 G4double prM2=mNeu2; // Proto of the nucleon sq mass 00587 if (!tD) // Quasi-elastic scattering of the proj QF nucleon 00588 { 00589 if(pD > 1) pD=1; 00590 if(!pN || (pZ && pA*G4UniformRand() < pZ) ) // proj QF proton 00591 { 00592 pPDG=2212; 00593 prM=mProt; 00594 prM2=mPro2; 00595 pC=1; 00596 } 00597 G4double Mom=0.; 00598 G4LorentzVector com4M = targ4M; // Proto of 4mom of compound 00599 G4double tgM=targ4M.e(); 00600 G4double rNM=0.; 00601 G4LorentzVector rNuc4M(0.,0.,0.,0.); 00602 if(pD==pB) 00603 { 00604 Mom=proj4M.rho(); 00605 com4M += proj4M; // Total 4-mom for scattering 00606 rNM=prM; 00607 } 00608 else // It is necessary to split the nucleon (with fermiM) 00609 { 00610 G4ThreeVector fm=mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); 00611 rNM=G4QPDGCode(2112).GetNuclMass(pZ-pC, pN, 0); 00612 G4double rNE=std::sqrt(fm*fm+rNM*rNM); 00613 rNuc4M=G4LorentzVector(fm,rNE); 00614 G4ThreeVector boostV=proj4M.vect()/proj4M.e(); 00615 rNuc4M.boost(boostV); 00616 G4LorentzVector qfN4M=proj4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS 00617 com4M += qfN4M; // Calculate Total 4Mom for NA scattering 00618 G4double pNE = qfN4M.e(); // Energy of the QF nucleon in LS 00619 if(rNE >= prM) Mom = std::sqrt(pNE*pNE-prM2); // Mom(s) fake value 00620 else break; // Break the while loop 00621 } 00622 xSec=0.; 00623 if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); 00624 else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); 00625 if( xSec <= 0. ) break; // Break the while loop 00626 G4double mint=0.; // Prototype of functional randomized -t 00627 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t 00628 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t 00629 G4double maxt=0.; // Prototype of maximum -t 00630 if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t 00631 else maxt=NCSmanager->GetHMaxT(); // maximum -t 00632 G4double cost=1.-mint/maxt; // cos(theta) in CMS 00633 if(cost>1. || cost<-1.) break; // Break the while loop 00634 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tgM); // 4mom of recoil target 00635 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N 00636 G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01); 00637 if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) 00638 { 00639 G4cout<<"G4QLE::Pt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl; 00640 break; // Break the while loop 00641 } 00642 G4Track* projSpect = 0; 00643 G4Track* aNucTrack = 0; 00644 if(pB > pD) // Fill the proj residual nucleus 00645 { 00646 G4int rZ=pZ-pC; 00647 G4int rA=pB-1; 00648 G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition 00649 if(rA==1) 00650 { 00651 if(!rZ) theDefinition = aNeutron; 00652 else theDefinition = aProton; 00653 } 00654 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ); 00655 G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M); 00656 projSpect = new G4Track(resN, localtime, position ); 00657 projSpect->SetWeight(weight); // weighted 00658 projSpect->SetTouchableHandle(trTouchable); 00659 #ifdef pdebug 00660 G4cout<<"G4QLowEn::PSDI:-->ProjQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl; 00661 #endif 00662 ++nSec; 00663 } 00664 G4ParticleDefinition* theDefinition = G4Neutron::Neutron(); 00665 if(pPDG==2212) theDefinition = G4Proton::Proton(); 00666 G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M); 00667 aNucTrack = new G4Track(scatN, localtime, position ); 00668 aNucTrack->SetWeight(weight); // weighted 00669 aNucTrack->SetTouchableHandle(trTouchable); 00670 #ifdef pdebug 00671 G4cout<<"G4QLowEn::PSDI:-->TgQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl; 00672 #endif 00673 ++nSec; 00674 G4Track* aFraTrack=0; 00675 if(tA==1) 00676 { 00677 if(!tZ) theDefinition = aNeutron; 00678 else theDefinition = aProton; 00679 } 00680 else if(tA==8 && tC==4) // Be8 decay 00681 { 00682 theDefinition = anAlpha; 00683 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha 00684 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha 00685 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) 00686 { 00687 G4cout<<"*G4QLE::TS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; 00688 } 00689 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 00690 aFraTrack = new G4Track(pAl, localtime, position ); 00691 aFraTrack->SetWeight(weight); // weighted 00692 aFraTrack->SetTouchableHandle(trTouchable); 00693 #ifdef pdebug 00694 G4cout<<"G4QLowEn::PSDI:-->TgRQFA4M="<<f4M<<G4endl; 00695 #endif 00696 ++nSec; 00697 reco4M=s4M; 00698 } 00699 else if(tA==5 && (tC==2 || tC==3)) // He5/Li5 decay 00700 { 00701 theDefinition = aProton; 00702 G4double mNuc = mProt; 00703 if(tC==2) 00704 { 00705 theDefinition = aNeutron; 00706 mNuc = mNeut; 00707 } 00708 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon 00709 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha 00710 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) 00711 { 00712 G4cout<<"*G4QLE::TS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; 00713 } 00714 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 00715 aFraTrack = new G4Track(pAl, localtime, position ); 00716 aFraTrack->SetWeight(weight); // weighted 00717 aFraTrack->SetTouchableHandle(trTouchable); 00718 #ifdef pdebug 00719 G4cout<<"G4QLowEn::PSDI:-->TgQFRN4M="<<f4M<<G4endl; 00720 #endif 00721 ++nSec; 00722 theDefinition = anAlpha; 00723 reco4M=s4M; 00724 } 00725 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(tZ,tB,0,tZ); 00726 ++nSec; 00727 #ifdef pdebug 00728 G4cout<<"G4QLowEn::PSDI:-->PQF_nSec="<<nSec<<G4endl; 00729 #endif 00730 aParticleChange.SetNumberOfSecondaries(nSec); 00731 if(projSpect) aParticleChange.AddSecondary(projSpect); 00732 if(aNucTrack) aParticleChange.AddSecondary(aNucTrack); 00733 if(aFraTrack) aParticleChange.AddSecondary(aFraTrack); 00734 G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M); 00735 G4Track* aResTrack = new G4Track(resA, localtime, position ); 00736 aResTrack->SetWeight(weight); // weighted 00737 aResTrack->SetTouchableHandle(trTouchable); 00738 #ifdef pdebug 00739 G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl; 00740 #endif 00741 aParticleChange.AddSecondary(aResTrack); 00742 } 00743 else // !pD : QF target Nucleon on the whole Projectile 00744 { 00745 if(tD > 1) tD=1; 00746 if(!tN || (tZ && tA*G4UniformRand() < tZ) ) // target QF proton 00747 { 00748 pPDG=2212; 00749 prM=mProt; 00750 prM2=mPro2; 00751 tC=1; 00752 } 00753 G4double Mom=0.; 00754 G4LorentzVector com4M=proj4M; // Proto of 4mom of compound 00755 prM=proj4M.m(); 00756 G4double rNM=0.; 00757 G4LorentzVector rNuc4M(0.,0.,0.,0.); 00758 if(tD==tB) 00759 { 00760 Mom=prM*proj4M.rho()/proj4M.m(); 00761 com4M += targ4M; // Total 4-mom for scattering 00762 rNM=prM; 00763 } 00764 else // It is necessary to split the nucleon (with fermiM) 00765 { 00766 G4ThreeVector fm=250.*std::pow(G4UniformRand(),third)*G4RandomDirection(); 00767 rNM=G4QPDGCode(2112).GetNuclMass(tZ-tC, tN, 0)/MeV; 00768 G4double rNE=std::sqrt(fm*fm+rNM*rNM); 00769 rNuc4M=G4LorentzVector(fm,rNE); 00770 G4LorentzVector qfN4M=targ4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS 00771 com4M += qfN4M; // Calculate Total 4Mom for NA scattering 00772 G4ThreeVector boostV=proj4M.vect()/proj4M.e(); 00773 qfN4M.boost(-boostV); 00774 G4double tNE = qfN4M.e(); // Energy of the QF nucleon in LS 00775 if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2); // Mom(s) fake value 00776 else break; // Break the while loop 00777 } 00778 xSec=0.; 00779 if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); 00780 else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); 00781 if( xSec <= 0. ) break; // Break the while loop 00782 G4double mint=0.; // Prototype of functional randomized -t 00783 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t 00784 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t 00785 G4double maxt=0.; // Prototype of maximum -t 00786 if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t 00787 else maxt=NCSmanager->GetHMaxT(); // maximum -t 00788 G4double cost=1.-mint/maxt; // cos(theta) in CMS 00789 if(cost>1. || cost<-1.) break; // Break the while loop 00790 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,prM); // 4mom of recoil target 00791 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N 00792 G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01); 00793 if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) 00794 { 00795 G4cout<<"G4QLE::Tt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl; 00796 break; // Break the while loop 00797 } 00798 G4Track* targSpect = 0; 00799 G4Track* aNucTrack = 0; 00800 if(tB > tD) // Fill the residual nucleus 00801 { 00802 G4int rZ=tZ-tC; 00803 G4int rA=tB-1; 00804 G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition 00805 if(rA==1) 00806 { 00807 if(!rZ) theDefinition = aNeutron; 00808 else theDefinition = aProton; 00809 } 00810 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ); 00811 G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M); 00812 targSpect = new G4Track(resN, localtime, position ); 00813 targSpect->SetWeight(weight); // weighted 00814 targSpect->SetTouchableHandle(trTouchable); 00815 #ifdef pdebug 00816 G4cout<<"G4QLowEn::PSDI:-->TargQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl; 00817 #endif 00818 ++nSec; 00819 } 00820 G4ParticleDefinition* theDefinition = G4Neutron::Neutron(); 00821 if(pPDG==2212) theDefinition = G4Proton::Proton(); 00822 G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M); 00823 aNucTrack = new G4Track(scatN, localtime, position ); 00824 aNucTrack->SetWeight(weight); // weighted 00825 aNucTrack->SetTouchableHandle(trTouchable); 00826 #ifdef pdebug 00827 G4cout<<"G4QLowEn::PSDI:-->PrQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl; 00828 #endif 00829 ++nSec; 00830 G4Track* aFraTrack=0; 00831 if(pA==1) 00832 { 00833 if(!pZ) theDefinition = aNeutron; 00834 else theDefinition = aProton; 00835 } 00836 else if(pA==8 && pC==4) // Be8 decay 00837 { 00838 theDefinition = anAlpha; 00839 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha 00840 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha 00841 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) 00842 { 00843 G4cout<<"*G4QLE::PS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; 00844 } 00845 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 00846 aFraTrack = new G4Track(pAl, localtime, position ); 00847 aFraTrack->SetWeight(weight); // weighted 00848 aFraTrack->SetTouchableHandle(trTouchable); 00849 #ifdef pdebug 00850 G4cout<<"G4QLowEn::PSDI:-->PrRQFA4M="<<f4M<<G4endl; 00851 #endif 00852 ++nSec; 00853 reco4M=s4M; 00854 } 00855 else if(pA==5 && (pC==2 || pC==3)) // He5/Li5 decay 00856 { 00857 theDefinition = aProton; 00858 G4double mNuc = mProt; 00859 if(pC==2) 00860 { 00861 theDefinition = aNeutron; 00862 mNuc = mNeut; 00863 } 00864 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon 00865 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha 00866 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) 00867 { 00868 G4cout<<"*G4QLE::PS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; 00869 } 00870 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 00871 aFraTrack = new G4Track(pAl, localtime, position ); 00872 aFraTrack->SetWeight(weight); // weighted 00873 aFraTrack->SetTouchableHandle(trTouchable); 00874 #ifdef pdebug 00875 G4cout<<"G4QLowEn::PSDI:-->PrQFRN4M="<<f4M<<G4endl; 00876 #endif 00877 ++nSec; 00878 theDefinition = anAlpha; 00879 reco4M=s4M; 00880 } 00881 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(pZ,pB,0,pZ); 00882 ++nSec; 00883 #ifdef pdebug 00884 G4cout<<"G4QLowEn::PSDI:-->TQF_nSec="<<nSec<<G4endl; 00885 #endif 00886 aParticleChange.SetNumberOfSecondaries(nSec); 00887 if(targSpect) aParticleChange.AddSecondary(targSpect); 00888 if(aNucTrack) aParticleChange.AddSecondary(aNucTrack); 00889 if(aFraTrack) aParticleChange.AddSecondary(aFraTrack); 00890 G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M); 00891 G4Track* aResTrack = new G4Track(resA, localtime, position ); 00892 aResTrack->SetWeight(weight); // weighted 00893 aResTrack->SetTouchableHandle(trTouchable); 00894 #ifdef pdebug 00895 G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl; 00896 #endif 00897 aParticleChange.AddSecondary( aResTrack ); 00898 } 00899 #ifdef debug 00900 G4cout<<"G4QLowEnergy::PostStepDoIt:***PostStepDoIt is done:Quasi-El***"<<G4endl; 00901 #endif 00902 return G4VDiscreteProcess::PostStepDoIt(track, step); // ===> Reaction is DONE 00903 } 00904 else // The cental region compound can be created 00905 { 00906 // First calculate the isotopic state of the parts of the compound 00907 if (!pZ) pC = 0; 00908 else if(!pN) pC = pD; 00909 else 00910 { 00911 #ifdef pdebug 00912 G4cout<<"G4QLowEn::PSDI: pD="<<pD<<", pZ="<<pZ<<", pA="<<pA<<G4endl; 00913 #endif 00914 G4double C=pD*pZ/pA; 00915 pC=static_cast<G4int>(C); 00916 if(G4UniformRand() < C-pC) ++pC; 00917 } 00918 if(!tZ) tC=0; 00919 else if(!tN) tC=tD; 00920 else 00921 { 00922 #ifdef pdebug 00923 G4cout<<"G4QLowEn::PSDI: tD="<<tD<<", tZ="<<tZ<<", tA="<<tA<<G4endl; 00924 #endif 00925 G4double C=tD*tZ/tA; 00926 tC=static_cast<G4int>(C); 00927 if(G4UniformRand() < C-tC) ++tC; 00928 } 00929 // calculate the transferred momentum 00930 G4ThreeVector pFM(0.,0.,0.); 00931 if(pD < pB) // The projectile nucleus must be splitted 00932 { 00933 G4int nc=pD; 00934 if(pD+pD>pB) nc=pB-pD; 00935 pFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); 00936 for(i=1; i < nc; ++i) 00937 pFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); 00938 } 00939 G4ThreeVector tFM(0.,0.,0.); 00940 if(tD<tB) // The projectile nucleus must be splitted 00941 { 00942 G4int nc=pD; 00943 if(tD+tD>tB) nc=tB-tD; 00944 tFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); 00945 for(i=1; i < nc; ++i) 00946 tFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); 00947 } 00948 #ifdef pdebug 00949 G4cout<<"G4QLE::PSDI:pC="<<pC<<", tC="<<tC<<", pFM="<<pFM<<", tFM="<<tFM<<G4endl; 00950 #endif 00951 // Split the projectile spectator 00952 G4int rpZ=pZ-pC; // number of protons in the projectile spectator 00953 G4int pF_value=pD-pC; // number of neutrons in the projectile part of comp 00954 G4int rpN=pN-pF_value; // number of neutrons in the projectile spectator 00955 G4double rpNM=G4QPDGCode(2112).GetNuclMass(rpZ, rpN, 0); // Mass of the spectator 00956 G4ThreeVector boostV=proj4M.vect()/proj4M.e(); // Antilab Boost Vector 00957 G4double rpE=std::sqrt(rpNM*rpNM+pFM.mag2()); 00958 G4LorentzVector rp4M(pFM,rpE); 00959 #ifdef pdebug 00960 G4cout<<"G4QLE::PSDI: boostV="<<boostV<<",rp4M="<<rp4M<<",pr4M="<<proj4M<<G4endl; 00961 #endif 00962 rp4M.boost(boostV); 00963 #ifdef pdebug 00964 G4cout<<"G4QLE::PSDI: After boost, rp4M="<<rp4M<<G4endl; 00965 #endif 00966 G4ParticleDefinition* theDefinition; // Prototype of projSpectatorNuclDefinition 00967 G4int rpA=rpZ+rpN; 00968 G4Track* aFraPTrack = 0; 00969 theDefinition = 0; 00970 if(rpA==1) 00971 { 00972 if(!rpZ) theDefinition = G4Neutron::Neutron(); 00973 else theDefinition = G4Proton::Proton(); 00974 #ifdef pdebug 00975 G4cout<<"G4QLE::PSDI: rpA=1, rpZ"<<rpZ<<G4endl; 00976 #endif 00977 } 00978 else if(rpA==2 && rpZ==0) // nn decay 00979 { 00980 theDefinition = aNeutron; 00981 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron 00982 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron 00983 #ifdef pdebug 00984 G4cout<<"G4QLE::CPS->n+n,nn="<<rp4M.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl; 00985 #endif 00986 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M)) 00987 { 00988 G4cout<<"*W*G4QLE::CPS->n+n,t="<<rp4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl; 00989 } 00990 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M); 00991 aFraPTrack = new G4Track(pNeu, localtime, position ); 00992 aFraPTrack->SetWeight(weight); // weighted 00993 aFraPTrack->SetTouchableHandle(trTouchable); 00994 tt4M-=f4M; 00995 #ifdef edebug 00996 totBaN-=2; 00997 tch4M -=f4M; 00998 G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 00999 #endif 01000 #ifdef pdebug 01001 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; 01002 #endif 01003 rp4M=s4M; 01004 } 01005 else if(rpA>2 && rpZ==0) // Z=0 decay 01006 { 01007 theDefinition = aNeutron; 01008 G4LorentzVector f4M=rp4M/rpA; // 4mom of 1st neutron 01009 #ifdef pdebug 01010 G4cout<<"G4QLE::CPS->Nn,M="<<rp4M.m()<<" >? N*MNeutron="<<rpA*mNeutron<<G4endl; 01011 #endif 01012 for(G4int it=1; it<rpA; ++it) // Fill (N-1) neutrons to output 01013 { 01014 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M); 01015 G4Track* aNTrack = new G4Track(pNeu, localtime, position ); 01016 aNTrack->SetWeight(weight); // weighted 01017 aNTrack->SetTouchableHandle(trTouchable); 01018 result.push_back(aNTrack); 01019 } 01020 G4int nesc = rpA-1; 01021 tt4M-=f4M*nesc; 01022 #ifdef edebug 01023 totBaN-=nesc; 01024 tch4M -=f4M*nesc; 01025 G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01026 #endif 01027 #ifdef pdebug 01028 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; 01029 #endif 01030 rp4M=f4M; 01031 } 01032 else if(rpA==8 && rpZ==4) // Be8 decay 01033 { 01034 theDefinition = anAlpha; 01035 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha 01036 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha 01037 #ifdef pdebug 01038 G4cout<<"G4QLE::CPS->A+A,mBe8="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; 01039 #endif 01040 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M)) 01041 { 01042 G4cout<<"*W*G4QLE::CPS->A+A,t="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; 01043 } 01044 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 01045 aFraPTrack = new G4Track(pAl, localtime, position ); 01046 aFraPTrack->SetWeight(weight); // weighted 01047 aFraPTrack->SetTouchableHandle(trTouchable); 01048 tt4M-=f4M; 01049 #ifdef edebug 01050 totChg-=2; 01051 totBaN-=4; 01052 tch4M -=f4M; 01053 G4cout<<">>G4QLEn::PSDI:1,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01054 #endif 01055 #ifdef pdebug 01056 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; 01057 #endif 01058 rp4M=s4M; 01059 } 01060 else if(rpA==5 && (rpZ==2 || rpZ==3)) // He5/Li5 decay 01061 { 01062 theDefinition = aProton; 01063 G4double mNuc = mProt; 01064 if(rpZ==2) 01065 { 01066 theDefinition = aNeutron; 01067 mNuc = mNeut; 01068 } 01069 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon 01070 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha 01071 #ifdef pdebug 01072 G4cout<<"G4QLowE::CPS->N+A, tM5="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; 01073 #endif 01074 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M)) 01075 { 01076 G4cout<<"*W*G4QLE::CPS->N+A,t="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; 01077 } 01078 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 01079 aFraPTrack = new G4Track(pAl, localtime, position ); 01080 aFraPTrack->SetWeight(weight); // weighted 01081 aFraPTrack->SetTouchableHandle(trTouchable); 01082 tt4M-=f4M; 01083 #ifdef edebug 01084 if(theDefinition == aProton) totChg-=1; 01085 totBaN-=1; 01086 tch4M -=f4M; 01087 G4cout<<">>G4QLEn::PSDI:2,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01088 #endif 01089 #ifdef pdebug 01090 G4cout<<"G4QLowEn::PSDI:-->ProjSpectN4M="<<f4M<<G4endl; 01091 #endif 01092 theDefinition = anAlpha; 01093 rp4M=s4M; 01094 } 01095 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rpZ,rpA,0,rpZ); 01096 if(!theDefinition) 01097 { 01098 // G4cout<<"-Warning-G4QLowEn::PSDI: pDef=0, Z="<<rpZ<<", A="<<rpA<<G4endl; 01099 // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer"); 01100 G4ExceptionDescription ed; 01101 ed << "Particle definition is a null pointer: pDef=0, Z= " << rpZ 01102 << ", A=" << rpA << G4endl; 01103 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000", 01104 JustWarning, ed); 01105 } 01106 #ifdef edebug 01107 if(theDefinition == anAlpha) 01108 { 01109 totChg-=2; 01110 totBaN-=4; 01111 } 01112 else 01113 { 01114 totChg-=rpZ; 01115 totBaN-=rpA; 01116 } 01117 tch4M -=rp4M; 01118 G4cout<<">>G4QLEn::PSDI:3, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl; 01119 #endif 01120 G4DynamicParticle* rpD = new G4DynamicParticle(theDefinition, rp4M); 01121 G4Track* aNewPTrack = new G4Track(rpD, localtime, position); 01122 aNewPTrack->SetWeight(weight);// weighted 01123 aNewPTrack->SetTouchableHandle(trTouchable); 01124 tt4M-=rp4M; 01125 #ifdef pdebug 01126 G4cout<<"G4QLowEn::PSDI:-->ProjSpectR4M="<<rp4M<<",Z="<<rpZ<<",A="<<rpA<<G4endl; 01127 #endif 01128 // 01129 // Split the target spectator 01130 G4int rtZ=tZ-tC; // number of protons in the target spectator 01131 G4int tF_value=tD-tC; // number of neutrons in the target part of comp 01132 G4int rtN=tN-tF_value; // number of neutrons in the target spectator 01133 G4double rtNM=G4QPDGCode(2112).GetNuclMass(rtZ, rtN, 0); // Mass of the spectator 01134 G4double rtE=std::sqrt(rtNM*rtNM+tFM.mag2()); 01135 G4LorentzVector rt4M(tFM,rtE); 01136 G4int rtA=rtZ+rtN; 01137 G4Track* aFraTTrack = 0; 01138 theDefinition = 0; 01139 if(rtA==1) 01140 { 01141 if(!rtZ) theDefinition = G4Neutron::Neutron(); 01142 else theDefinition = G4Proton::Proton(); 01143 #ifdef pdebug 01144 G4cout<<"G4QLE::PSDI: rtA=1, rtZ"<<rtZ<<G4endl; 01145 #endif 01146 } 01147 else if(rtA==2 && rtZ==0) // nn decay 01148 { 01149 theDefinition = aNeutron; 01150 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron 01151 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron 01152 #ifdef pdebug 01153 G4cout<<"G4QLE::CPS->n+n,nn="<<rptM.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl; 01154 #endif 01155 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M)) 01156 G4cout<<"*W*G4QLE::CPS->n+n,t="<<rt4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl; 01157 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M); 01158 aFraPTrack = new G4Track(pNeu, localtime, position ); 01159 aFraPTrack->SetWeight(weight); // weighted 01160 aFraPTrack->SetTouchableHandle(trTouchable); 01161 tt4M-=f4M; 01162 #ifdef edebug 01163 totBaN-=2; 01164 tch4M -=f4M; 01165 G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01166 #endif 01167 #ifdef pdebug 01168 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; 01169 #endif 01170 rt4M=s4M; 01171 } 01172 else if(rtA>2 && rtZ==0) // Z=0 decay 01173 { 01174 theDefinition = aNeutron; 01175 G4LorentzVector f4M=rt4M/rtA; // 4mom of 1st neutron 01176 #ifdef pdebug 01177 G4cout<<"G4QLE::CPS->Nn,M="<<rt4M.m()<<" >? N*MNeutron="<<rtA*mNeutron<<G4endl; 01178 #endif 01179 for(G4int it=1; it<rtA; ++it) // Fill (N-1) neutrons to output 01180 { 01181 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M); 01182 G4Track* aNTrack = new G4Track(pNeu, localtime, position ); 01183 aNTrack->SetWeight(weight); // weighted 01184 aNTrack->SetTouchableHandle(trTouchable); 01185 result.push_back(aNTrack); 01186 } 01187 G4int nesc = rtA-1; 01188 tt4M-=f4M*nesc; 01189 #ifdef edebug 01190 totBaN-=nesc; 01191 tch4M -=f4M*nesc; 01192 G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01193 #endif 01194 #ifdef pdebug 01195 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; 01196 #endif 01197 rt4M=f4M; 01198 } 01199 else if(rtA==8 && rtZ==4) // Be8 decay 01200 { 01201 theDefinition = anAlpha; 01202 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha 01203 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha 01204 #ifdef pdebug 01205 G4cout<<"G4QLE::CPS->A+A,mBe8="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; 01206 #endif 01207 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M)) 01208 { 01209 G4cout<<"*W*G4QLE::CTS->A+A,t="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; 01210 } 01211 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 01212 aFraTTrack = new G4Track(pAl, localtime, position ); 01213 aFraTTrack->SetWeight(weight); // weighted 01214 aFraTTrack->SetTouchableHandle(trTouchable); 01215 tt4M-=f4M; 01216 #ifdef edebug 01217 totChg-=2; 01218 totBaN-=4; 01219 tch4M -=f4M; 01220 G4cout<<">>G4QLEn::PSDI:4,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01221 #endif 01222 #ifdef pdebug 01223 G4cout<<"G4QLowEn::PSDI:-->TargSpectA4M="<<f4M<<G4endl; 01224 #endif 01225 rt4M=s4M; 01226 } 01227 else if(rtA==5 && (rtZ==2 || rtZ==3)) // He5/Li5 decay 01228 { 01229 theDefinition = aProton; 01230 G4double mNuc = mProt; 01231 if(rpZ==2) 01232 { 01233 theDefinition = aNeutron; 01234 mNuc = mNeut; 01235 } 01236 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon 01237 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha 01238 #ifdef pdebug 01239 G4cout<<"G4QLowE::CPS->N+A, tM5="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; 01240 #endif 01241 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M)) 01242 { 01243 G4cout<<"*W*G4QLE::CTS->N+A,t="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; 01244 } 01245 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); 01246 aFraTTrack = new G4Track(pAl, localtime, position ); 01247 aFraTTrack->SetWeight(weight); // weighted 01248 aFraTTrack->SetTouchableHandle(trTouchable); 01249 tt4M-=f4M; 01250 #ifdef edebug 01251 if(theDefinition == aProton) totChg-=1; 01252 totBaN-=1; 01253 tch4M -=f4M; 01254 G4cout<<">>G4QLEn::PSDI:5,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; 01255 #endif 01256 #ifdef pdebug 01257 G4cout<<"G4QLowEn::PSDI:-->TargSpectN4M="<<f4M<<G4endl; 01258 #endif 01259 theDefinition = anAlpha; 01260 rt4M=s4M; 01261 } 01262 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rtZ,rtA,0,rtZ); 01263 if(!theDefinition) 01264 { 01265 // G4cout<<"-Warning-G4QLowEn::PSDI: tDef=0, Z="<<rtZ<<", A="<<rtA<<G4endl; 01266 // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer"); 01267 G4ExceptionDescription ed; 01268 ed << "Particle definition is a null pointer: tDef=0, Z= " << rtZ 01269 << ", A=" << rtA << G4endl; 01270 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000", 01271 JustWarning, ed); 01272 } 01273 #ifdef edebug 01274 if(theDefinition == anAlpha) 01275 { 01276 totChg-=2; 01277 totBaN-=4; 01278 } 01279 else 01280 { 01281 totChg-=rtZ; 01282 totBaN-=rtA; 01283 } 01284 tch4M -=rt4M; 01285 G4cout<<">>G4QLEn::PSDI:6, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl; 01286 #endif 01287 G4DynamicParticle* rtD = new G4DynamicParticle(theDefinition, rt4M); 01288 G4Track* aNewTTrack = new G4Track(rtD, localtime, position); 01289 aNewTTrack->SetWeight(weight); // weighted 01290 aNewTTrack->SetTouchableHandle(trTouchable); 01291 tt4M-=rt4M; 01292 #ifdef pdebug 01293 G4cout<<"G4QLowEn::PSDI:-->TargSpectR4M="<<rt4M<<",Z="<<rtZ<<",A="<<rtA<<G4endl; 01294 #endif 01295 if(aFraPTrack) result.push_back(aFraPTrack); 01296 if(aNewPTrack) result.push_back(aNewPTrack); 01297 if(aFraTTrack) result.push_back(aFraTTrack); 01298 if(aNewTTrack) result.push_back(aNewTTrack); 01299 tcM = tt4M.m(); // total CMS mass of the compound (after reduction!) 01300 G4int sN=tF_value+pF_value; 01301 G4int sZ=tC+pC; 01302 tnM = targQPDG.GetNuclMass(sZ,sN,0); // GSM 01303 #ifdef pdebug 01304 G4cout<<"G4QLEn::PSDI:At#"<<nAt<<",totM="<<tcM<<",gsM="<<tnM<<",Z="<<sZ<<",N=" 01305 <<sN<<G4endl; 01306 #endif 01307 if(tcM > tnM) // !! The totalresidual reaction is modified ! 01308 { 01309 pZ=pC; 01310 pN=pF_value; 01311 tZ=tC; 01312 tN=tF_value; 01313 tot4M=tt4M; 01314 } 01315 //else 01316 //{ 01317 // G4cout<<"***G4QLowEnergy::PostStepDoIt: totM="<<tcM<<" <= GSM="<<tnM<<G4endl; 01318 // throw G4QException("G4QLowEnergy::PostStepDoIt: ResibualNucl below GSM shell"); 01319 //} 01320 } 01321 } // At least one is splitted 01322 else if(tD==tB || pD==pB) // Total absorption 01323 { 01324 tD=tZ+tN; 01325 pD=pZ+pN; 01326 tcM=tnM+1.; 01327 } 01328 } 01329 else // Total compound (define tD to go out of the while) 01330 { 01331 tD=tZ+tN; 01332 pD=pZ+pN; 01333 tcM=tnM+1.; 01334 } 01335 } // End of the interaction WHILE 01336 G4double totM=tot4M.m(); // total CMS mass of the reaction 01337 G4int totN=tN+pN; 01338 G4int totZ=tZ+pZ; 01339 #ifdef edebug 01340 G4cout<<">>>G4QLEn::PSDI: dZ="<<totChg-totZ<<", dB="<<totBaN-totN-totZ<<", d4M=" 01341 <<tch4M-tot4M<<",N="<<totN<<",Z="<<totZ<<G4endl; 01342 #endif 01343 // @@ Here mass[i] can be calculated if mass=0 01344 G4double mass[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 01345 0.,0.,0.,0.,0.,0.}; 01346 mass[0] = tM = targQPDG.GetNuclMass(totZ, totN, 0); // gamma+gamma and GSM update 01347 #ifdef pdebug 01348 G4cout<<"G4QLEn::PSDI:TotM="<<totM<<",NucM="<<tM<<",totZ="<<totZ<<",totN="<<totN<<G4endl; 01349 #endif 01350 if (totN>0 && totZ>0) 01351 { 01352 mass[1] = targQPDG.GetNuclMass(totZ,totN-1,0); // gamma+neutron 01353 mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0); // gamma+proton 01354 } 01355 if ( totZ > 1 && totN > 1 ) mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0); //g+d 01356 if ( totZ > 1 && totN > 2 ) mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0); //g+t 01357 if ( totZ > 2 && totN > 1 ) mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0); //g+3 01358 if ( totZ > 2 && totN > 2 ) mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0); //g+a 01359 if ( totZ > 0 && totN > 2 ) mass[7] = targQPDG.GetNuclMass(totZ ,totN-2,0); //n+n 01360 mass[ 8] = mass[3]; // neutron+proton (the same as a deuteron) 01361 if ( totZ > 2 ) mass[9] = targQPDG.GetNuclMass(totZ-2,totN ,0); //p+p 01362 mass[10] = mass[5]; // proton+deuteron (the same as He3) 01363 mass[11] = mass[4]; // neutron+deuteron (the same as t) 01364 mass[12] = mass[6]; // deuteron+deuteron (the same as alpha) 01365 mass[13] = mass[6]; // proton+tritium (the same as alpha) 01366 if ( totZ > 1 && totN > 3 ) mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0);//n+t 01367 if ( totZ > 3 && totN > 1 ) mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0);//He3+p 01368 mass[16] = mass[6]; // neutron+He3 (the same as alpha) 01369 if ( totZ > 3 && totN > 2 ) mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0);//pa 01370 if ( totZ > 2 && totN > 3 ) mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0);//na 01371 if(pL>0) // @@ Not debugged @@ 01372 { 01373 G4int pL1=pL-1; 01374 if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ ,totN ,pL1);// Lambda+gamma 01375 if( (totN > 0 && totZ > 0) || totZ > 1 ) 01376 mass[20]=targQPDG.GetNuclMass(totZ-1,totN ,pL1);//Lp 01377 if( (totN > 0 && totZ > 0) || totN > 0 ) 01378 mass[21]=targQPDG.GetNuclMass(totZ ,totN-1,pL1);//Ln 01379 if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) ) 01380 mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);//Ld 01381 if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) ) 01382 mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);//Lt 01383 if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) ) 01384 mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);//L3 01385 if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) ) 01386 mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);//La 01387 } 01388 #ifdef debug 01389 G4cout<<"G4QLowEn::PSDI: Residual masses are calculated"<<G4endl; 01390 #endif 01391 tA=tZ+tN; 01392 pA=pZ+pN; 01393 G4double prZ=pZ/pA+tZ/tA; 01394 G4double prN=pN/pA+tN/tA; 01395 G4double prD=prN*prZ; 01396 G4double prA=prD*prD; 01397 G4double prH=prD*prZ; 01398 G4double prT=prD*prN; 01399 G4double fhe3=6.*std::sqrt(tA); 01400 G4double prL=0.; 01401 if(pL>0) prL=pL/pA; 01402 G4double qval[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 01403 0.,0.,0.,0.,0.,0.}; 01404 qval[ 0] = (totM - mass[ 0])/137./137.; 01405 qval[ 1] = (totM - mass[ 1] - mNeut)*prN/137.; 01406 qval[ 2] = (totM - mass[ 2] - mProt)*prZ/137.; 01407 qval[ 3] = (totM - mass[ 3] - mDeut)*prD/3./137.; 01408 qval[ 4] = (totM - mass[ 4] - mTrit)*prT/6./137.; 01409 qval[ 5] = (totM - mass[ 5] - mHel3)*prH/fhe3/137.; 01410 qval[ 6] = (totM - mass[ 6] - mAlph)*prA/9./137.; 01411 qval[ 7] = (totM - mass[ 7] - mNeut - mNeut)*prN*prN; 01412 qval[ 8] = (totM - mass[ 8] - mNeut - mProt)*prD; 01413 qval[ 9] = (totM - mass[ 9] - mProt - mProt)*prZ*prZ; 01414 qval[10] = (totM - mass[10] - mProt - mDeut)*prH/3.; 01415 qval[11] = (totM - mass[11] - mNeut - mDeut)*prT/3.; 01416 qval[12] = (totM - mass[12] - mDeut - mDeut)*prA/3./3.; 01417 qval[13] = (totM - mass[13] - mProt - mTrit)*prA/6.; 01418 qval[14] = (totM - mass[14] - mNeut - mTrit)*prT*prN/6.; 01419 qval[15] = (totM - mass[15] - mProt - mHel3)*prH*prZ/fhe3; 01420 qval[16] = (totM - mass[16] - mNeut - mHel3)*prA/fhe3; 01421 qval[17] = (totM - mass[17] - mProt - mAlph)*prZ*prA/9.; 01422 qval[18] = (totM - mass[18] - mNeut - mAlph)*prN*prA/9.; 01423 if(pZ>0) 01424 { 01425 qval[19] = (totM - mass[19] - mLamb)*prL; 01426 qval[20] = (totM - mass[20] - mProt - mLamb)*prL*prZ; 01427 qval[21] = (totM - mass[21] - mNeut - mLamb)*prL*prN; 01428 qval[22] = (totM - mass[22] - mDeut - mLamb)*prL*prD/2.; 01429 qval[23] = (totM - mass[23] - mTrit - mLamb)*prL*prT/3.; 01430 qval[24] = (totM - mass[24] - mHel3 - mLamb)*prL*prH/fhe3; 01431 qval[25] = (totM - mass[25] - mAlph - mLamb)*prL*prA/4; 01432 } 01433 #ifdef debug 01434 G4cout<<"G4QLowEn::PSDI: Q-values are calculated, tgA="<<tA<<", prA="<<pA<<G4endl; 01435 #endif 01436 01437 G4double qv = 0.0; // Total sum of probabilities (q-values) 01438 for(i=0; i<nCh; ++i ) 01439 { 01440 #ifdef sdebug 01441 G4cout<<"G4QLowEn::PSDI: i="<<i<<", q="<<qval[i]<<",m="<<mass[i]<<G4endl; 01442 #endif 01443 if( mass[i] < 500.*MeV ) qval[i] = 0.0; // Close A/Z impossible channels (mesons) 01444 if( qval[i] < 0.0 ) qval[i] = 0.0; // Close the splitting impossible channels 01445 qv += qval[i]; 01446 } 01447 // Select the channel 01448 G4double qv1 = 0.0; 01449 G4double ran = qv*G4UniformRand(); 01450 G4int index = 0; 01451 for( index=0; index<nCh; ++index ) if( qval[index] > 0.0 ) 01452 { 01453 qv1 += qval[index]; 01454 if( ran <= qv1 ) break; 01455 } 01456 #ifdef debug 01457 G4cout<<"G4QLowEn::PSDI: index="<<index<<" < "<<nCh<<G4endl; 01458 #endif 01459 if(index == nCh) 01460 { 01461 G4cout<<"***G4QLowEnergy::PoStDI:Decay is impossible,totM="<<totM<<",GSM="<<tM<<G4endl; 01462 G4cout<<"G4QLowEn::PoStDI:Reaction "<<projPDG<<"+"<<targPDG<<", P="<<Momentum<<G4endl; 01463 G4int nRes=result.size(); 01464 if(nRes)G4cout<<"G4QLowEnergy::PoStDI: result exists with N="<<nRes<<" tracks"<<G4endl; 01465 // else throw G4QException("***G4QLowEnergy::PostStepDoIt: Can't decay ResidualCompound"); 01466 else { 01467 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000", 01468 FatalException, "Can't decay ResidualCompound"); 01469 } 01470 G4double minEx=1000000.; // Prototype of minimal excess 01471 G4bool found = false; // The solution is not found yet 01472 G4int cInd = 0; // Prototype of the best index 01473 G4int ctN = 0; // To 01474 G4int ctZ = 0; // avoid 01475 G4LorentzVector c4M=tot4M; // recalculation 01476 G4double ctM=0.; // of found 01477 for(G4int it=0; it<nRes; ++it) 01478 { 01479 G4Track* iTrack = result[it]; 01480 const G4DynamicParticle* iHadron = iTrack->GetDynamicParticle(); 01481 G4ParticleDefinition* iParticle=iHadron->GetDefinition(); 01482 G4int iPDG = iParticle->GetPDGEncoding(); 01483 G4LorentzVector ih4M = projHadron->Get4Momentum(); 01484 G4cout<<" G4QLowEn::PoStDI: H["<<it<<"] = [ "<<iPDG<<", "<<ih4M<<" ]"<<G4endl; 01485 G4int iB = iParticle->GetBaryonNumber(); // A of the secondary 01486 G4int iC = G4int(iParticle->GetPDGCharge()); // Z of the secondary 01487 G4LorentzVector new4M = tot4M + ih4M; // Make temporary tot4M 01488 G4int ntN=totN + iB-iC; // Make temporary totN 01489 G4int ntZ=totZ + iC; // Make temporary totZ 01490 G4double ntM = targQPDG.GetNuclMass(ntZ, ntN, 0); // Make temporary GSM 01491 G4double ttM = new4M.m(); 01492 if(ttM > ntM) // >>> This is at least one possible solution ! <<< 01493 { 01494 G4double nEx = ttM - ntM; 01495 if(found && nEx < minEx) // This one is better than the previous found 01496 { 01497 cInd = it; 01498 minEx= nEx; 01499 ctN = ntN; 01500 ctZ = ntZ; 01501 c4M = new4M; 01502 ctM = ttM; 01503 mass[0] = ntM; 01504 } 01505 found = true; 01506 } 01507 } 01508 if(found) 01509 { 01510 tot4M = c4M; 01511 totM = ctM; 01512 totN = ctN; 01513 totZ = ctZ; 01514 delete result[cInd]; 01515 G4int nR1 = nRes -1; 01516 if(cInd < nR1) result[cInd] = result[nR1]; 01517 result.pop_back(); 01518 //nRes=nR1; // @@ can be used for the two pt correction 01519 index = 0; // @@ emergency gamma+gamma decay is selected 01520 } 01521 else 01522 { 01523 G4cout<<"***G4QLowEnergy::PoStDI: One hadron coddection did not succeed***"<<G4endl; 01524 if(nRes>1) G4cout<<"***G4QLowEnergy::PoStDI:nRes's big enough to use 2PtCor"<<G4endl; 01525 // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't correct by one particle."); 01526 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000", 01527 FatalException, "Can't correct by one particle"); 01528 } 01529 } 01530 // @@ Convert it to G4QHadrons 01531 G4DynamicParticle* ResSec = new G4DynamicParticle; 01532 G4DynamicParticle* FstSec = new G4DynamicParticle; 01533 G4DynamicParticle* SecSec = new G4DynamicParticle; 01534 #ifdef debug 01535 G4cout<<"G4QLowEn::PSDI: Dynamic particles are created pL="<<pL<<G4endl; 01536 #endif 01537 01538 G4LorentzVector res4Mom(zeroMom,mass[index]*MeV); // The recoil nucleus prototype 01539 G4double mF=0.; 01540 G4double mS=0.; 01541 G4int rA=totZ+totN; 01542 G4int rZ=totZ; 01543 G4int rL=pL; 01544 G4int complete=3; 01545 G4ParticleDefinition* theDefinition; // Prototype for qfNucleon 01546 switch( index ) 01547 { 01548 case 0: 01549 if(!evaporate || rA<2) 01550 { 01551 if(!rZ) theDefinition=aNeutron; 01552 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01553 if(!theDefinition) 01554 { 01555 // G4cerr<<"-Warning-G4LE::PSDI: notDef(0), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01556 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01557 G4ExceptionDescription ed; 01558 ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ 01559 << ", A=" << rA << ", L=" << rL << G4endl; 01560 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000", 01561 JustWarning, ed); 01562 } 01563 ResSec->SetDefinition( theDefinition ); 01564 FstSec->SetDefinition( aGamma ); 01565 SecSec->SetDefinition( aGamma ); 01566 } 01567 else 01568 { 01569 delete ResSec; 01570 delete FstSec; 01571 delete SecSec; 01572 complete=0; 01573 } 01574 break; 01575 case 1: 01576 rA-=1; // gamma+n 01577 if(!evaporate || rA<2) 01578 { 01579 if(!rZ) theDefinition=aNeutron; 01580 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01581 if(!theDefinition) 01582 { 01583 // G4cerr<<"-Warning-G4LE::PSDI: notDef(1), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01584 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01585 G4ExceptionDescription ed; 01586 ed << "Particle definition is a null pointer: notDef(1), Z= " << rZ 01587 << ", A=" << rA << ", L=" << rL << G4endl; 01588 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0001", 01589 JustWarning, ed); 01590 } 01591 ResSec->SetDefinition( theDefinition ); 01592 SecSec->SetDefinition( aGamma ); 01593 } 01594 else 01595 { 01596 delete ResSec; 01597 delete SecSec; 01598 complete=1; 01599 } 01600 FstSec->SetDefinition( aNeutron ); 01601 mF=mNeut; // First hadron 4-momentum 01602 break; 01603 case 2: 01604 rA-=1; 01605 rZ-=1; // gamma+p 01606 if(!evaporate || rA<2) 01607 { 01608 if(!rZ) theDefinition=aNeutron; 01609 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01610 if(!theDefinition) 01611 { 01612 // G4cerr<<"-Warning-G4LE::PSDI: notDef(2), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01613 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01614 G4ExceptionDescription ed; 01615 ed << "Particle definition is a null pointer: notDef(2), Z= " << rZ 01616 << ", A=" << rA << ", L="<< rL << G4endl; 01617 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0002", 01618 JustWarning, ed); 01619 } 01620 ResSec->SetDefinition( theDefinition ); 01621 SecSec->SetDefinition( aGamma ); 01622 } 01623 else 01624 { 01625 delete ResSec; 01626 delete SecSec; 01627 complete=1; 01628 } 01629 FstSec->SetDefinition( aProton ); 01630 mF=mProt; // First hadron 4-momentum 01631 break; 01632 case 3: 01633 rA-=2; 01634 rZ-=1; // gamma+d 01635 if(!evaporate || rA<2) 01636 { 01637 if(!rZ) theDefinition=aNeutron; 01638 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01639 if(!theDefinition) 01640 { 01641 // G4cerr<<"-Warning-G4LE::PSDI: notDef(3), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01642 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01643 G4ExceptionDescription ed; 01644 ed << "Particle definition is a null pointer: notDef(3), Z= " << rZ 01645 << ", A=" << rA << ", L=" << rL << G4endl; 01646 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0003", 01647 JustWarning, ed); 01648 } 01649 ResSec->SetDefinition( theDefinition ); 01650 SecSec->SetDefinition( aGamma ); 01651 } 01652 else 01653 { 01654 delete ResSec; 01655 delete SecSec; 01656 complete=1; 01657 } 01658 FstSec->SetDefinition( aDeuteron ); 01659 mF=mDeut; // First hadron 4-momentum 01660 break; 01661 case 4: 01662 rA-=3; // gamma+t 01663 rZ-=1; 01664 if(!evaporate || rA<2) 01665 { 01666 if(!rZ) theDefinition=aNeutron; 01667 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01668 if(!theDefinition) 01669 { 01670 // G4cerr<<"-Warning-G4LE::PSDI: notDef(4), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01671 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01672 G4ExceptionDescription ed; 01673 ed << "Particle definition is a null pointer: notDef(4), Z= " << rZ 01674 << ", A=" << rA << ", L=" << rL << G4endl; 01675 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0004", 01676 JustWarning, ed); 01677 } 01678 ResSec->SetDefinition( theDefinition ); 01679 SecSec->SetDefinition( aGamma ); 01680 } 01681 else 01682 { 01683 delete ResSec; 01684 delete SecSec; 01685 complete=1; 01686 } 01687 FstSec->SetDefinition( aTriton ); 01688 mF=mTrit; // First hadron 4-momentum 01689 break; 01690 case 5: // gamma+He3 01691 rA-=3; 01692 rZ-=2; 01693 if(!evaporate || rA<2) 01694 { 01695 if(!rZ) theDefinition=aNeutron; 01696 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01697 if(!theDefinition) 01698 { 01699 // G4cerr<<"-Warning-G4LE::PSDI: notDef(5), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01700 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01701 G4ExceptionDescription ed; 01702 ed << "Particle definition is a null pointer: notDef(5), Z= " << rZ 01703 << ", A=" << rA << ", L=" << rL << G4endl; 01704 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0005", 01705 JustWarning, ed); 01706 } 01707 ResSec->SetDefinition( theDefinition ); 01708 SecSec->SetDefinition( aGamma ); 01709 } 01710 else 01711 { 01712 delete ResSec; 01713 delete SecSec; 01714 complete=1; 01715 } 01716 FstSec->SetDefinition( aHe3); 01717 mF=mHel3; // First hadron 4-momentum 01718 break; 01719 case 6: 01720 rA-=4; 01721 rZ-=2; // gamma+He4 01722 if(!evaporate || rA<2) 01723 { 01724 if(!rZ) theDefinition=aNeutron; 01725 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01726 if(!theDefinition) 01727 { 01728 // G4cerr<<"-Warning-G4LE::PSDI: notDef(6), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01729 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01730 G4ExceptionDescription ed; 01731 ed << "Particle definition is a null pointer: notDef(6), Z= " << rZ 01732 << ", A=" << rA << ", L=" << rL << G4endl; 01733 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0006", 01734 JustWarning, ed); 01735 } 01736 ResSec->SetDefinition( theDefinition ); 01737 SecSec->SetDefinition( aGamma ); 01738 } 01739 else 01740 { 01741 delete ResSec; 01742 delete SecSec; 01743 complete=1; 01744 } 01745 FstSec->SetDefinition( anAlpha ); 01746 mF=mAlph; // First hadron 4-momentum 01747 break; 01748 case 7: 01749 rA-=2; // n+n 01750 if(rA==1 && !rZ) theDefinition=aNeutron; 01751 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01752 if(!theDefinition) 01753 { 01754 // G4cerr<<"-Warning-G4LE::PSDI: notDef(7), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01755 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer"); 01756 G4ExceptionDescription ed; 01757 ed << "Particle definition is a null pointer: notDef(7), Z= " << rZ 01758 << ", A=" << rA << ", L=" << rL << G4endl; 01759 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0007", 01760 JustWarning, ed); 01761 } 01762 ResSec->SetDefinition( theDefinition ); 01763 FstSec->SetDefinition( aNeutron ); 01764 SecSec->SetDefinition( aNeutron ); 01765 mF=mNeut; // First hadron 4-momentum 01766 mS=mNeut; // Second hadron 4-momentum 01767 break; 01768 case 8: 01769 rZ-=1; 01770 rA-=2; // n+p 01771 if(rA==1 && !rZ) theDefinition=aNeutron; 01772 else if(rA==2 && !rZ) 01773 { 01774 index=7; 01775 ResSec->SetDefinition( aDeuteron); 01776 FstSec->SetDefinition( aNeutron ); 01777 SecSec->SetDefinition( aNeutron ); 01778 mF=mNeut; // First hadron 4-momentum 01779 mS=mNeut; // Second hadron 4-momentum 01780 break; 01781 } 01782 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01783 if(!theDefinition) 01784 { 01785 // G4cerr<<"-Warning-G4LE::PSDI: notDef(8), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01786 // throw G4QException("G4QLowEn::PostStepDoIt: Darticle Definition is a null pointer"); 01787 G4ExceptionDescription ed; 01788 ed << "Particle definition is a null pointer: notDef(8), Z= " << rZ 01789 << ", A=" << rA << ", L=" << rL << G4endl; 01790 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0008", 01791 JustWarning, ed); 01792 } 01793 ResSec->SetDefinition( theDefinition ); 01794 FstSec->SetDefinition( aNeutron ); 01795 SecSec->SetDefinition( aProton ); 01796 mF=mNeut; // First hadron 4-momentum 01797 mS=mProt; // Second hadron 4-momentum 01798 break; 01799 case 9: 01800 rZ-=2; 01801 rA-=2; // p+p 01802 if(rA==1 && !rZ) theDefinition=aNeutron; 01803 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01804 if(!theDefinition) 01805 { 01806 // G4cerr<<"-Warning-G4LE::PSDI: notDef(9), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01807 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01808 G4ExceptionDescription ed; 01809 ed << "Particle definition is a null pointer: notDef(9), Z= " << rZ 01810 << ", A=" << rA << ", L=" << rL << G4endl; 01811 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0009", 01812 JustWarning, ed); 01813 } 01814 ResSec->SetDefinition( theDefinition ); 01815 FstSec->SetDefinition( aProton ); 01816 SecSec->SetDefinition( aProton ); 01817 mF=mProt; // First hadron 4-momentum 01818 mS=mProt; // Second hadron 4-momentum 01819 break; 01820 case 10: 01821 rZ-=2; 01822 rA-=3; // p+d 01823 if(rA==1 && !rZ) theDefinition=aNeutron; 01824 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01825 if(!theDefinition) 01826 { 01827 // G4cerr<<"-Warning-G4LE::PSDI: notDef(10), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01828 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01829 G4ExceptionDescription ed; 01830 ed << "Particle definition is a null pointer: notDef(10), Z= " << rZ 01831 << ", A=" << rA << ", L=" << rL << G4endl; 01832 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0010", 01833 JustWarning, ed); 01834 } 01835 ResSec->SetDefinition( theDefinition ); 01836 FstSec->SetDefinition( aProton ); 01837 SecSec->SetDefinition( aDeuteron ); 01838 mF=mProt; // First hadron 4-momentum 01839 mS=mDeut; // Second hadron 4-momentum 01840 break; 01841 case 11: 01842 rZ-=1; 01843 rA-=3; // n+d 01844 if(rA==1 && !rZ) theDefinition=aNeutron; 01845 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01846 if(!theDefinition) 01847 { 01848 // G4cerr<<"-Warning-G4LE::PSDI: notDef(11), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01849 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01850 G4ExceptionDescription ed; 01851 ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ 01852 << ", A=" << rA << ", L=" << rL << G4endl; 01853 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0011", 01854 JustWarning, ed); 01855 } 01856 ResSec->SetDefinition( theDefinition ); 01857 FstSec->SetDefinition( aNeutron ); 01858 SecSec->SetDefinition( aDeuteron ); 01859 mF=mNeut; // First hadron 4-momentum 01860 mS=mDeut; // Second hadron 4-momentum 01861 break; 01862 case 12: 01863 rZ-=2; 01864 rA-=4; // d+d 01865 if(rA==1 && !rZ) theDefinition=aNeutron; 01866 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01867 if(!theDefinition) 01868 { 01869 // G4cerr<<"-Warning-G4LE::PSDI: notDef(12), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01870 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01871 G4ExceptionDescription ed; 01872 ed << "Particle definition is a null pointer: notDef(12), Z= " << rZ 01873 << ", A=" << rA << ", L=" << rL << G4endl; 01874 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0012", 01875 JustWarning, ed); 01876 } 01877 ResSec->SetDefinition( theDefinition ); 01878 FstSec->SetDefinition( aDeuteron ); 01879 SecSec->SetDefinition( aDeuteron ); 01880 mF=mDeut; // First hadron 4-momentum 01881 mS=mDeut; // Second hadron 4-momentum 01882 break; 01883 case 13: 01884 rZ-=2; 01885 rA-=4; // p+t 01886 if(rA==1 && !rZ) theDefinition=aNeutron; 01887 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01888 if(!theDefinition) 01889 { 01890 // G4cerr<<"-Warning-G4LE::PSDI: notDef(13), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01891 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01892 G4ExceptionDescription ed; 01893 ed << "Particle definition is a null pointer: notDef(13), Z= " << rZ 01894 << ", A=" << rA << ", L=" << rL << G4endl; 01895 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0013", 01896 JustWarning, ed); 01897 } 01898 ResSec->SetDefinition( theDefinition ); 01899 FstSec->SetDefinition( aProton ); 01900 SecSec->SetDefinition( aTriton ); 01901 mF=mProt; // First hadron 4-momentum 01902 mS=mTrit; // Second hadron 4-momentum 01903 break; 01904 case 14: 01905 rZ-=1; 01906 rA-=4; // n+t 01907 if(rA==1 && !rZ) theDefinition=aNeutron; 01908 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01909 if(!theDefinition) 01910 { 01911 // G4cerr<<"-Warning-G4LE::PSDI: notDef(14), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01912 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01913 G4ExceptionDescription ed; 01914 ed << "Particle definition is a null pointer: notDef(14), Z= " << rZ 01915 << ", A=" << rA << ", L=" << rL << G4endl; 01916 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0014", 01917 JustWarning, ed); 01918 } 01919 ResSec->SetDefinition( theDefinition ); 01920 FstSec->SetDefinition( aNeutron ); 01921 SecSec->SetDefinition( aTriton ); 01922 mF=mNeut; // First hadron 4-momentum 01923 mS=mTrit; // Second hadron 4-momentum 01924 break; 01925 case 15: 01926 rZ-=3; 01927 rA-=4; // p+He3 01928 if(rA==1 && !rZ) theDefinition=aNeutron; 01929 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01930 if(!theDefinition) 01931 { 01932 // G4cerr<<"-Warning-G4LE::PSDI: notDef(15), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01933 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01934 G4ExceptionDescription ed; 01935 ed << "Particle definition is a null pointer: notDef(15), Z= " << rZ 01936 << ", A=" << rA << ", L=" << rL << G4endl; 01937 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0015", 01938 JustWarning, ed); 01939 } 01940 ResSec->SetDefinition( theDefinition ); 01941 FstSec->SetDefinition( aProton); 01942 SecSec->SetDefinition( aHe3 ); 01943 mF=mProt; // First hadron 4-momentum 01944 mS=mHel3; // Second hadron 4-momentum 01945 break; 01946 case 16: 01947 rZ-=2; 01948 rA-=4; // n+He3 01949 if(rA==1 && !rZ) theDefinition=aNeutron; 01950 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01951 if(!theDefinition) 01952 { 01953 // G4cerr<<"-Warning-G4LE::PSDI: notDef(16), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01954 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01955 G4ExceptionDescription ed; 01956 ed << "Particle definition is a null pointer: notDef(16), Z= " << rZ 01957 << ", A=" << rA << ", L=" << rL << G4endl; 01958 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0016", 01959 JustWarning, ed); 01960 } 01961 ResSec->SetDefinition( theDefinition ); 01962 FstSec->SetDefinition( aNeutron ); 01963 SecSec->SetDefinition( aHe3 ); 01964 mF=mNeut; // First hadron 4-momentum 01965 mS=mHel3; // Second hadron 4-momentum 01966 break; 01967 case 17: 01968 rZ-=3; 01969 rA-=5; // p+alph 01970 if(rA==1 && !rZ) theDefinition=aNeutron; 01971 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01972 if(!theDefinition) 01973 { 01974 // G4cerr<<"-Warning-G4LE::PSDI: notDef(17), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01975 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01976 G4ExceptionDescription ed; 01977 ed << "Particle definition is a null pointer: notDef(17), Z= " << rZ 01978 << ", A=" << rA << ", L=" << rL << G4endl; 01979 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0017", 01980 JustWarning, ed); 01981 } 01982 ResSec->SetDefinition( theDefinition ); 01983 FstSec->SetDefinition( aProton ); 01984 SecSec->SetDefinition( anAlpha ); 01985 mF=mProt; // First hadron 4-momentum 01986 mS=mAlph; // Second hadron 4-momentum 01987 break; 01988 case 18: 01989 rZ-=2; 01990 rA-=5; // n+alph 01991 if(rA==1 && !rZ) theDefinition=aNeutron; 01992 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 01993 if(!theDefinition) 01994 { 01995 // G4cerr<<"-Warning-G4LE::PSDI: notDef(18), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 01996 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 01997 G4ExceptionDescription ed; 01998 ed << "Particle definition is a null pointer: notDef(18), Z= " << rZ 01999 << ", A=" << rA << ", L=" << rL << G4endl; 02000 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0018", 02001 JustWarning, ed); 02002 } 02003 ResSec->SetDefinition( theDefinition ); 02004 FstSec->SetDefinition( aNeutron ); 02005 SecSec->SetDefinition( anAlpha ); 02006 mF=mNeut; // First hadron 4-momentum 02007 mS=mAlph; // Second hadron 4-momentum 02008 break; 02009 case 19: 02010 rL-=1; // L+gamma (@@ rA inludes rL?) 02011 rA-=1; 02012 if(rA==1 && !rZ) theDefinition=aNeutron; 02013 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02014 if(!theDefinition) 02015 { 02016 // G4cerr<<"-Warning-G4LE::PSDI: notDef(19), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02017 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 02018 G4ExceptionDescription ed; 02019 ed << "Particle definition is a null pointer: notDef(19), Z= " << rZ 02020 << ", A=" << rA << ", L=" << rL << G4endl; 02021 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0019", 02022 JustWarning, ed); 02023 } 02024 ResSec->SetDefinition( theDefinition ); 02025 FstSec->SetDefinition( aLambda ); 02026 SecSec->SetDefinition( aGamma ); 02027 mF=mLamb; // First hadron 4-momentum 02028 break; 02029 case 20: 02030 rL-=1; // L+p (@@ rA inludes rL?) 02031 rZ-=1; 02032 rA-=2; 02033 if(rA==1 && !rZ) theDefinition=aNeutron; 02034 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02035 if(!theDefinition) 02036 { 02037 // G4cerr<<"-Warning-G4LE::PSDI: notDef(20), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02038 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 02039 G4ExceptionDescription ed; 02040 ed << "Particle definition is a null pointer: notDef(20), Z= " << rZ 02041 << ", A=" << rA << ", L=" << rL << G4endl; 02042 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0020", 02043 JustWarning, ed); 02044 } 02045 ResSec->SetDefinition( theDefinition ); 02046 FstSec->SetDefinition( aProton ); 02047 SecSec->SetDefinition( aLambda ); 02048 mF=mProt; // First hadron 4-momentum 02049 mS=mLamb; // Second hadron 4-momentum 02050 break; 02051 case 21: 02052 rL-=1; // L+n (@@ rA inludes rL?) 02053 rA-=2; 02054 if(rA==1 && !rZ) theDefinition=aNeutron; 02055 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02056 if(!theDefinition) 02057 { 02058 // G4cerr<<"-Warning-G4LE::PSDI: notDef(21), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02059 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Pefinition is a null pointer"); 02060 G4ExceptionDescription ed; 02061 ed << "Particle definition is a null pointer: notDef(21), Z= " << rZ 02062 << ", A=" << rA << ", L=" << rL << G4endl; 02063 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0021", 02064 JustWarning, ed); 02065 } 02066 ResSec->SetDefinition( theDefinition ); 02067 FstSec->SetDefinition( aNeutron ); 02068 SecSec->SetDefinition( aLambda ); 02069 mF=mNeut; // First hadron 4-momentum 02070 mS=mLamb; // Second hadron 4-momentum 02071 break; 02072 case 22: 02073 rL-=1; // L+d (@@ rA inludes rL?) 02074 rZ-=1; 02075 rA-=3; 02076 if(rA==1 && !rZ) theDefinition=aNeutron; 02077 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02078 if(!theDefinition) 02079 { 02080 // G4cerr<<"-Warning-G4LE::PSDI: notDef(22), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02081 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 02082 G4ExceptionDescription ed; 02083 ed << "Particle definition is a null pointer: notDef(22), Z= " << rZ 02084 << ", A=" << rA << ", L=" << rL << G4endl; 02085 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0022", 02086 JustWarning, ed); 02087 } 02088 ResSec->SetDefinition( theDefinition ); 02089 FstSec->SetDefinition( aDeuteron ); 02090 SecSec->SetDefinition( aLambda ); 02091 mF=mDeut; // First hadron 4-momentum 02092 mS=mLamb; // Second hadron 4-momentum 02093 break; 02094 case 23: 02095 rL-=1; // L+t (@@ rA inludes rL?) 02096 rZ-=1; 02097 rA-=4; 02098 if(rA==1 && !rZ) theDefinition=aNeutron; 02099 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02100 if(!theDefinition) 02101 { 02102 // G4cerr<<"-Warning-G4LE::PSDI: notDef(23), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02103 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 02104 G4ExceptionDescription ed; 02105 ed << "Particle definition is a null pointer: notDef(23), Z= " << rZ 02106 << ", A=" << rA << ", L=" << rL << G4endl; 02107 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0023", 02108 JustWarning, ed); 02109 } 02110 ResSec->SetDefinition( theDefinition ); 02111 FstSec->SetDefinition( aTriton ); 02112 SecSec->SetDefinition( aLambda ); 02113 mF=mTrit; // First hadron 4-momentum 02114 mS=mLamb; // Second hadron 4-momentum 02115 break; 02116 case 24: 02117 rL-=1; // L+He3 (@@ rA inludes rL?) 02118 rZ-=2; 02119 rA-=4; 02120 if(rA==1 && !rZ) theDefinition=aNeutron; 02121 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02122 if(!theDefinition) 02123 { 02124 // G4cerr<<"-Warning-G4LE::PSDI: notDef(24), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02125 // throw G4QException("G4QLowEn::PostStepDoIt: particle definition is a null pointer"); 02126 G4ExceptionDescription ed; 02127 ed << "Particle definition is a null pointer: notDef(24), Z= " << rZ 02128 << ", A=" << rA << ", L=" << rL << G4endl; 02129 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0024", 02130 JustWarning, ed); 02131 } 02132 ResSec->SetDefinition( theDefinition ); 02133 FstSec->SetDefinition( aHe3 ); 02134 SecSec->SetDefinition( aLambda ); 02135 mF=mHel3; // First hadron 4-momentum 02136 mS=mLamb; // Second hadron 4-momentum 02137 break; 02138 case 25: 02139 rL-=1; // L+alph (@@ rA inludes rL?) 02140 rZ-=2; 02141 rA-=5; 02142 if(rA==1 && !rZ) theDefinition=aNeutron; 02143 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); 02144 if(!theDefinition) 02145 { 02146 // G4cerr<<"-Warning-G4LE::PSDI: notDef(25), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; 02147 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer"); 02148 G4ExceptionDescription ed; 02149 ed << "Particle definition is a null pointer: notDef(25), Z= " << rZ 02150 << ", A=" << rA << ", L=" << rL << G4endl; 02151 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0025", 02152 JustWarning, ed); 02153 } 02154 ResSec->SetDefinition( theDefinition ); 02155 FstSec->SetDefinition( anAlpha ); 02156 SecSec->SetDefinition( aLambda ); 02157 mF=mAlph; // First hadron 4-momentum 02158 mS=mLamb; // Second hadron 4-momentum 02159 break; 02160 } 02161 #ifdef debug 02162 G4cout<<"G4QLowEn::PSDI:F="<<mF<<",S="<<mS<<",com="<<complete<<",ev="<<evaporate<<G4endl; 02163 #endif 02164 G4LorentzVector fst4Mom(zeroMom,mF); // Prototype of the first hadron 4-momentum 02165 G4LorentzVector snd4Mom(zeroMom,mS); // Prototype of the second hadron 4-momentum 02166 G4LorentzVector dir4Mom=tot4M; // Prototype of the resN decay direction 4-momentum 02167 dir4Mom.setE(tot4M.e()/2.); // Get half energy and total 3-momentum 02168 // @@ Can be repeated to take into account the Coulomb Barrier 02169 if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp)) 02170 { 02171 // G4cerr<<"**G4LowEnergy::PoStDoIt:i="<<index<<",tM="<<totM<<"->M1="<<res4Mom.m()<<"+M2=" 02172 // <<fst4Mom.m()<<"+M3="<<snd4Mom.m()<<"=="<<res4Mom.m()+fst4Mom.m()+snd4Mom.m()<<G4endl; 02173 // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound"); 02174 G4ExceptionDescription ed; 02175 ed << "Can't decay the Compound: i=" << index << ",tM=" << totM << "->M1=" 02176 << res4Mom.m() << "+M2=" << fst4Mom.m() << "+M3=" << snd4Mom.m() << "==" 02177 << res4Mom.m()+fst4Mom.m()+snd4Mom.m() << G4endl; 02178 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000", 02179 FatalException, ed); 02180 } 02181 #ifdef debug 02182 G4cout<<"G4QLowEn::PSDI:r4M="<<res4Mom<<",f4M="<<fst4Mom<<",s4M="<<snd4Mom<<G4endl; 02183 #endif 02184 G4Track* aNewTrack = 0; 02185 if(complete) 02186 { 02187 FstSec->Set4Momentum(fst4Mom); 02188 aNewTrack = new G4Track(FstSec, localtime, position ); 02189 aNewTrack->SetWeight(weight); // weighted 02190 aNewTrack->SetTouchableHandle(trTouchable); 02191 result.push_back( aNewTrack ); 02192 EnMomConservation-=fst4Mom; 02193 #ifdef debug 02194 G4cout<<"G4QLowEn::PSDI: ***Filled*** 1stH4M="<<fst4Mom 02195 <<", PDG="<<FstSec->GetDefinition()->GetPDGEncoding()<<G4endl; 02196 #endif 02197 if(complete>2) // Final solution 02198 { 02199 ResSec->Set4Momentum(res4Mom); 02200 aNewTrack = new G4Track(ResSec, localtime, position ); 02201 aNewTrack->SetWeight(weight); // weighted 02202 aNewTrack->SetTouchableHandle(trTouchable); 02203 result.push_back( aNewTrack ); 02204 EnMomConservation-=res4Mom; 02205 #ifdef debug 02206 G4cout<<"G4QLowEn::PSDI: ***Filled*** rA4M="<<res4Mom<<",rZ="<<rZ<<",rA="<<rA<<",rL=" 02207 <<rL<<G4endl; 02208 #endif 02209 SecSec->Set4Momentum(snd4Mom); 02210 aNewTrack = new G4Track(SecSec, localtime, position ); 02211 aNewTrack->SetWeight(weight); // weighted 02212 aNewTrack->SetTouchableHandle(trTouchable); 02213 result.push_back( aNewTrack ); 02214 EnMomConservation-=snd4Mom; 02215 #ifdef debug 02216 G4cout<<"G4QLowEn::PSDI: ***Filled*** 2ndH4M="<<snd4Mom 02217 <<", PDG="<<SecSec->GetDefinition()->GetPDGEncoding()<<G4endl; 02218 #endif 02219 } 02220 else res4Mom+=snd4Mom; 02221 } 02222 else res4Mom=tot4M; 02223 if(complete<3) // Evaporation of the residual must be done 02224 { 02225 G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom); // Input hadron-nucleus 02226 G4QHadronVector* evaHV = new G4QHadronVector; // Output vector of hadrons (delete!) 02227 Nuc.EvaporateNucleus(rHadron, evaHV); // here a pion can appear ! 02228 G4int nOut=evaHV->size(); 02229 for(i=0; i<nOut; i++) 02230 { 02231 G4QHadron* curH = (*evaHV)[i]; 02232 G4int hPDG=curH->GetPDGCode(); 02233 G4LorentzVector h4Mom=curH->Get4Momentum(); 02234 EnMomConservation-=h4Mom; 02235 #ifdef debug 02236 G4cout<<"G4QLowEn::PSDI: ***FillingCand#"<<i<<"*** evaH="<<hPDG<<h4Mom<<G4endl; 02237 #endif 02238 if (hPDG==90000001 || hPDG==2112) theDefinition = aNeutron; 02239 else if(hPDG==90001000 || hPDG==2212) theDefinition = aProton; 02240 else if(hPDG==91000000 || hPDG==3122) theDefinition = aLambda; 02241 else if(hPDG== 22 ) theDefinition = aGamma; 02242 else if(hPDG== 111) theDefinition = aPiZero; 02243 else if(hPDG== 211) theDefinition = aPiPlus; 02244 else if(hPDG== -211) theDefinition = aPiMinus; 02245 else 02246 { 02247 G4int hZ=curH->GetCharge(); 02248 G4int hA=curH->GetBaryonNumber(); 02249 G4int hS=curH->GetStrangeness(); 02250 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(hZ,hA,hS,0); // ion 02251 } 02252 if(theDefinition) 02253 { 02254 G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom); 02255 G4Track* evaQH = new G4Track(theEQH, localtime, position ); 02256 evaQH->SetWeight(weight); // weighted 02257 evaQH->SetTouchableHandle(trTouchable); 02258 result.push_back( evaQH ); 02259 } 02260 else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl; 02261 } 02262 } 02263 // algorithm implementation --- STOPS HERE 02264 G4int nres=result.size(); 02265 aParticleChange.SetNumberOfSecondaries(nres); 02266 for(i=0; i<nres; ++i) aParticleChange.AddSecondary(result[i]); 02267 #ifdef debug 02268 G4cout<<"G4QLowEnergy::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl; 02269 #endif 02270 return G4VDiscreteProcess::PostStepDoIt(track, step); 02271 }
void G4QLowEnergy::SwitchOffEvaporation | ( | ) | [inline] |
void G4QLowEnergy::SwitchOnEvaporation | ( | ) | [inline] |