// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: G4QCaptureAtRest.cc,v 1.19 2009/01/24 11:57:46 vnivanch Exp $ // GEANT4 tag $Name: geant4-09-02-ref-02 $ // // ---------------- G4QCaptureAtRest class ----------------- // by Mikhail Kossov, December 2003. // G4QCaptureAtRest class of the CHIPS Simulation Branch in GEANT4 // --------------------------------------------------------------- // **************************************************************************************** // ********** This CLASS is temporary moved from the photolepton_hadron directory ********* // ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ****** // **************************************************************************************** //#define debug //#define pdebug //#define tdebug #include "G4QCaptureAtRest.hh" #include "G4HadronicProcessStore.hh" G4QCaptureAtRest::G4QCaptureAtRest(const G4String& processName) : G4VRestProcess(processName, fHadronic), Time(0.), EnergyDeposition(0.) { SetProcessSubType(fCapture); #ifdef debug G4cout<<"G4QCaptureAtRest::Constructor is called"<0) G4cout << GetProcessName() << " is created "<< G4endl; //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture G4HadronicProcessStore::Instance()->RegisterExtraProcess(this); } G4bool G4QCaptureAtRest::manualFlag=false; // If false then standard parameters are used G4double G4QCaptureAtRest::Temperature=180.; // Critical Temperature (sensitive at High En) G4double G4QCaptureAtRest::SSin2Gluons=0.3; // Supression of s-quarks (in respect to u&d) G4double G4QCaptureAtRest::EtaEtaprime=0.3; // Supression of eta mesons (gg->qq/3g->qq) G4double G4QCaptureAtRest::freeNuc=0.5; // Percentage of free nucleons on the surface G4double G4QCaptureAtRest::freeDib=0.05; // Percentage of free diBaryons on the surface G4double G4QCaptureAtRest::clustProb=5.; // Nuclear clusterization parameter G4double G4QCaptureAtRest::mediRatio=10.; // medium/vacuum hadronization ratio G4int G4QCaptureAtRest::nPartCWorld=152; // The#of particles initialized in CHIPS World G4double G4QCaptureAtRest::SolidAngle=0.5; // Part of Solid Angle to capture (@@A-dep.) G4bool G4QCaptureAtRest::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon) G4double G4QCaptureAtRest::PiPrThresh=141.4; // Pion Production Threshold for gammas G4double G4QCaptureAtRest::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma G4double G4QCaptureAtRest::DiNuclMass=1880.; // DoubleNucleon Mass for VirtualNormalization void G4QCaptureAtRest::SetManual() {manualFlag=true;} void G4QCaptureAtRest::SetStandard() {manualFlag=false;} // Fill the private parameters void G4QCaptureAtRest::SetParameters(G4double temper, G4double ssin2g, G4double etaetap, G4double fN, G4double fD, G4double cP, G4double mR, G4int nParCW, G4double solAn, G4bool efFlag, G4double piThresh, G4double mpisq, G4double dinum) {// ============================================================================= Temperature=temper; SSin2Gluons=ssin2g; EtaEtaprime=etaetap; freeNuc=fN; freeDib=fD; clustProb=cP; mediRatio=mR; nPartCWorld = nParCW; EnergyFlux=efFlag; SolidAngle=solAn; PiPrThresh=piThresh; M2ShiftVir=mpisq; DiNuclMass=dinum; G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture } // Destructor G4QCaptureAtRest::~G4QCaptureAtRest() { // deregister in the store G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this); } G4LorentzVector G4QCaptureAtRest::GetEnegryMomentumConservation() { return EnMomConservation; } G4int G4QCaptureAtRest::GetNumberOfNeutronsInTarget() { return nOfNeutrons; } G4bool G4QCaptureAtRest::IsApplicable(const G4ParticleDefinition& particle) { if (particle == *( G4PionMinus::PionMinus() )) return true; else if (particle == *( G4KaonMinus::KaonMinus() )) return true; else if (particle == *( G4AntiProton::AntiProton() )) return true; else if (particle == *( G4MuonMinus::MuonMinus() )) return true; else if (particle == *( G4TauMinus::TauMinus() )) return true; else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true; else if (particle == *( G4XiMinus::XiMinus() )) return true; else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true; else if (particle == *( G4Neutron::Neutron() )) return true; else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true; else if (particle == *(G4AntiSigmaPlus::AntiSigmaPlus())) return true; #ifdef debug G4cout<<"***G4QCaptureAtRest::IsApplicable: PDG="<RegisterParticleForExtraProcess(this, &p); } void G4QCaptureAtRest::BuildPhysicsTable(const G4ParticleDefinition& p) { G4HadronicProcessStore::Instance()->PrintInfo(&p); } G4VParticleChange* G4QCaptureAtRest::AtRestDoIt(const G4Track& track, const G4Step& step) { static const G4double mNeut = G4QPDGCode(2112).GetMass(); static const G4double mNeut2= mNeut*mNeut; static const G4double dmNeut= mNeut+mNeut; static const G4double mProt = G4QPDGCode(2212).GetMass(); static const G4double dmProt= mProt+mProt; static const G4double mPi0 = G4QPDGCode(111).GetMass(); static const G4double mDeut = G4QPDGCode(2112).GetNuclMass(1,1,0); static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0); //static const G4double mPi = G4QPDGCode(211).GetMass(); static const G4double mMu = G4QPDGCode(13).GetMass(); //static const G4double mMu2 = mMu*mMu; //static const G4double dmMu = mMu+mMu; // The best parameters for mu->e+nu+nu decay static const G4double b1=2.784; static const G4double ga=.000015; static const G4double rb=1./b1; //static const G4double mTau = G4QPDGCode(15).GetMass(); static const G4double mEl = G4QPDGCode(11).GetMass(); static const G4double mEl2 = mEl*mEl; //------------------------------------------------------------------------------------- static G4bool CWinit = true; // CHIPS Warld needs to be initted if(CWinit) { CWinit=false; G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) } //------------------------------------------------------------------------------------- const G4DynamicParticle* stoppedHadron = track.GetDynamicParticle(); const G4ParticleDefinition* particle=stoppedHadron->GetDefinition(); Time=0.; EnergyDeposition=0.; #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt is called,EnDeposition="<GetElementVector(); G4int i=0; G4double sum=0.; G4int nE=material->GetNumberOfElements(); #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: "< sumfra; for(i=0; iGetFractionVector()[i]; if(projPDG==13||projPDG==15) { G4int cZ=static_cast((*theElementVector)[i]->GetZ()); frac*=cZ; if(cZ==9||cZ==35||cZ==53||cZ==85) frac*=.66; else if (cZ== 3) frac*=.50; else if (cZ==24||cZ==28) frac*=.90; else if (cZ== 5||cZ==17) frac*=.70; else if (cZ== 8) frac*=.56; } sum+=frac; sumfra.push_back(sum); // remember the summation steps } G4double rnd = sum*G4UniformRand(); for(i=0; i(pElement->GetZ()); if(Z<=0) { G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt:Element with Z="<GetIsotopeVector(); if(isoVector) isoSize=isoVector->size(); // Get real size of the isotopeVector if exists #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: isovectorLength="<GetIndex()+1; // Index of the non-trivial element is an order #ifdef debug G4cout<<"G4QCapAR::GetMFP: iE="<* pr= new std::pair(N,abund); #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt:pair#="<n+A,Z="<push_back(secnuc); // Fill recoil nucleus to the output G4QHadron* neutron = new G4QHadron(2112,n4Mom); // Create Hadron for the Neutron output->push_back(neutron); // Fill the neutron to the output #ifdef debug CV=27; G4cout<<"G4QCaptureAtRest::AtRestDoIt:ElasN="<push_back(secnuc); // Fill recharged nucleus to the output G4QHadron* proton = new G4QHadron(2212,p4Mom); // Create Hadron for the Proton output->push_back(proton) ; // Fill the proton to the output #ifdef debug CV=21; G4cout<<"G4QCaptureAtRest::AtRestDoIt:ChExP="<n+pi0)/p+pi-=>n+gamma) = 3/2 #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: Panofsky targPDG="<0.6) { pigamPDG=22; pigamM=0.; } G4LorentzVector g4Mom(0.,0.,0.,pigamM); // mass of the photon/Pi0 G4LorentzVector n4Mom(0.,0.,0.,mNeut); // mass of the secondary neutron if(!G4QHadron(totLV).DecayIn2(g4Mom,n4Mom)) { G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: H1+(pi-)=>n+"<push_back(pigamma); // Fill pi0 or Gamma in the output G4QHadron* neutron = new G4QHadron(2112,n4Mom); // Create Hadron for the Neutron output->push_back(neutron); // Fill the neutron to the output } // @@ For pi-,d reactions one can use just nn dedcay (see above) ? // @@ For K- Capture the quasifree n+Lambda+(A-n-p) reaction can be applyed as well... else if(projPDG==-211 && G4UniformRand()>1.&& Z>0&&N>0) // @@Quasi-Free PiCapture => tune { G4double mt=G4QPDGCode(targPDG).GetMass();// Mass of the target Nucleus G4LorentzVector totLV(0.,0.,0.,mp+mt); // 4-momentum of the (A+pi-) compound system if(Z==1 && N==1) // Quasi-Free process on Deuteron { G4LorentzVector f4Mom(0.,0.,0.,mNeut); // First neutron G4LorentzVector s4Mom(0.,0.,0.,mNeut); // Second neutron if(!G4QHadron(totLV).DecayIn2(f4Mom,s4Mom)) { G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: H2+(pi-)=>n+n"<push_back(neutr1); // Fill pi0 or Gamma in the output G4QHadron* neutr2 = new G4QHadron(2112,s4Mom); // Create Hadron for the 2nd Neutron output->push_back(neutr2); // Fill the neutron to the output } else { G4int rPDG=targPDG-1001; G4double mr=G4QPDGCode(rPDG).GetMass();// Mass of the residual Nucleus G4LorentzVector f4Mom(0.,0.,0.,mNeut); // First neutron G4LorentzVector s4Mom(0.,0.,0.,mNeut); // Second neutron G4LorentzVector r4Mom(0.,0.,0.,mr); // Residual nucleus if(!G4QHadron(totLV).DecayIn3(f4Mom,s4Mom,r4Mom)) { G4cerr<<"---Worning---G4QCaptureAtRest::AtRestDoIt: A+(pi-)=>n+n+(A-n-p)"<push_back(neutr1); // Fill the first neutron in the output G4QHadron* neutr2 = new G4QHadron(2112,s4Mom); // Create Hadron for the 2nd Neutron output->push_back(neutr2); // Fill the second neutron to the output G4QHadron* resnuc = new G4QHadron(rPDG,r4Mom); // Create Hadron for the ResidualNucl output->push_back(resnuc); // Fill the Residual Nucleus to the output } } else if((projPDG==13||projPDG==15) && !lepChan)//Normal BoundLepton->e+nu+anti_nu_e decay { G4double mbm=mp-EnergyDeposition; // Mass of the bounded muon G4LorentzVector totLV(0.,0.,0.,mbm); // 4-momentum of bounded muon/tau if(projPDG==13) // now for tau it is only energy deposition, for mu this EMCascade { MuCaptureEMCascade(Z, N, cascE); G4int nsec=cascE->size(); G4DynamicParticle* theSec = 0; // Prototype to fill particle in the G4ParticleChange for(G4int is=0; isoperator[](is); if(ener>0) theSec = new G4DynamicParticle(G4Electron::Electron(),RndmDir(),ener); else theSec = new G4DynamicParticle(G4Gamma::Gamma(),RndmDir(),-ener); totLV-=theSec->Get4Momentum(); G4Track* aNewTrack = new G4Track(theSec, localtime, position ); aNewTrack->SetWeight(weight); // weighted aNewTrack->SetTouchableHandle(trTouchable); cascT->push_back(aNewTrack); } #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: e+nu+nu decay 4M="<=mbm*(mmt-mbm/2)/mmt) { G4cout<<"-W-G4QCaptureAtRest::AtRestDoIt: Unrealistic E="<mu+A"<push_back(electron); // Fill the electron to the output G4QHadron* resnuc = new G4QHadron(targPDG,s4Mom); // Create Residual Nucleus output->push_back(resnuc); // Fill the Residual Nucleus to the output } else { G4double mr=std::sqrt(mmt*(mmt-ee-ee)+mEl2); // Mass of the M+nu+nu system if(mre+(A+2nu)"<A+NuMu+antiNuE"<.55) // Use mu decay instead of e-decay { deL=13; deN=-14; mdl=mMu; } G4LorentzVector e4Mom(0.,0.,0.,mdl); // mass of the electron G4LorentzVector n4Mom(0.,0.,0.,0.); // muon neutrino G4LorentzVector a4Mom(0.,0.,0.,0.); // electron anti-nutrino if(!G4QHadron(totLV).DecayIn3(e4Mom,n4Mom,a4Mom)) { G4cerr<<"--Worning--G4QCaptureAtRest::AtRestDoIt:Tau_b=>L+Nu_tau+anti_NuL"<n+nu_mu"<push_back(neutrino); // Fill pi0 or Gamma in the output G4QHadron* neutron = new G4QHadron(2112,n4Mom); // Create Hadron for the Neutron output->push_back(neutron); // Fill the neutron to the output } else if((projPDG==13||projPDG==15)&&lepChan&&targPDG==90001001)//LeptonCapture onDeuteron { G4LorentzVector totLV(0.,0.,0.,mp+mDeut-EnergyDeposition);// 4-mom of theCompoundSystem #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: CapOnDeutr decay 4M="<n+n+nu_mu"<push_back(neutrino); // Fill pi0 or Gamma in the output G4QHadron* neut1 = new G4QHadron(2112,n4Mom); // Create Hadron for the FirstNeutron output->push_back(neut1); // Fill the neutron to the output G4QHadron* neut2 = new G4QHadron(2112,s4Mom); // Create Hadron for the SecondNeutron output->push_back(neut2); // Fill the neutron to the output } // *** Closed *** else if((projPDG==13||projPDG==15)&&lepChan&&G4UniformRand()>1&&Z>0&&N>0)//@@QuasiFreeCap { G4double mt=G4QPDGCode(targPDG).GetMass();// Mass of the target Nucleus G4LorentzVector totLV(0.,0.,0.,mp+mt-EnergyDeposition);// 4-mom of the(A+mu-) compound #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: Quasi-Free decay 4M="<nu_mu+n+(A-p)"<push_back(neutrino); // Fill nutrino_mu in the output G4QHadron* neutron = new G4QHadron(2112,s4Mom);// Create Hadron for the 2nd Neutron output->push_back(neutron); // Fill the neutron to the output G4QHadron* resnuc = new G4QHadron(rPDG,r4Mom); // Create Hadron for the ResidualNucl output->push_back(resnuc); // Fill the Residual Nucleus to the output } else // *** THIS works for all particles ! *** //else if(1>2)// !! Immediately change back - Very temporary to avoid nuclear capture !! { G4QHadron* neutr = 0; // Create Neutrino if(projPDG==13||projPDG==15) mp-=EnergyDeposition;//TheEnergyDeposit is only for LepCap #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: CHIPS M="<size(); G4DynamicParticle* theSec = 0; // Prototype to fill particle in the G4ParticleChange for(G4int is=0; isoperator[](is); if(ener>0) theSec = new G4DynamicParticle(G4Electron::Electron(),RndmDir(),ener); else theSec = new G4DynamicParticle(G4Gamma::Gamma(),RndmDir(),-ener); projLV-=theSec->Get4Momentum(); G4Track* aNewTrack = new G4Track(theSec, localtime, position ); aNewTrack->SetWeight(weight); // weighted aNewTrack->SetTouchableHandle(trTouchable); cascT->push_back(aNewTrack); } // ----- Ericson mu to pi conversion ----- ????? ----- if(G4UniformRand()<.04) // .04 is a parameter ! { projPDG=-211; // Phase space decay of mu->nu+q+aq M=1 //G4LorentzVector f4Mom(0.,0.,0.,0.); // Muon neutrino //G4LorentzVector s4Mom(0.,0.,0.,0.); // Quark //G4LorentzVector r4Mom(0.,0.,0.,0.); // Anti-Quark //if(!G4QHadron(projLV).DecayIn3(f4Mom,s4Mom,r4Mom)) //{ // G4cerr<<"--Worning--G4QCaptureAtRest::AtRestDoIt: (mu-)=>q+aq+nu, M=1"<nu+q+aq with matrix element G4double mmu2=projLV.m2(); G4double mmu=std::sqrt(mmu2); G4double hmm=mmu/2.; G4double dmm=mmu+mmu; G4double qp=std::pow((std::pow(1.+ga*std::pow(hmm,b1),G4UniformRand())-1.)/ga,rb); G4double mqq=0.; if(qpq+aq+nu, M#1"<G4QCapAR:E/MCons, p="< DELETED---->---->----+ G4QHadronVector projHV; // | projHV.push_back(pH); // DESTROYED over 2 lines -+ | G4QEnvironment* pan= new G4QEnvironment(projHV,targPDG);// ---> DELETED --->-----+ | | std::for_each(projHV.begin(), projHV.end(), DeleteQHadron()); // <---<------<----+-+-+ projHV.clear(); // <------------<---------------<-------------------<------------+-+ #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: pPDG="< work for Hisaya @@ // The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body) G4QHadron* hadr=output->operator[](i); // Pointer to the output hadron if(hadr->GetNFragments()) // Intermediate hadron { #ifdef debug G4cout<<"G4QCaptureAtRest::AtRestDoIt: Intermediate particle is found i="<GetPDGCode(); #ifdef pdebug G4cout<<"G4QCaptureAtRest::AtRestDoIt:#"<.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S } else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus(); else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus(); else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus(); else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero(); else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus(); else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!) { G4int aZ = hadr->GetCharge(); G4int aA = hadr->GetBaryonNumber(); #ifdef pdebug G4cout<<"G4QCaptureAtRest::AtRestDoIt:Ion Z="<GetDefinition())==*(G4MuonMinus::MuonMinus()) || *(stoppedHadron->GetDefinition())==*(G4TauMinus::TauMinus()) ) return Time; else return 0.; } // Muon can decay or to be captured by the nucleus (Z,N): true=MuCapture, false=MuDecay G4bool G4QCaptureAtRest::RandomizeMuDecayOrCapture(G4int Z, G4int N) { static G4int mZ=0; // static memory about the last Z static G4int mN=0; // static memory about the last N static G4double mH=0.; // static memory about the last Huff(Z,N) static G4double mR=0.; // static memory about the last rate(Z,N) static const G4int nAZ=17; // total number of tabulated isotopes static const G4int nZm=10; // (maximumZ)+1 for which rates are tabulated static const G4int rin[nZm]={0,0,1,1,1,2,3,4,5,6}; // i=rin[Z]+N for the tabulated rate static const G4double rate[nAZ]={.00000045, .00000047, .00000215, .000000356, .00000468, .00000226, .00000610, .00002750, .000023500, .00003790, .00003500, .00006600, .00006200, .000102500, .00009500, .00008800, .00022900}; #ifdef debug G4cout<<"G4QCaptureAtRest::RandomizeMuDecayOrCapture is called"<9) Huff=1.-.000394*std::pow(z,2.19)/(1.+12.18*std::exp(z*.01373)); G4double pD=.00045516*Huff; // 1/MeanLifeTime of muon in atoms (in ns^-1)? G4double pC=1.e99; // Default 1/MeanLifeTime of muon NuclCapture(in ns^-1) // @@ Use Primakov correction (1-1.5625*N/(Z+N)) for isotopes. (M.K.) if(Z==mZ && N==mN) pC=mR; // Use already calculated value else if(Z>9) pC=.0000001256*std::pow(z,3.5)/(1.+.00429*std::exp(1.67*std::pow(z,.39))); else if(Z>0) pC=rate[rin[Z]+N]; // Tabulated light isotopes else G4cout<<"--Warning--G4QCaptureAtRest::RandomizeMuDecayOrCapture: NEG Z="<"<