Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QGluonString.cc @ 1007

Last change on this file since 1007 was 1007, checked in by garnier, 15 years ago
update to geant4.9.2
File size: 30.0 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: G4QGluonString.cc,v 1.3 2008/10/02 21:10:07 dennis Exp $
27	// GEANT4 tag $Name: geant4-09-02 $
28	//
29	// ---------------- G4QGluonString class -----------------
30	// by Mikhail Kossov, December 2003.
31	// G4QGluonString class of the CHIPS Simulation Branch in GEANT4
32	// ---------------------------------------------------------------
33	// ****************************************************************************************
34	// ******** This CLASS is temporary moved from the photolepton_hadron directory *******
35	// ****************************************************************************************
36
37	//#define debug
38	//#define pdebug
39	//#define ppdebug
40
41	#include "G4QGluonString.hh"
42
43	// Initialization of static vectors
44	std::vector<G4int> G4QGluonString::ElementZ; // Z of the element(i) in theLastCalc
45	std::vector<G4double> G4QGluonString::ElProbInMat; // SumProbabilityElements in Material
46	std::vector<std::vector<G4int>*> G4QGluonString::ElIsoN; // N of isotope(j) of Element(i)
47	std::vector<std::vector<G4double>*>G4QGluonString::IsoProbInEl;//SumProbabIsotopeInElementI
48
49	G4QGluonString::G4QGluonString(const G4String& processName):
50	G4VDiscreteProcess(processName, fHadronic)
51	{
52	#ifdef debug
53	G4cout<<"G4QGluonString::Constructor is called"<<G4endl;
54	#endif
55	if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created by CHIPS"<<G4endl;
56	SetProcessSubType(fHadronInelastic);
57	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
58	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
59	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
60	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
61	//@@ Initialize here other parameters
62	}
63
64	G4bool G4QGluonString::manualFlag=false; // If false then standard parameters are used
65	G4double G4QGluonString::Temperature=180.; // Critical Temperature (sensitive at High En)
66	G4double G4QGluonString::SSin2Gluons=0.3; // Supression of s-quarks (in respect to u&d)
67	G4double G4QGluonString::EtaEtaprime=0.3; // Supression of eta mesons (gg->qq/3g->qq)
68	G4double G4QGluonString::freeNuc=0.5; // Percentage of free nucleons on the surface
69	G4double G4QGluonString::freeDib=0.05; // Percentage of free diBaryons on the surface
70	G4double G4QGluonString::clustProb=5.; // Nuclear clusterization parameter
71	G4double G4QGluonString::mediRatio=10.; // medium/vacuum hadronization ratio
72	G4int G4QGluonString::nPartCWorld=152; // The#of particles initialized in CHIPS World
73	G4double G4QGluonString::SolidAngle=0.5; // Part of Solid Angle to capture (@@A-dep.)
74	G4bool G4QGluonString::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
75	G4double G4QGluonString::PiPrThresh=141.4; // Pion Production Threshold for gammas
76	G4double G4QGluonString::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
77	G4double G4QGluonString::DiNuclMass=1880.; // DoubleNucleon Mass for VirtualNormalization
78
79	void G4QGluonString::SetManual() {manualFlag=true;}
80	void G4QGluonString::SetStandard() {manualFlag=false;}
81
82	// Fill the private parameters
83	void G4QGluonString::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
84	G4double fN, G4double fD, G4double cP, G4double mR,
85	G4int nParCW, G4double solAn, G4bool efFlag,
86	G4double piThresh, G4double mpisq, G4double dinum)
87	{// =============================================================================
88	Temperature=temper;
89	SSin2Gluons=ssin2g;
90	EtaEtaprime=etaetap;
91	freeNuc=fN;
92	freeDib=fD;
93	clustProb=cP;
94	mediRatio=mR;
95	nPartCWorld = nParCW;
96	EnergyFlux=efFlag;
97	SolidAngle=solAn;
98	PiPrThresh=piThresh;
99	M2ShiftVir=mpisq;
100	DiNuclMass=dinum;
101	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
102	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
103	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
104	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
105	}
106
107	// Destructor
108
109	G4QGluonString::~G4QGluonString() {}
110
111	// Internal E/P conservation 4-Mom & the selected number of neutrons in the Element
112
113	G4LorentzVector G4QGluonString::GetEnegryMomentumConservation()
114	{
115	return EnMomConservation;
116	}
117
118	G4int G4QGluonString::GetNumberOfNeutronsInTarget()
119	{
120	return nOfNeutrons;
121	}
122
123	// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
124	// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
125	// ******** All CHIPS cross sections are calculated in the surface units **********
126	G4double G4QGluonString::GetMeanFreePath(const G4Track& aTrack,
127	G4double, G4ForceCondition* Fc)
128	{
129	#ifdef debug
130	G4cout<<"G4QGluonString::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
131	#endif
132	*Fc = NotForced;
133	#ifdef debug
134	G4cout<<"G4QGluonString::GetMeanFreePath: Before GetDynPart"<<G4endl;
135	#endif
136	const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
137	#ifdef debug
138	G4cout<<"G4QGluonString::GetMeanFreePath: Before GetDef"<<G4endl;
139	#endif
140	G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
141	if( !IsApplicable(*incidentParticleDefinition))
142	G4cout<<"-W-G4QGluonString::GetMeanFreePath called for NotImplementedParticle"<<G4endl;
143	// Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
144	G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
145	#ifdef debug
146	G4cout<<"G4QCollis::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
147	#endif
148	const G4Material* material = aTrack.GetMaterial(); // Get the current material
149	const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
150	const G4ElementVector* theElementVector = material->GetElementVector();
151	G4int nE=material->GetNumberOfElements();
152	#ifdef debug
153	G4cout<<"G4QGluonString::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
154	#endif
155	G4VQCrossSection* CSmanager=0;
156	G4int pPDG=0;
157	if(incidentParticleDefinition == G4Proton::Proton())
158	{
159	CSmanager=G4QProtonNuclearCrossSection::GetPointer();
160	pPDG=2212;
161	} //@@ Make cross-section mahnagers for other mesons & baryons
162	else
163	{
164	G4cerr<<"***G4QGluonString::GetMeanFreePath: Particle isn't implemented"<<G4endl;
165	G4Exception("G4QGluonString::PostStepDoIt:","72",FatalException,"BadProjectile");
166	}
167
168	G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
169	G4double sigma=0.; // Sums over elements for the material
170	G4int IPIE=IsoProbInEl.size(); // How many old elements?
171	if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
172	{
173	std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
174	SPI->clear();
175	delete SPI;
176	std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
177	IsN->clear();
178	delete IsN;
179	}
180	ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
181	ElementZ.clear(); // Clear the body vector for Z of Elements
182	IsoProbInEl.clear(); // Clear the body vector for SPI
183	ElIsoN.clear(); // Clear the body vector for N of Isotopes
184	for(G4int i=0; i<nE; ++i)
185	{
186	G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
187	G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
188	ElementZ.push_back(Z); // Remember Z of the Element
189	G4int isoSize=0; // The default for the isoVectorLength is 0
190	G4int indEl=0; // Index of non-trivial element or 0(default)
191	G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
192	if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
193	#ifdef debug
194	G4cout<<"G4QGluonString::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
195	#endif
196	if(isoSize) // The Element has non-trivial abumdance set
197	{
198	indEl=pElement->GetIndex(); // Index of the non-trivial element
199	if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
200	{
201	std::vector<std::pair<G4int,G4double>> newAbund =
202	new std::vector<std::pair<G4int,G4double>*>;
203	G4double* abuVector=pElement->GetRelativeAbundanceVector();
204	for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
205	{
206	G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
207	if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QCaptureAtRest::GetMeanFreePath"
208	<<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
209	G4double abund=abuVector[j];
210	std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
211	#ifdef debug
212	G4cout<<"G4QGluonString::GetMeanFreePath:p#"<<j<<",N="<<N<<",ab="<<abund<<G4endl;
213	#endif
214	newAbund->push_back(pr);
215	}
216	#ifdef debug
217	G4cout<<"G4QGluonString::PostStepDoIt:pairVectorLength="<<newAbund->size()<<G4endl;
218	#endif
219	indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
220	for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
221	delete newAbund; // Was "new" in the beginning of the name space
222	}
223	}
224	std::vector<std::pair<G4int,G4double>> cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
225	std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
226	IsoProbInEl.push_back(SPI);
227	std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
228	ElIsoN.push_back(IsN);
229	G4int nIs=cs->size(); // A#Of Isotopes in the Element
230	G4double susi=0.; // sum of CS over isotopes
231	if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
232	{
233	std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
234	G4int N=curIs->first; // #of Neuterons in the isotope j of El i
235	IsN->push_back(N); // Remember Min N for the Element
236	G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
237	#ifdef debug
238	G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
239	#endif
240	curIs->second = CSI; // Remenber the calculated cross-section
241	susi+=CSI; // Make a sum per isotopes
242	SPI->push_back(susi); // Remember summed cross-section
243	} // End of temporary initialization of the cross sections in the G4QIsotope singeltone
244	sigma+=Isotopes->GetMeanCrossSection(Z,indEl)NOfNucPerVolume[i];//SUM(MeanCSNOfNperV)
245	ElProbInMat.push_back(sigma);
246	} // End of LOOP over Elements
247	#ifdef debug
248	G4cout<<"G4QCol::GetMeanFrPa: S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
249	#endif
250	// Check that cross section is not zero and return the mean free path
251	if(sigma > 0.) return 1./sigma; // Mean path [distance]
252	return DBL_MAX; // If Sigma=0, return max value for PATH
253	}
254
255	// Check applicability of the process
256	G4bool G4QGluonString::IsApplicable(const G4ParticleDefinition& particle)
257	{
258	if (particle == *( G4Proton::Proton() )) return true;
259	//else if (particle == *( G4Neutron::Neutron() )) return true;
260	//else if (particle == *( G4PionMinus::PionMinus() )) return true;
261	//else if (particle == *( G4PionPlus::PionPlus() )) return true;
262	//else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
263	//else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
264	//else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
265	//else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
266	//else if (particle == *( G4Lambda::Lambda() )) return true;
267	//else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
268	//else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
269	//else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
270	//else if (particle == *( G4XiMinus::XiMinus() )) return true;
271	//else if (particle == *( G4XiZero::XiZero() )) return true;
272	//else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
273	//else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
274	//else if (particle == *( G4AntiProton::AntiProton() )) return true;
275	//else if (particle == *( G4AntiLambda::AntiLambda() )) return true;
276	//else if (particle == *( G4AntiSigmaPlus::AntiSigmaPlus() )) return true;
277	//else if (particle == *(G4AntiSigmaMinus::AntiSigmaMinus())) return true;
278	//else if (particle == *( G4AntiSigmaZero::AntiSigmaZero() )) return true;
279	//else if (particle == *( G4AntiXiMinus::AntiXiMinus() )) return true;
280	//else if (particle == *( G4AntiXiZero::AntiXiZero() )) return true;
281	//else if (particle == *(G4AntiOmegaMinus::AntiOmegaMinus())) return true;
282	else G4cerr<<"***G4QGluonString::IsApplicable: PDG?="<<particle.GetPDGEncoding()<<G4endl;
283	#ifdef debug
284	G4cout<<"***G4QGluonString::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
285	#endif
286	return false;
287	}
288
289	G4VParticleChange* G4QGluonString::PostStepDoIt(const G4Track& track, const G4Step& step)
290	{
291	//static const G4double dpi=M_PI+M_PI; // 2pi (for Phi distr.) changed to twopi*
292	//static const G4double mNeut= G4QPDGCode(2112).GetMass();
293	//static const G4double mNeut2= mNeut*mNeut; // Squared neutron mass
294	//static const G4double mProt= G4QPDGCode(2212).GetMass();
295	//static const G4double mProt2= mProt*mProt; // Squared Proton mass
296	//static const G4double dM=mProt+mNeut; // doubled nucleon mass
297	//static const G4double hdM=dM/2.; // M of the "nucleon"
298	//static const G4double hdM2=hdM*hdM; // M2 of the "nucleon"
299	//static const G4double mPi0 = G4QPDGCode(111).GetMass();
300	//static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
301	//static const G4double mPi = G4QPDGCode(211).GetMass();
302	//static const G4double tmPi = mPi+mPi; // DoubledMass of the charged pion
303	//static const G4double stmPi= tmPi*tmPi; // SquareDoubledMass of ChargedPion
304	//static const G4double mPPi = mPi+mProt; // Delta threshold
305	//static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2
306	//-------------------------------------------------------------------------------------
307	const G4DynamicParticle* projHadron = track.GetDynamicParticle();
308	const G4ParticleDefinition* particle=projHadron->GetDefinition();
309	#ifdef debug
310	G4cout<<"G4QGluonString::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
311	#endif
312	G4ForceCondition cond=NotForced;
313	GetMeanFreePath(track, 1., &cond); // Just to check that still sig>0
314	#ifdef debug
315	G4cout<<"G4QGluonString::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
316	#endif
317	G4LorentzVector proj4M=projHadron->Get4Momentum();
318	G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
319	G4double Momentum=proj4M.rho();
320	if(std::fabs(Momentum-momentum)>.001)
321	G4cerr<<"G4QGluonString::PostStepDoIt: P="<<Momentum<<"="<<momentum<<G4endl;
322	#ifdef debug
323	G4double mp=proj4M.m();
324	G4cout<<"G4QGluonString::PostStepDoIt is called, P="<<Momentum<<"="<<momentum<<G4endl;
325	#endif
326	if (!IsApplicable(*particle)) // Check applicability
327	{
328	G4cerr<<"G4QGluonString::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
329	<<G4endl;
330	return 0;
331	}
332	const G4Material* material = track.GetMaterial(); // Get the current material
333	G4int Z=0;
334	const G4ElementVector* theElementVector = material->GetElementVector();
335	G4int nE=material->GetNumberOfElements();
336	#ifdef debug
337	G4cout<<"G4QGluonString::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
338	#endif
339	G4int projPDG=0; // PDG Code prototype for the captured hadron
340	// Not all these particles are implemented yet (see Is Applicable)
341	if (particle == G4Proton::Proton() ) projPDG= 2212;
342	//else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
343	//else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
344	//else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
345	//else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 2112;
346	//else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
347	//else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
348	//else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310;
349	//else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
350	//else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
351	//else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
352	//else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
353	//else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
354	//else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
355	//else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
356	//else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
357	//else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
358	//else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
359	//else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
360	//else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
361	//else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
362	//else if (particle == G4AntiLambda::AntiLambda() ) projPDG=-3122;
363	//else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) projPDG=-3222;
364	//else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112;
365	//else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) projPDG=-3212;
366	//else if (particle == G4AntiXiMinus::AntiXiMinus() ) projPDG=-3312;
367	//else if (particle == G4AntiXiZero::AntiXiZero() ) projPDG=-3322;
368	//else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334;
369	//G4int aProjPDG=std::abs(projPDG);
370	#ifdef debug
371	G4int prPDG=particle->GetPDGEncoding();
372	G4cout<<"G4QGluonString::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
373	#endif
374	if(!projPDG)
375	{
376	G4cerr<<"--Warning--G4QGluonString::PostStepDoIt:Undefined interacting hadron"<<G4endl;
377	return 0;
378	}
379	G4int EPIM=ElProbInMat.size();
380	#ifdef debug
381	G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
382	#endif
383	G4int i=0;
384	if(EPIM>1)
385	{
386	G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
387	for(i=0; i<nE; ++i)
388	{
389	#ifdef debug
390	G4cout<<"G4QGluonString::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
391	#endif
392	if (rnd<ElProbInMat[i]) break;
393	}
394	if(i>=nE) i=nE-1; // Top limit for the Element
395	}
396	G4Element* pElement=(*theElementVector)[i];
397	Z=static_cast<G4int>(pElement->GetZ());
398	#ifdef debug
399	G4cout<<"G4QGluonString::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
400	#endif
401	if(Z<=0)
402	{
403	G4cerr<<"---Warning---G4QGluonString::PostStepDoIt: Element with Z="<<Z<<G4endl;
404	if(Z<0) return 0;
405	}
406	std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
407	std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
408	G4int nofIsot=SPI->size(); // #of isotopes in the element i
409	#ifdef debug
410	G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
411	#endif
412	G4int j=0;
413	if(nofIsot>1)
414	{
415	G4double rndI=(SPI)[nofIsot-1]G4UniformRand(); // Randomize the isotop of the Element
416	for(j=0; j<nofIsot; ++j)
417	{
418	#ifdef debug
419	G4cout<<"G4QGluonString::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
420	#endif
421	if(rndI < (*SPI)[j]) break;
422	}
423	if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
424	}
425	G4int N =(*IsN)[j]; ; // Randomized number of neutrons
426	#ifdef debug
427	G4cout<<"G4QGluonString::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<G4endl;
428	#endif
429	if(N<0)
430	{
431	G4cerr<<"-Warning-G4QGluonString::PostStepDoIt:Isotope with N="<<Z<<"<0,Z="<<Z<<G4endl;
432	return 0;
433	}
434	nOfNeutrons=N; // Remember it for the energy-momentum check
435	G4double dd=0.025;
436	G4double am=Z+N;
437	G4double sr=std::sqrt(am);
438	G4double dsr=0.01*(sr+sr);
439	if(dsr<dd)dsr=dd;
440	if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio);// ManualPa
441	//else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP CluPars
442	//else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.); //Pi- CluPars
443	#ifdef debug
444	G4cout<<"G4QGluonString::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
445	#endif
446	if(N<0)
447	{
448	G4cerr<<"---Warning---G4QGluonString::PostStepDoIt:Element with N="<<N<< G4endl;
449	return 0;
450	}
451	aParticleChange.Initialize(track);
452	G4double localtime = track.GetGlobalTime();
453	G4ThreeVector position = track.GetPosition();
454	G4TouchableHandle trTouchable = track.GetTouchableHandle();
455	//
456	G4int targPDG=90000000+Z*1000+N; // PDG Code of the target nucleus
457	// =========================
458	G4QPDGCode targQPDG(targPDG);
459	G4double tM=targQPDG.GetMass();
460	G4QHadronVector* output=new G4QHadronVector;// Prototype of EnvironOutput G4QHadronVector
461	G4double absMom = 0.; // Prototype of absorbed by nucleus Moment
462	G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVectorum
463	G4LorentzVector lead4M(0.,0.,0.,0.); // Prototype of LeadingQ 4-momentum
464	EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
465	if(absMom) EnMomConservation+=lead4M; // Add E/M of leading System
466	#ifdef debug
467	G4cout<<"G4QGluonString::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
468	#endif
469	//G4QHadron* pH = new G4QHadron(projPDG,proj4M);
470	// @@@@@@@@@@@@@@@@@@@@@@@@@@ Hrere the CHIPS_QGS must be implemented @@@@@@@@@@@@@@@
471	G4int qNH=leadhs->size();
472	if(absMom)
473	{
474	if(qNH) for(G4int iq=0; iq<qNH; iq++)
475	{
476	G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron
477	output->push_back(loh);
478	}
479	delete leadhs;
480	}
481	// ------------- From here the secondaries are filled -------------------------
482	G4int tNH = output->size(); // A#of hadrons in the output
483	aParticleChange.SetNumberOfSecondaries(tNH);
484	// Now add nuclear fragments
485	#ifdef debug
486	G4cout<<"G4QGluonString::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
487	#endif
488	#ifdef ppdebug
489	if(absMom)G4cout<<"G4QGluonString::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
490	#endif
491	G4int nOut=output->size(); // Real length of the output @@ Temporary
492	if(tNH==1) // @@ Temporary. Find out why it happened!
493	{
494	G4cout<<"-Warning-G4QGluonString::PostStepDoIt: 1 secondary! absMom="<<absMom;
495	if(absMom) G4cout<<", qNH="<<qNH;
496	G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
497	G4cout<<G4endl;
498	tNH=0;
499	delete output->operator[](0); // delete the creazy hadron
500	output->pop_back(); // clean up the output vector
501	}
502	if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QGluonString::PostStepDoIt:2 # "<<nOut<<G4endl;
503	// Deal with ParticleChange final state interface to GEANT4 output of the process
504	//if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
505	if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
506	{
507	// Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
508	// Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
509	// The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
510	G4QHadron* hadr=(*output)[i]; // Pointer to the output hadron
511	G4int PDGCode = hadr->GetPDGCode();
512	G4int nFrag = hadr->GetNFragments();
513	#ifdef pdebug
514	G4cout<<"G4QGluonString::AtRestDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag<<G4endl;
515	#endif
516	if(nFrag) // Skip intermediate (decayed) hadrons
517	{
518	#ifdef debug
519	G4cout<<"G4QGluonString::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
520	#endif
521	delete hadr;
522	continue;
523	}
524	G4DynamicParticle* theSec = new G4DynamicParticle;
525	G4ParticleDefinition* theDefinition;
526	if (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
527	else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
528	else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
529	else if(PDGCode==311 \|\| PDGCode==-311)
530	{
531	if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L
532	else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
533	}
534	else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
535	else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
536	else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
537	else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
538	else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
539	else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
540	{
541	G4int aZ = hadr->GetCharge();
542	G4int aA = hadr->GetBaryonNumber();
543	#ifdef pdebug
544	G4cout<<"G4QGluonString::AtRestDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
545	#endif
546	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
547	}
548	//else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
549	else
550	{
551	#ifdef pdebug
552	G4cout<<"G4QGluonString::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
553	#endif
554	theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
555	#ifdef pdebug
556	G4cout<<"G4QGluonString::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
557	#endif
558	}
559	if(!theDefinition)
560	{
561	#ifdef debug
562	G4cout<<"---Warning---G4QGluonString::PostStepDoIt: drop PDG="<<PDGCode<<G4endl;
563	#endif
564	delete hadr;
565	continue;
566	}
567	#ifdef pdebug
568	G4cout<<"G4QGluonString::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
569	#endif
570	theSec->SetDefinition(theDefinition);
571	G4LorentzVector h4M=hadr->Get4Momentum();
572	EnMomConservation-=h4M;
573	#ifdef tdebug
574	G4cout<<"G4QCollis::PSDI:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
575	#endif
576	#ifdef debug
577	G4cout<<"G4QGluonString::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
578	#endif
579	theSec->Set4Momentum(h4M); // ^
580	delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
581	#ifdef debug
582	G4ThreeVector curD=theSec->GetMomentumDirection(); // ^
583	G4double curM=theSec->GetMass(); // \|
584	G4double curE=theSec->GetKineticEnergy()+curM; // ^
585	G4cout<<"G4QCollis::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;// \|
586	#endif
587	G4Track* aNewTrack = new G4Track(theSec, localtime, position ); // ^
588	aNewTrack->SetTouchableHandle(trTouchable); // \|
589	aParticleChange.AddSecondary( aNewTrack ); // \|
590	#ifdef debug
591	G4cout<<"G4QGluonString::PostStepDoIt:#"<<i<<" is done"<<G4endl;// \|
592	#endif
593	} // \|
594	delete output; // instances of the G4QHadrons from the output are already deleted above +
595	#ifdef debug
596	G4cout<<"G4QGluonString::PostStDoIt: afterSt="<<aParticleChange.GetTrackStatus()<<G4endl;
597	#endif
598	aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the absorbed particle
599	#ifdef debug
600	G4cout<<"G4QGluonString::PostStepDoIt:* PostStepDoIt is done *, P="<<aProjPDG
601	<<", St="<<aParticleChange.GetTrackStatus()<<G4endl;
602	#endif
603	return G4VDiscreteProcess::PostStepDoIt(track, step);
604	}

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

Download in other formats: