00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4QLowEnergy.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4HadronicDeprecate.hh"
00048
00049
00050
00051
00052
00053 G4int G4QLowEnergy::nPartCWorld=85;
00054 std::vector<G4int> G4QLowEnergy::ElementZ;
00055 std::vector<G4double> G4QLowEnergy::ElProbInMat;
00056 std::vector<std::vector<G4int>*> G4QLowEnergy::ElIsoN;
00057 std::vector<std::vector<G4double>*>G4QLowEnergy::IsoProbInEl;
00058
00059
00060 G4QLowEnergy::G4QLowEnergy(const G4String& processName):
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);
00070 }
00071
00072
00073 G4QLowEnergy::~G4QLowEnergy() {}
00074
00075
00076 G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation(){return EnMomConservation;}
00077
00078 G4int G4QLowEnergy::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
00079
00080
00081
00082
00083 G4double G4QLowEnergy::GetMeanFreePath(const G4Track&Track, G4double, G4ForceCondition*F)
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
00091 G4double Momentum = incidentParticle->GetTotalMomentum();
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();
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();
00122 Momentum/=incidentParticleDefinition->GetBaryonNumber();
00123 G4QIsotope* Isotopes = G4QIsotope::Get();
00124 G4double sigma=0.;
00125 G4int IPIE=IsoProbInEl.size();
00126 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)
00127 {
00128 std::vector<G4double>* SPI=IsoProbInEl[ip];
00129 SPI->clear();
00130 delete SPI;
00131 std::vector<G4int>* IsN=ElIsoN[ip];
00132 IsN->clear();
00133 delete IsN;
00134 }
00135 ElProbInMat.clear();
00136 ElementZ.clear();
00137 IsoProbInEl.clear();
00138 ElIsoN.clear();
00139 for(G4int i=0; i<nE; ++i)
00140 {
00141 G4Element* pElement=(*theElementVector)[i];
00142 Z = static_cast<G4int>(pElement->GetZ());
00143 ElementZ.push_back(Z);
00144 G4int isoSize=0;
00145 G4int indEl=0;
00146 G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
00147 if(isoVector) isoSize=isoVector->size();
00148 #ifdef debug
00149 G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
00150 #endif
00151 if(isoSize)
00152 {
00153 indEl=pElement->GetIndex()+1;
00154 #ifdef debug
00155 G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
00156 #endif
00157 if(!Isotopes->IsDefined(Z,indEl))
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++)
00163 {
00164 G4int N=pElement->GetIsotope(j)->GetN()-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);
00178 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
00179 delete newAbund;
00180 }
00181 }
00182 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);
00183 std::vector<G4double>* SPI = new std::vector<G4double>;
00184 IsoProbInEl.push_back(SPI);
00185 std::vector<G4int>* IsN = new std::vector<G4int>;
00186 ElIsoN.push_back(IsN);
00187 G4int nIs=cs->size();
00188 #ifdef debug
00189 G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
00190 #endif
00191 G4double susi=0.;
00192 if(nIs) for(G4int j=0; j<nIs; j++)
00193 {
00194 std::pair<G4int,G4double>* curIs=(*cs)[j];
00195 G4int N=curIs->first;
00196 IsN->push_back(N);
00197 #ifdef debug
00198 G4cout<<"G4QLoE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl;
00199 #endif
00200 G4bool ccsf=true;
00201 #ifdef debug
00202 G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl;
00203 #endif
00204 G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);
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;
00211 SPI->push_back(susi);
00212 }
00213 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];
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 }
00220
00221 #ifdef debug
00222 G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
00223 #endif
00224 if(sigma > 0.000000001) return 1./sigma;
00225 return DBL_MAX;
00226 }
00227
00228 G4bool G4QLowEnergy::IsApplicable(const G4ParticleDefinition& particle)
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 }
00246
00247 G4VParticleChange* G4QLowEnergy::PostStepDoIt(const G4Track& track, const G4Step& step)
00248 {
00249 static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV;
00250 static const G4double mPro2= mProt*mProt;
00251 static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV;
00252 static const G4double mNeu2= mNeut*mNeut;
00253 static const G4double mLamb= G4QPDGCode(3122).GetMass()/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;
00273 static G4QNucleus Nuc;
00274
00275
00276 static G4bool CWinit = true;
00277 if(CWinit)
00278 {
00279 CWinit=false;
00280 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
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
00291
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;
00297 G4double momentum = projHadron->GetTotalMomentum()/MeV;
00298 G4double Momentum = proj4M.rho();
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))
00306 {
00307 G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl;
00308 return 0;
00309 }
00310 const G4Material* material = track.GetMaterial();
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;
00317
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;
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
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;
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];
00373 std::vector<G4int>* IsN = ElIsoN[i];
00374 G4int nofIsot=SPI->size();
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();
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;
00390 }
00391 G4int tN =(*IsN)[j]; ;
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;
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;
00424 G4QPDGCode targQPDG(targPDG);
00425 G4double tM=targQPDG.GetMass();
00426 G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3);
00427 G4int pZ=static_cast<G4int>(particle->GetPDGCharge());
00428 G4int pN=particle->GetBaryonNumber()-pZ-pL;
00429 G4double pM=targQPDG.GetNuclMass(pZ,pN,0);
00430 G4double cosp=-14*Momentum*(tM-pM)/tM/pM;
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;
00436 G4ParticleMomentum dir = projHadron->GetMomentumDirection();
00437 G4LorentzVector targ4M=G4LorentzVector(0.,0.,0.,tM);
00438 G4LorentzVector tot4M =proj4M+targ4M;
00439 #ifdef pdebug
00440 G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
00441 #endif
00442 EnMomConservation=tot4M;
00443
00444 G4double dER = tot4M.e() - tM - pM;
00445 G4double pA = pZ+pN;
00446 G4double tA = tZ+tN;
00447 G4double CBE = Nuc.CoulombBarGen(tZ, tA, pZ, pA);
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)
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
00459 aParticleChange.ProposeEnergy(kinEnergy);
00460 aParticleChange.ProposeLocalEnergyDeposit(0.);
00461 aParticleChange.ProposeMomentumDirection(dir) ;
00462 return G4VDiscreteProcess::PostStepDoIt(track,step);
00463 }
00464
00465
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);
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
00479 if(xSec <= 0.)
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
00486 aParticleChange.ProposeEnergy(kinEnergy);
00487 aParticleChange.ProposeLocalEnergyDeposit(0.);
00488 aParticleChange.ProposeMomentumDirection(dir) ;
00489 return G4VDiscreteProcess::PostStepDoIt(track,step);
00490 }
00491
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
00500 tA=tB;
00501 pA=pB;
00502 G4double tR=1.1;
00503 if(tB > 1) tR*=std::pow(tA,third);
00504 G4double pR=1.1;
00505 if(pB > 1) pR*=std::pow(pA,third);
00506 G4double R=tR+pR;
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;
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))
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
00547 G4double DpR3=DpR2*DpR;
00548
00549 G4double DpR4=DpR3*DpR;
00550 G4double tR2=tR*tR;
00551 G4double pR2=pR*pR;
00552 G4double tR3=tR2*tR;
00553
00554 G4double tR4=tR3*tR;
00555
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
00562
00563
00564 if(std::fabs(tF-4.) < 1.) tD=4;
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
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;
00577 G4int tC=0;
00578 if((tD || pD) && tD<tB && pD<pB)
00579 {
00580 if(!tD || !pD)
00581 {
00582 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
00583 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
00584 G4int pPDG=2112;
00585 G4double prM =mNeut;
00586 G4double prM2=mNeu2;
00587 if (!tD)
00588 {
00589 if(pD > 1) pD=1;
00590 if(!pN || (pZ && pA*G4UniformRand() < pZ) )
00591 {
00592 pPDG=2212;
00593 prM=mProt;
00594 prM2=mPro2;
00595 pC=1;
00596 }
00597 G4double Mom=0.;
00598 G4LorentzVector com4M = targ4M;
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;
00606 rNM=prM;
00607 }
00608 else
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;
00617 com4M += qfN4M;
00618 G4double pNE = qfN4M.e();
00619 if(rNE >= prM) Mom = std::sqrt(pNE*pNE-prM2);
00620 else break;
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;
00626 G4double mint=0.;
00627 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG);
00628 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG);
00629 G4double maxt=0.;
00630 if(pPDG==2212) maxt=PCSmanager->GetHMaxT();
00631 else maxt=NCSmanager->GetHMaxT();
00632 G4double cost=1.-mint/maxt;
00633 if(cost>1. || cost<-1.) break;
00634 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tgM);
00635 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM);
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;
00641 }
00642 G4Track* projSpect = 0;
00643 G4Track* aNucTrack = 0;
00644 if(pB > pD)
00645 {
00646 G4int rZ=pZ-pC;
00647 G4int rA=pB-1;
00648 G4ParticleDefinition* theDefinition;
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);
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);
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)
00681 {
00682 theDefinition = anAlpha;
00683 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph);
00684 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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))
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);
00709 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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);
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
00744 {
00745 if(tD > 1) tD=1;
00746 if(!tN || (tZ && tA*G4UniformRand() < tZ) )
00747 {
00748 pPDG=2212;
00749 prM=mProt;
00750 prM2=mPro2;
00751 tC=1;
00752 }
00753 G4double Mom=0.;
00754 G4LorentzVector com4M=proj4M;
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;
00762 rNM=prM;
00763 }
00764 else
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;
00771 com4M += qfN4M;
00772 G4ThreeVector boostV=proj4M.vect()/proj4M.e();
00773 qfN4M.boost(-boostV);
00774 G4double tNE = qfN4M.e();
00775 if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2);
00776 else break;
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;
00782 G4double mint=0.;
00783 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG);
00784 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG);
00785 G4double maxt=0.;
00786 if(pPDG==2212) maxt=PCSmanager->GetHMaxT();
00787 else maxt=NCSmanager->GetHMaxT();
00788 G4double cost=1.-mint/maxt;
00789 if(cost>1. || cost<-1.) break;
00790 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,prM);
00791 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM);
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;
00797 }
00798 G4Track* targSpect = 0;
00799 G4Track* aNucTrack = 0;
00800 if(tB > tD)
00801 {
00802 G4int rZ=tZ-tC;
00803 G4int rA=tB-1;
00804 G4ParticleDefinition* theDefinition;
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);
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);
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)
00837 {
00838 theDefinition = anAlpha;
00839 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph);
00840 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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))
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);
00865 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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);
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);
00903 }
00904 else
00905 {
00906
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
00930 G4ThreeVector pFM(0.,0.,0.);
00931 if(pD < pB)
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)
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
00952 G4int rpZ=pZ-pC;
00953 G4int pF_value=pD-pC;
00954 G4int rpN=pN-pF_value;
00955 G4double rpNM=G4QPDGCode(2112).GetNuclMass(rpZ, rpN, 0);
00956 G4ThreeVector boostV=proj4M.vect()/proj4M.e();
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;
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)
00979 {
00980 theDefinition = aNeutron;
00981 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut);
00982 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut);
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);
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)
01006 {
01007 theDefinition = aNeutron;
01008 G4LorentzVector f4M=rp4M/rpA;
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)
01013 {
01014 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
01015 G4Track* aNTrack = new G4Track(pNeu, localtime, position );
01016 aNTrack->SetWeight(weight);
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)
01033 {
01034 theDefinition = anAlpha;
01035 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph);
01036 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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))
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);
01070 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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
01099
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);
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
01130 G4int rtZ=tZ-tC;
01131 G4int tF_value=tD-tC;
01132 G4int rtN=tN-tF_value;
01133 G4double rtNM=G4QPDGCode(2112).GetNuclMass(rtZ, rtN, 0);
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)
01148 {
01149 theDefinition = aNeutron;
01150 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut);
01151 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut);
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);
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)
01173 {
01174 theDefinition = aNeutron;
01175 G4LorentzVector f4M=rt4M/rtA;
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)
01180 {
01181 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
01182 G4Track* aNTrack = new G4Track(pNeu, localtime, position );
01183 aNTrack->SetWeight(weight);
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)
01200 {
01201 theDefinition = anAlpha;
01202 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph);
01203 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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))
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);
01237 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph);
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);
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
01266
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);
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();
01300 G4int sN=tF_value+pF_value;
01301 G4int sZ=tC+pC;
01302 tnM = targQPDG.GetNuclMass(sZ,sN,0);
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)
01308 {
01309 pZ=pC;
01310 pN=pF_value;
01311 tZ=tC;
01312 tN=tF_value;
01313 tot4M=tt4M;
01314 }
01315
01316
01317
01318
01319
01320 }
01321 }
01322 else if(tD==tB || pD==pB)
01323 {
01324 tD=tZ+tN;
01325 pD=pZ+pN;
01326 tcM=tnM+1.;
01327 }
01328 }
01329 else
01330 {
01331 tD=tZ+tN;
01332 pD=pZ+pN;
01333 tcM=tnM+1.;
01334 }
01335 }
01336 G4double totM=tot4M.m();
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
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);
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);
01353 mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0);
01354 }
01355 if ( totZ > 1 && totN > 1 ) mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0);
01356 if ( totZ > 1 && totN > 2 ) mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0);
01357 if ( totZ > 2 && totN > 1 ) mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0);
01358 if ( totZ > 2 && totN > 2 ) mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0);
01359 if ( totZ > 0 && totN > 2 ) mass[7] = targQPDG.GetNuclMass(totZ ,totN-2,0);
01360 mass[ 8] = mass[3];
01361 if ( totZ > 2 ) mass[9] = targQPDG.GetNuclMass(totZ-2,totN ,0);
01362 mass[10] = mass[5];
01363 mass[11] = mass[4];
01364 mass[12] = mass[6];
01365 mass[13] = mass[6];
01366 if ( totZ > 1 && totN > 3 ) mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0);
01367 if ( totZ > 3 && totN > 1 ) mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0);
01368 mass[16] = mass[6];
01369 if ( totZ > 3 && totN > 2 ) mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0);
01370 if ( totZ > 2 && totN > 3 ) mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0);
01371 if(pL>0)
01372 {
01373 G4int pL1=pL-1;
01374 if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ ,totN ,pL1);
01375 if( (totN > 0 && totZ > 0) || totZ > 1 )
01376 mass[20]=targQPDG.GetNuclMass(totZ-1,totN ,pL1);
01377 if( (totN > 0 && totZ > 0) || totN > 0 )
01378 mass[21]=targQPDG.GetNuclMass(totZ ,totN-1,pL1);
01379 if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) )
01380 mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);
01381 if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) )
01382 mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);
01383 if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) )
01384 mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);
01385 if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) )
01386 mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);
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;
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;
01444 if( qval[i] < 0.0 ) qval[i] = 0.0;
01445 qv += qval[i];
01446 }
01447
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
01466 else {
01467 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01468 FatalException, "Can't decay ResidualCompound");
01469 }
01470 G4double minEx=1000000.;
01471 G4bool found = false;
01472 G4int cInd = 0;
01473 G4int ctN = 0;
01474 G4int ctZ = 0;
01475 G4LorentzVector c4M=tot4M;
01476 G4double ctM=0.;
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();
01486 G4int iC = G4int(iParticle->GetPDGCharge());
01487 G4LorentzVector new4M = tot4M + ih4M;
01488 G4int ntN=totN + iB-iC;
01489 G4int ntZ=totZ + iC;
01490 G4double ntM = targQPDG.GetNuclMass(ntZ, ntN, 0);
01491 G4double ttM = new4M.m();
01492 if(ttM > ntM)
01493 {
01494 G4double nEx = ttM - ntM;
01495 if(found && nEx < minEx)
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
01519 index = 0;
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
01526 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
01527 FatalException, "Can't correct by one particle");
01528 }
01529 }
01530
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);
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;
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
01556
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;
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
01584
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;
01602 break;
01603 case 2:
01604 rA-=1;
01605 rZ-=1;
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
01613
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;
01631 break;
01632 case 3:
01633 rA-=2;
01634 rZ-=1;
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
01642
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;
01660 break;
01661 case 4:
01662 rA-=3;
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
01671
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;
01689 break;
01690 case 5:
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
01700
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;
01718 break;
01719 case 6:
01720 rA-=4;
01721 rZ-=2;
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
01729
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;
01747 break;
01748 case 7:
01749 rA-=2;
01750 if(rA==1 && !rZ) theDefinition=aNeutron;
01751 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01752 if(!theDefinition)
01753 {
01754
01755
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;
01766 mS=mNeut;
01767 break;
01768 case 8:
01769 rZ-=1;
01770 rA-=2;
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;
01779 mS=mNeut;
01780 break;
01781 }
01782 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01783 if(!theDefinition)
01784 {
01785
01786
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;
01797 mS=mProt;
01798 break;
01799 case 9:
01800 rZ-=2;
01801 rA-=2;
01802 if(rA==1 && !rZ) theDefinition=aNeutron;
01803 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01804 if(!theDefinition)
01805 {
01806
01807
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;
01818 mS=mProt;
01819 break;
01820 case 10:
01821 rZ-=2;
01822 rA-=3;
01823 if(rA==1 && !rZ) theDefinition=aNeutron;
01824 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01825 if(!theDefinition)
01826 {
01827
01828
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;
01839 mS=mDeut;
01840 break;
01841 case 11:
01842 rZ-=1;
01843 rA-=3;
01844 if(rA==1 && !rZ) theDefinition=aNeutron;
01845 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01846 if(!theDefinition)
01847 {
01848
01849
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;
01860 mS=mDeut;
01861 break;
01862 case 12:
01863 rZ-=2;
01864 rA-=4;
01865 if(rA==1 && !rZ) theDefinition=aNeutron;
01866 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01867 if(!theDefinition)
01868 {
01869
01870
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;
01881 mS=mDeut;
01882 break;
01883 case 13:
01884 rZ-=2;
01885 rA-=4;
01886 if(rA==1 && !rZ) theDefinition=aNeutron;
01887 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01888 if(!theDefinition)
01889 {
01890
01891
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;
01902 mS=mTrit;
01903 break;
01904 case 14:
01905 rZ-=1;
01906 rA-=4;
01907 if(rA==1 && !rZ) theDefinition=aNeutron;
01908 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01909 if(!theDefinition)
01910 {
01911
01912
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;
01923 mS=mTrit;
01924 break;
01925 case 15:
01926 rZ-=3;
01927 rA-=4;
01928 if(rA==1 && !rZ) theDefinition=aNeutron;
01929 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01930 if(!theDefinition)
01931 {
01932
01933
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;
01944 mS=mHel3;
01945 break;
01946 case 16:
01947 rZ-=2;
01948 rA-=4;
01949 if(rA==1 && !rZ) theDefinition=aNeutron;
01950 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01951 if(!theDefinition)
01952 {
01953
01954
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;
01965 mS=mHel3;
01966 break;
01967 case 17:
01968 rZ-=3;
01969 rA-=5;
01970 if(rA==1 && !rZ) theDefinition=aNeutron;
01971 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01972 if(!theDefinition)
01973 {
01974
01975
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;
01986 mS=mAlph;
01987 break;
01988 case 18:
01989 rZ-=2;
01990 rA-=5;
01991 if(rA==1 && !rZ) theDefinition=aNeutron;
01992 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
01993 if(!theDefinition)
01994 {
01995
01996
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;
02007 mS=mAlph;
02008 break;
02009 case 19:
02010 rL-=1;
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
02017
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;
02028 break;
02029 case 20:
02030 rL-=1;
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
02038
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;
02049 mS=mLamb;
02050 break;
02051 case 21:
02052 rL-=1;
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
02059
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;
02070 mS=mLamb;
02071 break;
02072 case 22:
02073 rL-=1;
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
02081
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;
02092 mS=mLamb;
02093 break;
02094 case 23:
02095 rL-=1;
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
02103
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;
02114 mS=mLamb;
02115 break;
02116 case 24:
02117 rL-=1;
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
02125
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;
02136 mS=mLamb;
02137 break;
02138 case 25:
02139 rL-=1;
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
02147
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;
02158 mS=mLamb;
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);
02165 G4LorentzVector snd4Mom(zeroMom,mS);
02166 G4LorentzVector dir4Mom=tot4M;
02167 dir4Mom.setE(tot4M.e()/2.);
02168
02169 if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp))
02170 {
02171
02172
02173
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);
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)
02198 {
02199 ResSec->Set4Momentum(res4Mom);
02200 aNewTrack = new G4Track(ResSec, localtime, position );
02201 aNewTrack->SetWeight(weight);
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);
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)
02224 {
02225 G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom);
02226 G4QHadronVector* evaHV = new G4QHadronVector;
02227 Nuc.EvaporateNucleus(rHadron, evaHV);
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);
02251 }
02252 if(theDefinition)
02253 {
02254 G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom);
02255 G4Track* evaQH = new G4Track(theEQH, localtime, position );
02256 evaQH->SetWeight(weight);
02257 evaQH->SetTouchableHandle(trTouchable);
02258 result.push_back( evaQH );
02259 }
02260 else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl;
02261 }
02262 }
02263
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 }
02272
02273 G4double G4QLowEnergy::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG)
02274 {
02275 static G4bool first=true;
02276 static G4VQCrossSection* CSmanager;
02277 if(first)
02278 {
02279 CSmanager=G4QIonIonCrossSection::GetPointer();
02280 if(PDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer();
02281 first=false;
02282 }
02283 #ifdef debug
02284 G4cout<<"G4QLowE::CXS: *DONE* p="<<p<<",Z="<<Z<<",N="<<N<<",PDG="<<PDG<<G4endl;
02285 #endif
02286 return CSmanager->GetCrossSection(true, p, Z, N, PDG);
02287 }
