#include <G4QInelastic.hh>
Inheritance diagram for G4QInelastic:
Definition at line 122 of file G4QInelastic.hh.
| G4QInelastic::G4QInelastic | ( | const G4String & | processName = "CHIPS_Inelastic" |
) |
Definition at line 57 of file G4QInelastic.cc.
References G4cout, G4endl, G4HadronicDeprecate, G4VProcess::GetProcessName(), G4QEnvironment::SetParameters(), G4Quasmon::SetParameters(), G4QNucleus::SetParameters(), and G4VProcess::verboseLevel.
00057 : 00058 G4VDiscreteProcess(processName, fHadronic) 00059 { 00060 G4HadronicDeprecate("G4QInelastic"); 00061 00062 EnMomConservation = G4LorentzVector(0.,0.,0.,0.); 00063 nOfNeutrons = 0; 00064 #ifdef debug 00065 G4cout<<"G4QInelastic::Constructor is called"<<G4endl; 00066 #endif 00067 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl; 00068 //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPSWorld (234 part.max) 00069 G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's 00070 G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters 00071 G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture 00072 //@@ Initialize here the G4QuasmonString parameters 00073 }
| G4QInelastic::~G4QInelastic | ( | ) |
Definition at line 127 of file G4QInelastic.cc.
00128 { 00129 // The following is just a copy of what is done in PostStepDoIt every interaction ! 00130 // The correction is if(IPIE), so just for(...;ip<IPIE;...) does not work ! @@ 00131 G4int IPIE=IsoProbInEl.size(); // How many old elements? 00132 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) 00133 { 00134 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector 00135 SPI->clear(); 00136 delete SPI; 00137 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector 00138 IsN->clear(); 00139 delete IsN; 00140 } 00141 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) 00142 ElementZ.clear(); // Clear the body vector for Z of Elements 00143 IsoProbInEl.clear(); // Clear the body vector for SPI 00144 ElIsoN.clear(); // Clear the body vector for N of Isotopes 00145 }
| void G4QInelastic::BuildPhysicsTable | ( | const G4ParticleDefinition & | ) | [inline, virtual] |
| G4LorentzVector G4QInelastic::GetEnegryMomentumConservation | ( | ) |
| G4double G4QInelastic::GetMeanFreePath | ( | const G4Track & | aTrack, | |
| G4double | previousStepSize, | |||
| G4ForceCondition * | condition | |||
| ) | [virtual] |
Implements G4VDiscreteProcess.
Definition at line 161 of file G4QInelastic.cc.
References G4AntiKaonZero::AntiKaonZero(), G4AntiLambda::AntiLambda(), G4AntiNeutrinoE::AntiNeutrinoE(), G4AntiNeutrinoMu::AntiNeutrinoMu(), G4AntiNeutrinoTau::AntiNeutrinoTau(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), DBL_MAX, G4Electron::Electron(), G4cerr, G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4QIsotope::Get(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4VQCrossSection::GetCrossSection(), G4QIsotope::GetCSVector(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4Track::GetMaterial(), G4QIsotope::GetMeanCrossSection(), G4Material::GetNumberOfElements(), G4ParticleDefinition::GetPDGEncoding(), G4QANuENuclearCrossSection::GetPointer(), G4QNuENuclearCrossSection::GetPointer(), G4QANuANuNuclearCrossSection::GetPointer(), G4QANuMuNuclearCrossSection::GetPointer(), G4QNuNuNuclearCrossSection::GetPointer(), G4QNuMuNuclearCrossSection::GetPointer(), G4QTauNuclearCrossSection::GetPointer(), G4QElectronNuclearCrossSection::GetPointer(), G4QPhotonNuclearCrossSection::GetPointer(), G4QAntiBaryonPlusNuclearCrossSection::GetPointer(), G4QAntiBaryonNuclearCrossSection::GetPointer(), G4QMuonNuclearCrossSection::GetPointer(), G4QHyperonPlusNuclearCrossSection::GetPointer(), G4QHyperonNuclearCrossSection::GetPointer(), G4QKaonZeroNuclearCrossSection::GetPointer(), G4QKaonPlusNuclearCrossSection::GetPointer(), G4QKaonMinusNuclearCrossSection::GetPointer(), G4QPionPlusNuclearCrossSection::GetPointer(), G4QPionMinusNuclearCrossSection::GetPointer(), G4QProtonNuclearCrossSection::GetPointer(), G4QNeutronCaptureRatio::GetPointer(), G4QNeutronNuclearCrossSection::GetPointer(), G4QNeutronCaptureRatio::GetRatio(), G4DynamicParticle::GetTotalMomentum(), G4Material::GetVecNbOfAtomsPerVolume(), G4QIsotope::InitElement(), IsApplicable(), G4QIsotope::IsDefined(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZero::KaonZero(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), G4NeutrinoE::NeutrinoE(), G4NeutrinoMu::NeutrinoMu(), G4NeutrinoTau::NeutrinoTau(), G4Neutron::Neutron(), NotForced, G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4Positron::Positron(), G4Proton::Proton(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4TauMinus::TauMinus(), G4TauPlus::TauPlus(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().
Referenced by PostStepDoIt().
00162 { 00163 #ifdef debug 00164 G4cout<<"G4QInelastic::GetMeanFreePath: Called Fc="<<*Fc<<G4endl; 00165 #endif 00166 *Fc = NotForced; 00167 #ifdef debug 00168 G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDynPart"<<G4endl; 00169 #endif 00170 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle(); 00171 #ifdef debug 00172 G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDef"<<G4endl; 00173 #endif 00174 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition(); 00175 if( !IsApplicable(*incidentParticleDefinition)) 00176 { 00177 G4cout<<"-W-G4QInelastic::GetMeanFreePath called for not implemented particle"<<G4endl; 00178 return DBL_MAX; 00179 } 00180 // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material 00181 G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle 00182 #ifdef debug 00183 G4cout<<"G4QInelastic::GetMeanFreePath: BeforeGetMaterial"<<G4endl; 00184 #endif 00185 const G4Material* material = aTrack.GetMaterial(); // Get the current material 00186 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume(); 00187 const G4ElementVector* theElementVector = material->GetElementVector(); 00188 G4int nE=material->GetNumberOfElements(); 00189 #ifdef debug 00190 G4cout<<"G4QInelastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl; 00191 #endif 00192 //G4bool leptoNuc=false; // By default the reaction is not lepto-nuclear *Growing point* 00193 G4VQCrossSection* CSmanager=0; 00194 G4VQCrossSection* CSmanager2=0; 00195 G4QNeutronCaptureRatio* capMan=0; 00196 G4int pPDG=0; 00197 G4int pZ = incidentParticleDefinition->GetAtomicNumber(); 00198 G4int pA = incidentParticleDefinition->GetAtomicMass(); 00199 if(incidentParticleDefinition == G4Neutron::Neutron()) 00200 { 00201 CSmanager=G4QNeutronNuclearCrossSection::GetPointer(); 00202 capMan=G4QNeutronCaptureRatio::GetPointer(); 00203 #ifdef debug 00204 G4cout<<"G4QInelastic::GetMeanFreePath: CSmanager is defined for the neutron"<<G4endl; 00205 #endif 00206 pPDG=2112; 00207 } 00208 else if(incidentParticleDefinition == G4Proton::Proton()) 00209 { 00210 CSmanager=G4QProtonNuclearCrossSection::GetPointer(); 00211 pPDG=2212; 00212 } 00213 else if(incidentParticleDefinition == G4PionMinus::PionMinus()) 00214 { 00215 CSmanager=G4QPionMinusNuclearCrossSection::GetPointer(); 00216 pPDG=-211; 00217 } 00218 else if(incidentParticleDefinition == G4PionPlus::PionPlus()) 00219 { 00220 CSmanager=G4QPionPlusNuclearCrossSection::GetPointer(); 00221 pPDG=211; 00222 } 00223 else if(incidentParticleDefinition == G4KaonMinus::KaonMinus()) 00224 { 00225 CSmanager=G4QKaonMinusNuclearCrossSection::GetPointer(); 00226 pPDG=-321; 00227 } 00228 else if(incidentParticleDefinition == G4KaonPlus::KaonPlus()) 00229 { 00230 CSmanager=G4QKaonPlusNuclearCrossSection::GetPointer(); 00231 pPDG=321; 00232 } 00233 else if(incidentParticleDefinition == G4KaonZeroLong::KaonZeroLong() || 00234 incidentParticleDefinition == G4KaonZeroShort::KaonZeroShort() || 00235 incidentParticleDefinition == G4KaonZero::KaonZero() || 00236 incidentParticleDefinition == G4AntiKaonZero::AntiKaonZero() ) 00237 { 00238 CSmanager=G4QKaonZeroNuclearCrossSection::GetPointer(); 00239 if(G4UniformRand() > 0.5) pPDG= 311; 00240 else pPDG=-311; 00241 } 00242 else if(incidentParticleDefinition == G4Lambda::Lambda()) 00243 { 00244 CSmanager=G4QHyperonNuclearCrossSection::GetPointer(); 00245 pPDG=3122; 00246 } 00247 else if(pZ > 0 && pA > 1) // Ions (not implemented yet, should not be used) 00248 { 00249 G4cout<<"-Warning-G4QInelastic::GetMeanFreePath: G4QInelastic called for ions"<<G4endl; 00250 CSmanager=G4QProtonNuclearCrossSection::GetPointer(); 00251 pPDG=90000000+999*pZ+pA; 00252 } 00253 else if(incidentParticleDefinition == G4SigmaPlus::SigmaPlus()) 00254 { 00255 CSmanager=G4QHyperonPlusNuclearCrossSection::GetPointer(); 00256 pPDG=3222; 00257 } 00258 else if(incidentParticleDefinition == G4SigmaMinus::SigmaMinus()) 00259 { 00260 CSmanager=G4QHyperonNuclearCrossSection::GetPointer(); 00261 pPDG=3112; 00262 } 00263 else if(incidentParticleDefinition == G4SigmaZero::SigmaZero()) 00264 { 00265 CSmanager=G4QHyperonNuclearCrossSection::GetPointer(); 00266 pPDG=3212; 00267 } 00268 else if(incidentParticleDefinition == G4XiMinus::XiMinus()) 00269 { 00270 CSmanager=G4QHyperonNuclearCrossSection::GetPointer(); 00271 pPDG=3312; 00272 } 00273 else if(incidentParticleDefinition == G4XiZero::XiZero()) 00274 { 00275 CSmanager=G4QHyperonNuclearCrossSection::GetPointer(); 00276 pPDG=3322; 00277 } 00278 else if(incidentParticleDefinition == G4OmegaMinus::OmegaMinus()) 00279 { 00280 CSmanager=G4QHyperonNuclearCrossSection::GetPointer(); 00281 pPDG=3334; 00282 } 00283 else if(incidentParticleDefinition == G4MuonPlus::MuonPlus() || 00284 incidentParticleDefinition == G4MuonMinus::MuonMinus()) 00285 { 00286 CSmanager=G4QMuonNuclearCrossSection::GetPointer(); 00287 //leptoNuc=true; 00288 pPDG=13; 00289 } 00290 else if(incidentParticleDefinition == G4AntiNeutron::AntiNeutron()) 00291 { 00292 CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer(); 00293 pPDG=-2112; 00294 } 00295 else if(incidentParticleDefinition == G4AntiProton::AntiProton()) 00296 { 00297 CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer(); 00298 pPDG=-2212; 00299 } 00300 else if(incidentParticleDefinition == G4AntiLambda::AntiLambda()) 00301 { 00302 CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer(); 00303 pPDG=-3122; 00304 } 00305 else if(incidentParticleDefinition == G4AntiSigmaPlus::AntiSigmaPlus()) 00306 { 00307 CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer(); 00308 pPDG=-3222; 00309 } 00310 else if(incidentParticleDefinition == G4AntiSigmaMinus::AntiSigmaMinus()) 00311 { 00312 CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer(); 00313 pPDG=-3112; 00314 } 00315 else if(incidentParticleDefinition == G4AntiSigmaZero::AntiSigmaZero()) 00316 { 00317 CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer(); 00318 pPDG=-3212; 00319 } 00320 else if(incidentParticleDefinition == G4AntiXiMinus::AntiXiMinus()) 00321 { 00322 CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer(); 00323 pPDG=-3312; 00324 } 00325 else if(incidentParticleDefinition == G4AntiXiZero::AntiXiZero()) 00326 { 00327 CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer(); 00328 pPDG=-3322; 00329 } 00330 else if(incidentParticleDefinition == G4AntiOmegaMinus::AntiOmegaMinus()) 00331 { 00332 CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer(); 00333 pPDG=-3334; 00334 } 00335 else if(incidentParticleDefinition == G4Gamma::Gamma()) 00336 { 00337 CSmanager=G4QPhotonNuclearCrossSection::GetPointer(); 00338 pPDG=22; 00339 } 00340 else if(incidentParticleDefinition == G4Electron::Electron() || 00341 incidentParticleDefinition == G4Positron::Positron()) 00342 { 00343 CSmanager=G4QElectronNuclearCrossSection::GetPointer(); 00344 //leptoNuc=true; 00345 pPDG=11; 00346 } 00347 else if(incidentParticleDefinition == G4TauPlus::TauPlus() || 00348 incidentParticleDefinition == G4TauMinus::TauMinus()) 00349 { 00350 CSmanager=G4QTauNuclearCrossSection::GetPointer(); 00351 //leptoNuc=true; 00352 pPDG=15; 00353 } 00354 else if(incidentParticleDefinition == G4NeutrinoMu::NeutrinoMu() ) 00355 { 00356 CSmanager=G4QNuMuNuclearCrossSection::GetPointer(); 00357 CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); 00358 //leptoNuc=true; 00359 pPDG=14; 00360 } 00361 else if(incidentParticleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMu() ) 00362 { 00363 CSmanager=G4QANuMuNuclearCrossSection::GetPointer(); 00364 CSmanager2=G4QANuANuNuclearCrossSection::GetPointer(); 00365 //leptoNuc=true; 00366 pPDG=-14; 00367 } 00368 else if(incidentParticleDefinition == G4NeutrinoE::NeutrinoE() ) 00369 { 00370 CSmanager=G4QNuENuclearCrossSection::GetPointer(); 00371 CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); 00372 //leptoNuc=true; 00373 pPDG=12; 00374 } 00375 else if(incidentParticleDefinition == G4AntiNeutrinoE::AntiNeutrinoE() ) 00376 { 00377 CSmanager=G4QANuENuclearCrossSection::GetPointer(); 00378 CSmanager2=G4QANuANuNuclearCrossSection::GetPointer(); 00379 //leptoNuc=true; 00380 pPDG=-12; 00381 } 00382 else 00383 { 00384 G4cout<<"-Warning-G4QInelastic::GetMeanFreePath:Particle " 00385 <<incidentParticleDefinition->GetPDGEncoding()<<" isn't supported"<<G4endl; 00386 return DBL_MAX; 00387 } 00388 G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton 00389 G4double sigma=0.; // Sums over elements for the material 00390 G4int IPIE=IsoProbInEl.size(); // How many old elements? 00391 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) 00392 { 00393 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector 00394 SPI->clear(); 00395 delete SPI; 00396 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector 00397 IsN->clear(); 00398 delete IsN; 00399 } 00400 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) 00401 ElementZ.clear(); // Clear the body vector for Z of Elements 00402 IsoProbInEl.clear(); // Clear the body vector for SPI 00403 ElIsoN.clear(); // Clear the body vector for N of Isotopes 00404 for(G4int i=0; i<nE; ++i) 00405 { 00406 G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element 00407 G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element 00408 ElementZ.push_back(Z); // Remember Z of the Element 00409 G4int isoSize=0; // The default for the isoVectorLength is 0 00410 G4int indEl=0; // Index of non-trivial element or 0(default) 00411 G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect 00412 if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector 00413 #ifdef debug 00414 G4cout<<"G4QInelastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result 00415 #endif 00416 if(isoSize) // The Element has non-trivial abundance set 00417 { 00418 indEl=pElement->GetIndex()+1; // Index of the non-trivial element 00419 if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define 00420 { 00421 std::vector<std::pair<G4int,G4double>*>* newAbund = 00422 new std::vector<std::pair<G4int,G4double>*>; 00423 G4double* abuVector=pElement->GetRelativeAbundanceVector(); 00424 for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes 00425 { 00426 G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z ! 00427 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QInelastic::GetMeanFreePath" 00428 <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl; 00429 G4double abund=abuVector[j]; 00430 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund); 00431 #ifdef debug 00432 G4cout<<"G4QInelastic::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl; 00433 #endif 00434 newAbund->push_back(pr); 00435 } 00436 #ifdef debug 00437 G4cout<<"G4QInelastic::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl; 00438 #endif 00439 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd 00440 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary 00441 delete newAbund; // Was "new" in the beginning of the name space 00442 } 00443 } 00444 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer 00445 std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector 00446 IsoProbInEl.push_back(SPI); 00447 std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector 00448 ElIsoN.push_back(IsN); 00449 G4int nIs=cs->size(); // A#Of Isotopes in the Element 00450 G4double susi=0.; // sum of CS over isotopes 00451 #ifdef debug 00452 G4cout<<"G4QInelastic::GetMeanFreePath: Before Loop nIs="<<nIs<<G4endl; 00453 #endif 00454 if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El 00455 { 00456 std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice 00457 G4int N=curIs->first; // #of Neuterons in the isotope j of El i 00458 IsN->push_back(N); // Remember Min N for the Element 00459 #ifdef debug 00460 G4cout<<"G4QInel::GetMeanFrP: Before CS, P="<<Momentum<<",Z="<<Z<<",N="<<N<<G4endl; 00461 #endif 00462 if(!pPDG) G4cout<<"-Warning-G4QInelastic::GetMeanFrP: (1) projectile PDG=0"<<G4endl; 00463 G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope 00464 if(CSmanager2)CSI+=CSmanager2->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i)nu,nu 00465 if(capMan) CSI*=(1.-capMan->GetRatio(Momentum, Z, N)); 00466 #ifdef debug 00467 G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl; 00468 #endif 00469 curIs->second = CSI; 00470 susi+=CSI; // Make a sum per isotopes 00471 SPI->push_back(susi); // Remember summed cross-section 00472 } // End of temporary initialization of the cross sections in the G4QIsotope singeltone 00473 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV) 00474 ElProbInMat.push_back(sigma); 00475 } // End of LOOP over Elements 00476 #ifdef debug 00477 G4cout<<"G4QInel::GetMeanFrPa:S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl; 00478 #endif 00479 // Check that cross section is not zero and return the mean free path 00480 if(photNucBias!=1.) if(incidentParticleDefinition == G4Gamma::Gamma() || 00481 incidentParticleDefinition == G4MuonPlus::MuonPlus() || 00482 incidentParticleDefinition == G4MuonMinus::MuonMinus() || 00483 incidentParticleDefinition == G4Electron::Electron() || 00484 incidentParticleDefinition == G4Positron::Positron() || 00485 incidentParticleDefinition == G4TauMinus::TauMinus() || 00486 incidentParticleDefinition == G4TauPlus::TauPlus() ) 00487 sigma*=photNucBias; 00488 if(weakNucBias!=1.) if(incidentParticleDefinition==G4NeutrinoE::NeutrinoE() || 00489 incidentParticleDefinition==G4AntiNeutrinoE::AntiNeutrinoE() || 00490 incidentParticleDefinition==G4NeutrinoTau::NeutrinoTau() || 00491 incidentParticleDefinition==G4AntiNeutrinoTau::AntiNeutrinoTau()|| 00492 incidentParticleDefinition==G4NeutrinoMu::NeutrinoMu() || 00493 incidentParticleDefinition==G4AntiNeutrinoMu::AntiNeutrinoMu() ) 00494 sigma*=weakNucBias; 00495 if(sigma > 0.) return 1./sigma; // Mean path [distance] 00496 return DBL_MAX; 00497 }
| G4int G4QInelastic::GetNumberOfNeutronsInTarget | ( | ) |
| G4double G4QInelastic::GetPhotNucBias | ( | ) | [inline] |
| G4double G4QInelastic::GetWeakNucBias | ( | ) | [inline] |
| G4bool G4QInelastic::IsApplicable | ( | const G4ParticleDefinition & | particle | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 500 of file G4QInelastic.cc.
References G4AntiLambda::AntiLambda(), G4AntiNeutrinoE::AntiNeutrinoE(), G4AntiNeutrinoMu::AntiNeutrinoMu(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4Electron::Electron(), G4cout, G4endl, G4Gamma::Gamma(), G4GenericIon::GenericIon(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4ParticleDefinition::GetPDGEncoding(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), G4NeutrinoE::NeutrinoE(), G4NeutrinoMu::NeutrinoMu(), G4Neutron::Neutron(), G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4Positron::Positron(), G4Proton::Proton(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4TauMinus::TauMinus(), G4TauPlus::TauPlus(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().
Referenced by GetMeanFreePath(), and PostStepDoIt().
00501 { 00502 G4int Z=particle.GetAtomicNumber(); 00503 G4int A=particle.GetAtomicMass(); 00504 if (particle == *( G4Neutron::Neutron() )) return true; 00505 else if (particle == *( G4Proton::Proton() )) return true; 00506 else if (particle == *( G4MuonPlus::MuonPlus() )) return true; 00507 else if (particle == *( G4MuonMinus::MuonMinus() )) return true; 00508 else if (particle == *( G4Gamma::Gamma() )) return true; 00509 else if (particle == *( G4Electron::Electron() )) return true; 00510 else if (particle == *( G4Positron::Positron() )) return true; 00511 else if (particle == *( G4PionMinus::PionMinus() )) return true; 00512 else if (particle == *( G4PionPlus::PionPlus() )) return true; 00513 else if (particle == *( G4KaonPlus::KaonPlus() )) return true; 00514 else if (particle == *( G4KaonMinus::KaonMinus() )) return true; 00515 else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true; 00516 else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true; 00517 else if (particle == *( G4Lambda::Lambda() )) return true; 00518 else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true; 00519 else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true; 00520 else if (particle == *( G4SigmaZero::SigmaZero() )) return true; 00521 else if (particle == *( G4XiMinus::XiMinus() )) return true; 00522 else if (particle == *( G4XiZero::XiZero() )) return true; 00523 else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true; 00524 else if (particle == *(G4GenericIon::GenericIon()) || (Z > 0 && A > 1)) return true; 00525 else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true; 00526 else if (particle == *( G4AntiProton::AntiProton() )) return true; 00527 else if (particle == *( G4AntiLambda::AntiLambda() )) return true; 00528 else if (particle == *( G4AntiSigmaPlus::AntiSigmaPlus() )) return true; 00529 else if (particle == *(G4AntiSigmaMinus::AntiSigmaMinus())) return true; 00530 else if (particle == *( G4AntiSigmaZero::AntiSigmaZero() )) return true; 00531 else if (particle == *( G4AntiXiMinus::AntiXiMinus() )) return true; 00532 else if (particle == *( G4AntiXiZero::AntiXiZero() )) return true; 00533 else if (particle == *(G4AntiOmegaMinus::AntiOmegaMinus())) return true; 00534 else if (particle == *( G4TauPlus::TauPlus() )) return true; 00535 else if (particle == *( G4TauMinus::TauMinus() )) return true; 00536 else if (particle == *( G4AntiNeutrinoE::AntiNeutrinoE() )) return true; 00537 else if (particle == *( G4NeutrinoE::NeutrinoE() )) return true; 00538 else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true; 00539 else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true; 00540 //else if (particle == *(G4AntiNeutrinoTau::AntiNeutrinoTau())) return true; 00541 //else if (particle == *( G4NeutrinoTau::NeutrinoTau() )) return true; 00542 #ifdef debug 00543 G4cout<<"***G4QInelastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl; 00544 #endif 00545 return false; 00546 }
| G4VParticleChange * G4QInelastic::PostStepDoIt | ( | const G4Track & | aTrack, | |
| const G4Step & | aStep | |||
| ) | [virtual] |
Reimplemented from G4VDiscreteProcess.
Definition at line 548 of file G4QInelastic.cc.
References G4ParticleChange::AddSecondary(), G4Alpha::Alpha(), G4AntiKaonZero::AntiKaonZero(), G4AntiLambda::AntiLambda(), G4AntiNeutrinoE::AntiNeutrinoE(), G4AntiNeutrinoMu::AntiNeutrinoMu(), G4AntiNeutrinoTau::AntiNeutrinoTau(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4VProcess::aParticleChange, G4QuasiFreeRatios::ChExElCoef(), G4Deuteron::Deuteron(), G4Electron::Electron(), fAlive, FatalException, G4ParticleTable::FindIon(), G4QFragmentation::Fragment(), G4QIonIonCollision::Fragment(), G4Quasmon::Fragment(), fStopAndKill, fStopButAlive, G4cerr, G4cout, G4endl, G4Exception(), G4QThd(), G4RandomDirection(), G4UniformRand, G4Gamma::Gamma(), G4QPDGToG4Particle::Get(), G4QCHIPSWorld::Get(), G4DynamicParticle::Get4Momentum(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4VQCrossSection::GetCrossSection(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Material::GetElementVector(), G4ParticleChange::GetEnergy(), G4VQCrossSection::GetExchangeEnergy(), G4VQCrossSection::GetExchangePDGCode(), G4VQCrossSection::GetExchangeQ2(), G4Track::GetGlobalTime(), G4DynamicParticle::GetKineticEnergy(), G4VQCrossSection::GetLastQELCS(), G4VQCrossSection::GetLastTOTCS(), G4DynamicParticle::GetMass(), G4QPDGCode::GetMass(), G4Track::GetMaterial(), GetMeanFreePath(), G4ParticleChange::GetMomentumDirection(), G4DynamicParticle::GetMomentumDirection(), G4VQCrossSection::GetNQE_ExchangeQ2(), G4QPDGCode::GetNuclMass(), G4Material::GetNumberOfElements(), G4QPDGToG4Particle::GetParticleDefinition(), G4ParticleDefinition::GetParticleName(), G4QCHIPSWorld::GetParticles(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4QuasiFreeRatios::GetPointer(), G4QANuENuclearCrossSection::GetPointer(), G4QNuENuclearCrossSection::GetPointer(), G4QANuANuNuclearCrossSection::GetPointer(), G4QANuMuNuclearCrossSection::GetPointer(), G4QNuNuNuclearCrossSection::GetPointer(), G4QNuMuNuclearCrossSection::GetPointer(), G4QPhotonNuclearCrossSection::GetPointer(), G4QTauNuclearCrossSection::GetPointer(), G4QMuonNuclearCrossSection::GetPointer(), G4QElectronNuclearCrossSection::GetPointer(), G4Track::GetPosition(), G4VQCrossSection::GetQEL_ExchangeQ2(), G4QHadron::GetQPDG(), G4QuasiFreeRatios::GetRatios(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetTouchableHandle(), G4VParticleChange::GetTrackStatus(), G4VQCrossSection::GetVirtualFactor(), G4Track::GetWeight(), G4Element::GetZ(), G4He3::He3(), G4ParticleChange::Initialize(), IsApplicable(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZero::KaonZero(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), G4NeutrinoE::NeutrinoE(), G4NeutrinoMu::NeutrinoMu(), G4NeutrinoTau::NeutrinoTau(), G4Neutron::Neutron(), NotForced, G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), position, G4Positron::Positron(), G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4Proton::Proton(), G4QuasiFreeRatios::Scatter(), G4DynamicParticle::Set4Momentum(), G4DynamicParticle::SetDefinition(), G4VParticleChange::SetNumberOfSecondaries(), G4QNucleus::SetParameters(), G4Track::SetTouchableHandle(), G4Track::SetWeight(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4TauMinus::TauMinus(), G4TauPlus::TauPlus(), G4Triton::Triton(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().
00549 { 00550 static const G4double third = 1./3.; 00551 static const G4double me=G4Electron::Electron()->GetPDGMass(); // electron mass 00552 static const G4double me2=me*me; // squared electron mass 00553 static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass 00554 static const G4double mu2=mu*mu; // squared muon mass 00555 static const G4double mt=G4TauMinus::TauMinus()->GetPDGMass(); // tau mass 00556 static const G4double mt2=mt*mt; // squared tau mass 00557 //static const G4double dpi=M_PI+M_PI; // 2*pi (for Phi distr.) ***changed to twopi*** 00558 static const G4double mNeut= G4QPDGCode(2112).GetMass(); 00559 static const G4double mNeut2= mNeut*mNeut; 00560 static const G4double mProt= G4QPDGCode(2212).GetMass(); 00561 static const G4double mProt2= mProt*mProt; 00562 static const G4double dM=mProt+mNeut; // doubled nucleon mass 00563 static const G4double hdM=dM/2.; // M of the "nucleon" 00564 static const G4double hdM2=hdM*hdM; // M2 of the "nucleon" 00565 static const G4double mLamb= G4QPDGCode(3122).GetMass(); // Mass of Lambda/antiL 00566 static const G4double mPi0 = G4QPDGCode(111).GetMass(); 00567 static const G4double mPi0s= mPi0*mPi0; 00568 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron 00569 //static const G4double mDeut2=mDeut*mDeut; // SquaredMassOfDeuteron 00570 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium 00571 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3 00572 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha 00573 static const G4double mPi = G4QPDGCode(211).GetMass(); 00574 static const G4double tmPi = mPi+mPi; // Doubled mass of the charged pion 00575 static const G4double stmPi= tmPi*tmPi; // Squared Doubled mass of the charged pion 00576 static const G4double mPPi = mPi+mProt; // Delta threshold 00577 static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2 00578 //static const G4double mDel2= 1400*1400; // Delta up threshold for W2 (in MeV^2) 00579 // Static definitions for electrons (nu,e) ----------------------------------------- 00580 static const G4double meN = mNeut+me; 00581 static const G4double meN2= meN*meN; 00582 static const G4double fmeN= 4*mNeut2*me2; 00583 static const G4double mesN= mNeut2+me2; 00584 static const G4double meP = mProt+me; 00585 static const G4double meP2= meP*meP; 00586 static const G4double fmeP= 4*mProt2*me2; 00587 static const G4double mesP= mProt2+me2; 00588 static const G4double medM= me2/dM; // for x limit 00589 static const G4double meD = mPPi+me; // Multiperipheral threshold 00590 static const G4double meD2= meD*meD; 00591 // Static definitions for muons (nu,mu) ----------------------------------------- 00592 static const G4double muN = mNeut+mu; 00593 static const G4double muN2= muN*muN; 00594 static const G4double fmuN= 4*mNeut2*mu2; 00595 static const G4double musN= mNeut2+mu2; 00596 static const G4double muP = mProt+mu; 00597 static const G4double muP2= muP*muP; // + 00598 static const G4double fmuP= 4*mProt2*mu2; // + 00599 static const G4double musP= mProt2+mu2; 00600 static const G4double mudM= mu2/dM; // for x limit 00601 static const G4double muD = mPPi+mu; // Multiperipheral threshold 00602 static const G4double muD2= muD*muD; 00603 // Static definitions for muons (nu,nu) ----------------------------------------- 00604 //static const G4double nuN = mNeut; 00605 //static const G4double nuN2= mNeut2; 00606 //static const G4double fnuN= 0.; 00607 //static const G4double nusN= mNeut2; 00608 //static const G4double nuP = mProt; 00609 //static const G4double nuP2= mProt2; 00610 //static const G4double fnuP= 0.; 00611 //static const G4double nusP= mProt2; 00612 //static const G4double nudM= 0.; // for x limit 00613 //static const G4double nuD = mPPi; // Multiperipheral threshold 00614 //static const G4double nuD2= mPPi2; 00615 static const G4LorentzVector vacuum4M(0.,0.,0.,0.); 00616 //------------------------------------------------------------------------------------- 00617 static G4bool CWinit = true; // CHIPS Warld needs to be initted 00618 if(CWinit) 00619 { 00620 CWinit=false; 00621 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) 00622 } 00623 //------------------------------------------------------------------------------------- 00624 const G4DynamicParticle* projHadron = track.GetDynamicParticle(); 00625 const G4ParticleDefinition* particle=projHadron->GetDefinition(); 00626 #ifdef debug 00627 G4cout<<"G4QInelastic::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl; 00628 #endif 00629 G4ForceCondition cond=NotForced; 00630 GetMeanFreePath(track, 1., &cond); 00631 #ifdef debug 00632 G4cout<<"G4QInelastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl; 00633 #endif 00634 G4bool scat=false; // No CHEX in proj scattering 00635 G4int scatPDG=0; // Must be filled if true (CHEX) 00636 G4LorentzVector proj4M=projHadron->Get4Momentum(); // 4-momentum of the projectile (IU?) 00637 G4LorentzVector scat4M=proj4M; // Must be filled if true 00638 G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle 00639 G4double Momentum=proj4M.rho(); 00640 if(std::fabs(Momentum-momentum)>.001) 00641 G4cerr<<"*G4QInelastic::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl; 00642 #ifdef debug 00643 G4double mp=proj4M.m(); 00644 G4cout<<"-->G4QInel::PostStDoIt:*called*, 4M="<<proj4M<<", P="<<Momentum<<"="<<momentum 00645 <<",m="<<mp<<G4endl; 00646 #endif 00647 if (!IsApplicable(*particle)) // Check applicability 00648 { 00649 G4cerr<<"G4QInelastic::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented." 00650 <<G4endl; 00651 return 0; 00652 } 00653 const G4Material* material = track.GetMaterial(); // Get the current material 00654 G4int Z=0; 00655 const G4ElementVector* theElementVector = material->GetElementVector(); 00656 G4int nE=material->GetNumberOfElements(); 00657 #ifdef debug 00658 G4cout<<"G4QInelastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl; 00659 #endif 00660 G4int projPDG=0; // PDG Code prototype for the captured hadron 00661 G4int pZ=particle->GetAtomicNumber(); 00662 G4int pA=particle->GetAtomicMass(); 00663 if (particle == G4Neutron::Neutron() ) projPDG= 2112; 00664 else if (particle == G4Proton::Proton() ) projPDG= 2212; 00665 else if (particle == G4PionMinus::PionMinus() ) projPDG= -211; 00666 else if (particle == G4PionPlus::PionPlus() ) projPDG= 211; 00667 else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321; 00668 else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321; 00669 else if (particle == G4KaonZeroLong::KaonZeroLong() || 00670 particle == G4KaonZeroShort::KaonZeroShort() || 00671 particle == G4KaonZero::KaonZero() || 00672 particle == G4AntiKaonZero::AntiKaonZero() ) 00673 { 00674 if(G4UniformRand() > 0.5) projPDG= 311; 00675 else projPDG= -311; 00676 } 00677 else if ( pZ > 0 && pA > 1 ) projPDG = 90000000+999*pZ+pA; 00678 else if (particle == G4Lambda::Lambda() ) projPDG= 3122; 00679 else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222; 00680 else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112; 00681 else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212; 00682 else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312; 00683 else if (particle == G4XiZero::XiZero() ) projPDG= 3322; 00684 else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334; 00685 else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112; 00686 else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212; 00687 else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13; 00688 else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13; 00689 else if (particle == G4Gamma::Gamma() ) projPDG= 22; 00690 else if (particle == G4Electron::Electron() ) projPDG= 11; 00691 else if (particle == G4Positron::Positron() ) projPDG= -11; 00692 else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14; 00693 else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14; 00694 else if (particle == G4AntiLambda::AntiLambda() ) projPDG=-3122; 00695 else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) projPDG=-3222; 00696 else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112; 00697 else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) projPDG=-3212; 00698 else if (particle == G4AntiXiMinus::AntiXiMinus() ) projPDG=-3312; 00699 else if (particle == G4AntiXiZero::AntiXiZero() ) projPDG=-3322; 00700 else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334; 00701 else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12; 00702 else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12; 00703 else if (particle == G4TauPlus::TauPlus() ) projPDG= -15; 00704 else if (particle == G4TauMinus::TauMinus() ) projPDG= 15; 00705 else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16; 00706 else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16; 00707 G4int aProjPDG=std::abs(projPDG); 00708 #ifdef debug 00709 G4int prPDG=particle->GetPDGEncoding(); 00710 G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl; 00711 #endif 00712 if(!projPDG) 00713 { 00714 G4cerr<<"---Warning---G4QInelastic::PostStepDoIt:Undefined interacting hadron"<<G4endl; 00715 return 0; 00716 } 00717 G4int EPIM=ElProbInMat.size(); 00718 #ifdef debug 00719 G4cout<<"G4QInel::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl; 00720 #endif 00721 G4int i=0; 00722 if(EPIM>1) 00723 { 00724 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand(); 00725 for(i=0; i<nE; ++i) 00726 { 00727 #ifdef debug 00728 G4cout<<"G4QInelastic::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl; 00729 #endif 00730 if (rnd<ElProbInMat[i]) break; 00731 } 00732 if(i>=nE) i=nE-1; // Top limit for the Element 00733 } 00734 G4Element* pElement=(*theElementVector)[i]; 00735 Z=static_cast<G4int>(pElement->GetZ()); 00736 #ifdef debug 00737 G4cout<<"G4QInelastic::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl; 00738 #endif 00739 if(Z<=0) 00740 { 00741 G4cerr<<"---Warning---G4QInelastic::PostStepDoIt: Element with Z="<<Z<<G4endl; 00742 if(Z<0) return 0; 00743 } 00744 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes 00745 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i] 00746 G4int nofIsot=SPI->size(); // #of isotopes in the element i 00747 #ifdef debug 00748 G4cout<<"G4QInel::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl; 00749 #endif 00750 G4int j=0; 00751 if(nofIsot>1) 00752 { 00753 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element 00754 for(j=0; j<nofIsot; ++j) 00755 { 00756 #ifdef debug 00757 G4cout<<"G4QInelastic::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl; 00758 #endif 00759 if(rndI < (*SPI)[j]) break; 00760 } 00761 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope 00762 } 00763 G4int N =(*IsN)[j]; ; // Randomized number of neutrons 00764 #ifdef debug 00765 G4cout<<"G4QInelastic::PostStepDoIt: Z="<<Z<<", j="<<i<<", N(isotope)="<<N<<G4endl; 00766 #endif 00767 G4double kinEnergy= projHadron->GetKineticEnergy(); 00768 G4ParticleMomentum dir = projHadron->GetMomentumDirection(); 00769 if(projPDG==22 && Z==1 && !N && Momentum<150.) 00770 { 00771 #ifdef debug 00772 G4cerr<<"---LowPHOTON---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl; 00773 #endif 00774 //Do Nothing Action insead of the reaction 00775 aParticleChange.ProposeEnergy(kinEnergy); 00776 aParticleChange.ProposeLocalEnergyDeposit(0.); 00777 aParticleChange.ProposeMomentumDirection(dir); 00778 aParticleChange.ProposeTrackStatus(fAlive); 00779 return G4VDiscreteProcess::PostStepDoIt(track,step); 00780 } 00781 if(N<0) 00782 { 00783 G4cerr<<"-Warning-G4QInelastic::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl; 00784 return 0; 00785 } 00786 nOfNeutrons=N; // Remember it for the energy-momentum check 00787 //G4double am=Z+N; 00788 G4int am=Z+N; 00789 //G4double dd=0.025; 00790 //G4double sr=std::sqrt(am); 00791 //G4double dsr=0.01*(sr+sr); 00792 //if(dsr<dd)dsr=dd; 00793 //G4double medRA=mediRatio*G4QThd(am); 00794 G4double medRA=mediRatio; 00795 if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,medRA); // ManualP 00796 //else if(projPDG==-2212)G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,9.);//aP ClustPars 00797 //else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.); //Pi-ClustPars 00798 //else if(projPDG ==-2212 || projPDG ==-2112) 00799 // G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,1.);//aP ClustPars 00800 //else if(projPDG ==-211 ||projPDG == 211 ) 00801 // G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,1.); //Pi-ClustPars 00802 #ifdef debug 00803 G4cout<<"G4QInelastic::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl; 00804 #endif 00805 aParticleChange.Initialize(track); 00806 G4double weight = track.GetWeight(); 00807 #ifdef debug 00808 G4cout<<"G4QInelastic::PostStepDoIt: weight="<<weight<<G4endl; 00809 #endif 00810 if(photNucBias!=1.) weight/=photNucBias; 00811 else if(weakNucBias!=1.) weight/=weakNucBias; 00812 G4double localtime = track.GetGlobalTime(); 00813 #ifdef debug 00814 G4cout<<"G4QInelastic::PostStepDoIt: localtime="<<localtime<<G4endl; 00815 #endif 00816 G4ThreeVector position = track.GetPosition(); 00817 G4TouchableHandle trTouchable = track.GetTouchableHandle(); 00818 #ifdef debug 00819 G4cout<<"G4QInelastic::PostStepDoIt: position="<<position<<G4endl; 00820 #endif 00821 // 00822 G4int targPDG=90000000+Z*1000+N; // PDG Code of the target nucleus 00823 G4QPDGCode targQPDG(targPDG); 00824 G4double tgM=targQPDG.GetMass(); // Target mass 00825 G4double tM=tgM; // Target mass (copy to be changed) 00826 G4QHadronVector* output=new G4QHadronVector;// Prototype of Fragm Output G4QHadronVector 00827 G4double absMom = 0.; // Prototype of absorbed by nucleus Moment 00828 G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVector 00829 G4LorentzVector lead4M(0.,0.,0.,0.); // Prototype of LeadingQ 4-momentum 00830 #ifdef debug 00831 G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<aProjPDG<<", targPDG="<<targPDG<<G4endl; 00832 #endif 00833 // 00834 // Leptons with photonuclear 00835 // Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?) 00836 // 00837 if (aProjPDG == 11 || aProjPDG == 13 || aProjPDG == 15) 00838 { 00839 #ifdef debug 00840 G4cout<<"G4QInelastic::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl; 00841 #endif 00842 G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer(); 00843 G4double ml=me; 00844 G4double ml2=me2; 00845 if(aProjPDG== 13) 00846 { 00847 CSmanager=G4QMuonNuclearCrossSection::GetPointer(); 00848 ml=mu; 00849 ml2=mu2; 00850 } 00851 if(aProjPDG== 15) 00852 { 00853 CSmanager=G4QTauNuclearCrossSection::GetPointer(); 00854 ml=mt; 00855 ml2=mt2; 00856 } 00857 // @@ Probably this is not necessary any more (?) 00858 if(!aProjPDG) G4cout<<"-Warning-G4QInelast::PostStepDoIt:(2) projectile PDG=0"<<G4endl; 00859 G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,aProjPDG);// Recalculate XS 00860 // @@ check a possibility to separate p, n, or alpha (!) 00861 if(Z==1 && !N && Momentum<150.) xSec=0.; 00862 #ifdef debug 00863 G4cout<<"-Forse-G4QInel::PStDoIt:P="<<Momentum<<",PDG="<<projPDG<<",xS="<<xSec<<G4endl; 00864 #endif 00865 if(xSec <= 0.) // The cross-section is 0 -> Do Nothing 00866 { 00867 #ifdef debug 00868 G4cerr<<"---OUT---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl; 00869 #endif 00870 //Do Nothing Action insead of the reaction 00871 aParticleChange.ProposeEnergy(kinEnergy); 00872 aParticleChange.ProposeLocalEnergyDeposit(0.); 00873 aParticleChange.ProposeMomentumDirection(dir); 00874 aParticleChange.ProposeTrackStatus(fAlive); 00875 delete output; 00876 return G4VDiscreteProcess::PostStepDoIt(track,step); 00877 } 00878 G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart 00879 #ifdef debug 00880 G4cout<<"G4QInel::PStDoIt:kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl; 00881 #endif 00882 if( kinEnergy < photonEnergy || photonEnergy < 0.) 00883 { 00884 //Do Nothing Action insead of the reaction 00885 G4cerr<<"--G4QInelastic::PSDoIt: photE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl; 00886 aParticleChange.ProposeEnergy(kinEnergy); 00887 aParticleChange.ProposeLocalEnergyDeposit(0.); 00888 aParticleChange.ProposeMomentumDirection(dir); 00889 aParticleChange.ProposeTrackStatus(fAlive); 00890 delete output; 00891 return G4VDiscreteProcess::PostStepDoIt(track,step); 00892 } 00893 G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart 00894 G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon 00895 if(tM<999.) W-=mPi0+mPi0s/dM; // Pion production threshold for a nucleon target 00896 if(W<0.) 00897 { 00898 //Do Nothing Action insead of the reaction 00899 #ifdef debug 00900 G4cout<<"--G4QInelastic::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl; 00901 #endif 00902 aParticleChange.ProposeEnergy(kinEnergy); 00903 aParticleChange.ProposeLocalEnergyDeposit(0.); 00904 aParticleChange.ProposeMomentumDirection(dir); 00905 aParticleChange.ProposeTrackStatus(fAlive); 00906 delete output; 00907 return G4VDiscreteProcess::PostStepDoIt(track,step); 00908 } 00909 // Update G4VParticleChange for the scattered muon 00910 G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer(); 00911 G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy,Z,N,22);//Integrated XS 00912 G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N, 22); // Real XS 00913 G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2); 00914 if(sigNu*G4UniformRand()>sigK*rndFraction) 00915 { 00916 //NothingToDo Action insead of the reaction 00917 #ifdef debug 00918 G4cout<<"-DoNoth-G4QInelastic::PostStepDoIt: probab. correction - DoNothing"<<G4endl; 00919 #endif 00920 aParticleChange.ProposeEnergy(kinEnergy); 00921 aParticleChange.ProposeLocalEnergyDeposit(0.); 00922 aParticleChange.ProposeMomentumDirection(dir); 00923 aParticleChange.ProposeTrackStatus(fAlive); 00924 delete output; 00925 return G4VDiscreteProcess::PostStepDoIt(track,step); 00926 } 00927 G4double iniE=kinEnergy+ml; // Initial total energy of the lepton 00928 G4double finE=iniE-photonEnergy; // Final total energy of the lepton 00929 #ifdef pdebug 00930 G4cout<<"G4QInelastic::PoStDoIt:E="<<iniE<<",lE="<<finE<<"-"<<ml<<"="<<finE-ml<<G4endl; 00931 #endif 00932 aParticleChange.ProposeEnergy(finE-ml); 00933 if(finE<=ml) // Secondary lepton (e/mu/tau) at rest disappears 00934 { 00935 aParticleChange.ProposeEnergy(0.); 00936 if(aProjPDG== 11) aParticleChange.ProposeTrackStatus(fStopAndKill); 00937 else aParticleChange.ProposeTrackStatus(fStopButAlive); 00938 aParticleChange.ProposeMomentumDirection(dir); 00939 } 00940 else aParticleChange.ProposeTrackStatus(fAlive); 00941 G4double iniP=std::sqrt(iniE*iniE-ml2); // Initial momentum of the electron 00942 G4double finP=std::sqrt(finE*finE-ml2); // Final momentum of the electron 00943 G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton) 00944 #ifdef pdebug 00945 G4cout<<"G4QI::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl; 00946 #endif 00947 if(cost>1.) cost=1.; // To avoid the accuracy of calculation problem 00948 if(cost<-1.) cost=-1.; // To avoid the accuracy of calculation problem 00949 // 00950 // Scatter the lepton ( @@ make the same thing for real photons) 00951 // At this point we have photonEnergy and photonQ2 (with notDefinedPhi)->SelectProjPart 00952 G4double absEn = G4QThd(am)*GeV; // @@(b) Mean Energy Absorbed by a Nucleus 00953 //G4double absEn = std::pow(am,third)*GeV; // @@(b) Mean Energy Absorbed by a Nucleus 00954 if(am>1 && absEn < photonEnergy) // --> the absorption of energy can happen 00955 //if(absEn < photonEnergy) // --> the absorption of energy can happen 00956 { 00957 G4double abtEn = absEn+hdM; // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N) 00958 G4double abEn2 = abtEn*abtEn; // Squared absorbed Energy + MN 00959 G4double abMo2 = abEn2-hdM2; // Squared absorbed Momentum of compound system 00960 G4double phEn2 = photonEnergy*photonEnergy; 00961 G4double phMo2 = phEn2+photonQ2; // Squared momentum of primary virtual photon 00962 G4double phMo = std::sqrt(phMo2); // Momentum of the primary virtual photon 00963 absMom = std::sqrt(abMo2); // Absorbed Momentum 00964 if(absMom < phMo) // --> the absorption of momentum can happen 00965 { 00966 G4double dEn = photonEnergy - absEn; // Leading energy 00967 G4double dMo = phMo - absMom; // Leading momentum 00968 G4double sF = dEn*dEn - dMo*dMo;// s of leading particle 00969 #ifdef ppdebug 00970 G4cout<<"-PhotoAbsorption-G4QIn::PStDoIt: sF="<<sF<<",phEn="<<photonEnergy<<G4endl; 00971 #endif 00972 if(sF > stmPi) // --> Leading fragmentation is possible 00973 { 00974 photonEnergy = absEn; // New value of the photon energy 00975 photonQ2=abMo2-absEn*absEn; // New value of the photon Q2 00976 absEn = dEn; // Put energy of leading particle to absEn (!) 00977 } 00978 else absMom=0.; // Flag that nothing has happened 00979 } 00980 else absMom=0.; // Flag that nothing has happened 00981 } 00982 // ------------- End of ProjPart selection 00983 // 00984 // Scattering in respect to the derection of the incident muon is made impicitly: 00985 G4ThreeVector ort=dir.orthogonal(); // Not normed orthogonal vector (!) (to dir) 00986 G4ThreeVector ortx = ort.unit(); // First unit vector orthogonal to the direction 00987 G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction 00988 G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component 00989 G4double phi=twopi*G4UniformRand(); // phi of scattered electron 00990 G4double sinx=sint*std::sin(phi); // x perpendicular component 00991 G4double siny=sint*std::cos(phi); // y perpendicular component 00992 G4ThreeVector findir=cost*dir+sinx*ortx+siny*orty; 00993 aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton 00994 #ifdef pdebug 00995 G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-" 00996 <<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl; 00997 #endif 00998 G4ThreeVector photon3M=iniP*dir-finP*findir;// 3D total momentum of photon 00999 if(absMom) // Photon must be reduced & LeadingSyst fragmented 01000 { 01001 G4double ptm=photon3M.mag(); // 3M of the virtual photon 01002 #ifdef ppdebug 01003 G4cout<<"-Absorption-G4QInelastic::PostStepDoIt: ph3M="<<photon3M<<", eIn3M=" 01004 <<iniP*dir<<", eFin3M="<<finP*findir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl; 01005 #endif 01006 G4ThreeVector lead3M=photon3M*(ptm-absMom)/ptm; // Keep the direction for leading Q 01007 photon3M-=lead3M; // Reduced photon Momentum (photEn already = absEn) 01008 proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System 01009 #ifdef ppdebug 01010 G4cout<<"-->G4QI::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl; 01011 #endif 01012 lead4M=proj4M; // Remember 4-mom for the total 4-momentum 01013 G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+ 01014 try // | 01015 { // | 01016 if(leadhs) delete leadhs; // | 01017 G4QNucleus vac(90000000); // | 01018 leadhs=pan->Fragment(vac,1); // DELETED after it is copied to output vector | 01019 } // | 01020 catch (G4QException& error) // | 01021 { // | 01022 G4cerr<<"***G4QInelastic::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//| 01023 // G4Exception("G4QInelastic::PostStepDoIt:","72",FatalException,"QuasmonCrash"); //| 01024 G4Exception("G4QInelastic::PostStepDoIt()","HAD_CHPS_0072", 01025 FatalException, "QuasmonCrash"); 01026 } // | 01027 delete pan; // Delete the Nuclear Environment <----<---+ 01028 #ifdef ppdebug 01029 G4cout<<"G4QInel::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt=" 01030 <<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl; 01031 #endif 01032 } 01033 else 01034 { 01035 G4int qNH=0; 01036 if(leadhs) qNH=leadhs->size(); 01037 if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq]; 01038 if(leadhs) delete leadhs; 01039 leadhs=0; 01040 } 01041 projPDG=22; 01042 proj4M=G4LorentzVector(photon3M,photonEnergy); 01043 #ifdef debug 01044 G4cout<<"G4QInelastic::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m=" 01045 <<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl; 01046 #endif 01047 01048 // 01049 // neutrinoNuclear interactions (nu_e, nu_mu only) 01050 // 01051 } 01052 else if (aProjPDG == 12 || aProjPDG == 14) 01053 { 01054 kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino 01055 G4double dKinE=kinEnergy+kinEnergy; // doubled energy for s calculation 01056 #ifdef debug 01057 G4cout<<"G4QInelastic::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl; 01058 #endif 01059 dir = projHadron->GetMomentumDirection(); // unit vector 01060 G4double ml = mu; 01061 G4double ml2 = mu2; 01062 //G4double mlN = muN; 01063 G4double mlN2= muN2; 01064 G4double fmlN= fmuN; 01065 G4double mlsN= musN; 01066 //G4double mlP = muP; 01067 G4double mlP2= muP2; 01068 G4double fmlP= fmuP; 01069 G4double mlsP= musP; 01070 G4double mldM= mudM; 01071 //G4double mlD = muD; 01072 G4double mlD2= muD2; 01073 if(aProjPDG==12) 01074 { 01075 ml = me; 01076 ml2 = me2; 01077 //mlN = meN; 01078 mlN2= meN2; 01079 fmlN= fmeN; 01080 mlsN= mesN; 01081 //mlP = meP; 01082 mlP2= meP2; 01083 fmlP= fmeP; 01084 mlsP= mesP; 01085 mldM= medM; 01086 //mlD = meD; 01087 mlD2= meD2; 01088 } 01089 G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l) 01090 G4VQCrossSection* CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); // (nu,nu) 01091 proj4M=G4LorentzVector(dir*kinEnergy,kinEnergy); // temporary 01092 G4bool nuanu=true; 01093 scatPDG=13; // Prototype = secondary scattered mu- 01094 if(projPDG==-14) 01095 { 01096 nuanu=false; // Anti-neutrino 01097 CSmanager=G4QANuMuNuclearCrossSection::GetPointer(); // (anu,mu+) CC @@ open 01098 CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open 01099 scatPDG=-13; // secondary scattered mu+ 01100 } 01101 else if(projPDG==12) 01102 { 01103 CSmanager=G4QNuENuclearCrossSection::GetPointer(); // @@ open (only CC is changed) 01104 scatPDG=11; // secondary scattered e- 01105 } 01106 else if(projPDG==-12) 01107 { 01108 nuanu=false; // anti-neutrino 01109 CSmanager=G4QANuENuclearCrossSection::GetPointer(); // (anu,e+) CC @@ open 01110 CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open 01111 scatPDG=-11; // secondary scattered e+ 01112 } 01113 // @@ Probably this is not necessary any more 01114 if(!projPDG) G4cout<<"-Warning-G4QInelastic::PostStepDoIt:(3)projectile PDG=0"<<G4endl; 01115 G4double xSec1=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG); //Recalculate XS 01116 G4double xSec2=CSmanager2->GetCrossSection(false,Momentum,Z,N,projPDG);//Recalculate XS 01117 G4double xSec=xSec1+xSec2; 01118 // @@ check a possibility to separate p, n, or alpha (!) 01119 if(xSec <= 0.) // The cross-section = 0 -> Do Nothing 01120 { 01121 G4cerr<<"G4QInelastic::PSDoIt:nuE="<<kinEnergy<<",X1="<<xSec1<<",X2="<<xSec2<<G4endl; 01122 //Do Nothing Action insead of the reaction 01123 aParticleChange.ProposeEnergy(kinEnergy); 01124 aParticleChange.ProposeLocalEnergyDeposit(0.); 01125 aParticleChange.ProposeMomentumDirection(dir); 01126 aParticleChange.ProposeTrackStatus(fAlive); 01127 delete output; 01128 return G4VDiscreteProcess::PostStepDoIt(track,step); 01129 } 01130 G4bool secnu=false; 01131 if(xSec*G4UniformRand()>xSec1) // recover neutrino/antineutrino 01132 { 01133 if(scatPDG>0) scatPDG++; 01134 else scatPDG--; 01135 secnu=true; 01136 } 01137 scat=true; // event with changed scattered projectile 01138 G4double totCS1 = CSmanager->GetLastTOTCS(); // the last total cross section1(isotope?) 01139 G4double totCS2 = CSmanager2->GetLastTOTCS();// the last total cross section2(isotope?) 01140 G4double totCS = totCS1+totCS2; // the last total cross section (isotope?) 01141 if(std::fabs(xSec-totCS*millibarn)/xSec>.0001) 01142 G4cout<<"-Warning-G4QInelastic::PostStepDoIt: xS="<<xSec<<"# CS="<<totCS<<G4endl; 01143 G4double qelCS1 = CSmanager->GetLastQELCS(); // the last quasi-elastic cross section1 01144 G4double qelCS2 = CSmanager2->GetLastQELCS();// the last quasi-elastic cross section2 01145 G4double qelCS = qelCS1+qelCS2; // the last quasi-elastic cross section 01146 if(totCS - qelCS < 0.) // only at low energies 01147 { 01148 totCS = qelCS; 01149 totCS1 = qelCS1; 01150 totCS2 = qelCS2; 01151 } 01152 // make different definitions for neutrino and antineutrino 01153 G4double mIN=mProt; // Just a prototype (for anu, Z=1, N=0) 01154 G4double mOT=mNeut; 01155 G4double OT=mlN2; 01156 //G4double mOT2=mNeut2; // Other formula: sd=s-mlsOT -Not used?- 01157 G4double mlOT=fmlN; 01158 G4double mlsOT=mlsN; 01159 if(secnu) 01160 { 01161 if(am*G4UniformRand()>Z) // Neutron target 01162 { 01163 targPDG-=1; // subtract neutron 01164 projPDG=2112; // neutron is going out 01165 mIN =mNeut; 01166 OT =mNeut2; 01167 //mOT2=mNeut2; 01168 mlOT=0.; 01169 mlsOT=mNeut2; 01170 } 01171 else 01172 { 01173 targPDG-=1000; // subtract neutron 01174 projPDG=2212; // neutron is going out 01175 mOT =mProt; 01176 OT =mProt2; 01177 //mOT2 =mProt2; 01178 mlOT =0.; 01179 mlsOT=mProt2; 01180 } 01181 ml=0.; 01182 ml2=0.; 01183 mldM=0.; 01184 mlD2=mPPi2; 01185 G4QPDGCode temporary_targQPDG(targPDG); 01186 targQPDG = temporary_targQPDG; 01187 G4double rM=targQPDG.GetMass(); 01188 mIN=tM-rM; // bounded in-mass of the neutron 01189 tM=rM; 01190 } 01191 else if(nuanu) 01192 { 01193 targPDG-=1; // Neutrino -> subtract neutron 01194 G4QPDGCode temporary_targQPDG(targPDG); 01195 targQPDG = temporary_targQPDG; 01196 G4double rM=targQPDG.GetMass(); 01197 mIN=tM-rM; // bounded in-mass of the neutron 01198 tM=rM; 01199 mOT=mProt; 01200 OT=mlP2; 01201 //mOT2=mProt2; 01202 mlOT=fmlP; 01203 mlsOT=mlsP; 01204 projPDG=2212; // proton is going out 01205 } 01206 else 01207 { 01208 if(Z>1||N>0) // Calculate the splitted mass 01209 { 01210 targPDG-=1000; // Anti-Neutrino -> subtract proton 01211 G4QPDGCode temporary_targQPDG(targPDG); 01212 targQPDG = temporary_targQPDG; 01213 G4double rM=targQPDG.GetMass(); 01214 mIN=tM-rM; // bounded in-mass of the proton 01215 tM=rM; 01216 } 01217 else 01218 { 01219 targPDG=0; 01220 mIN=tM; 01221 tM=0.; 01222 } 01223 projPDG=2112; // neutron is going out 01224 } 01225 G4double s_value=mIN*(mIN+dKinE); // s_value=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu) 01226 #ifdef debug 01227 G4cout<<"G4QInelastic::PostStDoIt: s="<<s_value<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl; 01228 #endif 01229 if(s_value<=OT) // *** Do nothing solution *** 01230 { 01231 //Do NothingToDo Action insead of the reaction (@@ Can we make it common?) 01232 G4cout<<"G4QInelastic::PostStepDoIt: tooSmallFinalMassOfCompound: DoNothing"<<G4endl; 01233 aParticleChange.ProposeEnergy(kinEnergy); 01234 aParticleChange.ProposeLocalEnergyDeposit(0.); 01235 aParticleChange.ProposeMomentumDirection(dir); 01236 aParticleChange.ProposeTrackStatus(fAlive); 01237 delete output; 01238 return G4VDiscreteProcess::PostStepDoIt(track,step); 01239 } 01240 #ifdef debug 01241 G4cout<<"G4QInelastic::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl; 01242 #endif 01243 aParticleChange.ProposeEnergy(0.); 01244 aParticleChange.ProposeTrackStatus(fStopAndKill); // the initial neutrino is killed 01245 // There is no way back from here ! 01246 if ( ((secnu || !nuanu || N) && totCS*G4UniformRand() < qelCS) || s_value < mlD2 ) 01247 { // Quasi-Elastic interaction 01248 G4double Q2=0.; // Simulate transferred momentum, in MeV^2 01249 if(secnu) Q2=CSmanager2->GetQEL_ExchangeQ2(); 01250 else Q2=CSmanager->GetQEL_ExchangeQ2(); 01251 #ifdef debug 01252 G4cout<<"G4QInelastic::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s_value<<",Q2="<<Q2<<G4endl; 01253 #endif 01254 //G4double ds=s_value+s_value; // doubled s_value 01255 G4double sqs=std::sqrt(s_value); // M_cm 01256 G4double dsqs=sqs+sqs; // 2*M_cm 01257 G4double p_init=(s_value-mIN*mIN)/dsqs; // initial momentum in CMS (checked MK) 01258 G4double dpi=p_init+p_init; // doubled initial momentum in CMS 01259 G4double sd=s_value-mlsOT; // s_value-ml2-mOT2 (mlsOT=m^2_neut+m^2_lept) 01260 G4double qo2=(sd*sd-mlOT)/dsqs; // squared momentum of secondaries in CMS 01261 G4double qo=std::sqrt(qo2); // momentum of secondaries in CMS 01262 G4double cost=(dpi*std::sqrt(qo2+ml2)-Q2-ml2)/dpi/qo; // cos(theta) in CMS (chck MK) 01263 G4LorentzVector t4M(0.,0.,0.,mIN); // 4mom of the effective target 01264 G4LorentzVector c4M=t4M+proj4M; // 4mom of the compound system 01265 t4M.setT(mOT); // now it is 4mom of the outgoing nucleon 01266 scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scattered muon 01267 if(!G4QHadron(c4M).RelDecayIn2(scat4M, t4M, proj4M, cost, cost)) 01268 { 01269 G4cerr<<"G4QIn::PStD:c4M="<<c4M<<sqs<<",mM="<<ml<<",tM="<<mOT<<",c="<<cost<<G4endl; 01270 // throw G4QException("G4QInelastic::HadronizeQuasm: Can't dec QE nu,lept Compound"); 01271 G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0000", 01272 FatalException, "Hadronize quasmon: Can't dec QE nu,lept Compound"); 01273 } 01274 proj4M=t4M; // 4mom of the new projectile nucleon 01275 } 01276 else // ***** Non Quasi Elastic interaction 01277 { 01278 // Recover the target PDG Code if the quasi-elastic scattering did not happened 01279 if ( (secnu && projPDG == 2212) || (!secnu && projPDG == 2112) ) targPDG+=1; 01280 else if ( (secnu && projPDG == 2112) || (!secnu && projPDG == 2212) ) targPDG+=1000; 01281 G4double Q2=0; // Simulate transferred momentum, in MeV^2 01282 if(secnu) Q2=CSmanager->GetNQE_ExchangeQ2(); 01283 else Q2=CSmanager2->GetNQE_ExchangeQ2(); 01284 #ifdef debug 01285 G4cout<<"G4QInel::PStDoIt: MultiPeriferal s="<<s_value<<",Q2="<<Q2<<",T="<<targPDG<<G4endl; 01286 #endif 01287 if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma 01288 else projPDG=CSmanager->GetExchangePDGCode(); // PDG Code of the effective pion 01289 //@@ Temporary made only for direct interaction and for N=3 (good for small Q2) 01290 //@@ inFuture use N=GetNPartons and directFraction=GetDirectPart, @@ W2... 01291 G4double r=G4UniformRand(); 01292 G4double r1=0.5; // (1-x) 01293 if(r<0.5) r1=std::sqrt(r+r)*(.5+.1579*(r-.5)); 01294 else if(r>0.5) r1=1.-std::sqrt(2.-r-r)*(.5+.1579*(.5-r)); 01295 G4double xn=1.-mldM/Momentum; // Normalization of (1-x) [x>mldM/Mom] 01296 G4double x1=xn*r1; // (1-x) 01297 G4double x=1.-x1; // x=2k/M 01298 //G4double W2=(hdM2+Q2/x)*x1; // W2 candidate 01299 G4double mx=hdM*x; // Part of the target to interact with 01300 G4double we=Q2/(mx+mx); // transfered energy 01301 if(we>=kinEnergy-ml-.001) we=kinEnergy-ml-.0001; // safety to avoid nan=sqrt(neg) 01302 G4double pot=kinEnergy-we; // energy of the secondary lepton 01303 G4double mlQ2=ml2+Q2; 01304 G4double cost=(pot-mlQ2/dKinE)/std::sqrt(pot*pot-ml2); // LS cos(theta) 01305 if(std::fabs(cost)>1) 01306 { 01307 #ifdef debug 01308 G4cout<<"*G4QInelastic::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx=" 01309 <<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl; 01310 #endif 01311 if(cost>1.) cost=1.; 01312 else cost=-1.; 01313 pot=mlQ2/dKinE+dKinE*ml2/mlQ2; // extreme output momentum 01314 } 01315 G4double lEn=std::sqrt(pot*pot+ml2); // Lepton energy 01316 G4double lPl=pot*cost; // Lepton longitudinal momentum 01317 G4double lPt=pot*std::sqrt(1.-cost*cost); // Lepton transverse momentum 01318 std::pair<G4double,G4double> d2d=Random2DDirection(); // Randomize phi 01319 G4double lPx=lPt*d2d.first; 01320 G4double lPy=lPt*d2d.second; 01321 G4ThreeVector vdir=proj4M.vect(); // 3D momentum of the projectile 01322 G4ThreeVector vz= vdir.unit(); // Ort in the direction of the projectile 01323 G4ThreeVector vv= vz.orthogonal(); // Not normed orthogonal vector (!) 01324 G4ThreeVector vx= vv.unit(); // First ort orthogonal to the direction 01325 G4ThreeVector vy= vz.cross(vx); // Second ort orthoganal to the direction 01326 G4ThreeVector lP= lPl*vz+lPx*vx+lPy*vy; // 3D momentum of the scattered lepton 01327 scat4M=G4LorentzVector(lP,lEn); // 4mom of the scattered lepton 01328 proj4M-=scat4M; // 4mom of the W/Z (effective pion/gamma) 01329 #ifdef debug 01330 G4cout<<"G4QInelastic::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl; 01331 #endif 01332 // Check that the en/mom transfer is possible, if not -> elastic 01333 G4int fintPDG=targPDG; // Prototype for the compound nucleus 01334 if(!secnu) 01335 { 01336 if(projPDG<0) fintPDG-= 999; 01337 else fintPDG+= 999; 01338 } 01339 G4double fM=G4QPDGCode(fintPDG).GetMass();// compound nucleus Mass (MeV) 01340 G4double fM2=fM*fM; 01341 G4LorentzVector tg4M=G4LorentzVector(0.,0.,0.,tgM); 01342 G4LorentzVector c4M=tg4M+proj4M; 01343 #ifdef debug 01344 G4cout<<"G4QInelastic::PSDI:fM2="<<fM2<<"<? mc4M="<<c4M.m2()<<",dM="<<fM-tgM<<G4endl; 01345 #endif 01346 if(fM2>=c4M.m2()) // Elastic scattering should be done 01347 { 01348 G4LorentzVector tot4M=tg4M+proj4M+scat4M; // recover the total 4-momentum 01349 s_value=tot4M.m2(); 01350 G4double fs=s_value-fM2-ml2; 01351 G4double fMl=fM2*ml2; 01352 G4double hQ2max=(fs*fs/2-fMl-fMl)/s_value; // Maximum possible Q2/2 01353 cost=1.-Q2/hQ2max; // cos(theta) in CMS (use MultProd Q2) 01354 #ifdef debug 01355 G4cout<<"G4QI::PSDI:ct="<<cost<<",Q2="<<Q2<<",hQ2="<<hQ2max<<",4M="<<tot4M<<G4endl; 01356 #endif 01357 G4double acost=std::fabs(cost); 01358 if(acost>1.) 01359 { 01360 if(acost>1.001) G4cout<<"-Warning-G4QInelastic::PostStDoIt: cost="<<cost<<G4endl; 01361 if (cost> 1.) cost= 1.; 01362 else if(cost<-1.) cost=-1.; 01363 } 01364 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,fM); // 4mom of the recoilNucleus 01365 scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scatteredLepton 01366 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-ml)*.01); 01367 if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) 01368 { 01369 G4cerr<<"G4QI::PSDI:t4M="<<tot4M<<",lM="<<ml<<",rM="<<fM<<",cost="<<cost<<G4endl; 01370 //G4Exception("G4QInelastic::PostStepDoIt:","027",FatalException,"ElasticDecay"); 01371 } 01372 #ifdef debug 01373 G4cout<<"G4QIn::PStDoI:l4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl; 01374 #endif 01375 // ---------------------------------------------------- 01376 G4ParticleDefinition* theDefinition=0; // Prototype of a particle for E-Secondaries 01377 // Fill scattered lepton 01378 if (scatPDG==-11) theDefinition = G4Positron::Positron(); 01379 else if(scatPDG== 11) theDefinition = G4Electron::Electron(); 01380 else if(scatPDG== 13) theDefinition = G4MuonMinus::MuonMinus(); 01381 else if(scatPDG==-13) theDefinition = G4MuonPlus::MuonPlus(); 01382 //else if(scatPDG== 15) theDefinition = G4TauMinus::TauMinus(); 01383 //else if(scatPDG==-15) theDefinition = G4TauPlus::TauPlus(); 01384 if (scatPDG==-12) theDefinition = G4AntiNeutrinoE::AntiNeutrinoE(); 01385 else if(scatPDG== 12) theDefinition = G4NeutrinoE::NeutrinoE(); 01386 else if(scatPDG== 14) theDefinition = G4NeutrinoMu::NeutrinoMu(); 01387 else if(scatPDG==-14) theDefinition = G4AntiNeutrinoMu::AntiNeutrinoMu(); 01388 //else if(scatPDG== 16) theDefinition = G4NeutrinoTau::NeutrinoTau(); 01389 //else if(scatPDG==-16) theDefinition = G4AntiNeutrinoTau::AntiNeutrinoTau(); 01390 else G4cout<<"-Warning-G4QInelastic::PostStDoIt: UnknownLepton="<<scatPDG<<G4endl; 01391 G4DynamicParticle* theScL = new G4DynamicParticle(theDefinition,scat4M); 01392 G4Track* scatLep = new G4Track(theScL, localtime, position ); // scattered 01393 scatLep->SetWeight(weight); // weighted 01394 scatLep->SetTouchableHandle(trTouchable); // residual 01395 aParticleChange.AddSecondary(scatLep); // lepton 01396 // Fill residual nucleus 01397 if (fintPDG==90000001) theDefinition = G4Neutron::Neutron(); // neutron 01398 else if(fintPDG==90001000) theDefinition = G4Proton::Proton(); // proton 01399 else // ion 01400 { 01401 G4int fm=static_cast<G4int>(fintPDG/1000000); // Strange part 01402 G4int ZN=fintPDG-1000000*fm; 01403 G4int rZ=static_cast<G4int>(ZN/1000); 01404 G4int rA=ZN-999*rZ; 01405 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ); 01406 } 01407 G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,reco4M); 01408 G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered 01409 scatReN->SetWeight(weight); // weighted 01410 scatReN->SetTouchableHandle(trTouchable); // residual 01411 aParticleChange.AddSecondary(scatReN); // nucleus 01412 delete output; 01413 return G4VDiscreteProcess::PostStepDoIt(track, step); 01414 } 01415 } // End of non-quasi-elastic interaction 01416 // 01417 // quasi-elastic (& pickup process) for p+A(Z,N) 01418 // 01419 } 01420 //else if ((projPDG == 2212 || projPDG == 2112) && Z > 0 && N > 0) // Never (in Fragm!) 01421 else if(2>3) // Possibility is closed as quasi-free scattering is made in fragmentation 01422 { 01423 if(momentum<500. && projPDG == 2112) // @@ It's reasonable to add proton capture too ! 01424 { 01425 G4LorentzVector tot4M=G4LorentzVector(0.,0.,0.,tgM)+proj4M; 01426 G4double totM2=tot4M.m2(); 01427 G4int tZ=Z; 01428 G4int tN=N; 01429 if(projPDG == 2212) tZ++; 01430 else tN++; // @@ Now a neutron is the only alternative 01431 if(momentum<100.) // Try the production threshold (@@ Why 5 MeV?) 01432 { 01433 G4QPDGCode FakeP(2112); 01434 //G4int m1p=tZ-1; 01435 G4int m2n=tN-2; 01436 // @@ For the secondary proton the Coulomb barrier should be taken into account 01437 //G4double mRP= FakeP.GetNuclMass(m1p,tN,0); // Residual to the secondary proton 01438 G4double mR2N= FakeP.GetNuclMass(tZ,m2n,0); // Residual to two secondary neutrons 01439 //G4double mRAl= FakeP.GetNuclMass(tZ-2,m2n,0); // Residual to the secondary alpha 01440 //G4double mR2N= FakeP.GetNuclMass(m1p,tN-1,0); // Residual to the p+n secondaries 01441 G4double minM=mR2N+mNeut+mNeut; 01442 // Compare with other possibilities (sciped for acceleration) 01443 if(totM2 < minM*minM) // DoNothing (under reasonable threshold) 01444 { 01445 aParticleChange.ProposeEnergy(kinEnergy); 01446 aParticleChange.ProposeLocalEnergyDeposit(0.); 01447 aParticleChange.ProposeMomentumDirection(dir); 01448 aParticleChange.ProposeTrackStatus(fAlive); 01449 return G4VDiscreteProcess::PostStepDoIt(track,step); 01450 } 01451 } 01452 } 01453 // Now the quasi-free and PickUp processes should be done: 01454 G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer(); 01455 std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, projPDG, Z, N); 01456 G4double qepart=fief.first*fief.second; // @@ charge exchange is not included @@ 01457 // Keep QE part for the quasi-free nucleons 01458 const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing 01459 G4double clProb[lCl]={.65,.85,.95};// N/P,D,t/He3,He4, integroProbab for .65,.2,.1,.05 01460 if(qepart>0.45) qepart=.45; // Quasielastic part is too large - shrink 01461 //else qepart/=clProb[0]; // Add corresponding number of 2N, 3N, & 4N clusters 01462 qepart=qepart/clProb[0]-qepart;// Add QE for 2N, 3N, & 4N clusters (N is made in G4QF) 01463 G4double pickup=1.-qepart; // Estimate the rest of the cross-section 01464 G4double thresh=100.; 01465 //if(momentum >thresh) pickup*=50./momentum/std::pow(G4double(Z+N),third);// 50 is Par 01466 if(momentum > thresh) pickup*=50./momentum/G4QThd(Z+N);// 50 is Par 01467 // pickup = 0.; // To exclude the pickup process 01468 if (N) pickup+=qepart; 01469 else pickup =qepart; 01470 G4double rnd=G4UniformRand(); 01471 #ifdef ldebug 01472 G4cout<<"-->G4QInelastic::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<")]="<<qepart 01473 <<", pickup="<<pickup<<G4endl; 01474 #endif 01475 if(rnd<pickup) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim,CoulBar,PauliBl 01476 { 01477 G4double dmom=91.; // Fermi momentum (proto default for a deuteron) 01478 G4double fmm=287.; // @@ Can be A-dependent 01479 if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third); 01480 // Values to be redefined for quasi-elastic 01481 G4LorentzVector r4M(0.,0.,0.,0.); // Prototype of 4mom of the residual nucleus 01482 G4LorentzVector n4M(0.,0.,0.,0.); // Prototype of 4mom of the quasi-cluster 01483 G4int nPDG=90000001; // Prototype for quasiCluster mass calculation 01484 G4int restPDG=targPDG; // Prototype should be reduced by quasiCluster 01485 G4int rA=Z+N-1; // Prototype for the residualNucl definition 01486 G4int rZ=Z; // residZ: OK for the quasi-free neutron 01487 G4int nA=1; // Prototype for the quasi-cluster definition 01488 G4int nZ=0; // nA=1,nZ=0: OK for the quasi-free neutron 01489 G4double qM=mNeut; // Free mass of the quasi-free cluster 01490 G4int A = Z + N; // Baryon number of the nucleus 01491 if(rnd<qepart) 01492 { 01493 // ===> First decay a nucleus in a nucleon and a residual (A-1) nucleus 01494 // Calculate clusterProbabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters) 01495 G4double base=1.; // Base for randomization (can be reduced by totZ & totN) 01496 G4int max=lCl; // Number of boundaries (can be reduced by totZ & totN) 01497 // Take into account that at least one nucleon must be left ! 01498 if(Z<2 || N<2 || A < 6) base = clProb[--max]; // Alpha cluster is impossible 01499 if((Z > 1 && N < 2) || (Z < 2 && N > 1)) base=(clProb[max]+clProb[max-1])/2;// t/He3 01500 if ( (Z < 2 && N < 2) || A < 5) base=clProb[--max];// He3&t clusters are impossible 01501 if(A<3) base=clProb[--max]; // Deuteron cluster is impossible 01502 //G4int cln=0; // Cluster#0 (Default for the selectedNucleon) 01503 G4int cln=1; // Cluster#1 (Default for the selectedDeutron) 01504 //if(max) // Not only nucleons are possible 01505 if(max>1) // Not only deuterons are possible 01506 { 01507 base-=clProb[0]; // Exclude scattering on QF Nucleon 01508 G4double ran=+clProb[0]+base*G4UniformRand(); // Base can be reduced 01509 G4int ic=1; // Start from the smallest cluster boundary 01510 if(max>1) while(ic<max) if(ran>clProb[ic++]) cln=ic; 01511 } 01512 G4ParticleDefinition* theDefinition=0;// Prototype for qfNucleon 01513 G4bool cp1 = cln+2==A; // A=ClusterBN+1 condition 01514 if(!cln || cp1) // Split in nucleon + (A-1) with Fermi momentum 01515 { 01516 G4int nln=0; 01517 if(cln==2) nln=1; // @@ only for cp1: t/He3 choice from A=4 01518 // mass(A)=tM. Calculate masses of A-1 (rM) and mN (mNeut or mProt bounded mass) 01519 if ( ((!cln || cln == 2) && G4UniformRand()*(A-cln) > (N-nln)) || 01520 ((cln == 3 || cln == 1) && Z > N) ) 01521 { 01522 nPDG=90001000; // Update quasi-free nucleon PDGCode to P 01523 nZ=1; // Change charge of the quasiFree nucleon 01524 qM=mProt; // Update quasi-free nucleon mass 01525 rZ--; // Reduce the residual Z 01526 restPDG-=1000; // Reduce the residual PDGCode 01527 } 01528 else restPDG--; 01529 G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed 01530 G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus 01531 r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus 01532 G4double rM2=rM*rM; 01533 G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-nucleon 01534 n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon 01535 #ifdef qedebug 01536 G4cout<<"G4QInel::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl; 01537 #endif 01538 if(!G4QHadron(t4M).DecayIn2(r4M, n4M)) 01539 { 01540 //G4cerr<<"G4QInel::PostStDoIt:M="<<tM<<"<r="<<rM<<"+n="<<nM<<"="<<rM+nM<<G4endl; 01541 //G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0002", 01542 // FatalException,"Hadronize quasmon: Can't Dec totNuc->QENuc+ResNuc"); 01543 G4cerr<<"-W-G4QInel::PoStDoI:M="<<tM<<"<rM="<<rM<<"+nM="<<nM<<"="<<rM+nM<<G4endl; 01544 aParticleChange.ProposeEnergy(kinEnergy); 01545 aParticleChange.ProposeLocalEnergyDeposit(0.); 01546 aParticleChange.ProposeMomentumDirection(dir); 01547 aParticleChange.ProposeTrackStatus(fAlive); 01548 return G4VDiscreteProcess::PostStepDoIt(track,step); 01549 } 01550 #ifdef qedebug 01551 G4cout<<"G4QIn::PStDoIt:QE-N, RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl; 01552 #endif 01553 if(cp1 && cln) // Quasi-cluster case: swap the output 01554 { 01555 qM=rM; // Scattering will be made on a cluster 01556 nln=nPDG; 01557 nPDG=restPDG; 01558 restPDG=nln; 01559 t4M=n4M; 01560 n4M=r4M; 01561 r4M=t4M; 01562 nln=nZ; 01563 nZ=rZ; 01564 rZ=nln; 01565 nln=nA; 01566 nA=rA; 01567 rA=nln; 01568 } 01569 } 01570 else // Split the cluster (with or w/o "Fermi motion" and "Fermi decay") 01571 { 01572 if(cln==1) 01573 { 01574 nPDG=90001001; // Deuteron 01575 qM=mDeut; 01576 nA=2; 01577 nZ=1; 01578 restPDG-=1001; 01579 // Recalculate dmom 01580 G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) + 01581 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01582 dmom=sumv.mag(); 01583 } 01584 else if(cln==2) 01585 { 01586 nA=3; 01587 if(G4UniformRand()*(A-2)>(N-1)) // He3 01588 { 01589 nPDG=90002001; 01590 qM=mHel3; 01591 nZ=2; 01592 restPDG-=2001; 01593 } 01594 else // tritium 01595 { 01596 nPDG=90001002; 01597 qM=mTrit; 01598 nZ=1; 01599 restPDG-=1002; 01600 } 01601 // Recalculate dmom 01602 G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) + 01603 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+ 01604 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01605 dmom=sumv.mag(); 01606 } 01607 else 01608 { 01609 nPDG=90002002; // Alpha 01610 qM=mAlph; 01611 nA=4; 01612 nZ=2; 01613 restPDG-=2002; 01614 G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) + 01615 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+ 01616 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+ 01617 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01618 dmom=sumv.mag(); 01619 } 01620 rA=A-nA; 01621 rZ=Z-nZ; 01622 // This is a simple case of cluster at rest 01623 //G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus 01624 //r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus 01625 //n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster 01626 // --- End of the "simple case of cluster at rest" 01627 // Make a fake quasi-Fermi distribution for clusters (clusters are not at rest) 01628 G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed 01629 G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus 01630 r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus 01631 G4double rM2=rM*rM; 01632 G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-cluster 01633 n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon 01634 #ifdef qedebug 01635 G4cout<<"G4QInel::PStDoIt:QEC, p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl; 01636 #endif 01637 if(!G4QHadron(t4M).DecayIn2(r4M, n4M)) 01638 { 01639 //G4cerr<<"G4QInel::PostStDoIt:M="<<tM<<"<r="<<rM<<"+c="<<nM<<"="<<rM+nM<<G4endl; 01640 //G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0003", 01641 // FatalException,"Hadronize quasmon: Can't Dec totNuc->QEClu+ResNuc"); 01642 G4cerr<<"-W-G4QInel::PoStDoI:M="<<tM<<"<rM="<<rM<<"+cM="<<nM<<"="<<rM+nM<<G4endl; 01643 aParticleChange.ProposeEnergy(kinEnergy); 01644 aParticleChange.ProposeLocalEnergyDeposit(0.); 01645 aParticleChange.ProposeMomentumDirection(dir); 01646 aParticleChange.ProposeTrackStatus(fAlive); 01647 return G4VDiscreteProcess::PostStepDoIt(track,step); 01648 } 01649 // --- End of the moving cluster implementation --- 01650 #ifdef qedebug 01651 G4cout<<"G4QIn::PStDoIt:QEC, RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl; 01652 #endif 01653 } 01654 G4LorentzVector s4M=n4M+proj4M; // Tot 4-momentum for scattering 01655 G4double prjM2 = proj4M.m2(); // Squared mass of the projectile 01656 G4double prjM = std::sqrt(prjM2); // @@ Get from pPDG (?) 01657 G4double minM = prjM+qM; // Min mass sum for the final products 01658 G4double cmM2 =s4M.m2(); 01659 if(cmM2>minM*minM) 01660 { 01661 #ifdef qedebug 01662 G4cout<<"G4QInel::PStDoIt:***Enter***,cmM2="<<cmM2<<" > minM2="<<minM*minM<<G4endl; 01663 #endif 01664 // Estimate and randomize charge-exchange with a quasi-free cluster 01665 G4bool chex=false; // Flag of the charge exchange scattering 01666 G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true 01667 // This is reserved for the charge-exchange scattering (to help neutron spectras) 01668 //if(cln&&!cp1 &&(projPDG==2212&&rA>rZ || projPDG==2112&&rZ>1))// @@ Use proj chex 01669 if(2>3) // @@ charge exchange is not implemented yet @@ 01670 { 01671 #ifdef qedebug 01672 G4cout<<"G4QIn::PStDoIt:-Enter, P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl; 01673 #endif 01674 G4double tprM=mProt; 01675 G4double tprM2=mProt2; 01676 G4int tprPDG=2212; 01677 G4int tresPDG=restPDG+999; 01678 if(projPDG==2212) 01679 { 01680 projpt=G4Neutron::Neutron(); 01681 tprM=mNeut; 01682 tprM2=mNeut2; 01683 tprPDG=2112; 01684 tresPDG=restPDG-999; 01685 } 01686 minM=tprM+qM; 01687 G4double efE=(cmM2-tprM2-qM*qM)/(qM+qM); 01688 G4double efP=std::sqrt(efE*efE-tprM2); 01689 G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!) 01690 #ifdef qedebug 01691 G4cout<<"G4QInel::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl; 01692 #endif 01693 if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl)) // minM is redefined 01694 { 01695 projPDG=tprPDG; 01696 prjM=tprM; 01697 G4double rM=G4QPDGCode(tresPDG).GetMass();// Mass of the residual nucleus 01698 r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus 01699 n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster 01700 chex=true; // Confirm charge exchange scattering 01701 } 01702 } 01703 // 01704 std::pair<G4LorentzVector,G4LorentzVector> sctout=qfMan->Scatter(nPDG, n4M, 01705 projPDG, proj4M); 01706 #ifdef qedebug 01707 G4cout<<"G4QInelastic::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM 01708 <<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl; 01709 #endif 01710 aParticleChange.ProposeLocalEnergyDeposit(0.); // Everything is in particles 01711 // @@ @@ @@ Coulomb barriers must be checked !! @@ @@ @@ Skip if not 01712 if(chex) // ==> Projectile is changed: fill everything to secondaries 01713 { 01714 aParticleChange.ProposeEnergy(0.); // @@ ?? 01715 aParticleChange.ProposeTrackStatus(fStopAndKill); // projectile nucleon is killed 01716 aParticleChange.SetNumberOfSecondaries(3); 01717 G4DynamicParticle* thePrH = new G4DynamicParticle(projpt,sctout.second); 01718 G4Track* scatPrH = new G4Track(thePrH, localtime, position ); // scattered & chex 01719 scatPrH->SetWeight(weight); // weighted 01720 scatPrH->SetTouchableHandle(trTouchable); // projectile 01721 aParticleChange.AddSecondary(scatPrH); // hadron 01722 } 01723 else // ==> The leading particle is filled to the updated projectilee 01724 { 01725 aParticleChange.SetNumberOfSecondaries(2); // @@ if proj=leading 01726 G4double ldT=(sctout.second).e()-prjM; // kin Energy of scat project. 01727 aParticleChange.ProposeEnergy(ldT); // Change the kin Energy 01728 G4ThreeVector ldV=(sctout.second).vect(); // Change momentum direction 01729 aParticleChange.ProposeMomentumDirection(ldV/ldV.mag()); 01730 aParticleChange.ProposeTrackStatus(fAlive); 01731 } 01732 // --------------------------------------------------------- 01733 // Fill scattered quasi-free nucleon/fragment 01734 if (nPDG==90000001) theDefinition = G4Neutron::Neutron(); 01735 else if(nPDG==90001000) theDefinition = G4Proton::Proton(); 01736 else if(nZ>0 && nA>1) 01737 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ); 01738 #ifdef debug 01739 else G4cout<<"-Warning_G4QIn::PSD:scatqfPDG="<<nPDG<<", Z="<<nZ<<",A="<<nA<<G4endl; 01740 #endif 01741 if(nZ>0 && nA>0) 01742 { 01743 G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,sctout.first); 01744 G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // scattered 01745 scatQFN->SetWeight(weight); // weighted 01746 scatQFN->SetTouchableHandle(trTouchable); // quasi-free 01747 aParticleChange.AddSecondary(scatQFN); // nucleon/cluster 01748 } 01749 // ---------------------------------------------------- 01750 // Fill residual nucleus 01751 if (restPDG==90000001) theDefinition = G4Neutron::Neutron(); 01752 else if(restPDG==90001000) theDefinition = G4Proton::Proton(); 01753 else if(rZ>0 && rA>1) 01754 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ); 01755 #ifdef debug 01756 else G4cout<<"-Warning_G4QIn::PSD: resPDG="<<restPDG<<",Z="<<rZ<<",A="<<rA<<G4endl; 01757 #endif 01758 if(rZ>0 && rA>0) 01759 { 01760 G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M); 01761 G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered 01762 scatReN->SetWeight(weight); // weighted 01763 scatReN->SetTouchableHandle(trTouchable); // residual 01764 aParticleChange.AddSecondary(scatReN); // nucleus 01765 } 01766 delete output; 01767 return G4VDiscreteProcess::PostStepDoIt(track, step); 01768 } 01769 #ifdef qedebug 01770 else G4cout<<"G4QInel::PSD:OUT,M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl; 01771 #endif 01772 } // end of proton quasi-elastic (including QE on NF) 01773 else // @@ make cost-condition @@ Pickup process (pickup 1 or 2 n and make d or t) 01774 { 01775 if(projPDG==2212) restPDG--; // Res. nucleus PDG is one neutron less than targ. PDG 01776 else 01777 { 01778 restPDG-=1000; // Res. nucleus PDG is one proton less than targ. PDG 01779 rZ--; // Reduce the charge of the residual nucleus 01780 } 01781 G4double mPUF=mDeut; // Prototype of the mass of the created pickup fragment 01782 G4ParticleDefinition* theDefinition = G4Deuteron::Deuteron(); // Default: make d 01783 //theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion 01784 G4bool din=false; // Flag of picking up 2 neutrons to create t (tritium)(p) 01785 G4bool dip=false; // Flag of picking up 2 protons to create t (tritium) (n) 01786 G4bool pin=false; // Flag of picking up 1 n + 1 p to create He3/t (p/n) 01787 G4bool hin=false; // Flag of PickUp creation of alpha (He4) (p/n) 01788 G4double frn=G4UniformRand(); 01789 if(N>2 && frn > 0.95) // .95 is a parameter to pickup 2N (to cteate t or He3) 01790 { 01791 if(projPDG==2212) 01792 { 01793 if(N>1 && G4UniformRand()*(Z+.5*(N-1)) > Z) 01794 { 01795 theDefinition = G4Triton::Triton(); 01796 mPUF=mTrit; // The mass of the created pickup fragment is triton 01797 restPDG--; // Res. nucleus PDG is two neutrons less than target PDG 01798 din=true; 01799 } 01800 else 01801 { 01802 theDefinition = G4He3::He3(); 01803 mPUF=mHel3; // The mass of the created pickup fragment is Helium3 01804 restPDG-=1000; // Res. nucleus PDG is two neutrons less than target PDG 01805 rZ--; // Reduce the charge of the residual nucleus 01806 pin=true; 01807 } 01808 } 01809 else // Neutron projectile (2112) 01810 { 01811 if(Z>1 && G4UniformRand()*(N+.5*(Z-1)) > N) 01812 { 01813 theDefinition = G4He3::He3(); 01814 mPUF=mHel3; // The mass of the created pickup fragment is Helium3 01815 restPDG-=1000; // Res. nucleus PDG is two neutrons less than target PDG 01816 rZ--; // Reduce the charge of the residual nucleus 01817 dip=true; 01818 } 01819 else 01820 { 01821 theDefinition = G4Triton::Triton(); 01822 mPUF=mTrit; // The mass of the created pickup fragment is triton 01823 restPDG--; // Res. nucleus PDG is two neutrons less than target PDG 01824 pin=true; 01825 } 01826 } 01827 rA--; // Reduce the baryon number of the residual nucleus 01828 // Recalculate dmom 01829 G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) + 01830 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01831 dmom=sumv.mag(); 01832 if(Z>1 && frn > 0.99) // .99 is a parameter to pickup 2n1p || 1n2p & create alpha 01833 { 01834 theDefinition = G4Alpha::Alpha(); 01835 if((din && projPDG==2212) || (pin && projPDG==2112)) 01836 { 01837 restPDG-=1000; // Residual nucleus PDG is reduced to create alpha 01838 rZ--; // Reduce the charge of the residual nucleus 01839 } 01840 else if((pin && projPDG==2212) || (dip && projPDG==2112)) restPDG--; 01841 else G4cout<<"-Warning-G4QIn::PSD: PickUp logic error, proj="<<projPDG<<G4endl; 01842 hin=true; 01843 mPUF=mAlph; // The mass of the created pickup fragment (alpha) 01844 rA--; // Reduce the baryon number of the residual nucleus 01845 // Recalculate dmom 01846 sumv += fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01847 dmom=sumv.mag(); 01848 } 01849 } 01850 G4double mPUF2=mPUF*mPUF; 01851 G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed 01852 G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus 01853 G4double rM2=rM*rM; // Squared mass of the residual nucleus 01854 G4double nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom);// M2 of boundQF-nucleon(2N) 01855 if(nM2 < 0) nM2=0.; 01856 G4double den2=(dmom*dmom+nM2); // squared energy of bQF-nucleon 01857 G4double den=std::sqrt(den2); // energy of bQF-nucleon 01858 #ifdef qedebug 01859 G4double nM=std::sqrt(nM2); // M of bQF-nucleon 01860 G4cout<<"G4QInel::PStDoIt:PiU, p="<<dmom<<",tM="<<tM<<", R="<<rM<<",N="<<nM<<G4endl; 01861 #endif 01862 G4double qp=momentum*dmom; 01863 G4double srm=proj4M.e()*den; 01864 G4double mNucl2=mProt2; 01865 if(projPDG == 2112) mNucl2=mNeut2; 01866 G4double cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp); 01867 G4int ict=0; 01868 while(std::fabs(cost)>1. && ict<3) 01869 { 01870 dmom=91.; // Fermi momentum (proDefault for a deuteron) 01871 if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third); 01872 if(din || pin || dip) // Recalculate dmom for the final t/He3 01873 { 01874 G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) + 01875 fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01876 if(hin) 01877 sumv+=fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection(); 01878 dmom=sumv.mag(); 01879 } 01880 nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom); // M2 of bQF-nucleon 01881 if(nM2 < 0) nM2=0.; 01882 //nM=std::sqrt(nM2); // M of bQF-nucleon --Not used?-- 01883 den2=(dmom*dmom+nM2); // squared energy of bQF-nucleon 01884 den=std::sqrt(den2); // energy of bQF-nucleon 01885 qp=momentum*dmom; 01886 srm=proj4M.e()*den; 01887 cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp); 01888 ict++; 01889 #ifdef ldebug 01890 if(ict>2)G4cout<<"G4QInelast::PStDoIt:i="<<ict<<",d="<<dmom<<",ct="<<cost<<G4endl; 01891 #endif 01892 } 01893 if(std::fabs(cost)<=1.) 01894 {// ---- @@ From here can be a MF for QF nucleon extraction (if used by others) 01895 G4ThreeVector vdir = proj4M.vect(); // 3-Vector in the projectile direction 01896 G4ThreeVector vx(0.,0.,1.); // ProtoOrt in theDirection of the projectile 01897 G4ThreeVector vy(0.,1.,0.); // First ProtoOrt orthogonal to the direction 01898 G4ThreeVector vz(1.,0.,0.); // SecondProtoOrt orthoganal to the direction 01899 if(vdir.mag2() > 0.) // the projectile isn't at rest 01900 { 01901 vx = vdir.unit(); // Ort in the direction of the projectile 01902 G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!) 01903 vy = vv.unit(); // First ort orthogonal to the proj direction 01904 vz = vx.cross(vy); // Second ort orthoganal to the projDirection 01905 } 01906 // ---- @@ End of possible MF (Similar is in G4QEnvironment) 01907 G4double tdom=dmom*std::sqrt(1-cost*cost);// Transversal component of the Fermi-Mom 01908 G4double phi=twopi*G4UniformRand(); // Phi of the Fermi-Mom 01909 G4ThreeVector vedm=vx*dmom*cost+vy*tdom*std::sin(phi)+vz*tdom*std::cos(phi);// bQFN 01910 G4LorentzVector bqf(vedm,den); // 4-mom of the bQF nucleon 01911 r4M=t4M-bqf; // 4mom of the residual nucleus 01912 #ifdef debug 01913 if(std::fabs(r4M.m()-rM)>.001) G4cout<<"G4QIn::PSD: rM="<<rM<<"#"<<r4M.m()<<G4endl; 01914 #endif 01915 G4LorentzVector f4M=proj4M+bqf; // Prototype of 4-mom of Deuteron 01916 #ifdef debug 01917 if(std::fabs(f4M.m()-mPUF)>.001)G4cout<<"G4QI::PSD:m="<<mPUF<<"#"<<f4M.m()<<G4endl; 01918 #endif 01919 aParticleChange.ProposeEnergy(0.); // @@ ?? 01920 aParticleChange.ProposeTrackStatus(fStopAndKill);// projectile nucleon is killed 01921 aParticleChange.SetNumberOfSecondaries(2); 01922 // Fill pickup hadron 01923 G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,f4M); 01924 G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // pickup 01925 scatQFN->SetWeight(weight); // weighted 01926 scatQFN->SetTouchableHandle(trTouchable); // nuclear 01927 aParticleChange.AddSecondary(scatQFN); // cluster 01928 #ifdef pickupd 01929 G4cout<<"G4QInelastic::PStDoIt: f="<<theDefinition<<",f4M="<<f4M.m()<<f4M<<G4endl; 01930 #endif 01931 // ---------------------------------------------------- 01932 // Fill residual nucleus 01933 if (restPDG==90000001) theDefinition = G4Neutron::Neutron(); 01934 else if(restPDG==90001000) theDefinition = G4Proton::Proton(); 01935 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion 01936 G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M); 01937 G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered 01938 scatReN->SetWeight(weight); // weighted 01939 scatReN->SetTouchableHandle(trTouchable); // residual 01940 aParticleChange.AddSecondary(scatReN); // nucleus 01941 #ifdef pickupd 01942 G4cout<<"G4QInelastic::PStDoIt:rZ="<<rZ<<",rA="<<rA<<",r4M="<<r4M.m()<<r4M<<G4endl; 01943 #endif 01944 delete output; 01945 return G4VDiscreteProcess::PostStepDoIt(track, step); 01946 } 01947 } 01948 } // end of the quasi-elastic decoupled process for the neucleon projectile 01949 } // end lepto-nuclear, neutrino-nuclear, proton/neutron decoupled processes 01950 EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction 01951 if(absMom) EnMomConservation+=lead4M; // Add E/M of leading System 01952 #ifdef debug 01953 G4cout<<"G4QInelastic::PostStDoIt:before St="<<aParticleChange.GetTrackStatus()<<G4endl; 01954 #endif 01955 01956 // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin 01957 //G4bool elF=false; // Flag of the ellastic scattering is "false" by default 01958 //G4double eWei=1.; 01959 // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End 01960 #ifdef debug 01961 G4cout<<"^^G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl; 01962 #endif 01963 const G4QNucleus targNuc(Z,N); // Define the target nucleus 01964 if(projPDG>9000000) 01965 { 01966 delete output; // Before the output creation 01967 G4QNucleus projNuc(proj4M,projPDG); // Define the projectile nucleus 01968 G4QIonIonCollision IonIon(projNuc, targNuc); // Define deep-inelastic ion-ion reaction 01969 #ifdef debug 01970 G4cout<<"G4QInel::PStDoIt: ProjNuc="<<projPDG<<proj4M<<", TargNuc="<<targPDG<<G4endl; 01971 #endif 01972 try {output = IonIon.Fragment();} // DINR AA interaction products 01973 catch (G4QException& error) 01974 { 01975 G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl; 01976 // G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash"); 01977 G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0027", 01978 FatalException, "CHIPS hA crash"); 01979 } 01980 } 01981 else // --> A projectile hadron 01982 { 01983 // *** Let to produce elastic final states for acceleration *** 01984 #ifdef debug 01985 G4int maxCn=7; 01986 #endif 01987 G4int atCn=0; // Attempts counter 01988 //G4bool inel=false; 01989 //while (inel==false && ++atCn <= maxCn) // To evoid elastic final state 01990 //{ 01991 #ifdef debug 01992 G4cout<<"G4QInelastic::PStDoIt: atCn="<<atCn<<" <= maxCn="<<maxCn<<G4endl; 01993 #endif 01994 G4int outN=output->size(); 01995 if(outN) // The output is not empty 01996 { 01997 std::for_each(output->begin(), output->end(), DeleteQHadron()); 01998 output->clear(); 01999 } 02000 delete output; // Before the new output creation 02001 const G4QHadron projH(projPDG,proj4M); // Define the projectile particle 02002 #ifdef debug 02003 G4cout<<"G4QInelastic::PStDoIt: Before Fragmentation"<<G4endl; 02004 #endif 02005 //if(aProjPDG==22 && projPDG==22 && proj4M.m2()<.01 && proj4M.e()<150.) 02006 //{ 02007 // G4QHadron* tgH= new G4QHadron(targNuc.GetPDGCode(),targNuc.Get4Momentum()); 02008 // G4QHadron* prH= new G4QHadron(projPDG,proj4M); 02009 // output = new G4QHadronVector(); 02010 // output->push_back(prH); 02011 // output->push_back(tgH); 02012 //} 02013 //else 02014 //{ 02015 G4QFragmentation DINR(targNuc, projH); // Define deep-inel. h-A reaction 02016 #ifdef debug 02017 G4cout<<"G4QInelastic::PStDoIt:Proj="<<projPDG<<proj4M<<",TgNuc="<<targPDG<<G4endl; 02018 #endif 02019 try {output = DINR.Fragment();} // DINR hA fragmentation 02020 catch (G4QException& error) 02021 { 02022 G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl; 02023 // G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash"); 02024 G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0027", 02025 FatalException, "CHIPS hA crash"); 02026 } 02027 //} 02028 outN=output->size(); 02029 #ifdef debug 02030 G4cout<<"G4QInelastic::PostStepDoIt: At# "<<atCn<<", nSec="<<outN<<", fPDG=" 02031 <<(*output)[0]->GetPDGCode()<<", pPDG="<< projPDG<<G4endl; 02032 #endif 02033 //inel=true; // For while 02034 if(outN < 2) 02035 { 02036 G4cout<<"-Warning-G4QInelastic::PostStepDoIt: nSec="<<outN<<", At# "<<atCn<<G4endl; 02037 //inel=false; // For while 02038 } 02039 //else if(outN==2 && (*output)[0]->GetPDGCode() == projPDG) inel=false; // For while 02040 #ifdef debug 02041 if(atCn==maxCn)G4cout<<"-Warning-G4QI::PostStDoIt:mAt="<<atCn<<" is reached"<<G4endl; 02042 #endif 02043 //} 02044 } 02045 // 02046 // --- the scattered hadron with changed nature can be added here --- 02047 if(scat) 02048 { 02049 G4QHadron* scatHadron = new G4QHadron(scatPDG,scat4M); 02050 output->push_back(scatHadron); 02051 } 02052 G4int qNH=0; 02053 if(leadhs) qNH=leadhs->size(); 02054 if(absMom) 02055 { 02056 if(qNH) for(G4int iq=0; iq<qNH; iq++) 02057 { 02058 G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron 02059 output->push_back(loh); 02060 } 02061 if(leadhs) delete leadhs; 02062 leadhs=0; 02063 } 02064 else 02065 { 02066 if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq]; 02067 if(leadhs) delete leadhs; 02068 leadhs=0; 02069 } 02070 // ------------- From here the secondaries are filled ------------------------- 02071 G4int tNH = output->size(); // A#of hadrons in the output 02072 aParticleChange.SetNumberOfSecondaries(tNH); // @@ tNH can be changed (drop/decay!) 02073 // Now add nuclear fragments 02074 #ifdef debug 02075 G4cout<<"G4QInelastic::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl; 02076 #endif 02077 #ifdef ppdebug 02078 if(absMom)G4cout<<"G4QInelastic::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl; 02079 #endif 02080 G4int nOut=output->size(); // Real length of the output @@ Temporary 02081 if(tNH==1 && !scat) // @@ Temporary. Find out why it happened! 02082 { 02083 G4cout<<"-Warning-G4QInelastic::PostStepDoIt: 1 secondary! absMom="<<absMom; 02084 if(absMom) G4cout<<", qNH="<<qNH; 02085 G4cout<<", PDG0="<<(*output)[0]->GetPDGCode(); 02086 G4cout<<G4endl; 02087 tNH=0; 02088 delete output->operator[](0); // delete the creazy hadron 02089 output->pop_back(); // clean up the output vector 02090 } 02091 if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QInelastic::PostStepDoIt: 2 # "<<nOut<<G4endl; 02092 // Deal with ParticleChange final state interface to GEANT4 output of the process 02093 //if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH 02094 if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH 02095 { 02096 // Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside) 02097 // Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@ 02098 // The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body) 02099 G4QHadron* hadr=(*output)[i]; // Pointer to the output hadron 02100 G4int PDGCode = hadr->GetPDGCode(); 02101 G4int nFrag = hadr->GetNFragments(); 02102 #ifdef pdebug 02103 G4cout<<"G4QInelastic::PostStepDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag 02104 <<", 4Mom="<<hadr->Get4Momentum()<<G4endl; 02105 #endif 02106 if(nFrag) // Skip intermediate (decayed) hadrons 02107 { 02108 #ifdef debug 02109 G4cout<<"G4QInelastic::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl; 02110 #endif 02111 delete hadr; 02112 continue; 02113 } 02114 G4DynamicParticle* theSec = new G4DynamicParticle; 02115 G4ParticleDefinition* theDefinition; 02116 if (PDGCode==90000001) theDefinition = G4Neutron::Neutron(); 02117 else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions 02118 else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda(); 02119 else if(PDGCode==311 || PDGCode==-311) 02120 { 02121 if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L 02122 else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S 02123 } 02124 else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus(); 02125 else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus(); 02126 else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus(); 02127 else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero(); 02128 else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus(); 02129 else if(PDGCode==90000000) 02130 { 02131 #ifdef pdebug 02132 G4cout<<"-Warning-G4QInelastic::PostStepDoIt:Vacuum PDG="<<PDGCode 02133 <<", 4Mom="<<hadr->Get4Momentum()<<G4endl; 02134 #endif 02135 theDefinition = G4Gamma::Gamma(); 02136 G4double hM2=(hadr->Get4Momentum()).m2(); 02137 if(std::fabs(hM2)>.01) G4cout<<"-Warning-G4QInel::PoStDoIt:90000000M2="<<hM2<<G4endl; 02138 else if(hadr->Get4Momentum()==vacuum4M) 02139 { 02140 delete hadr; 02141 continue; 02142 } 02143 } 02144 else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!) 02145 { 02146 G4int aZ = hadr->GetCharge(); 02147 G4int aA = hadr->GetBaryonNumber(); 02148 #ifdef pdebug 02149 G4cout<<"G4QInelastic::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl; 02150 #endif 02151 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ); 02152 } 02153 //else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode); 02154 else if(PDGCode==89999998 || PDGCode==89998000 || PDGCode==88000000) // anti-di-baryon 02155 { 02156 G4double rM=mNeut; // Prototype of the baryon mass 02157 G4int rPDG=-2112; // Prototype of the baryon PDG 02158 G4ParticleDefinition* BarDef=G4Neutron::Neutron();// Prototype of theBaryonDefinition 02159 if (PDGCode==89998000) 02160 { 02161 rM=mProt; 02162 rPDG=-2212; 02163 BarDef=G4Proton::Proton(); 02164 } 02165 else if(PDGCode==88000000) 02166 { 02167 rM=mLamb; 02168 rPDG=-3122; 02169 BarDef=G4Lambda::Lambda(); 02170 } 02171 G4LorentzVector t4M=hadr->Get4Momentum(); // 4m of the di-baryon 02172 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 1st anti-baryon 02173 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 2nd anti-baryon 02174 #ifdef qedebug 02175 G4cout<<"G4QInel::PStDoIt:AntiDiBar,t4M="<<tM<<",m="<<rM<<",PDG="<<PDGCode<<G4endl; 02176 #endif 02177 if(!G4QHadron(t4M).DecayIn2(f4M, s4M)) 02178 { 02179 G4cerr<<"G4QIn::PostStDoIt: ADB, M="<<t4M.m()<<" < 2*rM="<<rM<<" = "<<2*rM<<G4endl; 02180 // throw G4QException("G4QInelastic::HadronizeQuasm:Can't decay anti-dibaryon"); 02181 G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0004", 02182 FatalException, "Hadronize quasmon: Can't decay anti-dibaryon"); 02183 } 02184 // --- End of the moving cluster implementation --- 02185 #ifdef qedebug 02186 G4cout<<"G4QInelastic::PStDoIt: ADB, H1="<<rM<<f4M<<", H2="<<rM<<s4M<<G4endl; 02187 #endif 02188 G4QHadron fH(rPDG,f4M); 02189 hadr->Set4Momentum(f4M); 02190 hadr->SetQPDG(fH.GetQPDG()); 02191 theDefinition=BarDef; 02192 #ifdef debug 02193 G4cout<<"G4QInel::PostStDoIt:Anti-DiBar, DecayIn2, h1="<<rPDG<<f4M<<G4endl; 02194 #endif 02195 G4QHadron* sH = new G4QHadron(rPDG,s4M); 02196 output->push_back(sH); 02197 ++tNH; 02198 #ifdef debug 02199 G4cout<<"G4QInel::PostStDoIt:Anti-DiBar, DecayIn2, h2="<<rPDG<<s4M<<G4endl; 02200 #endif 02201 } 02202 else if(PDGCode==90000997 || PDGCode==89997001) // anti-NDelta 02203 { 02204 G4double rM=mNeut; // Prototype of the baryon mass 02205 G4int rPDG=-2112; // Prototype of the baryon PDG 02206 G4double iM=mPi; // Prototype of the pion mass 02207 G4int iPDG= 211; // Prototype of the pion PDG 02208 G4ParticleDefinition* BarDef=G4Neutron::Neutron();// Prototype of theBaryonDefinition 02209 if(PDGCode==90000997) 02210 { 02211 rM=mProt; 02212 rPDG=-2212; 02213 iPDG=-211; 02214 BarDef=G4Proton::Proton(); 02215 } 02216 G4LorentzVector t4M=hadr->Get4Momentum(); // 4m of the di-baryon 02217 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 1st anti-baryon 02218 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 2nd anti-baryon 02219 G4LorentzVector u4M=G4LorentzVector(0.,0.,0.,iM); // 4mom of the pion 02220 #ifdef qedebug 02221 G4cout<<"G4QInel::PStDoIt:AntiNDelta, t4M="<<tM<<",m="<<rM<<",PDG="<<PDGCode<<G4endl; 02222 #endif 02223 if(!G4QHadron(t4M).DecayIn3(f4M, s4M, u4M)) 02224 { 02225 G4cerr<<"G4QIn::PostStDoIt: AND, tM="<<t4M.m()<<" < 2*mB+mPi="<<2*rM+iM<<G4endl; 02226 // throw G4QException("G4QInelastic::HadronizeQuasm:Can't decay anti-NDelta"); 02227 G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0005", 02228 FatalException, "Hadronize quasmon: Can't decay anti-NDelta"); 02229 } 02230 // --- End of the moving cluster implementation --- 02231 #ifdef qedebug 02232 G4cout<<"G4QInel::PStDoI:AND,B1="<<rM<<f4M<<",B2="<<rM<<s4M<<",Pi="<<iM<<u4M<<G4endl; 02233 #endif 02234 G4QHadron fH(rPDG,f4M); 02235 hadr->Set4Momentum(f4M); 02236 hadr->SetQPDG(fH.GetQPDG()); 02237 theDefinition=BarDef; 02238 #ifdef debug 02239 G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h1="<<rPDG<<f4M<<G4endl; 02240 #endif 02241 G4QHadron* sH = new G4QHadron(rPDG,s4M); 02242 output->push_back(sH); 02243 ++tNH; 02244 #ifdef debug 02245 G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h2="<<rPDG<<s4M<<G4endl; 02246 #endif 02247 G4QHadron* uH = new G4QHadron(iPDG,u4M); 02248 output->push_back(uH); 02249 ++tNH; 02250 #ifdef debug 02251 G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h2="<<rPDG<<s4M<<G4endl; 02252 #endif 02253 } 02254 else 02255 { 02256 #ifdef pdebug 02257 G4cout<<"G4QInelastic::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl; 02258 #endif 02259 theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode); 02260 #ifdef pdebug 02261 G4cout<<"G4QInelastic::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl; 02262 #endif 02263 } 02264 if(!theDefinition) 02265 { // !! Commenting of this print costs days of searching for E/mom nonconservation !! 02266 if(PDGCode!=90000000 || hadr->Get4Momentum()!=vacuum4M) 02267 G4cout<<"---Warning---G4QInelastic::PostStepDoIt: drop PDG="<<PDGCode<<", 4Mom=" 02268 <<hadr->Get4Momentum()<<G4endl; 02269 delete hadr; 02270 continue; 02271 } 02272 #ifdef pdebug 02273 G4cout<<"G4QInelastic::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl; 02274 #endif 02275 theSec->SetDefinition(theDefinition); 02276 G4LorentzVector h4M=hadr->Get4Momentum(); 02277 EnMomConservation-=h4M; 02278 #ifdef tdebug 02279 G4cout<<"G4QInelast::PSDI: "<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl; 02280 #endif 02281 #ifdef debug 02282 G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl; 02283 #endif 02284 theSec->Set4Momentum(h4M); // ^ 02285 delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+ 02286 #ifdef debug 02287 G4ThreeVector curD=theSec->GetMomentumDirection(); // ^ 02288 G4double curM=theSec->GetMass(); // | 02289 G4double curE=theSec->GetKineticEnergy()+curM; // ^ 02290 G4cout<<"G4QInelast::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;//| 02291 #endif 02292 G4Track* aNewTrack = new G4Track(theSec, localtime, position ); // ^ 02293 aNewTrack->SetWeight(weight); // weighted | 02294 aNewTrack->SetTouchableHandle(trTouchable); // | 02295 aParticleChange.AddSecondary( aNewTrack ); // | 02296 #ifdef debug 02297 G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<" is done."<<G4endl; // | 02298 #endif 02299 } // | 02300 delete output; // instances of the G4QHadrons from the output are already deleted above + 02301 if(leadhs) // To satisfy Valgrind ( How can that be?) 02302 { 02303 qNH=leadhs->size(); 02304 if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq]; 02305 delete leadhs; 02306 leadhs=0; 02307 } 02308 #ifdef debug 02309 G4cout<<"G4QInelastic::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl; 02310 #endif 02311 if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15) 02312 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the absorbed particle 02313 #ifdef pdebug 02314 G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy() 02315 <<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl; 02316 #endif 02317 #ifdef ldebug 02318 G4cout<<"G4QInelastic::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG<<", St=" 02319 <<aParticleChange.GetTrackStatus()<<G4endl; 02320 #endif 02321 return G4VDiscreteProcess::PostStepDoIt(track, step); 02322 }
| void G4QInelastic::PrintInfoDefinition | ( | ) | [inline] |
| void G4QInelastic::SetManual | ( | ) | [static] |
| void G4QInelastic::SetParameters | ( | G4double | temper = 180., |
|
| G4double | ssin2g = .1, |
|||
| G4double | etaetap = .3, |
|||
| G4double | fN = 0., |
|||
| G4double | fD = 0., |
|||
| G4double | cP = 1., |
|||
| G4double | mR = 1., |
|||
| G4int | npCHIPSWorld = 234, |
|||
| G4double | solAn = .5, |
|||
| G4bool | efFlag = false, |
|||
| G4double | piTh = 141.4, |
|||
| G4double | mpi2 = 20000., |
|||
| G4double | dinum = 1880. | |||
| ) | [static] |
Definition at line 98 of file G4QInelastic.cc.
References G4QCHIPSWorld::Get(), G4QCHIPSWorld::GetParticles(), G4QEnvironment::SetParameters(), G4Quasmon::SetParameters(), and G4QNucleus::SetParameters().
00102 { 00103 Temperature=temper; 00104 SSin2Gluons=ssin2g; 00105 EtaEtaprime=etaetap; 00106 freeNuc=fN; 00107 freeDib=fD; 00108 clustProb=cP; 00109 mediRatio=mR; 00110 nPartCWorld = nParCW; 00111 EnergyFlux=efFlag; 00112 SolidAngle=solAn; 00113 PiPrThresh=piThresh; 00114 M2ShiftVir=mpisq; 00115 DiNuclMass=dinum; 00116 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles 00117 G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's 00118 G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters 00119 G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture 00120 }
| void G4QInelastic::SetPhotNucBias | ( | G4double | phnB = 1. |
) | [static] |
Definition at line 122 of file G4QInelastic.cc.
Referenced by G4QPhotoNuclearPhysics::ConstructProcess().
| void G4QInelastic::SetStandard | ( | ) | [static] |
| void G4QInelastic::SetWeakNucBias | ( | G4double | ccnB = 1. |
) | [static] |
Definition at line 123 of file G4QInelastic.cc.
Referenced by G4QNeutrinoPhysics::ConstructProcess().
1.4.7