Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QDiffraction.cc @ 968

Last change on this file since 968 was 962, checked in by garnier, 15 years ago
update processes
File size: 28.5 KB

Line
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: G4QDiffraction.cc,v 1.5 2008/10/02 21:10:07 dennis Exp $
27	// GEANT4 tag $Name: geant4-09-02-ref-02 $
28	//
29	// ---------------- G4QDiffraction class -----------------
30	// by Mikhail Kossov, Aug 2007.
31	// G4QDiffraction class of the CHIPS Simulation Branch in GEANT4
32	// ---------------------------------------------------------------
33	// ****************************************************************************************
34	// ******** This CLASS is temporary moved from the "chips/interface" directory *******
35	// ****************************************************************************************
36
37	//#define debug
38	//#define pdebug
39	//#define tdebug
40	//#define nandebug
41	//#define ppdebug
42
43	#include "G4QDiffraction.hh"
44
45	// Initialization of static vectors
46	G4int G4QDiffraction::nPartCWorld=152; // #of particles initialized in CHIPS
47	std::vector<G4int> G4QDiffraction::ElementZ; // Z of element(i) in theLastCalc
48	std::vector<G4double> G4QDiffraction::ElProbInMat; // SumProbOfElem in Material
49	std::vector<std::vector<G4int>*> G4QDiffraction::ElIsoN;// N of isotope(j), E(i)
50	std::vector<std::vector<G4double>*>G4QDiffraction::IsoProbInEl;//SumProbIsotE(i)
51
52	// Constructor
53	G4QDiffraction::G4QDiffraction(const G4String& processName):
54	G4VDiscreteProcess(processName, fHadronic)
55	{
56	#ifdef debug
57	G4cout<<"G4QDiffraction::Constructor is called processName="<<processName<<G4endl;
58	#endif
59	if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
60	SetProcessSubType(fHadronInelastic);
61	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
62	}
63
64	// Destructor
65	G4QDiffraction::~G4QDiffraction() {}
66
67
68	G4LorentzVector G4QDiffraction::GetEnegryMomentumConservation(){return EnMomConservation;}
69
70	G4int G4QDiffraction::GetNumberOfNeutronsInTarget() {return nOfNeutrons;}
71
72	// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
73	// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
74	// ******** All CHIPS cross sections are calculated in the surface units **********
75	G4double G4QDiffraction::GetMeanFreePath(const G4Track&Track,G4double Q,G4ForceCondition*F)
76	{
77	*F = NotForced;
78	const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
79	G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
80	if( !IsApplicable(*incidentParticleDefinition))
81	G4cout<<"-Warning-G4QDiffraction::GetMeanFreePath for notImplemented Particle"<<G4endl;
82	// Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
83	G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
84	#ifdef debug
85	G4double KinEn = incidentParticle->GetKineticEnergy();
86	G4cout<<"G4QDiffraction::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
87	#endif
88	const G4Material* material = Track.GetMaterial(); // Get the current material
89	const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
90	const G4ElementVector* theElementVector = material->GetElementVector();
91	G4int nE=material->GetNumberOfElements();
92	#ifdef debug
93	G4cout<<"G4QDiffraction::GetMeanFreePath:"<<nE<<" Elems in Material="<<*material<<G4endl;
94	#endif
95	G4int pPDG=0;
96	// @@ At present it is made only for n & p, but can be extended if inXS are available
97	if (incidentParticleDefinition == G4Proton::Proton() ) pPDG=2212;
98	else if(incidentParticleDefinition == G4Neutron::Neutron()) pPDG=2112;
99	else G4cout<<"G4QDiffraction::GetMeanFreePath: only nA & pA are implemented"<<G4endl;
100
101	G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
102	G4double sigma=0.; // Sums over elements for the material
103	G4int IPIE=IsoProbInEl.size(); // How many old elements?
104	if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
105	{
106	std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
107	SPI->clear();
108	delete SPI;
109	std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
110	IsN->clear();
111	delete IsN;
112	}
113	ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
114	ElementZ.clear(); // Clear the body vector for Z of Elements
115	IsoProbInEl.clear(); // Clear the body vector for SPI
116	ElIsoN.clear(); // Clear the body vector for N of Isotopes
117	for(G4int i=0; i<nE; ++i)
118	{
119	G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
120	G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
121	ElementZ.push_back(Z); // Remember Z of the Element
122	G4int isoSize=0; // The default for the isoVectorLength is 0
123	G4int indEl=0; // Index of non-natural element or 0(default)
124	G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
125	if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
126	#ifdef debug
127	G4cout<<"G4QDiffraction::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
128	#endif
129	if(isoSize) // The Element has non-trivial abundance set
130	{
131	indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order
132	#ifdef debug
133	G4cout<<"G4QDiffr::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
134	#endif
135	if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
136	{
137	std::vector<std::pair<G4int,G4double>> newAbund =
138	new std::vector<std::pair<G4int,G4double>*>;
139	G4double* abuVector=pElement->GetRelativeAbundanceVector();
140	for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
141	{
142	G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
143	if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
144	<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
145	G4double abund=abuVector[j];
146	std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
147	#ifdef debug
148	G4cout<<"G4QDiffract::GetMeanFreePath:pair#"<<j<<",N="<<N<<",ab="<<abund<<G4endl;
149	#endif
150	newAbund->push_back(pr);
151	}
152	#ifdef debug
153	G4cout<<"G4QDiffract::GetMeanFreePath:pairVectorLength="<<newAbund->size()<<G4endl;
154	#endif
155	indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
156	for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
157	delete newAbund; // Was "new" in the beginning of the name space
158	}
159	}
160	std::vector<std::pair<G4int,G4double>> cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
161	std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
162	IsoProbInEl.push_back(SPI);
163	std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
164	ElIsoN.push_back(IsN);
165	G4int nIs=cs->size(); // A#Of Isotopes in the Element
166	#ifdef debug
167	G4cout<<"G4QDiffract::GetMFP:=***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
168	#endif
169	G4double susi=0.; // sum of CS over isotopes
170	if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
171	{
172	std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
173	G4int N=curIs->first; // #of Neuterons in the isotope j of El i
174	IsN->push_back(N); // Remember Min N for the Element
175	#ifdef debug
176	G4cout<<"G4QDiff::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
177	#endif
178	G4bool ccsf=true;
179	if(Q==-27.) ccsf=false;
180	#ifdef debug
181	G4cout<<"G4QDiffraction::GMFP: GetCS #1 j="<<j<<G4endl;
182	#endif
183	G4double CSI=CalculateXS(Momentum, Z, N, pPDG); // XS(j,i) for theIsotope
184
185	#ifdef debug
186	G4cout<<"G4QDiffraction::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
187	<<Momentu<<", XSec="<<CSI/millibarn<<G4endl;
188	#endif
189	curIs->second = CSI;
190	susi+=CSI; // Make a sum per isotopes
191	SPI->push_back(susi); // Remember summed cross-section
192	} // End of temporary initialization of the cross sections in the G4QIsotope singeltone
193	sigma+=Isotopes->GetMeanCrossSection(Z,indEl)NOfNucPerVolume[i];//SUM(MeanCSNOfNperV)
194	#ifdef debug
195	G4cout<<"G4QDiffraction::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
196	<<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
197	#endif
198	ElProbInMat.push_back(sigma);
199	} // End of LOOP over Elements
200	// Check that cross section is not zero and return the mean free path
201	#ifdef debug
202	G4cout<<"G4QDiffraction::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
203	#endif
204	if(sigma > 0.) return 1./sigma; // Mean path [distance]
205	return DBL_MAX;
206	}
207
208	G4bool G4QDiffraction::IsApplicable(const G4ParticleDefinition& particle)
209	{
210	if (particle == *( G4Proton::Proton() )) return true;
211	else if (particle == *( G4Neutron::Neutron() )) return true;
212	//else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
213	//else if (particle == *( G4TauPlus::TauPlus() )) return true;
214	//else if (particle == *( G4TauMinus::TauMinus() )) return true;
215	//else if (particle == *( G4Electron::Electron() )) return true;
216	//else if (particle == *( G4Positron::Positron() )) return true;
217	//else if (particle == *( G4Gamma::Gamma() )) return true;
218	//else if (particle == *( G4MuonPlus::MuonPlus() )) return true;
219	//else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
220	//else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
221	//else if (particle == *( G4PionMinus::PionMinus() )) return true;
222	//else if (particle == *( G4PionPlus::PionPlus() )) return true;
223	//else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
224	//else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
225	//else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
226	//else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
227	//else if (particle == *( G4Lambda::Lambda() )) return true;
228	//else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
229	//else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
230	//else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
231	//else if (particle == *( G4XiMinus::XiMinus() )) return true;
232	//else if (particle == *( G4XiZero::XiZero() )) return true;
233	//else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
234	//else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
235	//else if (particle == *( G4AntiProton::AntiProton() )) return true;
236	#ifdef debug
237	G4cout<<"***>>G4QDiffraction::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<G4endl;
238	#endif
239	return false;
240	}
241
242	G4VParticleChange* G4QDiffraction::PostStepDoIt(const G4Track& track, const G4Step& step)
243	{
244	static const G4double mProt= G4QPDGCode(2212).GetMass(); // CHIPS proton Mass in MeV
245	static const G4double mNeut= G4QPDGCode(2112).GetMass(); // CHIPS neutron Mass in MeV
246	static const G4double mPion= G4QPDGCode(111).GetMass(); // CHIPS Pi0 Mass in MeV
247	static G4QDiffractionRatio* diffRatio;
248	//
249	//-------------------------------------------------------------------------------------
250	static G4bool CWinit = true; // CHIPS Warld needs to be initted
251	if(CWinit)
252	{
253	CWinit=false;
254	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
255	diffRatio=G4QDiffractionRatio::GetPointer();
256	}
257	//-------------------------------------------------------------------------------------
258	const G4DynamicParticle* projHadron = track.GetDynamicParticle();
259	const G4ParticleDefinition* particle=projHadron->GetDefinition();
260	#ifdef debug
261	G4cout<<"G4QDiffraction::PostStepDoIt: Before the GetMeanFreePath is called In4M="
262	<<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
263	<<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
264	#endif
265	G4ForceCondition cond=NotForced;
266	GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters?
267	#ifdef debug
268	G4cout<<"G4QDiffraction::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
269	#endif
270	G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
271	G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
272	G4double Momentum = proj4M.rho(); // @@ Just for the test purposes
273	if(std::fabs(Momentum-momentum)>.000001)
274	G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
275	#ifdef pdebug
276	G4cout<<"G4QDiffraction::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum
277	<<",proj4M="<<proj4M<<", projM="<<proj4M.m()<<G4endl;
278	#endif
279	if (!IsApplicable(*particle)) // Check applicability
280	{
281	G4cerr<<"G4QDiffraction::PostStepDoIt: Only NA is implemented."<<G4endl;
282	return 0;
283	}
284	const G4Material* material = track.GetMaterial(); // Get the current material
285	G4int Z=0;
286	const G4ElementVector* theElementVector = material->GetElementVector();
287	G4int nE=material->GetNumberOfElements();
288	#ifdef debug
289	G4cout<<"G4QDiffraction::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
290	#endif
291	G4int projPDG=0; // PDG Code prototype for the captured hadron
292	// Not all these particles are implemented yet (see Is Applicable)
293	if (particle == G4Proton::Proton() ) projPDG= 2212;
294	else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
295	//else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
296	//else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
297	//else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321;
298	//else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
299	//else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
300	//else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310;
301	//else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
302	//else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
303	//else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14;
304	//else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14;
305	//else if (particle == G4Electron::Electron() ) projPDG= 11;
306	//else if (particle == G4Positron::Positron() ) projPDG= -11;
307	//else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12;
308	//else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12;
309	//else if (particle == G4Gamma::Gamma() ) projPDG= 22;
310	//else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
311	//else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
312	//else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
313	//else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
314	//else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
315	//else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
316	//else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
317	//else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
318	//else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
319	//else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
320	//else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
321	//else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
322	//else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
323	#ifdef debug
324	G4int prPDG=particle->GetPDGEncoding();
325	G4cout<<"G4QDiffraction::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
326	#endif
327	if(!projPDG)
328	{
329	G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
330	return 0;
331	}
332	//G4double pM2=proj4M.m2(); // in MeV^2
333	//G4double pM=std::sqrt(pM2); // in MeV
334	G4double pM=mNeut;
335	G4int fPDG=2112;
336	if(projPDG==2112)
337	{
338	pM=mProt;
339	fPDG=2212;
340	}
341	// Element treatment
342	G4int EPIM=ElProbInMat.size();
343	#ifdef debug
344	G4cout<<"G4QDiffra::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
345	#endif
346	G4int i=0;
347	if(EPIM>1)
348	{
349	G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
350	for(i=0; i<nE; ++i)
351	{
352	#ifdef debug
353	G4cout<<"G4QDiffra::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
354	#endif
355	if (rnd<ElProbInMat[i]) break;
356	}
357	if(i>=nE) i=nE-1; // Top limit for the Element
358	}
359	G4Element* pElement=(*theElementVector)[i];
360	Z=static_cast<G4int>(pElement->GetZ());
361	#ifdef debug
362	G4cout<<"G4QDiffraction::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
363	#endif
364	if(Z<=0)
365	{
366	G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Element with Z="<<Z<<G4endl;
367	if(Z<0) return 0;
368	}
369	std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
370	std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
371	G4int nofIsot=SPI->size(); // #of isotopes in the element i
372	#ifdef debug
373	G4cout<<"G4QDiffraction::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
374	#endif
375	G4int j=0;
376	if(nofIsot>1)
377	{
378	G4double rndI=(SPI)[nofIsot-1]G4UniformRand(); // Randomize the isotop of the Element
379	for(j=0; j<nofIsot; ++j)
380	{
381	#ifdef debug
382	G4cout<<"G4QDiffraction::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
383	#endif
384	if(rndI < (*SPI)[j]) break;
385	}
386	if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
387	}
388	G4int N =(*IsN)[j]; ; // Randomized number of neutrons
389	#ifdef debug
390	G4cout<<"G4QDiffraction::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<", MeV="<<MeV<<G4endl;
391	#endif
392	if(N<0)
393	{
394	G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Isotope Z="<<Z<<" has 0>N="<<N<<G4endl;
395	return 0;
396	}
397	nOfNeutrons=N; // Remember it for the energy-momentum check
398	#ifdef debug
399	G4cout<<"G4QDiffraction::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
400	#endif
401	if(N<0)
402	{
403	G4cerr<<"WarningG4QDiffraction::PostStepDoIt:Element with N="<<N<< G4endl;
404	return 0;
405	}
406	aParticleChange.Initialize(track);
407	#ifdef debug
408	G4cout<<"G4QDiffraction::PostStepDoIt: track is initialized"<<G4endl;
409	#endif
410	G4double weight = track.GetWeight();
411	G4double localtime = track.GetGlobalTime();
412	G4ThreeVector position = track.GetPosition();
413	#ifdef debug
414	G4cout<<"G4QDiffraction::PostStepDoIt: before Touchable extraction"<<G4endl;
415	#endif
416	G4TouchableHandle trTouchable = track.GetTouchableHandle();
417	#ifdef debug
418	G4cout<<"G4QDiffraction::PostStepDoIt: Touchable is extracted"<<G4endl;
419	#endif
420	G4int targPDG=90000000+Z*1000+N; // CHIPS PDG Code of the target nucleus
421	G4QPDGCode targQPDG(targPDG); // @@ use G4Ion and get rid of CHIPS World
422	G4double tM=targQPDG.GetMass(); // CHIPS final nucleus mass in MeV
423	G4double kinEnergy= projHadron->GetKineticEnergy()MeV; // Kin energy in MeV (Is MeV n?)
424	G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
425	G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
426	#ifdef debug
427	G4cout<<"G4QDiffraction::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
428	#endif
429	EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation
430	// @@ Probably this is not necessary any more
431	#ifdef debug
432	G4cout<<"G4QDiff::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl;
433	#endif
434	G4double xSec=CalculateXS(Momentum, Z, N, projPDG); // Recalculate CrossSection
435	#ifdef debug
436	G4cout<<"G4QDiffra::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",XS="<<xSec/millibarn<<G4endl;
437	#endif
438	#ifdef nandebug
439	if(xSec>0. \|\| xSec<0. \|\| xSec==0);
440	else G4cout<<"-Warning-G4QDiffraction::PostSDI: NAN xSec="<<xSec/millibarn<<G4endl;
441	#endif
442	// @@ check a possibility to separate p, n, or alpha (!)
443	if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
444	{
445	#ifdef pdebug
446	G4cerr<<"WarningG4QDiffraction::PSDoIt:Zero cross-section PDG="<<projPDG
447	<<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl;
448	#endif
449	//Do Nothing Action insead of the reaction
450	aParticleChange.ProposeEnergy(kinEnergy);
451	aParticleChange.ProposeLocalEnergyDeposit(0.);
452	aParticleChange.ProposeMomentumDirection(dir) ;
453	return G4VDiscreteProcess::PostStepDoIt(track,step);
454	}
455	G4double totCMMass=tot4M.m(); // Total CM mass, pM=projectileMass, tM=targetMass
456	if(totCMMass < mPion+pM+tM) // The diffraction reaction is impossible -> Do Nothing
457	{
458	#ifdef pdebug
459	G4cerr<<"WarningG4QDiffraction::PSDoIt:Below Diffraction Threshold cmM="<<totCMMass
460	<<">pM="<<pM<<"+tM="<<tM<<"+pi0="<<mPion<<"=="<<pM+tM+mPion<<G4endl;
461	#endif
462	//Do Nothing Action insead of the reaction
463	aParticleChange.ProposeEnergy(kinEnergy);
464	aParticleChange.ProposeLocalEnergyDeposit(0.);
465	aParticleChange.ProposeMomentumDirection(dir) ;
466	return G4VDiscreteProcess::PostStepDoIt(track,step);
467	}
468	// Kill interacting hadron
469	aParticleChange.ProposeTrackStatus(fStopAndKill);
470	G4QHadronVector* out=diffRatio->TargFragment(projPDG, proj4M, Z, N);
471	G4int nSec=out->size(); // #of secondaries in the diffraction reaction
472	G4DynamicParticle* theSec=0; // A prototype for secondary for the secondary
473	G4LorentzVector dif4M(0.,0.,0.,0.); // Prototype for the secondary 4-momentum
474	G4int difPDG=0; // PDG code of the secondary
475	G4QHadron* difQH=0; // Prototype for a Q-secondary
476	#ifdef pdebug
477	G4cout<<"G4QDiffraction::PostStepDoIt: =====found===== nSecondaries="<<nSec<<G4endl;
478	#endif
479	for(G4int i=0; i<nSec; i++)
480	{
481	difQH = (*out)[i];
482	difPDG= difQH->GetPDGCode();
483	G4ParticleDefinition* theDefinition=0;
484	if (difPDG==2212 \|\| difPDG==90001000) theDefinition=G4Proton::Proton();
485	else if(difPDG==2112 \|\| difPDG==90000001) theDefinition=G4Neutron::Neutron();
486	else if(difPDG== 22) theDefinition=G4Gamma::Gamma();
487	else if(difPDG== 111) theDefinition=G4PionZero::PionZero();
488	else if(difPDG==-211 \|\| difPDG==89999001) theDefinition=G4PionMinus::PionMinus();
489	else if(difPDG== 211 \|\| difPDG==90000999) theDefinition=G4PionPlus::PionPlus();
490	else if(difPDG== 321 \|\| difPDG==89001000) theDefinition=G4KaonPlus::KaonPlus();
491	else if(difPDG==-321 \|\| difPDG==90999000) theDefinition=G4KaonMinus::KaonMinus();
492	else if(difPDG== 130 \|\| difPDG==-311 \|\| difPDG==89000001)
493	theDefinition=G4KaonZeroLong::KaonZeroLong();
494	else if(difPDG== 310 \|\| difPDG== 311 \|\| difPDG==90999999)
495	theDefinition=G4KaonZeroShort::KaonZeroShort();
496	else if(difPDG==3122 \|\| difPDG==91000000) theDefinition=G4Lambda::Lambda();
497	else if(difPDG== 3222) theDefinition=G4SigmaPlus::SigmaPlus();
498	else if(difPDG== 3112) theDefinition=G4SigmaMinus::SigmaMinus();
499	else if(difPDG== 3212) theDefinition=G4SigmaZero::SigmaZero();
500	else if(difPDG== 3312) theDefinition=G4XiMinus::XiMinus();
501	else if(difPDG== 3322) theDefinition=G4XiZero::XiZero();
502	else if(difPDG== 3334) theDefinition=G4OmegaMinus::OmegaMinus();
503	else if(difPDG==-2112) theDefinition=G4AntiNeutron::AntiNeutron();
504	else if(difPDG==-2212) theDefinition=G4AntiProton::AntiProton();
505	else if(difPDG==-3122) theDefinition=G4AntiLambda::AntiLambda();
506	else if(difPDG==-3222) theDefinition=G4AntiSigmaPlus::AntiSigmaPlus();
507	else if(difPDG==-3112) theDefinition=G4AntiSigmaMinus::AntiSigmaMinus();
508	else if(difPDG==-3212) theDefinition=G4AntiSigmaZero::AntiSigmaZero();
509	else if(difPDG==-3312) theDefinition=G4AntiXiMinus::AntiXiMinus();
510	else if(difPDG==-3322) theDefinition=G4AntiXiZero::AntiXiZero();
511	else if(difPDG==-3334) theDefinition=G4AntiOmegaMinus::AntiOmegaMinus();
512	else if(difPDG== -11) theDefinition=G4Electron::Electron();
513	else if(difPDG== -13) theDefinition=G4MuonMinus::MuonMinus();
514	else if(difPDG== 11) theDefinition=G4Positron::Positron();
515	else if(difPDG== 13) theDefinition=G4MuonPlus::MuonPlus();
516	else
517	{
518	G4int Z = difQH->GetCharge();
519	G4int B = difQH->GetBaryonNumber();
520	G4int S = difQH->GetStrangeness();
521	if(S\|\|Z>B\|\|Z<0)G4cout<<"-Warning-G4QDif::PoStDoIt:Z="<<Z<<",A="<<B<<",S="<<S<<G4endl;
522	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(Z,B,0,0);
523	#ifdef pdebug
524	G4cout<<"G4QDiffraction::PoStDoIt:Ion,Z="<<Z<<",A="<<B<<",D="<<theDefinition<<G4endl;
525	#endif
526	}
527	if(theDefinition)
528	{
529	theSec = new G4DynamicParticle; // A secondary for the recoil hadron
530	theSec->SetDefinition(theDefinition);
531	dif4M = difQH->Get4Momentum();
532	EnMomConservation-=dif4M;
533	theSec->Set4Momentum(dif4M);
534	G4Track* aNewTrack = new G4Track(theSec, localtime, position );
535	aNewTrack->SetWeight(weight); // weighted
536	aNewTrack->SetTouchableHandle(trTouchable);
537	aParticleChange.AddSecondary( aNewTrack );
538	#ifdef pdebug
539	G4cout<<"G4QDiffraction::PostStepDoIt: Filled 4M="<<dif4M<<", PDG="<<difPDG<<G4endl;
540	#endif
541	}
542	else G4cout<<"-Warning-G4QDif::PSDI: Lost PDG="<<difPDG<<", Z="<<difQH->GetCharge()
543	<<", A="<<difQH->GetBaryonNumber()<<",S ="<<difQH->GetStrangeness()<<G4endl;
544	delete difQH; // Clean up the output QHadrons
545	}
546	delete out; // Delete the output QHadron-vector
547	#ifdef debug
548	G4cout<<"G4QDiffraction::PostStepDoIt:* PostStepDoIt is done *"<<G4endl;
549	#endif
550	return G4VDiscreteProcess::PostStepDoIt(track, step);
551	}
552
553	G4double G4QDiffraction::CalculateXS(G4double p, G4int Z, G4int N, G4int PDG)
554	{
555	static G4bool first=true;
556	static G4VQCrossSection* CSmanager;
557	static G4QDiffractionRatio* diffRatio;
558	if(first) // Connection with a singletone
559	{
560	CSmanager=G4QProtonNuclearCrossSection::GetPointer();
561	diffRatio=G4QDiffractionRatio::GetPointer();
562	first=false;
563	}
564	//G4double x=CSmanager->GetCrossSection(true, p, Z, N, PDG); // inelastic XS
565	//G4double pIU=p*GeV; // IndependentUnistMomentum
566	//G4double r=diffRatio->GetRatio(pIU, PDG, Z, N); // Proj. Diffraction Part
567	//G4double s=x*r; // XS for proj. diffraction
568	G4double s=diffRatio->GetTargSingDiffXS(p, PDG, Z, N); // XS for target diffraction
569	#ifdef debug
570	G4cout<<"G4QDiff::CXS:p="<<p<<",Z="<<Z<<",N="<<N<<",C="<<PDG<<",XS="<<s<<G4endl;
571	#endif
572	return s;
573	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: