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