ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VMultipleScattering.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VMultipleScattering.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 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4VMultipleScattering
33 //
34 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 //
36 // Creation date: 25.03.2003
37 //
38 // Modifications:
39 //
40 // 13.04.03 Change printout (V.Ivanchenko)
41 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
42 // 16-07-03 Use G4VMscModel interface (V.Ivanchenko)
43 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko)
44 // 04-11-03 Update PrintInfoDefinition (V.Ivanchenko)
45 // 01-03-04 SampleCosineTheta signature changed
46 // 22-04-04 SampleCosineTheta signature changed back to original
47 // 27-08-04 Add InitialiseForRun method (V.Ivanchneko)
48 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
49 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
50 // 15-04-05 optimize internal interface (V.Ivanchenko)
51 // 15-04-05 remove boundary flag (V.Ivanchenko)
52 // 27-10-05 introduce virtual function MscStepLimitation() (V.Ivanchenko)
53 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
54 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
55 // 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko)
56 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
57 // 04-06-13 Adoptation to MT mode (V.Ivanchenko)
58 //
59 // Class Description:
60 //
61 // It is the generic process of multiple scattering it includes common
62 // part of calculations for all charged particles
63 
64 // -------------------------------------------------------------------
65 //
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 
69 #include "G4VMultipleScattering.hh"
70 #include "G4PhysicalConstants.hh"
71 #include "G4SystemOfUnits.hh"
72 #include "G4LossTableManager.hh"
73 #include "G4MaterialCutsCouple.hh"
74 #include "G4Step.hh"
75 #include "G4ParticleDefinition.hh"
76 #include "G4VEmFluctuationModel.hh"
77 #include "G4UnitsTable.hh"
78 #include "G4ProductionCutsTable.hh"
79 #include "G4Electron.hh"
80 #include "G4GenericIon.hh"
82 #include "G4SafetyHelper.hh"
83 #include "G4ParticleTable.hh"
84 #include "G4ProcessVector.hh"
85 #include "G4ProcessManager.hh"
86 #include <iostream>
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
92  numberOfModels(0),
93  firstParticle(nullptr),
94  currParticle(nullptr),
95  stepLimit(fUseSafety),
96  facrange(0.04),
97  latDisplacement(true),
98  isIon(false),
99  fDispBeyondSafety(false),
100  fNewPosition(0.,0.,0.),
101  fNewDirection(0.,0.,1.)
102 {
104  SetVerboseLevel(1);
106  if("ionmsc" == name) { firstParticle = G4GenericIon::GenericIon(); }
107 
109 
111  fIonisation = nullptr;
112 
113  geomMin = 0.05*CLHEP::nm;
115 
117  safetyHelper = nullptr;
118  fPositionChanged = false;
119  isActive = false;
120 
121  currentModel = nullptr;
124  mscModels.reserve(2);
125  emManager->Register(this);
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
131 {
132  /*
133  if(1 < verboseLevel) {
134  G4cout << "G4VMultipleScattering destruct " << GetProcessName()
135  << G4endl;
136  }
137  */
138  delete modelManager;
139  emManager->DeRegister(this);
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143 
145  const G4Region* region)
146 {
147  if(!p) { return; }
148  G4VEmFluctuationModel* fm = nullptr;
149  modelManager->AddEmModel(order, p, fm, region);
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
156 {
157  for(auto & msc : mscModels) { if(msc == p) { return; } }
158  mscModels.push_back(p);
159 }
160 
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162 
164 {
165  return (index < mscModels.size()) ? mscModels[index] : nullptr;
166 }
167 
168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169 
170 void
172 {
173  if(1 < verboseLevel) {
174  G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
175  << GetProcessName()
176  << " and particle " << part.GetParticleName()
177  << G4endl;
178  }
179  G4bool master = emManager->IsMaster();
180 
181  if(!firstParticle) { firstParticle = &part; }
182  if(part.GetParticleType() == "nucleus") {
184  latDisplacement = false;
185  facrange = 0.2;
186  G4String pname = part.GetParticleName();
187  if(pname != "deuteron" && pname != "triton" &&
188  pname != "alpha+" && pname != "helium" &&
189  pname != "alpha" && pname != "He3" &&
190  pname != "hydrogen") {
191 
192  const G4ParticleDefinition* theGenericIon =
194  if(&part == theGenericIon) { isIon = true; }
195 
196  if(theGenericIon && firstParticle != theGenericIon) {
197  G4ProcessManager* pm = theGenericIon->GetProcessManager();
199  size_t n = v->size();
200  for(size_t j=0; j<n; ++j) {
201  if((*v)[j] == this) {
202  firstParticle = theGenericIon;
203  isIon = true;
204  break;
205  }
206  }
207  }
208  }
209  }
210 
211  emManager->PreparePhysicsTable(&part, this, master);
212  currParticle = nullptr;
213 
214  if(1 < verboseLevel) {
215  G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
216  << GetProcessName()
217  << " and particle " << part.GetParticleName()
218  << " local particle " << firstParticle->GetParticleName()
219  << " isIon: " << isIon << " isMaster: " << master
220  << G4endl;
221  }
222 
223  if(firstParticle == &part) {
224 
225  // initialise process
227 
228  // heavy particles and not ions
229  if(!isIon) {
230  if(part.GetPDGMass() > MeV) {
234  } else {
238  }
239  if(latDisplacement) {
241  }
242  }
243  if(master) { SetVerboseLevel(theParameters->Verbose()); }
245 
246  // initialisation of models
248  /*
249  G4cout << "### G4VMultipleScattering::PreparePhysicsTable() for "
250  << GetProcessName()
251  << " and particle " << part.GetParticleName()
252  << " Nmod= " << numberOfModels << " " << this
253  << G4endl;
254  */
255  for(G4int i=0; i<numberOfModels; ++i) {
256  G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i));
257  if(!msc) { continue; }
258  msc->SetIonisation(nullptr, firstParticle);
259  msc->SetMasterThread(master);
260  currentModel = msc;
262  G4double emax =
264  msc->SetHighEnergyLimit(emax);
265  }
266 
268  10.0, verboseLevel);
269 
270  if(!safetyHelper) {
272  ->GetSafetyHelper();
274  }
275  }
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
281 {
282  G4String num = part.GetParticleName();
283  G4bool master = emManager->IsMaster();
284  if(1 < verboseLevel) {
285  G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
286  << GetProcessName()
287  << " and particle " << num << " isIon: " << isIon
288  << " IsMaster: " << master << G4endl;
289  }
290  const G4VMultipleScattering* masterProcess =
291  static_cast<const G4VMultipleScattering*>(GetMasterProcess());
292 
293  if(firstParticle == &part) {
294  /*
295  G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
296  << GetProcessName()
297  << " and particle " << num
298  << " IsMaster= " << G4LossTableManager::Instance()->IsMaster()
299  << " " << this
300  << G4endl;
301  */
303 
304  if(!master) {
305  // initialisation of models
306  /*
307  G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
308  << GetProcessName()
309  << " and particle " << num
310  << " Nmod= " << numberOfModels << " " << this
311  << G4endl;
312  */
313  for(G4int i=0; i<numberOfModels; ++i) {
314  G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i));
315  if(!msc) { continue; }
316  G4VMscModel* msc0=
317  static_cast<G4VMscModel*>(masterProcess->GetModelByIndex(i));
318  msc->SetCrossSectionTable(msc0->GetCrossSectionTable(), false);
319  msc->InitialiseLocal(firstParticle, msc0);
320  }
321  }
322  }
323 
324  // explicitly defined printout by particle name
325  if(1 < verboseLevel ||
326  (0 < verboseLevel && (num == "e-" ||
327  num == "e+" || num == "mu+" ||
328  num == "mu-" || num == "proton"||
329  num == "pi+" || num == "pi-" ||
330  num == "kaon+" || num == "kaon-" ||
331  num == "alpha" || num == "anti_proton" ||
332  num == "GenericIon" || num == "alpha+" ||
333  num == "alpha++" )))
334  {
335  StreamInfo(G4cout, part);
336  }
337 
338  if(1 < verboseLevel) {
339  G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
340  << GetProcessName()
341  << " and particle " << num
342  << G4endl;
343  }
344 }
345 
346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
347 
348 void G4VMultipleScattering::StreamInfo(std::ostream& outFile,
349  const G4ParticleDefinition& part, G4bool rst) const
350 {
351  G4String indent = (rst ? " " : "");
352  outFile << G4endl << indent << GetProcessName() << ": ";
353  if (!rst) outFile << " for " << part.GetParticleName();
354  outFile << " SubType= " << GetProcessSubType() << G4endl;
355  StreamProcessInfo(outFile);
357 }
358 
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360 
362 {
363  G4VEnergyLossProcess* eloss = nullptr;
364  if(track->GetParticleDefinition() != currParticle) {
367  eloss = fIonisation;
368  }
369  /*
370  G4cout << "G4VMultipleScattering::StartTracking Nmod= " << numberOfModels
371  << " " << currParticle->GetParticleName()
372  << " E(MeV)= " << track->GetKineticEnergy()
373  << " Ion= " << eloss << " " << fIonisation << " IsMaster= "
374  << G4LossTableManager::Instance()->IsMaster()
375  << G4endl;
376  */
377  for(G4int i=0; i<numberOfModels; ++i) {
378  /*
379  G4cout << "Next model " << i << " " << msc
380  << " Emin= " << msc->LowEnergyLimit()
381  << " Emax= " << msc->HighEnergyLimit()
382  << " Eact= " << msc->LowEnergyActivationLimit() << G4endl;
383  */
384  G4VEmModel* msc = GetModelByIndex(i);
385  msc->StartTracking(track);
386  if(eloss) {
387  G4VMscModel* mscmod = static_cast<G4VMscModel*>(msc);
388  if(mscmod) { mscmod->SetIonisation(fIonisation, currParticle); }
389  }
390  }
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394 
396  const G4Track& track,
397  G4double,
398  G4double currentMinimalStep,
399  G4double&,
400  G4GPILSelection* selection)
401 {
402  // get Step limit proposed by the process
403  *selection = NotCandidateForSelection;
404  physStepLimit = gPathLength = tPathLength = currentMinimalStep;
405 
406  G4double ekin = track.GetKineticEnergy();
407  /*
408  G4cout << "MSC::AlongStepGPIL: Ekin= " << ekin
409  << " " << currParticle->GetParticleName()
410  << " currMod " << currentModel
411  << G4endl;
412  */
413  // isIon flag is used only to select a model
414  if(isIon) {
415  ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass();
416  }
417 
418  // select new model
419  if(1 < numberOfModels) {
420  currentModel = static_cast<G4VMscModel*>(
421  SelectModel(ekin,track.GetMaterialCutsCouple()->GetIndex()));
422  }
423  // msc is active is model is active, energy above the limit,
424  // and step size is above the limit;
425  // if it is active msc may limit the step
427  && ekin >= lowestKinEnergy) {
428  isActive = true;
429  tPathLength =
431  if (tPathLength < physStepLimit) {
432  *selection = CandidateForSelection;
433  }
434  } else { isActive = false; }
435 
436  //if(currParticle->GetPDGMass() > GeV)
437  /*
438  G4cout << "MSC::AlongStepGPIL: Ekin= " << ekin
439  << " " << currParticle->GetParticleName()
440  << " gPathLength= " << gPathLength
441  << " tPathLength= " << tPathLength
442  << " currentMinimalStep= " << currentMinimalStep
443  << " isActive " << isActive << G4endl;
444  */
445  return gPathLength;
446 }
447 
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
449 
450 G4double
453 {
454  *condition = NotForced;
455  return DBL_MAX;
456 }
457 
458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
459 
462 {
467  fPositionChanged = false;
468 
469  G4double geomLength = step.GetStepLength();
470 
471  // very small step - no msc
472  if(!isActive) {
473  tPathLength = geomLength;
474 
475  // sample msc
476  } else {
477  G4double range =
479  track.GetMaterialCutsCouple());
480 
482 
483  /*
484  if(currParticle->GetPDGMass() > 0.9*GeV)
485  G4cout << "G4VMsc::AlongStepDoIt: GeomLength= "
486  << geomLength
487  << " tPathLength= " << tPathLength
488  << " physStepLimit= " << physStepLimit
489  << " dr= " << range - tPathLength
490  << " ekin= " << track.GetKineticEnergy() << G4endl;
491  */
492  // protection against wrong t->g->t conversion
494 
495  // do not sample scattering at the last or at a small step
496  if(tPathLength < range && tPathLength > geomMin) {
497 
498  static const G4double minSafety = 1.20*CLHEP::nm;
499  static const G4double sFact = 0.99;
500 
502  step.GetPostStepPoint()->GetMomentumDirection(),minSafety);
503 
504  G4double r2 = displacement.mag2();
505  //G4cout << " R= " << sqrt(r2) << " Rmin= " << sqrt(minDisplacement2)
506  // << " flag= " << fDispBeyondSafety << G4endl;
507  if(r2 > minDisplacement2) {
508 
509  fPositionChanged = true;
510  G4double dispR = std::sqrt(r2);
511  G4double postSafety =
512  sFact*safetyHelper->ComputeSafety(fNewPosition, dispR);
513  //G4cout<<" R= "<< dispR<<" postSafety= "<<postSafety<<G4endl;
514 
515  // far away from geometry boundary
516  if(postSafety > 0.0 && dispR <= postSafety) {
517  fNewPosition += displacement;
518 
519  //near the boundary
520  } else {
521  // displaced point is definitely within the volume
522  //G4cout<<" R= "<<dispR<<" postSafety= "<<postSafety<<G4endl;
523  if(dispR < postSafety) {
524  fNewPosition += displacement;
525 
526  // optional extra mechanism is applied only if a particle
527  // is stopped by the boundary
528  } else if(fDispBeyondSafety && 0.0 == postSafety) {
529  fNewPosition += displacement;
530  G4double maxshift =
531  std::min(2.0*dispR, geomLength*(physStepLimit/tPathLength - 1.0));
532  G4double dist = 0.0;
533  G4double safety = postSafety + dispR;
535  /*
536  G4cout << "##MSC before Recheck maxshift= " << maxshift
537  << " postsafety= " << postSafety
538  << " Ekin= " << track.GetKineticEnergy()
539  << " " << track.GetDefinition()->GetParticleName()
540  << G4endl;
541  */
542  // check if it is possible to shift to the boundary
543  // and the shift is not large
545  fNewDirection, maxshift, &dist, &safety)
546  && std::abs(dist) < maxshift) {
547  /*
548  G4cout << "##MSC after Recheck dist= " << dist
549  << " postsafety= " << postSafety
550  << " t= " << tPathLength
551  << " g= " << geomLength
552  << " p= " << physStepLimit
553  << G4endl;
554  */
555  // shift is positive
556  if(dist >= 0.0) {
557  tPathLength *= (1.0 + dist/geomLength);
558  fNewPosition += dist*fNewDirection;
559 
560  // shift is negative cannot be larger than geomLength
561  } else {
562  maxshift = std::min(maxshift, geomLength);
563  if(0.0 < maxshift + dist) {
564  const G4ThreeVector& postpoint =
565  step.GetPostStepPoint()->GetPosition();
567  G4double R2 = (postpoint - point).mag2();
568  G4double newdist = dist;
569  // check not more than 10 extra boundaries
570  for(G4int i=0; i<10; ++i) {
571  dist = 0.0;
573  point, fNewDirection, maxshift, &dist, &safety)
574  && std::abs(newdist + dist) < maxshift) {
575  point += dist*fNewDirection;
576  G4double R2new = (postpoint - point).mag2();
577  //G4cout << "Backward i= " << i << " dist= " << dist
578  // << " R2= " << R2new << G4endl;
579  if(dist >= 0.0 || R2new > R2) { break; }
580  R2 = R2new;
582  newdist += dist;
583  } else {
584  break;
585  }
586  }
587  tPathLength *= (1.0 + newdist/geomLength);
588  // shift on boundary is not possible for negative disp
589  } else {
590  fNewPosition += displacement*(postSafety/dispR - 1.0);
591  }
592  }
593  // shift on boundary is not possible for any disp
594  } else {
595  fNewPosition += displacement*(postSafety/dispR - 1.0);
596  }
597  // reduced displacement
598  } else if(postSafety > geomMin) {
599  fNewPosition += displacement*(postSafety/dispR);
600 
601  // very small postSafety
602  } else {
603  fPositionChanged = false;
604  }
605  }
606  if(fPositionChanged) {
609  }
610  }
611  }
612  }
614  return &fParticleChange;
615 }
616 
617 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
618 
621 {
623  return &fParticleChange;
624 }
625 
626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
627 
629  const G4Track& track,
630  G4double previousStepSize,
631  G4double currentMinimalStep,
632  G4double& currentSafety)
633 {
635  G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize,
636  currentMinimalStep,
637  currentSafety,
638  &selection);
639  return x;
640 }
641 
642 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
643 
645  const G4Track& track,
646  G4double previousStepSize,
647  G4double currentMinimalStep,
648  G4double& currentSafety)
649 {
650  return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep,
651  currentSafety);
652 }
653 
654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
655 
658 {
659  *condition = Forced;
660  return DBL_MAX;
661 }
662 
663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
664 
665 G4bool
667  const G4String& directory,
668  G4bool ascii)
669 {
670  G4bool yes = true;
671  if(part != firstParticle) { return yes; }
672  const G4VMultipleScattering* masterProcess =
673  static_cast<const G4VMultipleScattering*>(GetMasterProcess());
674  if(masterProcess && masterProcess != this) { return yes; }
675 
677  static const G4String ss[4] = {"1","2","3","4"};
678  for(G4int i=0; i<nmod; ++i) {
679  G4VEmModel* msc = modelManager->GetModel(i);
680  yes = true;
681  G4PhysicsTable* table = msc->GetCrossSectionTable();
682  if (table) {
683  G4int j = std::min(i,3);
684  G4String name =
685  GetPhysicsTableFileName(part,directory,"LambdaMod"+ss[j],ascii);
686  yes = table->StorePhysicsTable(name,ascii);
687 
688  if ( yes ) {
689  if ( verboseLevel>0 ) {
690  G4cout << "Physics table are stored for "
691  << part->GetParticleName()
692  << " and process " << GetProcessName()
693  << " with a name <" << name << "> " << G4endl;
694  }
695  } else {
696  G4cout << "Fail to store Physics Table for "
697  << part->GetParticleName()
698  << " and process " << GetProcessName()
699  << " in the directory <" << directory
700  << "> " << G4endl;
701  }
702  }
703  }
704  return yes;
705 }
706 
707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
708 
709 G4bool
711  const G4String&,
712  G4bool)
713 {
714  return true;
715 }
716 
717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
718 
720 {
721  for(G4int i=0; i<numberOfModels; ++i) {
722  G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i, true));
723  if(msc) { msc->SetIonisation(p, firstParticle); }
724  }
725 }
726 
727 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
728 
729 void G4VMultipleScattering::ProcessDescription(std::ostream& outFile) const
730 {
731  if(firstParticle) {
732  StreamInfo(outFile, *firstParticle, true);
733  }
734 }
735 
736 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
737