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