Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4NeutrinoElectronProcess.cc
이 파일의 문서화 페이지로 가기
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: G4NeutrinoElectronProcess.cc 92396 2015-08-31 14:12:40Z gcosmo $
27 //
28 // Geant4 Hadron Elastic Scattering Process
29 //
30 // Created from G4HadronElasticProcess
31 //
32 // Modified:
33 //
34 // 2.2.18 V.Grichine - PostStepDoIt implementation
35 
36 #include <iostream>
37 #include <typeinfo>
38 
40 #include "G4SystemOfUnits.hh"
41 #include "G4Nucleus.hh"
42 #include "G4ProcessManager.hh"
44 #include "G4HadronElasticDataSet.hh" //???
45 #include "G4ProductionCutsTable.hh"
46 #include "G4HadronicException.hh"
47 #include "G4HadronicDeprecate.hh"
48 #include "G4HadronicInteraction.hh"
49 #include "G4VCrossSectionRatio.hh"
50 #include "G4VDiscreteProcess.hh"
51 
53 //#include "G4NeutrinoElectronCcModel.hh"
54 //#include "G4NeutrinoElectronNcModel.hh"
55 
56 #include "G4RotationMatrix.hh"
57 #include "G4ThreeVector.hh"
58 #include "G4AffineTransform.hh"
59 #include "G4DynamicParticle.hh"
60 #include "G4StepPoint.hh"
61 #include "G4VSolid.hh"
62 #include "G4LogicalVolume.hh"
63 #include "G4SafetyHelper.hh"
65 
67 
68 
70  : G4HadronicProcess( pName, fHadronElastic ), isInitialised(false), fBiased(true) // fHadronElastic???
71 {
72  // AddDataSet(new G4HadronElasticDataSet); //???
73  lowestEnergy = 1.*keV;
74  fEnvelope = nullptr;
75  fEnvelopeName = anEnvelopeName;
76  fTotXsc = nullptr; // new G4NeutrinoElectronTotXsc();
77  fNuEleCcBias=1.;
78  fNuEleNcBias=1.;
81 }
82 
84 {
85  // if( fTotXsc ) delete fTotXsc;
86 }
88 {
89  fNuEleCcBias=bfCc;
90  fNuEleNcBias=bfNc;
91 
93  fTotXsc->SetBiasingFactors(bfCc, bfNc);
94 }
95 
97 
98 void G4NeutrinoElectronProcess::ProcessDescription(std::ostream& outFile) const
99 {
100 
101  outFile << "G4NeutrinoElectronProcess handles the scattering of \n"
102  << "neutrino on electrons by invoking the following model(s) and \n"
103  << "cross section(s).\n";
104 
105 }
106 
108 
111 {
112  // if( track.GetVolume()->GetLogicalVolume() != fEnvelope )
113  if( track.GetVolume()->GetLogicalVolume()->GetName() != fEnvelopeName )
114  {
115  if( verboseLevel > 0 )
116  {
117  G4cout<<"Go out from G4NeutrinoElectronProcess::PostStepDoIt: wrong volume "<<G4endl;
118  }
119  return G4VDiscreteProcess::PostStepDoIt( track, step );
120  }
122  theTotalResult->Initialize(track);
123  G4double weight = track.GetWeight();
124  theTotalResult->ProposeWeight(weight);
125 
126  if( track.GetTrackStatus() != fAlive )
127  {
128  return theTotalResult;
129  }
130  // Next check for illegal track status
131  //
132  if (track.GetTrackStatus() != fAlive &&
133  track.GetTrackStatus() != fSuspend)
134  {
135  if (track.GetTrackStatus() == fStopAndKill ||
138  {
140  ed << "G4HadronicProcess: track in unusable state - "
141  << track.GetTrackStatus() << G4endl;
142  ed << "G4HadronicProcess: returning unchanged track " << G4endl;
143  DumpState(track,"PostStepDoIt",ed);
144  G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
145  }
146  // No warning for fStopButAlive which is a legal status here
147  return theTotalResult;
148  }
149 
150  // For elastic scattering, _any_ result is considered an interaction
152 
154  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
155  const G4ParticleDefinition* part = dynParticle->GetDefinition();
156 
157  // NOTE: Very low energy scatters were causing numerical (FPE) errors
158  // in earlier releases; these limits have not been changed since.
159 
160  if ( kineticEnergy <= lowestEnergy ) return theTotalResult;
161 
162  const G4Material* material = track.GetMaterial();
163  G4Nucleus* targNucleus = GetTargetNucleusPointer();
164 
166 
167  const G4StepPoint* pPostStepPoint = step.GetPostStepPoint();
168  const G4DynamicParticle* aParticle = track.GetDynamicParticle();
169  G4ThreeVector position = pPostStepPoint->GetPosition(), newPosition=position;
170  G4ParticleMomentum direction = aParticle->GetMomentumDirection();
171  G4double startTime = pPostStepPoint->GetGlobalTime();
172 
173 
174  if( fNuEleCcBias > 1.0 || fNuEleNcBias > 1.0) // = true, if fBiasingfactor != 1., i.e. xsc is biased
175  {
176  const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
177  G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
178  G4AffineTransform transform = G4AffineTransform(rotM,transl);
179  transform.Invert();
180 
181  G4ThreeVector localP = transform.TransformPoint(position);
182  G4ThreeVector localV = transform.TransformAxis(direction);
183 
184  G4double forward = track.GetVolume()->GetLogicalVolume()->GetSolid()->DistanceToOut(localP, localV);
185  G4double backward = track.GetVolume()->GetLogicalVolume()->GetSolid()->DistanceToOut(localP, -localV);
186 
187  G4double distance = forward+backward;
188 
189  // G4cout<<distance/cm<<", ";
190 
191  // uniform sampling of nu-e interaction point
192  // along neutrino direction in current volume
193 
194  G4double range = -backward+G4UniformRand()*distance;
195 
196  G4double delta = range - backward;
197 
198  startTime += delta/track.GetVelocity();
199 
200  newPosition = position + range*direction;
201 
202  safetyHelper->ReLocateWithinVolume(newPosition);
203 
204  theTotalResult->ProposePosition(newPosition); // G4Exception : GeomNav1002
205  // theTotalResult->ProposeGlobalTime(startTime); // time is updated for 'elastic' only
206  }
207  G4HadProjectile theProj( track );
208  G4HadronicInteraction* hadi = nullptr;
209  G4HadFinalState* result = nullptr;
210 
211  // Select element
212  const G4Element* elm = nullptr;
213  G4int ZZ=1;
214 
215  try
216  {
217  elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
218  *targNucleus);
219  }
220  catch( G4HadronicException & aR )
221  {
223  aR.Report(ed);
224  DumpState(track,"SampleZandA",ed);
225  ed << " PostStepDoIt failed on element selection" << G4endl;
226  G4Exception("G4NeutrinoElectronProcess::PostStepDoIt", "had003",
227  FatalException, ed);
228  }
229  if( elm ) ZZ = elm->GetZ();
230 
231  G4double xsc = fTotXsc->GetElementCrossSection(dynParticle, ZZ, material);
232  xsc *= 1.;
233  G4double ccTotRatio = fTotXsc->GetCcRatio();
234 
235  if( G4UniformRand() < ccTotRatio ) // Cc-model
236  {
237  // Initialize the hadronic projectile from the track
238  thePro.Initialise(track);
239 
240  hadi = (GetHadronicInteractionList())[0];
241 
242  result = hadi->ApplyYourself( thePro, *targNucleus);
243 
245 
247 
248  FillResult(result, track);
249  }
250  else // Nc-model, like 'elastic', 2->2 scattering
251  {
252 
253  hadi = (GetHadronicInteractionList())[1];
254 
255  size_t idx = track.GetMaterialCutsCouple()->GetIndex();
256 
258 
259  hadi->SetRecoilEnergyThreshold(tcut);
260 
261  if( verboseLevel > 1 )
262  {
263  G4cout << "G4NeutrinoElectronProcess::PostStepDoIt for "
264  << part->GetParticleName()
265  << " in " << material->GetName()
266  << " Target Z= " << targNucleus->GetZ_asInt()
267  << " A= " << targNucleus->GetA_asInt() << G4endl;
268  }
269  try
270  {
271  result = hadi->ApplyYourself( theProj, *targNucleus);
272  }
273  catch(G4HadronicException & aR)
274  {
276  aR.Report(ed);
277  ed << "Call for " << hadi->GetModelName() << G4endl;
278  ed << "Target element "<< elm->GetName()<<" Z= "
279  << targNucleus->GetZ_asInt()
280  << " A= " << targNucleus->GetA_asInt() << G4endl;
281  DumpState(track,"ApplyYourself",ed);
282  ed << " ApplyYourself failed" << G4endl;
283  G4Exception("G4NeutrinoElectronProcess::PostStepDoIt", "had006",
284  FatalException, ed);
285  }
286  // directions
287 
288  G4ThreeVector indir = track.GetMomentumDirection();
290  G4ThreeVector it(0., 0., 1.);
291  G4ThreeVector outdir = result->GetMomentumChange();
292 
293  if(verboseLevel>1)
294  {
295  G4cout << "Efin= " << result->GetEnergyChange()
296  << " de= " << result->GetLocalEnergyDeposit()
297  << " nsec= " << result->GetNumberOfSecondaries()
298  << " dir= " << outdir
299  << G4endl;
300  }
301  // energies
302 
303  G4double edep = result->GetLocalEnergyDeposit();
304  G4double efinal = result->GetEnergyChange();
305 
306  if(efinal < 0.0) { efinal = 0.0; }
307  if(edep < 0.0) { edep = 0.0; }
308 
309  // NOTE: Very low energy scatters were causing numerical (FPE) errors
310  // in earlier releases; these limits have not been changed since.
311 
312  if(efinal <= lowestEnergy)
313  {
314  edep += efinal;
315  efinal = 0.0;
316  }
317  // primary change
318 
319  theTotalResult->ProposeEnergy(efinal);
320 
321  G4TrackStatus status = track.GetTrackStatus();
322 
323  if(efinal > 0.0)
324  {
325  outdir.rotate(phi, it);
326  outdir.rotateUz(indir);
328  }
329  else
330  {
331  if( part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
332  {
333  status = fStopButAlive;
334  }
335  else
336  {
337  status = fStopAndKill;
338  }
340  }
341  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
342 
344 
345  // recoil
346 
347  if( result->GetNumberOfSecondaries() > 0 )
348  {
349  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
350 
351  if(p->GetKineticEnergy() > tcut)
352  {
355 
356  // G4cout << "recoil " << pdir << G4endl;
358 
359  pdir.rotate(phi, it);
360  pdir.rotateUz(indir);
361 
362  // G4cout << "recoil rotated " << pdir << G4endl;
363 
364  p->SetMomentumDirection(pdir);
365 
366  // in elastic scattering time and weight are not changed
367 
368  G4Track* t = new G4Track(p, track.GetGlobalTime(),
369  track.GetPosition());
370  t->SetWeight(weight);
373  }
374  else
375  {
376  edep += p->GetKineticEnergy();
377  delete p;
378  }
379  }
382  result->Clear();
383  }
384  return theTotalResult;
385 }
386 
387 void
389 {
390  if(!isInitialised) {
391  isInitialised = true;
392  if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
393  }
395 }
396 
397 void
399 {
400  lowestEnergy = val;
401 }
402 
G4double GetKineticEnergy() const
G4ParticleChange * theTotalResult
void SetTrafoToLab(const G4LorentzRotation &aT)
G4Nucleus * GetTargetNucleusPointer()
void ReLocateWithinVolume(const G4ThreeVector &pGlobalPoint)
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
const G4VTouchable * GetTouchable() const
G4SafetyHelper * GetSafetyHelper() const
G4LorentzRotation & GetTrafoToLab()
G4HadProjectile thePro
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
const G4ThreeVector & GetMomentumChange() const
G4LogicalVolume * GetLogicalVolume() const
void SetWeight(G4double aValue)
const G4String & GetName() const
Definition: G4Element.hh:127
static constexpr double keV
Definition: G4SIunits.hh:216
G4HadSecondary * GetSecondary(size_t i)
void ProposeWeight(G4double finalWeight)
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const char * p
Definition: xmltok.h:285
const G4String & GetParticleName() const
const G4TouchableHandle & GetTouchableHandle() const
virtual const G4ThreeVector & GetTranslation(G4int depth=0) const =0
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
virtual const G4RotationMatrix * GetRotation(G4int depth=0) const =0
void PreparePhysicsTable(const G4ParticleDefinition &) override
virtual void SetLowestEnergy(G4double)
void SetBiasingFactors(G4double bfCc, G4double bfNc)
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
void AddSecondary(G4Track *aSecondary)
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4double GetWeight() const
const G4String & GetName() const
Definition: G4Material.hh:179
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
void Initialise(const G4Track &aT)
double G4double
Definition: G4Types.hh:76
TString part[npart]
Definition: Style.C:32
G4ParticleDefinition * GetDefinition() const
#define position
Definition: xmlparse.cc:622
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
virtual void Initialize(const G4Track &)
const G4ParticleDefinition const G4Material *G4double range
const G4ThreeVector & GetPosition() const
void ProposePosition(G4double x, G4double y, G4double z)
G4NeutrinoElectronProcess(G4String anEnvelopeName, const G4String &procName="neutrino-electron")
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *)
#define G4UniformRand()
Definition: Randomize.hh:53
G4double GetGlobalTime() const
G4AffineTransform & Invert()
void SetTouchableHandle(const G4TouchableHandle &apValue)
const G4ThreeVector & GetPosition() const
static constexpr double eV
Definition: G4SIunits.hh:215
G4double GetGlobalTime() const
Definition: G4Step.hh:76
G4StepPoint * GetPostStepPoint() const
void ProposeEnergy(G4double finalEnergy)
G4Material * GetMaterial() const
void InitialiseHelper()
G4double GetVelocity() const
void FillResult(G4HadFinalState *aR, const G4Track &aT)
void SetRecoilEnergyThreshold(G4double val)
Double_t edep
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4int GetNumberOfSecondaries() const
G4double G4ParticleHPJENDLHEData::G4double result
static G4TransportationManager * GetTransportationManager()
virtual void ProcessDescription(std::ostream &outFile) const
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4double GetLocalEnergyDeposit() const
static G4ProductionCutsTable * GetProductionCutsTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
G4CrossSectionDataStore * GetCrossSectionDataStore()
double weight
Definition: plottest35.C:25
int G4int
Definition: G4Types.hh:78
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4ProcessManager * GetProcessManager() const
G4int verboseLevel
Definition: G4VProcess.hh:371
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4DynamicParticle * GetParticle()
G4double GetKineticEnergy() const
G4int size() const
G4GLOB_DLL std::ostream G4cout
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:450
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4VSolid * GetSolid() const
G4double GetZ() const
Definition: G4Element.hh:131
void Report(std::ostream &aS)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double GetEnergyChange() const
G4TrackStatus GetTrackStatus() const
const G4ThreeVector & GetMomentumDirection() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4NeutrinoElectronTotXsc * fTotXsc
void ProposeTrackStatus(G4TrackStatus status)
const G4DynamicParticle * GetDynamicParticle() const
void SetBiasingFactors(G4double bfCc, G4double bfNc)
std::vector< G4HadronicInteraction * > & GetHadronicInteractionList()
const G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4TrackStatus
const G4String & GetName() const
G4VPhysicalVolume * GetVolume() const
const G4String & GetModelName() const
void SetNumberOfSecondaries(G4int totSecondaries)