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
00046
00047
00048
00049 #include "G4QDiffraction.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4HadronicDeprecate.hh"
00052
00053
00054
00055
00056
00057 G4int G4QDiffraction::nPartCWorld=85;
00058 std::vector<G4int> G4QDiffraction::ElementZ;
00059 std::vector<G4double> G4QDiffraction::ElProbInMat;
00060 std::vector<std::vector<G4int>*> G4QDiffraction::ElIsoN;
00061 std::vector<std::vector<G4double>*>G4QDiffraction::IsoProbInEl;
00062
00063
00064 G4QDiffraction::G4QDiffraction(const G4String& processName):
00065 G4VDiscreteProcess(processName, fHadronic)
00066 {
00067 G4HadronicDeprecate("G4QDiffraction");
00068
00069 #ifdef debug
00070 G4cout<<"G4QDiffraction::Constructor is called processName="<<processName<<G4endl;
00071 #endif
00072 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
00073 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
00074 }
00075
00076
00077 G4QDiffraction::~G4QDiffraction() {}
00078
00079
00080 G4LorentzVector G4QDiffraction::GetEnegryMomentumConservation(){return EnMomConservation;}
00081
00082 G4int G4QDiffraction::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
00083
00084
00085
00086
00087 G4double G4QDiffraction::GetMeanFreePath(const G4Track& Track,G4double,G4ForceCondition *F)
00088 {
00089 *F = NotForced;
00090 const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
00091 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
00092 if( !IsApplicable(*incidentParticleDefinition))
00093 G4cout<<"-Warning-G4QDiffraction::GetMeanFreePath for notImplemented Particle"<<G4endl;
00094
00095 G4double Momentum = incidentParticle->GetTotalMomentum();
00096 #ifdef debug
00097 G4double KinEn = incidentParticle->GetKineticEnergy();
00098 G4cout<<"G4QDiffraction::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
00099 #endif
00100 const G4Material* material = Track.GetMaterial();
00101 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
00102 const G4ElementVector* theElementVector = material->GetElementVector();
00103 G4int nE=material->GetNumberOfElements();
00104 #ifdef debug
00105 G4cout<<"G4QDiffraction::GetMeanFreePath:"<<nE<<" Elems in Material="<<*material<<G4endl;
00106 #endif
00107 G4int pPDG=0;
00108
00109 if (incidentParticleDefinition == G4Proton::Proton() ) pPDG=2212;
00110 else if(incidentParticleDefinition == G4Neutron::Neutron()) pPDG=2112;
00111 else G4cout<<"G4QDiffraction::GetMeanFreePath: only nA & pA are implemented"<<G4endl;
00112
00113 G4QIsotope* Isotopes = G4QIsotope::Get();
00114 G4double sigma=0.;
00115 G4int IPIE=IsoProbInEl.size();
00116 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)
00117 {
00118 std::vector<G4double>* SPI=IsoProbInEl[ip];
00119 SPI->clear();
00120 delete SPI;
00121 std::vector<G4int>* IsN=ElIsoN[ip];
00122 IsN->clear();
00123 delete IsN;
00124 }
00125 ElProbInMat.clear();
00126 ElementZ.clear();
00127 IsoProbInEl.clear();
00128 ElIsoN.clear();
00129 for(G4int i=0; i<nE; ++i)
00130 {
00131 G4Element* pElement=(*theElementVector)[i];
00132 G4int Z = static_cast<G4int>(pElement->GetZ());
00133 ElementZ.push_back(Z);
00134 G4int isoSize=0;
00135 G4int indEl=0;
00136 G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
00137 if(isoVector) isoSize=isoVector->size();
00138 #ifdef debug
00139 G4cout<<"G4QDiffraction::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
00140 #endif
00141 if(isoSize)
00142 {
00143 indEl=pElement->GetIndex()+1;
00144 #ifdef debug
00145 G4cout<<"G4QDiffr::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
00146 #endif
00147 if(!Isotopes->IsDefined(Z,indEl))
00148 {
00149 std::vector<std::pair<G4int,G4double>*>* newAbund =
00150 new std::vector<std::pair<G4int,G4double>*>;
00151 G4double* abuVector=pElement->GetRelativeAbundanceVector();
00152 for(G4int j=0; j<isoSize; j++)
00153 {
00154 G4int N=pElement->GetIsotope(j)->GetN()-Z;
00155 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
00156 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00157 G4double abund=abuVector[j];
00158 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00159 #ifdef debug
00160 G4cout<<"G4QDiffract::GetMeanFreePath:pair#"<<j<<",N="<<N<<",ab="<<abund<<G4endl;
00161 #endif
00162 newAbund->push_back(pr);
00163 }
00164 #ifdef debug
00165 G4cout<<"G4QDiffract::GetMeanFreePath:pairVectorLength="<<newAbund->size()<<G4endl;
00166 #endif
00167 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund);
00168 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
00169 delete newAbund;
00170 }
00171 }
00172 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);
00173 std::vector<G4double>* SPI = new std::vector<G4double>;
00174 IsoProbInEl.push_back(SPI);
00175 std::vector<G4int>* IsN = new std::vector<G4int>;
00176 ElIsoN.push_back(IsN);
00177 G4int nIs=cs->size();
00178 #ifdef debug
00179 G4cout<<"G4QDiffract::GetMFP:=***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
00180 #endif
00181 G4double susi=0.;
00182 if(nIs) for(G4int j=0; j<nIs; j++)
00183 {
00184 std::pair<G4int,G4double>* curIs=(*cs)[j];
00185 G4int N=curIs->first;
00186 IsN->push_back(N);
00187 #ifdef debug
00188 G4cout<<"G4Q::GMFP:j="<<j<<",P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PD="<<pPDG<<G4endl;
00189 #endif
00190 G4double CSI=CalculateXS(Momentum, Z, N, pPDG);
00191
00192 #ifdef debug
00193 G4cout<<"G4QDiffraction::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
00194 <<Momentu<<", XSec="<<CSI/millibarn<<G4endl;
00195 #endif
00196 curIs->second = CSI;
00197 susi+=CSI;
00198 SPI->push_back(susi);
00199 }
00200 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];
00201 #ifdef debug
00202 G4cout<<"G4QDiffraction::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
00203 <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
00204 #endif
00205 ElProbInMat.push_back(sigma);
00206 }
00207
00208 #ifdef debug
00209 G4cout<<"G4QDiffraction::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
00210 #endif
00211 if(sigma > 0.) return 1./sigma;
00212 return DBL_MAX;
00213 }
00214
00215 G4bool G4QDiffraction::IsApplicable(const G4ParticleDefinition& particle)
00216 {
00217 if (particle == *( G4Proton::Proton() )) return true;
00218 else if (particle == *( G4Neutron::Neutron() )) return true;
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 #ifdef debug
00244 G4cout<<"***>>G4QDiffraction::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<G4endl;
00245 #endif
00246 return false;
00247 }
00248
00249 G4VParticleChange* G4QDiffraction::PostStepDoIt(const G4Track& track, const G4Step& step)
00250 {
00251 static const G4double mProt= G4QPDGCode(2212).GetMass();
00252 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00253 static const G4double mPion= G4QPDGCode(111).GetMass();
00254 static G4QDiffractionRatio* diffRatio;
00255
00256
00257 static G4bool CWinit = true;
00258 if(CWinit)
00259 {
00260 CWinit=false;
00261 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
00262 diffRatio=G4QDiffractionRatio::GetPointer();
00263 }
00264
00265 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
00266 const G4ParticleDefinition* particle=projHadron->GetDefinition();
00267 #ifdef debug
00268 G4cout<<"G4QDiffraction::PostStepDoIt: Before the GetMeanFreePath is called In4M="
00269 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
00270 <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
00271 #endif
00272 G4ForceCondition cond=NotForced;
00273 GetMeanFreePath(track, -27., &cond);
00274 #ifdef debug
00275 G4cout<<"G4QDiffraction::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
00276 #endif
00277 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV;
00278 G4double momentum = projHadron->GetTotalMomentum()/MeV;
00279 G4double Momentum = proj4M.rho();
00280 if(std::fabs(Momentum-momentum)>.000001)
00281 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
00282 #ifdef pdebug
00283 G4cout<<"G4QDiffraction::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum
00284 <<",proj4M="<<proj4M<<", projM="<<proj4M.m()<<G4endl;
00285 #endif
00286 if (!IsApplicable(*particle))
00287 {
00288 G4cerr<<"G4QDiffraction::PostStepDoIt: Only NA is implemented."<<G4endl;
00289 return 0;
00290 }
00291 const G4Material* material = track.GetMaterial();
00292 G4int Z=0;
00293 const G4ElementVector* theElementVector = material->GetElementVector();
00294 G4int nE=material->GetNumberOfElements();
00295 #ifdef debug
00296 G4cout<<"G4QDiffraction::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
00297 #endif
00298 G4int projPDG=0;
00299
00300 if (particle == G4Proton::Proton() ) projPDG= 2212;
00301 else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 #ifdef debug
00331 G4int prPDG=particle->GetPDGEncoding();
00332 G4cout<<"G4QDiffraction::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
00333 #endif
00334 if(!projPDG)
00335 {
00336 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
00337 return 0;
00338 }
00339
00340
00341 G4double pM=mNeut;
00342 if(projPDG==2112) pM=mProt;
00343
00344 G4int EPIM=ElProbInMat.size();
00345 #ifdef debug
00346 G4cout<<"G4QDiffra::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
00347 #endif
00348 G4int i=0;
00349 if(EPIM>1)
00350 {
00351 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
00352 for(i=0; i<nE; ++i)
00353 {
00354 #ifdef debug
00355 G4cout<<"G4QDiffra::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
00356 #endif
00357 if (rnd<ElProbInMat[i]) break;
00358 }
00359 if(i>=nE) i=nE-1;
00360 }
00361 G4Element* pElement=(*theElementVector)[i];
00362 Z=static_cast<G4int>(pElement->GetZ());
00363 #ifdef debug
00364 G4cout<<"G4QDiffraction::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
00365 #endif
00366 if(Z<=0)
00367 {
00368 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Element with Z="<<Z<<G4endl;
00369 if(Z<0) return 0;
00370 }
00371 std::vector<G4double>* SPI = IsoProbInEl[i];
00372 std::vector<G4int>* IsN = ElIsoN[i];
00373 G4int nofIsot=SPI->size();
00374 #ifdef debug
00375 G4cout<<"G4QDiffraction::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
00376 #endif
00377 G4int j=0;
00378 if(nofIsot>1)
00379 {
00380 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand();
00381 for(j=0; j<nofIsot; ++j)
00382 {
00383 #ifdef debug
00384 G4cout<<"G4QDiffraction::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
00385 #endif
00386 if(rndI < (*SPI)[j]) break;
00387 }
00388 if(j>=nofIsot) j=nofIsot-1;
00389 }
00390 G4int N =(*IsN)[j]; ;
00391 #ifdef debug
00392 G4cout<<"G4QDiffraction::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<", MeV="<<MeV<<G4endl;
00393 #endif
00394 if(N<0)
00395 {
00396 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Isotope Z="<<Z<<" has 0>N="<<N<<G4endl;
00397 return 0;
00398 }
00399 nOfNeutrons=N;
00400 #ifdef debug
00401 G4cout<<"G4QDiffraction::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
00402 #endif
00403 if(N<0)
00404 {
00405 G4cerr<<"*Warning*G4QDiffraction::PostStepDoIt:Element with N="<<N<< G4endl;
00406 return 0;
00407 }
00408 aParticleChange.Initialize(track);
00409 #ifdef debug
00410 G4cout<<"G4QDiffraction::PostStepDoIt: track is initialized"<<G4endl;
00411 #endif
00412 G4double weight = track.GetWeight();
00413 G4double localtime = track.GetGlobalTime();
00414 G4ThreeVector position = track.GetPosition();
00415 #ifdef debug
00416 G4cout<<"G4QDiffraction::PostStepDoIt: before Touchable extraction"<<G4endl;
00417 #endif
00418 G4TouchableHandle trTouchable = track.GetTouchableHandle();
00419 #ifdef debug
00420 G4cout<<"G4QDiffraction::PostStepDoIt: Touchable is extracted"<<G4endl;
00421 #endif
00422 G4int targPDG=90000000+Z*1000+N;
00423 G4QPDGCode targQPDG(targPDG);
00424 G4double tM=targQPDG.GetMass();
00425 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV;
00426 G4ParticleMomentum dir = projHadron->GetMomentumDirection();
00427 G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM);
00428 #ifdef debug
00429 G4cout<<"G4QDiffraction::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
00430 #endif
00431 EnMomConservation=tot4M;
00432
00433 #ifdef debug
00434 G4cout<<"G4QDiff::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl;
00435 #endif
00436 G4double xSec=CalculateXS(Momentum, Z, N, projPDG);
00437 #ifdef debug
00438 G4cout<<"G4QDiffra::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",XS="<<xSec/millibarn<<G4endl;
00439 #endif
00440 #ifdef nandebug
00441 if(xSec>0. || xSec<0. || xSec==0);
00442 else G4cout<<"-Warning-G4QDiffraction::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
00443 #endif
00444
00445 if(xSec <= 0.)
00446 {
00447 #ifdef pdebug
00448 G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Zero cross-section* PDG="<<projPDG
00449 <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl;
00450 #endif
00451
00452 aParticleChange.ProposeEnergy(kinEnergy);
00453 aParticleChange.ProposeLocalEnergyDeposit(0.);
00454 aParticleChange.ProposeMomentumDirection(dir) ;
00455 return G4VDiscreteProcess::PostStepDoIt(track,step);
00456 }
00457 G4double totCMMass=tot4M.m();
00458 if(totCMMass < mPion+pM+tM)
00459 {
00460 #ifdef pdebug
00461 G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Below Diffraction Threshold* cmM="<<totCMMass
00462 <<">pM="<<pM<<"+tM="<<tM<<"+pi0="<<mPion<<"=="<<pM+tM+mPion<<G4endl;
00463 #endif
00464
00465 aParticleChange.ProposeEnergy(kinEnergy);
00466 aParticleChange.ProposeLocalEnergyDeposit(0.);
00467 aParticleChange.ProposeMomentumDirection(dir) ;
00468 return G4VDiscreteProcess::PostStepDoIt(track,step);
00469 }
00470
00471 aParticleChange.ProposeTrackStatus(fStopAndKill);
00472 G4QHadronVector* out=diffRatio->TargFragment(projPDG, proj4M, Z, N);
00473 G4int nSec=out->size();
00474 G4DynamicParticle* theSec=0;
00475 G4LorentzVector dif4M(0.,0.,0.,0.);
00476 G4int difPDG=0;
00477 G4QHadron* difQH=0;
00478 #ifdef pdebug
00479 G4cout<<"G4QDiffraction::PostStepDoIt: =---=found=---= nSecondaries="<<nSec<<G4endl;
00480 #endif
00481 for(i=0; i<nSec; i++)
00482 {
00483 difQH = (*out)[i];
00484 difPDG= difQH->GetPDGCode();
00485 G4ParticleDefinition* theDefinition=0;
00486 if (difPDG==2212 || difPDG==90001000) theDefinition=G4Proton::Proton();
00487 else if(difPDG==2112 || difPDG==90000001) theDefinition=G4Neutron::Neutron();
00488 else if(difPDG== 22) theDefinition=G4Gamma::Gamma();
00489 else if(difPDG== 111) theDefinition=G4PionZero::PionZero();
00490 else if(difPDG==-211 || difPDG==89999001) theDefinition=G4PionMinus::PionMinus();
00491 else if(difPDG== 211 || difPDG==90000999) theDefinition=G4PionPlus::PionPlus();
00492 else if(difPDG== 321 || difPDG==89001000) theDefinition=G4KaonPlus::KaonPlus();
00493 else if(difPDG==-321 || difPDG==90999000) theDefinition=G4KaonMinus::KaonMinus();
00494 else if(difPDG== 130 || difPDG==-311 || difPDG==89000001)
00495 theDefinition=G4KaonZeroLong::KaonZeroLong();
00496 else if(difPDG== 310 || difPDG== 311 || difPDG==90999999)
00497 theDefinition=G4KaonZeroShort::KaonZeroShort();
00498 else if(difPDG==3122 || difPDG==91000000) theDefinition=G4Lambda::Lambda();
00499 else if(difPDG== 3222) theDefinition=G4SigmaPlus::SigmaPlus();
00500 else if(difPDG== 3112) theDefinition=G4SigmaMinus::SigmaMinus();
00501 else if(difPDG== 3212) theDefinition=G4SigmaZero::SigmaZero();
00502 else if(difPDG== 3312) theDefinition=G4XiMinus::XiMinus();
00503 else if(difPDG== 3322) theDefinition=G4XiZero::XiZero();
00504 else if(difPDG== 3334) theDefinition=G4OmegaMinus::OmegaMinus();
00505 else if(difPDG==-2112) theDefinition=G4AntiNeutron::AntiNeutron();
00506 else if(difPDG==-2212) theDefinition=G4AntiProton::AntiProton();
00507 else if(difPDG==-3122) theDefinition=G4AntiLambda::AntiLambda();
00508 else if(difPDG==-3222) theDefinition=G4AntiSigmaPlus::AntiSigmaPlus();
00509 else if(difPDG==-3112) theDefinition=G4AntiSigmaMinus::AntiSigmaMinus();
00510 else if(difPDG==-3212) theDefinition=G4AntiSigmaZero::AntiSigmaZero();
00511 else if(difPDG==-3312) theDefinition=G4AntiXiMinus::AntiXiMinus();
00512 else if(difPDG==-3322) theDefinition=G4AntiXiZero::AntiXiZero();
00513 else if(difPDG==-3334) theDefinition=G4AntiOmegaMinus::AntiOmegaMinus();
00514 else if(difPDG== -11) theDefinition=G4Electron::Electron();
00515 else if(difPDG== -13) theDefinition=G4MuonMinus::MuonMinus();
00516 else if(difPDG== 11) theDefinition=G4Positron::Positron();
00517 else if(difPDG== 13) theDefinition=G4MuonPlus::MuonPlus();
00518 else
00519 {
00520 Z = difQH->GetCharge();
00521 G4int B = difQH->GetBaryonNumber();
00522 G4int S = difQH->GetStrangeness();
00523 if(S||Z>B||Z<0)G4cout<<"-Warning-G4QDif::PoStDoIt:Z="<<Z<<",A="<<B<<",S="<<S<<G4endl;
00524 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(Z,B,0,0);
00525 #ifdef pdebug
00526 G4cout<<"G4QDiffraction::PoStDoIt:Ion,Z="<<Z<<",A="<<B<<",D="<<theDefinition<<G4endl;
00527 #endif
00528 }
00529 if(theDefinition)
00530 {
00531 theSec = new G4DynamicParticle;
00532 theSec->SetDefinition(theDefinition);
00533 dif4M = difQH->Get4Momentum();
00534 EnMomConservation-=dif4M;
00535 theSec->Set4Momentum(dif4M);
00536 G4Track* aNewTrack = new G4Track(theSec, localtime, position );
00537 aNewTrack->SetWeight(weight);
00538 aNewTrack->SetTouchableHandle(trTouchable);
00539 aParticleChange.AddSecondary( aNewTrack );
00540 #ifdef pdebug
00541 G4cout<<"G4QDiffraction::PostStepDoIt: Filled 4M="<<dif4M<<", PDG="<<difPDG<<G4endl;
00542 #endif
00543 }
00544 else G4cout<<"-Warning-G4QDif::PSDI: Lost PDG="<<difPDG<<", Z="<<difQH->GetCharge()
00545 <<", A="<<difQH->GetBaryonNumber()<<",S ="<<difQH->GetStrangeness()<<G4endl;
00546 delete difQH;
00547 }
00548 delete out;
00549 #ifdef debug
00550 G4cout<<"G4QDiffraction::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
00551 #endif
00552 return G4VDiscreteProcess::PostStepDoIt(track, step);
00553 }
00554
00555 G4double G4QDiffraction::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG)
00556 {
00557 static G4bool first=true;
00558
00559 static G4QDiffractionRatio* diffRatio;
00560 if(first)
00561 {
00562
00563 diffRatio=G4QDiffractionRatio::GetPointer();
00564 first=false;
00565 }
00566
00567
00568
00569
00570 G4double xs_value=diffRatio->GetTargSingDiffXS(p, PDG, Z, N);
00571 #ifdef debug
00572 G4cout<<"G4QDiff::CXS:p="<<p<<",Z="<<Z<<",N="<<N<<",C="<<PDG<<",XS="<<xs_value<<G4endl;
00573 #endif
00574 return xs_value;
00575 }