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

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago

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