[968] | 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 | // |
---|
[1055] | 26 | // $Id: G4QIonIonElastic.cc,v 1.4 2009/02/23 09:49:24 mkossov Exp $ |
---|
| 27 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
[968] | 28 | // |
---|
| 29 | // ---------------- G4QIonIonElastic class ----------------- |
---|
| 30 | // by Mikhail Kossov, December 2006. |
---|
| 31 | // G4QIonIonElastic class of the CHIPS Simulation Branch in GEANT4 |
---|
| 32 | // --------------------------------------------------------------- |
---|
| 33 | // **************************************************************************************** |
---|
| 34 | // ********** This CLASS is temporary moved from the photolepton_hadron directory ********* |
---|
| 35 | // **************************************************************************************** |
---|
[1055] | 36 | // Short description: a simple process for the Ion-Ion elastic scattering. |
---|
| 37 | // For heavy by heavy ions it can reach 50% of the total cross-section. |
---|
| 38 | // ----------------------------------------------------------------------- |
---|
[968] | 39 | |
---|
| 40 | //#define debug |
---|
| 41 | //#define pdebug |
---|
| 42 | //#define tdebug |
---|
| 43 | //#define nandebug |
---|
| 44 | //#define ppdebug |
---|
| 45 | |
---|
| 46 | #include "G4QIonIonElastic.hh" |
---|
| 47 | |
---|
| 48 | // Initialization of static vectors |
---|
| 49 | G4int G4QIonIonElastic::nPartCWorld=152; // The#of particles initialized in CHIPS World |
---|
| 50 | std::vector<G4int> G4QIonIonElastic::ElementZ; // Z of the element(i) in theLastCalc |
---|
| 51 | std::vector<G4double> G4QIonIonElastic::ElProbInMat; // SumProbabilityElements in Material |
---|
| 52 | std::vector<std::vector<G4int>*> G4QIonIonElastic::ElIsoN; // N of isotope(j) of Element(i) |
---|
| 53 | std::vector<std::vector<G4double>*>G4QIonIonElastic::IsoProbInEl;//SumProbabIsotopes in ElI |
---|
| 54 | |
---|
| 55 | // Constructor |
---|
| 56 | G4QIonIonElastic::G4QIonIonElastic(const G4String& processName): |
---|
| 57 | G4VDiscreteProcess(processName, fHadronic) |
---|
| 58 | { |
---|
| 59 | #ifdef debug |
---|
| 60 | G4cout<<"G4QIonIonElastic::Constructor is called processName="<<processName<<G4endl; |
---|
| 61 | #endif |
---|
| 62 | if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl; |
---|
| 63 | SetProcessSubType(fHadronElastic); |
---|
| 64 | //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max) |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | // Destructor |
---|
| 68 | G4QIonIonElastic::~G4QIonIonElastic() {} |
---|
| 69 | |
---|
| 70 | // output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)), |
---|
| 71 | // where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3) |
---|
| 72 | // ********** All CHIPS cross sections are calculated in the surface units ************ |
---|
| 73 | G4double G4QIonIonElastic::GetMeanFreePath(const G4Track& aTrack, G4double, |
---|
| 74 | G4ForceCondition* Fc) |
---|
| 75 | { |
---|
| 76 | *Fc = NotForced; |
---|
| 77 | const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle(); |
---|
| 78 | G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition(); |
---|
| 79 | if( !IsApplicable(*incidentParticleDefinition)) |
---|
| 80 | G4cout<<"-Warning-G4QIonIonElastic::GetMeanFreePath: notImplementedParticle"<<G4endl; |
---|
| 81 | // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material |
---|
| 82 | G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle |
---|
| 83 | #ifdef debug |
---|
| 84 | G4double KinEn = incidentParticle->GetKineticEnergy(); |
---|
| 85 | G4cout<<"G4QIonIonElastic::GetMeanFreePath: kinE="<<KinEn<<",Mom="<<Momentum<<G4endl; |
---|
| 86 | #endif |
---|
| 87 | const G4Material* material = aTrack.GetMaterial(); // Get the current material |
---|
| 88 | const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume(); |
---|
| 89 | const G4ElementVector* theElementVector = material->GetElementVector(); |
---|
| 90 | G4int nE=material->GetNumberOfElements(); |
---|
| 91 | #ifdef debug |
---|
| 92 | G4cout<<"G4QIonIonElastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl; |
---|
| 93 | #endif |
---|
| 94 | G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer(); |
---|
| 95 | G4int pPDG=0; |
---|
| 96 | // Probably enough: pPDG=incidentParticleDefinition->GetPDGEncoding(); |
---|
| 97 | if ( incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 100001002; |
---|
| 98 | else if ( incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 100002004; |
---|
| 99 | else if ( incidentParticleDefinition == G4Triton::Triton() ) pPDG = 100001003; |
---|
| 100 | else if ( incidentParticleDefinition == G4He3::He3() ) pPDG = 100002003; |
---|
| 101 | else if ( incidentParticleDefinition == G4GenericIon::GenericIon() ) |
---|
| 102 | { |
---|
| 103 | pPDG=incidentParticleDefinition->GetPDGEncoding(); |
---|
| 104 | #ifdef debug |
---|
| 105 | G4int B=incidentParticleDefinition->GetBaryonNumber(); |
---|
| 106 | G4int C=incidentParticleDefinition->GetPDGCharge(); |
---|
| 107 | prPDG=100000000+1000*C+B; |
---|
| 108 | G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl; |
---|
| 109 | #endif |
---|
| 110 | } |
---|
| 111 | else G4cout<<"-Warning-G4QIonIonElastic::GetMeanFreePath:Unknown projectile Ion"<<G4endl; |
---|
| 112 | Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A |
---|
| 113 | G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton |
---|
| 114 | G4double sigma=0.; // Sums over elements for the material |
---|
| 115 | G4int IPIE=IsoProbInEl.size(); // How many old elements? |
---|
| 116 | if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI) |
---|
| 117 | { |
---|
| 118 | std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector |
---|
| 119 | SPI->clear(); |
---|
| 120 | delete SPI; |
---|
| 121 | std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector |
---|
| 122 | IsN->clear(); |
---|
| 123 | delete IsN; |
---|
| 124 | } |
---|
| 125 | ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE) |
---|
| 126 | ElementZ.clear(); // Clear the body vector for Z of Elements |
---|
| 127 | IsoProbInEl.clear(); // Clear the body vector for SPI |
---|
| 128 | ElIsoN.clear(); // Clear the body vector for N of Isotopes |
---|
| 129 | for(G4int i=0; i<nE; ++i) |
---|
| 130 | { |
---|
| 131 | G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element |
---|
| 132 | G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element |
---|
| 133 | ElementZ.push_back(Z); // Remember Z of the Element |
---|
| 134 | G4int isoSize=0; // The default for the isoVectorLength is 0 |
---|
| 135 | G4int indEl=0; // Index of non-natural element or 0(default) |
---|
| 136 | G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect |
---|
| 137 | if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector |
---|
| 138 | #ifdef debug |
---|
| 139 | G4cout<<"G4QIonIonElastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; |
---|
| 140 | #endif |
---|
| 141 | if(isoSize) // The Element has non-trivial abundance set |
---|
| 142 | { |
---|
| 143 | indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order |
---|
| 144 | #ifdef debug |
---|
| 145 | G4cout<<"G4QIIEl::GetMFP:iE="<<indEl<<", def="<<Isotopes->IsDefined(Z,indEl)<<G4endl; |
---|
| 146 | #endif |
---|
| 147 | if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define |
---|
| 148 | { |
---|
| 149 | std::vector<std::pair<G4int,G4double>*>* newAbund = |
---|
| 150 | new std::vector<std::pair<G4int,G4double>*>; |
---|
| 151 | G4double* abuVector=pElement->GetRelativeAbundanceVector(); |
---|
| 152 | for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes |
---|
| 153 | { |
---|
| 154 | G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z ! |
---|
| 155 | if(pElement->GetIsotope(j)->GetZ()!=Z) G4cerr<<"G4QIonIonEl::GetMeanFreePath Z=" |
---|
| 156 | <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl; |
---|
| 157 | G4double abund=abuVector[j]; |
---|
[1055] | 158 | std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund); |
---|
[968] | 159 | #ifdef debug |
---|
| 160 | G4cout<<"G4QIonIonElastic::GetMeanFP:pair#="<<j<<",N="<<N<<",ab="<<abund<<G4endl; |
---|
| 161 | #endif |
---|
| 162 | newAbund->push_back(pr); |
---|
[1055] | 163 | } |
---|
[968] | 164 | #ifdef debug |
---|
| 165 | G4cout<<"G4QIonIonElastic::GetMeanFP: pairVectorLength="<<newAbund->size()<<G4endl; |
---|
| 166 | #endif |
---|
| 167 | indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd |
---|
| 168 | for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary |
---|
| 169 | delete newAbund; // Was "new" in the beginning of the name space |
---|
| 170 | } |
---|
| 171 | } |
---|
| 172 | std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer |
---|
| 173 | std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector |
---|
| 174 | IsoProbInEl.push_back(SPI); |
---|
| 175 | std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector |
---|
| 176 | ElIsoN.push_back(IsN); |
---|
| 177 | G4int nIs=cs->size(); // A#Of Isotopes in the Element |
---|
| 178 | #ifdef debug |
---|
| 179 | G4cout<<"G4QIonIonEl::GetMFP:=***=> #isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl; |
---|
| 180 | #endif |
---|
| 181 | G4double susi=0.; // sum of CS over isotopes |
---|
| 182 | if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El |
---|
| 183 | { |
---|
| 184 | std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice |
---|
| 185 | G4int N=curIs->first; // #of Neuterons in the isotope j of El i |
---|
| 186 | IsN->push_back(N); // Remember Min N for the Element |
---|
| 187 | #ifdef debug |
---|
| 188 | G4cout<<"G4QIIEl::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl; |
---|
| 189 | #endif |
---|
[1055] | 190 | G4bool ccsf=false; // Extract elastic Ion-Ion cross-section |
---|
[968] | 191 | #ifdef debug |
---|
| 192 | G4cout<<"G4QIonIonElastic::GMFP: GetCS #1 j="<<j<<G4endl; |
---|
| 193 | #endif |
---|
| 194 | G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope |
---|
| 195 | #ifdef debug |
---|
| 196 | G4cout<<"G4QIonIonElastic::GMFP: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="<<Momentum |
---|
| 197 | <<", XSec="<<CSI/millibarn<<G4endl; |
---|
| 198 | #endif |
---|
| 199 | curIs->second = CSI; |
---|
| 200 | susi+=CSI; // Make a sum per isotopes |
---|
| 201 | SPI->push_back(susi); // Remember summed cross-section |
---|
| 202 | } // End of temporary initialization of the cross sections in the G4QIsotope singeltone |
---|
| 203 | sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV) |
---|
| 204 | #ifdef debug |
---|
| 205 | G4cout<<"G4QIonIonEl::GMFP:<S>="<<Isotopes->GetMeanCrossSection(Z,indEl)<<",AddToSig=" |
---|
| 206 | <<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl; |
---|
| 207 | #endif |
---|
| 208 | ElProbInMat.push_back(sigma); |
---|
| 209 | } // End of LOOP over Elements |
---|
| 210 | // Check that cross section is not zero and return the mean free path |
---|
| 211 | #ifdef debug |
---|
| 212 | G4cout<<"G4QIonIonElastic::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl; |
---|
| 213 | #endif |
---|
| 214 | if(sigma > 0.) return 1./sigma; // Mean path [distance] |
---|
| 215 | return DBL_MAX; |
---|
| 216 | } |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | G4bool G4QIonIonElastic::IsApplicable(const G4ParticleDefinition& particle) |
---|
| 220 | { |
---|
| 221 | if (particle == *( G4Deuteron::Deuteron() )) return true; |
---|
| 222 | else if (particle == *( G4Alpha::Alpha() )) return true; |
---|
| 223 | else if (particle == *( G4Triton::Triton() )) return true; |
---|
| 224 | else if (particle == *( G4He3::He3() )) return true; |
---|
| 225 | else if (particle == *( G4GenericIon::GenericIon() )) return true; |
---|
| 226 | #ifdef debug |
---|
| 227 | G4cout<<"***>>G4QIonIonElastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl; |
---|
| 228 | #endif |
---|
| 229 | return false; |
---|
| 230 | } |
---|
| 231 | |
---|
| 232 | G4VParticleChange* G4QIonIonElastic::PostStepDoIt(const G4Track& track, const G4Step& step) |
---|
| 233 | { |
---|
| 234 | static const G4double fm2MeV2 = 3*38938./1.09; // (3/1.09)*(hc)^2 in fm^2*MeV^2 |
---|
| 235 | static G4bool CWinit = true; // CHIPS Warld needs to be initted |
---|
| 236 | if(CWinit) |
---|
[1055] | 237 | { |
---|
[968] | 238 | CWinit=false; |
---|
| 239 | G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max) |
---|
| 240 | } |
---|
| 241 | //------------------------------------------------------------------------------------- |
---|
| 242 | const G4DynamicParticle* projHadron = track.GetDynamicParticle(); |
---|
| 243 | const G4ParticleDefinition* particle=projHadron->GetDefinition(); |
---|
| 244 | #ifdef debug |
---|
| 245 | G4cout<<"G4QIonIonElastic::PostStepDoIt: Before the GetMeanFreePath is called In4M=" |
---|
| 246 | <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type=" |
---|
| 247 | <<particle->GetParticleType()<<", Subtp="<<particle->GetParticleSubType()<<G4endl; |
---|
| 248 | #endif |
---|
| 249 | //G4ForceCondition cond=NotForced; |
---|
| 250 | //GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters? |
---|
| 251 | #ifdef debug |
---|
| 252 | G4cout<<"G4QIonIonElastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl; |
---|
| 253 | #endif |
---|
| 254 | G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV! |
---|
| 255 | G4LorentzVector scat4M=proj4M; // @@ Must be filled (?) |
---|
| 256 | G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV |
---|
| 257 | G4double Momentum = proj4M.rho(); // @@ Just for the test purposes |
---|
| 258 | if(std::fabs(Momentum-momentum)>.000001) |
---|
| 259 | G4cerr<<"-Warn-G4QIonIonElastic::PostStepDoIt:P(IU)="<<Momentum<<"#"<<momentum<<G4endl; |
---|
| 260 | G4double pM2=proj4M.m2(); // in MeV^2 |
---|
| 261 | G4double pM=std::sqrt(pM2); // in MeV |
---|
| 262 | #ifdef pdebug |
---|
| 263 | G4cout<<"G4QIonIonElastic::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",pM="<<pM |
---|
| 264 | <<",scat4M="<<scat4M<<scat4M.m()<<G4endl; |
---|
| 265 | #endif |
---|
| 266 | if (!IsApplicable(*particle)) // Check applicability |
---|
| 267 | { |
---|
| 268 | G4cerr<<"G4QIonIonElastic::PostStepDoIt: Only NA elastic is implemented."<<G4endl; |
---|
| 269 | return 0; |
---|
| 270 | } |
---|
| 271 | const G4Material* material = track.GetMaterial(); // Get the current material |
---|
| 272 | G4int Z=0; |
---|
| 273 | const G4ElementVector* theElementVector = material->GetElementVector(); |
---|
| 274 | G4int nE=material->GetNumberOfElements(); |
---|
| 275 | #ifdef debug |
---|
| 276 | G4cout<<"G4QIonIonElastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl; |
---|
| 277 | #endif |
---|
| 278 | // Probably enough: projPDG=particle->GetPDGEncoding(); |
---|
| 279 | G4int projPDG=0; // CHIPS PDG Code for the captured hadron |
---|
| 280 | if (particle == G4Deuteron::Deuteron() ) projPDG= 100001002; |
---|
| 281 | else if (particle == G4Alpha::Alpha() ) projPDG= 100002004; |
---|
| 282 | else if (particle == G4Triton::Triton() ) projPDG= 100001003; |
---|
| 283 | else if (particle == G4He3::He3() ) projPDG= 100002003; |
---|
| 284 | else if (particle == G4GenericIon::GenericIon() ) |
---|
| 285 | { |
---|
| 286 | projPDG=particle->GetPDGEncoding(); |
---|
| 287 | #ifdef debug |
---|
| 288 | G4int B=particle->GetBaryonNumber(); |
---|
| 289 | G4int C=particle->GetPDGCharge(); |
---|
| 290 | prPDG=100000000+1000*C+B; |
---|
| 291 | G4cout<<"G4QIonIonElastic::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl; |
---|
| 292 | #endif |
---|
| 293 | } |
---|
| 294 | else G4cout<<"-Warning-G4QIonIonElastic::PostStepDoIt:Unknown projectile Ion"<<G4endl; |
---|
| 295 | #ifdef debug |
---|
| 296 | G4int prPDG=particle->GetPDGEncoding(); |
---|
[1055] | 297 | G4cout<<"G4QIonIonElastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl; |
---|
[968] | 298 | #endif |
---|
| 299 | if(!projPDG) |
---|
| 300 | { |
---|
| 301 | G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:Undefined interactingNucleus"<<G4endl; |
---|
| 302 | return 0; |
---|
| 303 | } |
---|
| 304 | G4double pA=particle->GetBaryonNumber(); // Projectile A |
---|
| 305 | G4double pZ=particle->GetPDGCharge(); // Projectile Z |
---|
| 306 | G4double pN=pA-pZ; // Projectile N |
---|
| 307 | G4int EPIM=ElProbInMat.size(); |
---|
| 308 | #ifdef debug |
---|
[1055] | 309 | G4cout<<"G4QIonIonElastic::PSDI:m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl; |
---|
[968] | 310 | #endif |
---|
| 311 | G4int i=0; |
---|
| 312 | if(EPIM>1) |
---|
| 313 | { |
---|
| 314 | G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand(); |
---|
| 315 | for(i=0; i<nE; ++i) |
---|
[1055] | 316 | { |
---|
[968] | 317 | #ifdef debug |
---|
[1055] | 318 | G4cout<<"G4QIonIonElastic::PSDI: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl; |
---|
[968] | 319 | #endif |
---|
| 320 | if (rnd<ElProbInMat[i]) break; |
---|
| 321 | } |
---|
| 322 | if(i>=nE) i=nE-1; // Top limit for the Element |
---|
| 323 | } |
---|
| 324 | G4Element* pElement=(*theElementVector)[i]; |
---|
| 325 | Z=static_cast<G4int>(pElement->GetZ()); |
---|
| 326 | #ifdef debug |
---|
[1055] | 327 | G4cout<<"G4QIonIonElastic::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl; |
---|
[968] | 328 | #endif |
---|
| 329 | if(Z<=0) |
---|
| 330 | { |
---|
| 331 | G4cerr<<"---Warning---G4QIonIonElastic::PostStepDoIt: Element with Z="<<Z<<G4endl; |
---|
| 332 | if(Z<0) return 0; |
---|
| 333 | } |
---|
| 334 | std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes |
---|
| 335 | std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i] |
---|
| 336 | G4int nofIsot=SPI->size(); // #of isotopes in the element i |
---|
| 337 | #ifdef debug |
---|
[1055] | 338 | G4cout<<"G4QIonIonElastic::PosStDoIt: nI="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl; |
---|
[968] | 339 | #endif |
---|
| 340 | G4int j=0; |
---|
| 341 | if(nofIsot>1) |
---|
| 342 | { |
---|
| 343 | G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element |
---|
| 344 | for(j=0; j<nofIsot; ++j) |
---|
| 345 | { |
---|
| 346 | #ifdef debug |
---|
[1055] | 347 | G4cout<<"G4QIonIonElastic::PostStDI: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl; |
---|
[968] | 348 | #endif |
---|
| 349 | if(rndI < (*SPI)[j]) break; |
---|
| 350 | } |
---|
| 351 | if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope |
---|
| 352 | } |
---|
| 353 | G4int N =(*IsN)[j]; ; // Randomized number of neutrons |
---|
| 354 | #ifdef debug |
---|
[1055] | 355 | G4cout<<"G4QIonIonElastic::PostStepDoIt:j="<<i<<",N(isotope)="<<N<<", MeV="<<MeV<<G4endl; |
---|
[968] | 356 | #endif |
---|
| 357 | if(N<0) |
---|
| 358 | { |
---|
| 359 | G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:IsotopeZ="<<Z<<" has 0>N="<<N<<G4endl; |
---|
| 360 | return 0; |
---|
| 361 | } |
---|
| 362 | nOfNeutrons=N; // Remember it for the energy-momentum check |
---|
| 363 | #ifdef debug |
---|
| 364 | G4cout<<"G4QIonIonElastic::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl; |
---|
| 365 | #endif |
---|
| 366 | if(N<0) |
---|
| 367 | { |
---|
| 368 | G4cerr<<"-Warning-G4QIonIonElastic::PostStepDoIt:Element with N="<<N<< G4endl; |
---|
| 369 | return 0; |
---|
| 370 | } |
---|
| 371 | aParticleChange.Initialize(track); |
---|
| 372 | #ifdef debug |
---|
| 373 | G4cout<<"G4QIonIonElastic::PostStepDoIt: track is initialized"<<G4endl; |
---|
| 374 | #endif |
---|
| 375 | G4double weight = track.GetWeight(); |
---|
| 376 | G4double localtime = track.GetGlobalTime(); |
---|
| 377 | G4ThreeVector position = track.GetPosition(); |
---|
| 378 | #ifdef debug |
---|
| 379 | G4cout<<"G4QIonIonElastic::PostStepDoIt: before Touchable extraction"<<G4endl; |
---|
| 380 | #endif |
---|
| 381 | G4TouchableHandle trTouchable = track.GetTouchableHandle(); |
---|
| 382 | #ifdef debug |
---|
| 383 | G4cout<<"G4QIonIonElastic::PostStepDoIt: Touchable is extracted"<<G4endl; |
---|
| 384 | #endif |
---|
| 385 | // |
---|
| 386 | G4double tA=Z+N; |
---|
| 387 | G4int targPDG=90000000+Z*1000+N; // CHIPS PDG Code of the target nucleus |
---|
| 388 | G4QPDGCode targQPDG(targPDG); // @@ one can use G4Ion & get rid of CHIPS World |
---|
| 389 | G4double tM=targQPDG.GetMass(); // CHIPS target mass in MeV |
---|
| 390 | G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?) |
---|
| 391 | G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector |
---|
| 392 | G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction |
---|
| 393 | #ifdef debug |
---|
| 394 | G4cout<<"G4QIonIonElastic::PostStDI: tM="<<tM<<", p4M="<<proj4M<<", t4M="<<tot4M<<G4endl; |
---|
| 395 | #endif |
---|
| 396 | EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation |
---|
| 397 | G4VQCrossSection* ELmanager=G4QElasticCrossSection::GetPointer(); |
---|
| 398 | G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer(); |
---|
| 399 | // @@ Probably this is not necessary any more |
---|
| 400 | #ifdef debug |
---|
| 401 | G4cout<<"G4QIIEl::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl; |
---|
| 402 | #endif |
---|
| 403 | // false means elastic cross-section |
---|
| 404 | G4double xSec=CSmanager->GetCrossSection(false, Momentum, Z, N, projPDG);// Rec.CrossSect |
---|
| 405 | #ifdef debug |
---|
| 406 | G4cout<<"G4QIIEl::PSDI: pPDG="<<projPDG<<",P="<<Momentum<<",CS="<<xSec/millibarn<<G4endl; |
---|
| 407 | #endif |
---|
| 408 | #ifdef nandebug |
---|
| 409 | if(xSec>0. || xSec<0. || xSec==0); |
---|
| 410 | else G4cout<<"-NaN-Warning-G4QIonIonElastic::PostStDoIt: xSec="<<xSec/millibarn<<G4endl; |
---|
| 411 | #endif |
---|
| 412 | // @@ check a possibility to separate p, n, or alpha (!) |
---|
| 413 | if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing |
---|
| 414 | { |
---|
| 415 | #ifdef pdebug |
---|
| 416 | G4cerr<<"-Warning-G4QIonIonElastic::PostStDoIt: *Zero cross-section* PDG="<<projPDG |
---|
| 417 | <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl; |
---|
| 418 | #endif |
---|
| 419 | //Do Nothing Action insead of the reaction |
---|
| 420 | aParticleChange.ProposeEnergy(kinEnergy); |
---|
| 421 | aParticleChange.ProposeLocalEnergyDeposit(0.); |
---|
| 422 | aParticleChange.ProposeMomentumDirection(dir) ; |
---|
| 423 | return G4VDiscreteProcess::PostStepDoIt(track,step); |
---|
| 424 | } |
---|
| 425 | G4double dtM=tM+tM; |
---|
| 426 | G4double PA=Momentum*pA; |
---|
| 427 | G4double PA2=PA*PA; |
---|
| 428 | G4double maxt=dtM*PA2/(std::sqrt(PA2+pM2)+tM/2+pM2/dtM); |
---|
| 429 | #ifdef pdebug |
---|
| 430 | G4cout<<"G4QIonIonElastic::PostStDoIt:pPDG="<<projPDG<<",tPDG="<<targPDG<<",P="<<Momentum |
---|
| 431 | <<",CS="<<xSec<<",maxt="<<maxt<<G4endl; |
---|
| 432 | #endif |
---|
| 433 | xSec=ELmanager->GetCrossSection(false, Momentum, 1, 0, 2212);// pp=nn |
---|
| 434 | G4double B1=ELmanager->GetSlope(1,0,2212); // slope for pp=nn |
---|
| 435 | xSec=ELmanager->GetCrossSection(false, Momentum, 1, 0, 2112);// np=pn |
---|
| 436 | G4double B2 =ELmanager->GetSlope(1,0,2112); // slope for np=pn |
---|
| 437 | G4double mB =((pZ*Z+pN*N)*B1+(pZ*N+pN*Z)*B2)/(pA+tA); |
---|
| 438 | G4double pR2=std::pow(pA+4.,.305)/fm2MeV2; |
---|
| 439 | G4double tR2=std::pow(tA+4.,.305)/fm2MeV2; |
---|
| 440 | G4double eB =mB+pR2+tR2; |
---|
| 441 | G4double mint=-std::log(1.-G4UniformRand()*(1.-std::exp(-eB*maxt)))/eB; |
---|
| 442 | #ifdef pdebug |
---|
| 443 | G4cout<<"G4QIonIonElastic::PostStDoIt:B1="<<B1<<",B2="<<B2<<",mB="<<mB |
---|
| 444 | <<",pR2="<<pR2<<",tR2="<<tR2<<",eB="<<eB<<",mint="<<mint<<G4endl; |
---|
| 445 | #endif |
---|
| 446 | #ifdef nandebug |
---|
| 447 | if(mint>-.0000001); |
---|
| 448 | else G4cout<<"-Warning-G4QIonIonElastic::PostStDoIt:-t="<<mint<<G4endl; |
---|
| 449 | #endif |
---|
| 450 | G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS |
---|
| 451 | // |
---|
| 452 | #ifdef ppdebug |
---|
| 453 | G4cout<<"G4QIonIonElastic::PoStDoI:t="<<mint<<",dpcm2="<<CSmanager->GetHMaxT()<<",Ek=" |
---|
| 454 | <<kinEnergy<<",tM="<<tM<<",pM="<<pM<<",cost="<<cost<<G4endl; |
---|
| 455 | #endif |
---|
| 456 | if(cost>1. || cost<-1. || !(cost>-1. || cost<1.)) |
---|
| 457 | { |
---|
| 458 | if(cost>1.000001 || cost<-1.000001 || !(cost>-1. || cost<1.)) |
---|
| 459 | { |
---|
| 460 | G4double tM2=tM*tM; // Squared target mass |
---|
| 461 | G4double pEn=pM+kinEnergy; // tot projectile Energy in MeV |
---|
| 462 | G4double sM=dtM*pEn+tM2+pM2; // Mondelstam s |
---|
| 463 | G4double twop2cm=(tM2+tM2)*(pEn*pEn-pM2)/sM;// Max_t/2 (2*p^2_cm) |
---|
| 464 | G4cout<<"-Warning-G4QIonIonElastic::PoStDI:cos="<<cost<<",t="<<mint<<",T="<<kinEnergy |
---|
| 465 | <<",tM="<<tM<<",tmax="<<2*kinEnergy*tM<<",p="<<projPDG<<",t="<<targPDG<<G4endl; |
---|
| 466 | G4cout<<"G4QIonIonElastic::PSDI:dpcm2="<<twop2cm<<"="<<CSmanager->GetHMaxT()<<G4endl; |
---|
| 467 | } |
---|
| 468 | if (cost>1.) cost=1.; |
---|
| 469 | else if(cost<-1.) cost=-1.; |
---|
| 470 | } |
---|
| 471 | G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tM); // 4mom of the recoil target |
---|
| 472 | G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-tM-pM)*.01); |
---|
| 473 | if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost)) |
---|
| 474 | { |
---|
| 475 | G4cerr<<"G4QIonIonE::PSDI:t4M="<<tot4M<<",pM="<<pM<<",tM="<<tM<<",cost="<<cost<<G4endl; |
---|
| 476 | } |
---|
| 477 | #ifdef debug |
---|
| 478 | G4cout<<"G4QIonIonElast::PSDI:s4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl; |
---|
| 479 | G4cout<<"G4QIonIonElastic::PSDI:scatE="<<scat4M.e()-pM<<",recoE="<<reco4M.e()-tM<<",d4M=" |
---|
| 480 | <<tot4M-scat4M-reco4M<<G4endl; |
---|
| 481 | #endif |
---|
| 482 | // Update G4VParticleChange for the scattered projectile |
---|
| 483 | G4double finE=scat4M.e()-pM; // Final kinetic energy of the scattered proton |
---|
| 484 | if(finE>0.0) aParticleChange.ProposeEnergy(finE); |
---|
| 485 | else |
---|
| 486 | { |
---|
| 487 | if(finE<-1.e-8 || !(finE>-1.||finE<1.)) // NAN or negative |
---|
| 488 | G4cerr<<"*Warning*G4QIonIonElastic::PostStDoIt: Zero or negative scattered E="<<finE |
---|
| 489 | <<", s4M="<<scat4M<<", r4M="<<reco4M<<", d4M="<<tot4M-scat4M-reco4M<<G4endl; |
---|
| 490 | aParticleChange.ProposeEnergy(0.) ; |
---|
| 491 | aParticleChange.ProposeTrackStatus(fStopButAlive); |
---|
| 492 | } |
---|
| 493 | G4ThreeVector findir=scat4M.vect()/scat4M.rho(); // Unit vector in new direction |
---|
| 494 | aParticleChange.ProposeMomentumDirection(findir); // new direction for the scattered part |
---|
| 495 | EnMomConservation-=scat4M; // It must be initialized by (pE+tM,pP) |
---|
| 496 | // This is how in general the secondary should be identified |
---|
[1055] | 497 | G4DynamicParticle* theSec = new G4DynamicParticle; // A secondary for the recoil hadron |
---|
[968] | 498 | G4int aA = Z+N; |
---|
| 499 | #ifdef pdebug |
---|
[1055] | 500 | G4cout<<"G4QIonIonElastic::PostStepDoIt: Ion Z="<<Z<<", A="<<aA<<G4endl; |
---|
[968] | 501 | #endif |
---|
| 502 | G4ParticleDefinition* theDefinition=G4ParticleTable::GetParticleTable() |
---|
| 503 | ->FindIon(Z,aA,0,Z); |
---|
| 504 | if(!theDefinition)G4cout<<"-Warning-G4QIonIonElastic::PoStDI:drop PDG="<<targPDG<<G4endl; |
---|
| 505 | #ifdef pdebug |
---|
| 506 | G4cout<<"G4QIonIonElastic::PoStDI:RecoilName="<<theDefinition->GetParticleName()<<G4endl; |
---|
| 507 | #endif |
---|
| 508 | theSec->SetDefinition(theDefinition); |
---|
| 509 | |
---|
| 510 | EnMomConservation-=reco4M; |
---|
| 511 | #ifdef tdebug |
---|
| 512 | G4cout<<"G4QIonIonElastic::PSD:"<<targPDG<<reco4M<<reco4M.m()<<EnMomConservation<<G4endl; |
---|
| 513 | #endif |
---|
| 514 | theSec->Set4Momentum(reco4M); |
---|
| 515 | #ifdef debug |
---|
| 516 | G4ThreeVector curD=theSec->GetMomentumDirection(); |
---|
| 517 | G4double curM=theSec->GetMass(); |
---|
| 518 | G4double curE=theSec->GetKineticEnergy()+curM; |
---|
| 519 | G4cout<<"G4QIonIonElastic::PSDI: p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl; |
---|
| 520 | #endif |
---|
| 521 | // Make a recoil nucleus |
---|
| 522 | G4Track* aNewTrack = new G4Track(theSec, localtime, position ); |
---|
| 523 | aNewTrack->SetWeight(weight); // weighted |
---|
| 524 | aNewTrack->SetTouchableHandle(trTouchable); |
---|
| 525 | aParticleChange.AddSecondary( aNewTrack ); |
---|
| 526 | #ifdef debug |
---|
| 527 | G4cout<<"G4QIonIonElastic::PostStepDoIt: **** PostStepDoIt is done ****"<<G4endl; |
---|
| 528 | #endif |
---|
| 529 | return G4VDiscreteProcess::PostStepDoIt(track, step); |
---|
| 530 | } |
---|