Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/processes/src/G4QInelastic.cc @ 1199

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 104.3 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: G4QInelastic.cc,v 1.2 2009/11/19 16:32:56 mkossov Exp $
27	// GEANT4 tag $Name: geant4-09-03-cand-01 $
28	//
29	// ---------------- G4QInelastic class -----------------
30	// by Mikhail Kossov, December 2003.
31	// G4QInelastic class of the CHIPS Simulation Branch in GEANT4
32	// ---------------------------------------------------------------
33	// ****************************************************************************************
34	// ********* This HEADER is a property of the CHIPS physics package (M. Kosov) ********
35	// ***** DO NOT MAKE ANY CHANGE YOURSELF! Send proposals to Mikhail.Kossov@cern.ch ****
36	// ****************************************************************************************
37	// Short description: This is a universal class for the incoherent (inelastic)
38	// nuclear interactions in the CHIPS model.
39	// ---------------------------------------------------------------------------
40	//#define debug
41	//#define pdebug
42	//#define pickupd
43	//#define ldebug
44	//#define ppdebug
45	//#define qedebug
46
47	#include "G4QInelastic.hh"
48
49	// Initialization of static vectors
50	std::vector<G4int> G4QInelastic::ElementZ; // Z of the element(i) in theLastCalc
51	std::vector<G4double> G4QInelastic::ElProbInMat; // SumProbabilityElements in Material
52	std::vector<std::vector<G4int>*> G4QInelastic::ElIsoN;// N of isotope(j) of Element(i)
53	std::vector<std::vector<G4double>*>G4QInelastic::IsoProbInEl;//SumProbabIsotopes inElementI
54
55	G4QInelastic::G4QInelastic(const G4String& processName):
56	G4VDiscreteProcess(processName, fHadronic)
57	{
58	#ifdef debug
59	G4cout<<"G4QInelastic::Constructor is called"<<G4endl;
60	#endif
61	if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
62	//G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPSWorld (234 part.max)
63	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
64	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
65	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
66	//@@ Initialize here the G4QuasmonString parameters
67	}
68
69	G4bool G4QInelastic::manualFlag=false; // If false then standard parameters are used
70	G4double G4QInelastic::Temperature=180.; // Critical Temperature (sensitive at High En)
71	G4double G4QInelastic::SSin2Gluons=0.3; // Supression of s-quarks (in respect to u&d)
72	G4double G4QInelastic::EtaEtaprime=0.3; // Supression of eta mesons (gg->qq/3g->qq)
73	G4double G4QInelastic::freeNuc=0.5; // Percentage of free nucleons on the surface
74	G4double G4QInelastic::freeDib=0.05; // Percentage of free diBaryons on the surface
75	G4double G4QInelastic::clustProb=5.; // Nuclear clusterization parameter
76	G4double G4QInelastic::mediRatio=10.; // medium/vacuum hadronization ratio
77	G4int G4QInelastic::nPartCWorld=152; // The#of particles initialized in CHIPS World
78	G4double G4QInelastic::SolidAngle=0.5; // Part of Solid Angle to capture (@@A-dep.)
79	G4bool G4QInelastic::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
80	G4double G4QInelastic::PiPrThresh=141.4; // Pion Production Threshold for gammas
81	G4double G4QInelastic::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
82	G4double G4QInelastic::DiNuclMass=1880.; // DoubleNucleon Mass for VirtualNormalization
83	G4double G4QInelastic::photNucBias=1.; // BiasingParameter for photo(e,mu,tau)Nuclear
84	G4double G4QInelastic::weakNucBias=1.; // BiasingParameter for ChargedCurrents(nu,mu)
85
86	void G4QInelastic::SetManual() {manualFlag=true;}
87	void G4QInelastic::SetStandard() {manualFlag=false;}
88
89	// Fill the private parameters
90	void G4QInelastic::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
91	G4double fN, G4double fD, G4double cP, G4double mR,
92	G4int nParCW, G4double solAn, G4bool efFlag,
93	G4double piThresh, G4double mpisq, G4double dinum)
94	{// =============================================================================
95	Temperature=temper;
96	SSin2Gluons=ssin2g;
97	EtaEtaprime=etaetap;
98	freeNuc=fN;
99	freeDib=fD;
100	clustProb=cP;
101	mediRatio=mR;
102	nPartCWorld = nParCW;
103	EnergyFlux=efFlag;
104	SolidAngle=solAn;
105	PiPrThresh=piThresh;
106	M2ShiftVir=mpisq;
107	DiNuclMass=dinum;
108	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
109	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
110	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
111	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
112	}
113
114	void G4QInelastic::SetPhotNucBias(G4double phnB) {photNucBias=phnB;}
115	void G4QInelastic::SetWeakNucBias(G4double ccnB) {weakNucBias=ccnB;}
116
117	// Destructor
118
119	G4QInelastic::~G4QInelastic()
120	{
121	ElementZ.clear();
122	ElProbInMat.clear();
123	for(unsigned i=0; i < ElIsoN.size(); ++i)
124	{
125	std::vector<G4int>* curEIN=ElIsoN[i];
126	curEIN->clear();
127	delete curEIN;
128	}
129	for(unsigned i=0; i < IsoProbInEl.size(); ++i)
130	{
131	std::vector<G4double>* curEIN=IsoProbInEl[i];
132	curEIN->clear();
133	delete curEIN;
134	}
135	}
136
137
138	G4LorentzVector G4QInelastic::GetEnegryMomentumConservation()
139	{
140	return EnMomConservation;
141	}
142
143	G4int G4QInelastic::GetNumberOfNeutronsInTarget()
144	{
145	return nOfNeutrons;
146	}
147
148	// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
149	// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
150	// ******** All CHIPS cross sections are calculated in the surface units **********
151	G4double G4QInelastic::GetMeanFreePath(const G4Track& aTrack,G4double,G4ForceCondition* Fc)
152	{
153	#ifdef debug
154	G4cout<<"G4QInelastic::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
155	#endif
156	*Fc = NotForced;
157	#ifdef debug
158	G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDynPart"<<G4endl;
159	#endif
160	const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
161	#ifdef debug
162	G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDef"<<G4endl;
163	#endif
164	G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
165	if( !IsApplicable(*incidentParticleDefinition))
166	G4cout<<"-W-G4QInelastic::GetMeanFreePath called for not implemented particle"<<G4endl;
167	// Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
168	G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
169	#ifdef debug
170	G4cout<<"G4QCollis::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
171	#endif
172	const G4Material* material = aTrack.GetMaterial(); // Get the current material
173	const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
174	const G4ElementVector* theElementVector = material->GetElementVector();
175	G4int nE=material->GetNumberOfElements();
176	#ifdef debug
177	G4cout<<"G4QInelastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
178	#endif
179	G4bool leptoNuc=false; // By default the reaction is not lepto-nuclear
180	G4VQCrossSection* CSmanager=0;
181	G4VQCrossSection* CSmanager2=0;
182	G4int pPDG=0;
183	G4int pZ = incidentParticleDefinition->GetAtomicNumber();
184	G4int pA = incidentParticleDefinition->GetAtomicMass();
185	if(incidentParticleDefinition == G4Neutron::Neutron()) // @@ make a switch
186	{
187	CSmanager=G4QNeutronNuclearCrossSection::GetPointer();
188	#ifdef debug
189	G4cout<<"G4QInelastic::GetMeanFreePath: CSmanager is defined for the neutron"<<G4endl;
190	#endif
191	pPDG=2112;
192	}
193	else if(incidentParticleDefinition == G4Proton::Proton())
194	{
195	CSmanager=G4QProtonNuclearCrossSection::GetPointer();
196	pPDG=2212;
197	}
198	else if(incidentParticleDefinition == G4PionMinus::PionMinus())
199	{
200	CSmanager=G4QPionMinusNuclearCrossSection::GetPointer();
201	pPDG=-211;
202	}
203	else if(incidentParticleDefinition == G4PionPlus::PionPlus())
204	{
205	CSmanager=G4QPionPlusNuclearCrossSection::GetPointer();
206	pPDG=211;
207	}
208	else if(incidentParticleDefinition == G4KaonMinus::KaonMinus())
209	{
210	CSmanager=G4QKaonMinusNuclearCrossSection::GetPointer();
211	pPDG=-321;
212	}
213	else if(incidentParticleDefinition == G4KaonPlus::KaonPlus())
214	{
215	CSmanager=G4QKaonPlusNuclearCrossSection::GetPointer();
216	pPDG=321;
217	}
218	else if(incidentParticleDefinition == G4KaonZeroLong::KaonZeroLong() \|\|
219	incidentParticleDefinition == G4KaonZeroShort::KaonZeroShort() \|\|
220	incidentParticleDefinition == G4KaonZero::KaonZero() \|\|
221	incidentParticleDefinition == G4AntiKaonZero::AntiKaonZero() )
222	{
223	CSmanager=G4QKaonZeroNuclearCrossSection::GetPointer();
224	if(G4UniformRand() > 0.5) pPDG= 311;
225	else pPDG=-311;
226	}
227	else if(incidentParticleDefinition == G4Lambda::Lambda())
228	{
229	CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
230	pPDG=3122;
231	}
232	else if(pZ > 0 && pA > 1) // Ions (not implemented yet (should not be used)
233	{
234	G4cout<<"-Warning-G4QInelastic::GetMeanFreePath: G4QInelastic called for ions"<<G4endl;
235	CSmanager=G4QProtonNuclearCrossSection::GetPointer();
236	pPDG=90000000+999*pZ+pA;
237	}
238	else if(incidentParticleDefinition == G4SigmaPlus::SigmaPlus())
239	{
240	CSmanager=G4QHyperonPlusNuclearCrossSection::GetPointer();
241	pPDG=3222;
242	}
243	else if(incidentParticleDefinition == G4SigmaMinus::SigmaMinus())
244	{
245	CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
246	pPDG=3112;
247	}
248	else if(incidentParticleDefinition == G4SigmaZero::SigmaZero())
249	{
250	CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
251	pPDG=3212;
252	}
253	else if(incidentParticleDefinition == G4XiMinus::XiMinus())
254	{
255	CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
256	pPDG=3312;
257	}
258	else if(incidentParticleDefinition == G4XiZero::XiZero())
259	{
260	CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
261	pPDG=3322;
262	}
263	else if(incidentParticleDefinition == G4OmegaMinus::OmegaMinus())
264	{
265	CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
266	pPDG=3334;
267	}
268	else if(incidentParticleDefinition == G4MuonPlus::MuonPlus() \|\|
269	incidentParticleDefinition == G4MuonMinus::MuonMinus())
270	{
271	CSmanager=G4QMuonNuclearCrossSection::GetPointer();
272	leptoNuc=true;
273	pPDG=13;
274	}
275	else if(incidentParticleDefinition == G4AntiNeutron::AntiNeutron())
276	{
277	CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
278	pPDG=-2112;
279	}
280	else if(incidentParticleDefinition == G4AntiProton::AntiProton())
281	{
282	CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
283	pPDG=-2212;
284	}
285	else if(incidentParticleDefinition == G4AntiLambda::AntiLambda())
286	{
287	CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
288	pPDG=-3122;
289	}
290	else if(incidentParticleDefinition == G4AntiSigmaPlus::AntiSigmaPlus())
291	{
292	CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
293	pPDG=-3222;
294	}
295	else if(incidentParticleDefinition == G4AntiSigmaMinus::AntiSigmaMinus())
296	{
297	CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
298	pPDG=-3112;
299	}
300	else if(incidentParticleDefinition == G4AntiSigmaZero::AntiSigmaZero())
301	{
302	CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
303	pPDG=-3212;
304	}
305	else if(incidentParticleDefinition == G4AntiXiMinus::AntiXiMinus())
306	{
307	CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
308	pPDG=-3312;
309	}
310	else if(incidentParticleDefinition == G4AntiXiZero::AntiXiZero())
311	{
312	CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
313	pPDG=-3322;
314	}
315	else if(incidentParticleDefinition == G4AntiOmegaMinus::AntiOmegaMinus())
316	{
317	CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
318	pPDG=-3334;
319	}
320	else if(incidentParticleDefinition == G4Gamma::Gamma())
321	{
322	CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
323	pPDG=22;
324	}
325	else if(incidentParticleDefinition == G4Electron::Electron() \|\|
326	incidentParticleDefinition == G4Positron::Positron())
327	{
328	CSmanager=G4QElectronNuclearCrossSection::GetPointer();
329	leptoNuc=true;
330	pPDG=11;
331	}
332	else if(incidentParticleDefinition == G4TauPlus::TauPlus() \|\|
333	incidentParticleDefinition == G4TauMinus::TauMinus())
334	{
335	CSmanager=G4QTauNuclearCrossSection::GetPointer();
336	leptoNuc=true;
337	pPDG=15;
338	}
339	else if(incidentParticleDefinition == G4NeutrinoMu::NeutrinoMu() )
340	{
341	CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
342	CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
343	leptoNuc=true;
344	pPDG=14;
345	}
346	else if(incidentParticleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMu() )
347	{
348	CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
349	CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
350	leptoNuc=true;
351	pPDG=-14;
352	}
353	else if(incidentParticleDefinition == G4NeutrinoE::NeutrinoE() )
354	{
355	CSmanager=G4QNuENuclearCrossSection::GetPointer();
356	CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
357	leptoNuc=true;
358	pPDG=12;
359	}
360	else if(incidentParticleDefinition == G4AntiNeutrinoE::AntiNeutrinoE() )
361	{
362	CSmanager=G4QANuENuclearCrossSection::GetPointer();
363	CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
364	leptoNuc=true;
365	pPDG=-12;
366	}
367	else G4cout<<"-Warning-G4QInelastic::GetMeanFreePath:Particle "
368	<<incidentParticleDefinition->GetPDGEncoding()
369	<<" isn't supported by CHIPS"<<G4endl;
370
371	G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
372	G4double sigma=0.; // Sums over elements for the material
373	G4int IPIE=IsoProbInEl.size(); // How many old elements?
374	if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
375	{
376	std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
377	SPI->clear();
378	delete SPI;
379	std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
380	IsN->clear();
381	delete IsN;
382	}
383	ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
384	ElementZ.clear(); // Clear the body vector for Z of Elements
385	IsoProbInEl.clear(); // Clear the body vector for SPI
386	ElIsoN.clear(); // Clear the body vector for N of Isotopes
387	for(G4int i=0; i<nE; ++i)
388	{
389	G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
390	G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
391	ElementZ.push_back(Z); // Remember Z of the Element
392	G4int isoSize=0; // The default for the isoVectorLength is 0
393	G4int indEl=0; // Index of non-trivial element or 0(default)
394	G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
395	if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
396	#ifdef debug
397	G4cout<<"G4QInelastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
398	#endif
399	if(isoSize) // The Element has non-trivial abundance set
400	{
401	indEl=pElement->GetIndex()+1; // Index of the non-trivial element
402	if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
403	{
404	std::vector<std::pair<G4int,G4double>> newAbund =
405	new std::vector<std::pair<G4int,G4double>*>;
406	G4double* abuVector=pElement->GetRelativeAbundanceVector();
407	for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
408	{
409	G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
410	if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QInelastic::GetMeanFreePath"
411	<<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
412	G4double abund=abuVector[j];
413	std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
414	#ifdef debug
415	G4cout<<"G4QInelastic::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
416	#endif
417	newAbund->push_back(pr);
418	}
419	#ifdef debug
420	G4cout<<"G4QInelastic::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
421	#endif
422	indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
423	for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
424	delete newAbund; // Was "new" in the beginning of the name space
425	}
426	}
427	std::vector<std::pair<G4int,G4double>> cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
428	std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
429	IsoProbInEl.push_back(SPI);
430	std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
431	ElIsoN.push_back(IsN);
432	G4int nIs=cs->size(); // A#Of Isotopes in the Element
433	G4double susi=0.; // sum of CS over isotopes
434	#ifdef debug
435	G4cout<<"G4QInelastic::GetMeanFreePath: Before Loop nIs="<<nIs<<G4endl;
436	#endif
437	if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
438	{
439	std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
440	G4int N=curIs->first; // #of Neuterons in the isotope j of El i
441	IsN->push_back(N); // Remember Min N for the Element
442	#ifdef debug
443	G4cout<<"G4QCollis::GetMeanFrP: Before CS, P="<<Momentum<<",Z="<<Z<<",N="<<N<<G4endl;
444	#endif
445	if(!pPDG) G4cout<<"-Warning-G4QCollis::GetMeanFrP: (1) projectile PDG=0"<<G4endl;
446	G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
447	if(CSmanager2)CSI+=CSmanager2->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i)nu,nu
448	#ifdef debug
449	G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
450	#endif
451	curIs->second = CSI;
452	susi+=CSI; // Make a sum per isotopes
453	SPI->push_back(susi); // Remember summed cross-section
454	} // End of temporary initialization of the cross sections in the G4QIsotope singeltone
455	sigma+=Isotopes->GetMeanCrossSection(Z,indEl)NOfNucPerVolume[i];//SUM(MeanCSNOfNperV)
456	ElProbInMat.push_back(sigma);
457	} // End of LOOP over Elements
458	#ifdef debug
459	G4cout<<"G4QCol::GetMeanFrPa: S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
460	#endif
461	// Check that cross section is not zero and return the mean free path
462	if(photNucBias!=1.) if(incidentParticleDefinition == G4Gamma::Gamma() \|\|
463	incidentParticleDefinition == G4MuonPlus::MuonPlus() \|\|
464	incidentParticleDefinition == G4MuonMinus::MuonMinus() \|\|
465	incidentParticleDefinition == G4Electron::Electron() \|\|
466	incidentParticleDefinition == G4Positron::Positron() \|\|
467	incidentParticleDefinition == G4TauMinus::TauMinus() \|\|
468	incidentParticleDefinition == G4TauPlus::TauPlus() )
469	sigma*=photNucBias;
470	if(weakNucBias!=1.) if(incidentParticleDefinition==G4NeutrinoE::NeutrinoE() \|\|
471	incidentParticleDefinition==G4AntiNeutrinoE::AntiNeutrinoE() \|\|
472	incidentParticleDefinition==G4NeutrinoTau::NeutrinoTau() \|\|
473	incidentParticleDefinition==G4AntiNeutrinoTau::AntiNeutrinoTau()\|\|
474	incidentParticleDefinition==G4NeutrinoMu::NeutrinoMu() \|\|
475	incidentParticleDefinition==G4AntiNeutrinoMu::AntiNeutrinoMu() )
476	sigma*=weakNucBias;
477	if(sigma > 0.) return 1./sigma; // Mean path [distance]
478	return DBL_MAX;
479	}
480
481
482	G4bool G4QInelastic::IsApplicable(const G4ParticleDefinition& particle)
483	{
484	G4int Z=particle.GetAtomicNumber();
485	G4int A=particle.GetAtomicMass();
486	if (particle == *( G4Neutron::Neutron() )) return true;
487	else if (particle == *( G4Proton::Proton() )) return true;
488	else if (particle == *( G4MuonPlus::MuonPlus() )) return true;
489	else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
490	else if (particle == *( G4Gamma::Gamma() )) return true;
491	else if (particle == *( G4Electron::Electron() )) return true;
492	else if (particle == *( G4Positron::Positron() )) return true;
493	else if (particle == *( G4PionMinus::PionMinus() )) return true;
494	else if (particle == *( G4PionPlus::PionPlus() )) return true;
495	else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
496	else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
497	else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
498	else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
499	else if (particle == *( G4Lambda::Lambda() )) return true;
500	else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
501	else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
502	else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
503	else if (particle == *( G4XiMinus::XiMinus() )) return true;
504	else if (particle == *( G4XiZero::XiZero() )) return true;
505	else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
506	else if (particle == *(G4GenericIon::GenericIon()) \|\| (Z > 0 && A > 1)) return true;
507	else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
508	else if (particle == *( G4AntiProton::AntiProton() )) return true;
509	else if (particle == *( G4AntiLambda::AntiLambda() )) return true;
510	else if (particle == *( G4AntiSigmaPlus::AntiSigmaPlus() )) return true;
511	else if (particle == *(G4AntiSigmaMinus::AntiSigmaMinus())) return true;
512	else if (particle == *( G4AntiSigmaZero::AntiSigmaZero() )) return true;
513	else if (particle == *( G4AntiXiMinus::AntiXiMinus() )) return true;
514	else if (particle == *( G4AntiXiZero::AntiXiZero() )) return true;
515	else if (particle == *(G4AntiOmegaMinus::AntiOmegaMinus())) return true;
516	else if (particle == *( G4TauPlus::TauPlus() )) return true;
517	else if (particle == *( G4TauMinus::TauMinus() )) return true;
518	else if (particle == *( G4AntiNeutrinoE::AntiNeutrinoE() )) return true;
519	else if (particle == *( G4NeutrinoE::NeutrinoE() )) return true;
520	else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
521	else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
522	//else if (particle == *(G4AntiNeutrinoTau::AntiNeutrinoTau())) return true;
523	//else if (particle == *( G4NeutrinoTau::NeutrinoTau() )) return true;
524	#ifdef debug
525	G4cout<<"***G4QInelastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
526	#endif
527	return false;
528	}
529
530	G4VParticleChange* G4QInelastic::PostStepDoIt(const G4Track& track, const G4Step& step)
531	{
532	static const G4double third = 1./3.;
533	static const G4double me=G4Electron::Electron()->GetPDGMass(); // electron mass
534	static const G4double me2=me*me; // squared electron mass
535	static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass
536	static const G4double mu2=mu*mu; // squared muon mass
537	static const G4double mt=G4TauMinus::TauMinus()->GetPDGMass(); // tau mass
538	static const G4double mt2=mt*mt; // squared tau mass
539	//static const G4double dpi=M_PI+M_PI; // 2pi (for Phi distr.) changed to twopi*
540	static const G4double mNeut= G4QPDGCode(2112).GetMass();
541	static const G4double mNeut2= mNeut*mNeut;
542	static const G4double mProt= G4QPDGCode(2212).GetMass();
543	static const G4double mProt2= mProt*mProt;
544	static const G4double dM=mProt+mNeut; // doubled nucleon mass
545	static const G4double hdM=dM/2.; // M of the "nucleon"
546	static const G4double hdM2=hdM*hdM; // M2 of the "nucleon"
547	static const G4double mPi0 = G4QPDGCode(111).GetMass();
548	static const G4double mPi0s= mPi0*mPi0;
549	static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
550	//static const G4double mDeut2=mDeut*mDeut; // SquaredMassOfDeuteron
551	static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
552	static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
553	static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
554	static const G4double mPi = G4QPDGCode(211).GetMass();
555	static const G4double tmPi = mPi+mPi; // Doubled mass of the charged pion
556	static const G4double stmPi= tmPi*tmPi; // Squared Doubled mass of the charged pion
557	static const G4double mPPi = mPi+mProt; // Delta threshold
558	static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2
559	//static const G4double mDel2= 1400*1400; // Delta up threshold for W2 (in MeV^2)
560	// Static definitions for electrons (nu,e) -----------------------------------------
561	static const G4double meN = mNeut+me;
562	static const G4double meN2= meN*meN;
563	static const G4double fmeN= 4mNeut2me2;
564	static const G4double mesN= mNeut2+me2;
565	static const G4double meP = mProt+me;
566	static const G4double meP2= meP*meP;
567	static const G4double fmeP= 4mProt2me2;
568	static const G4double mesP= mProt2+me2;
569	static const G4double medM= me2/dM; // for x limit
570	static const G4double meD = mPPi+me; // Multiperipheral threshold
571	static const G4double meD2= meD*meD;
572	// Static definitions for muons (nu,mu) -----------------------------------------
573	static const G4double muN = mNeut+mu;
574	static const G4double muN2= muN*muN;
575	static const G4double fmuN= 4mNeut2mu2;
576	static const G4double musN= mNeut2+mu2;
577	static const G4double muP = mProt+mu;
578	static const G4double muP2= muP*muP; // +
579	static const G4double fmuP= 4mProt2mu2; // +
580	static const G4double musP= mProt2+mu2;
581	static const G4double mudM= mu2/dM; // for x limit
582	static const G4double muD = mPPi+mu; // Multiperipheral threshold
583	static const G4double muD2= muD*muD;
584	// Static definitions for muons (nu,nu) -----------------------------------------
585	//static const G4double nuN = mNeut;
586	//static const G4double nuN2= mNeut2;
587	//static const G4double fnuN= 0.;
588	//static const G4double nusN= mNeut2;
589	//static const G4double nuP = mProt;
590	//static const G4double nuP2= mProt2;
591	//static const G4double fnuP= 0.;
592	//static const G4double nusP= mProt2;
593	//static const G4double nudM= 0.; // for x limit
594	//static const G4double nuD = mPPi; // Multiperipheral threshold
595	//static const G4double nuD2= mPPi2;
596	//-------------------------------------------------------------------------------------
597	static G4bool CWinit = true; // CHIPS Warld needs to be initted
598	if(CWinit)
599	{
600	CWinit=false;
601	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
602	}
603	//-------------------------------------------------------------------------------------
604	const G4DynamicParticle* projHadron = track.GetDynamicParticle();
605	const G4ParticleDefinition* particle=projHadron->GetDefinition();
606	#ifdef debug
607	G4cout<<"G4QInelastic::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
608	#endif
609	G4ForceCondition cond=NotForced;
610	GetMeanFreePath(track, 1., &cond);
611	#ifdef debug
612	G4cout<<"G4QInelastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
613	#endif
614	G4bool scat=false; // No CHEX in proj scattering
615	G4int scatPDG=0; // Must be filled if true (CHEX)
616	G4LorentzVector proj4M=projHadron->Get4Momentum(); // 4-momentum of the projectile (IU?)
617	G4LorentzVector scat4M=proj4M; // Must be filled if true
618	G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
619	G4double Momentum=proj4M.rho();
620	if(std::fabs(Momentum-momentum)>.001)
621	G4cerr<<"*G4QInelastic::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl;
622	#ifdef debug
623	G4double mp=proj4M.m();
624	G4cout<<"-->G4QCollis::PostStDoIt:called, 4M="<<proj4M<<", P="<<Momentum<<"="<<momentum
625	<<",m="<<mp<<G4endl;
626	#endif
627	if (!IsApplicable(*particle)) // Check applicability
628	{
629	G4cerr<<"G4QInelastic::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
630	<<G4endl;
631	return 0;
632	}
633	const G4Material* material = track.GetMaterial(); // Get the current material
634	G4int Z=0;
635	const G4ElementVector* theElementVector = material->GetElementVector();
636	G4int nE=material->GetNumberOfElements();
637	#ifdef debug
638	G4cout<<"G4QInelastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
639	#endif
640	G4int projPDG=0; // PDG Code prototype for the captured hadron
641	G4int pZ=particle->GetAtomicNumber();
642	G4int pA=particle->GetAtomicMass();
643	if (particle == G4Neutron::Neutron() ) projPDG= 2112;
644	else if (particle == G4Proton::Proton() ) projPDG= 2212;
645	else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
646	else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
647	else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321;
648	else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
649	else if (particle == G4KaonZeroLong::KaonZeroLong() \|\|
650	particle == G4KaonZeroShort::KaonZeroShort() \|\|
651	particle == G4KaonZero::KaonZero() \|\|
652	particle == G4AntiKaonZero::AntiKaonZero() )
653	{
654	if(G4UniformRand() > 0.5) projPDG= 311;
655	else projPDG= -311;
656	}
657	else if ( pZ > 0 && pA > 1 ) projPDG = 90000000+999*pZ+pA;
658	else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
659	else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
660	else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
661	else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
662	else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
663	else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
664	else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
665	else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
666	else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
667	else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
668	else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
669	else if (particle == G4Gamma::Gamma() ) projPDG= 22;
670	else if (particle == G4Electron::Electron() ) projPDG= 11;
671	else if (particle == G4Positron::Positron() ) projPDG= -11;
672	else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14;
673	else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14;
674	else if (particle == G4AntiLambda::AntiLambda() ) projPDG=-3122;
675	else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) projPDG=-3222;
676	else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112;
677	else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) projPDG=-3212;
678	else if (particle == G4AntiXiMinus::AntiXiMinus() ) projPDG=-3312;
679	else if (particle == G4AntiXiZero::AntiXiZero() ) projPDG=-3322;
680	else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334;
681	else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12;
682	else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12;
683	else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
684	else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
685	else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
686	else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
687	G4int aProjPDG=std::abs(projPDG);
688	#ifdef debug
689	G4int prPDG=particle->GetPDGEncoding();
690	G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
691	#endif
692	if(!projPDG)
693	{
694	G4cerr<<"---Warning---G4QInelastic::PostStepDoIt:Undefined interacting hadron"<<G4endl;
695	return 0;
696	}
697	G4int EPIM=ElProbInMat.size();
698	#ifdef debug
699	G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
700	#endif
701	G4int i=0;
702	if(EPIM>1)
703	{
704	G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
705	for(i=0; i<nE; ++i)
706	{
707	#ifdef debug
708	G4cout<<"G4QInelastic::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
709	#endif
710	if (rnd<ElProbInMat[i]) break;
711	}
712	if(i>=nE) i=nE-1; // Top limit for the Element
713	}
714	G4Element* pElement=(*theElementVector)[i];
715	Z=static_cast<G4int>(pElement->GetZ());
716	#ifdef debug
717	G4cout<<"G4QInelastic::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
718	#endif
719	if(Z<=0)
720	{
721	G4cerr<<"---Warning---G4QInelastic::PostStepDoIt: Element with Z="<<Z<<G4endl;
722	if(Z<0) return 0;
723	}
724	std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
725	std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
726	G4int nofIsot=SPI->size(); // #of isotopes in the element i
727	#ifdef debug
728	G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
729	#endif
730	G4int j=0;
731	if(nofIsot>1)
732	{
733	G4double rndI=(SPI)[nofIsot-1]G4UniformRand(); // Randomize the isotop of the Element
734	for(j=0; j<nofIsot; ++j)
735	{
736	#ifdef debug
737	G4cout<<"G4QInelastic::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
738	#endif
739	if(rndI < (*SPI)[j]) break;
740	}
741	if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
742	}
743	G4int N =(*IsN)[j]; ; // Randomized number of neutrons
744	#ifdef debug
745	G4cout<<"G4QInelastic::PostStepDoIt: Z="<<Z<<", j="<<i<<", N(isotope)="<<N<<G4endl;
746	#endif
747	if(N<0)
748	{
749	G4cerr<<"-Warning-G4QInelastic::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
750	return 0;
751	}
752	nOfNeutrons=N; // Remember it for the energy-momentum check
753	G4double dd=0.025;
754	G4double am=Z+N;
755	G4double sr=std::sqrt(am);
756	G4double dsr=0.01*(sr+sr);
757	if(dsr<dd)dsr=dd;
758	if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // ManualP
759	else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP ClustPars
760	else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.);//Pi- ClustPars
761	#ifdef debug
762	G4cout<<"G4QInelastic::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
763	#endif
764	if(N<0)
765	{
766	G4cerr<<"---Warning---G4QInelastic::PostStepDoIt:Element with N="<<N<< G4endl;
767	return 0;
768	}
769	aParticleChange.Initialize(track);
770	G4double weight = track.GetWeight();
771	#ifdef debug
772	G4cout<<"G4QInelastic::PostStepDoIt: weight="<<weight<<G4endl;
773	#endif
774	if(photNucBias!=1.) weight/=photNucBias;
775	else if(weakNucBias!=1.) weight/=weakNucBias;
776	G4double localtime = track.GetGlobalTime();
777	#ifdef debug
778	G4cout<<"G4QInelastic::PostStepDoIt: localtime="<<localtime<<G4endl;
779	#endif
780	G4ThreeVector position = track.GetPosition();
781	G4TouchableHandle trTouchable = track.GetTouchableHandle();
782	#ifdef debug
783	G4cout<<"G4QInelastic::PostStepDoIt: position="<<position<<G4endl;
784	#endif
785	//
786	G4int targPDG=90000000+Z*1000+N; // PDG Code of the target nucleus
787	G4QPDGCode targQPDG(targPDG);
788	G4double tgM=targQPDG.GetMass(); // Target mass
789	G4double tM=tgM; // Target mass (copy to be changed)
790	G4QHadronVector* output=new G4QHadronVector;// Prototype of Fragm Output G4QHadronVector
791	G4double absMom = 0.; // Prototype of absorbed by nucleus Moment
792	G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVector
793	G4LorentzVector lead4M(0.,0.,0.,0.); // Prototype of LeadingQ 4-momentum
794	#ifdef debug
795	G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<aProjPDG<<", targPDG="<<targPDG<<G4endl;
796	#endif
797	//
798	// Leptons with photonuclear
799	// Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?)
800	//
801	if (aProjPDG == 11 \|\| aProjPDG == 13 \|\| aProjPDG == 15)
802	{
803	#ifdef debug
804	G4cout<<"G4QInelastic::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl;
805	#endif
806	G4double kinEnergy= projHadron->GetKineticEnergy();
807	G4ParticleMomentum dir = projHadron->GetMomentumDirection();
808	G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
809	G4double ml=me;
810	G4double ml2=me2;
811	if(aProjPDG== 13)
812	{
813	CSmanager=G4QMuonNuclearCrossSection::GetPointer();
814	ml=mu;
815	ml2=mu2;
816	}
817	if(aProjPDG== 15)
818	{
819	CSmanager=G4QTauNuclearCrossSection::GetPointer();
820	ml=mt;
821	ml2=mt2;
822	}
823	// @@ Probably this is not necessary any more (?)
824	if(!aProjPDG) G4cout<<"-Warning-G4QCollis::PostStepDoIt: (2) projectile PDG=0"<<G4endl;
825	G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,aProjPDG);// Recalculate XS
826	// @@ check a possibility to separate p, n, or alpha (!)
827	if(xSec <= 0.) // The cross-section is 0 -> Do Nothing
828	{
829	#ifdef debug
830	G4cerr<<"---OUT---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl;
831	#endif
832	//Do Nothing Action insead of the reaction
833	aParticleChange.ProposeEnergy(kinEnergy);
834	aParticleChange.ProposeLocalEnergyDeposit(0.);
835	aParticleChange.ProposeMomentumDirection(dir);
836	aParticleChange.ProposeTrackStatus(fAlive);
837	delete output;
838	return G4VDiscreteProcess::PostStepDoIt(track,step);
839	}
840	G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
841	#ifdef debug
842	G4cout<<"G4QCol::PStDoIt: kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
843	#endif
844	if( kinEnergy < photonEnergy \|\| photonEnergy < 0.)
845	{
846	//Do Nothing Action insead of the reaction
847	G4cerr<<"--G4QInelastic::PSDoIt: photE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl;
848	aParticleChange.ProposeEnergy(kinEnergy);
849	aParticleChange.ProposeLocalEnergyDeposit(0.);
850	aParticleChange.ProposeMomentumDirection(dir);
851	aParticleChange.ProposeTrackStatus(fAlive);
852	delete output;
853	return G4VDiscreteProcess::PostStepDoIt(track,step);
854	}
855	G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart
856	G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon
857	if(tM<999.) W-=mPi0+mPi0s/dM; // Pion production threshold for a nucleon target
858	if(W<0.)
859	{
860	//Do Nothing Action insead of the reaction
861	#ifdef debug
862	G4cout<<"--G4QInelastic::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl;
863	#endif
864	aParticleChange.ProposeEnergy(kinEnergy);
865	aParticleChange.ProposeLocalEnergyDeposit(0.);
866	aParticleChange.ProposeMomentumDirection(dir);
867	aParticleChange.ProposeTrackStatus(fAlive);
868	delete output;
869	return G4VDiscreteProcess::PostStepDoIt(track,step);
870	}
871	// Update G4VParticleChange for the scattered muon
872	G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer();
873	G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy,Z,N,22);//Integrated XS
874	G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N, 22); // Real XS
875	G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2);
876	if(sigNuG4UniformRand()>sigKrndFraction)
877	{
878	//Do NothingToDo Action insead of the reaction
879	#ifdef debug
880	G4cout<<"-DoNoth-G4QInelastic::PostStepDoIt: probab. correction - DoNothing"<<G4endl;
881	#endif
882	aParticleChange.ProposeEnergy(kinEnergy);
883	aParticleChange.ProposeLocalEnergyDeposit(0.);
884	aParticleChange.ProposeMomentumDirection(dir);
885	aParticleChange.ProposeTrackStatus(fAlive);
886	delete output;
887	return G4VDiscreteProcess::PostStepDoIt(track,step);
888	}
889	G4double iniE=kinEnergy+ml; // Initial total energy of the lepton
890	G4double finE=iniE-photonEnergy; // Final total energy of the lepton
891	#ifdef pdebug
892	G4cout<<"G4QInelastic::PoStDoIt:E="<<iniE<<",lE="<<finE<<"-"<<ml<<"="<<finE-ml<<G4endl;
893	#endif
894	aParticleChange.ProposeEnergy(finE-ml);
895	if(finE<=ml) // Secondary lepton (e/mu/tau) at rest disappears
896	{
897	aParticleChange.ProposeEnergy(0.);
898	if(aProjPDG== 11) aParticleChange.ProposeTrackStatus(fStopAndKill);
899	else aParticleChange.ProposeTrackStatus(fStopButAlive);
900	aParticleChange.ProposeMomentumDirection(dir);
901	}
902	else aParticleChange.ProposeTrackStatus(fAlive);
903	G4double iniP=std::sqrt(iniE*iniE-ml2); // Initial momentum of the electron
904	G4double finP=std::sqrt(finE*finE-ml2); // Final momentum of the electron
905	G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton)
906	#ifdef pdebug
907	G4cout<<"G4QC::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
908	#endif
909	if(cost>1.) cost=1.; // To avoid the accuracy of calculation problem
910	if(cost<-1.) cost=-1.; // To avoid the accuracy of calculation problem
911	//
912	// Scatter the lepton ( @@ make the same thing for real photons)
913	// At this point we have photonEnergy and photonQ2 (with notDefinedPhi)->SelectProjPart
914	G4double absEn = std::pow(am,third)*GeV; // @@(b) Mean Energy Absorbed by a Nucleus
915	if(am>1 && absEn < photonEnergy) // --> the absorption of energy can happen
916	//if(absEn < photonEnergy) // --> the absorption of energy can happen
917	{
918	G4double abtEn = absEn+hdM; // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N)
919	G4double abEn2 = abtEn*abtEn; // Squared absorbed Energy + MN
920	G4double abMo2 = abEn2-hdM2; // Squared absorbed Momentum of compound system
921	G4double phEn2 = photonEnergy*photonEnergy;
922	G4double phMo2 = phEn2+photonQ2; // Squared momentum of primary virtual photon
923	G4double phMo = std::sqrt(phMo2); // Momentum of the primary virtual photon
924	absMom = std::sqrt(abMo2); // Absorbed Momentum
925	if(absMom < phMo) // --> the absorption of momentum can happen
926	{
927	G4double dEn = photonEnergy - absEn; // Leading energy
928	G4double dMo = phMo - absMom; // Leading momentum
929	G4double sF = dEndEn - dModMo;// s of leading particle
930	#ifdef ppdebug
931	G4cout<<"-PhotoAbsorption-G4QCol::PStDoIt:sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
932	#endif
933	if(sF > stmPi) // --> Leading fragmentation is possible
934	{
935	photonEnergy = absEn; // New value of the photon energy
936	photonQ2=abMo2-absEn*absEn; // New value of the photon Q2
937	absEn = dEn; // Put energy of leading particle to absEn (!)
938	}
939	else absMom=0.; // Flag that nothing has happened
940	}
941	else absMom=0.; // Flag that nothing has happened
942	}
943	// ------------- End of ProjPart selection
944	//
945	// Scattering in respect to the derection of the incident muon is made impicitly:
946	G4ThreeVector ort=dir.orthogonal(); // Not normed orthogonal vector (!) (to dir)
947	G4ThreeVector ortx = ort.unit(); // First unit vector orthogonal to the direction
948	G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction
949	G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
950	G4double phi=twopi*G4UniformRand(); // phi of scattered electron
951	G4double sinx=sint*std::sin(phi); // x perpendicular component
952	G4double siny=sint*std::cos(phi); // y perpendicular component
953	G4ThreeVector findir=costdir+sinxortx+siny*orty;
954	aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton
955	#ifdef pdebug
956	G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
957	<<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl;
958	#endif
959	G4ThreeVector photon3M=iniPdir-finPfindir;// 3D total momentum of photon
960	if(absMom) // Photon must be reduced & LeadingSyst fragmented
961	{
962	G4double ptm=photon3M.mag(); // 3M of the virtual photon
963	#ifdef ppdebug
964	G4cout<<"-Absorption-G4QInelastic::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
965	<<iniPdir<<", eFin3M="<<finPfindir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl;
966	#endif
967	G4ThreeVector lead3M=photon3M*(ptm-absMom)/ptm; // Keep the direction for leading Q
968	photon3M-=lead3M; // Reduced photon Momentum (photEn already = absEn)
969	proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System
970	#ifdef ppdebug
971	G4cout<<"-->G4QC::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
972	#endif
973	lead4M=proj4M; // Remember 4-mom for the total 4-momentum
974	G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+
975	try // \|
976	{ // \|
977	if(leadhs) delete leadhs; // \|
978	G4QNucleus vac(90000000); // \|
979	leadhs=pan->Fragment(vac,1); // DELETED after it is copied to output vector \|
980	} // \|
981	catch (G4QException& error) // \|
982	{ // \|
983	G4cerr<<"***G4QInelastic::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//\|
984	G4Exception("G4QInelastic::PostStepDoIt:","72",FatalException,"QuasmonCrash"); //\|
985	} // \|
986	delete pan; // Delete the Nuclear Environment <----<---+
987	#ifdef ppdebug
988	G4cout<<"G4QCol::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
989	<<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl;
990	#endif
991	}
992	else
993	{
994	G4int qNH=0;
995	if(leadhs) qNH=leadhs->size();
996	if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
997	if(leadhs) delete leadhs;
998	leadhs=0;
999	}
1000	projPDG=22;
1001	proj4M=G4LorentzVector(photon3M,photonEnergy);
1002	#ifdef debug
1003	G4cout<<"G4QInelastic::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
1004	<<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl;
1005	#endif
1006
1007	//
1008	// neutrinoNuclear interactions (nu_e, nu_mu only)
1009	//
1010	}
1011	else if (aProjPDG == 12 \|\| aProjPDG == 14)
1012	{
1013	G4double kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino
1014	G4double dKinE=kinEnergy+kinEnergy; // doubled energy for s calculation
1015	#ifdef debug
1016	G4cout<<"G4QInelastic::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
1017	#endif
1018	G4ParticleMomentum dir = projHadron->GetMomentumDirection(); // unit vector
1019	G4double ml = mu;
1020	G4double ml2 = mu2;
1021	//G4double mlN = muN;
1022	G4double mlN2= muN2;
1023	G4double fmlN= fmuN;
1024	G4double mlsN= musN;
1025	//G4double mlP = muP;
1026	G4double mlP2= muP2;
1027	G4double fmlP= fmuP;
1028	G4double mlsP= musP;
1029	G4double mldM= mudM;
1030	//G4double mlD = muD;
1031	G4double mlD2= muD2;
1032	if(aProjPDG==12)
1033	{
1034	ml = me;
1035	ml2 = me2;
1036	//mlN = meN;
1037	mlN2= meN2;
1038	fmlN= fmeN;
1039	mlsN= mesN;
1040	//mlP = meP;
1041	mlP2= meP2;
1042	fmlP= fmeP;
1043	mlsP= mesP;
1044	mldM= medM;
1045	//mlD = meD;
1046	mlD2= meD2;
1047	}
1048	G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l)
1049	G4VQCrossSection* CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); // (nu,nu)
1050	proj4M=G4LorentzVector(dir*kinEnergy,kinEnergy); // temporary
1051	G4bool nuanu=true;
1052	scatPDG=13; // Prototype = secondary scattered mu-
1053	if(projPDG==-14)
1054	{
1055	nuanu=false; // Anti-neutrino
1056	CSmanager=G4QANuMuNuclearCrossSection::GetPointer(); // (anu,mu+) CC @@ open
1057	CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
1058	scatPDG=-13; // secondary scattered mu+
1059	}
1060	else if(projPDG==12)
1061	{
1062	CSmanager=G4QNuENuclearCrossSection::GetPointer(); // @@ open (only CC is changed)
1063	scatPDG=11; // secondary scattered e-
1064	}
1065	else if(projPDG==-12)
1066	{
1067	nuanu=false; // anti-neutrino
1068	CSmanager=G4QANuENuclearCrossSection::GetPointer(); // (anu,e+) CC @@ open
1069	CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
1070	scatPDG=-11; // secondary scattered e+
1071	}
1072	// @@ Probably this is not necessary any more
1073	if(!projPDG) G4cout<<"-Warning-G4QCollis::PostStepDoIt: (3) projectile PDG=0"<<G4endl;
1074	G4double xSec1=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG); //Recalculate XS
1075	G4double xSec2=CSmanager2->GetCrossSection(false,Momentum,Z,N,projPDG);//Recalculate XS
1076	G4double xSec=xSec1+xSec2;
1077	// @@ check a possibility to separate p, n, or alpha (!)
1078	if(xSec <= 0.) // The cross-section = 0 -> Do Nothing
1079	{
1080	G4cerr<<"G4QInelastic::PSDoIt:nuE="<<kinEnergy<<",X1="<<xSec1<<",X2="<<xSec2<<G4endl;
1081	//Do Nothing Action insead of the reaction
1082	aParticleChange.ProposeEnergy(kinEnergy);
1083	aParticleChange.ProposeLocalEnergyDeposit(0.);
1084	aParticleChange.ProposeMomentumDirection(dir);
1085	aParticleChange.ProposeTrackStatus(fAlive);
1086	delete output;
1087	return G4VDiscreteProcess::PostStepDoIt(track,step);
1088	}
1089	G4bool secnu=false;
1090	if(xSec*G4UniformRand()>xSec1) // recover neutrino/antineutrino
1091	{
1092	if(scatPDG>0) scatPDG++;
1093	else scatPDG--;
1094	secnu=true;
1095	}
1096	scat=true; // event with changed scattered projectile
1097	G4double totCS1 = CSmanager->GetLastTOTCS(); // the last total cross section1(isotope?)
1098	G4double totCS2 = CSmanager2->GetLastTOTCS();// the last total cross section2(isotope?)
1099	G4double totCS = totCS1+totCS2; // the last total cross section (isotope?)
1100	if(std::fabs(xSec-totCS*millibarn)/xSec>.0001)
1101	G4cout<<"-Warning-G4QInelastic::PostStepDoIt: xS="<<xSec<<"# CS="<<totCS<<G4endl;
1102	G4double qelCS1 = CSmanager->GetLastQELCS(); // the last quasi-elastic cross section1
1103	G4double qelCS2 = CSmanager2->GetLastQELCS();// the last quasi-elastic cross section2
1104	G4double qelCS = qelCS1+qelCS2; // the last quasi-elastic cross section
1105	if(totCS - qelCS < 0.) // only at low energies
1106	{
1107	totCS = qelCS;
1108	totCS1 = qelCS1;
1109	totCS2 = qelCS2;
1110	}
1111	// make different definitions for neutrino and antineutrino
1112	G4double mIN=mProt; // Just a prototype (for anu, Z=1, N=0)
1113	G4double mOT=mNeut;
1114	G4double OT=mlN2;
1115	G4double mOT2=mNeut2;
1116	G4double mlOT=fmlN;
1117	G4double mlsOT=mlsN;
1118	if(secnu)
1119	{
1120	if(am*G4UniformRand()>Z) // Neutron target
1121	{
1122	targPDG-=1; // subtract neutron
1123	projPDG=2112; // neutron is going out
1124	mIN =mNeut;
1125	OT =mNeut2;
1126	mOT2=mNeut2;
1127	mlOT=0.;
1128	mlsOT=mNeut2;
1129	}
1130	else
1131	{
1132	targPDG-=1000; // subtract neutron
1133	projPDG=2212; // neutron is going out
1134	mOT =mProt;
1135	OT =mProt2;
1136	mOT2 =mProt2;
1137	mlOT =0.;
1138	mlsOT=mProt2;
1139	}
1140	ml=0.;
1141	ml2=0.;
1142	mldM=0.;
1143	mlD2=mPPi2;
1144	G4QPDGCode targQPDG(targPDG);
1145	G4double rM=targQPDG.GetMass();
1146	mIN=tM-rM; // bounded in-mass of the neutron
1147	tM=rM;
1148	}
1149	else if(nuanu)
1150	{
1151	targPDG-=1; // Neutrino -> subtract neutron
1152	G4QPDGCode targQPDG(targPDG);
1153	G4double rM=targQPDG.GetMass();
1154	mIN=tM-rM; // bounded in-mass of the neutron
1155	tM=rM;
1156	mOT=mProt;
1157	OT=mlP2;
1158	mOT2=mProt2;
1159	mlOT=fmlP;
1160	mlsOT=mlsP;
1161	projPDG=2212; // proton is going out
1162	}
1163	else
1164	{
1165	if(Z>1\|\|N>0) // Calculate the splitted mass
1166	{
1167	targPDG-=1000; // Anti-Neutrino -> subtract proton
1168	G4QPDGCode targQPDG(targPDG);
1169	G4double rM=targQPDG.GetMass();
1170	mIN=tM-rM; // bounded in-mass of the proton
1171	tM=rM;
1172	}
1173	else
1174	{
1175	targPDG=0;
1176	mIN=tM;
1177	tM=0.;
1178	}
1179	projPDG=2112; // neutron is going out
1180	}
1181	G4double s=mIN*(mIN+dKinE); // s=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu)
1182	#ifdef debug
1183	G4cout<<"G4QInelastic::PostStDoIt: s="<<s<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
1184	#endif
1185	if(s<=OT) // * Do nothing solution *
1186	{
1187	//Do NothingToDo Action insead of the reaction (@@ Can we make it common?)
1188	G4cout<<"G4QInelastic::PostStepDoIt: tooSmallFinalMassOfCompound: DoNothing"<<G4endl;
1189	aParticleChange.ProposeEnergy(kinEnergy);
1190	aParticleChange.ProposeLocalEnergyDeposit(0.);
1191	aParticleChange.ProposeMomentumDirection(dir);
1192	aParticleChange.ProposeTrackStatus(fAlive);
1193	delete output;
1194	return G4VDiscreteProcess::PostStepDoIt(track,step);
1195	}
1196	#ifdef debug
1197	G4cout<<"G4QInelastic::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
1198	#endif
1199	aParticleChange.ProposeEnergy(0.);
1200	aParticleChange.ProposeTrackStatus(fStopAndKill); // the initial neutrino is killed
1201	// There is no way back from here !
1202	if ( ((secnu \|\| !nuanu \|\| N) && totCS*G4UniformRand() < qelCS) \|\| s < mlD2 )
1203	{ // Quasi-Elastic interaction
1204	G4double Q2=0.; // Simulate transferred momentum, in MeV^2
1205	if(secnu) Q2=CSmanager2->GetQEL_ExchangeQ2();
1206	else Q2=CSmanager->GetQEL_ExchangeQ2();
1207	#ifdef debug
1208	G4cout<<"G4QInelastic::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s<<",Q2="<<Q2<<G4endl;
1209	#endif
1210	//G4double ds=s+s; // doubled s
1211	G4double sqs=std::sqrt(s); // M_cm
1212	G4double dsqs=sqs+sqs; // 2*M_cm
1213	G4double pi=(s-mIN*mIN)/dsqs; // initial momentum in CMS (checked MK)
1214	G4double dpi=pi+pi; // doubled initial momentum in CMS
1215	G4double sd=s-mlsOT; // s-ml2-mOT2 (mlsOT=m^2_neut+m^2_lept)
1216	G4double qo2=(sd*sd-mlOT)/dsqs; // squared momentum of secondaries in CMS
1217	G4double qo=std::sqrt(qo2); // momentum of secondaries in CMS
1218	G4double cost=(dpi*std::sqrt(qo2+ml2)-Q2-ml2)/dpi/qo; // cos(theta) in CMS (chck MK)
1219	G4LorentzVector t4M(0.,0.,0.,mIN); // 4mom of the effective target
1220	G4LorentzVector c4M=t4M+proj4M; // 4mom of the compound system
1221	t4M.setT(mOT); // now it is 4mom of the outgoing nucleon
1222	scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scattered muon
1223	if(!G4QHadron(c4M).RelDecayIn2(scat4M, t4M, proj4M, cost, cost))
1224	{
1225	G4cerr<<"G4QCol::PSD:c4M="<<c4M<<sqs<<",mM="<<ml<<",tM="<<mOT<<",c="<<cost<<G4endl;
1226	throw G4QException("G4QInelastic::HadronizeQuasm: Can't dec QE nu,lept Compound");
1227	}
1228	proj4M=t4M; // 4mom of the new projectile nucleon
1229	}
1230	else // ***** Non Quasi Elastic interaction
1231	{
1232	// Recover the target PDG Code if the quasi-elastic scattering did not happened
1233	if ( (secnu && projPDG == 2212) \|\| (!secnu && projPDG == 2112) ) targPDG+=1;
1234	else if ( (secnu && projPDG == 2112) \|\| (!secnu && projPDG == 2212) ) targPDG+=1000;
1235	G4double Q2=0; // Simulate transferred momentum, in MeV^2
1236	if(secnu) Q2=CSmanager->GetNQE_ExchangeQ2();
1237	else Q2=CSmanager2->GetNQE_ExchangeQ2();
1238	#ifdef debug
1239	G4cout<<"G4QColl::PStDoIt: MultiPeriferal s="<<s<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
1240	#endif
1241	if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma
1242	else projPDG=CSmanager->GetExchangePDGCode(); // PDG Code of the effective pion
1243	//@@ Temporary made only for direct interaction and for N=3 (good for small Q2)
1244	//@@ inFuture use N=GetNPartons and directFraction=GetDirectPart, @@ W2...
1245	G4double r=G4UniformRand();
1246	G4double r1=0.5; // (1-x)
1247	if(r<0.5) r1=std::sqrt(r+r)(.5+.1579(r-.5));
1248	else if(r>0.5) r1=1.-std::sqrt(2.-r-r)(.5+.1579(.5-r));
1249	G4double xn=1.-mldM/Momentum; // Normalization of (1-x) [x>mldM/Mom]
1250	G4double x1=xn*r1; // (1-x)
1251	G4double x=1.-x1; // x=2k/M
1252	//G4double W2=(hdM2+Q2/x)*x1; // W2 candidate
1253	G4double mx=hdM*x; // Part of the target to interact with
1254	G4double we=Q2/(mx+mx); // transfered energy
1255	if(we>=kinEnergy-ml-.001) we=kinEnergy-ml-.0001; // safety to avoid nan=sqrt(neg)
1256	G4double pot=kinEnergy-we; // energy of the secondary lepton
1257	G4double mlQ2=ml2+Q2;
1258	G4double cost=(pot-mlQ2/dKinE)/std::sqrt(pot*pot-ml2); // LS cos(theta)
1259	if(std::fabs(cost)>1)
1260	{
1261	#ifdef debug
1262	G4cout<<"*G4QInelastic::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
1263	<<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl;
1264	#endif
1265	if(cost>1.) cost=1.;
1266	else cost=-1.;
1267	pot=mlQ2/dKinE+dKinE*ml2/mlQ2; // extreme output momentum
1268	}
1269	G4double lEn=std::sqrt(pot*pot+ml2); // Lepton energy
1270	G4double lPl=pot*cost; // Lepton longitudinal momentum
1271	G4double lPt=potstd::sqrt(1.-costcost); // Lepton transverse momentum
1272	std::pair<G4double,G4double> d2d=Random2DDirection(); // Randomize phi
1273	G4double lPx=lPt*d2d.first;
1274	G4double lPy=lPt*d2d.second;
1275	G4ThreeVector vdir=proj4M.vect(); // 3D momentum of the projectile
1276	G4ThreeVector vz= vdir.unit(); // Ort in the direction of the projectile
1277	G4ThreeVector vv= vz.orthogonal(); // Not normed orthogonal vector (!)
1278	G4ThreeVector vx= vv.unit(); // First ort orthogonal to the direction
1279	G4ThreeVector vy= vz.cross(vx); // Second ort orthoganal to the direction
1280	G4ThreeVector lP= lPlvz+lPxvx+lPy*vy; // 3D momentum of the scattered lepton
1281	scat4M=G4LorentzVector(lP,lEn); // 4mom of the scattered lepton
1282	proj4M-=scat4M; // 4mom of the W/Z (effective pion/gamma)
1283	#ifdef debug
1284	G4cout<<"G4QInelastic::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
1285	#endif
1286	// Check that the en/mom transfer is possible, if not -> elastic
1287	G4int fintPDG=targPDG; // Prototype for the compound nucleus
1288	if(!secnu)
1289	{
1290	if(projPDG<0) fintPDG-= 999;
1291	else fintPDG+= 999;
1292	}
1293	G4double fM=G4QPDGCode(fintPDG).GetMass();// compound nucleus Mass (MeV)
1294	G4double fM2=fM*fM;
1295	G4LorentzVector tg4M=G4LorentzVector(0.,0.,0.,tgM);
1296	G4LorentzVector c4M=tg4M+proj4M;
1297	#ifdef debug
1298	G4cout<<"G4QCol::PSDI:fM2="<<fM2<<" <? mc4M="<<c4M.m2()<<",dM="<<fM-tgM<<G4endl;
1299	#endif
1300	if(fM2>=c4M.m2()) // Elastic scattering should be done
1301	{
1302	G4LorentzVector tot4M=tg4M+proj4M+scat4M; // recover the total 4-momentum
1303	s=tot4M.m2();
1304	G4double fs=s-fM2-ml2;
1305	G4double fMl=fM2*ml2;
1306	G4double hQ2max=(fs*fs/2-fMl-fMl)/s; // Maximum possible Q2/2
1307	G4double cost=1.-Q2/hQ2max; // cos(theta) in CMS (use MultProd Q2)
1308	#ifdef debug
1309	G4cout<<"G4QC::PSDI:ct="<<cost<<",Q2="<<Q2<<",hQ2="<<hQ2max<<",4M="<<tot4M<<G4endl;
1310	#endif
1311	G4double acost=std::fabs(cost);
1312	if(acost>1.)
1313	{
1314	if(acost>1.001) G4cout<<"-Warning-G4QInelastic::PostStDoIt: cost="<<cost<<G4endl;
1315	if (cost> 1.) cost= 1.;
1316	else if(cost<-1.) cost=-1.;
1317	}
1318	G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,fM); // 4mom of the recoilNucleus
1319	scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scatteredLepton
1320	G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-ml)*.01);
1321	if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
1322	{
1323	G4cerr<<"G4QC::PSDI:t4M="<<tot4M<<",lM="<<ml<<",rM="<<fM<<",cost="<<cost<<G4endl;
1324	//G4Exception("G4QInelastic::PostStepDoIt:","027",FatalException,"ElasticDecay");
1325	}
1326	#ifdef debug
1327	G4cout<<"G4QCol::PStDoI:l4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
1328	#endif
1329	// ----------------------------------------------------
1330	G4ParticleDefinition* theDefinition=0; // Prototype of a particle for E-Secondaries
1331	// Fill scattered lepton
1332	if (scatPDG==-11) theDefinition = G4Positron::Positron();
1333	else if(scatPDG== 11) theDefinition = G4Electron::Electron();
1334	else if(scatPDG== 13) theDefinition = G4MuonMinus::MuonMinus();
1335	else if(scatPDG==-13) theDefinition = G4MuonPlus::MuonPlus();
1336	//else if(scatPDG== 15) theDefinition = G4TauMinus::TauMinus();
1337	//else if(scatPDG==-15) theDefinition = G4TauPlus::TauPlus();
1338	if (scatPDG==-12) theDefinition = G4AntiNeutrinoE::AntiNeutrinoE();
1339	else if(scatPDG== 12) theDefinition = G4NeutrinoE::NeutrinoE();
1340	else if(scatPDG== 14) theDefinition = G4NeutrinoMu::NeutrinoMu();
1341	else if(scatPDG==-14) theDefinition = G4AntiNeutrinoMu::AntiNeutrinoMu();
1342	//else if(scatPDG== 16) theDefinition = G4NeutrinoTau::NeutrinoTau();
1343	//else if(scatPDG==-16) theDefinition = G4AntiNeutrinoTau::AntiNeutrinoTau();
1344	else G4cout<<"-Warning-G4QInelastic::PostStDoIt: UnknownLepton="<<scatPDG<<G4endl;
1345	G4DynamicParticle* theScL = new G4DynamicParticle(theDefinition,scat4M);
1346	G4Track* scatLep = new G4Track(theScL, localtime, position ); // scattered
1347	scatLep->SetWeight(weight); // weighted
1348	scatLep->SetTouchableHandle(trTouchable); // residual
1349	aParticleChange.AddSecondary(scatLep); // lepton
1350	// Fill residual nucleus
1351	if (fintPDG==90000001) theDefinition = G4Neutron::Neutron(); // neutron
1352	else if(fintPDG==90001000) theDefinition = G4Proton::Proton(); // proton
1353	else // ion
1354	{
1355	G4int fm=static_cast<G4int>(fintPDG/1000000); // Strange part
1356	G4int ZN=fintPDG-1000000*fm;
1357	G4int rZ=static_cast<G4int>(ZN/1000);
1358	G4int rA=ZN-999*rZ;
1359	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
1360	}
1361	G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,reco4M);
1362	G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1363	scatReN->SetWeight(weight); // weighted
1364	scatReN->SetTouchableHandle(trTouchable); // residual
1365	aParticleChange.AddSecondary(scatReN); // nucleus
1366	delete output;
1367	return G4VDiscreteProcess::PostStepDoIt(track, step);
1368	}
1369	}
1370	//
1371	// quasi-elastic (& pickup process) for p+A(Z,N)
1372	//
1373	}
1374	else if ((projPDG == 2212 \|\| projPDG == 2112) && Z > 0 && N > 0)
1375	//else if(2>3)
1376	{
1377	G4bool ncap = true; // As a prototype: the capture process have not happened
1378	// @@ in the same way the fission reaction can be added for heavy nuclei
1379	if(momentum<500. && projPDG == 2112) // @@ It's reasonable to add proton capture too !
1380	{
1381	G4QNeutronCaptureRatio* capMan=G4QNeutronCaptureRatio::GetPointer();
1382	if(G4UniformRand() <= capMan->GetRatio(momentum, Z, N)) ncap = false; // Make capture
1383	else if(momentum<100.) // Try the production threshold
1384	{
1385	G4int tZ=Z;
1386	G4int tN=N;
1387	if(projPDG == 2212) tZ++;
1388	else tN++; // @@ Now a neutron is the only alternative
1389	G4LorentzVector tot4M=G4LorentzVector(0.,0.,0.,tgM)+proj4M;
1390	G4double totM2=tot4M.m2();
1391	G4QPDGCode FakeP(2112);
1392	//G4int m1p=tZ-1;
1393	G4int m2n=tN-2;
1394	// @@ For the secondary proton the Coulomb barrier should be taken into account
1395	//G4double mRP= FakeP.GetNuclMass(m1p,tN,0); // Residual to the secondary proton
1396	G4double mR2N= FakeP.GetNuclMass(tZ,m2n,0); // Residual to two secondary neutrons
1397	//G4double mRAl= FakeP.GetNuclMass(tZ-2,m2n,0); // Residual to the secondary alpha
1398	//G4double mR2N= FakeP.GetNuclMass(m1p,tN-1,0); // Residual to the p+n secondaries
1399	G4double minM=mR2N+mNeut+mNeut;
1400	// Compare with other possibilities (sciped for acceleration)
1401	if(totM2 < minM*minM) ncap = false; // Make capture (n,gamma)
1402	if(!ncap)
1403	{
1404	G4int resPDG=targPDG;
1405	if(projPDG == 2212) resPDG+=1000;
1406	else resPDG++; // @@ Now neutron is theOnlyAlternative
1407	G4double rM=G4QPDGCode(resPDG).GetMass(); // Mass of the residual nucleus
1408	G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1409	G4LorentzVector g4M=G4LorentzVector(0.,0.,0.,0.); // 4mom of the gamma
1410	#ifdef qedebug
1411	G4cout<<"G4QCollis::PStDoIt: (n,gamma), tM="<<tM<<", R="<<rM<<G4endl;
1412	#endif
1413	if(!G4QHadron(tot4M).DecayIn2(r4M, g4M))
1414	{
1415	G4cerr<<"G4QCol::PostStDoIt: tM="<<std::sqrt(totM2)<<" < rM="<<rM<<G4endl;
1416	throw G4QException("G4QInelastic::HadronizeQuasm:Can'tDec TotNuc->ResNuc+gam");
1417	}
1418	#ifdef qedebug
1419	G4cout<<"G4QCol::PStDoIt:n->g,R="<<r4M.rho()<<r4M<<",G="<<g4M.rho()<<g4M<<G4endl;
1420	#endif
1421	aParticleChange.ProposeEnergy(0.); // @@ ??
1422	aParticleChange.ProposeTrackStatus(fStopAndKill);// projectile nucleon is killed
1423	aParticleChange.SetNumberOfSecondaries(2);
1424	// Fill gamma
1425	G4ParticleDefinition* theDefinition = G4Gamma::Gamma();
1426	G4DynamicParticle* theGam = new G4DynamicParticle(theDefinition,g4M);
1427	G4Track* capGamma = new G4Track(theGam, localtime, position );
1428	capGamma->SetWeight(weight);
1429	capGamma->SetTouchableHandle(trTouchable);
1430	aParticleChange.AddSecondary(capGamma);
1431	// ----------------------------------------------------
1432	// Fill residual nucleus
1433	G4int tA=tZ+tN;
1434	if (resPDG==90000001) theDefinition = G4Neutron::Neutron();
1435	else if(resPDG==90001000) theDefinition = G4Proton::Proton();
1436	else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(tZ,tA,0,tZ);
1437	G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1438	G4Track* scatReN = new G4Track(theReN, localtime, position ); // weighted
1439	scatReN->SetWeight(weight); // residual
1440	scatReN->SetTouchableHandle(trTouchable); // recoil
1441	aParticleChange.AddSecondary(scatReN); // nucleus
1442	delete output;
1443	return G4VDiscreteProcess::PostStepDoIt(track, step);
1444	}
1445	}
1446	}
1447	if(ncap)
1448	{
1449	// Now the quasi-free and PickUp processes:
1450	G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer();
1451	std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, projPDG, Z, N);
1452	G4double qepart=fief.first*fief.second; // @@ charge exchange is not included @@
1453	// Keep QE part for the quasi-free nucleons
1454	const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
1455	G4double clProb[lCl]={.65,.85,.95};// N/P,D,t/He3,He4, integroProbab for .65,.2,.1,.05
1456	if(qepart>0.45) qepart=.45; // Quasielastic part is too large - shrink
1457	//else qepart/=clProb[0]; // Add corresponding number of 2N, 3N, & 4N clusters
1458	qepart=qepart/clProb[0]-qepart;// Add QE for 2N, 3N, & 4N clusters (N is made in G4QF)
1459	G4double pickup=1.-qepart; // Estimate the rest of the cross-section
1460	G4double thresh=100.;
1461	if(momentum > thresh) pickup*=50./momentum/std::pow(G4double(Z+N),third);// 50 is Par
1462	// pickup = 0.; // To exclude the pickup process
1463	if (N) pickup+=qepart;
1464	else pickup =qepart;
1465	G4double rnd=G4UniformRand();
1466	#ifdef ldebug
1467	G4cout<<"-->G4QCol::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<")]="<<qepart
1468	<<", pickup="<<pickup<<G4endl;
1469	#endif
1470	if(rnd<pickup) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim,CoulBar,PauliBl
1471	{
1472	G4double dmom=91.; // Fermi momentum (proto default for a deuteron)
1473	G4double fmm=287.; // @@ Can be A-dependent
1474	if(Z>1\|\|N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
1475	// Values to be redefined for quasi-elastic
1476	G4LorentzVector r4M(0.,0.,0.,0.); // Prototype of 4mom of the residual nucleus
1477	G4LorentzVector n4M(0.,0.,0.,0.); // Prototype of 4mom of the quasi-cluster
1478	G4int nPDG=90000001; // Prototype for quasiCluster mass calculation
1479	G4int restPDG=targPDG; // Prototype should be reduced by quasiCluster
1480	G4int rA=Z+N-1; // Prototype for the residualNucl definition
1481	G4int rZ=Z; // residZ: OK for the quasi-free neutron
1482	G4int nA=1; // Prototype for the quasi-cluster definition
1483	G4int nZ=0; // nA=1,nZ=0: OK for the quasi-free neutron
1484	G4double qM=mNeut; // Free mass of the quasi-free cluster
1485	G4int A = Z + N; // Baryon number of the nucleus
1486	if(rnd<qepart)
1487	{
1488	// ===> First decay a nucleus in a nucleon and a residual (A-1) nucleus
1489	// Calculate clusterProbabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters)
1490	G4double base=1.; // Base for randomization (can be reduced by totZ & totN)
1491	G4int max=lCl; // Number of boundaries (can be reduced by totZ & totN)
1492	// Take into account that at least one nucleon must be left !
1493	if (Z<2 \|\| N<2 \|\| A < 6) base = clProb[--max]; // Alpha cluster is impossible
1494	if ( (Z > 1 && N < 2) \|\| (Z < 2 && N > 1) )
1495	base=(clProb[max]+clProb[max-1])/2; // t or He3 is impossible
1496	if ( (Z < 2 && N < 2) \|\| A < 5) base=clProb[--max];// He3&t clusters are impossible
1497	if(A<3) base=clProb[--max]; // Deuteron cluster is impossible
1498	//G4int cln=0; // Cluster#0 (Default for the selectedNucleon)
1499	G4int cln=1; // Cluster#1 (Default for the selectedDeutron)
1500	//if(max) // Not only nucleons are possible
1501	if(max>1) // Not only deuterons are possible
1502	{
1503	base-=clProb[0]; // Exclude scattering on QF Nucleon
1504	G4double ran=+clProb[0]+base*G4UniformRand(); // Base can be reduced
1505	G4int ic=1; // Start from the smallest cluster boundary
1506	if(max>1) while(ic<max) if(ran>clProb[ic++]) cln=ic;
1507	}
1508	G4ParticleDefinition* theDefinition=0;// Prototype for qfNucleon
1509	G4bool cp1 = cln+2==A; // A=ClusterBN+1 condition
1510	if(!cln \|\| cp1) // Split in nucleon + (A-1) with Fermi momentum
1511	{
1512	G4int nln=0;
1513	if(cln==2) nln=1; // @@ only for cp1: t/He3 choice from A=4
1514	// mass(A)=tM. Calculate masses of A-1 (rM) and mN (mNeut or mProt bounded mass)
1515	if ( ((!cln \|\| cln == 2) && G4UniformRand()*(A-cln) > (N-nln)) \|\|
1516	((cln == 3 \|\| cln == 1) && Z > N) )
1517	{
1518	nPDG=90001000; // Update quasi-free nucleon PDGCode to P
1519	nZ=1; // Change charge of the quasiFree nucleon
1520	qM=mProt; // Update quasi-free nucleon mass
1521	rZ--; // Reduce the residual Z
1522	restPDG-=1000; // Reduce the residual PDGCode
1523	}
1524	else restPDG--;
1525	G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1526	G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1527	r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1528	G4double rM2=rM*rM;
1529	G4double nM=std::sqrt(rM2+tMtM-(tM+tM)std::sqrt(rM2+dmom*dmom));// M of q-nucleon
1530	n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon
1531	#ifdef qedebug
1532	G4cout<<"G4QCollis::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1533	#endif
1534	if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1535	{
1536	G4cerr<<"G4QCol::PostStDoIt: M="<<tM<<"<rM="<<rM<<"+nM="<<nM<<"="<<rM+nM<<G4endl;
1537	throw G4QException("G4QInelastic::HadronizeQuasm:Can'tDec totNuc->QENuc+ResNuc");
1538	}
1539	#ifdef qedebug
1540	G4cout<<"G4QCol::PStDoIt:QE-N,RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
1541	#endif
1542	if(cp1 && cln) // Quasi-cluster case: swap the output
1543	{
1544	qM=rM; // Scattering will be made on a cluster
1545	nln=nPDG;
1546	nPDG=restPDG;
1547	restPDG=nln;
1548	t4M=n4M;
1549	n4M=r4M;
1550	r4M=t4M;
1551	nln=nZ;
1552	nZ=rZ;
1553	rZ=nln;
1554	nln=nA;
1555	nA=rA;
1556	rA=nln;
1557	}
1558	}
1559	else // Split the cluster (with or w/o "Fermi motion" and "Fermi decay")
1560	{
1561	if(cln==1)
1562	{
1563	nPDG=90001001; // Deuteron
1564	qM=mDeut;
1565	nA=2;
1566	nZ=1;
1567	restPDG-=1001;
1568	// Recalculate dmom
1569	G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1570	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1571	dmom=sumv.mag();
1572	}
1573	else if(cln==2)
1574	{
1575	nA=3;
1576	if(G4UniformRand()*(A-2)>(N-1)) // He3
1577	{
1578	nPDG=90002001;
1579	qM=mHel3;
1580	nZ=2;
1581	restPDG-=2001;
1582	}
1583	else // tritium
1584	{
1585	nPDG=90001002;
1586	qM=mTrit;
1587	nZ=1;
1588	restPDG-=1002;
1589	}
1590	// Recalculate dmom
1591	G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1592	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection()+
1593	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1594	dmom=sumv.mag();
1595	}
1596	else
1597	{
1598	nPDG=90002002; // Alpha
1599	qM=mAlph;
1600	nA=4;
1601	nZ=2;
1602	restPDG-=2002;
1603	G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1604	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection()+
1605	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection()+
1606	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1607	dmom=sumv.mag();
1608	}
1609	rA=A-nA;
1610	rZ=Z-nZ;
1611	// This is a simple case of cluster at rest
1612	//G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1613	//r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1614	//n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster
1615	// --- End of the "simple case of cluster at rest"
1616	// Make a fake quasi-Fermi distribution for clusters (clusters are not at rest)
1617	G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1618	G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1619	r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1620	G4double rM2=rM*rM;
1621	G4double nM=std::sqrt(rM2+tMtM-(tM+tM)std::sqrt(rM2+dmom*dmom));// M of q-cluster
1622	n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon
1623	#ifdef qedebug
1624	G4cout<<"G4QCollis::PStDoIt:QEC,p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1625	#endif
1626	if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1627	{
1628	G4cerr<<"G4QCol::PostStDoIt: M="<<tM<<"<rM="<<rM<<"+cM="<<nM<<"="<<rM+nM<<G4endl;
1629	throw G4QException("G4QInelastic::HadronizeQuasm:Can'tDec totNuc->QEClu+ResNuc");
1630	}
1631	// --- End of the moving cluster implementation ---
1632	#ifdef qedebug
1633	G4cout<<"G4QCol::PStDoIt:QEC,RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
1634	#endif
1635	}
1636	G4LorentzVector s4M=n4M+proj4M; // Tot 4-momentum for scattering
1637	G4double prjM2 = proj4M.m2(); // Squared mass of the projectile
1638	G4double prjM = std::sqrt(prjM2); // @@ Get from pPDG (?)
1639	G4double minM = prjM+qM; // Min mass sum for the final products
1640	G4double cmM2 =s4M.m2();
1641	if(cmM2>minM*minM)
1642	{
1643	#ifdef qedebug
1644	G4cout<<"G4QCol::PStDoIt:*Enter,cmM2="<<cmM2<<" > minM2="<<minMminM<<G4endl;
1645	#endif
1646	// Estimate and randomize charge-exchange with a quasi-free cluster
1647	G4bool chex=false; // Flag of the charge exchange scattering
1648	G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true
1649	// This is reserved for the charge-exchange scattering (to help neutron spectras)
1650	//if(cln&&!cp1 &&(projPDG==2212&&rA>rZ \|\| projPDG==2112&&rZ>1))// @@ Use proj chex
1651	if(2>3) // @@ charge exchange is not implemented yet @@
1652	{
1653	#ifdef qedebug
1654	G4cout<<"G4QCol::PStDoIt:-Enter,P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
1655	#endif
1656	G4double tprM=mProt;
1657	G4double tprM2=mProt2;
1658	G4int tprPDG=2212;
1659	G4int tresPDG=restPDG+999;
1660	if(projPDG==2212)
1661	{
1662	projpt=G4Neutron::Neutron();
1663	tprM=mNeut;
1664	tprM2=mNeut2;
1665	tprPDG=2112;
1666	tresPDG=restPDG-999;
1667	}
1668	minM=tprM+qM;
1669	G4double efE=(cmM2-tprM2-qM*qM)/(qM+qM);
1670	G4double efP=std::sqrt(efE*efE-tprM2);
1671	G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!)
1672	#ifdef qedebug
1673	G4cout<<"G4QCol::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
1674	#endif
1675	if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl)) // minM is redefined
1676	{
1677	projPDG=tprPDG;
1678	prjM=tprM;
1679	G4double rM=G4QPDGCode(tresPDG).GetMass();// Mass of the residual nucleus
1680	r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1681	n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster
1682	chex=true; // Confirm charge exchange scattering
1683	}
1684	}
1685	//
1686	std::pair<G4LorentzVector,G4LorentzVector> sctout=qfMan->Scatter(nPDG, n4M,
1687	projPDG, proj4M);
1688	#ifdef qedebug
1689	G4cout<<"G4QCollis::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
1690	<<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl;
1691	#endif
1692	aParticleChange.ProposeLocalEnergyDeposit(0.); // Everything is in particles
1693	// @@ @@ @@ Coulomb barriers must be checked !! @@ @@ @@ Skip if not
1694	if(chex) // ==> Projectile is changed: fill everything to secondaries
1695	{
1696	aParticleChange.ProposeEnergy(0.); // @@ ??
1697	aParticleChange.ProposeTrackStatus(fStopAndKill); // projectile nucleon is killed
1698	aParticleChange.SetNumberOfSecondaries(3);
1699	G4DynamicParticle* thePrH = new G4DynamicParticle(projpt,sctout.second);
1700	G4Track* scatPrH = new G4Track(thePrH, localtime, position ); // scattered & chex
1701	scatPrH->SetWeight(weight); // weighted
1702	scatPrH->SetTouchableHandle(trTouchable); // projectile
1703	aParticleChange.AddSecondary(scatPrH); // hadron
1704	}
1705	else // ==> The leading particle is filled to the updated projectilee
1706	{
1707	aParticleChange.SetNumberOfSecondaries(2); // @@ if proj=leading
1708	G4double ldT=(sctout.second).e()-prjM; // kin Energy of scat project.
1709	aParticleChange.ProposeEnergy(ldT); // Change the kin Energy
1710	G4ThreeVector ldV=(sctout.second).vect(); // Change momentum direction
1711	aParticleChange.ProposeMomentumDirection(ldV/ldV.mag());
1712	aParticleChange.ProposeTrackStatus(fAlive);
1713	}
1714	// ---------------------------------------------------------
1715	// Fill scattered quasi-free nucleon/fragment
1716	if (nPDG==90000001) theDefinition = G4Neutron::Neutron();
1717	else if(nPDG==90001000) theDefinition = G4Proton::Proton();
1718	else if(nZ>0 && nA>1)
1719	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);
1720	#ifdef debug
1721	else G4cout<<"-Warning_G4QCol::PSD:scatqfPDG="<<nPDG<<",Z="<<nZ<<",A="<<nA<<G4endl;
1722	#endif
1723	if(nZ>0 && nA>0)
1724	{
1725	G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,sctout.first);
1726	G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // scattered
1727	scatQFN->SetWeight(weight); // weighted
1728	scatQFN->SetTouchableHandle(trTouchable); // quasi-free
1729	aParticleChange.AddSecondary(scatQFN); // nucleon/cluster
1730	}
1731	// ----------------------------------------------------
1732	// Fill residual nucleus
1733	if (restPDG==90000001) theDefinition = G4Neutron::Neutron();
1734	else if(restPDG==90001000) theDefinition = G4Proton::Proton();
1735	else if(rZ>0 && rA>1)
1736	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
1737	#ifdef debug
1738	else G4cout<<"-Warning_G4QCol::PSD:resPDG="<<restPDG<<",Z="<<rZ<<",A="<<rA<<G4endl;
1739	#endif
1740	if(rZ>0 && rA>0)
1741	{
1742	G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1743	G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1744	scatReN->SetWeight(weight); // weighted
1745	scatReN->SetTouchableHandle(trTouchable); // residual
1746	aParticleChange.AddSecondary(scatReN); // nucleus
1747	}
1748	delete output;
1749	return G4VDiscreteProcess::PostStepDoIt(track, step);
1750	}
1751	#ifdef qedebug
1752	else G4cout<<"G4QCol::PSD:OUT, M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl;
1753	#endif
1754	} // end of proton quasi-elastic (including QE on NF)
1755	else // @@ make cost-condition @@ Pickup process (pickup 1 or 2 n and make d or t)
1756	{
1757	if(projPDG==2212) restPDG--; // Res. nucleus PDG is one neutron less than targ. PDG
1758	else
1759	{
1760	restPDG-=1000; // Res. nucleus PDG is one proton less than targ. PDG
1761	rZ--; // Reduce the charge of the residual nucleus
1762	}
1763	G4double mPUF=mDeut; // Prototype of the mass of the created pickup fragment
1764	G4ParticleDefinition* theDefinition = G4Deuteron::Deuteron(); // Default: make d
1765	//theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion
1766	G4bool din=false; // Flag of picking up 2 neutrons to create t (tritium)(p)
1767	G4bool dip=false; // Flag of picking up 2 protons to create t (tritium) (n)
1768	G4bool pin=false; // Flag of picking up 1 n + 1 p to create He3/t (p/n)
1769	G4bool hin=false; // Flag of PickUp creation of alpha (He4) (p/n)
1770	G4double frn=G4UniformRand();
1771	if(N>2 && frn > 0.95) // .95 is a parameter to pickup 2N (to cteate t or He3)
1772	{
1773	if(projPDG==2212)
1774	{
1775	if(N>1 && G4UniformRand()(Z+.5(N-1)) > Z)
1776	{
1777	theDefinition = G4Triton::Triton();
1778	mPUF=mTrit; // The mass of the created pickup fragment is triton
1779	restPDG--; // Res. nucleus PDG is two neutrons less than target PDG
1780	din=true;
1781	}
1782	else
1783	{
1784	theDefinition = G4He3::He3();
1785	mPUF=mHel3; // The mass of the created pickup fragment is Helium3
1786	restPDG-=1000; // Res. nucleus PDG is two neutrons less than target PDG
1787	rZ--; // Reduce the charge of the residual nucleus
1788	pin=true;
1789	}
1790	}
1791	else // Neutron projectile (2112)
1792	{
1793	if(Z>1 && G4UniformRand()(N+.5(Z-1)) > N)
1794	{
1795	theDefinition = G4He3::He3();
1796	mPUF=mHel3; // The mass of the created pickup fragment is Helium3
1797	restPDG-=1000; // Res. nucleus PDG is two neutrons less than target PDG
1798	rZ--; // Reduce the charge of the residual nucleus
1799	dip=true;
1800	}
1801	else
1802	{
1803	theDefinition = G4Triton::Triton();
1804	mPUF=mTrit; // The mass of the created pickup fragment is triton
1805	restPDG--; // Res. nucleus PDG is two neutrons less than target PDG
1806	pin=true;
1807	}
1808	}
1809	rA--; // Reduce the baryon number of the residual nucleus
1810	// Recalculate dmom
1811	G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1812	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1813	dmom=sumv.mag();
1814	if(Z>1 && frn > 0.99) // .99 is a parameter to pickup 2n1p \|\| 1n2p & create alpha
1815	{
1816	theDefinition = G4Alpha::Alpha();
1817	if((din && projPDG==2212) \|\| (pin && projPDG==2112))
1818	{
1819	restPDG-=1000; // Residual nucleus PDG is reduced to create alpha
1820	rZ--; // Reduce the charge of the residual nucleus
1821	}
1822	else if((pin && projPDG==2212) \|\| (dip && projPDG==2112)) restPDG--;
1823	else G4cout<<"-Warning-G4QCol::PSD: PickUp logic error, proj="<<projPDG<<G4endl;
1824	hin=true;
1825	mPUF=mAlph; // The mass of the created pickup fragment (alpha)
1826	rA--; // Reduce the baryon number of the residual nucleus
1827	// Recalculate dmom
1828	sumv += fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1829	dmom=sumv.mag();
1830	}
1831	}
1832	G4double mPUF2=mPUF*mPUF;
1833	G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1834	G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1835	G4double rM2=rM*rM; // Squared mass of the residual nucleus
1836	G4double nM2=rM2+tMtM-(tM+tM)std::sqrt(rM2+dmom*dmom);// M2 of boundQF-nucleon(2N)
1837	if(nM2 < 0) nM2=0.;
1838	G4double nM=std::sqrt(nM2); // M of bQF-nucleon
1839	G4double den2=(dmom*dmom+nM2); // squared energy of bQF-nucleon
1840	G4double den=std::sqrt(den2); // energy of bQF-nucleon
1841	#ifdef qedebug
1842	G4cout<<"G4QCollis::PStDoIt:PiU,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1843	#endif
1844	G4double qp=momentum*dmom;
1845	G4double srm=proj4M.e()*den;
1846	G4double mNucl2=mProt2;
1847	if(projPDG == 2112) mNucl2=mNeut2;
1848	G4double cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp);
1849	G4int ict=0;
1850	while(std::fabs(cost)>1. && ict<3)
1851	{
1852	dmom=91.; // Fermi momentum (proDefault for a deuteron)
1853	if(Z>1\|\|N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
1854	if(din \|\| pin \|\| dip) // Recalculate dmom for the final t/He3
1855	{
1856	G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1857	fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1858	if(hin)
1859	sumv+=fmmstd::pow(-std::log(G4UniformRand()),third)G4RandomDirection();
1860	dmom=sumv.mag();
1861	}
1862	nM2=rM2+tMtM-(tM+tM)std::sqrt(rM2+dmom*dmom); // M2 of bQF-nucleon
1863	if(nM2 < 0) nM2=0.;
1864	nM=std::sqrt(nM2); // M of bQF-nucleon
1865	den2=(dmom*dmom+nM2); // squared energy of bQF-nucleon
1866	den=std::sqrt(den2); // energy of bQF-nucleon
1867	qp=momentum*dmom;
1868	srm=proj4M.e()*den;
1869	cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp);
1870	ict++;
1871	#ifdef ldebug
1872	if(ict>2)G4cout<<"G4QCollis::PStDoIt:i="<<ict<<",d="<<dmom<<",ct="<<cost<<G4endl;
1873	#endif
1874	}
1875	if(std::fabs(cost)<=1.)
1876	{// ---- @@ From here can be a MF for QF nucleon extraction (if used by others)
1877	G4ThreeVector vdir = proj4M.vect(); // 3-Vector in the projectile direction
1878	G4ThreeVector vx(0.,0.,1.); // ProtoOrt in theDirection of the projectile
1879	G4ThreeVector vy(0.,1.,0.); // First ProtoOrt orthogonal to the direction
1880	G4ThreeVector vz(1.,0.,0.); // SecondProtoOrt orthoganal to the direction
1881	if(vdir.mag2() > 0.) // the projectile isn't at rest
1882	{
1883	vx = vdir.unit(); // Ort in the direction of the projectile
1884	G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
1885	vy = vv.unit(); // First ort orthogonal to the proj direction
1886	vz = vx.cross(vy); // Second ort orthoganal to the projDirection
1887	}
1888	// ---- @@ End of possible MF (Similar is in G4QEnvironment)
1889	G4double tdom=dmomstd::sqrt(1-costcost);// Transversal component of the Fermi-Mom
1890	G4double phi=twopi*G4UniformRand(); // Phi of the Fermi-Mom
1891	G4ThreeVector vedm=vxdmomcost+vytdomstd::sin(phi)+vztdomstd::cos(phi);// bQFN
1892	G4LorentzVector bqf(vedm,den); // 4-mom of the bQF nucleon
1893	r4M=t4M-bqf; // 4mom of the residual nucleus
1894	if(std::fabs(r4M.m()-rM)>.001) G4cout<<"G4QCol::PSD:rM="<<rM<<"#"<<r4M.m()<<G4endl;
1895	G4LorentzVector f4M=proj4M+bqf; // Prototype of 4-mom of Deuteron
1896	if(std::fabs(f4M.m()-mPUF)>.001)
1897	G4cout<<"G4QCol::PSD:mDeut="<<mPUF<<" # "<<f4M.m()<<G4endl;
1898	aParticleChange.ProposeEnergy(0.); // @@ ??
1899	aParticleChange.ProposeTrackStatus(fStopAndKill);// projectile nucleon is killed
1900	aParticleChange.SetNumberOfSecondaries(2);
1901	// Fill pickup hadron
1902	G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,f4M);
1903	G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // pickup
1904	scatQFN->SetWeight(weight); // weighted
1905	scatQFN->SetTouchableHandle(trTouchable); // nuclear
1906	aParticleChange.AddSecondary(scatQFN); // cluster
1907	#ifdef pickupd
1908	G4cout<<"G4QCol::PStDoIt: f="<<theDefinition<<",f4M="<<f4M.m()<<f4M<<G4endl;
1909	#endif
1910	// ----------------------------------------------------
1911	// Fill residual nucleus
1912	if (restPDG==90000001) theDefinition = G4Neutron::Neutron();
1913	else if(restPDG==90001000) theDefinition = G4Proton::Proton();
1914	else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion
1915	G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1916	G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1917	scatReN->SetWeight(weight); // weighted
1918	scatReN->SetTouchableHandle(trTouchable); // residual
1919	aParticleChange.AddSecondary(scatReN); // nucleus
1920	#ifdef pickupd
1921	G4cout<<"G4QCol::PStDoIt: rZ="<<rZ<<", rA="<<rA<<",r4M="<<r4M.m()<<r4M<<G4endl;
1922	#endif
1923	delete output;
1924	return G4VDiscreteProcess::PostStepDoIt(track, step);
1925	}
1926	}
1927	} // end of the quasi-elastic decoupled process for the neucleon projectile
1928	} // end the PickUp and quasi-elastic decoupled processes (capture is made before them)
1929	} // end lepto-nuclear, neutrino-nuclear, proton/neutron decoupled processes
1930	EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
1931	if(absMom) EnMomConservation+=lead4M; // Add E/M of leading System
1932	#ifdef debug
1933	G4cout<<"G4QInelastic::PostStDoIt:before St="<<aParticleChange.GetTrackStatus()<<G4endl;
1934	#endif
1935
1936	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin
1937	//G4bool elF=false; // Flag of the ellastic scattering is "false" by default
1938	//G4double eWei=1.;
1939	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- End
1940	#ifdef ldebug
1941	G4cout<<"^^G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
1942	#endif
1943	const G4QNucleus targNuc(Z,N); // Define the target nucleus
1944	if(projPDG>9000000)
1945	{
1946	delete output; // Before the output creation
1947	G4QNucleus projNuc(proj4M,projPDG); // Define the projectile nucleus
1948	G4QIonIonCollision IonIon(projNuc, targNuc); // Define deep-inelastic ion-ion reaction
1949	#ifdef debug
1950	G4cout<<"G4QCol::PStDoIt: ProjNuc="<<projPDG<<proj4M<<", TargNuc="<<targPDG<<G4endl;
1951	#endif
1952	try {output = IonIon.Fragment();} // DINR AA interaction products
1953	catch (G4QException& error)
1954	{
1955	G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl;
1956	G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash");
1957	}
1958	}
1959	else // --> A projectile hadron
1960	{
1961	G4int maxCn=7;
1962	G4int atCn=0; // Attempts counter
1963	G4bool inel=false;
1964	while (inel==false && ++atCn <= maxCn) // To evoid elastic final state
1965	{
1966	G4int outN=output->size();
1967	if(outN) // The output is not empty
1968	{
1969	for_each(output->begin(), output->end(), DeleteQHadron());
1970	output->clear();
1971	}
1972	delete output; // Before the new output creation
1973	const G4QHadron projH(projPDG,proj4M); // Define the projectile particle
1974	G4QFragmentation DINR(targNuc, projH); // Define deep-inel. h-A reaction
1975	#ifdef debug
1976	G4cout<<"G4QInelastic::PStDoIt:Proj="<<projPDG<<proj4M<<",TargNuc="<<targPDG<<G4endl;
1977	#endif
1978	try {output = DINR.Fragment();} // DINR hA fragmentation
1979	catch (G4QException& error)
1980	{
1981	G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl;
1982	G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash");
1983	}
1984	outN=output->size();
1985	#ifdef debug
1986	G4cout<<"G4QInelastic::PostStepDoIt: At# "<<atCn<<", nSec="<<outN<<", fPDG="
1987	<<(*output)[0]->GetPDGCode()<<", pPDG="<< projPDG<<G4endl;
1988	#endif
1989	inel=true;
1990	if(outN < 2)
1991	{
1992	G4cout<<"-Warning-G4QInelastic::PostStepDoIt: nSec="<<outN<<", At# "<<atCn<<G4endl;
1993	inel=false;
1994	}
1995	else if(outN==2 && (*output)[0]->GetPDGCode() == projPDG) inel=false;
1996	#ifdef debug
1997	if(atCn==maxCn) G4cout<<"-Warn-G4QCol::PostStDoIt:mAt="<<atCn<<" is reached"<<G4endl;
1998	#endif
1999	}
2000	}
2001	//
2002	// --- the scattered hadron with changed nature can be added here ---
2003	if(scat)
2004	{
2005	G4QHadron* scatHadron = new G4QHadron(scatPDG,scat4M);
2006	output->push_back(scatHadron);
2007	}
2008	G4int qNH=0;
2009	if(leadhs) qNH=leadhs->size();
2010	if(absMom)
2011	{
2012	if(qNH) for(G4int iq=0; iq<qNH; iq++)
2013	{
2014	G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron
2015	output->push_back(loh);
2016	}
2017	if(leadhs) delete leadhs;
2018	leadhs=0;
2019	}
2020	else
2021	{
2022	if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
2023	if(leadhs) delete leadhs;
2024	leadhs=0;
2025	}
2026	// ------------- From here the secondaries are filled -------------------------
2027	G4int tNH = output->size(); // A#of hadrons in the output
2028	aParticleChange.SetNumberOfSecondaries(tNH);
2029	// Now add nuclear fragments
2030	#ifdef debug
2031	G4cout<<"G4QInelastic::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
2032	#endif
2033	#ifdef ppdebug
2034	if(absMom)G4cout<<"G4QInelastic::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
2035	#endif
2036	G4int nOut=output->size(); // Real length of the output @@ Temporary
2037	if(tNH==1 && !scat) // @@ Temporary. Find out why it happened!
2038	{
2039	G4cout<<"-Warning-G4QInelastic::PostStepDoIt: 1 secondary! absMom="<<absMom;
2040	if(absMom) G4cout<<", qNH="<<qNH;
2041	G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
2042	G4cout<<G4endl;
2043	tNH=0;
2044	delete output->operator[](0); // delete the creazy hadron
2045	output->pop_back(); // clean up the output vector
2046	}
2047	if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QInelastic::PostStepDoIt: 2 # "<<nOut<<G4endl;
2048	// Deal with ParticleChange final state interface to GEANT4 output of the process
2049	//if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
2050	if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
2051	{
2052	// Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
2053	// Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
2054	// The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
2055	G4QHadron* hadr=(*output)[i]; // Pointer to the output hadron
2056	G4int PDGCode = hadr->GetPDGCode();
2057	G4int nFrag = hadr->GetNFragments();
2058	#ifdef pdebug
2059	G4cout<<"G4QInelastic::PostStepDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag
2060	<<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
2061	#endif
2062	if(nFrag) // Skip intermediate (decayed) hadrons
2063	{
2064	#ifdef debug
2065	G4cout<<"G4QInelastic::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
2066	#endif
2067	delete hadr;
2068	continue;
2069	}
2070	G4DynamicParticle* theSec = new G4DynamicParticle;
2071	G4ParticleDefinition* theDefinition;
2072	if (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
2073	else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
2074	else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
2075	else if(PDGCode==311 \|\| PDGCode==-311)
2076	{
2077	if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L
2078	else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
2079	}
2080	else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
2081	else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
2082	else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
2083	else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
2084	else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
2085	else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
2086	{
2087	G4int aZ = hadr->GetCharge();
2088	G4int aA = hadr->GetBaryonNumber();
2089	#ifdef pdebug
2090	G4cout<<"G4QInelastic::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
2091	#endif
2092	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
2093	}
2094	//else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
2095	else
2096	{
2097	#ifdef pdebug
2098	G4cout<<"G4QInelastic::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
2099	#endif
2100	theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
2101	#ifdef pdebug
2102	G4cout<<"G4QInelastic::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
2103	#endif
2104	}
2105	if(!theDefinition)
2106	{
2107	#ifdef debug
2108	G4cout<<"---Warning---G4QInelastic::PostStepDoIt: drop PDG="<<PDGCode<<G4endl;
2109	#endif
2110	delete hadr;
2111	continue;
2112	}
2113	#ifdef pdebug
2114	G4cout<<"G4QInelastic::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
2115	#endif
2116	theSec->SetDefinition(theDefinition);
2117	G4LorentzVector h4M=hadr->Get4Momentum();
2118	EnMomConservation-=h4M;
2119	#ifdef tdebug
2120	G4cout<<"G4QCollis::PSDI:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
2121	#endif
2122	#ifdef debug
2123	G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
2124	#endif
2125	theSec->Set4Momentum(h4M); // ^
2126	delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
2127	#ifdef debug
2128	G4ThreeVector curD=theSec->GetMomentumDirection(); // ^
2129	G4double curM=theSec->GetMass(); // \|
2130	G4double curE=theSec->GetKineticEnergy()+curM; // ^
2131	G4cout<<"G4QCollis::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;// \|
2132	#endif
2133	G4Track* aNewTrack = new G4Track(theSec, localtime, position ); // ^
2134	aNewTrack->SetWeight(weight); // weighted \|
2135	aNewTrack->SetTouchableHandle(trTouchable); // \|
2136	aParticleChange.AddSecondary( aNewTrack ); // \|
2137	#ifdef debug
2138	G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<" is done."<<G4endl; // \|
2139	#endif
2140	} // \|
2141	delete output; // instances of the G4QHadrons from the output are already deleted above +
2142	if(leadhs) // To satisfy Valgrind ( How can that be?)
2143	{
2144	G4int qNH=leadhs->size();
2145	if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
2146	delete leadhs;
2147	leadhs=0;
2148	}
2149	#ifdef debug
2150	G4cout<<"G4QInelastic::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
2151	#endif
2152	if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15)
2153	aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the absorbed particle
2154	#ifdef pdebug
2155	G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()
2156	<<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl;
2157	#endif
2158	#ifdef ldebug
2159	G4cout<<"G4QInelastic::PostStepDoIt:* PostStepDoIt is done *, P="<<aProjPDG<<", St="
2160	<<aParticleChange.GetTrackStatus()<<G4endl;
2161	#endif
2162	return G4VDiscreteProcess::PostStepDoIt(track, step);
2163	}
2164
2165	std::pair<G4double,G4double> G4QInelastic::Random2DDirection()
2166	{
2167	G4double sp=0; // sin(phi)
2168	G4double cp=1.; // cos(phi)
2169	G4double r2=2.; // to enter the loop
2170	while(r2>1. \|\| r2<.0001) // pi/4 efficiency
2171	{
2172	G4double s=G4UniformRand();
2173	G4double c=G4UniformRand();
2174	sp=1.-s-s;
2175	cp=1.-c-c;
2176	r2=spsp+cpcp;
2177	}
2178	G4double norm=std::sqrt(r2);
2179	return std::make_pair(sp/norm,cp/norm);
2180	}

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

Download in other formats: