[1199] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
[1337] | 26 | // $Id: G4QLowEnergy.cc,v 1.5 2010/06/14 16:11:27 mkossov Exp $ |
---|
[1340] | 27 | // GEANT4 tag $Name: hadr-chips-V09-03-08 $ |
---|
[1199] | 28 | // |
---|
| 29 | // ---------------- G4QLowEnergy class ----------------- |
---|
| 30 | // by Mikhail Kossov, Aug 2007. |
---|
| 31 | // G4QLowEnergy class of the CHIPS Simulation Branch in GEANT4 |
---|
| 32 | // --------------------------------------------------------------- |
---|
| 33 | // Short description: This is a fast low energy algorithm for the |
---|
| 34 | // inelastic interactions of nucleons and nuclei (ions) with nuclei. |
---|
| 35 | // This is a fase-space algorithm, but not quark level. Provides |
---|
| 36 | // nuclear fragments upto alpha only. Never was tumed (but can be). |
---|
| 37 | // --------------------------------------------------------------- |
---|
| 38 | |
---|
| 39 | //#define debug |
---|
| 40 | //#define pdebug |
---|
| 41 | //#define edebug |
---|
| 42 | //#define tdebug |
---|
| 43 | //#define nandebug |
---|
| 44 | //#define ppdebug |
---|
| 45 | |
---|
| 46 | #include "G4QLowEnergy.hh" |
---|
| 47 | |
---|
| 48 | // Initialization of static vectors |
---|
| 49 | G4int G4QLowEnergy::nPartCWorld=152; // #of particles initialized in CHIPS |
---|
| 50 | std::vector<G4int> G4QLowEnergy::ElementZ; // Z of element(i) in theLastCalc |
---|
| 51 | std::vector<G4double> G4QLowEnergy::ElProbInMat; // SumProbOfElem in Material |
---|
| 52 | std::vector<std::vector<G4int>*> G4QLowEnergy::ElIsoN;// N of isotope(j), E(i) |
---|
| 53 | std::vector<std::vector<G4double>*>G4QLowEnergy::IsoProbInEl;//SumProbIsotE(i) |
---|
| 54 | |
---|
| 55 | // Constructor |
---|
| 56 | G4QLowEnergy::G4QLowEnergy(const G4String& processName): |
---|
| 57 | G4VDiscreteProcess(processName, fHadronic), evaporate(true) |
---|
| 58 | { |
---|
| 59 | #ifdef debug |
---|
| 60 | G4cout<<"G4QLowEnergy::Constructor is called processName="<<processName<<G4endl; |
---|
| 61 | #endif |
---|
| 62 | if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created "<<G4endl; |
---|
| 63 | G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max) |
---|
| 64 | } |
---|
| 65 | |
---|
| 66 | // Destructor |
---|
| 67 | G4QLowEnergy::~G4QLowEnergy() {} |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation(){return EnMomConservation;} |
---|
| 71 | |
---|
| 72 | G4int G4QLowEnergy::GetNumberOfNeutronsInTarget() {return nOfNeutrons;} |
---|
| 73 | |
---|
| 74 | // output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)), |
---|
| 75 | // where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3) |
---|
| 76 | // ********** All CHIPS cross sections are calculated in the surface units ************ |
---|
| 77 | G4double G4QLowEnergy::GetMeanFreePath(const G4Track&Track, G4double, G4ForceCondition*F) |
---|
| 78 | { |
---|
| 79 | *F = NotForced; |
---|
| 80 | const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle(); |
---|
| 81 | G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition(); |
---|
| 82 | if( !IsApplicable(*incidentParticleDefinition)) |
---|
| 83 | G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath for notImplemented Particle"<<G4endl; |
---|
| 84 | // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material |
---|
| 85 | G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle |
---|
| 86 | #ifdef debug |
---|
| 87 | G4double KinEn = incidentParticle->GetKineticEnergy(); |
---|
| 88 | G4cout<<"G4QLowEnergy::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl; |
---|
| 89 | #endif |
---|
| 90 | const G4Material* material = Track.GetMaterial(); // Get the current material |
---|
| 91 | const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume(); |
---|
| 92 | const G4ElementVector* theElementVector = material->GetElementVector(); |
---|
| 93 | G4int nE=material->GetNumberOfElements(); |
---|
| 94 | #ifdef debug |
---|
| 95 | G4cout<<"G4QLowEnergy::GetMeanFreePath:"<<nE<<" Elems"<<G4endl; |
---|
| 96 | #endif |
---|
| 97 | G4int pPDG=0; |
---|
| 98 | G4int Z=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge()); |
---|
| 99 | G4int A=incidentParticleDefinition->GetBaryonNumber(); |
---|
| 100 | if ( incidentParticleDefinition == G4Proton::Proton() ) pPDG = 2212; |
---|
| 101 | else if ( incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020; |
---|
| 102 | else if ( incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 1000020040; |
---|
| 103 | else if ( incidentParticleDefinition == G4Triton::Triton() ) pPDG = 1000010030; |
---|
| 104 | else if ( incidentParticleDefinition == G4He3::He3() ) pPDG = 1000020030; |
---|
| 105 | else if ( incidentParticleDefinition == G4GenericIon::GenericIon() || (Z > 0 && A > 0)) |
---|
| 106 | { |
---|
| 107 | pPDG=incidentParticleDefinition->GetPDGEncoding(); |
---|
| 108 | #ifdef debug |
---|
| 109 | G4int prPDG=1000000000+10000*A+10*Z; |
---|
| 110 | G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl; |
---|
| 111 | #endif |
---|
| 112 | } |
---|
| 113 | else G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath: only AA & pA implemented"<<G4endl; |
---|
| 114 | G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer(); |
---|
| 115 | if(pPDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer(); // just to test |
---|
| 116 | Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A |
---|
| 117 | G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton |
---|
| 118 | G4double sigma=0.; // Sums over elements for the material |
---|
| 119 | G4int IPIE=IsoProbInEl.size(); // How many old elements? |
---|
| 120 | if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) |
---|
| 121 | { |
---|
| 122 | std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector |
---|
| 123 | SPI->clear(); |
---|
| 124 | delete SPI; |
---|
| 125 | std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector |
---|
| 126 | IsN->clear(); |
---|
| 127 | delete IsN; |
---|
| 128 | } |
---|
| 129 | ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) |
---|
| 130 | ElementZ.clear(); // Clear the body vector for Z of Elements |
---|
| 131 | IsoProbInEl.clear(); // Clear the body vector for SPI |
---|
| 132 | ElIsoN.clear(); // Clear the body vector for N of Isotopes |
---|
| 133 | for(G4int i=0; i<nE; ++i) |
---|
| 134 | { |
---|
| 135 | G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element |
---|
| 136 | G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element |
---|
| 137 | ElementZ.push_back(Z); // Remember Z of the Element |
---|
| 138 | G4int isoSize=0; // The default for the isoVectorLength is 0 |
---|
| 139 | G4int indEl=0; // Index of non-natural element or 0(default) |
---|
| 140 | G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect |
---|
| 141 | if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector |
---|
| 142 | #ifdef debug |
---|
| 143 | G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl; |
---|
| 144 | #endif |
---|
| 145 | if(isoSize) // The Element has non-trivial abundance set |
---|
| 146 | { |
---|
| 147 | indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order |
---|
| 148 | #ifdef debug |
---|
| 149 | G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl; |
---|
| 150 | #endif |
---|
| 151 | if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define |
---|
| 152 | { |
---|
| 153 | std::vector<std::pair<G4int,G4double>*>* newAbund = |
---|
| 154 | new std::vector<std::pair<G4int,G4double>*>; |
---|
| 155 | G4double* abuVector=pElement->GetRelativeAbundanceVector(); |
---|
| 156 | for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes |
---|
| 157 | { |
---|
| 158 | G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z ! |
---|
| 159 | if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z=" |
---|
| 160 | <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl; |
---|
| 161 | G4double abund=abuVector[j]; |
---|
| 162 | std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund); |
---|
| 163 | #ifdef debug |
---|
| 164 | G4cout<<"G4QLowEnergy::GetMeanFreePath:pair#"<<j<<",N="<<N<<",a="<<abund<<G4endl; |
---|
| 165 | #endif |
---|
| 166 | newAbund->push_back(pr); |
---|
| 167 | } |
---|
| 168 | #ifdef debug |
---|
| 169 | G4cout<<"G4QLowEnergy::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl; |
---|
| 170 | #endif |
---|
| 171 | indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd |
---|
| 172 | for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary |
---|
| 173 | delete newAbund; // Was "new" in the beginning of the name space |
---|
| 174 | } |
---|
| 175 | } |
---|
| 176 | std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer |
---|
| 177 | std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector |
---|
| 178 | IsoProbInEl.push_back(SPI); |
---|
| 179 | std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector |
---|
| 180 | ElIsoN.push_back(IsN); |
---|
| 181 | G4int nIs=cs->size(); // A#Of Isotopes in the Element |
---|
| 182 | #ifdef debug |
---|
| 183 | G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl; |
---|
| 184 | #endif |
---|
| 185 | G4double susi=0.; // sum of CS over isotopes |
---|
| 186 | if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El |
---|
| 187 | { |
---|
| 188 | std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice |
---|
| 189 | G4int N=curIs->first; // #of Neuterons in the isotope j of El i |
---|
| 190 | IsN->push_back(N); // Remember Min N for the Element |
---|
| 191 | #ifdef debug |
---|
| 192 | G4cout<<"G4QLoE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl; |
---|
| 193 | #endif |
---|
| 194 | G4bool ccsf=true; // Extract inelastic Ion-Ion cross-section |
---|
| 195 | #ifdef debug |
---|
| 196 | G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl; |
---|
| 197 | #endif |
---|
| 198 | G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope |
---|
| 199 | #ifdef debug |
---|
| 200 | G4cout<<"G4QLowEnergy::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom=" |
---|
| 201 | <<Momentum<<", XSec="<<CSI/millibarn<<G4endl; |
---|
| 202 | #endif |
---|
| 203 | curIs->second = CSI; |
---|
| 204 | susi+=CSI; // Make a sum per isotopes |
---|
| 205 | SPI->push_back(susi); // Remember summed cross-section |
---|
| 206 | } // End of temporary initialization of the cross sections in the G4QIsotope singeltone |
---|
| 207 | sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV) |
---|
| 208 | #ifdef debug |
---|
| 209 | G4cout<<"G4QLowEnergy::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl) |
---|
| 210 | <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl; |
---|
| 211 | #endif |
---|
| 212 | ElProbInMat.push_back(sigma); |
---|
| 213 | } // End of LOOP over Elements |
---|
| 214 | // Check that cross section is not zero and return the mean free path |
---|
| 215 | #ifdef debug |
---|
| 216 | G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl; |
---|
| 217 | #endif |
---|
| 218 | if(sigma > 0.000000001) return 1./sigma; // Mean path [distance] |
---|
| 219 | return DBL_MAX; |
---|
| 220 | } |
---|
| 221 | |
---|
| 222 | G4bool G4QLowEnergy::IsApplicable(const G4ParticleDefinition& particle) |
---|
| 223 | { |
---|
| 224 | G4int Z=static_cast<G4int>(particle.GetPDGCharge()); |
---|
| 225 | G4int A=particle.GetBaryonNumber(); |
---|
| 226 | if (particle == *( G4Proton::Proton() )) return true; |
---|
| 227 | else if (particle == *( G4Neutron::Neutron() )) return true; |
---|
| 228 | else if (particle == *( G4Deuteron::Deuteron() )) return true; |
---|
| 229 | else if (particle == *( G4Alpha::Alpha() )) return true; |
---|
| 230 | else if (particle == *( G4Triton::Triton() )) return true; |
---|
| 231 | else if (particle == *( G4He3::He3() )) return true; |
---|
| 232 | else if (particle == *( G4GenericIon::GenericIon() )) return true; |
---|
| 233 | else if (Z > 0 && A > 0) return true; |
---|
| 234 | #ifdef debug |
---|
| 235 | G4cout<<"***>>G4QLowEnergy::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<", A=" |
---|
| 236 | <<A<<", Z="<<Z<<G4endl; |
---|
| 237 | #endif |
---|
| 238 | return false; |
---|
| 239 | } |
---|
| 240 | |
---|
| 241 | G4VParticleChange* G4QLowEnergy::PostStepDoIt(const G4Track& track, const G4Step& step) |
---|
| 242 | { |
---|
| 243 | static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV; // CHIPS proton Mass in MeV |
---|
| 244 | static const G4double mPro2= mProt*mProt; // CHIPS sq proton Mass |
---|
| 245 | static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV; // CHIPS neutron Mass in MeV |
---|
| 246 | static const G4double mNeu2= mNeut*mNeut; // CHIPS sq neutron Mass |
---|
| 247 | static const G4double mLamb= G4QPDGCode(3122).GetMass()/MeV; // CHIPS Lambda Mass in MeV |
---|
| 248 | static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0)/MeV; |
---|
| 249 | static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0)/MeV; |
---|
| 250 | static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0)/MeV; |
---|
| 251 | static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0)/MeV; |
---|
| 252 | static const G4double mFm= 250*MeV; |
---|
| 253 | static const G4double third= 1./3.; |
---|
| 254 | static const G4ThreeVector zeroMom(0.,0.,0.); |
---|
| 255 | static G4ParticleDefinition* aGamma = G4Gamma::Gamma(); |
---|
[1337] | 256 | static G4ParticleDefinition* aPiZero = G4PionZero::PionZero(); |
---|
| 257 | static G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus(); |
---|
| 258 | static G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus(); |
---|
[1199] | 259 | static G4ParticleDefinition* aProton = G4Proton::Proton(); |
---|
| 260 | static G4ParticleDefinition* aNeutron = G4Neutron::Neutron(); |
---|
| 261 | static G4ParticleDefinition* aLambda = G4Lambda::Lambda(); |
---|
| 262 | static G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron(); |
---|
| 263 | static G4ParticleDefinition* aTriton = G4Triton::Triton(); |
---|
| 264 | static G4ParticleDefinition* aHe3 = G4He3::He3(); |
---|
| 265 | static G4ParticleDefinition* anAlpha = G4Alpha::Alpha(); |
---|
| 266 | static const G4int nCh=26; // #of combinations |
---|
| 267 | static G4QNucleus Nuc; // A fake nucleus to call Evaporation |
---|
| 268 | // |
---|
| 269 | //------------------------------------------------------------------------------------- |
---|
| 270 | static G4bool CWinit = true; // CHIPS Warld needs to be initted |
---|
| 271 | if(CWinit) |
---|
| 272 | { |
---|
| 273 | CWinit=false; |
---|
| 274 | G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) |
---|
| 275 | } |
---|
| 276 | //------------------------------------------------------------------------------------- |
---|
| 277 | const G4DynamicParticle* projHadron = track.GetDynamicParticle(); |
---|
| 278 | const G4ParticleDefinition* particle=projHadron->GetDefinition(); |
---|
| 279 | #ifdef pdebug |
---|
| 280 | G4cout<<"G4QLowEnergy::PostStepDoIt: *** Called *** In4M=" |
---|
| 281 | <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type=" |
---|
| 282 | <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl; |
---|
| 283 | #endif |
---|
| 284 | //G4ForceCondition cond=NotForced; |
---|
| 285 | //GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters? |
---|
| 286 | #ifdef debug |
---|
| 287 | G4cout<<"G4QLowEnergy::PostStepDoIt: After GetMeanFreePath is called"<<G4endl; |
---|
| 288 | #endif |
---|
| 289 | std::vector<G4Track*> result; |
---|
| 290 | G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV! |
---|
| 291 | G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV |
---|
| 292 | G4double Momentum = proj4M.rho(); // @@ Just for the test purposes |
---|
| 293 | if(std::fabs(Momentum-momentum)>.000001) |
---|
| 294 | G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl; |
---|
| 295 | #ifdef debug |
---|
| 296 | G4cout<<"G4QLowEnergy::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",proj4M,m=" |
---|
| 297 | <<proj4M<<proj4M.m()<<G4endl; |
---|
| 298 | #endif |
---|
| 299 | if (!IsApplicable(*particle)) // Check applicability |
---|
| 300 | { |
---|
| 301 | G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl; |
---|
| 302 | return 0; |
---|
| 303 | } |
---|
| 304 | const G4Material* material = track.GetMaterial(); // Get the current material |
---|
| 305 | const G4ElementVector* theElementVector = material->GetElementVector(); |
---|
| 306 | G4int nE=material->GetNumberOfElements(); |
---|
| 307 | #ifdef debug |
---|
| 308 | G4cout<<"G4QLowEnergy::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl; |
---|
| 309 | #endif |
---|
| 310 | G4int projPDG=0; // PDG Code prototype for the captured hadron |
---|
| 311 | // Not all these particles are implemented yet (see Is Applicable) |
---|
| 312 | G4int Z=static_cast<G4int>(particle->GetPDGCharge()); |
---|
| 313 | G4int A=particle->GetBaryonNumber(); |
---|
| 314 | if (particle == G4Proton::Proton() ) projPDG= 2212; |
---|
| 315 | else if (particle == G4Neutron::Neutron() ) projPDG= 2112; |
---|
| 316 | else if (particle == G4Deuteron::Deuteron() ) projPDG= 1000010020; |
---|
| 317 | else if (particle == G4Alpha::Alpha() ) projPDG= 1000020040; |
---|
| 318 | else if (particle == G4Triton::Triton() ) projPDG= 1000010030; |
---|
| 319 | else if (particle == G4He3::He3() ) projPDG= 1000020030; |
---|
| 320 | else if (particle == G4GenericIon::GenericIon() || (Z > 0 && A > 0)) |
---|
| 321 | { |
---|
| 322 | projPDG=particle->GetPDGEncoding(); |
---|
| 323 | #ifdef debug |
---|
| 324 | G4int prPDG=1000000000+10000*Z+10*A; |
---|
| 325 | G4cout<<"G4QLowEnergy::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl; |
---|
| 326 | #endif |
---|
| 327 | } |
---|
| 328 | else G4cout<<"-Warning-G4QLowEnergy::PostStepDoIt:Unknown projectile Ion"<<G4endl; |
---|
| 329 | #ifdef pdebug |
---|
| 330 | G4int prPDG=particle->GetPDGEncoding(); |
---|
| 331 | G4cout<<"G4QLowEnergy::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl; |
---|
| 332 | #endif |
---|
| 333 | if(!projPDG) |
---|
| 334 | { |
---|
| 335 | G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl; |
---|
| 336 | return 0; |
---|
| 337 | } |
---|
| 338 | // Element treatment |
---|
| 339 | G4int EPIM=ElProbInMat.size(); |
---|
| 340 | #ifdef debug |
---|
| 341 | G4cout<<"G4QLowEn::PostStDoIt: m="<<EPIM<<", n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl; |
---|
| 342 | #endif |
---|
| 343 | G4int i=0; |
---|
| 344 | if(EPIM>1) |
---|
| 345 | { |
---|
| 346 | G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand(); |
---|
| 347 | for(i=0; i<nE; ++i) |
---|
| 348 | { |
---|
| 349 | #ifdef debug |
---|
| 350 | G4cout<<"G4QLowEn::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl; |
---|
| 351 | #endif |
---|
| 352 | if (rnd<ElProbInMat[i]) break; |
---|
| 353 | } |
---|
| 354 | if(i>=nE) i=nE-1; // Top limit for the Element |
---|
| 355 | } |
---|
| 356 | G4Element* pElement=(*theElementVector)[i]; |
---|
| 357 | G4int tZ=static_cast<G4int>(pElement->GetZ()); |
---|
| 358 | #ifdef debug |
---|
| 359 | G4cout<<"G4QLowEnergy::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl; |
---|
| 360 | #endif |
---|
| 361 | if(tZ<=0) |
---|
| 362 | { |
---|
| 363 | G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Element with Z="<<tZ<<G4endl; |
---|
| 364 | if(tZ<0) return 0; |
---|
| 365 | } |
---|
| 366 | std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes |
---|
| 367 | std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i] |
---|
| 368 | G4int nofIsot=SPI->size(); // #of isotopes in the element i |
---|
| 369 | #ifdef debug |
---|
| 370 | G4cout<<"G4QLowEnergy::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl; |
---|
| 371 | #endif |
---|
| 372 | G4int j=0; |
---|
| 373 | if(nofIsot>1) |
---|
| 374 | { |
---|
| 375 | G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element |
---|
| 376 | for(j=0; j<nofIsot; ++j) |
---|
| 377 | { |
---|
| 378 | #ifdef debug |
---|
| 379 | G4cout<<"G4QLowEnergy::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl; |
---|
| 380 | #endif |
---|
| 381 | if(rndI < (*SPI)[j]) break; |
---|
| 382 | } |
---|
| 383 | if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope |
---|
| 384 | } |
---|
| 385 | G4int tN =(*IsN)[j]; ; // Randomized number of neutrons |
---|
| 386 | #ifdef debug |
---|
| 387 | G4cout<<"G4QLowEnergy::PostStepDoIt: j="<<i<<", N(isotope)="<<tN<<", MeV="<<MeV<<G4endl; |
---|
| 388 | #endif |
---|
| 389 | if(tN<0) |
---|
| 390 | { |
---|
| 391 | G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Isotope Z="<<tZ<<" has 0>N="<<tN<<G4endl; |
---|
| 392 | return 0; |
---|
| 393 | } |
---|
| 394 | nOfNeutrons=tN; // Remember it for the energy-momentum check |
---|
| 395 | #ifdef debug |
---|
| 396 | G4cout<<"G4QLowEnergy::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl; |
---|
| 397 | #endif |
---|
| 398 | if(tN<0) |
---|
| 399 | { |
---|
| 400 | G4cerr<<"***Warning***G4QLowEnergy::PostStepDoIt: Element with N="<<tN<< G4endl; |
---|
| 401 | return 0; |
---|
| 402 | } |
---|
| 403 | aParticleChange.Initialize(track); |
---|
| 404 | #ifdef debug |
---|
| 405 | G4cout<<"G4QLowEnergy::PostStepDoIt: track is initialized"<<G4endl; |
---|
| 406 | #endif |
---|
| 407 | G4double weight = track.GetWeight(); |
---|
| 408 | G4double localtime = track.GetGlobalTime(); |
---|
| 409 | G4ThreeVector position = track.GetPosition(); |
---|
| 410 | #ifdef debug |
---|
| 411 | G4cout<<"G4QLowEnergy::PostStepDoIt: before Touchable extraction"<<G4endl; |
---|
| 412 | #endif |
---|
| 413 | G4TouchableHandle trTouchable = track.GetTouchableHandle(); |
---|
| 414 | #ifdef debug |
---|
| 415 | G4cout<<"G4QLowEnergy::PostStepDoIt: Touchable is extracted"<<G4endl; |
---|
| 416 | #endif |
---|
| 417 | G4QPDGCode targQPDG(90000000+tZ*1000+tN); // @@ use G4Ion and get rid of CHIPS World |
---|
| 418 | G4double tM=targQPDG.GetMass(); // CHIPS target nucleus mass in MeV |
---|
| 419 | G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3); // Strangeness |
---|
| 420 | G4int pZ=static_cast<G4int>(particle->GetPDGCharge()); // Charge of the projectile |
---|
| 421 | G4int pN=particle->GetBaryonNumber()-pZ-pL;// #of neutrons in projectile |
---|
| 422 | G4double pM=targQPDG.GetNuclMass(pZ,pN,0); // CHIPS projectile nucleus mass in MeV |
---|
| 423 | G4double cosp=-14*Momentum*(tM-pM)/tM/pM; // Asymmetry power for angular distribution |
---|
| 424 | #ifdef debug |
---|
| 425 | G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ |
---|
| 426 | <<","<<tN<<"), cosp="<<cosp<<G4endl; |
---|
| 427 | #endif |
---|
| 428 | G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?) |
---|
| 429 | G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector |
---|
| 430 | G4LorentzVector targ4M=G4LorentzVector(0.,0.,0.,tM); // Target's 4-mom |
---|
| 431 | G4LorentzVector tot4M =proj4M+targ4M; // Total 4-mom of the reaction |
---|
| 432 | #ifdef pdebug |
---|
| 433 | G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl; |
---|
| 434 | #endif |
---|
| 435 | EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation |
---|
| 436 | // @@ Probably this is not necessary any more |
---|
| 437 | #ifdef debug |
---|
| 438 | G4cout<<"G4QLE::PSDI:false,P="<<Momentum<<",Z="<<pZ<<",N="<<pN<<",PDG="<<projPDG<<G4endl; |
---|
| 439 | #endif |
---|
| 440 | G4double xSec=CalculateXS(Momentum, tZ, tN, projPDG); // Recalculate CrossSection |
---|
| 441 | #ifdef pdebug |
---|
| 442 | G4cout<<"G4QLowEn::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",tZ="<<tZ<<",N="<<tN<<",XS=" |
---|
| 443 | <<xSec/millibarn<<G4endl; |
---|
| 444 | #endif |
---|
| 445 | #ifdef nandebug |
---|
| 446 | if(xSec>0. || xSec<0. || xSec==0); |
---|
| 447 | else G4cout<<"-Warning-G4QLowEnergy::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl; |
---|
| 448 | #endif |
---|
| 449 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
| 450 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
| 451 | { |
---|
| 452 | #ifdef debug |
---|
| 453 | G4cerr<<"-Warning-G4QLowEnergy::PSDoIt:*Zero cross-section* PDG="<<projPDG |
---|
| 454 | <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl; |
---|
| 455 | #endif |
---|
| 456 | //Do Nothing Action insead of the reaction |
---|
| 457 | aParticleChange.ProposeEnergy(kinEnergy); |
---|
| 458 | aParticleChange.ProposeLocalEnergyDeposit(0.); |
---|
| 459 | aParticleChange.ProposeMomentumDirection(dir) ; |
---|
| 460 | return G4VDiscreteProcess::PostStepDoIt(track,step); |
---|
| 461 | } |
---|
| 462 | // Kill interacting hadron |
---|
| 463 | aParticleChange.ProposeEnergy(0.); |
---|
| 464 | aParticleChange.ProposeTrackStatus(fStopAndKill); |
---|
| 465 | G4int tB=tZ+tN; |
---|
| 466 | G4int pB=pZ+pN; |
---|
| 467 | #ifdef pdebug |
---|
| 468 | G4cout<<"G4QLowEn::PSDI: Projectile track is killed"<<", tA="<<tB<<", pA="<<pB<<G4endl; |
---|
| 469 | #endif |
---|
| 470 | // algorithm implementation STARTS HERE --- All calculations are in IU -------- |
---|
| 471 | G4double tA=tB; |
---|
| 472 | G4double pA=pB; |
---|
| 473 | G4double tR=1.1; // target nucleus R in fm |
---|
| 474 | if(tB > 1) tR*=std::pow(tA,third); // in fm |
---|
| 475 | G4double pR=1.1; // projectile nucleus R in fm |
---|
| 476 | if(pB > 1) pR*=std::pow(pA,third); // in fm |
---|
| 477 | G4double R=tR+pR; // total radius |
---|
| 478 | G4double R2=R*R; |
---|
| 479 | G4int tD=0; |
---|
| 480 | G4int pD=0; |
---|
| 481 | G4int nAt=0; |
---|
| 482 | G4int nAtM=27; |
---|
| 483 | G4int nSec = 0; |
---|
| 484 | G4double tcM=0.; |
---|
| 485 | G4double tnM=1.; |
---|
| 486 | #ifdef edebug |
---|
| 487 | G4int totChg=0; |
---|
| 488 | G4int totBaN=0; |
---|
| 489 | G4LorentzVector tch4M =tot4M; // Total 4-mom of the reaction |
---|
| 490 | #endif |
---|
| 491 | while((!tD && !pD && ++nAt<nAtM) || tcM<tnM) |
---|
| 492 | { |
---|
| 493 | #ifdef edebug |
---|
| 494 | totChg=tZ+pZ; |
---|
| 495 | totBaN=tB+pB; |
---|
| 496 | tch4M =tot4M; |
---|
| 497 | G4cout<<">G4QLEn::PSDI:#"<<nAt<<",tC="<<totChg<<",tA="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 498 | #endif |
---|
| 499 | G4LorentzVector tt4M=tot4M; |
---|
| 500 | G4int resN=result.size(); |
---|
| 501 | if(resN) |
---|
| 502 | { |
---|
| 503 | for(G4int i=0; i<resN; ++i) delete result[i]; |
---|
| 504 | result.clear(); |
---|
| 505 | } |
---|
| 506 | G4double D=std::sqrt(R2*G4UniformRand()); |
---|
| 507 | #ifdef pdebug |
---|
| 508 | G4cout<<"G4QLowEn::PSDI: D="<<D<<", tR="<<tR<<", pR="<<pR<<G4endl; |
---|
| 509 | #endif |
---|
| 510 | if(D > std::fabs(tR-pR)) // leading parts can be separated |
---|
| 511 | { |
---|
| 512 | nSec = 0; |
---|
| 513 | G4double DtR=D-tR; |
---|
| 514 | G4double DpR=D-pR; |
---|
| 515 | G4double DtR2=DtR*DtR; |
---|
| 516 | G4double DpR2=DpR*DpR; |
---|
| 517 | //G4double DtR3=DtR2*DtR; |
---|
| 518 | G4double DpR3=DpR2*DpR; |
---|
| 519 | //G4double DtR4=DtR3*DtR; |
---|
| 520 | G4double DpR4=DpR3*DpR; |
---|
| 521 | G4double tR2=tR*tR; |
---|
| 522 | G4double pR2=pR*pR; |
---|
| 523 | G4double tR3=tR2*tR; |
---|
| 524 | //G4double pR3=pR2*pR; |
---|
| 525 | G4double tR4=tR3*tR; |
---|
| 526 | //G4double pR4=pR3*pR; |
---|
| 527 | G4double HD=16.*D; |
---|
| 528 | G4double tF=tA*(6*tR2*(pR2-DtR2)-4*D*(tR3-DpR3)+3*(tR4-DpR4))/HD/tR3; |
---|
| 529 | G4double pF=tF; |
---|
| 530 | tD=static_cast<G4int>(tF); |
---|
| 531 | pD=static_cast<G4int>(pF); |
---|
| 532 | //if(G4UniformRand() < tF-tD) ++tD; // Simple solution |
---|
| 533 | //if(G4UniformRand() < pF-pD) ++pD; |
---|
| 534 | // Enhance alphas solution |
---|
| 535 | if(std::fabs(tF-4.)<1.) tD=4; // @@ 1. is the enhancement parameter |
---|
| 536 | else if(G4UniformRand() < tF-tD) ++tD; |
---|
| 537 | if(std::fabs(pF-4.)<1.) pD=4; |
---|
| 538 | else if(G4UniformRand() < pF-pD) ++pD; |
---|
| 539 | if(tD > tB) tD=tB; |
---|
| 540 | if(pD > pB) pD=tB; |
---|
| 541 | // @@ Quasi Free is not debugged @@ The following close it |
---|
| 542 | if(tD < 1) tD=1; |
---|
| 543 | if(pD < 1) pD=1; |
---|
| 544 | #ifdef pdebug |
---|
| 545 | G4cout<<"G4QLE::PSDI:#"<<nAt<<",pF="<<pF<<",tF="<<tF<<",pD="<<pD<<",tD="<<tD<<G4endl; |
---|
| 546 | #endif |
---|
| 547 | G4int pC=0; // charge of the projectile fraction |
---|
| 548 | G4int tC=0; // charge of the target fraction |
---|
| 549 | if((tD || pD) && tD<tB && pD<pB)// Periferal interaction |
---|
| 550 | { |
---|
| 551 | if(!tD || !pD) // Quasi-Elastic: B+A->(B-1)+N+A || ->B+N+(A-1) |
---|
| 552 | { |
---|
[1315] | 553 | G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); |
---|
| 554 | G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); |
---|
[1199] | 555 | G4int pPDG=2112; // Proto of the nucleon PDG (proton) |
---|
| 556 | G4double prM =mNeut; // Proto of the nucleon mass |
---|
| 557 | G4double prM2=mNeu2; // Proto of the nucleon sq mass |
---|
| 558 | if (!tD) // Quasi-elastic scattering of the proj QF nucleon |
---|
| 559 | { |
---|
| 560 | if(pD > 1) pD=1; |
---|
| 561 | if(!pN || (pZ && pA*G4UniformRand() < pZ) ) // proj QF proton |
---|
| 562 | { |
---|
| 563 | pPDG=2212; |
---|
| 564 | prM=mProt; |
---|
| 565 | prM2=mPro2; |
---|
| 566 | pC=1; |
---|
| 567 | } |
---|
| 568 | G4double Mom=0.; |
---|
| 569 | G4LorentzVector com4M = targ4M; // Proto of 4mom of compound |
---|
| 570 | G4double tgM=targ4M.e(); |
---|
| 571 | G4double rNM=0.; |
---|
| 572 | G4LorentzVector rNuc4M(0.,0.,0.,0.); |
---|
| 573 | if(pD==pB) |
---|
| 574 | { |
---|
| 575 | Mom=proj4M.rho(); |
---|
| 576 | com4M += proj4M; // Total 4-mom for scattering |
---|
| 577 | rNM=prM; |
---|
| 578 | } |
---|
| 579 | else // It is necessary to split the nucleon (with fermiM) |
---|
| 580 | { |
---|
| 581 | G4ThreeVector fm=mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); |
---|
| 582 | rNM=G4QPDGCode(2112).GetNuclMass(pZ-pC, pN, 0); |
---|
| 583 | G4double rNE=std::sqrt(fm*fm+rNM*rNM); |
---|
| 584 | rNuc4M=G4LorentzVector(fm,rNE); |
---|
| 585 | G4ThreeVector boostV=proj4M.vect()/proj4M.e(); |
---|
| 586 | rNuc4M.boost(boostV); |
---|
| 587 | G4LorentzVector qfN4M=proj4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS |
---|
| 588 | com4M += qfN4M; // Calculate Total 4Mom for NA scattering |
---|
| 589 | G4double pNE = qfN4M.e(); // Energy of the QF nucleon in LS |
---|
| 590 | if(rNE >= prM) Mom = std::sqrt(pNE*pNE-prM2); // Mom(s) fake value |
---|
| 591 | else break; // Break the while loop |
---|
| 592 | } |
---|
[1315] | 593 | G4double xSec=0.; |
---|
| 594 | if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); |
---|
| 595 | else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); |
---|
[1199] | 596 | if( xSec <= 0. ) break; // Break the while loop |
---|
[1315] | 597 | G4double mint=0.; // Prototype of functional randomized -t |
---|
| 598 | if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t |
---|
| 599 | else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t |
---|
| 600 | G4double maxt=0.; // Prototype of maximum -t |
---|
| 601 | if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t |
---|
| 602 | else maxt=NCSmanager->GetHMaxT(); // maximum -t |
---|
| 603 | G4double cost=1.-mint/maxt; // cos(theta) in CMS |
---|
[1199] | 604 | if(cost>1. || cost<-1.) break; // Break the while loop |
---|
| 605 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tgM); // 4mom of recoil target |
---|
| 606 | G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N |
---|
| 607 | G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01); |
---|
| 608 | if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) |
---|
| 609 | { |
---|
| 610 | G4cout<<"G4QLE::Pt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl; |
---|
| 611 | break; // Break the while loop |
---|
| 612 | } |
---|
| 613 | G4Track* projSpect = 0; |
---|
| 614 | G4Track* aNucTrack = 0; |
---|
| 615 | if(pB > pD) // Fill the proj residual nucleus |
---|
| 616 | { |
---|
| 617 | G4int rZ=pZ-pC; |
---|
| 618 | G4int rA=pB-1; |
---|
| 619 | G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition |
---|
| 620 | if(rA==1) |
---|
| 621 | { |
---|
| 622 | if(!rZ) theDefinition = aNeutron; |
---|
| 623 | else theDefinition = aProton; |
---|
| 624 | } |
---|
| 625 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ); |
---|
| 626 | G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M); |
---|
| 627 | projSpect = new G4Track(resN, localtime, position ); |
---|
| 628 | projSpect->SetWeight(weight); // weighted |
---|
| 629 | projSpect->SetTouchableHandle(trTouchable); |
---|
| 630 | #ifdef pdebug |
---|
| 631 | G4cout<<"G4QLowEn::PSDI:-->ProjQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl; |
---|
| 632 | #endif |
---|
| 633 | ++nSec; |
---|
| 634 | } |
---|
| 635 | G4ParticleDefinition* theDefinition = G4Neutron::Neutron(); |
---|
| 636 | if(pPDG==2212) theDefinition = G4Proton::Proton(); |
---|
| 637 | G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M); |
---|
| 638 | aNucTrack = new G4Track(scatN, localtime, position ); |
---|
| 639 | aNucTrack->SetWeight(weight); // weighted |
---|
| 640 | aNucTrack->SetTouchableHandle(trTouchable); |
---|
| 641 | #ifdef pdebug |
---|
| 642 | G4cout<<"G4QLowEn::PSDI:-->TgQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl; |
---|
| 643 | #endif |
---|
| 644 | ++nSec; |
---|
| 645 | G4Track* aFraTrack=0; |
---|
| 646 | if(tA==1) |
---|
| 647 | { |
---|
| 648 | if(!tZ) theDefinition = aNeutron; |
---|
| 649 | else theDefinition = aProton; |
---|
| 650 | } |
---|
| 651 | else if(tA==8 && tC==4) // Be8 decay |
---|
| 652 | { |
---|
| 653 | theDefinition = anAlpha; |
---|
| 654 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha |
---|
| 655 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha |
---|
| 656 | if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) |
---|
| 657 | { |
---|
| 658 | G4cout<<"*G4QLE::TS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; |
---|
| 659 | } |
---|
| 660 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 661 | aFraTrack = new G4Track(pAl, localtime, position ); |
---|
| 662 | aFraTrack->SetWeight(weight); // weighted |
---|
| 663 | aFraTrack->SetTouchableHandle(trTouchable); |
---|
| 664 | #ifdef pdebug |
---|
| 665 | G4cout<<"G4QLowEn::PSDI:-->TgRQFA4M="<<f4M<<G4endl; |
---|
| 666 | #endif |
---|
| 667 | ++nSec; |
---|
| 668 | reco4M=s4M; |
---|
| 669 | } |
---|
| 670 | else if(tA==5 && (tC==2 || tC==3)) // He5/Li5 decay |
---|
| 671 | { |
---|
| 672 | theDefinition = aProton; |
---|
| 673 | G4double mNuc = mProt; |
---|
| 674 | if(tC==2) |
---|
| 675 | { |
---|
| 676 | theDefinition = aNeutron; |
---|
| 677 | mNuc = mNeut; |
---|
| 678 | } |
---|
| 679 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon |
---|
| 680 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha |
---|
| 681 | if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) |
---|
| 682 | { |
---|
| 683 | G4cout<<"*G4QLE::TS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; |
---|
| 684 | } |
---|
| 685 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 686 | aFraTrack = new G4Track(pAl, localtime, position ); |
---|
| 687 | aFraTrack->SetWeight(weight); // weighted |
---|
| 688 | aFraTrack->SetTouchableHandle(trTouchable); |
---|
| 689 | #ifdef pdebug |
---|
| 690 | G4cout<<"G4QLowEn::PSDI:-->TgQFRN4M="<<f4M<<G4endl; |
---|
| 691 | #endif |
---|
| 692 | ++nSec; |
---|
| 693 | theDefinition = anAlpha; |
---|
| 694 | reco4M=s4M; |
---|
| 695 | } |
---|
| 696 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(tZ,tB,0,tZ); |
---|
| 697 | ++nSec; |
---|
| 698 | #ifdef pdebug |
---|
| 699 | G4cout<<"G4QLowEn::PSDI:-->PQF_nSec="<<nSec<<G4endl; |
---|
| 700 | #endif |
---|
| 701 | aParticleChange.SetNumberOfSecondaries(nSec); |
---|
| 702 | if(projSpect) aParticleChange.AddSecondary(projSpect); |
---|
| 703 | if(aNucTrack) aParticleChange.AddSecondary(aNucTrack); |
---|
| 704 | if(aFraTrack) aParticleChange.AddSecondary(aFraTrack); |
---|
| 705 | G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M); |
---|
| 706 | G4Track* aResTrack = new G4Track(resA, localtime, position ); |
---|
| 707 | aResTrack->SetWeight(weight); // weighted |
---|
| 708 | aResTrack->SetTouchableHandle(trTouchable); |
---|
| 709 | #ifdef pdebug |
---|
| 710 | G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl; |
---|
| 711 | #endif |
---|
| 712 | aParticleChange.AddSecondary(aResTrack); |
---|
| 713 | } |
---|
| 714 | else // !pD : QF target Nucleon on the whole Projectile |
---|
| 715 | { |
---|
| 716 | if(tD > 1) tD=1; |
---|
| 717 | if(!tN || (tZ && tA*G4UniformRand() < tZ) ) // target QF proton |
---|
| 718 | { |
---|
| 719 | pPDG=2212; |
---|
| 720 | prM=mProt; |
---|
| 721 | prM2=mPro2; |
---|
| 722 | tC=1; |
---|
| 723 | } |
---|
| 724 | G4double Mom=0.; |
---|
| 725 | G4LorentzVector com4M=proj4M; // Proto of 4mom of compound |
---|
| 726 | G4double prM=proj4M.m(); |
---|
| 727 | G4double rNM=0.; |
---|
| 728 | G4LorentzVector rNuc4M(0.,0.,0.,0.); |
---|
| 729 | if(tD==tB) |
---|
| 730 | { |
---|
| 731 | Mom=prM*proj4M.rho()/proj4M.m(); |
---|
| 732 | com4M += targ4M; // Total 4-mom for scattering |
---|
| 733 | rNM=prM; |
---|
| 734 | } |
---|
| 735 | else // It is necessary to split the nucleon (with fermiM) |
---|
| 736 | { |
---|
| 737 | G4ThreeVector fm=250.*std::pow(G4UniformRand(),third)*G4RandomDirection(); |
---|
| 738 | rNM=G4QPDGCode(2112).GetNuclMass(tZ-tC, tN, 0)/MeV; |
---|
| 739 | G4double rNE=std::sqrt(fm*fm+rNM*rNM); |
---|
| 740 | rNuc4M=G4LorentzVector(fm,rNE); |
---|
| 741 | G4LorentzVector qfN4M=targ4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS |
---|
[1315] | 742 | com4M += qfN4M; // Calculate Total 4Mom for NA scattering |
---|
[1199] | 743 | G4ThreeVector boostV=proj4M.vect()/proj4M.e(); |
---|
| 744 | qfN4M.boost(-boostV); |
---|
[1315] | 745 | G4double tNE = qfN4M.e(); // Energy of the QF nucleon in LS |
---|
[1199] | 746 | if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2); // Mom(s) fake value |
---|
[1315] | 747 | else break; // Break the while loop |
---|
[1199] | 748 | } |
---|
[1315] | 749 | G4double xSec=0.; |
---|
| 750 | if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); |
---|
| 751 | else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG); |
---|
| 752 | if( xSec <= 0. ) break; // Break the while loop |
---|
| 753 | G4double mint=0.; // Prototype of functional randomized -t |
---|
| 754 | if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t |
---|
| 755 | else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t |
---|
| 756 | G4double maxt=0.; // Prototype of maximum -t |
---|
| 757 | if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t |
---|
| 758 | else maxt=NCSmanager->GetHMaxT(); // maximum -t |
---|
| 759 | G4double cost=1.-mint/maxt; // cos(theta) in CMS |
---|
| 760 | if(cost>1. || cost<-1.) break; // Break the while loop |
---|
[1199] | 761 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,prM); // 4mom of recoil target |
---|
| 762 | G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N |
---|
| 763 | G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01); |
---|
| 764 | if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) |
---|
| 765 | { |
---|
| 766 | G4cout<<"G4QLE::Tt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl; |
---|
[1315] | 767 | break; // Break the while loop |
---|
[1199] | 768 | } |
---|
| 769 | G4Track* targSpect = 0; |
---|
| 770 | G4Track* aNucTrack = 0; |
---|
[1315] | 771 | if(tB > tD) // Fill the residual nucleus |
---|
[1199] | 772 | { |
---|
| 773 | G4int rZ=tZ-tC; |
---|
| 774 | G4int rA=tB-1; |
---|
| 775 | G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition |
---|
| 776 | if(rA==1) |
---|
| 777 | { |
---|
| 778 | if(!rZ) theDefinition = aNeutron; |
---|
| 779 | else theDefinition = aProton; |
---|
| 780 | } |
---|
| 781 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ); |
---|
| 782 | G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M); |
---|
| 783 | targSpect = new G4Track(resN, localtime, position ); |
---|
| 784 | targSpect->SetWeight(weight); // weighted |
---|
| 785 | targSpect->SetTouchableHandle(trTouchable); |
---|
| 786 | #ifdef pdebug |
---|
| 787 | G4cout<<"G4QLowEn::PSDI:-->TargQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl; |
---|
| 788 | #endif |
---|
| 789 | ++nSec; |
---|
| 790 | } |
---|
| 791 | G4ParticleDefinition* theDefinition = G4Neutron::Neutron(); |
---|
| 792 | if(pPDG==2212) theDefinition = G4Proton::Proton(); |
---|
| 793 | G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M); |
---|
| 794 | aNucTrack = new G4Track(scatN, localtime, position ); |
---|
| 795 | aNucTrack->SetWeight(weight); // weighted |
---|
| 796 | aNucTrack->SetTouchableHandle(trTouchable); |
---|
| 797 | #ifdef pdebug |
---|
| 798 | G4cout<<"G4QLowEn::PSDI:-->PrQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl; |
---|
| 799 | #endif |
---|
| 800 | ++nSec; |
---|
| 801 | G4Track* aFraTrack=0; |
---|
| 802 | if(pA==1) |
---|
| 803 | { |
---|
| 804 | if(!pZ) theDefinition = aNeutron; |
---|
| 805 | else theDefinition = aProton; |
---|
| 806 | } |
---|
| 807 | else if(pA==8 && pC==4) // Be8 decay |
---|
| 808 | { |
---|
| 809 | theDefinition = anAlpha; |
---|
| 810 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha |
---|
| 811 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha |
---|
| 812 | if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) |
---|
| 813 | { |
---|
| 814 | G4cout<<"*G4QLE::PS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; |
---|
| 815 | } |
---|
| 816 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 817 | aFraTrack = new G4Track(pAl, localtime, position ); |
---|
| 818 | aFraTrack->SetWeight(weight); // weighted |
---|
| 819 | aFraTrack->SetTouchableHandle(trTouchable); |
---|
| 820 | #ifdef pdebug |
---|
| 821 | G4cout<<"G4QLowEn::PSDI:-->PrRQFA4M="<<f4M<<G4endl; |
---|
| 822 | #endif |
---|
| 823 | ++nSec; |
---|
| 824 | reco4M=s4M; |
---|
| 825 | } |
---|
| 826 | else if(pA==5 && (pC==2 || pC==3)) // He5/Li5 decay |
---|
| 827 | { |
---|
| 828 | theDefinition = aProton; |
---|
| 829 | G4double mNuc = mProt; |
---|
| 830 | if(pC==2) |
---|
| 831 | { |
---|
| 832 | theDefinition = aNeutron; |
---|
| 833 | mNuc = mNeut; |
---|
| 834 | } |
---|
| 835 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon |
---|
| 836 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha |
---|
| 837 | if(!G4QHadron(reco4M).DecayIn2(f4M, s4M)) |
---|
| 838 | { |
---|
| 839 | G4cout<<"*G4QLE::PS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; |
---|
| 840 | } |
---|
| 841 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 842 | aFraTrack = new G4Track(pAl, localtime, position ); |
---|
| 843 | aFraTrack->SetWeight(weight); // weighted |
---|
| 844 | aFraTrack->SetTouchableHandle(trTouchable); |
---|
| 845 | #ifdef pdebug |
---|
| 846 | G4cout<<"G4QLowEn::PSDI:-->PrQFRN4M="<<f4M<<G4endl; |
---|
| 847 | #endif |
---|
| 848 | ++nSec; |
---|
| 849 | theDefinition = anAlpha; |
---|
| 850 | reco4M=s4M; |
---|
| 851 | } |
---|
| 852 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(pZ,pB,0,pZ); |
---|
| 853 | ++nSec; |
---|
| 854 | #ifdef pdebug |
---|
| 855 | G4cout<<"G4QLowEn::PSDI:-->TQF_nSec="<<nSec<<G4endl; |
---|
| 856 | #endif |
---|
| 857 | aParticleChange.SetNumberOfSecondaries(nSec); |
---|
| 858 | if(targSpect) aParticleChange.AddSecondary(targSpect); |
---|
| 859 | if(aNucTrack) aParticleChange.AddSecondary(aNucTrack); |
---|
| 860 | if(aFraTrack) aParticleChange.AddSecondary(aFraTrack); |
---|
| 861 | G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M); |
---|
| 862 | G4Track* aResTrack = new G4Track(resA, localtime, position ); |
---|
| 863 | aResTrack->SetWeight(weight); // weighted |
---|
| 864 | aResTrack->SetTouchableHandle(trTouchable); |
---|
| 865 | #ifdef pdebug |
---|
| 866 | G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl; |
---|
| 867 | #endif |
---|
| 868 | aParticleChange.AddSecondary( aResTrack ); |
---|
| 869 | } |
---|
| 870 | #ifdef debug |
---|
| 871 | G4cout<<"G4QLowEnergy::PostStepDoIt:***PostStepDoIt is done:Quasi-El***"<<G4endl; |
---|
| 872 | #endif |
---|
| 873 | return G4VDiscreteProcess::PostStepDoIt(track, step); |
---|
| 874 | } |
---|
| 875 | else // The cental region compound can be created |
---|
| 876 | { |
---|
| 877 | // First calculate the isotopic state of the parts of the compound |
---|
| 878 | if(!pZ) pC=0; |
---|
| 879 | else if(!pN) pC=pD; |
---|
| 880 | else |
---|
| 881 | { |
---|
| 882 | #ifdef pdebug |
---|
| 883 | G4cout<<"G4QLowEn::PSDI: pD="<<pD<<", pZ="<<pZ<<", pA="<<pA<<G4endl; |
---|
| 884 | #endif |
---|
| 885 | G4double C=pD*pZ/pA; |
---|
| 886 | pC=static_cast<G4int>(C); |
---|
| 887 | if(G4UniformRand() < C-pC) ++pC; |
---|
| 888 | } |
---|
| 889 | if(!tZ) tC=0; |
---|
| 890 | else if(!tN) tC=tD; |
---|
| 891 | else |
---|
| 892 | { |
---|
| 893 | #ifdef pdebug |
---|
| 894 | G4cout<<"G4QLowEn::PSDI: tD="<<tD<<", tZ="<<tZ<<", tA="<<tA<<G4endl; |
---|
| 895 | #endif |
---|
| 896 | G4double C=tD*tZ/tA; |
---|
| 897 | tC=static_cast<G4int>(C); |
---|
| 898 | if(G4UniformRand() < C-tC) ++tC; |
---|
| 899 | } |
---|
| 900 | // calculate the transferred momentum |
---|
| 901 | G4ThreeVector pFM(0.,0.,0.); |
---|
| 902 | if(pD<pB) // The projectile nucleus must be splitted |
---|
| 903 | { |
---|
| 904 | G4int nc=pD; |
---|
| 905 | if(pD+pD>pB) nc=pB-pD; |
---|
| 906 | pFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); |
---|
| 907 | for(G4int i=1; i < nc; ++i) |
---|
| 908 | pFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); |
---|
| 909 | } |
---|
| 910 | G4ThreeVector tFM(0.,0.,0.); |
---|
| 911 | if(tD<tB) // The projectile nucleus must be splitted |
---|
| 912 | { |
---|
| 913 | G4int nc=pD; |
---|
| 914 | if(tD+tD>tB) nc=tB-tD; |
---|
| 915 | tFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); |
---|
| 916 | for(G4int i=1; i < nc; ++i) |
---|
| 917 | tFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection(); |
---|
| 918 | } |
---|
| 919 | #ifdef pdebug |
---|
| 920 | G4cout<<"G4QLE::PSDI:pC="<<pC<<", tC="<<tC<<", pFM="<<pFM<<", tFM="<<tFM<<G4endl; |
---|
| 921 | #endif |
---|
| 922 | // Split the projectile spectator |
---|
| 923 | G4int rpZ=pZ-pC; // number of protons in the projectile spectator |
---|
| 924 | G4int pF=pD-pC; // number of neutrons in the projectile part of comp |
---|
| 925 | G4int rpN=pN-pF; // number of neutrons in the projectile spectator |
---|
| 926 | G4double rpNM=G4QPDGCode(2112).GetNuclMass(rpZ, rpN, 0); // Mass of the spectator |
---|
| 927 | G4ThreeVector boostV=proj4M.vect()/proj4M.e(); // Antilab Boost Vector |
---|
| 928 | G4double rpE=std::sqrt(rpNM*rpNM+pFM.mag2()); |
---|
| 929 | G4LorentzVector rp4M(pFM,rpE); |
---|
| 930 | #ifdef pdebug |
---|
| 931 | G4cout<<"G4QLE::PSDI: boostV="<<boostV<<",rp4M="<<rp4M<<",pr4M="<<proj4M<<G4endl; |
---|
| 932 | #endif |
---|
| 933 | rp4M.boost(boostV); |
---|
| 934 | #ifdef pdebug |
---|
| 935 | G4cout<<"G4QLE::PSDI: After boost, rp4M="<<rp4M<<G4endl; |
---|
| 936 | #endif |
---|
| 937 | G4ParticleDefinition* theDefinition; // Prototype of projSpectatorNuclDefinition |
---|
| 938 | G4int rpA=rpZ+rpN; |
---|
| 939 | G4Track* aFraPTrack = 0; |
---|
| 940 | theDefinition = 0; |
---|
| 941 | if(rpA==1) |
---|
| 942 | { |
---|
| 943 | if(!rpZ) theDefinition = G4Neutron::Neutron(); |
---|
| 944 | else theDefinition = G4Proton::Proton(); |
---|
| 945 | #ifdef pdebug |
---|
| 946 | G4cout<<"G4QLE::PSDI: rpA=1, rpZ"<<rpZ<<G4endl; |
---|
| 947 | #endif |
---|
| 948 | } |
---|
[1315] | 949 | else if(rpA==2 && rpZ==0) // nn decay |
---|
| 950 | { |
---|
| 951 | theDefinition = aNeutron; |
---|
| 952 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron |
---|
| 953 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron |
---|
| 954 | #ifdef pdebug |
---|
| 955 | G4cout<<"G4QLE::CPS->n+n,nn="<<rp4M.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl; |
---|
| 956 | #endif |
---|
| 957 | if(!G4QHadron(rp4M).DecayIn2(f4M, s4M)) |
---|
| 958 | { |
---|
| 959 | G4cout<<"*W*G4QLE::CPS->n+n,t="<<rp4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl; |
---|
| 960 | } |
---|
| 961 | G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M); |
---|
| 962 | aFraPTrack = new G4Track(pNeu, localtime, position ); |
---|
| 963 | aFraPTrack->SetWeight(weight); // weighted |
---|
| 964 | aFraPTrack->SetTouchableHandle(trTouchable); |
---|
| 965 | tt4M-=f4M; |
---|
| 966 | #ifdef edebug |
---|
| 967 | totBaN-=2; |
---|
| 968 | tch4M -=f4M; |
---|
| 969 | G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 970 | #endif |
---|
| 971 | #ifdef pdebug |
---|
| 972 | G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; |
---|
| 973 | #endif |
---|
| 974 | ++nSec; |
---|
| 975 | rp4M=s4M; |
---|
| 976 | } |
---|
| 977 | else if(rpA>2 && rpZ==0) // Z=0 decay |
---|
| 978 | { |
---|
| 979 | theDefinition = aNeutron; |
---|
| 980 | G4LorentzVector f4M=rp4M/rpA; // 4mom of 1st neutron |
---|
| 981 | #ifdef pdebug |
---|
| 982 | G4cout<<"G4QLE::CPS->Nn,M="<<rp4M.m()<<" >? N*MNeutron="<<rpA*mNeutron<<G4endl; |
---|
| 983 | #endif |
---|
| 984 | for(G4int it=1; it<rpA; ++it) // Fill (N-1) neutrons to output |
---|
| 985 | { |
---|
| 986 | G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M); |
---|
| 987 | G4Track* aNTrack = new G4Track(pNeu, localtime, position ); |
---|
| 988 | aNTrack->SetWeight(weight); // weighted |
---|
| 989 | aNTrack->SetTouchableHandle(trTouchable); |
---|
| 990 | result.push_back(aNTrack); |
---|
| 991 | } |
---|
| 992 | G4int nesc = rpA-1; |
---|
| 993 | tt4M-=f4M*nesc; |
---|
| 994 | #ifdef edebug |
---|
| 995 | totBaN-=nesc; |
---|
| 996 | tch4M -=f4M*nesc; |
---|
| 997 | G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 998 | #endif |
---|
| 999 | #ifdef pdebug |
---|
| 1000 | G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; |
---|
| 1001 | #endif |
---|
| 1002 | nSec+=nesc; |
---|
| 1003 | rp4M=f4M; |
---|
| 1004 | } |
---|
[1199] | 1005 | else if(rpA==8 && rpZ==4) // Be8 decay |
---|
| 1006 | { |
---|
| 1007 | theDefinition = anAlpha; |
---|
| 1008 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha |
---|
| 1009 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha |
---|
| 1010 | #ifdef pdebug |
---|
| 1011 | G4cout<<"G4QLE::CPS->A+A,mBe8="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; |
---|
| 1012 | #endif |
---|
| 1013 | if(!G4QHadron(rp4M).DecayIn2(f4M, s4M)) |
---|
| 1014 | { |
---|
| 1015 | G4cout<<"*W*G4QLE::CPS->A+A,t="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; |
---|
| 1016 | } |
---|
| 1017 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 1018 | aFraPTrack = new G4Track(pAl, localtime, position ); |
---|
| 1019 | aFraPTrack->SetWeight(weight); // weighted |
---|
| 1020 | aFraPTrack->SetTouchableHandle(trTouchable); |
---|
| 1021 | tt4M-=f4M; |
---|
| 1022 | #ifdef edebug |
---|
| 1023 | totChg-=2; |
---|
| 1024 | totBaN-=4; |
---|
| 1025 | tch4M -=f4M; |
---|
| 1026 | G4cout<<">>G4QLEn::PSDI:1,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 1027 | #endif |
---|
| 1028 | #ifdef pdebug |
---|
| 1029 | G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl; |
---|
| 1030 | #endif |
---|
| 1031 | ++nSec; |
---|
| 1032 | rp4M=s4M; |
---|
| 1033 | } |
---|
| 1034 | else if(rpA==5 && (rpZ==2 || rpZ==3)) // He5/Li5 decay |
---|
| 1035 | { |
---|
| 1036 | theDefinition = aProton; |
---|
| 1037 | G4double mNuc = mProt; |
---|
| 1038 | if(rpZ==2) |
---|
| 1039 | { |
---|
| 1040 | theDefinition = aNeutron; |
---|
| 1041 | mNuc = mNeut; |
---|
| 1042 | } |
---|
| 1043 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon |
---|
| 1044 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha |
---|
| 1045 | #ifdef pdebug |
---|
| 1046 | G4cout<<"G4QLowE::CPS->N+A, tM5="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; |
---|
| 1047 | #endif |
---|
| 1048 | if(!G4QHadron(rp4M).DecayIn2(f4M, s4M)) |
---|
| 1049 | { |
---|
| 1050 | G4cout<<"*W*G4QLE::CPS->N+A,t="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; |
---|
| 1051 | } |
---|
| 1052 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 1053 | aFraPTrack = new G4Track(pAl, localtime, position ); |
---|
| 1054 | aFraPTrack->SetWeight(weight); // weighted |
---|
| 1055 | aFraPTrack->SetTouchableHandle(trTouchable); |
---|
| 1056 | tt4M-=f4M; |
---|
| 1057 | #ifdef edebug |
---|
| 1058 | if(theDefinition == aProton) totChg-=1; |
---|
| 1059 | totBaN-=1; |
---|
| 1060 | tch4M -=f4M; |
---|
| 1061 | G4cout<<">>G4QLEn::PSDI:2,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 1062 | #endif |
---|
| 1063 | #ifdef pdebug |
---|
| 1064 | G4cout<<"G4QLowEn::PSDI:-->ProjSpectN4M="<<f4M<<G4endl; |
---|
| 1065 | #endif |
---|
| 1066 | ++nSec; |
---|
| 1067 | theDefinition = anAlpha; |
---|
| 1068 | rp4M=s4M; |
---|
| 1069 | } |
---|
| 1070 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rpZ,rpA,0,rpZ); |
---|
| 1071 | if(!theDefinition) |
---|
| 1072 | G4cout<<"-Warning-G4QLowEn::PSDI: pDef=0, Z="<<rpZ<<", A="<<rpA<<G4endl; |
---|
| 1073 | #ifdef edebug |
---|
| 1074 | if(theDefinition == anAlpha) |
---|
| 1075 | { |
---|
| 1076 | totChg-=2; |
---|
| 1077 | totBaN-=4; |
---|
| 1078 | } |
---|
| 1079 | else |
---|
| 1080 | { |
---|
| 1081 | totChg-=rpZ; |
---|
| 1082 | totBaN-=rpA; |
---|
| 1083 | } |
---|
| 1084 | tch4M -=rp4M; |
---|
| 1085 | G4cout<<">>G4QLEn::PSDI:3, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl; |
---|
| 1086 | #endif |
---|
| 1087 | G4DynamicParticle* rpD = new G4DynamicParticle(theDefinition, rp4M); |
---|
| 1088 | G4Track* aNewPTrack = new G4Track(rpD, localtime, position); |
---|
| 1089 | aNewPTrack->SetWeight(weight);// weighted |
---|
| 1090 | aNewPTrack->SetTouchableHandle(trTouchable); |
---|
| 1091 | tt4M-=rp4M; |
---|
| 1092 | #ifdef pdebug |
---|
| 1093 | G4cout<<"G4QLowEn::PSDI:-->ProjSpectR4M="<<rp4M<<",Z="<<rpZ<<",A="<<rpA<<G4endl; |
---|
| 1094 | #endif |
---|
| 1095 | ++nSec; |
---|
| 1096 | // |
---|
| 1097 | // Split the target spectator |
---|
| 1098 | G4int rtZ=tZ-tC; // number of protons in the target spectator |
---|
| 1099 | G4int tF=tD-tC; // number of neutrons in the target part of comp |
---|
| 1100 | G4int rtN=tN-tF; // number of neutrons in the target spectator |
---|
| 1101 | G4double rtNM=G4QPDGCode(2112).GetNuclMass(rtZ, rtN, 0); // Mass of the spectator |
---|
| 1102 | G4double rtE=std::sqrt(rtNM*rtNM+tFM.mag2()); |
---|
| 1103 | G4LorentzVector rt4M(tFM,rtE); |
---|
| 1104 | G4int rtA=rtZ+rtN; |
---|
| 1105 | G4Track* aFraTTrack = 0; |
---|
| 1106 | theDefinition = 0; |
---|
| 1107 | if(rtA==1) |
---|
| 1108 | { |
---|
| 1109 | if(!rtZ) theDefinition = G4Neutron::Neutron(); |
---|
| 1110 | else theDefinition = G4Proton::Proton(); |
---|
| 1111 | #ifdef pdebug |
---|
| 1112 | G4cout<<"G4QLE::PSDI: rtA=1, rtZ"<<rtZ<<G4endl; |
---|
| 1113 | #endif |
---|
| 1114 | } |
---|
| 1115 | else if(rtA==8 && rtZ==4) // Be8 decay |
---|
| 1116 | { |
---|
| 1117 | theDefinition = anAlpha; |
---|
| 1118 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha |
---|
| 1119 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha |
---|
| 1120 | #ifdef pdebug |
---|
| 1121 | G4cout<<"G4QLE::CPS->A+A,mBe8="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; |
---|
| 1122 | #endif |
---|
| 1123 | if(!G4QHadron(rt4M).DecayIn2(f4M, s4M)) |
---|
| 1124 | { |
---|
| 1125 | G4cout<<"*W*G4QLE::CTS->A+A,t="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl; |
---|
| 1126 | } |
---|
| 1127 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 1128 | aFraTTrack = new G4Track(pAl, localtime, position ); |
---|
| 1129 | aFraTTrack->SetWeight(weight); // weighted |
---|
| 1130 | aFraTTrack->SetTouchableHandle(trTouchable); |
---|
| 1131 | tt4M-=f4M; |
---|
| 1132 | #ifdef edebug |
---|
| 1133 | totChg-=2; |
---|
| 1134 | totBaN-=4; |
---|
| 1135 | tch4M -=f4M; |
---|
| 1136 | G4cout<<">>G4QLEn::PSDI:4,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 1137 | #endif |
---|
| 1138 | #ifdef pdebug |
---|
| 1139 | G4cout<<"G4QLowEn::PSDI:-->TargSpectA4M="<<f4M<<G4endl; |
---|
| 1140 | #endif |
---|
| 1141 | ++nSec; |
---|
| 1142 | rt4M=s4M; |
---|
| 1143 | } |
---|
| 1144 | else if(rtA==5 && (rtZ==2 || rtZ==3)) // He5/Li5 decay |
---|
| 1145 | { |
---|
| 1146 | theDefinition = aProton; |
---|
| 1147 | G4double mNuc = mProt; |
---|
| 1148 | if(rpZ==2) |
---|
| 1149 | { |
---|
| 1150 | theDefinition = aNeutron; |
---|
| 1151 | mNuc = mNeut; |
---|
| 1152 | } |
---|
| 1153 | G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon |
---|
| 1154 | G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha |
---|
| 1155 | #ifdef pdebug |
---|
| 1156 | G4cout<<"G4QLowE::CPS->N+A, tM5="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; |
---|
| 1157 | #endif |
---|
| 1158 | if(!G4QHadron(rt4M).DecayIn2(f4M, s4M)) |
---|
| 1159 | { |
---|
| 1160 | G4cout<<"*W*G4QLE::CTS->N+A,t="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl; |
---|
| 1161 | } |
---|
| 1162 | G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M); |
---|
| 1163 | aFraTTrack = new G4Track(pAl, localtime, position ); |
---|
| 1164 | aFraTTrack->SetWeight(weight); // weighted |
---|
| 1165 | aFraTTrack->SetTouchableHandle(trTouchable); |
---|
| 1166 | tt4M-=f4M; |
---|
| 1167 | #ifdef edebug |
---|
| 1168 | if(theDefinition == aProton) totChg-=1; |
---|
| 1169 | totBaN-=1; |
---|
| 1170 | tch4M -=f4M; |
---|
| 1171 | G4cout<<">>G4QLEn::PSDI:5,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl; |
---|
| 1172 | #endif |
---|
| 1173 | #ifdef pdebug |
---|
| 1174 | G4cout<<"G4QLowEn::PSDI:-->TargSpectN4M="<<f4M<<G4endl; |
---|
| 1175 | #endif |
---|
| 1176 | ++nSec; |
---|
| 1177 | theDefinition = anAlpha; |
---|
| 1178 | rt4M=s4M; |
---|
| 1179 | } |
---|
| 1180 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rtZ,rtA,0,rtZ); |
---|
| 1181 | if(!theDefinition) |
---|
| 1182 | G4cout<<"-Warning-G4QLowEn::PSDI: tDef=0, Z="<<rtZ<<", A="<<rtA<<G4endl; |
---|
| 1183 | #ifdef edebug |
---|
| 1184 | if(theDefinition == anAlpha) |
---|
| 1185 | { |
---|
| 1186 | totChg-=2; |
---|
| 1187 | totBaN-=4; |
---|
| 1188 | } |
---|
| 1189 | else |
---|
| 1190 | { |
---|
| 1191 | totChg-=rtZ; |
---|
| 1192 | totBaN-=rtA; |
---|
| 1193 | } |
---|
| 1194 | tch4M -=rt4M; |
---|
| 1195 | G4cout<<">>G4QLEn::PSDI:6, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl; |
---|
| 1196 | #endif |
---|
| 1197 | G4DynamicParticle* rtD = new G4DynamicParticle(theDefinition, rt4M); |
---|
| 1198 | G4Track* aNewTTrack = new G4Track(rtD, localtime, position); |
---|
| 1199 | aNewTTrack->SetWeight(weight); // weighted |
---|
| 1200 | aNewTTrack->SetTouchableHandle(trTouchable); |
---|
| 1201 | tt4M-=rt4M; |
---|
| 1202 | #ifdef pdebug |
---|
| 1203 | G4cout<<"G4QLowEn::PSDI:-->TargSpectR4M="<<rt4M<<",Z="<<rtZ<<",A="<<rtA<<G4endl; |
---|
| 1204 | #endif |
---|
| 1205 | ++nSec; |
---|
| 1206 | nSec+=3; |
---|
| 1207 | aParticleChange.SetNumberOfSecondaries(nSec); |
---|
| 1208 | if(aFraPTrack) result.push_back(aFraPTrack); |
---|
| 1209 | if(aNewPTrack) result.push_back(aNewPTrack); |
---|
| 1210 | if(aFraTTrack) result.push_back(aFraTTrack); |
---|
| 1211 | if(aNewTTrack) result.push_back(aNewTTrack); |
---|
| 1212 | tcM = tt4M.m(); // total CMS mass of the compound |
---|
| 1213 | G4int sN=tF+pF; |
---|
| 1214 | G4int sZ=tC+pC; |
---|
| 1215 | tnM = targQPDG.GetNuclMass(sZ,sN,0); // GSM |
---|
| 1216 | #ifdef pdebug |
---|
| 1217 | G4cout<<"G4QLEn::PSDI:At#"<<nAt<<",totM="<<tcM<<",gsM="<<tnM<<",Z="<<sZ<<",N=" |
---|
| 1218 | <<sN<<G4endl; |
---|
| 1219 | #endif |
---|
| 1220 | if(tcM>tnM) |
---|
| 1221 | { |
---|
| 1222 | pZ=pC; |
---|
| 1223 | pN=pF; |
---|
| 1224 | tZ=tC; |
---|
| 1225 | tN=tF; |
---|
| 1226 | tot4M=tt4M; |
---|
| 1227 | } |
---|
| 1228 | } |
---|
| 1229 | } // At least one is splitted |
---|
| 1230 | else if(tD==tB || pD==pB) // Total absorption |
---|
| 1231 | { |
---|
| 1232 | tD=tZ+tN; |
---|
| 1233 | pD=pZ+pN; |
---|
| 1234 | tcM=tnM+1.; |
---|
| 1235 | } |
---|
| 1236 | } |
---|
| 1237 | else // Total compound (define tD to go out of the while |
---|
| 1238 | { |
---|
| 1239 | tD=tZ+tN; |
---|
| 1240 | pD=pZ+pN; |
---|
| 1241 | tcM=tnM+1.; |
---|
| 1242 | } |
---|
| 1243 | } // End of the interaction WHILE |
---|
| 1244 | G4double totM=tot4M.m(); // total CMS mass of the reaction |
---|
| 1245 | G4int totN=tN+pN; |
---|
| 1246 | G4int totZ=tZ+pZ; |
---|
| 1247 | #ifdef edebug |
---|
| 1248 | G4cout<<">>>G4QLEn::PSDI: dZ="<<totChg-totZ<<", dB="<<totBaN-totN-totZ<<", d4M=" |
---|
| 1249 | <<tch4M-tot4M<<",N="<<totN<<",Z="<<totZ<<G4endl; |
---|
| 1250 | #endif |
---|
| 1251 | // @@ Here mass[i] can be calculated if mass=0 |
---|
| 1252 | G4double mass[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., |
---|
| 1253 | 0.,0.,0.,0.,0.,0.}; |
---|
| 1254 | mass[0] = tM = targQPDG.GetNuclMass(totZ,totN,0); // gamma+gamma |
---|
| 1255 | #ifdef pdebug |
---|
| 1256 | G4cout<<"G4QLEn::PSDI:TotM="<<totM<<",NucM="<<tM<<",totZ="<<totZ<<",totN="<<totN<<G4endl; |
---|
| 1257 | #endif |
---|
| 1258 | if (totN>0 && totZ>0) |
---|
| 1259 | { |
---|
| 1260 | mass[1] = targQPDG.GetNuclMass(totZ,totN-1,0); // gamma+neutron |
---|
| 1261 | mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0); // gamma+proton |
---|
| 1262 | } |
---|
| 1263 | if ( totZ > 1 && totN > 1 ) mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0); //g+d |
---|
| 1264 | if ( totZ > 1 && totN > 2 ) mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0); //g+t |
---|
| 1265 | if ( totZ > 2 && totN > 1 ) mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0); //g+3 |
---|
| 1266 | if ( totZ > 2 && totN > 2 ) mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0); //g+a |
---|
| 1267 | if ( totZ > 0 && totN > 2 ) mass[7] = targQPDG.GetNuclMass(totZ ,totN-2,0); //n+n |
---|
| 1268 | mass[ 8] = mass[3]; // neutron+proton (the same as a deuteron) |
---|
| 1269 | if ( totZ > 2 ) mass[9] = targQPDG.GetNuclMass(totZ-2,totN ,0); //p+p |
---|
| 1270 | mass[10] = mass[5]; // proton+deuteron (the same as He3) |
---|
| 1271 | mass[11] = mass[4]; // neutron+deuteron (the same as t) |
---|
| 1272 | mass[12] = mass[6]; // deuteron+deuteron (the same as alpha) |
---|
| 1273 | mass[13] = mass[6]; // proton+tritium (the same as alpha) |
---|
| 1274 | if ( totZ > 1 && totN > 3 ) mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0);//n+t |
---|
| 1275 | if ( totZ > 3 && totN > 1 ) mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0);//He3+p |
---|
| 1276 | mass[16] = mass[6]; // neutron+He3 (the same as alpha) |
---|
| 1277 | if ( totZ > 3 && totN > 2 ) mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0);//pa |
---|
| 1278 | if ( totZ > 2 && totN > 3 ) mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0);//na |
---|
| 1279 | if(pL>0) // @@ Not debugged @@ |
---|
| 1280 | { |
---|
| 1281 | G4int pL1=pL-1; |
---|
| 1282 | if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ ,totN ,pL1);// Lambda+gamma |
---|
| 1283 | if( (totN > 0 && totZ > 0) || totZ > 1 ) |
---|
| 1284 | mass[20]=targQPDG.GetNuclMass(totZ-1,totN ,pL1);//Lp |
---|
| 1285 | if( (totN > 0 && totZ > 0) || totN > 0 ) |
---|
| 1286 | mass[21]=targQPDG.GetNuclMass(totZ ,totN-1,pL1);//Ln |
---|
| 1287 | if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) ) |
---|
| 1288 | mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);//Ld |
---|
| 1289 | if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) ) |
---|
| 1290 | mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);//Lt |
---|
| 1291 | if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) ) |
---|
| 1292 | mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);//L3 |
---|
| 1293 | if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) ) |
---|
| 1294 | mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);//La |
---|
| 1295 | } |
---|
| 1296 | #ifdef debug |
---|
| 1297 | G4cout<<"G4QLowEn::PSDI: Residual masses are calculated"<<G4endl; |
---|
| 1298 | #endif |
---|
| 1299 | tA=tZ+tN; |
---|
| 1300 | pA=pZ+pN; |
---|
| 1301 | G4double prZ=pZ/pA+tZ/tA; |
---|
| 1302 | G4double prN=pN/pA+tN/tA; |
---|
| 1303 | G4double prD=prN*prZ; |
---|
| 1304 | G4double prA=prD*prD; |
---|
| 1305 | G4double prH=prD*prZ; |
---|
| 1306 | G4double prT=prD*prN; |
---|
| 1307 | G4double fhe3=6.*std::sqrt(tA); |
---|
| 1308 | G4double prL=0.; |
---|
| 1309 | if(pL>0) prL=pL/pA; |
---|
| 1310 | G4double qval[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., |
---|
| 1311 | 0.,0.,0.,0.,0.,0.}; |
---|
| 1312 | qval[ 0] = (totM - mass[ 0])/137./137.; |
---|
| 1313 | qval[ 1] = (totM - mass[ 1] - mNeut)*prN/137.; |
---|
| 1314 | qval[ 2] = (totM - mass[ 2] - mProt)*prZ/137.; |
---|
| 1315 | qval[ 3] = (totM - mass[ 3] - mDeut)*prD/3./137.; |
---|
| 1316 | qval[ 4] = (totM - mass[ 4] - mTrit)*prT/6./137.; |
---|
| 1317 | qval[ 5] = (totM - mass[ 5] - mHel3)*prH/fhe3/137.; |
---|
| 1318 | qval[ 6] = (totM - mass[ 6] - mAlph)*prA/9./137.; |
---|
| 1319 | qval[ 7] = (totM - mass[ 7] - mNeut - mNeut)*prN*prN; |
---|
| 1320 | qval[ 8] = (totM - mass[ 8] - mNeut - mProt)*prD; |
---|
| 1321 | qval[ 9] = (totM - mass[ 9] - mProt - mProt)*prZ*prZ; |
---|
| 1322 | qval[10] = (totM - mass[10] - mProt - mDeut)*prH/3.; |
---|
| 1323 | qval[11] = (totM - mass[11] - mNeut - mDeut)*prT/3.; |
---|
| 1324 | qval[12] = (totM - mass[12] - mDeut - mDeut)*prA/3./3.; |
---|
| 1325 | qval[13] = (totM - mass[13] - mProt - mTrit)*prA/6.; |
---|
| 1326 | qval[14] = (totM - mass[14] - mNeut - mTrit)*prT*prN/6.; |
---|
| 1327 | qval[15] = (totM - mass[15] - mProt - mHel3)*prH*prZ/fhe3; |
---|
| 1328 | qval[16] = (totM - mass[16] - mNeut - mHel3)*prA/fhe3; |
---|
| 1329 | qval[17] = (totM - mass[17] - mProt - mAlph)*prZ*prA/9.; |
---|
| 1330 | qval[18] = (totM - mass[18] - mNeut - mAlph)*prN*prA/9.; |
---|
| 1331 | if(pZ>0) |
---|
| 1332 | { |
---|
| 1333 | qval[19] = (totM - mass[19] - mLamb)*prL; |
---|
| 1334 | qval[20] = (totM - mass[20] - mProt - mLamb)*prL*prZ; |
---|
| 1335 | qval[21] = (totM - mass[21] - mNeut - mLamb)*prL*prN; |
---|
| 1336 | qval[22] = (totM - mass[22] - mDeut - mLamb)*prL*prD/2.; |
---|
| 1337 | qval[23] = (totM - mass[23] - mTrit - mLamb)*prL*prT/3.; |
---|
| 1338 | qval[24] = (totM - mass[24] - mHel3 - mLamb)*prL*prH/fhe3; |
---|
| 1339 | qval[25] = (totM - mass[25] - mAlph - mLamb)*prL*prA/4; |
---|
| 1340 | } |
---|
| 1341 | #ifdef debug |
---|
| 1342 | G4cout<<"G4QLowEn::PSDI: Q-values are calculated, tgA="<<tA<<", prA="<<pA<<G4endl; |
---|
| 1343 | #endif |
---|
| 1344 | |
---|
| 1345 | G4double qv = 0.0; // Total sum of probabilities (q-values) |
---|
| 1346 | for(G4int i=0; i<nCh; ++i ) |
---|
| 1347 | { |
---|
| 1348 | #ifdef sdebug |
---|
| 1349 | G4cout<<"G4QLowEn::PSDI: i="<<i<<", q="<<qval[i]<<",m="<<mass[i]<<G4endl; |
---|
| 1350 | #endif |
---|
| 1351 | if( mass[i] < 500.*MeV ) qval[i] = 0.0; // Close A/Z impossible channels |
---|
| 1352 | if( qval[i] < 0.0 ) qval[i] = 0.0; // Close the splitting impossible channels |
---|
| 1353 | qv += qval[i]; |
---|
| 1354 | } |
---|
| 1355 | // Select the channel |
---|
| 1356 | G4double qv1 = 0.0; |
---|
| 1357 | G4double ran = G4UniformRand(); |
---|
| 1358 | G4int index = 0; |
---|
| 1359 | for( index=0; index<nCh; ++index ) |
---|
| 1360 | { |
---|
| 1361 | if( qval[index] > 0.0 ) |
---|
| 1362 | { |
---|
| 1363 | qv1 += qval[index]/qv; |
---|
| 1364 | if( ran <= qv1 ) break; |
---|
| 1365 | } |
---|
| 1366 | } |
---|
| 1367 | #ifdef debug |
---|
| 1368 | G4cout<<"G4QLowEn::PSDI: index="<<index<<" < "<<nCh<<G4endl; |
---|
| 1369 | #endif |
---|
| 1370 | if(index == nCh) |
---|
| 1371 | { |
---|
| 1372 | G4cout<<"***G4QLowEnergy::PoStDI:Decay is impossible,totM="<<totM<<",GSM="<<tM<<G4endl; |
---|
| 1373 | throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound"); |
---|
| 1374 | } |
---|
| 1375 | // @@ Convert it to G4QHadrons |
---|
| 1376 | G4DynamicParticle* ResSec = new G4DynamicParticle; |
---|
| 1377 | G4DynamicParticle* FstSec = new G4DynamicParticle; |
---|
| 1378 | G4DynamicParticle* SecSec = new G4DynamicParticle; |
---|
| 1379 | #ifdef debug |
---|
| 1380 | G4cout<<"G4QLowEn::PSDI: Dynamic particles are created pL="<<pL<<G4endl; |
---|
| 1381 | #endif |
---|
| 1382 | |
---|
| 1383 | G4LorentzVector res4Mom(zeroMom,mass[index]*MeV); // The recoil nucleus prototype |
---|
| 1384 | G4double mF=0.; |
---|
| 1385 | G4double mS=0.; |
---|
| 1386 | G4int rA=totZ+totN; |
---|
| 1387 | G4int rZ=totZ; |
---|
| 1388 | G4int rL=pL; |
---|
| 1389 | G4int complete=3; |
---|
| 1390 | G4ParticleDefinition* theDefinition; // Prototype for qfNucleon |
---|
| 1391 | switch( index ) |
---|
| 1392 | { |
---|
| 1393 | case 0: |
---|
| 1394 | if(!evaporate || rA<2) |
---|
| 1395 | { |
---|
| 1396 | if(!rZ) theDefinition=aNeutron; |
---|
| 1397 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1398 | if(!theDefinition) |
---|
| 1399 | G4cerr<<"-Warning-G4LE::PSDI: notDef(1), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1400 | ResSec->SetDefinition( theDefinition ); |
---|
| 1401 | FstSec->SetDefinition( aGamma ); |
---|
| 1402 | SecSec->SetDefinition( aGamma ); |
---|
| 1403 | } |
---|
| 1404 | else |
---|
| 1405 | { |
---|
| 1406 | delete ResSec; |
---|
| 1407 | delete FstSec; |
---|
| 1408 | delete SecSec; |
---|
| 1409 | complete=0; |
---|
| 1410 | } |
---|
| 1411 | break; |
---|
| 1412 | case 1: |
---|
| 1413 | rA-=1; // gamma+n |
---|
| 1414 | if(!evaporate || rA<2) |
---|
| 1415 | { |
---|
| 1416 | if(!rZ) theDefinition=aNeutron; |
---|
| 1417 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1418 | if(!theDefinition) |
---|
| 1419 | G4cerr<<"-Warning-G4LE::PSDI: notDef(2), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1420 | ResSec->SetDefinition( theDefinition ); |
---|
| 1421 | SecSec->SetDefinition( aGamma ); |
---|
| 1422 | } |
---|
| 1423 | else |
---|
| 1424 | { |
---|
| 1425 | delete ResSec; |
---|
| 1426 | delete SecSec; |
---|
| 1427 | complete=1; |
---|
| 1428 | } |
---|
| 1429 | FstSec->SetDefinition( aNeutron ); |
---|
| 1430 | mF=mNeut; // First hadron 4-momentum |
---|
| 1431 | break; |
---|
| 1432 | case 2: |
---|
| 1433 | rA-=1; |
---|
| 1434 | rZ-=1; // gamma+p |
---|
| 1435 | if(!evaporate || rA<2) |
---|
| 1436 | { |
---|
| 1437 | if(!rZ) theDefinition=aNeutron; |
---|
| 1438 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1439 | if(!theDefinition) |
---|
| 1440 | G4cerr<<"-Warning-G4LE::PSDI: notDef(3), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1441 | ResSec->SetDefinition( theDefinition ); |
---|
| 1442 | SecSec->SetDefinition( aGamma ); |
---|
| 1443 | } |
---|
| 1444 | else |
---|
| 1445 | { |
---|
| 1446 | delete ResSec; |
---|
| 1447 | delete SecSec; |
---|
| 1448 | complete=1; |
---|
| 1449 | } |
---|
| 1450 | FstSec->SetDefinition( aProton ); |
---|
| 1451 | mF=mProt; // First hadron 4-momentum |
---|
| 1452 | break; |
---|
| 1453 | case 3: |
---|
| 1454 | rA-=2; |
---|
| 1455 | rZ-=1; // gamma+d |
---|
| 1456 | if(!evaporate || rA<2) |
---|
| 1457 | { |
---|
| 1458 | if(!rZ) theDefinition=aNeutron; |
---|
| 1459 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1460 | if(!theDefinition) |
---|
| 1461 | G4cerr<<"-Warning-G4LE::PSDI: notDef(4), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1462 | ResSec->SetDefinition( theDefinition ); |
---|
| 1463 | SecSec->SetDefinition( aGamma ); |
---|
| 1464 | } |
---|
| 1465 | else |
---|
| 1466 | { |
---|
| 1467 | delete ResSec; |
---|
| 1468 | delete SecSec; |
---|
| 1469 | complete=1; |
---|
| 1470 | } |
---|
| 1471 | FstSec->SetDefinition( aDeuteron ); |
---|
| 1472 | mF=mDeut; // First hadron 4-momentum |
---|
| 1473 | break; |
---|
| 1474 | case 4: |
---|
| 1475 | rA-=3; // gamma+t |
---|
| 1476 | rZ-=1; |
---|
| 1477 | if(!evaporate || rA<2) |
---|
| 1478 | { |
---|
| 1479 | if(!rZ) theDefinition=aNeutron; |
---|
| 1480 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1481 | if(!theDefinition) |
---|
| 1482 | G4cerr<<"-Warning-G4LE::PSDI: notDef(5), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1483 | ResSec->SetDefinition( theDefinition ); |
---|
| 1484 | SecSec->SetDefinition( aGamma ); |
---|
| 1485 | } |
---|
| 1486 | else |
---|
| 1487 | { |
---|
| 1488 | delete ResSec; |
---|
| 1489 | delete SecSec; |
---|
| 1490 | complete=1; |
---|
| 1491 | } |
---|
| 1492 | FstSec->SetDefinition( aTriton ); |
---|
| 1493 | mF=mTrit; // First hadron 4-momentum |
---|
| 1494 | break; |
---|
| 1495 | case 5: // gamma+He3 |
---|
| 1496 | rA-=3; |
---|
| 1497 | rZ-=2; |
---|
| 1498 | if(!evaporate || rA<2) |
---|
| 1499 | { |
---|
| 1500 | if(!rZ) theDefinition=aNeutron; |
---|
| 1501 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1502 | if(!theDefinition) |
---|
| 1503 | G4cerr<<"-Warning-G4LE::PSDI: notDef(6), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1504 | ResSec->SetDefinition( theDefinition ); |
---|
| 1505 | SecSec->SetDefinition( aGamma ); |
---|
| 1506 | } |
---|
| 1507 | else |
---|
| 1508 | { |
---|
| 1509 | delete ResSec; |
---|
| 1510 | delete SecSec; |
---|
| 1511 | complete=1; |
---|
| 1512 | } |
---|
| 1513 | FstSec->SetDefinition( aHe3); |
---|
| 1514 | mF=mHel3; // First hadron 4-momentum |
---|
| 1515 | break; |
---|
| 1516 | case 6: |
---|
| 1517 | rA-=4; |
---|
| 1518 | rZ-=2; // gamma+He4 |
---|
| 1519 | if(!evaporate || rA<2) |
---|
| 1520 | { |
---|
| 1521 | if(!rZ) theDefinition=aNeutron; |
---|
| 1522 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1523 | if(!theDefinition) |
---|
| 1524 | G4cerr<<"-Warning-G4LE::PSDI: notDef(7), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1525 | ResSec->SetDefinition( theDefinition ); |
---|
| 1526 | SecSec->SetDefinition( aGamma ); |
---|
| 1527 | } |
---|
| 1528 | else |
---|
| 1529 | { |
---|
| 1530 | delete ResSec; |
---|
| 1531 | delete SecSec; |
---|
| 1532 | complete=1; |
---|
| 1533 | } |
---|
| 1534 | FstSec->SetDefinition( anAlpha ); |
---|
| 1535 | mF=mAlph; // First hadron 4-momentum |
---|
| 1536 | break; |
---|
| 1537 | case 7: |
---|
| 1538 | rA-=2; // n+n |
---|
| 1539 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1540 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1541 | if(!theDefinition) |
---|
| 1542 | G4cerr<<"-Warning-G4LE::PSDI: notDef(8), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1543 | ResSec->SetDefinition( theDefinition ); |
---|
| 1544 | FstSec->SetDefinition( aNeutron ); |
---|
| 1545 | SecSec->SetDefinition( aNeutron ); |
---|
| 1546 | mF=mNeut; // First hadron 4-momentum |
---|
| 1547 | mS=mNeut; // Second hadron 4-momentum |
---|
| 1548 | break; |
---|
| 1549 | case 8: |
---|
| 1550 | rZ-=1; |
---|
| 1551 | rA-=2; // n+p |
---|
| 1552 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1553 | else if(rA==2 && !rZ) |
---|
| 1554 | { |
---|
| 1555 | index=7; |
---|
| 1556 | ResSec->SetDefinition( aDeuteron); |
---|
| 1557 | FstSec->SetDefinition( aNeutron ); |
---|
| 1558 | SecSec->SetDefinition( aNeutron ); |
---|
| 1559 | mF=mNeut; // First hadron 4-momentum |
---|
| 1560 | mS=mNeut; // Second hadron 4-momentum |
---|
| 1561 | break; |
---|
| 1562 | } |
---|
| 1563 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1564 | if(!theDefinition) |
---|
| 1565 | G4cerr<<"-Warning-G4LE::PSDI: notDef(9), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1566 | ResSec->SetDefinition( theDefinition ); |
---|
| 1567 | FstSec->SetDefinition( aNeutron ); |
---|
| 1568 | SecSec->SetDefinition( aProton ); |
---|
| 1569 | mF=mNeut; // First hadron 4-momentum |
---|
| 1570 | mS=mProt; // Second hadron 4-momentum |
---|
| 1571 | break; |
---|
| 1572 | case 9: |
---|
| 1573 | rZ-=2; |
---|
| 1574 | rA-=2; // p+p |
---|
| 1575 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1576 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1577 | if(!theDefinition) |
---|
| 1578 | G4cerr<<"-Warning-G4LE::PSDI: notDef(10), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1579 | ResSec->SetDefinition( theDefinition ); |
---|
| 1580 | FstSec->SetDefinition( aProton ); |
---|
| 1581 | SecSec->SetDefinition( aProton ); |
---|
| 1582 | mF=mProt; // First hadron 4-momentum |
---|
| 1583 | mS=mProt; // Second hadron 4-momentum |
---|
| 1584 | break; |
---|
| 1585 | case 10: |
---|
| 1586 | rZ-=2; |
---|
| 1587 | rA-=3; // p+d |
---|
| 1588 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1589 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1590 | if(!theDefinition) |
---|
| 1591 | G4cerr<<"-Warning-G4LE::PSDI: notDef(11), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1592 | ResSec->SetDefinition( theDefinition ); |
---|
| 1593 | FstSec->SetDefinition( aProton ); |
---|
| 1594 | SecSec->SetDefinition( aDeuteron ); |
---|
| 1595 | mF=mProt; // First hadron 4-momentum |
---|
| 1596 | mS=mDeut; // Second hadron 4-momentum |
---|
| 1597 | break; |
---|
| 1598 | case 11: |
---|
| 1599 | rZ-=1; |
---|
| 1600 | rA-=3; // n+d |
---|
| 1601 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1602 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1603 | if(!theDefinition) |
---|
| 1604 | G4cerr<<"-Warning-G4LE::PSDI: notDef(12), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1605 | ResSec->SetDefinition( theDefinition ); |
---|
| 1606 | FstSec->SetDefinition( aNeutron ); |
---|
| 1607 | SecSec->SetDefinition( aDeuteron ); |
---|
| 1608 | mF=mNeut; // First hadron 4-momentum |
---|
| 1609 | mS=mDeut; // Second hadron 4-momentum |
---|
| 1610 | break; |
---|
| 1611 | case 12: |
---|
| 1612 | rZ-=2; |
---|
| 1613 | rA-=4; // d+d |
---|
| 1614 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1615 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1616 | if(!theDefinition) |
---|
| 1617 | G4cerr<<"-Warning-G4LE::PSDI: notDef(13), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1618 | ResSec->SetDefinition( theDefinition ); |
---|
| 1619 | FstSec->SetDefinition( aDeuteron ); |
---|
| 1620 | SecSec->SetDefinition( aDeuteron ); |
---|
| 1621 | mF=mDeut; // First hadron 4-momentum |
---|
| 1622 | mS=mDeut; // Second hadron 4-momentum |
---|
| 1623 | break; |
---|
| 1624 | case 13: |
---|
| 1625 | rZ-=2; |
---|
| 1626 | rA-=4; // p+t |
---|
| 1627 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1628 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1629 | if(!theDefinition) |
---|
| 1630 | G4cerr<<"-Warning-G4LE::PSDI: notDef(14), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1631 | ResSec->SetDefinition( theDefinition ); |
---|
| 1632 | FstSec->SetDefinition( aProton ); |
---|
| 1633 | SecSec->SetDefinition( aTriton ); |
---|
| 1634 | mF=mProt; // First hadron 4-momentum |
---|
| 1635 | mS=mTrit; // Second hadron 4-momentum |
---|
| 1636 | break; |
---|
| 1637 | case 14: |
---|
| 1638 | rZ-=1; |
---|
| 1639 | rA-=4; // n+t |
---|
| 1640 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1641 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1642 | if(!theDefinition) |
---|
| 1643 | G4cerr<<"-Warning-G4LE::PSDI: notDef(15), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1644 | ResSec->SetDefinition( theDefinition ); |
---|
| 1645 | FstSec->SetDefinition( aNeutron ); |
---|
| 1646 | SecSec->SetDefinition( aTriton ); |
---|
| 1647 | mF=mNeut; // First hadron 4-momentum |
---|
| 1648 | mS=mTrit; // Second hadron 4-momentum |
---|
| 1649 | break; |
---|
| 1650 | case 15: |
---|
| 1651 | rZ-=3; |
---|
| 1652 | rA-=4; // p+He3 |
---|
| 1653 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1654 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1655 | if(!theDefinition) |
---|
| 1656 | G4cerr<<"-Warning-G4LE::PSDI: notDef(16), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1657 | ResSec->SetDefinition( theDefinition ); |
---|
| 1658 | FstSec->SetDefinition( aProton); |
---|
| 1659 | SecSec->SetDefinition( aHe3 ); |
---|
| 1660 | mF=mProt; // First hadron 4-momentum |
---|
| 1661 | mS=mHel3; // Second hadron 4-momentum |
---|
| 1662 | break; |
---|
| 1663 | case 16: |
---|
| 1664 | rZ-=2; |
---|
| 1665 | rA-=4; // n+He3 |
---|
| 1666 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1667 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1668 | if(!theDefinition) |
---|
| 1669 | G4cerr<<"-Warning-G4LE::PSDI: notDef(17), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1670 | ResSec->SetDefinition( theDefinition ); |
---|
| 1671 | FstSec->SetDefinition( aNeutron ); |
---|
| 1672 | SecSec->SetDefinition( aHe3 ); |
---|
| 1673 | mF=mNeut; // First hadron 4-momentum |
---|
| 1674 | mS=mHel3; // Second hadron 4-momentum |
---|
| 1675 | break; |
---|
| 1676 | case 17: |
---|
| 1677 | rZ-=3; |
---|
| 1678 | rA-=5; // p+alph |
---|
| 1679 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1680 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1681 | if(!theDefinition) |
---|
| 1682 | G4cerr<<"-Warning-G4LE::PSDI: notDef(18), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1683 | ResSec->SetDefinition( theDefinition ); |
---|
| 1684 | FstSec->SetDefinition( aProton ); |
---|
| 1685 | SecSec->SetDefinition( anAlpha ); |
---|
| 1686 | mF=mProt; // First hadron 4-momentum |
---|
| 1687 | mS=mAlph; // Second hadron 4-momentum |
---|
| 1688 | break; |
---|
| 1689 | case 18: |
---|
| 1690 | rZ-=2; |
---|
| 1691 | rA-=5; // n+alph |
---|
| 1692 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1693 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1694 | if(!theDefinition) |
---|
| 1695 | G4cerr<<"-Warning-G4LE::PSDI: notDef(19), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1696 | ResSec->SetDefinition( theDefinition ); |
---|
| 1697 | FstSec->SetDefinition( aNeutron ); |
---|
| 1698 | SecSec->SetDefinition( anAlpha ); |
---|
| 1699 | mF=mNeut; // First hadron 4-momentum |
---|
| 1700 | mS=mAlph; // Second hadron 4-momentum |
---|
| 1701 | break; |
---|
| 1702 | case 19: |
---|
| 1703 | rL-=1; // L+gamma (@@ rA inludes rL?) |
---|
| 1704 | rA-=1; |
---|
| 1705 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1706 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1707 | if(!theDefinition) |
---|
| 1708 | G4cerr<<"-Warning-G4LE::PSDI: notDef(20), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1709 | ResSec->SetDefinition( theDefinition ); |
---|
| 1710 | FstSec->SetDefinition( aLambda ); |
---|
| 1711 | SecSec->SetDefinition( aGamma ); |
---|
| 1712 | mF=mLamb; // First hadron 4-momentum |
---|
| 1713 | break; |
---|
| 1714 | case 20: |
---|
| 1715 | rL-=1; // L+p (@@ rA inludes rL?) |
---|
| 1716 | rZ-=1; |
---|
| 1717 | rA-=2; |
---|
| 1718 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1719 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1720 | if(!theDefinition) |
---|
| 1721 | G4cerr<<"-Warning-G4LE::PSDI: notDef(21), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1722 | ResSec->SetDefinition( theDefinition ); |
---|
| 1723 | FstSec->SetDefinition( aProton ); |
---|
| 1724 | SecSec->SetDefinition( aLambda ); |
---|
| 1725 | mF=mProt; // First hadron 4-momentum |
---|
| 1726 | mS=mLamb; // Second hadron 4-momentum |
---|
| 1727 | break; |
---|
| 1728 | case 21: |
---|
| 1729 | rL-=1; // L+n (@@ rA inludes rL?) |
---|
| 1730 | rA-=2; |
---|
| 1731 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1732 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1733 | if(!theDefinition) |
---|
| 1734 | G4cerr<<"-Warning-G4LE::PSDI: notDef(22), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1735 | ResSec->SetDefinition( theDefinition ); |
---|
| 1736 | FstSec->SetDefinition( aNeutron ); |
---|
| 1737 | SecSec->SetDefinition( aLambda ); |
---|
| 1738 | mF=mNeut; // First hadron 4-momentum |
---|
| 1739 | mS=mLamb; // Second hadron 4-momentum |
---|
| 1740 | break; |
---|
| 1741 | case 22: |
---|
| 1742 | rL-=1; // L+d (@@ rA inludes rL?) |
---|
| 1743 | rZ-=1; |
---|
| 1744 | rA-=3; |
---|
| 1745 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1746 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1747 | if(!theDefinition) |
---|
| 1748 | G4cerr<<"-Warning-G4LE::PSDI: notDef(23), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1749 | ResSec->SetDefinition( theDefinition ); |
---|
| 1750 | FstSec->SetDefinition( aDeuteron ); |
---|
| 1751 | SecSec->SetDefinition( aLambda ); |
---|
| 1752 | mF=mDeut; // First hadron 4-momentum |
---|
| 1753 | mS=mLamb; // Second hadron 4-momentum |
---|
| 1754 | break; |
---|
| 1755 | case 23: |
---|
| 1756 | rL-=1; // L+t (@@ rA inludes rL?) |
---|
| 1757 | rZ-=1; |
---|
| 1758 | rA-=4; |
---|
| 1759 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1760 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1761 | if(!theDefinition) |
---|
| 1762 | G4cerr<<"-Warning-G4LE::PSDI: notDef(24), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1763 | ResSec->SetDefinition( theDefinition ); |
---|
| 1764 | FstSec->SetDefinition( aTriton ); |
---|
| 1765 | SecSec->SetDefinition( aLambda ); |
---|
| 1766 | mF=mTrit; // First hadron 4-momentum |
---|
| 1767 | mS=mLamb; // Second hadron 4-momentum |
---|
| 1768 | break; |
---|
| 1769 | case 24: |
---|
| 1770 | rL-=1; // L+He3 (@@ rA inludes rL?) |
---|
| 1771 | rZ-=2; |
---|
| 1772 | rA-=4; |
---|
| 1773 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1774 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1775 | if(!theDefinition) |
---|
| 1776 | G4cerr<<"-Warning-G4LE::PSDI: notDef(25), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1777 | ResSec->SetDefinition( theDefinition ); |
---|
| 1778 | FstSec->SetDefinition( aHe3 ); |
---|
| 1779 | SecSec->SetDefinition( aLambda ); |
---|
| 1780 | mF=mHel3; // First hadron 4-momentum |
---|
| 1781 | mS=mLamb; // Second hadron 4-momentum |
---|
| 1782 | break; |
---|
| 1783 | case 25: |
---|
| 1784 | rL-=1; // L+alph (@@ rA inludes rL?) |
---|
| 1785 | rZ-=2; |
---|
| 1786 | rA-=5; |
---|
| 1787 | if(rA==1 && !rZ) theDefinition=aNeutron; |
---|
| 1788 | else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ); |
---|
| 1789 | if(!theDefinition) |
---|
| 1790 | G4cerr<<"-Warning-G4LE::PSDI: notDef(26), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl; |
---|
| 1791 | ResSec->SetDefinition( theDefinition ); |
---|
| 1792 | FstSec->SetDefinition( anAlpha ); |
---|
| 1793 | SecSec->SetDefinition( aLambda ); |
---|
| 1794 | mF=mAlph; // First hadron 4-momentum |
---|
| 1795 | mS=mLamb; // Second hadron 4-momentum |
---|
| 1796 | break; |
---|
| 1797 | } |
---|
| 1798 | #ifdef debug |
---|
| 1799 | G4cout<<"G4QLowEn::PSDI:F="<<mF<<",S="<<mS<<",com="<<complete<<",ev="<<evaporate<<G4endl; |
---|
| 1800 | #endif |
---|
| 1801 | G4LorentzVector fst4Mom(zeroMom,mF); // Prototype of the first hadron 4-momentum |
---|
| 1802 | G4LorentzVector snd4Mom(zeroMom,mS); // Prototype of the second hadron 4-momentum |
---|
| 1803 | G4LorentzVector dir4Mom=tot4M; // Prototype of the resN decay direction 4-momentum |
---|
| 1804 | dir4Mom.setE(tot4M.e()/2.); // Get half energy and total 3-momentum |
---|
| 1805 | // @@ Can be repeated to take into account the Coulomb Barrier |
---|
| 1806 | if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp)) |
---|
| 1807 | { // |
---|
| 1808 | G4cerr<<"**G4LowEnergy::PoStDoIt:i="<<index<<",tM="<<totM<<"->M1="<<res4Mom.m()<<"+M2=" |
---|
| 1809 | <<fst4Mom.m()<<"+M3="<<snd4Mom.m()<<"=="<<res4Mom.m()+fst4Mom.m()+snd4Mom.m()<<G4endl; |
---|
| 1810 | throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound"); |
---|
| 1811 | } // |
---|
| 1812 | #ifdef debug |
---|
| 1813 | G4cout<<"G4QLowEn::PSDI:r4M="<<res4Mom<<",f4M="<<fst4Mom<<",s4M="<<snd4Mom<<G4endl; |
---|
| 1814 | #endif |
---|
| 1815 | G4Track* aNewTrack = 0; |
---|
| 1816 | if(complete) |
---|
| 1817 | { |
---|
| 1818 | FstSec->Set4Momentum(fst4Mom); |
---|
| 1819 | aNewTrack = new G4Track(FstSec, localtime, position ); |
---|
| 1820 | aNewTrack->SetWeight(weight); // weighted |
---|
| 1821 | aNewTrack->SetTouchableHandle(trTouchable); |
---|
| 1822 | result.push_back( aNewTrack ); |
---|
| 1823 | EnMomConservation-=fst4Mom; |
---|
| 1824 | #ifdef debug |
---|
| 1825 | G4cout<<"G4QLowEn::PSDI: ***Filled*** 1stH4M="<<fst4Mom |
---|
| 1826 | <<", PDG="<<FstSec->GetDefinition()->GetPDGEncoding()<<G4endl; |
---|
| 1827 | #endif |
---|
| 1828 | if(complete>2) // Final solution |
---|
| 1829 | { |
---|
| 1830 | ResSec->Set4Momentum(res4Mom); |
---|
| 1831 | aNewTrack = new G4Track(ResSec, localtime, position ); |
---|
| 1832 | aNewTrack->SetWeight(weight); // weighted |
---|
| 1833 | aNewTrack->SetTouchableHandle(trTouchable); |
---|
| 1834 | result.push_back( aNewTrack ); |
---|
| 1835 | EnMomConservation-=res4Mom; |
---|
| 1836 | #ifdef debug |
---|
| 1837 | G4cout<<"G4QLowEn::PSDI: ***Filled*** rA4M="<<res4Mom<<",rZ="<<rZ<<",rA="<<rA<<",rL=" |
---|
| 1838 | <<rL<<G4endl; |
---|
| 1839 | #endif |
---|
| 1840 | SecSec->Set4Momentum(snd4Mom); |
---|
| 1841 | aNewTrack = new G4Track(SecSec, localtime, position ); |
---|
| 1842 | aNewTrack->SetWeight(weight); // weighted |
---|
| 1843 | aNewTrack->SetTouchableHandle(trTouchable); |
---|
| 1844 | result.push_back( aNewTrack ); |
---|
| 1845 | EnMomConservation-=snd4Mom; |
---|
| 1846 | #ifdef debug |
---|
| 1847 | G4cout<<"G4QLowEn::PSDI: ***Filled*** 2ndH4M="<<snd4Mom |
---|
| 1848 | <<", PDG="<<SecSec->GetDefinition()->GetPDGEncoding()<<G4endl; |
---|
| 1849 | #endif |
---|
| 1850 | } |
---|
| 1851 | else res4Mom+=snd4Mom; |
---|
| 1852 | } |
---|
| 1853 | else res4Mom=tot4M; |
---|
| 1854 | if(complete<3) // Evaporation of the residual must be done |
---|
| 1855 | { |
---|
| 1856 | G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom); // Input hadron-nucleus |
---|
| 1857 | G4QHadronVector* evaHV = new G4QHadronVector; // Output vector of hadrons (delete!) |
---|
[1337] | 1858 | Nuc.EvaporateNucleus(rHadron, evaHV); // here a pion can appear ! |
---|
[1199] | 1859 | G4int nOut=evaHV->size(); |
---|
| 1860 | for(G4int i=0; i<nOut; i++) |
---|
| 1861 | { |
---|
| 1862 | G4QHadron* curH = (*evaHV)[i]; |
---|
| 1863 | G4int hPDG=curH->GetPDGCode(); |
---|
| 1864 | G4LorentzVector h4Mom=curH->Get4Momentum(); |
---|
| 1865 | EnMomConservation-=h4Mom; |
---|
| 1866 | #ifdef debug |
---|
| 1867 | G4cout<<"G4QLowEn::PSDI: ***FillingCand#"<<i<<"*** evaH="<<hPDG<<h4Mom<<G4endl; |
---|
| 1868 | #endif |
---|
| 1869 | if (hPDG==90000001 || hPDG==2112) theDefinition = aNeutron; |
---|
| 1870 | else if(hPDG==90001000 || hPDG==2212) theDefinition = aProton; |
---|
| 1871 | else if(hPDG==91000000 || hPDG==3122) theDefinition = aLambda; |
---|
| 1872 | else if(hPDG== 22 ) theDefinition = aGamma; |
---|
[1337] | 1873 | else if(hPDG== 111) theDefinition = aPiZero; |
---|
| 1874 | else if(hPDG== 211) theDefinition = aPiPlus; |
---|
| 1875 | else if(hPDG== -211) theDefinition = aPiMinus; |
---|
[1199] | 1876 | else |
---|
| 1877 | { |
---|
| 1878 | G4int hZ=curH->GetCharge(); |
---|
| 1879 | G4int hA=curH->GetBaryonNumber(); |
---|
| 1880 | G4int hS=curH->GetStrangeness(); |
---|
| 1881 | theDefinition = G4ParticleTable::GetParticleTable()->FindIon(hZ,hA,hS,0); // ion |
---|
| 1882 | } |
---|
| 1883 | if(theDefinition) |
---|
| 1884 | { |
---|
| 1885 | G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom); |
---|
| 1886 | G4Track* evaQH = new G4Track(theEQH, localtime, position ); |
---|
| 1887 | evaQH->SetWeight(weight); // weighted |
---|
| 1888 | evaQH->SetTouchableHandle(trTouchable); |
---|
| 1889 | result.push_back( evaQH ); |
---|
| 1890 | } |
---|
| 1891 | else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl; |
---|
| 1892 | } |
---|
| 1893 | } |
---|
| 1894 | // algorithm implementation --- STOPS HERE |
---|
| 1895 | G4int nres=result.size(); |
---|
| 1896 | aParticleChange.SetNumberOfSecondaries(nres); |
---|
| 1897 | for(G4int i=0; i<nres; ++i) aParticleChange.AddSecondary(result[i]); |
---|
| 1898 | #ifdef debug |
---|
| 1899 | G4cout<<"G4QLowEnergy::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl; |
---|
| 1900 | #endif |
---|
| 1901 | return G4VDiscreteProcess::PostStepDoIt(track, step); |
---|
| 1902 | } |
---|
| 1903 | |
---|
| 1904 | G4double G4QLowEnergy::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG) |
---|
| 1905 | { |
---|
| 1906 | static G4bool first=true; |
---|
| 1907 | static G4VQCrossSection* CSmanager; |
---|
| 1908 | if(first) // Connection with a singletone |
---|
| 1909 | { |
---|
| 1910 | CSmanager=G4QIonIonCrossSection::GetPointer(); |
---|
| 1911 | if(PDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer(); |
---|
| 1912 | first=false; |
---|
| 1913 | } |
---|
| 1914 | #ifdef debug |
---|
| 1915 | G4cout<<"G4QLowE::CXS: *DONE* p="<<p<<",Z="<<Z<<",N="<<N<<",PDG="<<PDG<<G4endl; |
---|
| 1916 | #endif |
---|
| 1917 | return CSmanager->GetCrossSection(true, p, Z, N, PDG); |
---|
| 1918 | } |
---|