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: G4QLowEnergy.cc,v 1.5 2010/06/14 16:11:27 mkossov Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-04-beta-01 $ |
---|
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(); |
---|
256 | static G4ParticleDefinition* aPiZero = G4PionZero::PionZero(); |
---|
257 | static G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus(); |
---|
258 | static G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus(); |
---|
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 | { |
---|
553 | G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer(); |
---|
554 | G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer(); |
---|
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 | } |
---|
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); |
---|
596 | if( xSec <= 0. ) break; // Break the while loop |
---|
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 |
---|
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 |
---|
742 | com4M += qfN4M; // Calculate Total 4Mom for NA scattering |
---|
743 | G4ThreeVector boostV=proj4M.vect()/proj4M.e(); |
---|
744 | qfN4M.boost(-boostV); |
---|
745 | G4double tNE = qfN4M.e(); // Energy of the QF nucleon in LS |
---|
746 | if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2); // Mom(s) fake value |
---|
747 | else break; // Break the while loop |
---|
748 | } |
---|
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 |
---|
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; |
---|
767 | break; // Break the while loop |
---|
768 | } |
---|
769 | G4Track* targSpect = 0; |
---|
770 | G4Track* aNucTrack = 0; |
---|
771 | if(tB > tD) // Fill the residual nucleus |
---|
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 | } |
---|
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 | } |
---|
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!) |
---|
1858 | Nuc.EvaporateNucleus(rHadron, evaHV); // here a pion can appear ! |
---|
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; |
---|
1873 | else if(hPDG== 111) theDefinition = aPiZero; |
---|
1874 | else if(hPDG== 211) theDefinition = aPiPlus; |
---|
1875 | else if(hPDG== -211) theDefinition = aPiMinus; |
---|
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 | } |
---|