Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4NeutrinoElectronProcess.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//
27// Geant4 Hadron Elastic Scattering Process
28//
29// Created from G4HadronElasticProcess
30//
31// Modified:
32//
33// 2.2.18 V.Grichine - PostStepDoIt implementation
34// 03.10.18 V. Grichine - G4Region name and optionally total cross section biased in the region only.
35#include <iostream>
36#include <typeinfo>
37
39#include "G4SystemOfUnits.hh"
40#include "G4Nucleus.hh"
41#include "G4ProcessManager.hh"
47#include "G4VDiscreteProcess.hh"
48
50
51#include "G4RotationMatrix.hh"
52#include "G4ThreeVector.hh"
53#include "G4AffineTransform.hh"
54#include "G4DynamicParticle.hh"
55#include "G4StepPoint.hh"
56#include "G4VSolid.hh"
57#include "G4LogicalVolume.hh"
58#include "G4SafetyHelper.hh"
60
61///////////////////////////////////////////////////////////////////////////////
62
63
66{
67 lowestEnergy = 1.*keV;
68 fEnvelopeName = anEnvelopeName;
69 fTotXsc = new G4NeutrinoElectronTotXsc();
71 safetyHelper->InitialiseHelper();
72}
73
74///////////////////////////////////////////////////////
75
77{
78 fNuEleTotXscBias = bf;
79}
80
81///////////////////////////////////////////////////////
82
84{
85 fNuEleCcBias = bfCc;
86 fNuEleNcBias = bfNc;
87 fNuEleTotXscBias = std::max( fNuEleCcBias, fNuEleNcBias );
88}
89
90///////////////////////////////////////////////
91
99
100//////////////////////////////////////////////////
101
104{
105 //G4cout << "GetMeanFreePath " << aTrack.GetDefinition()->GetParticleName()
106 // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
108 G4double totxsc =
110 aTrack.GetMaterial());
111
112 if ( rName == fEnvelopeName )
113 {
114 totxsc *= fNuEleTotXscBias;
115 }
116 G4double res = (totxsc>0.0) ? 1.0/totxsc : DBL_MAX;
117 //G4cout << " xsection= " << totxsc << G4endl;
118 return res;
119}
120
121///////////////////////////////////////////////////
122
123void G4NeutrinoElectronProcess::ProcessDescription(std::ostream& outFile) const
124{
125 outFile << "G4NeutrinoElectronProcess handles the scattering of \n"
126 << "neutrino on electrons by invoking the following model(s) and \n"
127 << "cross section(s).\n";
128}
129
130///////////////////////////////////////////////////////////////////////
131
134{
136
137 if( rName != fEnvelopeName )
138 {
139 if( verboseLevel > 0 )
140 {
141 G4cout<<"Go out from G4NeutrinoElectronProcess::PostStepDoIt: wrong volume "<<G4endl;
142 }
143 return G4VDiscreteProcess::PostStepDoIt( track, step );
144 }
147 G4double weight = track.GetWeight();
149
150 if( track.GetTrackStatus() != fAlive )
151 {
152 return theTotalResult;
153 }
154 // Next check for illegal track status
155 //
156 if (track.GetTrackStatus() != fAlive &&
157 track.GetTrackStatus() != fSuspend)
158 {
159 if (track.GetTrackStatus() == fStopAndKill ||
162 {
164 ed << "G4HadronicProcess: track in unusable state - "
165 << track.GetTrackStatus() << G4endl;
166 ed << "G4HadronicProcess: returning unchanged track " << G4endl;
167 DumpState(track,"PostStepDoIt",ed);
168 G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
169 }
170 // No warning for fStopButAlive which is a legal status here
171 return theTotalResult;
172 }
173
174 // For elastic scattering, _any_ result is considered an interaction
176
177 G4double kineticEnergy = track.GetKineticEnergy();
178 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
179 const G4ParticleDefinition* part = dynParticle->GetDefinition();
180
181 // NOTE: Very low energy scatters were causing numerical (FPE) errors
182 // in earlier releases; these limits have not been changed since.
183
184 if ( kineticEnergy <= lowestEnergy ) return theTotalResult;
185
186 const G4Material* material = track.GetMaterial();
187 G4Nucleus* targNucleus = GetTargetNucleusPointer();
188
189 //////////////// uniform random spread of the neutrino interaction point ////////////
190
191 const G4StepPoint* pPostStepPoint = step.GetPostStepPoint();
192 const G4DynamicParticle* aParticle = track.GetDynamicParticle();
193 G4ThreeVector position = pPostStepPoint->GetPosition(), newPosition=position;
194 G4ParticleMomentum direction = aParticle->GetMomentumDirection();
195
196 if( fNuEleCcBias > 1.0 || fNuEleNcBias > 1.0) // = true, if fBiasingfactor != 1., i.e. xsc is biased
197 {
198 const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
199 G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
200 G4AffineTransform transform = G4AffineTransform(rotM,transl);
201 transform.Invert();
202
203 G4ThreeVector localP = transform.TransformPoint(position);
204 G4ThreeVector localV = transform.TransformAxis(direction);
205
206 G4double forward = track.GetVolume()->GetLogicalVolume()->GetSolid()->DistanceToOut(localP, localV);
207 G4double backward = track.GetVolume()->GetLogicalVolume()->GetSolid()->DistanceToOut(localP, -localV);
208
209 G4double distance = forward+backward;
210
211 // G4cout<<distance/cm<<", ";
212
213 // uniform sampling of nu-e interaction point
214 // along neutrino direction in current volume
215
216 G4double range = -backward+G4UniformRand()*distance;
217
218 newPosition = position + range*direction;
219
220 safetyHelper->ReLocateWithinVolume(newPosition);
221
222 theTotalResult->ProposePosition(newPosition); // G4Exception : GeomNav1002
223 }
224 G4HadProjectile theProj( track );
225 G4HadronicInteraction* hadi = nullptr;
226 G4HadFinalState* result = nullptr;
227
228 // Select element
229 const G4Element* elm = nullptr;
230
231 try
232 {
233 elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
234 *targNucleus);
235 }
236 catch( G4HadronicException & aR )
237 {
239 aR.Report(ed);
240 DumpState(track,"SampleZandA",ed);
241 ed << " PostStepDoIt failed on element selection" << G4endl;
242 G4Exception("G4NeutrinoElectronProcess::PostStepDoIt", "had003",
243 FatalException, ed);
244 }
245
246 G4double ccTotRatio = fTotXsc->GetCcRatio();
247
248 if( G4UniformRand() < ccTotRatio ) // Cc-model
249 {
250 // Initialize the hadronic projectile from the track
251 thePro.Initialise(track);
252
253 hadi = (GetHadronicInteractionList())[0];
254
255 result = hadi->ApplyYourself( thePro, *targNucleus);
256
258
260
261 FillResult(result, track);
262 }
263 else // Nc-model, like 'elastic', 2->2 scattering
264 {
265
266 hadi = (GetHadronicInteractionList())[1];
267
268 size_t idx = track.GetMaterialCutsCouple()->GetIndex();
269
271
272 hadi->SetRecoilEnergyThreshold(tcut);
273
274 if( verboseLevel > 1 )
275 {
276 G4cout << "G4NeutrinoElectronProcess::PostStepDoIt for "
277 << part->GetParticleName()
278 << " in " << material->GetName()
279 << " Target Z= " << targNucleus->GetZ_asInt()
280 << " A= " << targNucleus->GetA_asInt() << G4endl;
281 }
282 try
283 {
284 result = hadi->ApplyYourself( theProj, *targNucleus);
285 }
286 catch(G4HadronicException & aR)
287 {
289 aR.Report(ed);
290 ed << "Call for " << hadi->GetModelName() << G4endl;
291 ed << "Target element "<< elm->GetName()<<" Z= "
292 << targNucleus->GetZ_asInt()
293 << " A= " << targNucleus->GetA_asInt() << G4endl;
294 DumpState(track,"ApplyYourself",ed);
295 ed << " ApplyYourself failed" << G4endl;
296 G4Exception("G4NeutrinoElectronProcess::PostStepDoIt", "had006",
297 FatalException, ed);
298 }
299 // directions
300
301 G4ThreeVector indir = track.GetMomentumDirection();
302 G4double phi = CLHEP::twopi*G4UniformRand();
303 G4ThreeVector it(0., 0., 1.);
304 G4ThreeVector outdir = result->GetMomentumChange();
305
306 if(verboseLevel>1)
307 {
308 G4cout << "Efin= " << result->GetEnergyChange()
309 << " de= " << result->GetLocalEnergyDeposit()
310 << " nsec= " << result->GetNumberOfSecondaries()
311 << " dir= " << outdir
312 << G4endl;
313 }
314 // energies
315
316 G4double edep = result->GetLocalEnergyDeposit();
317 G4double efinal = result->GetEnergyChange();
318
319 if(efinal < 0.0) { efinal = 0.0; }
320 if(edep < 0.0) { edep = 0.0; }
321
322 // NOTE: Very low energy scatters were causing numerical (FPE) errors
323 // in earlier releases; these limits have not been changed since.
324
325 if(efinal <= lowestEnergy)
326 {
327 edep += efinal;
328 efinal = 0.0;
329 }
330 // primary change
331
333
334 G4TrackStatus status = track.GetTrackStatus();
335
336 if(efinal > 0.0)
337 {
338 outdir.rotate(phi, it);
339 outdir.rotateUz(indir);
341 }
342 else
343 {
344 if( part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
345 {
346 status = fStopButAlive;
347 }
348 else
349 {
350 status = fStopAndKill;
351 }
353 }
354 //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
355
357
358 // recoil
359
360 if( result->GetNumberOfSecondaries() > 0 )
361 {
362 G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
363
364 if(p->GetKineticEnergy() > tcut)
365 {
368
369 // G4cout << "recoil " << pdir << G4endl;
370 //!! is not needed for models inheriting G4NeutrinoElectron
371
372 pdir.rotate(phi, it);
373 pdir.rotateUz(indir);
374
375 // G4cout << "recoil rotated " << pdir << G4endl;
376
377 p->SetMomentumDirection(pdir);
378
379 // in elastic scattering time and weight are not changed
380
381 G4Track* t = new G4Track(p, track.GetGlobalTime(),
382 track.GetPosition());
383 t->SetWeight(weight);
386 }
387 else
388 {
389 edep += p->GetKineticEnergy();
390 delete p;
391 }
392 }
395 result->Clear();
396 }
397 return theTotalResult;
398}
399
400void
402{
403 lowestEnergy = val;
404}
405
G4double condition(const G4ErrorSymMatrix &m)
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4ForceCondition
G4TrackStatus
@ fKillTrackAndSecondaries
@ fSuspend
@ fAlive
@ fStopAndKill
@ fStopButAlive
@ fPostponeToNextEvent
double G4double
Definition G4Types.hh:83
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition Randomize.hh:52
Hep3Vector & rotateUz(const Hep3Vector &)
Hep3Vector & rotate(double, const Hep3Vector &)
G4AffineTransform & Invert()
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
G4double ComputeCrossSection(const G4DynamicParticle *, const G4Material *)
const G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
void SetMomentumDirection(const G4ThreeVector &aDirection)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4String & GetName() const
Definition G4Element.hh:115
G4double GetEnergyChange() const
void SetTrafoToLab(const G4LorentzRotation &aT)
G4double GetLocalEnergyDeposit() const
const G4ThreeVector & GetMomentumChange() const
std::size_t GetNumberOfSecondaries() const
G4HadSecondary * GetSecondary(size_t i)
void Initialise(const G4Track &aT)
G4LorentzRotation & GetTrafoToLab()
G4DynamicParticle * GetParticle()
void Report(std::ostream &aS) const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
const G4String & GetModelName() const
void SetRecoilEnergyThreshold(G4double val)
void FillResult(G4HadFinalState *aR, const G4Track &aT)
G4HadProjectile thePro
G4Nucleus * GetTargetNucleusPointer()
G4ParticleChange * theTotalResult
std::vector< G4HadronicInteraction * > & GetHadronicInteractionList()
G4CrossSectionDataStore * GetCrossSectionDataStore()
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
G4VSolid * GetSolid() const
G4Region * GetRegion() const
const G4String & GetName() const
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
void SetBiasingFactors(G4double bfCc, G4double bfNc)
G4NeutrinoElectronProcess(const G4String &anEnvelopeName, const G4String &procName="nuElectron")
void ProcessDescription(std::ostream &outFile) const override
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *) override
G4int GetA_asInt() const
Definition G4Nucleus.hh:99
G4int GetZ_asInt() const
Definition G4Nucleus.hh:105
void AddSecondary(G4Track *aSecondary)
void ProposePosition(G4double x, G4double y, G4double z)
void Initialize(const G4Track &) override
void ProposeEnergy(G4double finalEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4ProcessManager * GetProcessManager() const
const G4String & GetParticleName() const
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
std::size_t size() const
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
static G4ProductionCutsTable * GetProductionCutsTable()
const G4String & GetName() const
void ReLocateWithinVolume(const G4ThreeVector &pGlobalPoint)
const G4VTouchable * GetTouchable() const
const G4ThreeVector & GetPosition() const
G4VPhysicalVolume * GetPhysicalVolume() const
G4StepPoint * GetPreStepPoint() const
G4StepPoint * GetPostStepPoint() const
virtual const G4RotationMatrix * GetRotation(G4int depth=0) const
virtual const G4ThreeVector & GetTranslation(G4int depth=0) const
G4TrackStatus GetTrackStatus() const
G4VPhysicalVolume * GetVolume() const
G4double GetWeight() const
void SetWeight(G4double aValue)
const G4ThreeVector & GetPosition() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4double GetGlobalTime() const
G4Material * GetMaterial() const
const G4DynamicParticle * GetDynamicParticle() const
const G4TouchableHandle & GetTouchableHandle() const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const G4Step * GetStep() const
static G4TransportationManager * GetTransportationManager()
G4SafetyHelper * GetSafetyHelper() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
void ProposeWeight(G4double finalWeight)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetNumberOfSecondaries(G4int totSecondaries)
G4LogicalVolume * GetLogicalVolume() const
void ClearNumberOfInteractionLengthLeft()
G4int verboseLevel
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const =0
#define DBL_MAX
Definition templates.hh:62