Context Navigation

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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago
update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 105.0 KB

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

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

Download in other formats: