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

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

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