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 #include "G4QCaptureAtRest.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4HadronicDeprecate.hh"
00052
00053
00054 G4QCaptureAtRest::G4QCaptureAtRest(const G4String& processName)
00055 : G4VRestProcess(processName, fHadronic), Time(0.), EnergyDeposition(0.)
00056 {
00057
00058 G4HadronicDeprecate("G4QCaptureAtRest");
00059
00060 #ifdef debug
00061 G4cout<<"G4QCaptureAtRest::Constructor is called"<<G4endl;
00062 #endif
00063 if (verboseLevel>0) G4cout << GetProcessName() << " is created "<< G4endl;
00064
00065
00066 G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio);
00067 G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);
00068 G4QEnvironment::SetParameters(SolidAngle);
00069 }
00070
00071 G4bool G4QCaptureAtRest::manualFlag=false;
00072 G4double G4QCaptureAtRest::Temperature=180.;
00073 G4double G4QCaptureAtRest::SSin2Gluons=0.3;
00074 G4double G4QCaptureAtRest::EtaEtaprime=0.3;
00075 G4double G4QCaptureAtRest::freeNuc=0.04;
00076 G4double G4QCaptureAtRest::freeDib=0.08;
00077 G4double G4QCaptureAtRest::clustProb=4.;
00078 G4double G4QCaptureAtRest::mediRatio=10.;
00079
00080 G4int G4QCaptureAtRest::nPartCWorld=85;
00081 G4double G4QCaptureAtRest::SolidAngle=0.5;
00082 G4bool G4QCaptureAtRest::EnergyFlux=false;
00083 G4double G4QCaptureAtRest::PiPrThresh=141.4;
00084 G4double G4QCaptureAtRest::M2ShiftVir=20000.;
00085 G4double G4QCaptureAtRest::DiNuclMass=1880.;
00086
00087 void G4QCaptureAtRest::SetManual() {manualFlag=true;}
00088 void G4QCaptureAtRest::SetStandard() {manualFlag=false;}
00089
00090
00091 void G4QCaptureAtRest::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
00092 G4double fN, G4double fD, G4double cP, G4double mR,
00093 G4int nParCW, G4double solAn, G4bool efFlag,
00094 G4double piThresh, G4double mpisq, G4double dinum)
00095 {
00096 Temperature=temper;
00097 SSin2Gluons=ssin2g;
00098 EtaEtaprime=etaetap;
00099 freeNuc=fN;
00100 freeDib=fD;
00101 clustProb=cP;
00102 mediRatio=mR;
00103 nPartCWorld = nParCW;
00104 EnergyFlux=efFlag;
00105 SolidAngle=solAn;
00106 PiPrThresh=piThresh;
00107 M2ShiftVir=mpisq;
00108 DiNuclMass=dinum;
00109 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
00110 G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio);
00111 G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);
00112 G4QEnvironment::SetParameters(SolidAngle);
00113 }
00114
00115
00116
00117 G4QCaptureAtRest::~G4QCaptureAtRest()
00118 {}
00119
00120 G4LorentzVector G4QCaptureAtRest::GetEnegryMomentumConservation()
00121 {
00122 return EnMomConservation;
00123 }
00124
00125 G4int G4QCaptureAtRest::GetNumberOfNeutronsInTarget()
00126 {
00127 return nOfNeutrons;
00128 }
00129
00130 G4bool G4QCaptureAtRest::IsApplicable(const G4ParticleDefinition& particle)
00131 {
00132 if (particle == *( G4PionMinus::PionMinus() )) return true;
00133 else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
00134 else if (particle == *( G4AntiProton::AntiProton() )) return true;
00135 else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
00136 else if (particle == *( G4TauMinus::TauMinus() )) return true;
00137 else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
00138 else if (particle == *( G4XiMinus::XiMinus() )) return true;
00139 else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
00140 else if (particle == *( G4Neutron::Neutron() )) return true;
00141 else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
00142 else if (particle == *(G4AntiSigmaPlus::AntiSigmaPlus())) return true;
00143 #ifdef debug
00144 G4cout<<"***G4QCaptureAtRest::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
00145 #endif
00146 return false;
00147 }
00148
00149 G4VParticleChange* G4QCaptureAtRest::AtRestDoIt(const G4Track& track, const G4Step& step)
00150 {
00151 static const G4double mNeut = G4QPDGCode(2112).GetMass();
00152 static const G4double mNeut2= mNeut*mNeut;
00153 static const G4double dmNeut= mNeut+mNeut;
00154 static const G4double mProt = G4QPDGCode(2212).GetMass();
00155 static const G4double dmProt= mProt+mProt;
00156 static const G4double mPi0 = G4QPDGCode(111).GetMass();
00157 static const G4double mDeut = G4QPDGCode(2112).GetNuclMass(1,1,0);
00158 static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0);
00159
00160 static const G4double mMu = G4QPDGCode(13).GetMass();
00161
00162
00163
00164 static const G4double b1=2.784;
00165 static const G4double ga=.000015;
00166 static const G4double rb=1./b1;
00167
00168 static const G4double mEl = G4QPDGCode(11).GetMass();
00169 static const G4double mEl2 = mEl*mEl;
00170
00171 static G4bool CWinit = true;
00172 if(CWinit)
00173 {
00174 CWinit=false;
00175 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld);
00176 }
00177
00178 const G4DynamicParticle* stoppedHadron = track.GetDynamicParticle();
00179 const G4ParticleDefinition* particle=stoppedHadron->GetDefinition();
00180 Time=0.;
00181 EnergyDeposition=0.;
00182 #ifdef debug
00183 G4cout<<"G4QCaptureAtRest::AtRestDoIt is called,EnDeposition="<<EnergyDeposition<<G4endl;
00184 #endif
00185 if (!IsApplicable(*particle))
00186 {
00187 G4cerr<<"G4QCaptureAtRest::AtRestDoIt: Only mu-,pi-,K-,S-,X-,O-,aP,aN,aS+."<< G4endl;
00188 return 0;
00189 }
00190 const G4Material* material = track.GetMaterial();
00191 G4int Z=0;
00192 const G4ElementVector* theElementVector = material->GetElementVector();
00193 G4int i=0;
00194 G4double sum=0.;
00195 G4int nE=material->GetNumberOfElements();
00196 #ifdef debug
00197 G4cout<<"G4QCaptureAtRest::AtRestDoIt: "<<nE<<" elements in the material."<<G4endl;
00198 #endif
00199
00200
00201 G4double primEnergy = 0.0;
00202 G4double secEnergy = 0.0;
00203 G4bool productOK = true;
00204 const G4double limEnergy = 2*GeV;
00205 G4int counter = 0;
00206 do {
00207 secEnergy = 0.0;
00208 ++counter;
00209 G4int countA = particle->GetBaryonNumber();
00210
00211
00212
00213 G4int projPDG=0;
00214 if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
00215 else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
00216 else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
00217 else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
00218 else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
00219 else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
00220 else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
00221 else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
00222 else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
00223 else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
00224 else if (particle == G4AntiSigmaPlus::AntiSigmaPlus()) projPDG=-3222;
00225 #ifdef debug
00226 G4cout<<"G4QCaptureAtRest::AtRestDoIt: projPDG="<<projPDG<<G4endl;
00227 #endif
00228 if(!projPDG)
00229 {
00230 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: Undefined captured hadron"<<G4endl;
00231 return 0;
00232 }
00233 std::vector<G4double> sumfra;
00234 for(i=0; i<nE; ++i)
00235 {
00236 G4double frac=material->GetFractionVector()[i];
00237 if(projPDG==13||projPDG==15) frac*=(*theElementVector)[i]->GetZ();
00238 sum+=frac;
00239 sumfra.push_back(sum);
00240 }
00241 G4double rnd = sum*G4UniformRand();
00242 for(i=0; i<nE; ++i) if (rnd<sumfra[i]) break;
00243 G4Element* pElement=(*theElementVector)[i];
00244 Z=static_cast<G4int>(pElement->GetZ());
00245 if(Z<=0)
00246 {
00247 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt:Element with Z="<<Z<< G4endl;
00248 if(Z<0) return 0;
00249 }
00250 G4QIsotope* Isotopes = G4QIsotope::Get();
00251 G4int N = Z;
00252 G4int isoSize=0;
00253 G4int indEl=0;
00254 G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
00255 if(isoVector) isoSize=isoVector->size();
00256 #ifdef debug
00257 G4cout<<"G4QCaptureAtRest::AtRestDoIt: isovectorLength="<<isoSize<<G4endl;
00258 #endif
00259 if(isoSize)
00260 {
00261 indEl=pElement->GetIndex()+1;
00262 #ifdef debug
00263 G4cout<<"G4QCapAR::GetMFP: iE="<<indEl<<", def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
00264 #endif
00265 if(!Isotopes->IsDefined(Z,indEl))
00266 {
00267 std::vector<std::pair<G4int,G4double>*>* newAbund =
00268 new std::vector<std::pair<G4int,G4double>*>;
00269 G4double* abuVector=pElement->GetRelativeAbundanceVector();
00270 for(G4int j=0; j<isoSize; j++)
00271 {
00272 N=pElement->GetIsotope(j)->GetN()-Z;
00273 if(pElement->GetIsotope(j)->GetZ()!=Z) G4cerr<<"*G4QCaptureAtRest::AtRestDoIt: Z="
00274 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
00275 G4double abund=abuVector[j];
00276 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
00277 #ifdef debug
00278 G4cout<<"G4QCaptureAtRest::AtRestDoIt:pair#="<<j<<", N="<<N<<",ab="<<abund<<G4endl;
00279 #endif
00280 newAbund->push_back(pr);
00281 }
00282 #ifdef debug
00283 G4cout<<"G4QCaptureAtRest::AtRestDoIt: pairVectorLength="<<newAbund->size()<<G4endl;
00284 #endif
00285 indEl=Isotopes->InitElement(Z,indEl,newAbund);
00286 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
00287 delete newAbund;
00288 }
00289
00290 N = Isotopes->GetNeutrons(Z,indEl);
00291 }
00292 else N = Isotopes->GetNeutrons(Z);
00293 nOfNeutrons=N;
00294 G4double dd=0.025;
00295 G4double am=Z+N;
00296 G4double value_sr=std::sqrt(am);
00297 G4double dsr=0.01*(value_sr+value_sr);
00298 if(dsr<dd)dsr=dd;
00299 if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio);
00300 else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);
00301 else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.);
00302 #ifdef debug
00303 G4cout<<"G4QCaptureAtRest::AtRestDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
00304 #endif
00305 if(N<0)
00306 {
00307 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt:Element with N="<<N<< G4endl;
00308 return 0;
00309 }
00310 G4bool lepChan=true;
00311 if(projPDG==13)
00312 {
00313 CalculateEnergyDepositionOfMuCapture(Z);
00314 lepChan=RandomizeMuDecayOrCapture(Z, N);
00315 }
00316 else if(projPDG==15)
00317 {
00318 CalculateEnergyDepositionOfTauCapture(Z);
00319 lepChan=RandomizeTauDecayOrCapture(Z,N);
00320 }
00321 G4double mp=G4QPDGCode(projPDG).GetMass();
00322 G4int targPDG=90000000+Z*1000+N;
00323
00324
00325 primEnergy = mp+G4QPDGCode(targPDG).GetMass();
00326
00327 G4QHadronVector* output=new G4QHadronVector;
00328 #ifdef debug
00329 G4cout<<"G4QCaptureAtRest::AtRestDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
00330 #endif
00331 G4double weight = track.GetWeight();
00332 G4double localtime = track.GetGlobalTime();
00333 G4ThreeVector position = track.GetPosition();
00334 #ifdef debug
00335 G4cout<<"G4QCaptureAtRest::AtRestDoIt: t="<<localtime<<", p="<<position<<G4endl;
00336 #endif
00337 G4TouchableHandle trTouchable = track.GetTouchableHandle();
00338 #ifdef debug
00339 G4cout<<"G4QCaptureAtRest::AtRestDoIt: touch="<<trTouchable<<G4endl;
00340 #endif
00341 localtime += Time;
00342 G4bool neutronElastic = false;
00343 G4bool chargExElastic = false;
00344 G4LorentzVector proj4M;
00345 G4double mAT=0.;
00346 G4double totNE=0.;
00347 G4double mAP=0.;
00348 if(projPDG==2112)
00349 {
00350 proj4M=stoppedHadron->Get4Momentum();
00351 totNE=proj4M.e();
00352
00353 G4QPDGCode QPDGbase;
00354 mAT=QPDGbase.GetNuclMass(Z,N,0);
00355 G4double mAR=std::sqrt(mAT*(mAT+totNE+totNE)+mNeut2);
00356 G4double mAN=QPDGbase.GetNuclMass(Z,N+1,0);
00357 mAP=QPDGbase.GetNuclMass(Z-1,N+1,0);
00358 G4double mAA=1000000.;
00359 if(Z>=2 && N>=1) mAA=QPDGbase.GetNuclMass(Z-2,N-1,0);
00360 G4double eProt=mAR-mAP-mProt;
00361 if(mAR<mAN && eProt>0.)
00362 {
00363 #ifdef debug
00364 G4cout<<"G4QCaptureAtRest::AtRestDoIt: n-Capture isn't possible mC="<<mAR<<" < mGS="
00365 <<mAN<<", Ep="<<eProt<<G4endl;
00366 #endif
00367 G4double eNeut=totNE-mNeut;
00368 if(eNeut<0.) eNeut=0.;
00369 if(eNeut<.0001) chargExElastic=true;
00370 else
00371 {
00372 G4double probP=std::sqrt(eProt*(dmProt+eProt));
00373 G4double probN=std::sqrt(eNeut*(dmNeut+eNeut));
00374 if((probN+probP)*G4UniformRand()<probN) neutronElastic=true;
00375 else chargExElastic=true;
00376 }
00377 }
00378 else if(mAR<=mAN||(mAR<=mAP+mProt&&mAR<=mAA+mAlph))
00379 {
00380 #ifdef debug
00381 G4cout<<"G4QCaptureAtRest::AtRestDoIt: n-Capture only elastic is possible"<<G4endl;
00382 #endif
00383 neutronElastic=true;
00384 }
00385 #ifdef debug
00386 G4cout<<"G4QCaptureAtRest::AtRestDoIt: n-Capture El="<<neutronElastic<<", Ex="
00387 <<chargExElastic<<G4endl;
00388 #endif
00389 }
00390 G4int nuPDG=14;
00391 if(projPDG==15) nuPDG=16;
00392 #ifdef debug
00393 G4int CV=0;
00394 G4cout<<"G4QCaptureAtRest::AtRestDoIt:DecayIf is reached CV="<<CV<<G4endl;
00395 #endif
00396 if(projPDG==2112 && neutronElastic)
00397 {
00398 #ifdef debug
00399 G4cout<<"G4QCaptureAtRest::AtRestDoIt:nElast, 4M="<<proj4M<<",Z="<<Z<<",N="<<N<<G4endl;
00400 #endif
00401 G4LorentzVector a4Mom(0.,0.,0.,mAT);
00402 G4LorentzVector totLV=proj4M+a4Mom;
00403 G4LorentzVector n4Mom(0.,0.,0.,mNeut);
00404 if(!G4QHadron(totLV).DecayIn2(a4Mom,n4Mom))
00405 {
00406 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt:n+A=>n+A,Z="<<Z<<",N="<<N<<G4endl;
00407 return 0;
00408 }
00409 G4QHadron* secnuc = new G4QHadron(targPDG,a4Mom);
00410 output->push_back(secnuc);
00411 G4QHadron* neutron = new G4QHadron(2112,n4Mom);
00412 output->push_back(neutron);
00413 #ifdef debug
00414 CV=27;
00415 G4cout<<"G4QCaptureAtRest::AtRestDoIt:ElasN="<<n4Mom<<",A="<<a4Mom<<",CV="<<CV<<G4endl;
00416 #endif
00417 }
00418 else if(projPDG==2112 && chargExElastic)
00419 {
00420 #ifdef debug
00421 G4cout<<"G4QCaptureAtRest::AtRestDoIt:npChEx, 4M="<<proj4M<<",Z="<<Z<<",N="<<N<<G4endl;
00422 #endif
00423 G4LorentzVector totLV(0.,0.,0.,mAT);
00424 totLV+=proj4M;
00425 G4LorentzVector a4Mom(0.,0.,0.,mAP);
00426 G4LorentzVector p4Mom(0.,0.,0.,mProt);
00427 if(!G4QHadron(totLV).DecayIn2(a4Mom,p4Mom))
00428 {
00429 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt:n+A=p+A',Z="<<Z<<",N="<<N<<G4endl;
00430 return 0;
00431 }
00432 G4QHadron* secnuc = new G4QHadron(targPDG-999,a4Mom);
00433 output->push_back(secnuc);
00434 G4QHadron* proton = new G4QHadron(2212,p4Mom);
00435 output->push_back(proton) ;
00436 #ifdef debug
00437 CV=21;
00438 G4cout<<"G4QCaptureAtRest::AtRestDoIt:ChExP="<<p4Mom<<",A="<<a4Mom<<",CV="<<CV<<G4endl;
00439 #endif
00440 }
00441 else if(projPDG==-211 && targPDG==90001000)
00442 {
00443 #ifdef debug
00444 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Panofsky targPDG="<<targPDG<<G4endl;
00445 #endif
00446 G4LorentzVector totLV(0.,0.,0.,mp+mProt);
00447 G4int pigamPDG=111;
00448 G4double pigamM=mPi0;
00449 if(G4UniformRand()>0.6)
00450 {
00451 pigamPDG=22;
00452 pigamM=0.;
00453 }
00454 G4LorentzVector g4Mom(0.,0.,0.,pigamM);
00455 G4LorentzVector n4Mom(0.,0.,0.,mNeut);
00456 if(!G4QHadron(totLV).DecayIn2(g4Mom,n4Mom))
00457 {
00458 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: H1+(pi-)=>n+"<<pigamPDG<<G4endl;
00459 return 0;
00460 }
00461 G4QHadron* pigamma = new G4QHadron(pigamPDG,g4Mom);
00462 output->push_back(pigamma);
00463 G4QHadron* neutron = new G4QHadron(2112,n4Mom);
00464 output->push_back(neutron);
00465 }
00466
00467
00468 else if(projPDG==-211 && G4UniformRand()>1.&& Z>0&&N>0)
00469 {
00470 G4double mt=G4QPDGCode(targPDG).GetMass();
00471 G4LorentzVector totLV(0.,0.,0.,mp+mt);
00472 if(Z==1 && N==1)
00473 {
00474 G4LorentzVector f4Mom(0.,0.,0.,mNeut);
00475 G4LorentzVector s4Mom(0.,0.,0.,mNeut);
00476 if(!G4QHadron(totLV).DecayIn2(f4Mom,s4Mom))
00477 {
00478 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: H2+(pi-)=>n+n"<<G4endl;
00479 return 0;
00480 }
00481 G4QHadron* neutr1 = new G4QHadron(2112,f4Mom);
00482 output->push_back(neutr1);
00483 G4QHadron* neutr2 = new G4QHadron(2112,s4Mom);
00484 output->push_back(neutr2);
00485 }
00486 else
00487 {
00488 G4int rPDG=targPDG-1001;
00489 G4double mr=G4QPDGCode(rPDG).GetMass();
00490 G4LorentzVector f4Mom(0.,0.,0.,mNeut);
00491 G4LorentzVector s4Mom(0.,0.,0.,mNeut);
00492 G4LorentzVector r4Mom(0.,0.,0.,mr);
00493 if(!G4QHadron(totLV).DecayIn3(f4Mom,s4Mom,r4Mom))
00494 {
00495 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: A+(pi-)=>n+n+(A-n-p)"<<G4endl;
00496 return 0;
00497 }
00498 G4QHadron* neutr1 = new G4QHadron(2112,f4Mom);
00499 output->push_back(neutr1);
00500 G4QHadron* neutr2 = new G4QHadron(2112,s4Mom);
00501 output->push_back(neutr2);
00502 G4QHadron* resnuc = new G4QHadron(rPDG,r4Mom);
00503 output->push_back(resnuc);
00504 }
00505 }
00506 else if((projPDG==13||projPDG==15) && !lepChan)
00507 {
00508 G4double mbm=mp-EnergyDeposition;
00509 G4LorentzVector totLV(0.,0.,0.,mbm);
00510 if(projPDG==13)
00511 {
00512 #ifdef debug
00513 G4cout<<"G4QCaptureAtRest::AtRestDoIt: e+nu+nu decay 4M="<<totLV<<totLV.m()<<G4endl;
00514 #endif
00515 G4double mt=G4QPDGCode(targPDG).GetMass();
00516 G4double ee=mEl+RandomizeDecayElectron(Z);
00517 totLV+=G4LorentzVector(0.,0.,0.,mt);
00518 G4double mmt=totLV.m();
00519 if(ee>=mbm*(mmt-mbm/2)/mmt)
00520 {
00521 G4cout<<"-W-G4QCaptureAtRest::AtRestDoIt: Unrealistic E="<<ee<<", m="<<mMu<<G4endl;
00522 G4LorentzVector f4Mom(0.,0.,0.,mEl);
00523 G4LorentzVector s4Mom(0.,0.,0.,mt);
00524 if(!G4QHadron(totLV).DecayIn2(f4Mom,s4Mom))
00525 {
00526 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: (mubound+A)->mu+A"<<G4endl;
00527 return 0;
00528 }
00529 G4QHadron* electron = new G4QHadron(11,f4Mom);
00530 output->push_back(electron);
00531 G4QHadron* resnuc = new G4QHadron(targPDG,s4Mom);
00532 output->push_back(resnuc);
00533 }
00534 else
00535 {
00536 G4double mr=std::sqrt(mmt*(mmt-ee-ee)+mEl2);
00537 if(mr<mt) mr=mt+.000001;
00538 G4LorentzVector f4Mom(0.,0.,0.,mEl);
00539 G4LorentzVector s4Mom(0.,0.,0.,mr);
00540 if(!G4QHadron(totLV).DecayIn2(f4Mom,s4Mom))
00541 {
00542 G4cerr<<"---Warning---G4QCaptureAtRest::AtRestDoIt: (mu+A)=>e+(A+2nu)"<<G4endl;
00543 return 0;
00544 }
00545 #ifdef debug
00546 G4double fe=f4Mom.e();
00547 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Ei="<<ee<<",Ef="<<fe<<",de="<<fe-ee<<G4endl;
00548 #endif
00549 G4QHadron* electron = new G4QHadron(11,f4Mom);
00550 output->push_back(electron);
00551 G4LorentzVector r4Mom(0.,0.,0.,mt);
00552 G4LorentzVector n4Mom(0.,0.,0.,0.);
00553 G4LorentzVector a4Mom(0.,0.,0.,0.);
00554 if(!G4QHadron(s4Mom).DecayIn3(r4Mom,n4Mom,a4Mom))
00555 {
00556 G4cerr<<"-Warning-G4QCaptureAtRest::AtRestDoIt: (A+2nu)=>A+NuMu+antiNuE"<<G4endl;
00557 return 0;
00558 }
00559 #ifdef debug
00560 G4cout<<"G4QCaptureAtRest::AtRestDoIt: (A+2nu) Decay is successful - 2"<<G4endl;
00561 #endif
00562 G4QHadron* resnuc = new G4QHadron(targPDG,r4Mom);
00563 #ifdef debug
00564 G4cout<<"G4QCaptureAtRest::AtRestDoIt: ResNuc 4M="<<mt<<r4Mom<<r4Mom.m()<<G4endl;
00565 #endif
00566 output->push_back(resnuc);
00567 #ifdef debug
00568 G4cout<<"G4QCaptureAtRest::AtRestDoIt: ResNuc is filled nu="<<n4Mom<<nuPDG<<G4endl;
00569 #endif
00570 G4QHadron* numu = new G4QHadron(nuPDG,n4Mom);
00571 output->push_back(numu);
00572 #ifdef debug
00573 G4cout<<"G4QCaptureAtRest::AtRestDoIt:Nu is filled,anu="<<a4Mom<<a4Mom.m()<<G4endl;
00574 #endif
00575 G4QHadron* anue = new G4QHadron(-12,a4Mom);
00576 output->push_back(anue);
00577 #ifdef debug
00578 G4cout<<"G4QCaptureAtRest::AtRestDoIt: anu is filled == Success of Mu-cap"<<G4endl;
00579 #endif
00580 }
00581 }
00582 else
00583 {
00584 G4int deL=11;
00585 G4int deN=-12;
00586 G4double mdl=mEl;
00587 if(G4UniformRand()>.55)
00588 {
00589 deL=13;
00590 deN=-14;
00591 mdl=mMu;
00592 }
00593 G4LorentzVector e4Mom(0.,0.,0.,mdl);
00594 G4LorentzVector n4Mom(0.,0.,0.,0.);
00595 G4LorentzVector a4Mom(0.,0.,0.,0.);
00596 if(!G4QHadron(totLV).DecayIn3(e4Mom,n4Mom,a4Mom))
00597 {
00598 G4cerr<<"--Worning--G4QCaptureAtRest::AtRestDoIt:Tau_b=>L+Nu_tau+anti_NuL"<<G4endl;
00599 return 0;
00600 }
00601 #ifdef debug
00602 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Tau Decay is successful"<<G4endl;
00603 #endif
00604 G4QHadron* lept = new G4QHadron(deL,e4Mom);
00605 #ifdef debug
00606 G4cout<<"G4QCaptureAtRest::AtRestDoIt: lepton 4M="<<e4Mom<<e4Mom.m()<<G4endl;
00607 #endif
00608 output->push_back(lept);
00609 #ifdef debug
00610 G4cout<<"G4QCaptureAtRest::AtRestDoIt: lepton is filled nu="<<n4Mom<<nuPDG<<G4endl;
00611 #endif
00612 G4QHadron* nuta = new G4QHadron(nuPDG,n4Mom);
00613 #ifdef debug
00614 G4cout<<"G4QCaptureAtRest::AtRestDoIt: nu 4M="<<n4Mom<<n4Mom.m()<<G4endl;
00615 #endif
00616 output->push_back(nuta);
00617 G4QHadron* anul = new G4QHadron(deN,a4Mom);
00618 #ifdef debug
00619 G4cout<<"G4QCaptureAtRest::AtRestDoIt: antiNu 4M="<<a4Mom<<a4Mom.m()<<G4endl;
00620 #endif
00621 output->push_back(anul);
00622 }
00623 }
00624 else if((projPDG==13||projPDG==15)&&lepChan&&targPDG==90001000)
00625 {
00626 G4LorentzVector totLV(0.,0.,0.,mp+mProt-EnergyDeposition);
00627 #ifdef debug
00628 G4cout<<"G4QCaptureAtRest::AtRestDoIt:CapOnProton decay 4M="<<totLV<<totLV.m()<<G4endl;
00629 #endif
00630 G4LorentzVector g4Mom(0.,0.,0.,0.);
00631 G4LorentzVector n4Mom(0.,0.,0.,mNeut);
00632 if(!G4QHadron(totLV).DecayIn2(g4Mom,n4Mom))
00633 {
00634 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: H1+(mu-)=>n+nu_mu"<<G4endl;
00635 return 0;
00636 }
00637 G4QHadron* neutrino = new G4QHadron(nuPDG,g4Mom);
00638 output->push_back(neutrino);
00639 G4QHadron* neutron = new G4QHadron(2112,n4Mom);
00640 output->push_back(neutron);
00641 }
00642 else if((projPDG==13||projPDG==15)&&lepChan&&targPDG==90001001)
00643 {
00644 G4LorentzVector totLV(0.,0.,0.,mp+mDeut-EnergyDeposition);
00645 #ifdef debug
00646 G4cout<<"G4QCaptureAtRest::AtRestDoIt: CapOnDeutr decay 4M="<<totLV<<totLV.m()<<G4endl;
00647 #endif
00648 G4LorentzVector g4Mom(0.,0.,0.,0.);
00649 G4LorentzVector n4Mom(0.,0.,0.,mNeut);
00650 G4LorentzVector s4Mom(0.,0.,0.,mNeut);
00651 if(!G4QHadron(totLV).DecayIn3(g4Mom,n4Mom,s4Mom))
00652 {
00653 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: D+(mu-)=>n+n+nu_mu"<<G4endl;
00654 return 0;
00655 }
00656 G4QHadron* neutrino = new G4QHadron(nuPDG,g4Mom);
00657 output->push_back(neutrino);
00658 G4QHadron* neut1 = new G4QHadron(2112,n4Mom);
00659 output->push_back(neut1);
00660 G4QHadron* neut2 = new G4QHadron(2112,s4Mom);
00661 output->push_back(neut2);
00662 }
00663
00664 else if((projPDG==13||projPDG==15)&&lepChan&&G4UniformRand()>1&&Z>0&&N>0)
00665 {
00666 G4double mt=G4QPDGCode(targPDG).GetMass();
00667 G4LorentzVector totLV(0.,0.,0.,mp+mt-EnergyDeposition);
00668 #ifdef debug
00669 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Quasi-Free decay 4M="<<totLV<<totLV.m()<<G4endl;
00670 #endif
00671 G4int rPDG=targPDG-1000;
00672 G4double mr=G4QPDGCode(rPDG).GetMass();
00673 G4LorentzVector f4Mom(0.,0.,0.,0.);
00674 G4LorentzVector s4Mom(0.,0.,0.,mNeut);
00675 G4LorentzVector r4Mom(0.,0.,0.,mr);
00676 if(!G4QHadron(totLV).DecayIn3(f4Mom,s4Mom,r4Mom))
00677 {
00678 G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: A+(mu-)=>nu_mu+n+(A-p)"<<G4endl;
00679 return 0;
00680 }
00681 G4QHadron* neutrino = new G4QHadron(nuPDG,f4Mom);
00682 output->push_back(neutrino);
00683 G4QHadron* neutron = new G4QHadron(2112,s4Mom);
00684 output->push_back(neutron);
00685 G4QHadron* resnuc = new G4QHadron(rPDG,r4Mom);
00686 output->push_back(resnuc);
00687 }
00688 else
00689
00690 {
00691 G4QHadron* neutr = 0;
00692 if(projPDG==13||projPDG==15) mp-=EnergyDeposition;
00693 #ifdef debug
00694 G4cout<<"G4QCaptureAtRest::AtRestDoIt: CHIPS M="<<mp<<",dE="<<EnergyDeposition<<G4endl;
00695 #endif
00696 G4LorentzVector projLV(0.,0.,0.,mp);
00697 if(projPDG==13)
00698 {
00699 if(G4UniformRand()<.04)
00700 {
00701 projPDG=-211;
00702
00703 G4double mmu2=projLV.m2();
00704 G4double mmu=std::sqrt(mmu2);
00705 G4double hmm=mmu/2.;
00706 G4double dmm=mmu+mmu;
00707 G4double qp=std::pow((std::pow(1.+ga*std::pow(hmm,b1),G4UniformRand())-1.)/ga,rb);
00708 G4double mqq=0.;
00709 if(qp<hmm) mqq=std::sqrt(mmu2-dmm*qp);
00710 G4LorentzVector f4Mom(0.,0.,0.,0.);
00711 G4LorentzVector s4Mom(0.,0.,0.,mqq);
00712 if(!G4QHadron(projLV).DecayIn2(f4Mom,s4Mom))
00713 {
00714 G4cerr<<"---Warning---G4QCaptureAtRest::AtRestDoIt: (mu-)=>q+aq+nu, M#1"<<G4endl;
00715 return 0;
00716 }
00717 neutr = new G4QHadron(nuPDG,f4Mom);
00718 projLV-=f4Mom;
00719
00720 }
00721 }
00722 G4QPDGCode targQPDG(targPDG);
00723 G4double tM=mp+targQPDG.GetMass();
00724 EnMomConservation=G4LorentzVector(0.,0.,0.,tM);
00725
00726 #ifdef tdebug
00727 G4cout<<"=--=>G4QCapAR:E/MCons, p="<<mp<<","<<projPDG<<",t="<<tM<<","<<targPDG<<",t4M="
00728 <<EnMomConservation<<G4endl;
00729 #endif
00730 G4QHadron* pH = new G4QHadron(projPDG,projLV);
00731 G4QHadronVector projHV;
00732 projHV.push_back(pH);
00733 G4QEnvironment* pan= new G4QEnvironment(projHV,targPDG);
00734 std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00735 projHV.clear();
00736 #ifdef debug
00737 G4cout<<"G4QCaptureAtRest::AtRestDoIt: pPDG="<<projPDG<<", m="<<mp<<G4endl;
00738 #endif
00739 try
00740 {
00741 delete output;
00742 output = pan->Fragment();
00743 }
00744 catch (G4QException& error)
00745 {
00746
00747 G4cerr<<"***G4QCaptureAtRest::AtRestDoIt: Exception is catched"<<G4endl;
00748
00749 G4Exception("G4QCaptureAtRest::AtRestDoIt()","HAD_CHPS_0027",
00750 FatalException," General CHIPS Exception");
00751 }
00752 delete pan;
00753 #ifdef debug
00754 G4cout<<"G4QCaptureAtRest::AtRestDoIt: CHIPS fragmentation is done, CV="<<CV<<G4endl;
00755 #endif
00756 if(neutr) output->push_back(neutr);
00757 }
00758 aParticleChange.Initialize(track);
00759 G4int tNH = output->size();
00760 aParticleChange.SetNumberOfSecondaries(tNH);
00761
00762 #ifdef debug
00763 G4cout<<"G4QCaptureAtRest::AtRestDoIt: "<<tNH<<" particles are generated"<<G4endl;
00764 #endif
00765
00766
00767 for(i=0; i<tNH; i++)
00768 {
00769
00770
00771
00772 G4QHadron* hadr=output->operator[](i);
00773 if(hadr->GetNFragments())
00774 {
00775 #ifdef debug
00776 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Intermediate particle is found i="<<i<<G4endl;
00777 #endif
00778 delete hadr;
00779 continue;
00780 }
00781
00782 countA -= hadr->GetBaryonNumber();
00783
00784
00785 G4DynamicParticle* theSec = new G4DynamicParticle;
00786 G4int PDGCode = hadr->GetPDGCode();
00787 #ifdef pdebug
00788 G4cout<<"G4QCaptureAtRest::AtRestDoIt:#"<<i<<",PDG="<<PDGCode<<G4endl;
00789 #endif
00790 G4ParticleDefinition* theDefinition=0;
00791 if (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
00792 else if(PDGCode==90001000) theDefinition = G4Proton::Proton();
00793 else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
00794 else if(PDGCode==311 || PDGCode==-311)
00795 {
00796 if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong();
00797 else theDefinition = G4KaonZeroShort::KaonZeroShort();
00798 }
00799 else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
00800 else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
00801 else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
00802 else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
00803 else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
00804 else if(PDGCode >80000000)
00805 {
00806 G4int aZ = hadr->GetCharge();
00807 G4int aA = hadr->GetBaryonNumber();
00808 #ifdef pdebug
00809 G4cout<<"G4QCaptureAtRest::AtRestDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
00810 #endif
00811
00812
00813
00814
00815
00816 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
00817 }
00818 else
00819 {
00820 #ifdef pdebug
00821 G4cout<<"G4QCaptureAtRest::AtRestDoIt:Define particle with PDG="<<PDGCode<<G4endl;
00822 #endif
00823 theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
00824 #ifdef pdebug
00825 G4cout<<"G4QCaptureAtRest::AtRestDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
00826 #endif
00827 }
00828 if(!theDefinition)
00829 {
00830 G4cout<<"---Worning---G4QCaptureAtRest::AtRestDoIt: drop PDG="<<PDGCode<<G4endl;
00831 delete hadr;
00832 continue;
00833 }
00834 #ifdef pdebug
00835 G4cout<<"G4QCaptureAtRest::AtRestDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
00836 #endif
00837 theSec->SetDefinition(theDefinition);
00838 G4LorentzVector h4M=hadr->Get4Momentum();
00839 EnMomConservation-=h4M;
00840
00841
00842 secEnergy += h4M.e();
00843 if(h4M.e() > 200*GeV) {
00844 G4cout<<"G4QCaptureAtRest: wrong product"<<i<<" imax= "<<tNH
00845 <<" " << theDefinition->GetParticleName() << " 4-mom= "
00846 << h4M <<G4endl;
00847 }
00848
00849
00850 #ifdef tdebug
00851 G4cout<<"G4QCapAR::ARDoIt:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
00852 #endif
00853 #ifdef debug
00854 G4cout<<"G4QCaptureAtRest::AtRestDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
00855 #endif
00856 theSec->Set4Momentum(h4M);
00857 delete hadr;
00858 #ifdef debug
00859 G4ThreeVector curD=theSec->GetMomentumDirection();
00860 G4double curM=theSec->GetMass();
00861 G4double curE=theSec->GetKineticEnergy()+curM;
00862 G4cout<<"G4QCapAtRest::AtRDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;
00863 #endif
00864 G4Track* aNewTrack = new G4Track(theSec, localtime, position );
00865 aNewTrack->SetWeight(weight);
00866 aNewTrack->SetTouchableHandle(trTouchable);
00867 aParticleChange.AddSecondary( aNewTrack );
00868 #ifdef debug
00869 G4cout<<"G4QCaptureAtRest::AtRestDoIt:#"<<i<<" is done."<<G4endl;
00870 #endif
00871 }
00872 delete output;
00873
00874 productOK = true;
00875 countA += Z + N;
00876 secEnergy += neutron_mass_c2*countA;
00877 if(EnergyDeposition > limEnergy || std::fabs(secEnergy - primEnergy) > limEnergy) {
00878 productOK = false;
00879 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Big energy non-concervation we need redo sampling"<<G4endl;
00880 G4cout << " Z= " << Z << " N= " << N << G4endl;
00881 G4cout<<"G4QCaptureAtRest::AtRestDoIt: the EnergyDeposition(GeV)="<<EnergyDeposition/GeV<<G4endl;
00882 G4cout<<" primEnergy(GeV)= " <<primEnergy/GeV << " secEnergy(GeV)="<< secEnergy/GeV <<G4endl;
00883 for(i=0; i<tNH; i++)
00884 {
00885 delete aParticleChange.GetSecondary(i);
00886 }
00887 aParticleChange.Clear();
00888 if(counter >= 100) {
00889 G4cout<<"G4QCaptureAtRest::AtRestDoIt: Cannot sample final state after "
00890 << counter << " iterations" <<G4endl;
00891 G4Exception("G4QCaptureAtRest::AtRestDoIt:","VI",FatalException,"Cannot sample final state");
00892 }
00893 }
00894 } while(!productOK);
00895 #ifdef debug
00896 G4cout<<"G4QCaptureAtRest::AtRestDoIt: the EnergyDeposition="<<EnergyDeposition<<G4endl;
00897 #endif
00898 aParticleChange.ProposeLocalEnergyDeposit(EnergyDeposition);
00899 aParticleChange.ProposeTrackStatus(fStopAndKill);
00900
00901 return G4VRestProcess::AtRestDoIt(track, step);
00902 }
00903
00904
00905 G4double G4QCaptureAtRest::GetMeanLifeTime(const G4Track& aTrack, G4ForceCondition*)
00906 {
00907 const G4DynamicParticle* stoppedHadron = aTrack.GetDynamicParticle();
00908 #ifdef debug
00909 G4cout<<"G4QCaptureAtRest::GetMeanLifeTime is called"<<G4endl;
00910 #endif
00911 if (*(stoppedHadron->GetDefinition())==*(G4MuonMinus::MuonMinus()) ||
00912 *(stoppedHadron->GetDefinition())==*(G4TauMinus::TauMinus()) ) return Time;
00913 else return 0.;
00914 }
00915
00916
00917 G4bool G4QCaptureAtRest::RandomizeMuDecayOrCapture(G4int Z, G4int N)
00918 {
00919 static G4int mZ=0;
00920 static G4int mN=0;
00921 static G4double mH=0.;
00922 static G4double mR=0.;
00923 static const G4int nAZ=17;
00924 static const G4int nZm=10;
00925 static const G4int rin[nZm]={0,0,1,1,1,2,3,4,5,6};
00926 static const G4double rate[nAZ]={.00000045, .00000047, .00000215, .000000356, .00000468,
00927 .00000226, .00000610, .00002750, .000023500, .00003790,
00928 .00003500, .00006600, .00006200, .000102500, .00009500,
00929 .00008800, .00022900};
00930 #ifdef debug
00931 G4cout<<"G4QCaptureAtRest::RandomizeMuDecayOrCapture is called"<<G4endl;
00932 #endif
00933 G4double z=Z;
00934 G4double Huff=1.;
00935 if(Z==mZ && N==mN) Huff=mH;
00936 else if(Z==8 || Z==9) Huff=.998;
00937
00938 else if(Z>9) Huff=1.-.000394*std::pow(z,2.19)/(1.+12.18*std::exp(z*.01373));
00939 G4double pD=.00045516*Huff;
00940 G4double pC=1.e99;
00941
00942 if(Z==mZ && N==mN) pC=mR;
00943 else if(Z>9) pC=.0000001256*std::pow(z,3.5)/(1.+.00429*std::exp(1.67*std::pow(z,.39)));
00944 else if(Z>0) pC=rate[rin[Z]+N];
00945 else G4cout<<"--Warning--G4QCaptureAtRest::RandomizeMuDecayOrCapture: NEG Z="<<Z<<G4endl;
00946 mZ=Z; mN=N; mH=Huff; mR=pC;
00947
00948
00949
00950
00951
00952 #ifdef debug
00953
00954 #endif
00955
00956
00957
00958
00959
00960 #ifdef debug
00961
00962 #endif
00963
00964
00965 if((pD+pC)*G4UniformRand()>pD)
00966 {
00967 Time=-std::log(G4UniformRand())/pC;
00968 return true;
00969 }
00970 else
00971 {
00972 Time=-std::log(G4UniformRand())/pD;
00973 return false;
00974 }
00975 }
00976
00977
00978 void G4QCaptureAtRest::CalculateEnergyDepositionOfMuCapture(G4int Z)
00979 {
00980 EnergyDeposition = Z*Z/(3556.+403.4/Z/(2.15+.0039*Z));
00981 #ifdef debug
00982 G4cout<<"G4QCaptureAtR::CalculateEnergyDepositionOfMuCapture="<<EnergyDeposition<<G4endl;
00983 #endif
00984 }
00985
00986
00987 G4bool G4QCaptureAtRest::RandomizeTauDecayOrCapture(G4int Z, G4int N)
00988 {
00989 #ifdef debug
00990 G4cout<<"G4QCaptureAtRest::RandomizeMuDecayOrCapture is called"<<G4endl;
00991 #endif
00992 G4double Z27 =0.002727*Z;
00993 G4double Z227=Z27*Z27;
00994 G4double Z427=Z227*Z227;
00995 G4double Zeff=(Z-0.13782)*(1.2162-(0.09118-Z427)*std::sqrt((G4double)Z));
00996 G4double Ze2=Zeff*Zeff;
00997 G4double pD=3436.*(1.-Ze2*.00014658);
00998 G4double pC=227.*Ze2*Ze2/(33.563+N);
00999 if(Z==1&&N==0) pC=10.;
01000 if(Z==1&&N==1) pC=.2;
01001 G4double DLifeT=-std::log(G4UniformRand())/pD;
01002 G4double CLifeT=-std::log(G4UniformRand())/pC;
01003 if(DLifeT<CLifeT)
01004 {
01005 Time=DLifeT;
01006 #ifdef debug
01007 G4cout<<"G4QCaptureAtRest::RandomizeTauDecayOrCapture: DecayLifeTime="<<Time<<G4endl;
01008 #endif
01009 return false;
01010 }
01011 else
01012 {
01013 Time=CLifeT;
01014 #ifdef debug
01015 G4cout<<"G4QCaptureAtRest::RandomizeTauDecayOrCapture: CaptureLifeTime="<<Time<<G4endl;
01016 #endif
01017 return true;
01018 }
01019 }
01020
01021
01022 void G4QCaptureAtRest::CalculateEnergyDepositionOfTauCapture(G4int Z)
01023 {
01024 EnergyDeposition = Z*Z/(21.144+40.38/Z/(2.15+.0039*Z));
01025 #ifdef debug
01026 G4cout<<"G4QCapAtRest::CalculateEnergyDepositionOfTauCapture="<<EnergyDeposition<<G4endl;
01027 #endif
01028 }
01029
01030
01031 G4double G4QCaptureAtRest::RandomizeDecayElectron(G4int z)
01032 {
01033 static const G4int nZ=12;
01034 static const G4int nE=200;
01035 static const G4int nEl=nE-1;
01036 static const G4int nEb=nE-2;
01037 static G4double tZ[nZ]={1.,2.,4.,6.,8.,11.,15.,20.,30.,45.,65.,92.};
01038
01039 static G4double P0[nE]={
01040 0.00000,7.28739,9.21005,10.5820,11.6892,12.6350,13.4701,14.2238,14.9144,15.5546,
01041 16.1533,16.7171,17.2511,17.7593,18.2450,18.7106,19.1584,19.5901,20.0073,20.4112,
01042 20.8031,21.1839,21.5544,21.9153,22.2675,22.6114,22.9476,23.2766,23.5988,23.9146,
01043 24.2243,24.5284,24.8270,25.1205,25.4091,25.6931,25.9726,26.2479,26.5192,26.7865,
01044 27.0502,27.3103,27.5670,27.8204,28.0706,28.3178,28.5621,28.8035,29.0422,29.2783,
01045 29.5118,29.7428,29.9715,30.1978,30.4219,30.6439,30.8637,31.0815,31.2972,31.5111,
01046 31.7231,31.9332,32.1416,32.3482,32.5531,32.7564,32.9581,33.1583,33.3569,33.5540,
01047 33.7497,33.9439,34.1368,34.3283,34.5185,34.7074,34.8951,35.0815,35.2667,35.4507,
01048 35.6335,35.8152,35.9958,36.1753,36.3537,36.5311,36.7075,36.8828,37.0572,37.2306,
01049 37.4030,37.5745,37.7451,37.9147,38.0835,38.2515,38.4185,38.5847,38.7501,38.9147,
01050 39.0785,39.2415,39.4038,39.5652,39.7260,39.8859,40.0452,40.2038,40.3616,40.5188,
01051 40.6753,40.8311,40.9862,41.1407,41.2946,41.4478,41.6004,41.7524,41.9038,42.0546,
01052 42.2048,42.3545,42.5035,42.6520,42.8000,42.9474,43.0942,43.2405,43.3863,43.5316,
01053 43.6764,43.8206,43.9644,44.1076,44.2504,44.3927,44.5345,44.6759,44.8168,44.9572,
01054 45.0972,45.2367,45.3758,45.5145,45.6527,45.7905,45.9279,46.0648,46.2014,46.3375,
01055 46.4733,46.6086,46.7436,46.8781,47.0123,47.1461,47.2796,47.4126,47.5453,47.6776,
01056 47.8096,47.9412,48.0724,48.2034,48.3339,48.4641,48.5940,48.7236,48.8528,48.9817,
01057 49.1103,49.2385,49.3665,49.4941,49.6214,49.7484,49.8751,50.0015,50.1276,50.2534,
01058 50.3789,50.5041,50.6290,50.7536,50.8780,51.0021,51.1259,51.2494,51.3727,51.4957,
01059 51.6184,51.7409,51.8633,51.9856,52.1080,52.2311,52.3561,52.4857,52.6273,52.8065};
01060
01061 static G4double P1[nE]={
01062 0.00000,7.27695,9.19653,10.5662,11.6715,12.6156,13.4493,14.2016,14.8910,15.5300,
01063 16.1275,16.6903,17.2233,17.7306,18.2153,18.6800,19.1269,19.5578,19.9742,20.3774,
01064 20.7685,21.1485,21.5182,21.8785,22.2299,22.5732,22.9087,23.2370,23.5585,23.8737,
01065 24.1828,24.4862,24.7843,25.0772,25.3652,25.6486,25.9275,26.2023,26.4729,26.7398,
01066 27.0029,27.2625,27.5186,27.7715,28.0212,28.2679,28.5116,28.7525,28.9907,29.2263,
01067 29.4593,29.6899,29.9180,30.1439,30.3675,30.5890,30.8083,31.0256,31.2410,31.4544,
01068 31.6659,31.8756,32.0835,32.2897,32.4942,32.6970,32.8983,33.0980,33.2962,33.4929,
01069 33.6881,33.8820,34.0744,34.2655,34.4553,34.6438,34.8310,35.0170,35.2018,35.3854,
01070 35.5678,35.7492,35.9293,36.1085,36.2865,36.4635,36.6395,36.8144,36.9884,37.1614,
01071 37.3334,37.5046,37.6748,37.8441,38.0125,38.1800,38.3467,38.5126,38.6776,38.8418,
01072 39.0053,39.1679,39.3298,39.4909,39.6512,39.8109,39.9698,40.1280,40.2855,40.4423,
01073 40.5984,40.7539,40.9087,41.0629,41.2164,41.3693,41.5215,41.6732,41.8242,41.9747,
01074 42.1246,42.2739,42.4226,42.5707,42.7184,42.8654,43.0119,43.1579,43.3034,43.4483,
01075 43.5928,43.7367,43.8801,44.0231,44.1655,44.3075,44.4490,44.5900,44.7306,44.8707,
01076 45.0104,45.1496,45.2884,45.4267,45.5646,45.7021,45.8392,45.9758,46.1121,46.2479,
01077 46.3833,46.5184,46.6530,46.7873,46.9211,47.0546,47.1878,47.3205,47.4529,47.5849,
01078 47.7166,47.8479,47.9788,48.1094,48.2397,48.3696,48.4992,48.6285,48.7574,48.8860,
01079 49.0143,49.1422,49.2699,49.3972,49.5242,49.6509,49.7773,49.9034,50.0292,50.1548,
01080 50.2800,50.4050,50.5297,50.6541,50.7783,50.9023,51.0261,51.1498,51.2735,51.3973,
01081 51.5215,51.6462,51.7721,51.8999,52.0308,52.1668,52.3115,52.4720,52.6642,52.9396};
01082
01083 static G4double P2[nE]={
01084 0.00000,7.23479,9.14334,10.5051,11.6039,12.5425,13.3712,14.1191,14.8044,15.4396,
01085 16.0336,16.5930,17.1228,17.6270,18.1088,18.5707,19.0150,19.4432,19.8571,20.2579,
01086 20.6466,21.0243,21.3918,21.7499,22.0992,22.4404,22.7739,23.1002,23.4198,23.7330,
01087 24.0403,24.3419,24.6381,24.9292,25.2155,25.4971,25.7744,26.0474,26.3165,26.5817,
01088 26.8432,27.1012,27.3558,27.6071,27.8553,28.1005,28.3428,28.5823,28.8190,29.0531,
01089 29.2847,29.5139,29.7407,29.9652,30.1874,30.4076,30.6256,30.8416,31.0556,31.2677,
01090 31.4779,31.6863,31.8930,32.0979,32.3012,32.5028,32.7029,32.9014,33.0984,33.2939,
01091 33.4880,33.6806,33.8719,34.0619,34.2505,34.4379,34.6240,34.8088,34.9925,35.1750,
01092 35.3563,35.5366,35.7157,35.8937,36.0707,36.2466,36.4215,36.5954,36.7683,36.9403,
01093 37.1113,37.2814,37.4506,37.6189,37.7863,37.9528,38.1185,38.2834,38.4475,38.6107,
01094 38.7731,38.9348,39.0957,39.2559,39.4153,39.5740,39.7319,39.8892,40.0457,40.2016,
01095 40.3568,40.5114,40.6652,40.8185,40.9711,41.1231,41.2744,41.4252,41.5754,41.7249,
01096 41.8739,42.0223,42.1702,42.3174,42.4642,42.6104,42.7560,42.9012,43.0458,43.1899,
01097 43.3334,43.4765,43.6191,43.7612,43.9028,44.0440,44.1846,44.3248,44.4646,44.6039,
01098 44.7427,44.8811,45.0191,45.1566,45.2937,45.4304,45.5666,45.7025,45.8379,45.9730,
01099 46.1076,46.2419,46.3757,46.5092,46.6423,46.7750,46.9074,47.0394,47.1710,47.3023,
01100 47.4332,47.5637,47.6939,47.8238,47.9533,48.0826,48.2114,48.3400,48.4683,48.5962,
01101 48.7239,48.8512,48.9783,49.1052,49.2318,49.3581,49.4843,49.6103,49.7362,49.8620,
01102 49.9878,50.1137,50.2397,50.3660,50.4927,50.6200,50.7481,50.8774,51.0083,51.1412,
01103 51.2769,51.4163,51.5607,51.7120,51.8729,52.0476,52.2435,52.4742,52.7705,53.2313};
01104
01105 static G4double P3[nE]={
01106 0.00000,7.18201,9.07814,10.4312,11.5230,12.4557,13.2793,14.0225,14.7036,15.3349,
01107 15.9253,16.4814,17.0080,17.5092,17.9882,18.4474,18.8890,19.3148,19.7263,20.1247,
01108 20.5112,20.8867,21.2521,21.6082,21.9555,22.2948,22.6264,22.9509,23.2687,23.5802,
01109 23.8858,24.1857,24.4803,24.7699,25.0546,25.3347,25.6105,25.8821,26.1497,26.4135,
01110 26.6737,26.9303,27.1836,27.4336,27.6805,27.9244,28.1654,28.4037,28.6392,28.8721,
01111 29.1026,29.3306,29.5562,29.7796,30.0007,30.2198,30.4367,30.6516,30.8646,31.0756,
01112 31.2848,31.4923,31.6979,31.9019,32.1041,32.3048,32.5039,32.7015,32.8975,33.0921,
01113 33.2853,33.4770,33.6674,33.8565,34.0443,34.2308,34.4160,34.6000,34.7829,34.9645,
01114 35.1451,35.3245,35.5028,35.6800,35.8562,36.0313,36.2055,36.3786,36.5508,36.7220,
01115 36.8923,37.0616,37.2301,37.3976,37.5643,37.7302,37.8952,38.0593,38.2227,38.3852,
01116 38.5470,38.7080,38.8683,39.0277,39.1865,39.3445,39.5018,39.6585,39.8144,39.9696,
01117 40.1242,40.2781,40.4314,40.5840,40.7360,40.8874,41.0382,41.1884,41.3379,41.4869,
01118 41.6353,41.7832,41.9305,42.0772,42.2234,42.3690,42.5141,42.6587,42.8028,42.9463,
01119 43.0894,43.2319,43.3740,43.5156,43.6567,43.7973,43.9375,44.0772,44.2164,44.3552,
01120 44.4936,44.6315,44.7690,44.9060,45.0426,45.1789,45.3147,45.4501,45.5851,45.7197,
01121 45.8539,45.9877,46.1211,46.2542,46.3869,46.5193,46.6512,46.7829,46.9142,47.0451,
01122 47.1758,47.3061,47.4361,47.5658,47.6952,47.8243,47.9532,48.0819,48.2103,48.3385,
01123 48.4666,48.5945,48.7223,48.8501,48.9778,49.1056,49.2335,49.3616,49.4899,49.6187,
01124 49.7479,49.8779,50.0087,50.1406,50.2738,50.4087,50.5457,50.6854,50.8282,50.9752,
01125 51.1272,51.2858,51.4528,51.6310,51.8244,52.0390,52.2854,52.5832,52.9766,53.6078};
01126
01127 static G4double P4[nE]={
01128 0.00000,7.11944,9.00194,10.3456,11.4301,12.3566,13.1748,13.9133,14.5901,15.2176,
01129 15.8044,16.3570,16.8806,17.3789,17.8550,18.3116,18.7507,19.1741,19.5833,19.9796,
01130 20.3640,20.7375,21.1010,21.4552,21.8007,22.1382,22.4681,22.7910,23.1072,23.4172,
01131 23.7213,24.0198,24.3130,24.6011,24.8845,25.1633,25.4378,25.7082,25.9746,26.2372,
01132 26.4962,26.7517,27.0039,27.2528,27.4987,27.7415,27.9816,28.2188,28.4534,28.6854,
01133 28.9149,29.1420,29.3667,29.5892,29.8095,30.0277,30.2438,30.4579,30.6701,30.8803,
01134 31.0888,31.2954,31.5004,31.7036,31.9052,32.1051,32.3035,32.5004,32.6958,32.8898,
01135 33.0823,33.2734,33.4632,33.6517,33.8388,34.0247,34.2094,34.3929,34.5751,34.7563,
01136 34.9363,35.1151,35.2929,35.4696,35.6453,35.8200,35.9936,36.1663,36.3380,36.5088,
01137 36.6786,36.8475,37.0155,37.1827,37.3489,37.5144,37.6789,37.8427,38.0057,38.1678,
01138 38.3292,38.4899,38.6497,38.8089,38.9673,39.1250,39.2819,39.4382,39.5938,39.7488,
01139 39.9030,40.0566,40.2096,40.3619,40.5136,40.6647,40.8152,40.9651,41.1144,41.2632,
01140 41.4113,41.5589,41.7059,41.8524,41.9984,42.1438,42.2887,42.4330,42.5769,42.7202,
01141 42.8631,43.0054,43.1473,43.2887,43.4296,43.5701,43.7101,43.8496,43.9887,44.1274,
01142 44.2656,44.4034,44.5407,44.6777,44.8142,44.9503,45.0861,45.2214,45.3564,45.4910,
01143 45.6252,45.7590,45.8925,46.0256,46.1584,46.2909,46.4231,46.5549,46.6865,46.8178,
01144 46.9488,47.0796,47.2101,47.3404,47.4706,47.6005,47.7304,47.8601,47.9898,48.1194,
01145 48.2491,48.3788,48.5087,48.6388,48.7692,48.8999,49.0311,49.1629,49.2955,49.4289,
01146 49.5635,49.6993,49.8367,49.9760,50.1175,50.2617,50.4091,50.5605,50.7166,50.8785,
01147 51.0475,51.2256,51.4151,51.6195,51.8439,52.0962,52.3897,52.7496,53.2325,54.0205};
01148
01149 static G4double P5[nE]={
01150 0.00000,7.00695,8.86633,10.1943,11.2667,12.1832,12.9928,13.7237,14.3938,15.0151,
01151 15.5963,16.1438,16.6625,17.1563,17.6283,18.0810,18.5164,18.9362,19.3421,19.7352,
01152 20.1165,20.4872,20.8479,21.1994,21.5424,21.8775,22.2051,22.5257,22.8397,23.1476,
01153 23.4497,23.7462,24.0375,24.3239,24.6055,24.8826,25.1555,25.4243,25.6892,25.9503,
01154 26.2078,26.4619,26.7127,26.9604,27.2050,27.4466,27.6854,27.9215,28.1550,28.3859,
01155 28.6143,28.8404,29.0641,29.2857,29.5050,29.7223,29.9375,30.1508,30.3621,30.5716,
01156 30.7793,30.9852,31.1894,31.3919,31.5928,31.7921,31.9899,32.1862,32.3810,32.5743,
01157 32.7663,32.9569,33.1462,33.3341,33.5208,33.7063,33.8905,34.0735,34.2554,34.4361,
01158 34.6157,34.7942,34.9717,35.1481,35.3235,35.4978,35.6712,35.8436,36.0150,36.1856,
01159 36.3552,36.5239,36.6917,36.8586,37.0247,37.1900,37.3544,37.5181,37.6809,37.8430,
01160 38.0043,38.1648,38.3246,38.4837,38.6420,38.7997,38.9566,39.1129,39.2685,39.4234,
01161 39.5777,39.7314,39.8844,40.0368,40.1885,40.3397,40.4903,40.6403,40.7897,40.9385,
01162 41.0868,41.2345,41.3817,41.5284,41.6745,41.8201,41.9652,42.1097,42.2538,42.3974,
01163 42.5405,42.6831,42.8253,42.9670,43.1082,43.2490,43.3893,43.5292,43.6687,43.8078,
01164 43.9464,44.0847,44.2225,44.3600,44.4971,44.6338,44.7702,44.9062,45.0418,45.1772,
01165 45.3122,45.4469,45.5814,45.7155,45.8494,45.9831,46.1165,46.2497,46.3828,46.5156,
01166 46.6484,46.7810,46.9136,47.0462,47.1787,47.3113,47.4439,47.5768,47.7098,47.8430,
01167 47.9767,48.1107,48.2452,48.3804,48.5163,48.6530,48.7908,48.9297,49.0700,49.2119,
01168 49.3556,49.5015,49.6499,49.8012,49.9559,50.1146,50.2780,50.4470,50.6226,50.8063,
01169 50.9997,51.2054,51.4263,51.6671,51.9343,52.2382,52.5959,53.0400,53.6438,54.6445};
01170
01171 static G4double P6[nE]={
01172 0.00000,6.82422,8.64751,9.95140,11.0052,11.9065,12.7032,13.4228,14.0829,14.6952,
01173 15.2682,15.8082,16.3200,16.8074,17.2734,17.7205,18.1506,18.5656,18.9668,19.3555,
01174 19.7327,20.0995,20.4564,20.8044,21.1441,21.4759,21.8005,22.1182,22.4294,22.7347,
01175 23.0342,23.3283,23.6173,23.9014,24.1809,24.4560,24.7270,24.9939,25.2570,25.5164,
01176 25.7723,26.0249,26.2742,26.5204,26.7637,27.0040,27.2416,27.4765,27.7088,27.9387,
01177 28.1661,28.3912,28.6141,28.8347,29.0533,29.2698,29.4843,29.6969,29.9076,30.1165,
01178 30.3236,30.5290,30.7327,30.9348,31.1353,31.3342,31.5317,31.7277,31.9222,32.1153,
01179 32.3071,32.4975,32.6867,32.8745,33.0611,33.2465,33.4308,33.6138,33.7957,33.9765,
01180 34.1562,34.3349,34.5125,34.6890,34.8646,35.0392,35.2128,35.3855,35.5573,35.7281,
01181 35.8980,36.0671,36.2353,36.4027,36.5692,36.7350,36.8999,37.0640,37.2274,37.3900,
01182 37.5519,37.7130,37.8734,38.0331,38.1921,38.3505,38.5081,38.6651,38.8214,38.9771,
01183 39.1322,39.2866,39.4404,39.5937,39.7463,39.8983,40.0498,40.2007,40.3511,40.5009,
01184 40.6502,40.7989,40.9471,41.0948,41.2420,41.3887,41.5349,41.6806,41.8259,41.9707,
01185 42.1150,42.2589,42.4024,42.5454,42.6880,42.8302,42.9720,43.1134,43.2545,43.3951,
01186 43.5354,43.6754,43.8150,43.9542,44.0932,44.2319,44.3702,44.5084,44.6462,44.7838,
01187 44.9212,45.0584,45.1954,45.3323,45.4690,45.6057,45.7422,45.8787,46.0152,46.1517,
01188 46.2882,46.4249,46.5617,46.6987,46.8360,46.9735,47.1115,47.2499,47.3888,47.5284,
01189 47.6687,47.8098,47.9518,48.0950,48.2394,48.3852,48.5326,48.6819,48.8332,48.9869,
01190 49.1432,49.3027,49.4657,49.6327,49.8044,49.9815,50.1649,50.3557,50.5552,50.7653,
01191 50.9882,51.2267,51.4850,51.7686,52.0858,52.4496,52.8817,53.4232,54.1672,55.4144};
01192
01193 static G4double P7[nE]={
01194 0.00000,6.55304,8.32352,9.59220,10.6190,11.4982,12.2762,12.9795,13.6251,14.2244,
01195 14.7856,15.3149,15.8168,16.2951,16.7526,17.1918,17.6146,18.0227,18.4174,18.8000,
01196 19.1714,19.5327,19.8846,20.2277,20.5627,20.8902,21.2106,21.5244,21.8320,22.1337,
01197 22.4298,22.7208,23.0067,23.2880,23.5648,23.8373,24.1058,24.3704,24.6313,24.8886,
01198 25.1426,25.3933,25.6409,25.8855,26.1272,26.3661,26.6023,26.8360,27.0671,27.2959,
01199 27.5223,27.7465,27.9685,28.1884,28.4062,28.6221,28.8360,29.0481,29.2584,29.4669,
01200 29.6737,29.8789,30.0824,30.2844,30.4848,30.6838,30.8813,31.0774,31.2721,31.4654,
01201 31.6575,31.8482,32.0377,32.2260,32.4131,32.5990,32.7838,32.9675,33.1501,33.3316,
01202 33.5120,33.6915,33.8699,34.0474,34.2238,34.3994,34.5740,34.7478,34.9206,35.0926,
01203 35.2637,35.4339,35.6034,35.7720,35.9399,36.1070,36.2733,36.4389,36.6037,36.7678,
01204 36.9312,37.0939,37.2560,37.4173,37.5780,37.7380,37.8974,38.0562,38.2144,38.3719,
01205 38.5289,38.6853,38.8411,38.9963,39.1510,39.3051,39.4587,39.6117,39.7643,39.9163,
01206 40.0679,40.2189,40.3695,40.5196,40.6692,40.8184,40.9671,41.1154,41.2633,41.4108,
01207 41.5579,41.7045,41.8508,41.9967,42.1423,42.2875,42.4324,42.5769,42.7212,42.8651,
01208 43.0087,43.1521,43.2953,43.4381,43.5808,43.7233,43.8655,44.0077,44.1496,44.2915,
01209 44.4332,44.5749,44.7166,44.8582,44.9999,45.1416,45.2834,45.4253,45.5675,45.7098,
01210 45.8524,45.9953,46.1386,46.2824,46.4267,46.5716,46.7172,46.8635,47.0107,47.1589,
01211 47.3083,47.4588,47.6108,47.7643,47.9195,48.0768,48.2362,48.3981,48.5628,48.7306,
01212 48.9019,49.0772,49.2570,49.4420,49.6329,49.8306,50.0362,50.2510,50.4768,50.7155,
01213 50.9699,51.2436,51.5414,51.8702,52.2401,52.6667,53.1765,53.8195,54.7092,56.213};
01214
01215 static G4double P8[nE]={
01216 0.00000,5.93631,7.58588,8.77313,9.73699,10.5643,11.2979,11.9624,12.5734,13.1415,
01217 13.6742,14.1774,14.6552,15.1111,15.5477,15.9674,16.3719,16.7627,17.1412,17.5084,
01218 17.8654,18.2129,18.5517,18.8825,19.2058,19.5221,19.8319,20.1356,20.4335,20.7260,
01219 21.0135,21.2961,21.5742,21.8479,22.1175,22.3833,22.6453,22.9037,23.1588,23.4106,
01220 23.6593,23.9050,24.1479,24.3880,24.6255,24.8605,25.0931,25.3232,25.5512,25.7769,
01221 26.0005,26.2221,26.4417,26.6594,26.8753,27.0893,27.3017,27.5123,27.7214,27.9288,
01222 28.1348,28.3392,28.5422,28.7437,28.9439,29.1428,29.3404,29.5367,29.7318,29.9257,
01223 30.1184,30.3100,30.5005,30.6899,30.8782,31.0655,31.2518,31.4371,31.6215,31.8049,
01224 31.9874,32.1690,32.3497,32.5296,32.7086,32.8869,33.0643,33.2409,33.4168,33.5919,
01225 33.7663,33.9399,34.1129,34.2851,34.4567,34.6277,34.7979,34.9676,35.1366,35.3050,
01226 35.4728,35.6400,35.8067,35.9727,36.1383,36.3033,36.4677,36.6317,36.7951,36.9580,
01227 37.1205,37.2825,37.4440,37.6050,37.7656,37.9258,38.0855,38.2449,38.4038,38.5623,
01228 38.7205,38.8782,39.0356,39.1927,39.3494,39.5058,39.6618,39.8176,39.9730,40.1281,
01229 40.2830,40.4376,40.5920,40.7461,40.9000,41.0537,41.2072,41.3606,41.5137,41.6667,
01230 41.8196,41.9724,42.1251,42.2777,42.4303,42.5829,42.7354,42.8880,43.0407,43.1934,
01231 43.3462,43.4992,43.6524,43.8058,43.9595,44.1135,44.2679,44.4226,44.5779,44.7336,
01232 44.8900,45.0470,45.2047,45.3632,45.5227,45.6831,45.8447,46.0075,46.1716,46.3372,
01233 46.5045,46.6735,46.8446,47.0179,47.1936,47.3720,47.5535,47.7382,47.9267,48.1193,
01234 48.3166,48.5191,48.7275,48.9426,49.1653,49.3968,49.6383,49.8917,50.1588,50.4424,
01235 50.7458,51.0735,51.4316,51.8286,52.2771,52.7968,53.4209,54.2121,55.3131,57.1875};
01236
01237 static G4double P9[nE]={
01238 0.00000,5.00701,6.48187,7.55018,8.42075,9.17008,9.83596,10.4402,10.9967,11.5149,
01239 12.0017,12.4619,12.8995,13.3176,13.7185,14.1042,14.4765,14.8365,15.1856,15.5248,
01240 15.8548,16.1764,16.4903,16.7972,17.0974,17.3914,17.6797,17.9627,18.2406,18.5137,
01241 18.7825,19.0470,19.3075,19.5643,19.8176,20.0674,20.3141,20.5576,20.7983,21.0362,
01242 21.2714,21.5041,21.7344,21.9623,22.1880,22.4116,22.6331,22.8527,23.0703,23.2862,
01243 23.5003,23.7127,23.9234,24.1327,24.3404,24.5466,24.7514,24.9549,25.1571,25.3580,
01244 25.5577,25.7562,25.9535,26.1497,26.3449,26.5390,26.7322,26.9243,27.1155,27.3058,
01245 27.4952,27.6837,27.8714,28.0583,28.2445,28.4298,28.6144,28.7984,28.9816,29.1641,
01246 29.3460,29.5273,29.7079,29.8880,30.0675,30.2464,30.4248,30.6026,30.7799,30.9568,
01247 31.1332,31.3091,31.4845,31.6595,31.8341,32.0083,32.1821,32.3555,32.5286,32.7013,
01248 32.8736,33.0456,33.2173,33.3887,33.5598,33.7306,33.9011,34.0713,34.2414,34.4111,
01249 34.5806,34.7500,34.9191,35.0880,35.2567,35.4252,35.5936,35.7618,35.9299,36.0978,
01250 36.2656,36.4333,36.6008,36.7683,36.9357,37.1031,37.2704,37.4376,37.6048,37.7720,
01251 37.9391,38.1063,38.2735,38.4407,38.6080,38.7754,38.9428,39.1103,39.2779,39.4457,
01252 39.6136,39.7817,39.9500,40.1184,40.2872,40.4562,40.6255,40.7951,40.9650,41.1354,
01253 41.3061,41.4773,41.6490,41.8213,41.9941,42.1675,42.3417,42.5165,42.6922,42.8687,
01254 43.0462,43.2246,43.4042,43.5849,43.7669,43.9504,44.1353,44.3218,44.5102,44.7004,
01255 44.8928,45.0875,45.2847,45.4846,45.6876,45.8939,46.1038,46.3178,46.5362,46.7596,
01256 46.9885,47.2236,47.4656,47.7155,47.9743,48.2434,48.5242,48.8187,49.1292,49.4588,
01257 49.8113,50.1919,50.6075,51.0680,51.5878,52.1895,52.9113,53.8254,55.0955,57.2547};
01258
01259 static G4double PA[nE]={
01260 0.00000,3.72463,4.98461,5.91463,6.68054,7.34437,7.93719,8.47714,8.97585,9.44132,
01261 9.87929,10.2941,10.6890,11.0666,11.4290,11.7780,12.1150,12.4412,12.7576,13.0651,
01262 13.3644,13.6563,13.9413,14.2199,14.4925,14.7597,15.0217,15.2789,15.5315,15.7800,
01263 16.0244,16.2651,16.5023,16.7361,16.9667,17.1943,17.4190,17.6410,17.8604,18.0774,
01264 18.2920,18.5044,18.7147,18.9229,19.1292,19.3336,19.5362,19.7371,19.9364,20.1341,
01265 20.3303,20.5251,20.7185,20.9106,21.1014,21.2910,21.4794,21.6667,21.8529,22.0381,
01266 22.2223,22.4055,22.5879,22.7693,22.9500,23.1298,23.3088,23.4871,23.6647,23.8416,
01267 24.0179,24.1935,24.3685,24.5430,24.7170,24.8904,25.0633,25.2358,25.4078,25.5794,
01268 25.7506,25.9214,26.0919,26.2620,26.4319,26.6014,26.7706,26.9396,27.1084,27.2770,
01269 27.4453,27.6135,27.7815,27.9494,28.1171,28.2847,28.4522,28.6197,28.7871,28.9545,
01270 29.1218,29.2891,29.4564,29.6237,29.7911,29.9586,30.1260,30.2936,30.4613,30.6291,
01271 30.7970,30.9651,31.1333,31.3018,31.4704,31.6392,31.8082,31.9775,32.1471,32.3169,
01272 32.4870,32.6574,32.8281,32.9992,33.1706,33.3425,33.5147,33.6873,33.8603,34.0338,
01273 34.2078,34.3822,34.5572,34.7326,34.9087,35.0852,35.2624,35.4402,35.6187,35.7978,
01274 35.9775,36.1580,36.3393,36.5213,36.7041,36.8877,37.0722,37.2576,37.4439,37.6312,
01275 37.8196,38.0089,38.1994,38.3910,38.5838,38.7778,38.9731,39.1698,39.3680,39.5676,
01276 39.7688,39.9716,40.1762,40.3827,40.5910,40.8014,41.0140,41.2289,41.4463,41.6662,
01277 41.8890,42.1147,42.3437,42.5761,42.8122,43.0524,43.2969,43.5462,43.8008,44.0610,
01278 44.3276,44.6012,44.8826,45.1727,45.4727,45.7838,46.1077,46.4465,46.8024,47.1788,
01279 47.5795,48.0100,48.4774,48.9920,49.5687,50.2309,51.0179,52.0041,53.3576,55.6238};
01280
01281
01282 static G4double PB[nE]={
01283 0.00000,1.93933,2.80877,3.48857,4.06872,4.58456,5.05440,5.48917,5.89607,6.28010,
01284 6.64494,6.99335,7.32751,7.64913,7.95964,8.26019,8.55175,8.83516,9.11111,9.38023,
01285 9.64304,9.90002,10.1516,10.3981,10.6399,10.8772,11.1105,11.3398,11.5654,11.7875,
01286 12.0063,12.2220,12.4347,12.6446,12.8518,13.0564,13.2586,13.4585,13.6561,13.8516,
01287 14.0450,14.2364,14.4260,14.6137,14.7998,14.9841,15.1668,15.3480,15.5278,15.7060,
01288 15.8829,16.0585,16.2328,16.4059,16.5778,16.7485,16.9181,17.0867,17.2543,17.4208,
01289 17.5864,17.7511,17.9149,18.0779,18.2400,18.4013,18.5619,18.7218,18.8809,19.0394,
01290 19.1972,19.3544,19.5110,19.6670,19.8225,19.9774,20.1318,20.2858,20.4393,20.5924,
01291 20.7450,20.8973,21.0492,21.2007,21.3519,21.5028,21.6534,21.8038,21.9539,22.1037,
01292 22.2534,22.4029,22.5522,22.7013,22.8503,22.9992,23.1481,23.2968,23.4455,23.5942,
01293 23.7428,23.8915,24.0402,24.1889,24.3377,24.4866,24.6355,24.7846,24.9339,25.0833,
01294 25.2329,25.3828,25.5328,25.6832,25.8337,25.9846,26.1359,26.2874,26.4394,26.5917,
01295 26.7444,26.8976,27.0513,27.2055,27.3602,27.5155,27.6713,27.8278,27.9849,28.1427,
01296 28.3012,28.4604,28.6204,28.7812,28.9429,29.1055,29.2690,29.4334,29.5989,29.7655,
01297 29.9331,30.1019,30.2719,30.4431,30.6156,30.7895,30.9648,31.1416,31.3200,31.4999,
01298 31.6816,31.8650,32.0502,32.2374,32.4266,32.6178,32.8113,33.0071,33.2053,33.4060,
01299 33.6094,33.8156,34.0248,34.2370,34.4525,34.6714,34.8939,35.1203,35.3508,35.5855,
01300 35.8249,36.0691,36.3185,36.5734,36.8343,37.1015,37.3756,37.6571,37.9466,38.2448,
01301 38.5525,38.8705,39.2000,39.5422,39.8984,40.2703,40.6601,41.0701,41.5034,41.9638,
01302 42.4562,42.9870,43.5648,44.2017,44.9151,45.7322,46.6981,47.8972,49.5188,52.1676};
01303 static G4double* PP[nZ]= {P0,P1,P2,P3,P4,P5,P6,P7,P8,P9,PA,PB};
01304 static const G4int mZ1=93;
01305 static const G4int mZ=mZ1-1;
01306 static G4double* P[mZ1]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
01307 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
01308 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
01309
01310 if(z<1 || z>mZ)
01311 {
01312 G4cout<<"-W-G4QCapAtRest::RandomizeDecayElectron: <=0 or big(>"<<mZ<<") Z="<<z<<G4endl;
01313 return 0.;
01314 }
01315 G4double Z = z;
01316 G4double* nP = 0;
01317 if(!P[z])
01318 {
01319 for(G4int i=0; i<nZ; i++)
01320 {
01321 G4double fZ=tZ[i];
01322 if(Z<=fZ)
01323 {
01324 #ifdef debug
01325 G4cout<<"G4QCapAtRest::RandomizeDecayElectron: Z="<<Z<<", fZ="<<fZ<<G4endl;
01326 #endif
01327 nP = new G4double[nE];
01328 G4double* fP=PP[i];
01329 if(Z==fZ) for(G4int j=0; j<nE; j++) nP[j]=fP[j];
01330 else
01331 {
01332 G4int i1=i-1;
01333 G4double iZ=tZ[i1];
01334 G4double* iP=PP[i1];
01335 G4double rZ=(Z-iZ)/(fZ-iZ);
01336 for(G4int j=0; j<nE; j++) nP[j]=iP[j]+(fP[j]-iP[j])*rZ;
01337 }
01338 #ifdef debug
01339 for(G4int k=0; k<nE; k++)G4cout<<"G4QCAR::RandDecEle:P["<<k<<"]="<<nP[k]<<G4endl;
01340 #endif
01341 P[z]=nP;
01342 break;
01343 }
01344 }
01345 }
01346 else nP=P[z];
01347
01348 G4double R=G4UniformRand()*nE;
01349 G4int iR=static_cast<int>(R);
01350 if(iR > nEl) iR=nEl;
01351 else if(iR<0) iR=0;
01352 G4double nPi=nP[iR];
01353 G4double nPf=0.;
01354 if(iR<nEl) nPf=nP[iR+1];
01355 else nPf=nP[nEl]+(nP[nEl]-nP[nEb])/2;
01356 #ifdef debug
01357 G4cout<<"G4QCapAtR::RaDEl:R="<<R<<",Ei="<<nPi<<",E="<<MeV*(nPi+(R-iR)*(nPf-nPi))<<G4endl;
01358 #endif
01359 return MeV*(nPi+(R-iR)*(nPf-nPi));
01360 }