ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutrinoElectronProcess.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
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"
43 #include "G4HadronElasticDataSet.hh" //???
44 #include "G4ProductionCutsTable.hh"
45 #include "G4HadronicException.hh"
46 #include "G4HadronicDeprecate.hh"
47 #include "G4HadronicInteraction.hh"
48 #include "G4VCrossSectionRatio.hh"
49 #include "G4VDiscreteProcess.hh"
50 
52 //#include "G4NeutrinoElectronCcModel.hh"
53 //#include "G4NeutrinoElectronNcModel.hh"
54 
55 #include "G4RotationMatrix.hh"
56 #include "G4ThreeVector.hh"
57 #include "G4AffineTransform.hh"
58 #include "G4DynamicParticle.hh"
59 #include "G4StepPoint.hh"
60 #include "G4VSolid.hh"
61 #include "G4LogicalVolume.hh"
62 #include "G4SafetyHelper.hh"
64 
66 
67 
69  : G4HadronicProcess( pName, fHadronElastic ), isInitialised(false), fBiased(true) // fHadronElastic???
70 {
71  // AddDataSet(new G4HadronElasticDataSet); //???
72  lowestEnergy = 1.*keV;
73  fEnvelope = nullptr;
74  fEnvelopeName = anEnvelopeName;
75  fTotXsc = nullptr; // new G4NeutrinoElectronTotXsc();
76  fNuEleCcBias=1.;
77  fNuEleNcBias=1.;
81 }
82 
84 {
85  // if( fTotXsc ) delete fTotXsc;
86 }
87 
89 
91 {
92  fNuEleTotXscBias = bf;
93 
95  // fTotXsc->SetBiasingFactor(bf);
96 }
97 
99 
101 {
102  fNuEleCcBias=bfCc;
103  fNuEleNcBias=bfNc;
104 
106  fTotXsc->SetBiasingFactors(bfCc, bfNc);
107 }
108 
110 
113 {
114  //G4cout << "GetMeanFreePath " << aTrack.GetDefinition()->GetParticleName()
115  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
117  G4double totxsc(0.);
118  try
119  {
120  if( rName == fEnvelopeName && fNuEleTotXscBias > 1.)
121  {
122  totxsc = fNuEleTotXscBias*
124  aTrack.GetMaterial());
125  }
126  else
127  {
129  aTrack.GetMaterial());
130  }
131  }
132  catch(G4HadronicException & aR)
133  {
135  aR.Report(ed);
136  DumpState(aTrack,"GetMeanFreePath",ed);
137  ed << " Cross section is not available" << G4endl;
138  G4Exception("G4NeutrinoElectronProcess::GetMeanFreePath", "had002", FatalException,
139  ed);
140  }
141  G4double res = (totxsc>0.0) ? 1.0/totxsc : DBL_MAX;
142  //G4cout << " xsection= " << totxsc << G4endl;
143  return res;
144 }
145 
147 
148 void G4NeutrinoElectronProcess::ProcessDescription(std::ostream& outFile) const
149 {
150 
151  outFile << "G4NeutrinoElectronProcess handles the scattering of \n"
152  << "neutrino on electrons by invoking the following model(s) and \n"
153  << "cross section(s).\n";
154 
155 }
156 
158 
161 {
162  // track.GetVolume()->GetLogicalVolume()->GetName()
163  // if( track.GetVolume()->GetLogicalVolume() != fEnvelope )
164 
166 
167  if( rName != fEnvelopeName )
168  {
169  if( verboseLevel > 0 )
170  {
171  G4cout<<"Go out from G4NeutrinoElectronProcess::PostStepDoIt: wrong volume "<<G4endl;
172  }
173  return G4VDiscreteProcess::PostStepDoIt( track, step );
174  }
176  theTotalResult->Initialize(track);
177  G4double weight = track.GetWeight();
178  theTotalResult->ProposeWeight(weight);
179 
180  if( track.GetTrackStatus() != fAlive )
181  {
182  return theTotalResult;
183  }
184  // Next check for illegal track status
185  //
186  if (track.GetTrackStatus() != fAlive &&
187  track.GetTrackStatus() != fSuspend)
188  {
189  if (track.GetTrackStatus() == fStopAndKill ||
192  {
194  ed << "G4HadronicProcess: track in unusable state - "
195  << track.GetTrackStatus() << G4endl;
196  ed << "G4HadronicProcess: returning unchanged track " << G4endl;
197  DumpState(track,"PostStepDoIt",ed);
198  G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
199  }
200  // No warning for fStopButAlive which is a legal status here
201  return theTotalResult;
202  }
203 
204  // For elastic scattering, _any_ result is considered an interaction
206 
207  G4double kineticEnergy = track.GetKineticEnergy();
208  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
209  const G4ParticleDefinition* part = dynParticle->GetDefinition();
210 
211  // NOTE: Very low energy scatters were causing numerical (FPE) errors
212  // in earlier releases; these limits have not been changed since.
213 
214  if ( kineticEnergy <= lowestEnergy ) return theTotalResult;
215 
216  const G4Material* material = track.GetMaterial();
217  G4Nucleus* targNucleus = GetTargetNucleusPointer();
218 
220 
221  const G4StepPoint* pPostStepPoint = step.GetPostStepPoint();
222  const G4DynamicParticle* aParticle = track.GetDynamicParticle();
223  G4ThreeVector position = pPostStepPoint->GetPosition(), newPosition=position;
224  G4ParticleMomentum direction = aParticle->GetMomentumDirection();
225  G4double startTime = pPostStepPoint->GetGlobalTime();
226 
227 
228  if( fNuEleCcBias > 1.0 || fNuEleNcBias > 1.0) // = true, if fBiasingfactor != 1., i.e. xsc is biased
229  {
230  const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
231  G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
233  transform.Invert();
234 
235  G4ThreeVector localP = transform.TransformPoint(position);
236  G4ThreeVector localV = transform.TransformAxis(direction);
237 
238  G4double forward = track.GetVolume()->GetLogicalVolume()->GetSolid()->DistanceToOut(localP, localV);
239  G4double backward = track.GetVolume()->GetLogicalVolume()->GetSolid()->DistanceToOut(localP, -localV);
240 
241  G4double distance = forward+backward;
242 
243  // G4cout<<distance/cm<<", ";
244 
245  // uniform sampling of nu-e interaction point
246  // along neutrino direction in current volume
247 
248  G4double range = -backward+G4UniformRand()*distance;
249 
250  G4double delta = range - backward;
251 
252  startTime += delta/track.GetVelocity();
253 
254  newPosition = position + range*direction;
255 
256  safetyHelper->ReLocateWithinVolume(newPosition);
257 
258  theTotalResult->ProposePosition(newPosition); // G4Exception : GeomNav1002
259  // theTotalResult->ProposeGlobalTime(startTime); // time is updated for 'elastic' only
260  }
261  G4HadProjectile theProj( track );
262  G4HadronicInteraction* hadi = nullptr;
263  G4HadFinalState* result = nullptr;
264 
265  // Select element
266  const G4Element* elm = nullptr;
267  G4int ZZ=1;
268 
269  try
270  {
271  elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
272  *targNucleus);
273  }
274  catch( G4HadronicException & aR )
275  {
277  aR.Report(ed);
278  DumpState(track,"SampleZandA",ed);
279  ed << " PostStepDoIt failed on element selection" << G4endl;
280  G4Exception("G4NeutrinoElectronProcess::PostStepDoIt", "had003",
281  FatalException, ed);
282  }
283  if( elm ) ZZ = elm->GetZ();
284 
285  G4double xsc = fTotXsc->GetElementCrossSection(dynParticle, ZZ, material);
286  xsc *= 1.;
287  G4double ccTotRatio = fTotXsc->GetCcRatio();
288 
289  if( G4UniformRand() < ccTotRatio ) // Cc-model
290  {
291  // Initialize the hadronic projectile from the track
292  thePro.Initialise(track);
293 
294  hadi = (GetHadronicInteractionList())[0];
295 
296  result = hadi->ApplyYourself( thePro, *targNucleus);
297 
299 
301 
302  FillResult(result, track);
303  }
304  else // Nc-model, like 'elastic', 2->2 scattering
305  {
306 
307  hadi = (GetHadronicInteractionList())[1];
308 
309  size_t idx = track.GetMaterialCutsCouple()->GetIndex();
310 
312 
313  hadi->SetRecoilEnergyThreshold(tcut);
314 
315  if( verboseLevel > 1 )
316  {
317  G4cout << "G4NeutrinoElectronProcess::PostStepDoIt for "
318  << part->GetParticleName()
319  << " in " << material->GetName()
320  << " Target Z= " << targNucleus->GetZ_asInt()
321  << " A= " << targNucleus->GetA_asInt() << G4endl;
322  }
323  try
324  {
325  result = hadi->ApplyYourself( theProj, *targNucleus);
326  }
327  catch(G4HadronicException & aR)
328  {
330  aR.Report(ed);
331  ed << "Call for " << hadi->GetModelName() << G4endl;
332  ed << "Target element "<< elm->GetName()<<" Z= "
333  << targNucleus->GetZ_asInt()
334  << " A= " << targNucleus->GetA_asInt() << G4endl;
335  DumpState(track,"ApplyYourself",ed);
336  ed << " ApplyYourself failed" << G4endl;
337  G4Exception("G4NeutrinoElectronProcess::PostStepDoIt", "had006",
338  FatalException, ed);
339  }
340  // directions
341 
342  G4ThreeVector indir = track.GetMomentumDirection();
344  G4ThreeVector it(0., 0., 1.);
345  G4ThreeVector outdir = result->GetMomentumChange();
346 
347  if(verboseLevel>1)
348  {
349  G4cout << "Efin= " << result->GetEnergyChange()
350  << " de= " << result->GetLocalEnergyDeposit()
351  << " nsec= " << result->GetNumberOfSecondaries()
352  << " dir= " << outdir
353  << G4endl;
354  }
355  // energies
356 
357  G4double edep = result->GetLocalEnergyDeposit();
358  G4double efinal = result->GetEnergyChange();
359 
360  if(efinal < 0.0) { efinal = 0.0; }
361  if(edep < 0.0) { edep = 0.0; }
362 
363  // NOTE: Very low energy scatters were causing numerical (FPE) errors
364  // in earlier releases; these limits have not been changed since.
365 
366  if(efinal <= lowestEnergy)
367  {
368  edep += efinal;
369  efinal = 0.0;
370  }
371  // primary change
372 
373  theTotalResult->ProposeEnergy(efinal);
374 
375  G4TrackStatus status = track.GetTrackStatus();
376 
377  if(efinal > 0.0)
378  {
379  outdir.rotate(phi, it);
380  outdir.rotateUz(indir);
382  }
383  else
384  {
385  if( part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
386  {
387  status = fStopButAlive;
388  }
389  else
390  {
391  status = fStopAndKill;
392  }
394  }
395  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
396 
398 
399  // recoil
400 
401  if( result->GetNumberOfSecondaries() > 0 )
402  {
403  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
404 
405  if(p->GetKineticEnergy() > tcut)
406  {
409 
410  // G4cout << "recoil " << pdir << G4endl;
412 
413  pdir.rotate(phi, it);
414  pdir.rotateUz(indir);
415 
416  // G4cout << "recoil rotated " << pdir << G4endl;
417 
418  p->SetMomentumDirection(pdir);
419 
420  // in elastic scattering time and weight are not changed
421 
422  G4Track* t = new G4Track(p, track.GetGlobalTime(),
423  track.GetPosition());
424  t->SetWeight(weight);
427  }
428  else
429  {
430  edep += p->GetKineticEnergy();
431  delete p;
432  }
433  }
436  result->Clear();
437  }
438  return theTotalResult;
439 }
440 
441 void
443 {
444  if(!isInitialised) {
445  isInitialised = true;
446  if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
447  }
449 }
450 
451 void
453 {
454  lowestEnergy = val;
455 }
456