Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QAtomicElectronScattering.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago
update processes
File size: 33.2 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: G4QAtomicElectronScattering.cc,v 1.4 2008/10/02 21:10:07 dennis Exp $
27	// GEANT4 tag $Name: geant4-09-02-ref-02 $
28	//
29	// ---------------- G4QAtomicElectronScattering class -----------------
30	// by Mikhail Kossov, December 2003.
31	// G4QAtomicElectronScattering 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
40	#include "G4QAtomicElectronScattering.hh"
41
42	G4QAtomicElectronScattering::G4QAtomicElectronScattering(const G4String& processName):
43	G4VDiscreteProcess(processName, fElectromagnetic)
44	{
45	#ifdef debug
46	G4cout<<"G4QAtomicElectronScattering::Constructor is called"<<G4endl;
47	#endif
48	if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
49
50	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
51	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
52	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
53	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
54	//@@ Initialize here the G4QuasmonString parameters
55	}
56
57	G4bool G4QAtomicElectronScattering::manualFlag=false; // If false:use standard parameters
58	G4double G4QAtomicElectronScattering::Temperature=180.; // Critical Temperature (High Ener)
59	G4double G4QAtomicElectronScattering::SSin2Gluons=0.3; // Supression of s-quarks (to u&d)
60	G4double G4QAtomicElectronScattering::EtaEtaprime=0.3; // Supression of eta(gg->qq/3g->qq)
61	G4double G4QAtomicElectronScattering::freeNuc=0.5; // % of free nucleons on a surface
62	G4double G4QAtomicElectronScattering::freeDib=0.05; // % of free diBaryons on a surface
63	G4double G4QAtomicElectronScattering::clustProb=5.; // Nuclear clusterization parameter
64	G4double G4QAtomicElectronScattering::mediRatio=10.; // medium/vacuum hadronizationRatio
65	G4int G4QAtomicElectronScattering::nPartCWorld=152; // #of particles in the CHIPS World
66	G4double G4QAtomicElectronScattering::SolidAngle=0.5; // A part of Solid Angle to capture
67	G4bool G4QAtomicElectronScattering::EnergyFlux=false; // Flag to use EnergyFlux or MultyQ
68	G4double G4QAtomicElectronScattering::PiPrThresh=141.4; // PiProductionThreshold for gammas
69	G4double G4QAtomicElectronScattering::M2ShiftVir=20000.;// M2=-Q2=m_pi^2 shift of virtGamma
70	G4double G4QAtomicElectronScattering::DiNuclMass=1880.; // Double Nucleon Mass for VirtNorm
71
72	void G4QAtomicElectronScattering::SetManual() {manualFlag=true;}
73	void G4QAtomicElectronScattering::SetStandard() {manualFlag=false;}
74
75	// Fill the private parameters
76	void G4QAtomicElectronScattering::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
77	G4double fN, G4double fD, G4double cP, G4double mR,
78	G4int nParCW, G4double solAn, G4bool efFlag,
79	G4double piThresh, G4double mpisq, G4double dinum)
80	{// =============================================================================
81	Temperature=temper;
82	SSin2Gluons=ssin2g;
83	EtaEtaprime=etaetap;
84	freeNuc=fN;
85	freeDib=fD;
86	clustProb=cP;
87	mediRatio=mR;
88	nPartCWorld = nParCW;
89	EnergyFlux=efFlag;
90	SolidAngle=solAn;
91	PiPrThresh=piThresh;
92	M2ShiftVir=mpisq;
93	DiNuclMass=dinum;
94	G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
95	G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
96	G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
97	G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
98	}
99
100	// Destructor
101
102	G4QAtomicElectronScattering::~G4QAtomicElectronScattering() {}
103
104
105	G4LorentzVector G4QAtomicElectronScattering::GetEnegryMomentumConservation()
106	{
107	return EnMomConservation;
108	}
109
110	G4int G4QAtomicElectronScattering::GetNumberOfNeutronsInTarget()
111	{
112	return nOfNeutrons;
113	}
114
115	G4double G4QAtomicElectronScattering::GetMeanFreePath(const G4Track& aTrack,
116	G4double,G4ForceCondition* Fc)
117	{
118	*Fc = NotForced;
119	const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
120	G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
121	if( !IsApplicable(*incidentParticleDefinition))
122	G4cout<<"-Wa-G4QAtElScat::GetMeanFreePath called for not implemented particle"<<G4endl;
123	// Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
124	G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
125	const G4Material* material = aTrack.GetMaterial(); // Get the current material
126	const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
127	const G4ElementVector* theElementVector = material->GetElementVector();
128	G4int nE=material->GetNumberOfElements();
129	#ifdef debug
130	G4cout<<"G4QAtomElectScattering::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
131	#endif
132	G4bool leptoNuc=false; // By default the reaction is not lepto-nuclear
133	G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
134	if(incidentParticleDefinition == G4Electron::Electron())
135	{
136	CSmanager=G4QElectronNuclearCrossSection::GetPointer();
137	leptoNuc=true;
138	}
139	else G4cout<<"G4QAtomEScattering::GetMeanFreePath:Particle isn't known in CHIPS"<<G4endl;
140
141	G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singelton
142	G4double sigma=0.;
143	for(G4int i=0; i<nE; ++i)
144	{
145	G4int Z = static_cast<G4int>((*theElementVector)[i]->GetZ()); // Z of the Element
146	std::vector<std::pair<G4int,G4double>> cs= Isotopes->GetCSVector(Z); // Pointer to CS
147	G4int nIs=cs->size(); // A#Of Isotopes in the Element
148	if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
149	{
150	std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
151	G4int N=curIs->first; // #ofNeuterons in the isotope
152	curIs->second = CSmanager->GetCrossSection(true,Momentum,Z,N,13); // CS calculation
153	} // End of temporary initialization of the cross sections in the G4QIsotope singeltone
154	sigma+=Isotopes->GetMeanCrossSection(Z)NOfNucPerVolume[i]; // SUM(MeanCSNOFNperV)
155	} // End of LOOP over Elements
156
157	// Check that cross section is not zero and return the mean free path
158	if(sigma > 0.) return 1./sigma; // Mean path [distance]
159	return DBL_MAX;
160	}
161
162
163	G4bool G4QAtomicElectronScattering::IsApplicable(const G4ParticleDefinition& particle)
164	{
165	if (particle == *( G4MuonPlus::MuonPlus() )) return true;
166	else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
167	else if (particle == *( G4TauPlus::TauPlus() )) return true;
168	else if (particle == *( G4TauMinus::TauMinus() )) return true;
169	else if (particle == *( G4Electron::Electron() )) return true;
170	else if (particle == *( G4Positron::Positron() )) return true;
171	else if (particle == *( G4Gamma::Gamma() )) return true;
172	else if (particle == *( G4Proton::Proton() )) return true;
173	//else if (particle == *( G4Neutron::Neutron() )) return true;
174	//else if (particle == *( G4PionMinus::PionMinus() )) return true;
175	//else if (particle == *( G4PionPlus::PionPlus() )) return true;
176	//else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
177	//else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
178	//else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
179	//else if (particle == *(G4KaonZeroShort::KaonZeroShort())) return true;
180	//else if (particle == *( G4Lambda::Lambda() )) return true;
181	//else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
182	//else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
183	//else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
184	//else if (particle == *( G4XiMinus::XiMinus() )) return true;
185	//else if (particle == *( G4XiZero::XiZero() )) return true;
186	//else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
187	//else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
188	//else if (particle == *( G4AntiProton::AntiProton() )) return true;
189	#ifdef debug
190	G4cout<<"***G4QAtomElScattering::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
191	#endif
192	return false;
193	}
194
195	G4VParticleChange* G4QAtomicElectronScattering::PostStepDoIt(const G4Track& track,
196	const G4Step& step)
197	{
198	static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass
199	static const G4double mu2=mu*mu; // squared muon mass
200	//static const G4double dpi=M_PI+M_PI; // 2pi (for Phi distr.) changed to twopi*
201	static const G4double mNeut= G4QPDGCode(2112).GetMass();
202	static const G4double mProt= G4QPDGCode(2212).GetMass();
203	static const G4double dM=mProt+mNeut; // doubled nucleon mass
204	//static const G4double mPi0 = G4QPDGCode(111).GetMass();
205	//static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
206	//static const G4double mPi = G4QPDGCode(211).GetMass();
207	//static const G4double mMu = G4QPDGCode(13).GetMass();
208	//static const G4double mTau = G4QPDGCode(15).GetMass();
209	//static const G4double mEl = G4QPDGCode(11).GetMass();
210	//
211	const G4DynamicParticle* projHadron = track.GetDynamicParticle();
212	const G4ParticleDefinition* particle=projHadron->GetDefinition();
213	G4LorentzVector proj4M=projHadron->Get4Momentum();
214	G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
215	G4double Momentum=proj4M.rho();
216	if(std::fabs(Momentum-momentum)>.001) G4cerr<<"G4QAtElScat::PSDI P="<<Momentum<<"="
217	<<momentum<<G4endl;
218	#ifdef debug
219	G4double mp=proj4M.m();
220	G4cout<<"G4QAtomElScattering::PostStepDoIt called, P="<<Momentum<<"="<<momentum<<G4endl;
221	#endif
222	if (!IsApplicable(*particle)) // Check applicability
223	{
224	G4cerr<<"G4QAtomElectScat::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented"
225	<<G4endl;
226	return 0;
227	}
228	const G4Material* material = track.GetMaterial(); // Get the current material
229	G4int Z=0;
230	const G4ElementVector* theElementVector = material->GetElementVector();
231	G4int i=0;
232	G4double sum=0.;
233	G4int nE=material->GetNumberOfElements();
234	#ifdef debug
235	G4cout<<"G4QAtomElectronScat::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
236	#endif
237	G4int projPDG=0; // PDG Code prototype for the captured hadron
238	// Not all these particles are implemented yet (see Is Applicable)
239	if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
240	else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
241	else if (particle == G4Electron::Electron() ) projPDG= 11;
242	else if (particle == G4Positron::Positron() ) projPDG= -11;
243	else if (particle == G4Gamma::Gamma() ) projPDG= 22;
244	else if (particle == G4Proton::Proton() ) projPDG= 2212;
245	else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
246	else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
247	else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
248	else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 2112;
249	else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
250	else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
251	else if (particle == G4KaonZeroShort::KaonZeroShort()) projPDG= 310;
252	else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
253	else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
254	else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
255	else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
256	else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
257	else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
258	else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
259	else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
260	else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
261	else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
262	else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
263	#ifdef debug
264	G4int prPDG=particle->GetPDGEncoding();
265	G4cout<<"G4QAtomElScat::PostStepRestDoIt: projPDG="<<projPDG<<",stPDG="<<prPDG<<G4endl;
266	#endif
267	if(!projPDG)
268	{
269	G4cerr<<"-Warning-G4QAtomElScattering::PostStepDoIt:Undefined captured hadron"<<G4endl;
270	return 0;
271	}
272	// @@ It's a standard randomization procedure, which can be placed in G4QMaterial class
273	std::vector<G4double> sumfra;
274	for(i=0; i<nE; ++i)
275	{
276	G4double frac=material->GetFractionVector()[i];
277	sum+=frac;
278	sumfra.push_back(sum); // remember the summation steps
279	}
280	G4double rnd = sum*G4UniformRand();
281	for(i=0; i<nE; ++i) if (rnd<sumfra[i]) break;
282	G4Element* pElement=(*theElementVector)[i];
283	Z=static_cast<G4int>(pElement->GetZ());
284	if(Z<=0)
285	{
286	G4cerr<<"-Warning-G4QAtomicElectronScattering::PostStepDoIt: Element's Z="<<Z<<G4endl;
287	if(Z<0) return 0;
288	}
289	G4int N = Z;
290	G4int isoSize=0; // The default for the isoVectorLength is 0
291	G4IsotopeVector* isoVector=pElement->GetIsotopeVector();
292	if(isoVector) isoSize=isoVector->size(); // Get real size of the isotopeVector if exists
293	#ifdef debug
294	G4cout<<"G4QAtomicElectronScattering::PostStepDoIt: isovectorLength="<<isoSize<<G4endl;
295	#endif
296	if(isoSize) // The Element has not trivial abumdance set
297	{
298	// @@ the following solution is temporary till G4Element can contain the QIsotopIndex
299	G4int curInd=G4QIsotope::Get()->GetLastIndex(Z);
300	if(!curInd) // The new artificial element must be defined
301	{
302	std::vector<std::pair<G4int,G4double>> newAbund =
303	new std::vector<std::pair<G4int,G4double>*>;
304	G4double* abuVector=pElement->GetRelativeAbundanceVector();
305	for(G4int j=0; j<isoSize; j++)
306	{
307	N=pElement->GetIsotope(j)->GetN()-Z;
308	if(pElement->GetIsotope(j)->GetZ()!=Z) G4cerr<<"*G4QCaptureAtRest::AtRestDoIt: Z="
309	<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
310	G4double abund=abuVector[j];
311	std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
312	#ifdef debug
313	G4cout<<"G4QAtomElScat::PostStepDoIt:pair#="<<j<<", N="<<N<<",ab="<<abund<<G4endl;
314	#endif
315	newAbund->push_back(pr);
316	}
317	#ifdef debug
318	G4cout<<"G4QAtomElectScat::PostStepDoIt:pairVectorLength="<<newAbund->size()<<G4endl;
319	#endif
320	curInd=G4QIsotope::Get()->InitElement(Z,1,newAbund);
321	for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];
322	delete newAbund;
323	}
324	// @@ ^^^^^^^^^^ End of the temporary solution ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
325	N = G4QIsotope::Get()->GetNeutrons(Z,curInd);
326	}
327	else N = G4QIsotope::Get()->GetNeutrons(Z);
328	nOfNeutrons=N; // Remember it for energy-mom. check
329	G4double dd=0.025;
330	G4double am=Z+N;
331	G4double sr=std::sqrt(am);
332	G4double dsr=0.01*(sr+sr);
333	if(dsr<dd)dsr=dd;
334	if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // ManualP
335	else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP ClustPars
336	else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.);//Pi- ClustPars
337	#ifdef debug
338	G4cout<<"G4QAtomElectScattering::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
339	#endif
340	if(N<0)
341	{
342	G4cerr<<"---Warning---G4QAtomElectScat::PostStepDoIt:Element with N="<<N<< G4endl;
343	return 0;
344	}
345	if(projPDG==11\|\|projPDG==-11\|\|projPDG==13\|\|projPDG==-13\|\|projPDG==15\|\|projPDG==-15)
346	{ // Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + neutrino & QE
347	G4double kinEnergy= projHadron->GetKineticEnergy();
348	G4ParticleMomentum dir = projHadron->GetMomentumDirection();
349	G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
350	G4int aProjPDG=std::abs(projPDG);
351	if(aProjPDG==13) CSmanager=G4QMuonNuclearCrossSection::GetPointer();
352	if(aProjPDG==15) CSmanager=G4QTauNuclearCrossSection::GetPointer();
353	G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,13);//Recalculate CrossSect
354	// @@ check a possibility to separate p, n, or alpha (!)
355	if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
356	{
357	//Do Nothing Action insead of the reaction
358	aParticleChange.ProposeEnergy(kinEnergy);
359	aParticleChange.ProposeLocalEnergyDeposit(0.);
360	aParticleChange.ProposeMomentumDirection(dir) ;
361	return G4VDiscreteProcess::PostStepDoIt(track,step);
362	}
363	G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
364	if( kinEnergy < photonEnergy )
365	{
366	//Do Nothing Action insead of the reaction
367	G4cerr<<"G4QAtomElectScat::PSDoIt: phE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl;
368	aParticleChange.ProposeEnergy(kinEnergy);
369	aParticleChange.ProposeLocalEnergyDeposit(0.);
370	aParticleChange.ProposeMomentumDirection(dir) ;
371	return G4VDiscreteProcess::PostStepDoIt(track,step);
372	}
373	G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart
374	G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon
375	if(W<0.)
376	{
377	//Do Nothing Action insead of the reaction
378	G4cout<<"G4QAtomElScat::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl;
379	aParticleChange.ProposeEnergy(kinEnergy);
380	aParticleChange.ProposeLocalEnergyDeposit(0.);
381	aParticleChange.ProposeMomentumDirection(dir) ;
382	return G4VDiscreteProcess::PostStepDoIt(track,step);
383	}
384	// Update G4VParticleChange for the scattered muon
385	G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer();
386	G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy, Z, N);// Integrated CS
387	G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N); // Real CrosSect
388	G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2);
389	if(sigNuG4UniformRand()>sigKrndFraction)
390	{
391	//Do NothingToDo Action insead of the reaction
392	G4cout<<"G4QAtomElectScat::PostStepDoIt: probability correction - DoNothing"<<G4endl;
393	aParticleChange.ProposeEnergy(kinEnergy);
394	aParticleChange.ProposeLocalEnergyDeposit(0.);
395	aParticleChange.ProposeMomentumDirection(dir) ;
396	return G4VDiscreteProcess::PostStepDoIt(track,step);
397	}
398	G4double iniE=kinEnergy+mu; // Initial total energy of the muon
399	G4double finE=iniE-photonEnergy; // Final total energy of the muon
400	if(finE>0) aParticleChange.ProposeEnergy(finE) ;
401	else
402	{
403	aParticleChange.ProposeEnergy(0.) ;
404	aParticleChange.ProposeTrackStatus(fStopAndKill);
405	}
406	// Scatter the muon
407	G4double EEm=iniEfinE-mu2; // Just an intermediate value to avoid "2"
408	G4double iniP=std::sqrt(iniE*iniE-mu2); // Initial momentum of the electron
409	G4double finP=std::sqrt(finE*finE-mu2); // Final momentum of the electron
410	G4double cost=(EEm+EEm-photonQ2)/iniP/finP; // cos(theta) for the electron scattering
411	if(cost>1.) cost=1.; // To avoid the accuracy of calculation problem
412	//else if(cost>1.001) // @@ error report can be done, but not necessary
413	if(cost<-1.) cost=-1.; // To avoid the accuracy of calculation problem
414	//else if(cost<-1.001) // @@ error report can be done, but not necessary
415	// --- Example from electromagnetic physics --
416	//G4ThreeVector newMuonDirection(dirx,diry,dirz);
417	//newMuonDirection.rotateUz(dir);
418	//aParticleChange.ProposeMomentumDirection(newMuonDirection1) ;
419	// The scattering in respect to the derection of the incident muon is made impicitly:
420	G4ThreeVector ort=dir.orthogonal(); // Not normed orthogonal vector (!) (to dir)
421	G4ThreeVector ortx = ort.unit(); // First unit vector orthogonal to the direction
422	G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction
423	G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
424	G4double phi=twopi*G4UniformRand(); // phi of scattered electron
425	G4double sinx=sint*std::sin(phi); // x-component
426	G4double siny=sint*std::cos(phi); // y-component
427	G4ThreeVector findir=costdir+sinxortx+siny*orty;
428	aParticleChange.ProposeMomentumDirection(findir); // new direction for the muon
429	const G4ThreeVector photon3M=iniPdir-finPfindir;
430	projPDG=22;
431	proj4M=G4LorentzVector(photon3M,photon3M.mag());
432	}
433	G4int targPDG=90000000+Z*1000+N; // PDG Code of the target nucleus
434	G4QPDGCode targQPDG(targPDG);
435	G4double tM=targQPDG.GetMass();
436	EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
437	G4QHadronVector* output=new G4QHadronVector; // Prototype of the output G4QHadronVector
438	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin
439	//G4bool elF=false; // Flag of the ellastic scattering is "false" by default
440	//G4double eWei=1.;
441	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- End
442	#ifdef debug
443	G4cout<<"G4QAtomElScat::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
444	#endif
445	G4QHadron* pH = new G4QHadron(projPDG,proj4M); // ---> DELETED -->-- -+
446	if(momentum<1000.) // Condition for using G4QEnvironment (not G4QuasmonString) \|
447	{ // \|
448	G4QHadronVector projHV; // \|
449	projHV.push_back(pH); // DESTROYED over 2 lines -+ \|
450	G4QEnvironment* pan= new G4QEnvironment(projHV,targPDG);// ---> DELETED --->-----+ \| \|
451	std::for_each(projHV.begin(), projHV.end(), DeleteQHadron()); // <---<------<----+-+-+
452	projHV.clear(); // <------------<---------------<-------------------<------------+-+ .
453	#ifdef debug
454	G4cout<<"G4QAtomElectScat::PostStepDoIt: pPDG="<<projPDG<<", mp="<<mp<<G4endl;// \| .
455	#endif
456	try // \| .
457	{ // \| .
458	delete output; // \| .
459	output = pan->Fragment();// DESTROYED in the end of the LOOP work space \| .
460	} // \| .
461	catch (G4QException& error)// \| .
462	{ // \| .
463	//#ifdef pdebug
464	G4cerr<<"**G4QAtomElectScat::PostStepDoIt:G4QE Exception is catched"<<G4endl;//\| .
465	//#endif
466	G4Exception("G4QAtomElScat::PostStepDoIt:","27",FatalException,"CHIPScrash");//\| .
467	} // \| .
468	delete pan; // Delete the Nuclear Environment <--<--+ .
469	} // .
470	else // Use G4QuasmonString .
471	{ // ^
472	G4QuasmonString* pan= new G4QuasmonString(pH,false,targPDG,false);//-> DELETED --+ \|
473	delete pH; // --------<-------+---+
474	#ifdef debug
475	G4double mp=G4QPDGCode(projPDG).GetMass(); // Mass of the projectile particle \|
476	G4cout<<"G4QAtomElectScat::PostStepDoIt: pPDG="<<projPDG<<", pM="<<mp<<G4endl; //\|
477	#endif
478	//G4int tNH=0; // Prototype of the number of secondaries inOut\|
479	try // \|
480	{ // \|
481	delete output; // \|
482	output = pan->Fragment();// DESTROYED in the end of the LOOP work space \|
483	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin \|
484	//tNH=pan->GetNOfHadrons(); // For the test purposes of the String \|
485	//if(tNH==2) // At least 2 hadrons are in the Constr.Output \|
486	//{// \|
487	// elF=true; // Just put a flag for the ellastic Scattering \|
488	// delete output; // Delete a prototype of dummy G4QHadronVector \|
489	// output = pan->GetHadrons(); // DESTROYED in the end of the LOOP work space \|
490	//}// \|
491	//eWei=pan->GetWeight(); // Just an example for the weight of the event \|
492	#ifdef debug
493	//G4cout<<"=====>>G4QAtomElScat::PostStepDoIt:elF="<<elF<<",n="<<tNH<<G4endl;//\|
494	#endif
495	// @@@@@@@@@@@@@@ Temporary for the testing purposes --- End \|
496	} // \|
497	catch (G4QException& error)// \|
498	{ // \|
499	//#ifdef pdebug
500	G4cerr<<"**G4QAtomElectScat::PostStepDoIt: GEN Exception is catched"<<G4endl;//\|
501	//#endif
502	G4Exception("G4QAtomElSct::AtRestDoIt:","27",FatalException,"QString Excep");//\|
503	} // \|
504	delete pan; // Delete the Nuclear Environment ---<--+
505	}
506	aParticleChange.Initialize(track);
507	G4double localtime = track.GetGlobalTime();
508	G4ThreeVector position = track.GetPosition();
509	G4TouchableHandle trTouchable = track.GetTouchableHandle();
510	// ------------- From here the secondaries are filled -------------------------
511	G4int tNH = output->size(); // A#of hadrons in the output
512	aParticleChange.SetNumberOfSecondaries(tNH);
513	// Now add nuclear fragments
514	#ifdef debug
515	G4cout<<"G4QAtomElectronScat::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
516	#endif
517	G4int nOut=output->size(); // Real length of the output @@ Temporary
518	if(tNH==1) tNH=0; // @@ Temporary
519	if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QAtomElScat::PostStepDoIt: 2 # "<<nOut<<G4endl;
520	// Deal with ParticleChange final state interface to GEANT4 output of the process
521	//if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
522	if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
523	{
524	// Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
525	// Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
526	// The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
527	G4QHadron* hadr=output->operator[](i); // Pointer to the output hadron
528	G4int PDGCode = hadr->GetPDGCode();
529	G4int nFrag = hadr->GetNFragments();
530	#ifdef pdebug
531	G4cout<<"G4QAtomElectScat::AtRestDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag<<G4endl;
532	#endif
533	if(nFrag) // Skip intermediate (decayed) hadrons
534	{
535	#ifdef debug
536	G4cout<<"G4QAtomElScat::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
537	#endif
538	delete hadr;
539	continue;
540	}
541	G4DynamicParticle* theSec = new G4DynamicParticle;
542	G4ParticleDefinition* theDefinition;
543	if (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
544	else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
545	else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
546	else if(PDGCode==311 \|\| PDGCode==-311)
547	{
548	if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L
549	else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
550	}
551	else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
552	else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
553	else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
554	else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
555	else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
556	else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
557	{
558	G4int aZ = hadr->GetCharge();
559	G4int aA = hadr->GetBaryonNumber();
560	#ifdef pdebug
561	G4cout<<"G4QAtomicElectronScattering::AtRestDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
562	#endif
563	theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
564	}
565	//else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
566	else
567	{
568	#ifdef pdebug
569	G4cout<<"G4QAtomElectScat::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
570	#endif
571	theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
572	#ifdef pdebug
573	G4cout<<"G4QAtomElScat::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
574	#endif
575	}
576	if(!theDefinition)
577	{
578	G4cout<<"---Warning---G4QAtomElScattering::PostStepDoIt: drop PDG="<<PDGCode<<G4endl;
579	delete hadr;
580	continue;
581	}
582	#ifdef pdebug
583	G4cout<<"G4QAtomElScat::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
584	#endif
585	theSec->SetDefinition(theDefinition);
586	G4LorentzVector h4M=hadr->Get4Momentum();
587	EnMomConservation-=h4M;
588	#ifdef tdebug
589	G4cout<<"G4QCollis::PSDI:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
590	#endif
591	#ifdef debug
592	G4cout<<"G4QAtomElectScat::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
593	#endif
594	theSec->Set4Momentum(h4M); // ^
595	delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
596	#ifdef debug
597	G4ThreeVector curD=theSec->GetMomentumDirection(); // ^
598	G4double curM=theSec->GetMass(); // \|
599	G4double curE=theSec->GetKineticEnergy()+curM; // ^
600	G4cout<<"G4QCollis::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;// \|
601	#endif
602	G4Track* aNewTrack = new G4Track(theSec, localtime, position ); // ^
603	aNewTrack->SetTouchableHandle(trTouchable); // \|
604	aParticleChange.AddSecondary( aNewTrack ); // \|
605	#ifdef debug
606	G4cout<<"G4QAtomicElectronScattering::PostStepDoIt:#"<<i<<" is done."<<G4endl; // \|
607	#endif
608	} // \|
609	delete output; // instances of the G4QHadrons from the output are already deleted above +
610	aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the absorbed particle
611	//return &aParticleChange; // This is not enough (ClearILL)
612	#ifdef debug
613	G4cout<<"G4QAtomicElectronScattering::PostStepDoIt:**PostStepDoIt done**"<<G4endl;
614	#endif
615	return G4VDiscreteProcess::PostStepDoIt(track, step);
616	}

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

Download in other formats: