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 #include "G4QElastic.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4HadronicDeprecate.hh"
00044
00045
00046
00047
00048
00049 G4int G4QElastic::nPartCWorld=85;
00050 std::vector<G4int> G4QElastic::ElementZ;
00051 std::vector<G4double> G4QElastic::ElProbInMat;
00052 std::vector<std::vector<G4int>*> G4QElastic::ElIsoN;
00053 std::vector<std::vector<G4double>*>G4QElastic::IsoProbInEl;
00054
00055
00056 G4QElastic::G4QElastic(const G4String& processName):
00057 G4VDiscreteProcess(processName, fHadronic)
00058 {
00059 G4HadronicDeprecate("G4QElastic");
00060
00061 #ifdef debug
00062 G4cout<<"G4QElastic::Constructor is called processName="<<processName<<G4endl;
00063 #endif
00064 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
00065
00066 }
00067
00068
00069 G4QElastic::~G4QElastic() {}
00070
00071
00072 G4LorentzVector G4QElastic::GetEnegryMomentumConservation() {return EnMomConservation;}
00073
00074 G4int G4QElastic::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
00075
00076
00077
00078
00079 G4double G4QElastic::GetMeanFreePath(const G4Track& aTrack,G4double Q,G4ForceCondition* Fc)
00080 {
00081 *Fc = NotForced;
00082 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
00083 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
00084 if( !IsApplicable(*incidentParticleDefinition))
00085 G4cout<<"*W*G4QElastic::GetMeanFreePath: is called for notImplementedParticle"<<G4endl;
00086
00087 G4double Momentum = incidentParticle->GetTotalMomentum();
00088 #ifdef debug
00089 G4double KinEn = incidentParticle->GetKineticEnergy();
00090 G4cout<<"G4QElastic::GetMeanFreePath: kinE="<<KinEn<<",Mom="<<Momentum<<G4endl;
00091 #endif
00092 const G4Material* material = aTrack.GetMaterial();
00093 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
00094 const G4ElementVector* theElementVector = material->GetElementVector();
00095 G4int nE=material->GetNumberOfElements();
00096 G4int pPDG=0;
00097 if (incidentParticleDefinition == G4Proton::Proton() ) pPDG=2212;
00098 else if(incidentParticleDefinition == G4Neutron::Neutron() ) pPDG=2112;
00099 else if(incidentParticleDefinition == G4PionPlus::PionPlus() ) pPDG= 211;
00100 else if(incidentParticleDefinition == G4PionMinus::PionMinus() ) pPDG=-211;
00101 else if(incidentParticleDefinition == G4KaonPlus::KaonPlus() ) pPDG= 321;
00102 else if(incidentParticleDefinition == G4KaonMinus::KaonMinus() ) pPDG=-321;
00103 else if(incidentParticleDefinition == G4KaonZeroLong::KaonZeroLong() ) pPDG= 130;
00104 else if(incidentParticleDefinition == G4KaonZeroShort::KaonZeroShort() ) pPDG= 310;
00105 else if(incidentParticleDefinition == G4Lambda::Lambda() ) pPDG= 3122;
00106 else if(incidentParticleDefinition == G4SigmaPlus::SigmaPlus() ) pPDG= 3222;
00107 else if(incidentParticleDefinition == G4SigmaMinus::SigmaMinus() ) pPDG= 3112;
00108 else if(incidentParticleDefinition == G4SigmaZero::SigmaZero() ) pPDG= 3212;
00109 else if(incidentParticleDefinition == G4XiMinus::XiMinus() ) pPDG= 3312;
00110 else if(incidentParticleDefinition == G4XiZero::XiZero() ) pPDG= 3322;
00111 else if(incidentParticleDefinition == G4OmegaMinus::OmegaMinus() ) pPDG= 3334;
00112 else if(incidentParticleDefinition == G4AntiNeutron::AntiNeutron() ) pPDG=-2112;
00113 else if(incidentParticleDefinition == G4AntiProton::AntiProton() ) pPDG=-2212;
00114 else if(incidentParticleDefinition == G4AntiLambda::AntiLambda() ) pPDG=-3122;
00115 else if(incidentParticleDefinition == G4AntiSigmaPlus::AntiSigmaPlus() ) pPDG=-3222;
00116 else if(incidentParticleDefinition == G4AntiSigmaMinus::AntiSigmaMinus() ) pPDG=-3112;
00117 else if(incidentParticleDefinition == G4AntiSigmaZero::AntiSigmaZero() ) pPDG=-3212;
00118 else if(incidentParticleDefinition == G4AntiXiMinus::AntiXiMinus() ) pPDG=-3312;
00119 else if(incidentParticleDefinition == G4AntiXiZero::AntiXiZero() ) pPDG=-3322;
00120 else if(incidentParticleDefinition == G4AntiOmegaMinus::AntiOmegaMinus() ) pPDG=-3334;
00121 else G4cout<<"Warning-G4QElastic::GetMFP:"<<incidentParticleDefinition->GetParticleName()
00122 <<" isn't implemented (only SU(3) particles)"<<G4endl;
00123 #ifdef pdebug
00124 G4cout<<"G4QElastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial,proj="<<pPDG<<G4endl;
00125 #endif
00126 G4VQCrossSection* CSmanager=0;
00127 G4VQCrossSection* CSmanager2=0;
00128 if (pPDG==2212) CSmanager=G4QProtonElasticCrossSection::GetPointer();
00129 else if(pPDG==2112) CSmanager=G4QNeutronElasticCrossSection::GetPointer();
00130 else if(pPDG== 211) CSmanager=G4QPionPlusElasticCrossSection::GetPointer();
00131 else if(pPDG==-211) CSmanager=G4QPionMinusElasticCrossSection::GetPointer();
00132 else if(pPDG== 321) CSmanager=G4QKaonPlusElasticCrossSection::GetPointer();
00133 else if(pPDG==-321) CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00134 else if(pPDG== 130 || pPDG== 310)
00135 {
00136 CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00137 CSmanager2=G4QKaonPlusElasticCrossSection::GetPointer();
00138 }
00139 else if(pPDG==3222) CSmanager=G4QHyperonPlusElasticCrossSection::GetPointer();
00140 else if(pPDG>3110 && pPDG<3335) CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00141 else if(pPDG>-3335&&pPDG<-1110) CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00142 else G4cout<<"*Warning*G4QElastic::GetMeanFreePath: wrong PDG="<<pPDG<<G4endl;
00143 G4QIsotope* Isotopes = G4QIsotope::Get();
00144 G4double sigma=0.;
00145 G4int IPIE=IsoProbInEl.size();
00146 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)
00147 {
00148 std::vector<G4double>* SPI=IsoProbInEl[ip];
00149 SPI->clear();
00150 delete SPI;
00151 std::vector<G4int>* IsN=ElIsoN[ip];
00152 IsN->clear();
00153 delete IsN;
00154 }
00155 ElProbInMat.clear();
00156 ElementZ.clear();
00157 IsoProbInEl.clear();
00158 ElIsoN.clear();
00159 for(G4int i=0; i<nE; ++i)
00160 {
00161 G4Element* pElement=(*theElementVector)[i];
00162 G4int Z = static_cast<G4int>(pElement->GetZ());
00163 ElementZ.push_back(Z);
00164 G4int isoSize=0;
00165 G4int indEl=0;
00166 G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
00167 if(isoVector) isoSize=isoVector->size();
00168 #ifdef debug
00169 G4cout<<"G4QElastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl;
00170 #endif
00171 if(isoSize)
00172 {
00173 indEl=pElement->GetIndex()+1;
00174 #ifdef debug
00175 G4cout<<"G4QEl::GetMFP: iE="<<indEl<<", def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
00176 #endif
00177 if(!Isotopes->IsDefined(Z,indEl))
00178 {
00179 std::vector<std::pair<G4int,G4double>*>* newAbund =
00180 new std::vector<std::pair<G4int,G4double>*>;
00181 G4double* abuVector=pElement->GetRelativeAbundanceVector();
00182 for(G4int j=0; j<isoSize; j++)
00183 {
00184 G4int N=pElement->GetIsotope(j)->GetN()-Z;
00185 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QElastic::GetMeanFreePath: Z="
00186 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00187 G4double abund=abuVector[j];
00188 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00189 #ifdef debug
00190 G4cout<<"G4QElastic::GetMeanFreePath:pair#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
00191 #endif
00192 newAbund->push_back(pr);
00193 }
00194 #ifdef debug
00195 G4cout<<"G4QElastic::GetMeanFreePath: pairVectorLength="<<newAbund->size()<<G4endl;
00196 #endif
00197 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund);
00198 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
00199 delete newAbund;
00200 }
00201 }
00202 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);
00203 std::vector<G4double>* SPI = new std::vector<G4double>;
00204 IsoProbInEl.push_back(SPI);
00205 std::vector<G4int>* IsN = new std::vector<G4int>;
00206 ElIsoN.push_back(IsN);
00207 G4int nIs=cs->size();
00208 #ifdef debug
00209 G4cout<<"G4QEl::GMFP:=***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
00210 #endif
00211 G4double susi=0.;
00212 if(nIs) for(G4int j=0; j<nIs; j++)
00213 {
00214 std::pair<G4int,G4double>* curIs=(*cs)[j];
00215 G4int N=curIs->first;
00216 IsN->push_back(N);
00217 #ifdef debug
00218 G4cout<<"G4QEl::GMFP:*true*,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
00219 #endif
00220 G4bool ccsf=true;
00221 if(Q==-27.) ccsf=false;
00222 #ifdef debug
00223 G4cout<<"G4QEl::GMFP: GetCS #1 j="<<j<<G4endl;
00224 #endif
00225 G4double CSI=0.;
00226 if(!CSmanager2) CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);
00227 else CSI=(CSmanager->GetCrossSection(ccsf,Momentum,Z,N,-321)+
00228 CSmanager2->GetCrossSection(ccsf,Momentum,Z,N,321) )/2.;
00229 #ifdef debug
00230 G4cout<<"G4QEl::GMFP: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="<<Momentum<<", XSec="
00231 <<CSI/millibarn<<G4endl;
00232 #endif
00233 curIs->second = CSI;
00234 susi+=CSI;
00235 SPI->push_back(susi);
00236 }
00237 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];
00238 #ifdef debug
00239 G4cout<<"G4QEl::GMFP: <S>="<<Isotopes->GetMeanCrossSection(Z,indEl)<<", AddToSigma="
00240 <<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
00241 #endif
00242 ElProbInMat.push_back(sigma);
00243 }
00244
00245 #ifdef debug
00246 G4cout<<"G4QElastic::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
00247 #endif
00248 if(sigma > 0.) return 1./sigma;
00249 return DBL_MAX;
00250 }
00251
00252
00253 G4bool G4QElastic::IsApplicable(const G4ParticleDefinition& particle)
00254 {
00255 if (particle == *( G4Proton::Proton() )) return true;
00256 else if (particle == *( G4Neutron::Neutron() )) return true;
00257 else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
00258 else if (particle == *( G4TauPlus::TauPlus() )) return true;
00259 else if (particle == *( G4TauMinus::TauMinus() )) return true;
00260 else if (particle == *( G4Electron::Electron() )) return true;
00261 else if (particle == *( G4Positron::Positron() )) return true;
00262 else if (particle == *( G4Gamma::Gamma() )) return true;
00263 else if (particle == *( G4MuonPlus::MuonPlus() )) return true;
00264 else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
00265 else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
00266 else if (particle == *( G4PionMinus::PionMinus() )) return true;
00267 else if (particle == *( G4PionPlus::PionPlus() )) return true;
00268 else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
00269 else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
00270 else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
00271 else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
00272 else if (particle == *( G4Lambda::Lambda() )) return true;
00273 else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
00274 else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
00275 else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
00276 else if (particle == *( G4XiMinus::XiMinus() )) return true;
00277 else if (particle == *( G4XiZero::XiZero() )) return true;
00278 else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
00279 else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
00280 else if (particle == *( G4AntiProton::AntiProton() )) return true;
00281 #ifdef debug
00282 G4cout<<"***>>G4QElastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
00283 #endif
00284 return false;
00285 }
00286
00287 G4VParticleChange* G4QElastic::PostStepDoIt(const G4Track& track, const G4Step& step)
00288 {
00289
00290
00291
00292
00293
00294 static G4bool CWinit = true;
00295 if(CWinit)
00296 {
00297 CWinit=false;
00298 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
00299 }
00300
00301 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
00302 const G4ParticleDefinition* particle=projHadron->GetDefinition();
00303 #ifdef debug
00304 G4cout<<"G4QElastic::PostStepDoIt: Before the GetMeanFreePath is called In4M="
00305 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
00306 <<particle->GetParticleType()<<", Subtp="<<particle->GetParticleSubType()<<G4endl;
00307 #endif
00308 G4ForceCondition cond=NotForced;
00309 GetMeanFreePath(track, -27., &cond);
00310 #ifdef debug
00311 G4cout<<"G4QElastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
00312 #endif
00313 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV;
00314 G4LorentzVector scat4M=proj4M;
00315 G4double momentum = projHadron->GetTotalMomentum()/MeV;
00316 G4double Momentum = proj4M.rho();
00317 if(std::fabs(Momentum-momentum)>.000001)
00318 G4cerr<<"*War*G4QElastic::PostStepDoIt:P(IU)="<<Momentum<<"#"<<momentum<<G4endl;
00319 G4double pM2=proj4M.m2();
00320 G4double pM=std::sqrt(pM2);
00321 #ifdef debug
00322 G4cout<<"G4QElastic::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",pM="<<pM
00323 <<",scat4M="<<scat4M<<scat4M.m()<<G4endl;
00324 #endif
00325 if (!IsApplicable(*particle))
00326 {
00327 G4cerr<<"G4QElastic::PostStepDoIt: Only NA elastic is implemented."<<G4endl;
00328 return 0;
00329 }
00330 const G4Material* material = track.GetMaterial();
00331 const G4ElementVector* theElementVector = material->GetElementVector();
00332 G4int nE=material->GetNumberOfElements();
00333 #ifdef debug
00334 G4cout<<"G4QElastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
00335 #endif
00336 G4int projPDG=0;
00337
00338 if (particle == G4Proton::Proton() ) projPDG= 2212;
00339 else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
00340 else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
00341 else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
00342 else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 2112;
00343 else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
00344 else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
00345 else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310;
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359 else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
00360 else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
00361 else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
00362 else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
00363 else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
00364 else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
00365 else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
00366 else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
00367 else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
00368 else if (particle == G4AntiLambda::AntiLambda() ) projPDG=-3122;
00369 else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) projPDG=-3222;
00370 else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112;
00371 else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) projPDG=-3212;
00372 else if (particle == G4AntiXiMinus::AntiXiMinus() ) projPDG=-3312;
00373 else if (particle == G4AntiXiZero::AntiXiZero() ) projPDG=-3322;
00374 else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334;
00375 #ifdef pdebug
00376 G4int prPDG=particle->GetPDGEncoding();
00377 G4cout<<"G4QElastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
00378 #endif
00379 if(!projPDG)
00380 {
00381 G4cerr<<"*Warning*G4QElastic::PostStepDoIt: Undefined interacting hadron"<<G4endl;
00382 return 0;
00383 }
00384 G4int EPIM=ElProbInMat.size();
00385 #ifdef debug
00386 G4cout<<"G4QElastic::PostStDoIt:m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
00387 #endif
00388 G4int i=0;
00389 if(EPIM>1)
00390 {
00391 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
00392 for(i=0; i<nE; ++i)
00393 {
00394 #ifdef debug
00395 G4cout<<"G4QElastic::PostStepDoIt:EPM["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
00396 #endif
00397 if (rnd<ElProbInMat[i]) break;
00398 }
00399 if(i>=nE) i=nE-1;
00400 }
00401 G4Element* pElement=(*theElementVector)[i];
00402 G4int Z=static_cast<G4int>(pElement->GetZ());
00403 #ifdef debug
00404 G4cout<<"G4QElastic::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
00405 #endif
00406 if(Z<=0)
00407 {
00408 G4cerr<<"---Warning---G4QElastic::PostStepDoIt: Element with Z="<<Z<<G4endl;
00409 if(Z<0) return 0;
00410 }
00411 std::vector<G4double>* SPI = IsoProbInEl[i];
00412 std::vector<G4int>* IsN = ElIsoN[i];
00413 G4int nofIsot=SPI->size();
00414 #ifdef debug
00415 G4cout<<"G4QElastic::PosStDoIt: nI="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
00416 #endif
00417 G4int j=0;
00418 if(nofIsot>1)
00419 {
00420 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand();
00421 for(j=0; j<nofIsot; ++j)
00422 {
00423 #ifdef debug
00424 G4cout<<"G4QElastic::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
00425 #endif
00426 if(rndI < (*SPI)[j]) break;
00427 }
00428 if(j>=nofIsot) j=nofIsot-1;
00429 }
00430 G4int N =(*IsN)[j]; ;
00431 #ifdef debug
00432 G4cout<<"G4QElastic::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<", MeV="<<MeV<<G4endl;
00433 #endif
00434 if(N<0)
00435 {
00436 G4cerr<<"*Warning*G4QElastic::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
00437 return 0;
00438 }
00439 nOfNeutrons=N;
00440 #ifdef debug
00441 G4cout<<"G4QElastic::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
00442 #endif
00443 aParticleChange.Initialize(track);
00444 #ifdef debug
00445 G4cout<<"G4QElastic::PostStepDoIt: track is initialized"<<G4endl;
00446 #endif
00447 G4double weight = track.GetWeight();
00448 G4double localtime = track.GetGlobalTime();
00449 G4ThreeVector position = track.GetPosition();
00450 #ifdef debug
00451 G4cout<<"G4QElastic::PostStepDoIt: before Touchable extraction"<<G4endl;
00452 #endif
00453 G4TouchableHandle trTouchable = track.GetTouchableHandle();
00454 #ifdef debug
00455 G4cout<<"G4QElastic::PostStepDoIt: Touchable is extracted"<<G4endl;
00456 #endif
00457
00458 G4int targPDG=90000000+Z*1000+N;
00459 G4QPDGCode targQPDG(targPDG);
00460 G4double tM=targQPDG.GetMass();
00461 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV;
00462 G4ParticleMomentum dir = projHadron->GetMomentumDirection();
00463 G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM);
00464 #ifdef debug
00465 G4cout<<"G4QElastic::PostStepDoIt: tM="<<tM<<", p4M="<<proj4M<<", t4M="<<tot4M<<G4endl;
00466 #endif
00467 EnMomConservation=tot4M;
00468 G4VQCrossSection* CSmanager=0;
00469 G4VQCrossSection* CSmanager2=0;
00470 if (projPDG==2212) CSmanager=G4QProtonElasticCrossSection::GetPointer();
00471 else if(projPDG==2112) CSmanager=G4QNeutronElasticCrossSection::GetPointer();
00472 else if(projPDG== 211) CSmanager=G4QPionPlusElasticCrossSection::GetPointer();
00473 else if(projPDG==-211) CSmanager=G4QPionMinusElasticCrossSection::GetPointer();
00474 else if(projPDG== 321) CSmanager=G4QKaonPlusElasticCrossSection::GetPointer();
00475 else if(projPDG==-321) CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00476 else if(projPDG== 130 || projPDG== 310)
00477 {
00478 CSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00479 CSmanager2=G4QKaonPlusElasticCrossSection::GetPointer();
00480 }
00481 else if(projPDG==3222) CSmanager=G4QHyperonPlusElasticCrossSection::GetPointer();
00482 else if(projPDG>3110&&projPDG<3335)CSmanager=G4QHyperonElasticCrossSection::GetPointer();
00483 else if(projPDG>-3335 && projPDG<-1110)
00484 CSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00485 else G4cout<<"*Warning*G4QElastic::PostStepDoIt: wrong PDG="<<projPDG<<G4endl;
00486 #ifdef debug
00487 G4cout<<"G4QElas::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl;
00488 #endif
00489 G4double xSec=0.;
00490 if(!CSmanager2) xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG);
00491 else xSec=(CSmanager->GetCrossSection(false,Momentum,Z,N,-321)+
00492 CSmanager2->GetCrossSection(false,Momentum,Z,N,321) )/2.;
00493 #ifdef debug
00494 G4cout<<"G4QElast::PSDI:pPDG="<<projPDG<<",P="<<Momentum<<",CS="<<xSec/millibarn<<G4endl;
00495 #endif
00496 #ifdef nandebug
00497 if(xSec>0. || xSec<0. || xSec==0);
00498 else G4cout<<"*Warning*G4QElastic::PostStDoIt: xSec="<<xSec/millibarn<<G4endl;
00499 #endif
00500
00501 if(xSec <= 0.)
00502 {
00503 #ifdef debug
00504 G4cerr<<"*Warning*G4QElastic::PostStDoIt: *Zero cross-section* PDG="<<projPDG<<",tPDG="
00505 <<targPDG<<",P="<<Momentum<<G4endl;
00506 #endif
00507
00508 aParticleChange.ProposeEnergy(kinEnergy);
00509 aParticleChange.ProposeLocalEnergyDeposit(0.);
00510 aParticleChange.ProposeMomentumDirection(dir) ;
00511 return G4VDiscreteProcess::PostStepDoIt(track,step);
00512 }
00513 G4double mint=CSmanager->GetExchangeT(Z, N, projPDG);
00514 #ifdef debug
00515 G4cout<<"G4QElast::PSDI:pPDG="<<projPDG<<",tPDG="<<targPDG<<",P="<<Momentum<<",CS="
00516 <<xSec<<",-t="<<mint<<G4endl;
00517 #endif
00518 #ifdef nandebug
00519 if(mint>-.0000001);
00520 else G4cout<<"*Warning*G4QElastic::PostStDoIt:-t="<<mint<<G4endl;
00521 #endif
00522 G4double maxt=CSmanager->GetHMaxT();
00523 G4double cost=1.-mint/maxt;
00524
00525 #ifdef ppdebug
00526 G4cout<<"G4QElastic::PoStDoI:t="<<mint<<", maxt="<<maxt<<",Ek="<<kinEnergy<<",tM="<<tM
00527 <<",pM="<<pM<<",cost="<<cost<<G4endl;
00528 #endif
00529 if(cost>1. || cost<-1. || !(cost>-1. || cost<1.))
00530 {
00531 if(cost>1.000001 || cost<-1.000001 || !(cost>-1. || cost<1.))
00532 {
00533 G4double tM2=tM*tM;
00534 G4double pEn=pM+kinEnergy;
00535 G4double sM=(tM+tM)*pEn+tM2+pM2;
00536 G4double twop2cm=(tM2+tM2)*(pEn*pEn-pM2)/sM;
00537 G4cout<<"*Warning*G4QElastic::PostStepDoIt:cos="<<cost<<",t="<<mint<<",T="<<kinEnergy
00538 <<",tM="<<tM<<",tmax="<<2*kinEnergy*tM<<",p="<<projPDG<<",t="<<targPDG<<G4endl;
00539 G4cout<<"..G4QElastic::PoStDoI: dpcm2="<<twop2cm<<"="<<maxt<<G4endl;
00540 }
00541 if (cost>1.) cost=1.;
00542 else if(cost<-1.) cost=-1.;
00543 }
00544 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tM);
00545 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-tM-pM)*.01);
00546 if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
00547 {
00548 G4cerr<<"G4QElastic::PSD:t4M="<<tot4M<<",pM="<<pM<<",tM="<<tM<<",cost="<<cost<<G4endl;
00549
00550 }
00551 #ifdef debug
00552 G4cout<<"G4QElastic::PoStDoIt:s4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
00553 G4cout<<"G4QElastic::PoStDoIt: scatE="<<scat4M.e()-pM<<", recoE="<<reco4M.e()-tM<<",d4M="
00554 <<tot4M-scat4M-reco4M<<G4endl;
00555 #endif
00556
00557 G4double finE=scat4M.e()-pM;
00558 if(finE>0.0) aParticleChange.ProposeEnergy(finE);
00559 else
00560 {
00561 if(finE<-1.e-8 || !(finE>-1.||finE<1.))
00562 G4cerr<<"*Warning*G4QElastic::PostStDoIt: Zero or negative scattered E="<<finE
00563 <<", s4M="<<scat4M<<", r4M="<<reco4M<<", d4M="<<tot4M-scat4M-reco4M<<G4endl;
00564
00565 aParticleChange.ProposeEnergy(0.) ;
00566 aParticleChange.ProposeTrackStatus(fStopButAlive);
00567 }
00568 G4ThreeVector findir=scat4M.vect()/scat4M.rho();
00569 aParticleChange.ProposeMomentumDirection(findir);
00570 EnMomConservation-=scat4M;
00571
00572 G4DynamicParticle* theSec = new G4DynamicParticle;
00573 G4int aA = Z+N;
00574 #ifdef debug
00575 G4cout<<"G4QElastic::PostStepDoIt: Ion Z="<<Z<<", A="<<aA<<G4endl;
00576 #endif
00577 G4ParticleDefinition* theDefinition=G4ParticleTable::GetParticleTable()
00578 ->FindIon(Z,aA,0,Z);
00579 if(!theDefinition)G4cout<<"*Warning*G4QElastic::PostStepDoIt:drop PDG="<<targPDG<<G4endl;
00580 #ifdef debug
00581 G4cout<<"G4QElastic::PostStepDoIt:RecoilName="<<theDefinition->GetParticleName()<<G4endl;
00582 #endif
00583 theSec->SetDefinition(theDefinition);
00584
00585 EnMomConservation-=reco4M;
00586 #ifdef tdebug
00587 G4cout<<"G4QElastic::PostSDoIt:"<<targPDG<<reco4M<<reco4M.m()<<EnMomConservation<<G4endl;
00588 #endif
00589 theSec->Set4Momentum(reco4M);
00590 #ifdef debug
00591 G4ThreeVector curD=theSec->GetMomentumDirection();
00592 G4double curM=theSec->GetMass();
00593 G4double curE=theSec->GetKineticEnergy()+curM;
00594 G4cout<<"G4QElastic::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;
00595 #endif
00596
00597 G4Track* aNewTrack = new G4Track(theSec, localtime, position );
00598 aNewTrack->SetWeight(weight);
00599 aNewTrack->SetTouchableHandle(trTouchable);
00600 aParticleChange.AddSecondary( aNewTrack );
00601 #ifdef debug
00602 G4cout<<"G4QElastic::PostStepDoIt: **** PostStepDoIt is done ****"<<G4endl;
00603 #endif
00604 return G4VDiscreteProcess::PostStepDoIt(track, step);
00605 }
00606
00607
00608 G4double G4QElastic::CalculateXSt(G4bool oxs, G4bool xst, G4double p, G4int Z, G4int N,
00609 G4int pPDG)
00610 {
00611 static G4bool init=false;
00612 static G4bool first=true;
00613 static G4VQCrossSection* PCSmanager;
00614 static G4VQCrossSection* NCSmanager;
00615 static G4VQCrossSection* PiPCSmanager;
00616 static G4VQCrossSection* PiMCSmanager;
00617 static G4VQCrossSection* KaPCSmanager;
00618 static G4VQCrossSection* KaMCSmanager;
00619 static G4VQCrossSection* HypCSmanager;
00620 static G4VQCrossSection* HyPCSmanager;
00621 static G4VQCrossSection* aBaCSmanager;
00622 if(first)
00623 {
00624 PCSmanager =G4QProtonElasticCrossSection::GetPointer();
00625 NCSmanager =G4QNeutronElasticCrossSection::GetPointer();
00626 PiPCSmanager=G4QPionPlusElasticCrossSection::GetPointer();
00627 PiMCSmanager=G4QPionMinusElasticCrossSection::GetPointer();
00628 KaPCSmanager=G4QKaonPlusElasticCrossSection::GetPointer();
00629 KaMCSmanager=G4QKaonMinusElasticCrossSection::GetPointer();
00630 HypCSmanager=G4QHyperonElasticCrossSection::GetPointer();
00631 HyPCSmanager=G4QHyperonPlusElasticCrossSection::GetPointer();
00632 aBaCSmanager=G4QAntiBaryonElasticCrossSection::GetPointer();
00633 first=false;
00634 }
00635 G4double res=0.;
00636 if(oxs && xst)
00637 {
00638 if (pPDG==2212) res=PCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00639 else if(pPDG==2112) res=NCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00640 else if(pPDG== 211) res=PiPCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00641 else if(pPDG==-211) res=PiMCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00642 else if(pPDG== 321) res=KaPCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00643 else if(pPDG==-321) res=KaMCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00644 else if(pPDG==310||pPDG==130) res=(KaMCSmanager->GetCrossSection(true, p, Z, N, pPDG)+
00645 KaPCSmanager->GetCrossSection(true, p, Z, N, pPDG))/2.;
00646 else if(pPDG==3222) res=HyPCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00647 else if(pPDG>3110 && pPDG<3335) res=HypCSmanager->GetCrossSection(true, p, Z, N, pPDG);
00648 else if(pPDG>-3335 && pPDG<-1110) res=aBaCSmanager->GetCrossSection(true,p,Z,N,pPDG);
00649 else G4cout<<"*Warning*G4QElastic::CalculateXSt: (o) wrong PDG="<<pPDG<<G4endl;
00650 }
00651 else if(!oxs && xst)
00652 {
00653 if (pPDG==2212) res=PCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00654 else if(pPDG==2112) res=NCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00655 else if(pPDG== 211) res=PiPCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00656 else if(pPDG==-211) res=PiMCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00657 else if(pPDG== 321) res=KaPCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00658 else if(pPDG==-321) res=KaMCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00659 else if(pPDG==310||pPDG==130) res=(KaMCSmanager->GetCrossSection(false, p, Z, N, pPDG)+
00660 KaPCSmanager->GetCrossSection(false, p, Z, N, pPDG))/2.;
00661 else if(pPDG==3222) res=HyPCSmanager->GetCrossSection(false, p, Z, N, pPDG);
00662 else if(pPDG>3110 && pPDG<3335) res=HypCSmanager->GetCrossSection(false, p, Z,N, pPDG);
00663 else if(pPDG>-3335 && pPDG<-1110) res=aBaCSmanager->GetCrossSection(false,p,Z,N,pPDG);
00664 else G4cout<<"*Warning*G4QElastic::CalculateXSt: (x) wrong PDG="<<pPDG<<G4endl;
00665
00666 init=true;
00667 }
00668 else if(init)
00669 {
00670 if (pPDG==2212)
00671 {
00672 if(oxs) res=PCSmanager->GetHMaxT();
00673 else res=PCSmanager->GetExchangeT(Z, N, pPDG);
00674 }
00675 else if(pPDG==2112)
00676 {
00677 if(oxs) res=NCSmanager->GetHMaxT();
00678 else res=NCSmanager->GetExchangeT(Z, N, pPDG);
00679 }
00680 else if(pPDG== 211)
00681 {
00682 if(oxs) res=PiPCSmanager->GetHMaxT();
00683 else res=PiPCSmanager->GetExchangeT(Z, N, pPDG);
00684 }
00685 else if(pPDG==-211)
00686 {
00687 if(oxs) res=PiMCSmanager->GetHMaxT();
00688 else res=PiMCSmanager->GetExchangeT(Z, N, pPDG);
00689 }
00690 else if(pPDG==321 || pPDG==310 || pPDG==130)
00691 {
00692 if(oxs) res=KaPCSmanager->GetHMaxT();
00693 else res=KaPCSmanager->GetExchangeT(Z, N, pPDG);
00694 }
00695 else if(pPDG==-321)
00696 {
00697 if(oxs) res=KaMCSmanager->GetHMaxT();
00698 else res=KaMCSmanager->GetExchangeT(Z, N, pPDG);
00699 }
00700 else if(pPDG==3222)
00701 {
00702 if(oxs) res=HyPCSmanager->GetHMaxT();
00703 else res=HyPCSmanager->GetExchangeT(Z, N, pPDG);
00704 }
00705 else if(pPDG<3335 && pPDG>1110)
00706 {
00707 if(oxs) res=HypCSmanager->GetHMaxT();
00708 else res=HypCSmanager->GetExchangeT(Z, N, pPDG);
00709 }
00710 else if(pPDG>-3335 && pPDG<-1110)
00711 {
00712 if(oxs) res=aBaCSmanager->GetHMaxT();
00713 else res=aBaCSmanager->GetExchangeT(Z, N, pPDG);
00714 }
00715 else G4cout<<"*Warning*G4QElastic::CalculateXSt: (i) wrong PDG="<<pPDG<<G4endl;
00716 }
00717 else G4cout<<"*Warning*G4QElastic::CalculateXSt:*NotInitiatedScattering*"<<G4endl;
00718 return res;
00719 }