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 "G4QIonIonElastic.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4HadronicDeprecate.hh"
00048
00049
00050
00051
00052
00053 G4int G4QIonIonElastic::nPartCWorld=85;
00054 std::vector<G4int> G4QIonIonElastic::ElementZ;
00055 std::vector<G4double> G4QIonIonElastic::ElProbInMat;
00056 std::vector<std::vector<G4int>*> G4QIonIonElastic::ElIsoN;
00057 std::vector<std::vector<G4double>*>G4QIonIonElastic::IsoProbInEl;
00058
00059
00060 G4QIonIonElastic::G4QIonIonElastic(const G4String& processName):
00061 G4VDiscreteProcess(processName, fHadronic)
00062 {
00063 G4HadronicDeprecate("G4QIonIonElastic");
00064
00065 #ifdef debug
00066 G4cout<<"G4QIonIonElastic::Constructor is called processName="<<processName<<G4endl;
00067 #endif
00068 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
00069
00070 }
00071
00072
00073 G4QIonIonElastic::~G4QIonIonElastic() {}
00074
00075
00076
00077
00078 G4double G4QIonIonElastic::GetMeanFreePath(const G4Track& aTrack, G4double,
00079 G4ForceCondition* Fc)
00080 {
00081 static const G4double mProt = G4QPDGCode(2212).GetMass()/MeV;
00082 *Fc = NotForced;
00083 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
00084 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
00085 if( !IsApplicable(*incidentParticleDefinition))
00086 G4cout<<"-Warning-G4QIonIonElastic::GetMeanFreePath: notImplementedParticle"<<G4endl;
00087
00088 G4double Momentum = incidentParticle->GetTotalMomentum();
00089 #ifdef debug
00090 G4double KinEn = incidentParticle->GetKineticEnergy();
00091 G4cout<<"G4QIonIonElastic::GetMeanFreePath: kinE="<<KinEn<<",Mom="<<Momentum<<G4endl;
00092 #endif
00093 const G4Material* material = aTrack.GetMaterial();
00094 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
00095 const G4ElementVector* theElementVector = material->GetElementVector();
00096 G4int nE=material->GetNumberOfElements();
00097 #ifdef debug
00098 G4cout<<"G4QIonIonElastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
00099 #endif
00100 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer();
00101 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
00102 G4int pPDG=0;
00103
00104 G4int pZ=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge());
00105 G4int pA=incidentParticleDefinition->GetBaryonNumber();
00106 if (incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020;
00107 else if (incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 1000020040;
00108 else if (incidentParticleDefinition == G4Triton::Triton() ) pPDG = 1000010030;
00109 else if (incidentParticleDefinition == G4He3::He3() ) pPDG = 1000020030;
00110 else if (incidentParticleDefinition == G4GenericIon::GenericIon() || (pZ > 0 && pA > 0))
00111 {
00112 pPDG=incidentParticleDefinition->GetPDGEncoding();
00113 #ifdef debug
00114 G4int prPDG=1000000000+10000*pZ+10*pA;
00115 G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl;
00116 #endif
00117 }
00118 else G4cout<<"-Warning-G4QIonIonElastic::GetMeanFreePath:Unknown projectile Ion"<<G4endl;
00119 Momentum/=incidentParticleDefinition->GetBaryonNumber();
00120 G4QIsotope* Isotopes = G4QIsotope::Get();
00121 G4double sigma=0.;
00122 G4int IPIE=IsoProbInEl.size();
00123 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)
00124 {
00125 std::vector<G4double>* SPI=IsoProbInEl[ip];
00126 SPI->clear();
00127 delete SPI;
00128 std::vector<G4int>* IsN=ElIsoN[ip];
00129 IsN->clear();
00130 delete IsN;
00131 }
00132 ElProbInMat.clear();
00133 ElementZ.clear();
00134 IsoProbInEl.clear();
00135 ElIsoN.clear();
00136 for(G4int i=0; i<nE; ++i)
00137 {
00138 G4Element* pElement=(*theElementVector)[i];
00139 G4int Z = static_cast<G4int>(pElement->GetZ());
00140 ElementZ.push_back(Z);
00141 G4int isoSize=0;
00142 G4int indEl=0;
00143 G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
00144 if(isoVector) isoSize=isoVector->size();
00145 #ifdef debug
00146 G4cout<<"G4QIonIonElastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl;
00147 #endif
00148 if(isoSize)
00149 {
00150 indEl=pElement->GetIndex()+1;
00151 #ifdef debug
00152 G4cout<<"G4QIIEl::GetMFP:iE="<<indEl<<", def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
00153 #endif
00154 if(!Isotopes->IsDefined(Z,indEl))
00155 {
00156 std::vector<std::pair<G4int,G4double>*>* newAbund =
00157 new std::vector<std::pair<G4int,G4double>*>;
00158 G4double* abuVector=pElement->GetRelativeAbundanceVector();
00159 for(G4int j=0; j<isoSize; j++)
00160 {
00161 G4int N=pElement->GetIsotope(j)->GetN()-Z;
00162 if(pElement->GetIsotope(j)->GetZ()!=Z) G4cerr<<"G4QIonIonEl::GetMeanFreePath Z="
00163 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00164 G4double abund=abuVector[j];
00165 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00166 #ifdef debug
00167 G4cout<<"G4QIonIonElastic::GetMeanFP:pair#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
00168 #endif
00169 newAbund->push_back(pr);
00170 }
00171 #ifdef debug
00172 G4cout<<"G4QIonIonElastic::GetMeanFP: pairVectorLength="<<newAbund->size()<<G4endl;
00173 #endif
00174 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund);
00175 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
00176 delete newAbund;
00177 }
00178 }
00179 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);
00180 std::vector<G4double>* SPI = new std::vector<G4double>;
00181 IsoProbInEl.push_back(SPI);
00182 std::vector<G4int>* IsN = new std::vector<G4int>;
00183 ElIsoN.push_back(IsN);
00184 G4int nIs=cs->size();
00185 #ifdef debug
00186 G4cout<<"G4QIonIonEl::GetMFP:=***=> #isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
00187 #endif
00188 G4double susi=0.;
00189 if(nIs) for(G4int j=0; j<nIs; j++)
00190 {
00191 std::pair<G4int,G4double>* curIs=(*cs)[j];
00192 G4int N=curIs->first;
00193 IsN->push_back(N);
00194 #ifdef debug
00195 G4cout<<"G4QIIE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl;
00196 #endif
00197 G4bool ccsf=false;
00198 #ifdef debug
00199 G4cout<<"G4QIonIonElastic::GMFP: GetCS #1 j="<<j<<G4endl;
00200 #endif
00201 G4double CSI=0.;
00202 if(Z==1 && !N)
00203 {
00204 G4double projM=G4QPDGCode(2212).GetNuclMass(pZ,pA-pZ,0);
00205 CSI=PCSmanager->GetCrossSection(true, Momentum*mProt/projM, Z, N, 2212);
00206 }
00207 else CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);
00208 #ifdef debug
00209 G4cout<<"G4QIonIonElastic::GMFP: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="<<Momentum
00210 <<", XSec="<<CSI/millibarn<<G4endl;
00211 #endif
00212 curIs->second = CSI;
00213 susi+=CSI;
00214 SPI->push_back(susi);
00215 }
00216 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];
00217 #ifdef debug
00218 G4cout<<"G4QIonIonEl::GMFP:<S>="<<Isotopes->GetMeanCrossSection(Z,indEl)<<",AddToSig="
00219 <<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
00220 #endif
00221 ElProbInMat.push_back(sigma);
00222 }
00223
00224 #ifdef debug
00225 G4cout<<"G4QIonIonElastic::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
00226 #endif
00227 if(sigma > 0.00000000001) return 1./sigma;
00228 return DBL_MAX;
00229 }
00230
00231
00232 G4bool G4QIonIonElastic::IsApplicable(const G4ParticleDefinition& particle)
00233 {
00234 G4int Z=static_cast<G4int>(particle.GetPDGCharge());
00235 G4int A=particle.GetBaryonNumber();
00236 if (particle == *( G4Deuteron::Deuteron() )) return true;
00237 else if (particle == *( G4Alpha::Alpha() )) return true;
00238 else if (particle == *( G4Triton::Triton() )) return true;
00239 else if (particle == *( G4He3::He3() )) return true;
00240 else if (particle == *( G4GenericIon::GenericIon() )) return true;
00241 else if (Z > 0 && A > 0) return true;
00242 #ifdef debug
00243 G4cout<<"***>>G4QIonIonElastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<", A="
00244 <<A<<", Z="<<Z<<G4endl;
00245 #endif
00246 return false;
00247 }
00248
00249 G4VParticleChange* G4QIonIonElastic::PostStepDoIt(const G4Track& track, const G4Step& step)
00250 {
00251 static const G4double mProt= G4QPDGCode(2212).GetMass();
00252 static const G4double fm2MeV2 = 3*38938./1.09;
00253 static G4bool CWinit = true;
00254 if(CWinit)
00255 {
00256 CWinit=false;
00257 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
00258 }
00259
00260 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
00261 const G4ParticleDefinition* particle=projHadron->GetDefinition();
00262 #ifdef debug
00263 G4cout<<"G4QIonIonElastic::PostStepDoIt: Before the GetMeanFreePath is called In4M="
00264 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
00265 <<particle->GetParticleType()<<", Subtp="<<particle->GetParticleSubType()<<G4endl;
00266 #endif
00267 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV;
00268 G4double momentum = projHadron->GetTotalMomentum()/MeV;
00269 G4double Momentum = proj4M.rho();
00270 if(std::fabs(Momentum-momentum)>.000001)
00271 G4cerr<<"-Warn-G4QIonIonElastic::PostStepDoIt:P(IU)="<<Momentum<<"#"<<momentum<<G4endl;
00272 G4double pM2=proj4M.m2();
00273 G4double pM=std::sqrt(pM2);
00274 #ifdef pdebug
00275 G4cout<<"G4QIonIonEl::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",pM="<<pM<<G4endl;
00276 #endif
00277 if (!IsApplicable(*particle))
00278 {
00279 G4cerr<<"G4QIonIonElastic::PostStepDoIt: Only NA elastic is implemented."<<G4endl;
00280 return 0;
00281 }
00282 const G4Material* material = track.GetMaterial();
00283 const G4ElementVector* theElementVector = material->GetElementVector();
00284 G4int nE=material->GetNumberOfElements();
00285 #ifdef debug
00286 G4cout<<"G4QIonIonElastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
00287 #endif
00288
00289 G4int projPDG=0;
00290 G4int pZ=static_cast<G4int>(particle->GetPDGCharge());
00291 G4int pA=particle->GetBaryonNumber();
00292 if (particle == G4Deuteron::Deuteron() ) projPDG= 1000010020;
00293 else if (particle == G4Alpha::Alpha() ) projPDG= 1000020040;
00294 else if (particle == G4Triton::Triton() ) projPDG= 1000010030;
00295 else if (particle == G4He3::He3() ) projPDG= 1000020030;
00296 else if (particle == G4GenericIon::GenericIon() || (pZ > 0 && pA > 0))
00297 {
00298 projPDG=particle->GetPDGEncoding();
00299 #ifdef debug
00300 G4int prPDG=1000000000+10000*pZ+10*pA;
00301 G4cout<<"G4QIonIonElastic::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl;
00302 #endif
00303 }
00304 else G4cout<<"-Warning-G4QIonIonElastic::PostStepDoIt:Unknown projectile Ion"<<G4endl;
00305 #ifdef debug
00306 G4int prPDG=particle->GetPDGEncoding();
00307 G4cout<<"G4QIonIonElastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
00308 #endif
00309 if(!projPDG)
00310 {
00311 G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:Undefined interactingNucleus"<<G4endl;
00312 return 0;
00313 }
00314 G4int pN=pA-pZ;
00315 G4int EPIM=ElProbInMat.size();
00316 #ifdef debug
00317 G4cout<<"G4QIonIonElastic::PSDI:m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
00318 #endif
00319 G4int i=0;
00320 if(EPIM>1)
00321 {
00322 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
00323 for(i=0; i<nE; ++i)
00324 {
00325 #ifdef debug
00326 G4cout<<"G4QIonIonElastic::PSDI: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl;
00327 #endif
00328 if (rnd<ElProbInMat[i]) break;
00329 }
00330 if(i>=nE) i=nE-1;
00331 }
00332 G4Element* pElement=(*theElementVector)[i];
00333 G4int tZ=static_cast<G4int>(pElement->GetZ());
00334 #ifdef debug
00335 G4cout<<"G4QIonIonElastic::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl;
00336 #endif
00337 if(tZ<=0)
00338 {
00339 G4cerr<<"---Warning---G4QIonIonElastic::PostStepDoIt: Element with Z="<<tZ<<G4endl;
00340 if(tZ<0) return 0;
00341 }
00342 std::vector<G4double>* SPI = IsoProbInEl[i];
00343 std::vector<G4int>* IsN = ElIsoN[i];
00344 G4int nofIsot=SPI->size();
00345 #ifdef debug
00346 G4cout<<"G4QIonIonElastic::PosStDoIt: nI="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
00347 #endif
00348 G4int j=0;
00349 if(nofIsot>1)
00350 {
00351 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand();
00352 for(j=0; j<nofIsot; ++j)
00353 {
00354 #ifdef debug
00355 G4cout<<"G4QIonIonElastic::PostStDI: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
00356 #endif
00357 if(rndI < (*SPI)[j]) break;
00358 }
00359 if(j>=nofIsot) j=nofIsot-1;
00360 }
00361 G4int tN =(*IsN)[j]; ;
00362 #ifdef debug
00363 G4cout<<"G4QIonIonElastic::PostStepDoIt:j="<<i<<",N(isotope)="<<tN<<",MeV="<<MeV<<G4endl;
00364 #endif
00365 if(tN<0)
00366 {
00367 G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:IsotopeZ="<<tZ<<" & 0>N="<<tN<<G4endl;
00368 return 0;
00369 }
00370 nOfNeutrons=tN;
00371 #ifdef debug
00372 G4cout<<"G4QIonIonElastic::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl;
00373 #endif
00374 aParticleChange.Initialize(track);
00375 #ifdef debug
00376 G4cout<<"G4QIonIonElastic::PostStepDoIt: track is initialized"<<G4endl;
00377 #endif
00378 G4double weight = track.GetWeight();
00379 G4double localtime = track.GetGlobalTime();
00380 G4ThreeVector position = track.GetPosition();
00381 #ifdef debug
00382 G4cout<<"G4QIonIonElastic::PostStepDoIt: before Touchable extraction"<<G4endl;
00383 #endif
00384 G4TouchableHandle trTouchable = track.GetTouchableHandle();
00385 #ifdef debug
00386 G4cout<<"G4QIonIonElastic::PostStepDoIt: Touchable is extracted"<<G4endl;
00387 #endif
00388
00389 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV;
00390 G4ParticleMomentum dir = projHadron->GetMomentumDirection();
00391
00392 G4bool revkin=false;
00393 G4ThreeVector bvel(0.,0.,0.);
00394 G4int tA=tZ+tN;
00395 G4int targPDG=90000000+tZ*1000+tN;
00396 G4QPDGCode targQPDG(targPDG);
00397 G4double tM=targQPDG.GetMass();
00398 G4LorentzVector targ4M(0.,0.,0.,tM);
00399 if(tZ==1 && !tN)
00400 {
00401 G4ForceCondition cond=NotForced;
00402 GetMeanFreePath(track, -27., &cond);
00403 #ifdef debug
00404 G4cout<<"G4QIonIonElastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
00405 #endif
00406 revkin=true;
00407 tZ=pZ;
00408 tN=pN;
00409 tA=tZ+tN;
00410 tM=pM;
00411 pZ=1;
00412 pN=0;
00413 pA=1;
00414 pM=mProt;
00415 bvel=proj4M.vect()/proj4M.e();
00416 proj4M=targ4M.boost(-bvel);
00417 targ4M=G4LorentzVector(0.,0.,0.,tM);
00418 Momentum = proj4M.rho();
00419 }
00420 G4LorentzVector tot4M=proj4M+targ4M;
00421 #ifdef debug
00422 G4cout<<"G4QIonIonElastic::PostStDI: tM="<<tM<<", p4M="<<proj4M<<", t4M="<<tot4M<<G4endl;
00423 #endif
00424 EnMomConservation=tot4M;
00425 G4VQCrossSection* PELmanager=G4QProtonElasticCrossSection::GetPointer();
00426 G4VQCrossSection* NELmanager=G4QNeutronElasticCrossSection::GetPointer();
00427 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer();
00428
00429 #ifdef debug
00430 G4cout<<"G4QIIEl::PSDI:f, P="<<Momentum<<",Z="<<tZ<<",N="<<tN<<",tPDG="<<projPDG<<G4endl;
00431 #endif
00432
00433 G4double xSec=0.;
00434 if(revkin) xSec=PELmanager->GetCrossSection(false, Momentum, tZ, tN, 2212);
00435 else xSec=CSmanager->GetCrossSection(false, Momentum, tZ, tN, projPDG);
00436 #ifdef debug
00437 G4cout<<"G4QIIEl::PSDI: pPDG="<<projPDG<<",P="<<Momentum<<",CS="<<xSec/millibarn<<G4endl;
00438 #endif
00439 #ifdef nandebug
00440 if(xSec>0. || xSec<0. || xSec==0);
00441 else G4cout<<"-NaN-Warning-G4QIonIonElastic::PostStDoIt: xSec="<<xSec/millibarn<<G4endl;
00442 #endif
00443
00444 if(xSec <= 0.)
00445 {
00446 #ifdef pdebug
00447 G4cerr<<"-Warning-G4QIonIonElastic::PostStDoIt: *Zero cross-section* PDG="<<projPDG
00448 <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl;
00449 #endif
00450
00451 aParticleChange.ProposeEnergy(kinEnergy);
00452 aParticleChange.ProposeLocalEnergyDeposit(0.);
00453 aParticleChange.ProposeMomentumDirection(dir) ;
00454 return G4VDiscreteProcess::PostStepDoIt(track,step);
00455 }
00456 G4double mint=0;
00457 G4double maxt=0;
00458 G4double dtM=tM+tM;
00459 if(revkin)
00460 {
00461 mint=PELmanager->GetExchangeT(tZ,tN,2212);
00462 maxt=PELmanager->GetHMaxT();
00463 }
00464 else
00465 {
00466 G4double PA=Momentum*pA;
00467 G4double PA2=PA*PA;
00468 maxt=dtM*PA2/(std::sqrt(PA2+pM2)+tM/2+pM2/dtM);
00469 #ifdef pdebug
00470 G4cout<<"G4QIonIonElastic::PostStDoIt:pPDG="<<projPDG<<",tPDG="<<targPDG<<",P="
00471 <<Momentum<<",CS="<<xSec<<",maxt="<<maxt<<G4endl;
00472 #endif
00473 xSec=PELmanager->GetCrossSection(false, Momentum, 1, 0, 2212);
00474 G4double B1=PELmanager->GetSlope(1,0,2212);
00475 xSec=NELmanager->GetCrossSection(false, Momentum, 1, 0, 2112);
00476 G4double B2 =NELmanager->GetSlope(1,0,2112);
00477 G4double mB =((pZ*tZ+pN*tN)*B1+(pZ*tN+pN*tZ)*B2)/(pA+tA);
00478 G4double pR2=std::pow(pA+4.,.305)/fm2MeV2;
00479 G4double tR2=std::pow(tA+4.,.305)/fm2MeV2;
00480 G4double eB =mB+pR2+tR2;
00481 mint=-std::log(1.-G4UniformRand()*(1.-std::exp(-eB*maxt)))/eB;
00482 mint+=mint;
00483 #ifdef pdebug
00484 G4cout<<"G4QIonIonElastic::PostStDoIt:B1="<<B1<<",B2="<<B2<<",mB="<<mB
00485 <<",pR2="<<pR2<<",tR2="<<tR2<<",eB="<<eB<<",mint="<<mint<<G4endl;
00486 #endif
00487 }
00488 #ifdef nandebug
00489 if(mint>-.0000001);
00490 else G4cout<<"-Warning-G4QIonIonElastic::PostStDoIt:-t="<<mint<<G4endl;
00491 #endif
00492 G4double cost=1.-mint/maxt;
00493
00494 #ifdef ppdebug
00495 G4cout<<"G4QIonIonElastic::PoStDoI:t="<<mint<<",dpcm2="<<CSmanager->GetHMaxT()<<",Ek="
00496 <<kinEnergy<<",tM="<<tM<<",pM="<<pM<<",cost="<<cost<<G4endl;
00497 #endif
00498 if(cost>1. || cost<-1. || !(cost>-1. || cost<1.))
00499 {
00500 if(cost>1.000001 || cost<-1.000001 || !(cost>-1. || cost<1.))
00501 {
00502 G4double tM2=tM*tM;
00503 G4double pEn=pM+kinEnergy;
00504 G4double sM=dtM*pEn+tM2+pM2;
00505 G4double twop2cm=(tM2+tM2)*(pEn*pEn-pM2)/sM;
00506 G4cout<<"-Warning-G4QIonIonElastic::PoStDI:cos="<<cost<<",t="<<mint<<",T="<<kinEnergy
00507 <<",tM="<<tM<<",tmax="<<2*kinEnergy*tM<<",p="<<projPDG<<",t="<<targPDG<<G4endl;
00508 G4cout<<"G4QIonIonElastic::PSDI:dpcm2="<<twop2cm<<"="<<CSmanager->GetHMaxT()<<G4endl;
00509 }
00510 if (cost>1.) cost=1.;
00511 else if(cost<-1.) cost=-1.;
00512 }
00513 G4LorentzVector scat4M(0.,0.,0.,pM);
00514 G4LorentzVector reco4M(0.,0.,0.,tM);
00515 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-tM-pM)*.01);
00516 if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
00517 {
00518 G4cout<<"-Warning-G4QIonIonE::PSDI:t4M="<<tot4M<<",pM="<<pM<<",tM="<<tM<<",cost="
00519 <<cost<<G4endl;
00520 }
00521 if(revkin)
00522 {
00523 G4LorentzVector tmpLV=scat4M.boost(bvel);
00524 scat4M=reco4M.boost(bvel);
00525 reco4M=tmpLV;
00526 pM=tM;
00527 tM=mProt;
00528 }
00529 #ifdef debug
00530 G4cout<<"G4QIonIonElast::PSDI:s4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
00531 G4cout<<"G4QIonIonElastic::PSDI:scatE="<<scat4M.e()-pM<<",recoE="<<reco4M.e()-tM<<",d4M="
00532 <<tot4M-scat4M-reco4M<<G4endl;
00533 #endif
00534
00535 G4double finE=scat4M.e()-pM;
00536 if(finE>0.0) aParticleChange.ProposeEnergy(finE);
00537 else
00538 {
00539 if(finE<-1.e-8 || !(finE>-1.||finE<1.))
00540 G4cerr<<"*Warning*G4QIonIonElastic::PostStDoIt: Zero or negative scattered E="<<finE
00541 <<", s4M="<<scat4M<<", r4M="<<reco4M<<", d4M="<<tot4M-scat4M-reco4M<<G4endl;
00542 aParticleChange.ProposeEnergy(0.) ;
00543 aParticleChange.ProposeTrackStatus(fStopButAlive);
00544 }
00545 G4ThreeVector findir=scat4M.vect()/scat4M.rho();
00546 aParticleChange.ProposeMomentumDirection(findir);
00547 EnMomConservation-=scat4M;
00548
00549 G4DynamicParticle* theSec = new G4DynamicParticle;
00550 #ifdef pdebug
00551 G4cout<<"G4QIonIonElastic::PostStepDoIt: Ion tZ="<<tZ<<", tA="<<tA<<G4endl;
00552 #endif
00553 G4ParticleDefinition* theDefinition=0;
00554 if(revkin) theDefinition = G4Proton::Proton();
00555 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(tZ,tA,0,tZ);
00556 if(!theDefinition)G4cout<<"-Warning-G4QIonIonElastic::PoStDI:drop PDG="<<targPDG<<G4endl;
00557 #ifdef pdebug
00558 G4cout<<"G4QIonIonElastic::PoStDI:RecoilName="<<theDefinition->GetParticleName()<<G4endl;
00559 #endif
00560 theSec->SetDefinition(theDefinition);
00561 EnMomConservation-=reco4M;
00562 #ifdef tdebug
00563 G4cout<<"G4QIonIonElastic::PSD:"<<targPDG<<reco4M<<reco4M.m()<<EnMomConservation<<G4endl;
00564 #endif
00565 theSec->Set4Momentum(reco4M);
00566 #ifdef debug
00567 G4ThreeVector curD=theSec->GetMomentumDirection();
00568 G4double curM=theSec->GetMass();
00569 G4double curE=theSec->GetKineticEnergy()+curM;
00570 G4cout<<"G4QIonIonElastic::PSDI: p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;
00571 #endif
00572
00573 G4Track* aNewTrack = new G4Track(theSec, localtime, position );
00574 aNewTrack->SetWeight(weight);
00575 aNewTrack->SetTouchableHandle(trTouchable);
00576 aParticleChange.AddSecondary( aNewTrack );
00577 #ifdef debug
00578 G4cout<<"G4QIonIonElastic::PostStepDoIt: **** PostStepDoIt is done ****"<<G4endl;
00579 #endif
00580 return G4VDiscreteProcess::PostStepDoIt(track, step);
00581 }