Context Navigation

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

Last change on this file since 968 was 962, checked in by garnier, 15 years ago
update processes
File size: 82.6 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: G4QCollision.cc,v 1.28 2008/10/02 21:10:07 dennis Exp $
27	// GEANT4 tag $Name: geant4-09-02-ref-02 $
28	//
29	// ---------------- G4QCollision class -----------------
30	// by Mikhail Kossov, December 2003.
31	// G4QCollision class of the CHIPS Simulation Branch in GEANT4
32	// ---------------------------------------------------------------
33	// ****************************************************************************************
34	// ******** This CLASS is temporary moved from the photolepton_hadron directory *******
35	// ****************************************************************************************
36
37	//#define debug
38	//#define pdebug
39	//#define ppdebug
40	//#define qedebug
41
42	#include "G4QCollision.hh"
43
44	// Initialization of static vectors
45	std::vector<G4int> G4QCollision::ElementZ; // Z of the element(i) in theLastCalc
46	std::vector<G4double> G4QCollision::ElProbInMat; // SumProbabilityElements in Material
47	std::vector<std::vector<G4int>*> G4QCollision::ElIsoN; // N of isotope(j) of Element(i)
48	std::vector<std::vector<G4double>*>G4QCollision::IsoProbInEl;//SumProbabIsotopes inElementI
49
50	G4QCollision::G4QCollision(const G4String& processName):
51	G4VDiscreteProcess(processName, fHadronic)
52	{
53	#ifdef debug
54	G4cout<<"G4QCollision::Constructor is called"<<G4endl;
55	#endif
56	if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
57	SetProcessSubType(fHadronInelastic);
58	//G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPSWorld (234 part.max)
59	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
60	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
61	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
62	//@@ Initialize here the G4QuasmonString parameters
63	}
64
65	G4bool G4QCollision::manualFlag=false; // If false then standard parameters are used
66	G4double G4QCollision::Temperature=180.; // Critical Temperature (sensitive at High En)
67	G4double G4QCollision::SSin2Gluons=0.3; // Supression of s-quarks (in respect to u&d)
68	G4double G4QCollision::EtaEtaprime=0.3; // Supression of eta mesons (gg->qq/3g->qq)
69	G4double G4QCollision::freeNuc=0.5; // Percentage of free nucleons on the surface
70	G4double G4QCollision::freeDib=0.05; // Percentage of free diBaryons on the surface
71	G4double G4QCollision::clustProb=5.; // Nuclear clusterization parameter
72	G4double G4QCollision::mediRatio=10.; // medium/vacuum hadronization ratio
73	G4int G4QCollision::nPartCWorld=152; // The#of particles initialized in CHIPS World
74	G4double G4QCollision::SolidAngle=0.5; // Part of Solid Angle to capture (@@A-dep.)
75	G4bool G4QCollision::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
76	G4double G4QCollision::PiPrThresh=141.4; // Pion Production Threshold for gammas
77	G4double G4QCollision::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
78	G4double G4QCollision::DiNuclMass=1880.; // DoubleNucleon Mass for VirtualNormalization
79	G4double G4QCollision::photNucBias=1.; // BiasingParameter for photo(e,mu,tau)Nuclear
80	G4double G4QCollision::weakNucBias=1.; // BiasingParameter for ChargedCurrents(nu,mu)
81
82	void G4QCollision::SetManual() {manualFlag=true;}
83	void G4QCollision::SetStandard() {manualFlag=false;}
84
85	// Fill the private parameters
86	void G4QCollision::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
87	G4double fN, G4double fD, G4double cP, G4double mR,
88	G4int nParCW, G4double solAn, G4bool efFlag,
89	G4double piThresh, G4double mpisq, G4double dinum)
90	{// =============================================================================
91	Temperature=temper;
92	SSin2Gluons=ssin2g;
93	EtaEtaprime=etaetap;
94	freeNuc=fN;
95	freeDib=fD;
96	clustProb=cP;
97	mediRatio=mR;
98	nPartCWorld = nParCW;
99	EnergyFlux=efFlag;
100	SolidAngle=solAn;
101	PiPrThresh=piThresh;
102	M2ShiftVir=mpisq;
103	DiNuclMass=dinum;
104	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
105	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
106	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
107	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
108	}
109
110	void G4QCollision::SetPhotNucBias(G4double phnB) {photNucBias=phnB;}
111	void G4QCollision::SetWeakNucBias(G4double ccnB) {weakNucBias=ccnB;}
112
113	// Destructor
114
115	G4QCollision::~G4QCollision() {}
116
117
118	G4LorentzVector G4QCollision::GetEnegryMomentumConservation()
119	{
120	return EnMomConservation;
121	}
122
123	G4int G4QCollision::GetNumberOfNeutronsInTarget()
124	{
125	return nOfNeutrons;
126	}
127
128	// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
129	// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
130	// ******** All CHIPS cross sections are calculated in the surface units **********
131	G4double G4QCollision::GetMeanFreePath(const G4Track& aTrack,G4double,G4ForceCondition* Fc)
132	{
133	#ifdef debug
134	G4cout<<"G4QCollision::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
135	#endif
136	*Fc = NotForced;
137	#ifdef debug
138	G4cout<<"G4QCollision::GetMeanFreePath: Before GetDynPart"<<G4endl;
139	#endif
140	const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
141	#ifdef debug
142	G4cout<<"G4QCollision::GetMeanFreePath: Before GetDef"<<G4endl;
143	#endif
144	G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
145	if( !IsApplicable(*incidentParticleDefinition))
146	G4cout<<"-W-G4QCollision::GetMeanFreePath called for not implemented particle"<<G4endl;
147	// Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
148	G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
149	#ifdef debug
150	G4cout<<"G4QCollis::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
151	#endif
152	const G4Material* material = aTrack.GetMaterial(); // Get the current material
153	const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
154	const G4ElementVector* theElementVector = material->GetElementVector();
155	G4int nE=material->GetNumberOfElements();
156	#ifdef debug
157	G4cout<<"G4QCollision::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
158	#endif
159	G4bool leptoNuc=false; // By default the reaction is not lepto-nuclear
160	G4VQCrossSection* CSmanager=0;
161	G4VQCrossSection* CSmanager2=0;
162	G4int pPDG=0;
163	if(incidentParticleDefinition == G4Proton::Proton())
164	{
165	CSmanager=G4QProtonNuclearCrossSection::GetPointer();
166	pPDG=2212;
167	}
168	else if(incidentParticleDefinition == G4Gamma::Gamma())
169	{
170	CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
171	pPDG=22;
172	}
173	else if(incidentParticleDefinition == G4Electron::Electron() \|\|
174	incidentParticleDefinition == G4Positron::Positron())
175	{
176	CSmanager=G4QElectronNuclearCrossSection::GetPointer();
177	leptoNuc=true;
178	pPDG=11;
179	}
180	else if(incidentParticleDefinition == G4MuonPlus::MuonPlus() \|\|
181	incidentParticleDefinition == G4MuonMinus::MuonMinus())
182	{
183	CSmanager=G4QMuonNuclearCrossSection::GetPointer();
184	leptoNuc=true;
185	pPDG=13;
186	}
187	else if(incidentParticleDefinition == G4TauPlus::TauPlus() \|\|
188	incidentParticleDefinition == G4TauMinus::TauMinus())
189	{
190	CSmanager=G4QTauNuclearCrossSection::GetPointer();
191	leptoNuc=true;
192	pPDG=15;
193	}
194	else if(incidentParticleDefinition == G4NeutrinoMu::NeutrinoMu() )
195	{
196	CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
197	CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
198	leptoNuc=true;
199	pPDG=14;
200	}
201	else if(incidentParticleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMu() )
202	{
203	CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
204	CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
205	leptoNuc=true;
206	pPDG=-14;
207	}
208	else if(incidentParticleDefinition == G4NeutrinoE::NeutrinoE() )
209	{
210	CSmanager=G4QNuENuclearCrossSection::GetPointer();
211	CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
212	leptoNuc=true;
213	pPDG=12;
214	}
215	else if(incidentParticleDefinition == G4AntiNeutrinoE::AntiNeutrinoE() )
216	{
217	CSmanager=G4QANuENuclearCrossSection::GetPointer();
218	CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
219	leptoNuc=true;
220	pPDG=-12;
221	}
222	else G4cout<<"G4QCollision::GetMeanFreePath:Particle isn't implemented in CHIPS"<<G4endl;
223
224	G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
225	G4double sigma=0.; // Sums over elements for the material
226	G4int IPIE=IsoProbInEl.size(); // How many old elements?
227	if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
228	{
229	std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
230	SPI->clear();
231	delete SPI;
232	std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
233	IsN->clear();
234	delete IsN;
235	}
236	ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
237	ElementZ.clear(); // Clear the body vector for Z of Elements
238	IsoProbInEl.clear(); // Clear the body vector for SPI
239	ElIsoN.clear(); // Clear the body vector for N of Isotopes
240	for(G4int i=0; i<nE; ++i)
241	{
242	G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
243	G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
244	ElementZ.push_back(Z); // Remember Z of the Element
245	G4int isoSize=0; // The default for the isoVectorLength is 0
246	G4int indEl=0; // Index of non-trivial element or 0(default)
247	G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
248	if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
249	#ifdef debug
250	G4cout<<"G4QCollision::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
251	#endif
252	if(isoSize) // The Element has non-trivial abumdance set
253	{
254	indEl=pElement->GetIndex()+1; // Index of the non-trivial element
255	if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
256	{
257	std::vector<std::pair<G4int,G4double>> newAbund =
258	new std::vector<std::pair<G4int,G4double>*>;
259	G4double* abuVector=pElement->GetRelativeAbundanceVector();
260	for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
261	{
262	G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
263	if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QCollision::GetMeanFreePath"
264	<<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
265	G4double abund=abuVector[j];
266	std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
267	#ifdef debug
268	G4cout<<"G4QCollision::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
269	#endif
270	newAbund->push_back(pr);
271	}
272	#ifdef debug
273	G4cout<<"G4QCollision::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
274	#endif
275	indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
276	for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
277	delete newAbund; // Was "new" in the beginning of the name space
278	}
279	}
280	std::vector<std::pair<G4int,G4double>> cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
281	std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
282	IsoProbInEl.push_back(SPI);
283	std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
284	ElIsoN.push_back(IsN);
285	G4int nIs=cs->size(); // A#Of Isotopes in the Element
286	G4double susi=0.; // sum of CS over isotopes
287	if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
288	{
289	std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
290	G4int N=curIs->first; // #of Neuterons in the isotope j of El i
291	IsN->push_back(N); // Remember Min N for the Element
292	G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
293	if(CSmanager2)CSI+=CSmanager2->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i)nu,nu
294	#ifdef debug
295	G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
296	#endif
297	curIs->second = CSI;
298	susi+=CSI; // Make a sum per isotopes
299	SPI->push_back(susi); // Remember summed cross-section
300	} // End of temporary initialization of the cross sections in the G4QIsotope singeltone
301	sigma+=Isotopes->GetMeanCrossSection(Z,indEl)NOfNucPerVolume[i];//SUM(MeanCSNOfNperV)
302	ElProbInMat.push_back(sigma);
303	} // End of LOOP over Elements
304	#ifdef debug
305	G4cout<<"G4QCol::GetMeanFrPa: S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
306	#endif
307	// Check that cross section is not zero and return the mean free path
308	if(photNucBias!=1.) if(incidentParticleDefinition == G4Gamma::Gamma() \|\|
309	incidentParticleDefinition == G4MuonPlus::MuonPlus() \|\|
310	incidentParticleDefinition == G4MuonMinus::MuonMinus() \|\|
311	incidentParticleDefinition == G4Electron::Electron() \|\|
312	incidentParticleDefinition == G4Positron::Positron() \|\|
313	incidentParticleDefinition == G4TauMinus::TauMinus() \|\|
314	incidentParticleDefinition == G4TauPlus::TauPlus() )
315	sigma*=photNucBias;
316	if(weakNucBias!=1.) if(incidentParticleDefinition==G4NeutrinoE::NeutrinoE() \|\|
317	incidentParticleDefinition==G4AntiNeutrinoE::AntiNeutrinoE() \|\|
318	incidentParticleDefinition==G4NeutrinoTau::NeutrinoTau() \|\|
319	incidentParticleDefinition==G4AntiNeutrinoTau::AntiNeutrinoTau()\|\|
320	incidentParticleDefinition==G4NeutrinoMu::NeutrinoMu() \|\|
321	incidentParticleDefinition==G4AntiNeutrinoMu::AntiNeutrinoMu() )
322	sigma*=weakNucBias;
323	if(sigma > 0.) return 1./sigma; // Mean path [distance]
324	return DBL_MAX;
325	}
326
327
328	G4bool G4QCollision::IsApplicable(const G4ParticleDefinition& particle)
329	{
330	if (particle == *( G4MuonPlus::MuonPlus() )) return true;
331	else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
332	else if (particle == *( G4TauPlus::TauPlus() )) return true;
333	else if (particle == *( G4TauMinus::TauMinus() )) return true;
334	else if (particle == *( G4Electron::Electron() )) return true;
335	else if (particle == *( G4Positron::Positron() )) return true;
336	else if (particle == *( G4Gamma::Gamma() )) return true;
337	else if (particle == *( G4Proton::Proton() )) return true;
338	else if (particle == *( G4AntiNeutrinoE::AntiNeutrinoE() )) return true;
339	else if (particle == *( G4NeutrinoE::NeutrinoE() )) return true;
340	else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
341	else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
342	//else if (particle == *( G4Neutron::Neutron() )) return true;
343	//else if (particle == *( G4PionMinus::PionMinus() )) return true;
344	//else if (particle == *( G4PionPlus::PionPlus() )) return true;
345	//else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
346	//else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
347	//else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
348	//else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
349	//else if (particle == *( G4Lambda::Lambda() )) return true;
350	//else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
351	//else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
352	//else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
353	//else if (particle == *( G4XiMinus::XiMinus() )) return true;
354	//else if (particle == *( G4XiZero::XiZero() )) return true;
355	//else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
356	//else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
357	//else if (particle == *( G4AntiProton::AntiProton() )) return true;
358	//else if (particle == *(G4AntiNeutrinoTau::AntiNeutrinoTau())) return true;
359	//else if (particle == *( G4NeutrinoTau::NeutrinoTau() )) return true;
360	#ifdef debug
361	G4cout<<"***G4QCollision::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
362	#endif
363	return false;
364	}
365
366	G4VParticleChange* G4QCollision::PostStepDoIt(const G4Track& track, const G4Step& step)
367	{
368	static const G4double third = 1./3.;
369	static const G4double me=G4Electron::Electron()->GetPDGMass(); // electron mass
370	static const G4double me2=me*me; // squared electron mass
371	static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass
372	static const G4double mu2=mu*mu; // squared muon mass
373	static const G4double mt=G4TauMinus::TauMinus()->GetPDGMass(); // tau mass
374	static const G4double mt2=mt*mt; // squared tau mass
375	//static const G4double dpi=M_PI+M_PI; // 2pi (for Phi distr.) changed to twopi*
376	static const G4double mNeut= G4QPDGCode(2112).GetMass();
377	static const G4double mNeut2= mNeut*mNeut;
378	static const G4double mProt= G4QPDGCode(2212).GetMass();
379	static const G4double mProt2= mProt*mProt;
380	static const G4double dM=mProt+mNeut; // doubled nucleon mass
381	static const G4double hdM=dM/2.; // M of the "nucleon"
382	static const G4double hdM2=hdM*hdM; // M2 of the "nucleon"
383	static const G4double mPi0 = G4QPDGCode(111).GetMass();
384	static const G4double mPi0s= mPi0*mPi0;
385	static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
386	static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
387	static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
388	static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
389	static const G4double mPi = G4QPDGCode(211).GetMass();
390	static const G4double tmPi = mPi+mPi; // Doubled mass of the charged pion
391	static const G4double stmPi= tmPi*tmPi; // Squared Doubled mass of the charged pion
392	static const G4double mPPi = mPi+mProt; // Delta threshold
393	static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2
394	//static const G4double mDel2= 1400*1400; // Delta up threshold for W2 (in MeV^2)
395	// Static definitions for electrons (nu,e) -----------------------------------------
396	static const G4double meN = mNeut+me;
397	static const G4double meN2= meN*meN;
398	static const G4double fmeN= 4mNeut2me2;
399	static const G4double mesN= mNeut2+me2;
400	static const G4double meP = mProt+me;
401	static const G4double meP2= meP*meP;
402	static const G4double fmeP= 4mProt2me2;
403	static const G4double mesP= mProt2+me2;
404	static const G4double medM= me2/dM; // for x limit
405	static const G4double meD = mPPi+me; // Multiperipheral threshold
406	static const G4double meD2= meD*meD;
407	// Static definitions for muons (nu,mu) -----------------------------------------
408	static const G4double muN = mNeut+mu;
409	static const G4double muN2= muN*muN;
410	static const G4double fmuN= 4mNeut2mu2;
411	static const G4double musN= mNeut2+mu2;
412	static const G4double muP = mProt+mu;
413	static const G4double muP2= muP*muP; // +
414	static const G4double fmuP= 4mProt2mu2; // +
415	static const G4double musP= mProt2+mu2;
416	static const G4double mudM= mu2/dM; // for x limit
417	static const G4double muD = mPPi+mu; // Multiperipheral threshold
418	static const G4double muD2= muD*muD;
419	// Static definitions for muons (nu,nu) -----------------------------------------
420	//static const G4double nuN = mNeut;
421	//static const G4double nuN2= mNeut2;
422	//static const G4double fnuN= 0.;
423	//static const G4double nusN= mNeut2;
424	//static const G4double nuP = mProt;
425	//static const G4double nuP2= mProt2;
426	//static const G4double fnuP= 0.;
427	//static const G4double nusP= mProt2;
428	//static const G4double nudM= 0.; // for x limit
429	//static const G4double nuD = mPPi; // Multiperipheral threshold
430	//static const G4double nuD2= mPPi2;
431	//-------------------------------------------------------------------------------------
432	static G4bool CWinit = true; // CHIPS Warld needs to be initted
433	if(CWinit)
434	{
435	CWinit=false;
436	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
437	}
438	//-------------------------------------------------------------------------------------
439	const G4DynamicParticle* projHadron = track.GetDynamicParticle();
440	const G4ParticleDefinition* particle=projHadron->GetDefinition();
441	#ifdef debug
442	G4cout<<"G4QCollision::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
443	#endif
444	G4ForceCondition cond=NotForced;
445	GetMeanFreePath(track, 1., &cond);
446	#ifdef debug
447	G4cout<<"G4QCollision::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
448	#endif
449	G4bool scat=false; // No CHEX in proj scattering
450	G4int scatPDG=0; // Must be filled if true (CHEX)
451	G4LorentzVector proj4M=projHadron->Get4Momentum(); // 4-momentum of the projectile (IU?)
452	G4LorentzVector scat4M=proj4M; // Must be filled if true
453	G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
454	G4double Momentum=proj4M.rho();
455	if(std::fabs(Momentum-momentum)>.001)
456	G4cerr<<"*G4QCollision::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl;
457	#ifdef debug
458	G4double mp=proj4M.m();
459	G4cout<<"G4QCollis::PostStepDoIt:called, P="<<Momentum<<"="<<momentum<<",m="<<mp<<G4endl;
460	#endif
461	if (!IsApplicable(*particle)) // Check applicability
462	{
463	G4cerr<<"G4QCollision::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
464	<<G4endl;
465	return 0;
466	}
467	const G4Material* material = track.GetMaterial(); // Get the current material
468	G4int Z=0;
469	const G4ElementVector* theElementVector = material->GetElementVector();
470	G4int nE=material->GetNumberOfElements();
471	#ifdef debug
472	G4cout<<"G4QCollision::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
473	#endif
474	G4int projPDG=0; // PDG Code prototype for the captured hadron
475	// Not all these particles are implemented yet (see Is Applicable)
476	if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
477	else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
478	else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14;
479	else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14;
480	else if (particle == G4Electron::Electron() ) projPDG= 11;
481	else if (particle == G4Positron::Positron() ) projPDG= -11;
482	else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12;
483	else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12;
484	else if (particle == G4Gamma::Gamma() ) projPDG= 22;
485	else if (particle == G4Proton::Proton() ) projPDG= 2212;
486	else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
487	else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
488	else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
489	else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 2112;
490	else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
491	else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
492	else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310;
493	else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
494	else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
495	else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
496	else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
497	else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
498	else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
499	else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
500	else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
501	else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
502	else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
503	else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
504	else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
505	else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
506	G4int aProjPDG=std::abs(projPDG);
507	#ifdef debug
508	G4int prPDG=particle->GetPDGEncoding();
509	G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
510	#endif
511	if(!projPDG)
512	{
513	G4cerr<<"---Warning---G4QCollision::PostStepDoIt:Undefined interacting hadron"<<G4endl;
514	return 0;
515	}
516	G4int EPIM=ElProbInMat.size();
517	#ifdef debug
518	G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
519	#endif
520	G4int i=0;
521	if(EPIM>1)
522	{
523	G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
524	for(i=0; i<nE; ++i)
525	{
526	#ifdef debug
527	G4cout<<"G4QCollision::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
528	#endif
529	if (rnd<ElProbInMat[i]) break;
530	}
531	if(i>=nE) i=nE-1; // Top limit for the Element
532	}
533	G4Element* pElement=(*theElementVector)[i];
534	Z=static_cast<G4int>(pElement->GetZ());
535	#ifdef debug
536	G4cout<<"G4QCollision::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
537	#endif
538	if(Z<=0)
539	{
540	G4cerr<<"---Warning---G4QCollision::PostStepDoIt: Element with Z="<<Z<<G4endl;
541	if(Z<0) return 0;
542	}
543	std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
544	std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
545	G4int nofIsot=SPI->size(); // #of isotopes in the element i
546	#ifdef debug
547	G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
548	#endif
549	G4int j=0;
550	if(nofIsot>1)
551	{
552	G4double rndI=(SPI)[nofIsot-1]G4UniformRand(); // Randomize the isotop of the Element
553	for(j=0; j<nofIsot; ++j)
554	{
555	#ifdef debug
556	G4cout<<"G4QCollision::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
557	#endif
558	if(rndI < (*SPI)[j]) break;
559	}
560	if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
561	}
562	G4int N =(*IsN)[j]; ; // Randomized number of neutrons
563	#ifdef debug
564	G4cout<<"G4QCollision::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<G4endl;
565	#endif
566	if(N<0)
567	{
568	G4cerr<<"-Warning-G4QCollision::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
569	return 0;
570	}
571	nOfNeutrons=N; // Remember it for the energy-momentum check
572	G4double dd=0.025;
573	G4double am=Z+N;
574	G4double sr=std::sqrt(am);
575	G4double dsr=0.01*(sr+sr);
576	if(dsr<dd)dsr=dd;
577	if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // ManualP
578	else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP ClustPars
579	else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.);//Pi- ClustPars
580	#ifdef debug
581	G4cout<<"G4QCollision::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
582	#endif
583	if(N<0)
584	{
585	G4cerr<<"---Warning---G4QCollision::PostStepDoIt:Element with N="<<N<< G4endl;
586	return 0;
587	}
588	aParticleChange.Initialize(track);
589	G4double weight = track.GetWeight();
590	if(photNucBias!=1.) weight/=photNucBias;
591	else if(weakNucBias!=1.) weight/=weakNucBias;
592	G4double localtime = track.GetGlobalTime();
593	G4ThreeVector position = track.GetPosition();
594	G4TouchableHandle trTouchable = track.GetTouchableHandle();
595	//
596	G4int targPDG=90000000+Z*1000+N; // PDG Code of the target nucleus
597	G4QPDGCode targQPDG(targPDG);
598	G4double tgM=targQPDG.GetMass(); // Target mass
599	G4double tM=tgM; // Target mass (copy to be changed)
600	G4QHadronVector* output=new G4QHadronVector;// Prototype of EnvironOutput G4QHadronVector
601	G4double absMom = 0.; // Prototype of absorbed by nucleus Moment
602	G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVectorum
603	G4LorentzVector lead4M(0.,0.,0.,0.); // Prototype of LeadingQ 4-momentum
604
605	//
606	// Leptons with photonuclear
607	// Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?)
608	//
609	if (aProjPDG == 11 \|\| aProjPDG == 13 \|\| aProjPDG == 15) {
610
611	#ifdef debug
612	G4cout<<"G4QCollision::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl;
613	#endif
614	G4double kinEnergy= projHadron->GetKineticEnergy();
615	G4ParticleMomentum dir = projHadron->GetMomentumDirection();
616	G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
617	G4double ml=me;
618	G4double ml2=me2;
619	if(aProjPDG== 13)
620	{
621	CSmanager=G4QMuonNuclearCrossSection::GetPointer();
622	ml=mu;
623	ml2=mu2;
624	}
625	if(aProjPDG== 15)
626	{
627	CSmanager=G4QTauNuclearCrossSection::GetPointer();
628	ml=mt;
629	ml2=mt2;
630	}
631	// @@ Probably this is not necessary any more (?)
632	G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,aProjPDG);// Recalculate XS
633	// @@ check a possibility to separate p, n, or alpha (!)
634	if(xSec <= 0.) // The cross-section is 0 -> Do Nothing
635	{
636	#ifdef debug
637	G4cerr<<"---OUT---G4QCollision::PSDoIt: Called for zero Cross-section"<<G4endl;
638	#endif
639	//Do Nothing Action insead of the reaction
640	aParticleChange.ProposeEnergy(kinEnergy);
641	aParticleChange.ProposeLocalEnergyDeposit(0.);
642	aParticleChange.ProposeMomentumDirection(dir);
643	aParticleChange.ProposeTrackStatus(fAlive);
644	return G4VDiscreteProcess::PostStepDoIt(track,step);
645	}
646	G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
647	#ifdef debug
648	G4cout<<"G4QCol::PStDoIt: kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
649	#endif
650	if( kinEnergy < photonEnergy \|\| photonEnergy < 0.)
651	{
652	//Do Nothing Action insead of the reaction
653	G4cerr<<"--G4QCollision::PSDoIt: photE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl;
654	aParticleChange.ProposeEnergy(kinEnergy);
655	aParticleChange.ProposeLocalEnergyDeposit(0.);
656	aParticleChange.ProposeMomentumDirection(dir);
657	aParticleChange.ProposeTrackStatus(fAlive);
658	return G4VDiscreteProcess::PostStepDoIt(track,step);
659	}
660	G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart
661	G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon
662	if(tM<999.) W-=mPi0+mPi0s/dM; // Pion production threshold for a nucleon target
663	if(W<0.)
664	{
665	//Do Nothing Action insead of the reaction
666	#ifdef debug
667	G4cout<<"--G4QCollision::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl;
668	#endif
669	aParticleChange.ProposeEnergy(kinEnergy);
670	aParticleChange.ProposeLocalEnergyDeposit(0.);
671	aParticleChange.ProposeMomentumDirection(dir);
672	aParticleChange.ProposeTrackStatus(fAlive);
673	return G4VDiscreteProcess::PostStepDoIt(track,step);
674	}
675	// Update G4VParticleChange for the scattered muon
676	G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer();
677	G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy,Z,N,22);//Integrated XS
678	G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N, 22); // Real XS
679	G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2);
680	if(sigNuG4UniformRand()>sigKrndFraction)
681	{
682	//Do NothingToDo Action insead of the reaction
683	#ifdef debug
684	G4cout<<"-DoNoth-G4QCollision::PostStepDoIt: probab. correction - DoNothing"<<G4endl;
685	#endif
686	aParticleChange.ProposeEnergy(kinEnergy);
687	aParticleChange.ProposeLocalEnergyDeposit(0.);
688	aParticleChange.ProposeMomentumDirection(dir);
689	aParticleChange.ProposeTrackStatus(fAlive);
690	return G4VDiscreteProcess::PostStepDoIt(track,step);
691	}
692	G4double iniE=kinEnergy+ml; // Initial total energy of the lepton
693	G4double finE=iniE-photonEnergy; // Final total energy of the lepton
694	#ifdef pdebug
695	G4cout<<"G4QCollision::PoStDoIt:E="<<iniE<<",lE="<<finE<<"-"<<ml<<"="<<finE-ml<<G4endl;
696	#endif
697	aParticleChange.ProposeEnergy(finE-ml);
698	if(finE<=ml) // Secondary lepton (e/mu/tau) at rest disappears
699	{
700	aParticleChange.ProposeEnergy(0.);
701	if(aProjPDG== 11) aParticleChange.ProposeTrackStatus(fStopAndKill);
702	else aParticleChange.ProposeTrackStatus(fStopButAlive);
703	aParticleChange.ProposeMomentumDirection(dir);
704	}
705	else aParticleChange.ProposeTrackStatus(fAlive);
706	G4double iniP=std::sqrt(iniE*iniE-ml2); // Initial momentum of the electron
707	G4double finP=std::sqrt(finE*finE-ml2); // Final momentum of the electron
708	G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton)
709	#ifdef pdebug
710	G4cout<<"G4QC::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
711	#endif
712	if(cost>1.) cost=1.; // To avoid the accuracy of calculation problem
713	if(cost<-1.) cost=-1.; // To avoid the accuracy of calculation problem
714	//
715	// Scatter the lepton ( @@ make the same thing for real photons)
716	// At this point we have photonEnergy and photonQ2 (with notDefinedPhi)->SelectProjPart
717	G4double absEn = std::pow(am,third)*GeV; // @@(b) Mean Energy Absorbed by a Nucleus
718	//if(am>1 && absEn < photonEnergy) // --> the absorption of energy can happen
719	if(absEn < photonEnergy) // --> the absorption of energy can happen
720	{
721	G4double abtEn = absEn+hdM; // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N)
722	G4double abEn2 = abtEn*abtEn; // Squared absorbed Energy + MN
723	G4double abMo2 = abEn2-hdM2; // Squared absorbed Momentum of compound system
724	G4double phEn2 = photonEnergy*photonEnergy;
725	G4double phMo2 = phEn2+photonQ2; // Squared momentum of primary virtual photon
726	G4double phMo = std::sqrt(phMo2); // Momentum of the primary virtual photon
727	absMom = std::sqrt(abMo2); // Absorbed Momentum
728	if(absMom < phMo) // --> the absorption of momentum can happen
729	{
730	G4double dEn = photonEnergy - absEn; // Leading energy
731	G4double dMo = phMo - absMom; // Leading momentum
732	G4double sF = dEndEn - dModMo;// s of leading particle
733	#ifdef ppdebug
734	G4cout<<"-PhotoAbsorption-G4QCol::PStDoIt:sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
735	#endif
736	if(sF > stmPi) // --> Leading fragmentation is possible
737	{
738	photonEnergy = absEn; // New value of the photon energy
739	photonQ2=abMo2-absEn*absEn; // New value of the photon Q2
740	absEn = dEn; // Put energy of leading particle to absEn (!)
741	}
742	else absMom=0.; // Flag that nothing has happened
743	}
744	else absMom=0.; // Flag that nothing has happened
745	}
746	// ------------- End of ProjPart selection
747	//
748	// Scattering in respect to the derection of the incident muon is made impicitly:
749	G4ThreeVector ort=dir.orthogonal(); // Not normed orthogonal vector (!) (to dir)
750	G4ThreeVector ortx = ort.unit(); // First unit vector orthogonal to the direction
751	G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction
752	G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
753	G4double phi=twopi*G4UniformRand(); // phi of scattered electron
754	G4double sinx=sint*std::sin(phi); // x perpendicular component
755	G4double siny=sint*std::cos(phi); // y perpendicular component
756	G4ThreeVector findir=costdir+sinxortx+siny*orty;
757	aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton
758	#ifdef pdebug
759	G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
760	<<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl;
761	#endif
762	G4ThreeVector photon3M=iniPdir-finPfindir;// 3D total momentum of photon
763	if(absMom) // Photon must be reduced & LeadingSyst fragmented
764	{
765	G4double ptm=photon3M.mag(); // 3M of the virtual photon
766	#ifdef ppdebug
767	G4cout<<"-Absorption-G4QCollision::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
768	<<iniPdir<<", eFin3M="<<finPfindir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl;
769	#endif
770	G4ThreeVector lead3M=photon3M*(ptm-absMom)/ptm; // Keep the direction for leading Q
771	photon3M-=lead3M; // Reduced photon Momentum (photEn already = absEn)
772	proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System
773	#ifdef ppdebug
774	G4cout<<"-->G4QC::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
775	#endif
776	lead4M=proj4M; // Remember 4-mom for the total 4-momentum
777	G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+
778	try // \|
779	{ // \|
780	delete leadhs; // \|
781	G4QNucleus vac(90000000); // \|
782	leadhs=pan->Fragment(vac,1); // DELETED after it is copied to output vector \|
783	} // \|
784	catch (G4QException& error) // \|
785	{ // \|
786	G4cerr<<"***G4QCollision::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//\|
787	G4Exception("G4QCollision::PostStepDoIt:","72",FatalException,"QuasmonCrash"); //\|
788	} // \|
789	delete pan; // Delete the Nuclear Environment <----<---+
790	#ifdef ppdebug
791	G4cout<<"G4QCol::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
792	<<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl;
793	#endif
794	}
795	projPDG=22;
796	proj4M=G4LorentzVector(photon3M,photonEnergy);
797	#ifdef debug
798	G4cout<<"G4QCollision::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
799	<<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl;
800	#endif
801
802	//
803	// neutrinoNuclear interactions (nu_e, nu_mu only)
804	//
805	} else if (aProjPDG == 12 \|\| aProjPDG == 14) {
806
807	G4double kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino
808	G4double dKinE=kinEnergy+kinEnergy; // doubled energy for s calculation
809	#ifdef debug
810	G4cout<<"G4QCollision::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
811	#endif
812	G4ParticleMomentum dir = projHadron->GetMomentumDirection(); // unit vector
813	G4double ml = mu;
814	G4double ml2 = mu2;
815	//G4double mlN = muN;
816	G4double mlN2= muN2;
817	G4double fmlN= fmuN;
818	G4double mlsN= musN;
819	//G4double mlP = muP;
820	G4double mlP2= muP2;
821	G4double fmlP= fmuP;
822	G4double mlsP= musP;
823	G4double mldM= mudM;
824	//G4double mlD = muD;
825	G4double mlD2= muD2;
826	if(aProjPDG==12)
827	{
828	ml = me;
829	ml2 = me2;
830	//mlN = meN;
831	mlN2= meN2;
832	fmlN= fmeN;
833	mlsN= mesN;
834	//mlP = meP;
835	mlP2= meP2;
836	fmlP= fmeP;
837	mlsP= mesP;
838	mldM= medM;
839	//mlD = meD;
840	mlD2= meD2;
841	}
842	G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l)
843	G4VQCrossSection* CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); // (nu,nu)
844	proj4M=G4LorentzVector(dir*kinEnergy,kinEnergy); // temporary
845	G4bool nuanu=true;
846	scatPDG=13; // Prototype = secondary scattered mu-
847	if(projPDG==-14)
848	{
849	nuanu=false; // Anti-neutrino
850	CSmanager=G4QANuMuNuclearCrossSection::GetPointer(); // (anu,mu+) CC @@ open
851	CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
852	scatPDG=-13; // secondary scattered mu+
853	}
854	else if(projPDG==12)
855	{
856	CSmanager=G4QNuENuclearCrossSection::GetPointer(); // @@ open (only CC is changed)
857	scatPDG=11; // secondary scattered e-
858	}
859	else if(projPDG==-12)
860	{
861	nuanu=false; // anti-neutrino
862	CSmanager=G4QANuENuclearCrossSection::GetPointer(); // (anu,e+) CC @@ open
863	CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
864	scatPDG=-11; // secondary scattered e+
865	}
866	// @@ Probably this is not necessary any more
867	G4double xSec1=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG); //Recalculate XS
868	G4double xSec2=CSmanager2->GetCrossSection(false,Momentum,Z,N,projPDG);//Recalculate XS
869	G4double xSec=xSec1+xSec2;
870	// @@ check a possibility to separate p, n, or alpha (!)
871	if(xSec <= 0.) // The cross-section = 0 -> Do Nothing
872	{
873	G4cerr<<"G4QCollision::PSDoIt:nuE="<<kinEnergy<<",X1="<<xSec1<<",X2="<<xSec2<<G4endl;
874	//Do Nothing Action insead of the reaction
875	aParticleChange.ProposeEnergy(kinEnergy);
876	aParticleChange.ProposeLocalEnergyDeposit(0.);
877	aParticleChange.ProposeMomentumDirection(dir);
878	aParticleChange.ProposeTrackStatus(fAlive);
879	return G4VDiscreteProcess::PostStepDoIt(track,step);
880	}
881	G4bool secnu=false;
882	if(xSec*G4UniformRand()>xSec1) // recover neutrino/antineutrino
883	{
884	if(scatPDG>0) scatPDG++;
885	else scatPDG--;
886	secnu=true;
887	}
888	scat=true; // event with changed scattered projectile
889	G4double totCS1 = CSmanager->GetLastTOTCS(); // the last total cross section1(isotope?)
890	G4double totCS2 = CSmanager2->GetLastTOTCS();// the last total cross section2(isotope?)
891	G4double totCS = totCS1+totCS2; // the last total cross section (isotope?)
892	if(std::fabs(xSec-totCS*millibarn)/xSec>.0001)
893	G4cout<<"-Warning-G4QCollision::PostStepDoIt: xS="<<xSec<<"# CS="<<totCS<<G4endl;
894	G4double qelCS1 = CSmanager->GetLastQELCS(); // the last quasi-elastic cross section1
895	G4double qelCS2 = CSmanager2->GetLastQELCS();// the last quasi-elastic cross section2
896	G4double qelCS = qelCS1+qelCS2; // the last quasi-elastic cross section
897	if(totCS - qelCS < 0.) // only at low energies
898	{
899	totCS = qelCS;
900	totCS1 = qelCS1;
901	totCS2 = qelCS2;
902	}
903	// make different definitions for neutrino and antineutrino
904	G4double mIN=mProt; // Just a prototype (for anu, Z=1, N=0)
905	G4double mOT=mNeut;
906	G4double OT=mlN2;
907	G4double mOT2=mNeut2;
908	G4double mlOT=fmlN;
909	G4double mlsOT=mlsN;
910	if(secnu)
911	{
912	if(am*G4UniformRand()>Z) // Neutron target
913	{
914	targPDG-=1; // subtract neutron
915	projPDG=2112; // neutron is going out
916	mIN =mNeut;
917	OT =mNeut2;
918	mOT2=mNeut2;
919	mlOT=0.;
920	mlsOT=mNeut2;
921	}
922	else
923	{
924	targPDG-=1000; // subtract neutron
925	projPDG=2212; // neutron is going out
926	mOT =mProt;
927	OT =mProt2;
928	mOT2 =mProt2;
929	mlOT =0.;
930	mlsOT=mProt2;
931	}
932	ml=0.;
933	ml2=0.;
934	mldM=0.;
935	mlD2=mPPi2;
936	G4QPDGCode targQPDG(targPDG);
937	G4double rM=targQPDG.GetMass();
938	mIN=tM-rM; // bounded in-mass of the neutron
939	tM=rM;
940	}
941	else if(nuanu)
942	{
943	targPDG-=1; // Neutrino -> subtract neutron
944	G4QPDGCode targQPDG(targPDG);
945	G4double rM=targQPDG.GetMass();
946	mIN=tM-rM; // bounded in-mass of the neutron
947	tM=rM;
948	mOT=mProt;
949	OT=mlP2;
950	mOT2=mProt2;
951	mlOT=fmlP;
952	mlsOT=mlsP;
953	projPDG=2212; // proton is going out
954	}
955	else
956	{
957	if(Z>1\|\|N>0) // Calculate the splitted mass
958	{
959	targPDG-=1000; // Anti-Neutrino -> subtract proton
960	G4QPDGCode targQPDG(targPDG);
961	G4double rM=targQPDG.GetMass();
962	mIN=tM-rM; // bounded in-mass of the proton
963	tM=rM;
964	}
965	else
966	{
967	targPDG=0;
968	mIN=tM;
969	tM=0.;
970	}
971	projPDG=2112; // neutron is going out
972	}
973	G4double s=mIN*(mIN+dKinE); // s=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu)
974	#ifdef debug
975	G4cout<<"G4QCollision::PostStDoIt: s="<<s<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
976	#endif
977	if(s<=OT) // * Do nothing solution *
978	{
979	//Do NothingToDo Action insead of the reaction (@@ Can we make it common?)
980	G4cout<<"G4QCollision::PostStepDoIt: tooSmallFinalMassOfCompound: DoNothing"<<G4endl;
981	aParticleChange.ProposeEnergy(kinEnergy);
982	aParticleChange.ProposeLocalEnergyDeposit(0.);
983	aParticleChange.ProposeMomentumDirection(dir);
984	aParticleChange.ProposeTrackStatus(fAlive);
985	return G4VDiscreteProcess::PostStepDoIt(track,step);
986	}
987	#ifdef debug
988	G4cout<<"G4QCollision::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
989	#endif
990	aParticleChange.ProposeEnergy(0.);
991	aParticleChange.ProposeTrackStatus(fStopAndKill); // the initial neutrino is killed
992	// There is no way back from here
993
994	if ( ((secnu \|\| !nuanu \|\| N) && totCS*G4UniformRand() < qelCS) \|\| s < mlD2 )
995	{ // Quasi-Elastic interaction
996	G4double Q2=0.; // Simulate transferred momentum, in MeV^2
997	if(secnu) Q2=CSmanager2->GetQEL_ExchangeQ2();
998	else Q2=CSmanager->GetQEL_ExchangeQ2();
999	#ifdef debug
1000	G4cout<<"G4QCollision::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s<<",Q2="<<Q2<<G4endl;
1001	#endif
1002	//G4double ds=s+s; // doubled s
1003	G4double sqs=std::sqrt(s); // M_cm
1004	G4double dsqs=sqs+sqs; // 2*M_cm
1005	G4double pi=(s-mIN*mIN)/dsqs; // initial momentum in CMS (checked MK)
1006	G4double dpi=pi+pi; // doubled initial momentum in CMS
1007	G4double sd=s-mlsOT; // s-ml2-mOT2 (mlsOT=m^2_neut+m^2_lept)
1008	G4double qo2=(sd*sd-mlOT)/dsqs; // squared momentum of secondaries in CMS
1009	G4double qo=std::sqrt(qo2); // momentum of secondaries in CMS
1010	G4double cost=(dpi*std::sqrt(qo2+ml2)-Q2-ml2)/dpi/qo; // cos(theta) in CMS (chck MK)
1011	G4LorentzVector t4M(0.,0.,0.,mIN); // 4mom of the effective target
1012	G4LorentzVector c4M=t4M+proj4M; // 4mom of the compound system
1013	t4M.setT(mOT); // now it is 4mom of the outgoing nucleon
1014	scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scattered muon
1015	if(!G4QHadron(c4M).RelDecayIn2(scat4M, t4M, proj4M, cost, cost))
1016	{
1017	G4cerr<<"G4QCol::PSD:c4M="<<c4M<<sqs<<",mM="<<ml<<",tM="<<mOT<<",c="<<cost<<G4endl;
1018	throw G4QException("G4QCollision::HadronizeQuasm: Can't dec QE nu,lept Compound");
1019	}
1020	proj4M=t4M; // 4mom of the new projectile nucleon
1021	}
1022	else // ***** Non Quasi Elastic interaction
1023	{
1024
1025	if ( (secnu && projPDG == 2212) \|\| (!secnu && projPDG == 2112) ) {
1026	targPDG+=1; // Recover target PDG,
1027	} else if ( (secnu && projPDG == 2112) \|\| (!secnu && projPDG == 2212) ) {
1028	targPDG+=1000; // if not quasiEl
1029	}
1030	G4double Q2=0; // Simulate transferred momentum, in MeV^2
1031	if(secnu) Q2=CSmanager->GetNQE_ExchangeQ2();
1032	else Q2=CSmanager2->GetNQE_ExchangeQ2();
1033	#ifdef debug
1034	G4cout<<"G4QColl::PStDoIt: MultiPeriferal s="<<s<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
1035	#endif
1036	if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma
1037	else projPDG=CSmanager->GetExchangePDGCode(); // PDG Code of the effective pion
1038	//@@ Temporary made only for direct interaction and for N=3 (good for small Q2)
1039	//@@ inFuture use N=GetNPartons and directFraction=GetDirectPart, @@ W2...
1040	G4double r=G4UniformRand();
1041	G4double r1=0.5; // (1-x)
1042	if(r<0.5) r1=std::sqrt(r+r)(.5+.1579(r-.5));
1043	else if(r>0.5) r1=1.-std::sqrt(2.-r-r)(.5+.1579(.5-r));
1044	G4double xn=1.-mldM/Momentum; // Normalization of (1-x) [x>mldM/Mom]
1045	G4double x1=xn*r1; // (1-x)
1046	G4double x=1.-x1; // x=2k/M
1047	//G4double W2=(hdM2+Q2/x)*x1; // W2 candidate
1048	G4double mx=hdM*x; // Part of the target to interact with
1049	G4double we=Q2/(mx+mx); // transfered energy
1050	if(we>=kinEnergy-ml-.001) we=kinEnergy-ml-.0001; // safety to avoid nan=sqrt(neg)
1051	G4double pot=kinEnergy-we; // energy of the secondary lepton
1052	G4double mlQ2=ml2+Q2;
1053	G4double cost=(pot-mlQ2/dKinE)/std::sqrt(pot*pot-ml2); // LS cos(theta)
1054	if(std::fabs(cost)>1)
1055	{
1056	#ifdef debug
1057	G4cout<<"*G4QCollision::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
1058	<<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl;
1059	#endif
1060	if(cost>1.) cost=1.;
1061	else cost=-1.;
1062	pot=mlQ2/dKinE+dKinE*ml2/mlQ2; // extreme output momentum
1063	}
1064	G4double lEn=std::sqrt(pot*pot+ml2); // Lepton energy
1065	G4double lPl=pot*cost; // Lepton longitudinal momentum
1066	G4double lPt=potstd::sqrt(1.-costcost); // Lepton transverse momentum
1067	std::pair<G4double,G4double> d2d=Random2DDirection(); // Randomize phi
1068	G4double lPx=lPt*d2d.first;
1069	G4double lPy=lPt*d2d.second;
1070	G4ThreeVector vdir=proj4M.vect(); // 3D momentum of the projectile
1071	G4ThreeVector vz= vdir.unit(); // Ort in the direction of the projectile
1072	G4ThreeVector vv= vz.orthogonal(); // Not normed orthogonal vector (!)
1073	G4ThreeVector vx= vv.unit(); // First ort orthogonal to the direction
1074	G4ThreeVector vy= vz.cross(vx); // Second ort orthoganal to the direction
1075	G4ThreeVector lP= lPlvz+lPxvx+lPy*vy; // 3D momentum of the scattered lepton
1076	scat4M=G4LorentzVector(lP,lEn); // 4mom of the scattered lepton
1077	proj4M-=scat4M; // 4mom of the W/Z (effective pion/gamma)
1078	#ifdef debug
1079	G4cout<<"G4QCollision::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
1080	#endif
1081	// Check that the en/mom transfer is possible, if not -> elastic
1082	G4int fintPDG=targPDG; // Prototype for the compound nucleus
1083	if(!secnu)
1084	{
1085	if(projPDG<0) fintPDG-= 999;
1086	else fintPDG+= 999;
1087	}
1088	G4double fM=G4QPDGCode(fintPDG).GetMass();// compound nucleus Mass (MeV)
1089	G4double fM2=fM*fM;
1090	G4LorentzVector tg4M=G4LorentzVector(0.,0.,0.,tgM);
1091	G4LorentzVector c4M=tg4M+proj4M;
1092	#ifdef debug
1093	G4cout<<"G4QCol::PSDI:fM2="<<fM2<<" <? mc4M="<<c4M.m2()<<",dM="<<fM-tgM<<G4endl;
1094	#endif
1095	if(fM2>=c4M.m2()) // Elastic scattering should be done
1096	{
1097	G4LorentzVector tot4M=tg4M+proj4M+scat4M; // recover the total 4-momentum
1098	s=tot4M.m2();
1099	G4double fs=s-fM2-ml2;
1100	G4double fMl=fM2*ml2;
1101	G4double hQ2max=(fs*fs/2-fMl-fMl)/s; // Maximum possible Q2/2
1102	G4double cost=1.-Q2/hQ2max; // cos(theta) in CMS (use MultProd Q2)
1103	#ifdef debug
1104	G4cout<<"G4QC::PSDI:ct="<<cost<<",Q2="<<Q2<<",hQ2="<<hQ2max<<",4M="<<tot4M<<G4endl;
1105	#endif
1106	G4double acost=std::fabs(cost);
1107	if(acost>1.)
1108	{
1109	if(acost>1.001) G4cout<<"-Warning-G4QCollision::PostStDoIt: cost="<<cost<<G4endl;
1110	if (cost> 1.) cost= 1.;
1111	else if(cost<-1.) cost=-1.;
1112	}
1113	G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,fM); // 4mom of the recoilNucleus
1114	scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scatteredLepton
1115	G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-ml)*.01);
1116	if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
1117	{
1118	G4cerr<<"G4QC::PSDI:t4M="<<tot4M<<",lM="<<ml<<",rM="<<fM<<",cost="<<cost<<G4endl;
1119	//G4Exception("G4QCollision::PostStepDoIt:","027",FatalException,"ElasticDecay");
1120	}
1121	#ifdef debug
1122	G4cout<<"G4QCol::PStDoI:l4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
1123	#endif
1124	// ----------------------------------------------------
1125	G4ParticleDefinition* theDefinition=0; // Prototype of a particle for E-Secondaries
1126	// Fill scattered lepton
1127	if (scatPDG==-11) theDefinition = G4Positron::Positron();
1128	else if(scatPDG== 11) theDefinition = G4Electron::Electron();
1129	else if(scatPDG== 13) theDefinition = G4MuonMinus::MuonMinus();
1130	else if(scatPDG==-13) theDefinition = G4MuonPlus::MuonPlus();
1131	//else if(scatPDG== 15) theDefinition = G4TauMinus::TauMinus();
1132	//else if(scatPDG==-15) theDefinition = G4TauPlus::TauPlus();
1133	if (scatPDG==-12) theDefinition = G4AntiNeutrinoE::AntiNeutrinoE();
1134	else if(scatPDG== 12) theDefinition = G4NeutrinoE::NeutrinoE();
1135	else if(scatPDG== 14) theDefinition = G4NeutrinoMu::NeutrinoMu();
1136	else if(scatPDG==-14) theDefinition = G4AntiNeutrinoMu::AntiNeutrinoMu();
1137	//else if(scatPDG== 16) theDefinition = G4NeutrinoTau::NeutrinoTau();
1138	//else if(scatPDG==-16) theDefinition = G4AntiNeutrinoTau::AntiNeutrinoTau();
1139	else G4cout<<"-Warning-G4QCollision::PostStDoIt: UnknownLepton="<<scatPDG<<G4endl;
1140	G4DynamicParticle* theScL = new G4DynamicParticle(theDefinition,scat4M);
1141	G4Track* scatLep = new G4Track(theScL, localtime, position ); // scattered
1142	scatLep->SetWeight(weight); // weighted
1143	scatLep->SetTouchableHandle(trTouchable); // residual
1144	aParticleChange.AddSecondary(scatLep); // lepton
1145	// Fill residual nucleus
1146	if (fintPDG==90000001) theDefinition = G4Neutron::Neutron(); // neutron
1147	else if(fintPDG==90001000) theDefinition = G4Proton::Proton(); // proton
1148	else // ion
1149	{
1150	G4int fm=static_cast<G4int>(fintPDG/1000000); // Strange part
1151	G4int ZN=fintPDG-1000000*fm;
1152	G4int rZ=static_cast<G4int>(ZN/1000);
1153	G4int rA=ZN-999*rZ;
1154	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
1155	}
1156	G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,reco4M);
1157	G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1158	scatReN->SetWeight(weight); // weighted
1159	scatReN->SetTouchableHandle(trTouchable); // residual
1160	aParticleChange.AddSecondary(scatReN); // nucleus
1161	return G4VDiscreteProcess::PostStepDoIt(track, step);
1162	}
1163	}
1164
1165	//
1166	// quasi-elastic for p+A(Z,N)
1167	//
1168	} else if (aProjPDG == 2212 && Z > 0 && N > 0) {
1169	//else if(2>3)
1170
1171	G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer();
1172	std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, aProjPDG, Z, N);
1173	G4double qepart=fief.first*fief.second;
1174	#ifdef qedebug
1175	G4cout<<"G4QCol::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<",)="<<qepart<<G4endl;
1176	#endif
1177	if(G4UniformRand()<qepart) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim
1178	{
1179	// First decay a nucleus in a nucleon and a residual (A-1) nucleus
1180	G4double dmom=91.; // Fermi momentum (proto default for a deuteron)
1181	if(Z>1\|\|N>1) dmom=286.2*std::pow(-std::log(G4UniformRand()),third);// p_max=250 MeV/c
1182
1183	// Calculate cluster probabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters)
1184	const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
1185	G4double clProb[lCl]={0.6,0.7,0.8}; // N/P,D,t/He3,Alpha, integrated prob for .6,.1,.1,.2
1186	G4double base=1.; // Base for randomization (can be reduced by totZ & totN)
1187	G4int max=lCl; // Number of boundaries (can be reduced by totZ & totN)
1188
1189	// Take into account that at least one nucleon must be left !
1190	// Change max-- to --max - DHW 05/08
1191	G4int A = Z + N; // Baryon number of the nucleus
1192	if (Z<2 \|\| N<2 \|\| A<5) base = clProb[--max]; // Alpha cluster is impossible
1193	if ( (Z > 1 && N < 2) \|\| (Z < 2 && N > 1) )
1194	base=(clProb[max]+clProb[max-1])/2; // t or He3 is impossible
1195
1196	if ( (Z < 2 && N < 2) \|\| A < 4) base=clProb[--max]; // Both He3 and t clusters are impossible
1197
1198	if(A<3) base=clProb[--max]; // Deuteron cluster is impossible
1199	G4int cln=0; // Cluster#0 (Default for the selected nucleon)
1200	if(max) // Not only nucleons are possible
1201	//if(2>3)
1202	{
1203	G4double ran=base*G4UniformRand(); // Base can be reduced
1204	G4int ic=0; // Start from the smallest cluster boundary
1205	while(ic<max) if(ran>clProb[ic++]) cln=ic;
1206	}
1207	G4ParticleDefinition* theDefinition; // Prototype for qfNucleon
1208	G4bool cp1 = cln+2==A; // A=ClusterBN+1 condition
1209	// Values to be defined in the following IF/ELSE
1210	G4LorentzVector r4M(0.,0.,0.,0.); // Prototype of 4mom of the residual nucleus
1211	G4LorentzVector n4M(0.,0.,0.,0.); // Prototype of 4mom of the quasi-cluster
1212	G4int nPDG=90000001; // Prototype for quasi-cluster mass calculation
1213	G4int restPDG=targPDG; // Prototype should be reduced by quasi-cluster
1214	G4int rA=Z+N-1; // Prototype for the residualNucl definition
1215	G4int rZ=Z; // residZ: OK for the quasi-free neutron
1216	G4int nA=1; // Prototype for the quasi-cluster definition
1217	G4int nZ=0; // nA=1,nZ=0: OK for the quasi-free neutron
1218	G4double qM=mNeut; // Free mass of the quasi-free cluster
1219	if(!cln \|\| cp1) // Split in nucleon + (A-1) with Fermi momentum
1220	{
1221	G4int nln=0;
1222	if(cln==2) nln=1; // @@ only for cp1: t/He3 choice from A=4
1223	// mass(A)=tM. Calculate masses of A-1 (rM) and mN (mNeut or mProt bounded mass)
1224	if ( ((!cln \|\| cln == 2) && G4UniformRand()*(A-cln) > (N-nln)) \|\|
1225	((cln == 3 \|\| cln == 1) && Z > N) )
1226	{
1227	nPDG=90001000; // Update quasi-free nucleon PDGCode to P
1228	nZ=1; // Change charge of the quasiFree nucleon
1229	qM=mProt; // Update quasi-free nucleon mass
1230	rZ--; // Reduce the residual Z
1231	restPDG-=1000; // Reduce the residual PDGCode
1232	}
1233	else restPDG--;
1234	G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1235	G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1236	r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1237	G4double rM2=rM*rM;
1238	G4double nM=std::sqrt(rM2+tMtM-(tM+tM)std::sqrt(rM2+dmom*dmom));// M of q-nucleon
1239	n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon
1240	#ifdef qedebug
1241	G4cout<<"G4QCollis::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1242	#endif
1243	if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1244	{
1245	G4cerr<<"G4QCol::PostStDoIt: M="<<tM<<"<rM="<<rM<<"+nM="<<nM<<"="<<rM+nM<<G4endl;
1246	throw G4QException("G4QCollision::HadronizeQuasm:Can'tDec totNuc->QENuc+ResNuc");
1247	}
1248	#ifdef qedebug
1249	G4cout<<"G4QCol::PStDoIt:QE-N,RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
1250	#endif
1251	if(cp1 && cln) // Quasi-cluster case: swap the output
1252	{
1253	qM=rM; // Scattering will be made on a cluster
1254	nln=nPDG;
1255	nPDG=restPDG;
1256	restPDG=nln;
1257	t4M=n4M;
1258	n4M=r4M;
1259	r4M=t4M;
1260	nln=nZ;
1261	nZ=rZ;
1262	rZ=nln;
1263	nln=nA;
1264	nA=rA;
1265	rA=nln;
1266	}
1267	}
1268	else // Split a cluster (w or w/o "Fermi motion" and "Fermi decay")
1269	{
1270	if(cln==1)
1271	{
1272	nPDG=90001001; // Deuteron
1273	qM=mDeut;
1274	nA=2;
1275	nZ=1;
1276	restPDG-=1001;
1277	}
1278	else if(cln==2)
1279	{
1280	nA=3;
1281	if(G4UniformRand()*(A-2)>(N-1)) // He3
1282	{
1283	nPDG=90002001;
1284	qM=mHel3;
1285	nZ=2;
1286	restPDG-=2001;
1287	}
1288	else // tritium
1289	{
1290	nPDG=90001002;
1291	qM=mTrit;
1292	nZ=1;
1293	restPDG-=1002;
1294	}
1295	}
1296	else
1297	{
1298	nPDG=90002002; // Alpha
1299	qM=mAlph;
1300	nA=4;
1301	nZ=2;
1302	restPDG-=2002;
1303	}
1304	rA=A-nA;
1305	rZ=Z-nZ;
1306	// This is a simple case of cluster at rest
1307	//G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1308	//r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1309	//n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster
1310	// --- End of the "simple case of cluster at rest"
1311	// Make a fake quasi-Fermi distribution for clusters (clusters are not at rest)
1312	G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1313	G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1314	r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1315	G4double rM2=rM*rM;
1316	G4double nM=std::sqrt(rM2+tMtM-(tM+tM)std::sqrt(rM2+dmom*dmom));// M of q-cluster
1317	n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon
1318	#ifdef qedebug
1319	G4cout<<"G4QCollis::PStDoIt:QEC,p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1320	#endif
1321	if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1322	{
1323	G4cerr<<"G4QCol::PostStDoIt: M="<<tM<<"<rM="<<rM<<"+cM="<<nM<<"="<<rM+nM<<G4endl;
1324	throw G4QException("G4QCollision::HadronizeQuasm:Can'tDec totNuc->QEClu+ResNuc");
1325	}
1326	// --- End of the moving cluster implementation ---
1327	#ifdef qedebug
1328	G4cout<<"G4QCol::PStDoIt:QEC,RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
1329	#endif
1330	}
1331	G4LorentzVector s4M=n4M+proj4M; // Tot 4-momentum for scattering
1332	G4double prjM2 = proj4M.m2();
1333	G4double prjM = std::sqrt(prjM2); // @@ Get from pPDG (?)
1334	G4double minM = prjM+qM; // Min mass sum for the final products
1335	G4double cmM2 =s4M.m2();
1336	if(cmM2>minM*minM)
1337	{
1338	#ifdef qedebug
1339	G4cout<<"G4QCol::PStDoIt:*Enter,cmM2="<<cmM2<<" > minM2="<<minMminM<<G4endl;
1340	#endif
1341	// Estimate and randomize charge-exchange with quasi-free cluster
1342	G4bool chex=false; // Flag of the charge exchange scattering
1343	G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true
1344	//if(cln&&!cp1 &&(projPDG==2212&&rA>rZ \|\| projPDG==2112&&rZ>1))// @@ Use proj chex
1345	if(2>3)
1346	{
1347	#ifdef qedebug
1348	G4cout<<"G4QCol::PStDoIt:-Enter,P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
1349	#endif
1350	G4double tprM=mProt;
1351	G4double tprM2=mProt2;
1352	G4int tprPDG=2212;
1353	G4int tresPDG=restPDG+999;
1354	if(projPDG==2212)
1355	{
1356	projpt=G4Neutron::Neutron();
1357	tprM=mNeut;
1358	tprM2=mNeut2;
1359	tprPDG=2112;
1360	tresPDG=restPDG-999;
1361	}
1362	minM=tprM+qM;
1363	G4double efE=(cmM2-tprM2-qM*qM)/(qM+qM);
1364	G4double efP=std::sqrt(efE*efE-tprM2);
1365	G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!)
1366	#ifdef qedebug
1367	G4cout<<"G4QCol::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
1368	#endif
1369	if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl)) // minM is redefined
1370	{
1371	projPDG=tprPDG;
1372	prjM=tprM;
1373	G4double rM=G4QPDGCode(tresPDG).GetMass();// Mass of the residual nucleus
1374	r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1375	n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster
1376	chex=true; // Confirm charge exchange scattering
1377	}
1378	}
1379	//
1380	std::pair<G4LorentzVector,G4LorentzVector> sctout=qfMan->Scatter(nPDG, n4M,
1381	projPDG, proj4M);
1382	#ifdef qedebug
1383	G4cout<<"G4QCollis::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
1384	<<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl;
1385	#endif
1386	aParticleChange.ProposeLocalEnergyDeposit(0.); // Everything is in particles
1387	// @@ @@ @@ Coulomb barriers must be checked !! @@ @@ @@ Skip if not
1388	if(chex) // ==> Projectile is changed: fill everything to secondaries
1389	{
1390	aParticleChange.ProposeEnergy(0.); // @@ ??
1391	aParticleChange.ProposeTrackStatus(fStopAndKill); // projectile nucleon is killed
1392	aParticleChange.SetNumberOfSecondaries(3);
1393	G4DynamicParticle* thePrH = new G4DynamicParticle(projpt,sctout.second);
1394	G4Track* scatPrH = new G4Track(thePrH, localtime, position ); // scattered & chex
1395	scatPrH->SetWeight(weight); // weighted
1396	scatPrH->SetTouchableHandle(trTouchable); // projectile
1397	aParticleChange.AddSecondary(scatPrH); // hadron
1398	}
1399	else // ==> The leading particle is filled to the updated projectilee
1400	{
1401	aParticleChange.SetNumberOfSecondaries(2); // @@ if proj=leading
1402	G4double ldT=(sctout.second).e()-prjM; // kin Energy of scat project.
1403	aParticleChange.ProposeEnergy(ldT); // Change the kin Energy
1404	G4ThreeVector ldV=(sctout.second).vect(); // Change momentum direction
1405	aParticleChange.ProposeMomentumDirection(ldV/ldV.mag());
1406	aParticleChange.ProposeTrackStatus(fAlive);
1407	}
1408	// ---------------------------------------------------------
1409	// Fill scattered quasi-free nucleon
1410	if (nPDG==90000001) theDefinition = G4Neutron::Neutron();
1411	else if(nPDG==90001000) theDefinition = G4Proton::Proton();
1412	else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion
1413	G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,sctout.first);
1414	G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // scattered
1415	scatQFN->SetWeight(weight); // weighted
1416	scatQFN->SetTouchableHandle(trTouchable); // quasi-free
1417	aParticleChange.AddSecondary(scatQFN); // nucleon/cluster
1418	// ----------------------------------------------------
1419	// Fill residual nucleus
1420	if (restPDG==90000001) theDefinition = G4Neutron::Neutron();
1421	else if(restPDG==90001000) theDefinition = G4Proton::Proton();
1422	else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion
1423	G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1424	G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1425	scatReN->SetWeight(weight); // weighted
1426	scatReN->SetTouchableHandle(trTouchable); // residual
1427	aParticleChange.AddSecondary(scatReN); // nucleus
1428	return G4VDiscreteProcess::PostStepDoIt(track, step);
1429	}
1430	#ifdef qedebug
1431	else G4cout<<"G4QCol::PSD: OUT, M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl;
1432	#endif
1433	}
1434	} // end lepto-nuclear, neutrino-nuclear, proton quasi-elastic
1435
1436	EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
1437	if(absMom) EnMomConservation+=lead4M; // Add E/M of leading System
1438	#ifdef debug
1439	G4cout<<"G4QCollision::PostStDoIt:before St="<<aParticleChange.GetTrackStatus()<<G4endl;
1440	#endif
1441
1442	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin
1443	//G4bool elF=false; // Flag of the ellastic scattering is "false" by default
1444	//G4double eWei=1.;
1445	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- End
1446	#ifdef debug
1447	G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
1448	#endif
1449	G4QHadron* pH = new G4QHadron(projPDG,proj4M); // ---> DELETED -->-- -+
1450	//if(momentum<1000.) // Condition for using G4QEnvironment (not G4QuasmonString) \|
1451	{ // \|
1452	G4QHadronVector projHV; // \|
1453	projHV.push_back(pH); // DESTROYED over 2 lines-+ \|
1454	G4QEnvironment* pan= new G4QEnvironment(projHV,targPDG);// ---> DELETED --->----+ \| \|
1455	std::for_each(projHV.begin(), projHV.end(), DeleteQHadron()); // <---<------<---+-+-+
1456	projHV.clear(); // <------------<---------------<-------------------<-----------+-+ .
1457	#ifdef debug
1458	G4cout<<"G4QCol::PStDoIt:Proj="<<projPDG<<proj4M<<",Targ="<<targPDG<<G4endl; // \| .
1459	#endif
1460	try // \| .
1461	{ // \| .
1462	delete output; // \| .
1463	output = pan->Fragment();// DESTROYED in the end of the LOOP work space \| .
1464	} // \| .
1465	catch (G4QException& error)// \| .
1466	{ // \| .
1467	//#ifdef pdebug
1468	G4cerr<<"***G4QCollision::PostStepDoIt: G4QE Exception is catched"<<G4endl;// \| .
1469	//#endif
1470	G4Exception("G4QCollision::PostStepDoIt:","27",FatalException,"CHIPSCrash");//\| .
1471	} // \| .
1472	delete pan; // Delete the Nuclear Environment <-<--+ .
1473	} // .
1474	//else // Use G4QuasmonString .
1475	//{ // ^
1476	// G4QuasmonString* pan= new G4QuasmonString(pH,false,targPDG,false);//-> DELETED --+ \|
1477	// delete pH; // --------<-------+--+
1478	#ifdef debug
1479	// G4double mp=G4QPDGCode(projPDG).GetMass(); // Mass of the projectile particle \|
1480	// G4cout<<"G4QCollision::PostStepDoIt: pPDG="<<projPDG<<", pM="<<mp<<G4endl; // \|
1481	#endif
1482	// //G4int tNH=0; // Prototype of the number of secondaries inOut\|
1483	// try // \|
1484	// { // \|
1485	// delete output; // \|
1486	// output = pan->Fragment();// DESTROYED in the end of the LOOP work space \|
1487	// // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin \|
1488	// //tNH=pan->GetNOfHadrons(); // For the test purposes of the String \|
1489	// //if(tNH==2) // At least 2 hadrons are in the Constr.Output \|
1490	// //{// \|
1491	// // elF=true; // Just put a flag for the ellastic Scattering \|
1492	// // delete output; // Delete a prototype of dummy G4QHadronVector \|
1493	// // output = pan->GetHadrons(); // DESTROYED in the end of the LOOP work space \|
1494	// //}// \|
1495	// //eWei=pan->GetWeight(); // Just an example for the weight of the event \|
1496	#ifdef debug
1497	// //G4cout<<"=====>>G4QCollision::PostStepDoIt: elF="<<elF<<",n="<<tNH<<G4endl;//\|
1498	#endif
1499	// // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End \|
1500	// } // \|
1501	// catch (G4QException& error)// \|
1502	// { // \|
1503	// //#ifdef pdebug
1504	// G4cerr<<"***G4QCollision::PostStepDoIt: GEN Exception is catched"<<G4endl; // \|
1505	// //#endif
1506	// G4Exception("G4QCollision::PostStDoIt:","27",FatalException,"QString Excep");//\|
1507	// } // \|
1508	// delete pan; // Delete the Nuclear Environment ---<--+
1509	//}
1510	// --- the scattered hadron with changed nature can be added here ---
1511	if(scat)
1512	{
1513	G4QHadron* scatHadron = new G4QHadron(scatPDG,scat4M);
1514	output->push_back(scatHadron);
1515	}
1516	G4int qNH=leadhs->size();
1517	if(absMom)
1518	{
1519	if(qNH) for(G4int iq=0; iq<qNH; iq++)
1520	{
1521	G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron
1522	output->push_back(loh);
1523	}
1524	delete leadhs;
1525	}
1526
1527
1528	// ------------- From here the secondaries are filled -------------------------
1529	G4int tNH = output->size(); // A#of hadrons in the output
1530	aParticleChange.SetNumberOfSecondaries(tNH);
1531	// Now add nuclear fragments
1532	#ifdef debug
1533	G4cout<<"G4QCollision::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
1534	#endif
1535	#ifdef ppdebug
1536	if(absMom)G4cout<<"G4QCollision::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
1537	#endif
1538	G4int nOut=output->size(); // Real length of the output @@ Temporary
1539	if(tNH==1 && !scat) // @@ Temporary. Find out why it happened!
1540	{
1541	G4cout<<"-Warning-G4QCollision::PostStepDoIt: 1 secondary! absMom="<<absMom;
1542	if(absMom) G4cout<<", qNH="<<qNH;
1543	G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
1544	G4cout<<G4endl;
1545	tNH=0;
1546	delete output->operator[](0); // delete the creazy hadron
1547	output->pop_back(); // clean up the output vector
1548	}
1549	if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QCollision::PostStepDoIt: 2 # "<<nOut<<G4endl;
1550	// Deal with ParticleChange final state interface to GEANT4 output of the process
1551	//if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
1552	if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
1553	{
1554	// Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
1555	// Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
1556	// The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
1557	G4QHadron* hadr=(*output)[i]; // Pointer to the output hadron
1558	G4int PDGCode = hadr->GetPDGCode();
1559	G4int nFrag = hadr->GetNFragments();
1560	#ifdef pdebug
1561	G4cout<<"G4QCollision::PostStepDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag
1562	<<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
1563	#endif
1564	if(nFrag) // Skip intermediate (decayed) hadrons
1565	{
1566	#ifdef debug
1567	G4cout<<"G4QCollision::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
1568	#endif
1569	delete hadr;
1570	continue;
1571	}
1572	G4DynamicParticle* theSec = new G4DynamicParticle;
1573	G4ParticleDefinition* theDefinition;
1574	if (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
1575	else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
1576	else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
1577	else if(PDGCode==311 \|\| PDGCode==-311)
1578	{
1579	if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L
1580	else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
1581	}
1582	else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
1583	else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
1584	else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
1585	else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
1586	else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
1587	else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
1588	{
1589	G4int aZ = hadr->GetCharge();
1590	G4int aA = hadr->GetBaryonNumber();
1591	#ifdef pdebug
1592	G4cout<<"G4QCollision::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
1593	#endif
1594	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
1595	}
1596	//else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
1597	else
1598	{
1599	#ifdef pdebug
1600	G4cout<<"G4QCollision::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
1601	#endif
1602	theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
1603	#ifdef pdebug
1604	G4cout<<"G4QCollision::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
1605	#endif
1606	}
1607	if(!theDefinition)
1608	{
1609	#ifdef debug
1610	G4cout<<"---Warning---G4QCollision::PostStepDoIt: drop PDG="<<PDGCode<<G4endl;
1611	#endif
1612	delete hadr;
1613	continue;
1614	}
1615	#ifdef pdebug
1616	G4cout<<"G4QCollision::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
1617	#endif
1618	theSec->SetDefinition(theDefinition);
1619	G4LorentzVector h4M=hadr->Get4Momentum();
1620	EnMomConservation-=h4M;
1621	#ifdef tdebug
1622	G4cout<<"G4QCollis::PSDI:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
1623	#endif
1624	#ifdef debug
1625	G4cout<<"G4QCollision::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
1626	#endif
1627	theSec->Set4Momentum(h4M); // ^
1628	delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
1629	#ifdef debug
1630	G4ThreeVector curD=theSec->GetMomentumDirection(); // ^
1631	G4double curM=theSec->GetMass(); // \|
1632	G4double curE=theSec->GetKineticEnergy()+curM; // ^
1633	G4cout<<"G4QCollis::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;// \|
1634	#endif
1635	G4Track* aNewTrack = new G4Track(theSec, localtime, position ); // ^
1636	aNewTrack->SetWeight(weight); // weighted \|
1637	aNewTrack->SetTouchableHandle(trTouchable); // \|
1638	aParticleChange.AddSecondary( aNewTrack ); // \|
1639	#ifdef debug
1640	G4cout<<"G4QCollision::PostStepDoIt:#"<<i<<" is done."<<G4endl; // \|
1641	#endif
1642	} // \|
1643	delete output; // instances of the G4QHadrons from the output are already deleted above +
1644	#ifdef debug
1645	G4cout<<"G4QCollision::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
1646	#endif
1647	if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15)
1648	aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the absorbed particle
1649	//return &aParticleChange; // This is not enough (ClearILL)
1650	#ifdef pdebug
1651	G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()
1652	<<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl;
1653	#endif
1654	#ifdef debug
1655	G4cout<<"G4QCollision::PostStepDoIt:* PostStepDoIt is done *, P="<<aProjPDG<<", St="
1656	<<aParticleChange.GetTrackStatus()<<G4endl;
1657	#endif
1658	return G4VDiscreteProcess::PostStepDoIt(track, step);
1659	}
1660
1661	std::pair<G4double,G4double> G4QCollision::Random2DDirection()
1662	{
1663	G4double sp=0; // sin(phi)
1664	G4double cp=1.; // cos(phi)
1665	G4double r2=2.; // to enter the loop
1666	while(r2>1. \|\| r2<.0001) // pi/4 efficiency
1667	{
1668	G4double s=G4UniformRand();
1669	G4double c=G4UniformRand();
1670	sp=1.-s-s;
1671	cp=1.-c-c;
1672	r2=spsp+cpcp;
1673	}
1674	G4double norm=std::sqrt(r2);
1675	return std::make_pair(sp/norm,cp/norm);
1676	}

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

Download in other formats: