Context Navigation

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

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

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

Download in other formats: