ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LBE.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LBE.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 // For information related to this code contact: Alex Howard
30 // e-mail: alexander.howard@cern.ch
31 // --------------------------------------------------------------
32 // Comments
33 //
34 // Underground Advanced
35 //
36 // This physics list is taken from the underground_physics example with small
37 // modifications. It is an example of a "flat" physics list with no dependence
38 // on builders. The physics covered would be suitable for a low background
39 // experiment including the neutron_hp package
40 //
41 //
42 //
43 // PhysicsList program
44 //
45 // Modified:
46 //
47 // 14-02-03 Fix bugs in msc and hIon instanciation + cut per region
48 // 16-08-10 Remove inclusion of obsolete class of G4ParticleWithCuts
49 // 20-10-10 Migrate LowEnergy process to Livermore models, LP
50 // 28-03-13 Replace LEP/HEP with FTFP+BERT (A.R.)
51 // 15-01-20 Updated cross sections (A.R.)
52 // --------------------------------------------------------------
53 
54 #include <iomanip>
55 
56 #include "globals.hh"
57 #include <CLHEP/Units/SystemOfUnits.h>
58 #include "G4ios.hh"
59 #include "G4ProcessManager.hh"
60 #include "G4ProcessVector.hh"
61 
62 #include "G4ParticleTypes.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4ProductionCutsTable.hh"
65 
66 #include "G4UserLimits.hh"
67 #include "G4WarnPLStatus.hh"
68 
69 // Builder for all stopping processes
70 #include "G4StoppingPhysics.hh"
71 
72 #include "G4HadronicParameters.hh"
73 
74 #include "LBE.hh"
75 
76 // Constructor /////////////////////////////////////////////////////////////
77 LBE::LBE(G4int ver)
78 {
79 
80  G4cout << "You are using the simulation engine: LBE"<<G4endl;
81  G4cout <<G4endl<<G4endl;
82  defaultCutValue = 1.0*CLHEP::micrometer; //
83  cutForGamma = defaultCutValue;
84 // cutForElectron = 1.0*CLHEP::nanometer;
85  cutForElectron = 1.0*CLHEP::micrometer;
86  cutForPositron = defaultCutValue;
87  //not used:
88  // cutForProton = defaultCutValue;
89  // cutForAlpha = 1.0*CLHEP::nanometer;
90  // cutForGenericIon = 1.0*CLHEP::nanometer;
91 
92  stoppingPhysics = new G4StoppingPhysics;
93 
94  VerboseLevel = ver;
95  OpVerbLevel = 0;
96 
97  SetVerboseLevel(VerboseLevel);
98 }
99 
100 
101 // Destructor //////////////////////////////////////////////////////////////
102 LBE::~LBE()
103 {
104  delete stoppingPhysics;
105 }
106 
107 
108 // Construct Particles /////////////////////////////////////////////////////
109  void LBE::ConstructParticle()
110 {
111 
112  // In this method, static member functions should be called
113  // for all particles which you want to use.
114  // This ensures that objects of these particle types will be
115  // created in the program.
116 
117  ConstructMyBosons();
118  ConstructMyLeptons();
119  ConstructMyMesons();
120  ConstructMyBaryons();
121  ConstructMyIons();
122  ConstructMyShortLiveds();
123  stoppingPhysics->ConstructParticle(); // Anything not included above
124 }
125 
126 
127 // construct Bosons://///////////////////////////////////////////////////
128  void LBE::ConstructMyBosons()
129 {
130  // pseudo-particles
131  G4Geantino::GeantinoDefinition();
132  G4ChargedGeantino::ChargedGeantinoDefinition();
133 
134  // gamma
135  G4Gamma::GammaDefinition();
136 
137  //OpticalPhotons
138  G4OpticalPhoton::OpticalPhotonDefinition();
139 }
140 
141 
142 // construct Leptons://///////////////////////////////////////////////////
143  void LBE::ConstructMyLeptons()
144 {
145  // leptons
146  G4Electron::ElectronDefinition();
147  G4Positron::PositronDefinition();
148  G4MuonPlus::MuonPlusDefinition();
149  G4MuonMinus::MuonMinusDefinition();
150 
151  G4NeutrinoE::NeutrinoEDefinition();
152  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
153  G4NeutrinoMu::NeutrinoMuDefinition();
154  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
155 }
156 
157 #include "G4MesonConstructor.hh"
158 #include "G4BaryonConstructor.hh"
159 #include "G4IonConstructor.hh"
160 
161 
162 // construct Mesons://///////////////////////////////////////////////////
163  void LBE::ConstructMyMesons()
164 {
165  // mesons
166  G4MesonConstructor mConstructor;
167  mConstructor.ConstructParticle();
168 
169 }
170 
171 
172 // construct Baryons://///////////////////////////////////////////////////
173  void LBE::ConstructMyBaryons()
174 {
175  // baryons
176  G4BaryonConstructor bConstructor;
177  bConstructor.ConstructParticle();
178 
179 }
180 
181 
182 // construct Ions://///////////////////////////////////////////////////
183  void LBE::ConstructMyIons()
184 {
185  // ions
186  G4IonConstructor iConstructor;
187  iConstructor.ConstructParticle();
188 
189 }
190 
191 // construct Shortliveds://///////////////////////////////////////////////////
192  void LBE::ConstructMyShortLiveds()
193 {
194  // ShortLiveds
195  ;
196 }
197 
198 
199 
200 
201 // Construct Processes //////////////////////////////////////////////////////
202  void LBE::ConstructProcess()
203 {
204  AddTransportation();
205  ConstructEM();
206  ConstructOp();
207  ConstructHad();
208  ConstructGeneral();
209 }
210 
211 
212 // Transportation ///////////////////////////////////////////////////////////
213 #include "G4MaxTimeCuts.hh"
214 #include "G4MinEkineCuts.hh"
215 
216  void LBE::AddTransportation() {
217 
218  G4VUserPhysicsList::AddTransportation();
219 
220  auto myParticleIterator=G4ParticleTable::GetParticleTable()->GetIterator();
221  myParticleIterator->reset();
222  while( (*(myParticleIterator))() ){
223  G4ParticleDefinition* particle = myParticleIterator->value();
224  G4ProcessManager* pmanager = particle->GetProcessManager();
225  G4String particleName = particle->GetParticleName();
226  // time cuts for ONLY neutrons:
227  if(particleName == "neutron")
228  pmanager->AddDiscreteProcess(new G4MaxTimeCuts());
229  // Energy cuts to kill charged (embedded in method) particles:
230  pmanager->AddDiscreteProcess(new G4MinEkineCuts());
231  }
232 }
233 
234 
235 // Electromagnetic Processes ////////////////////////////////////////////////
236 // all charged particles
237 
238 #include "G4eMultipleScattering.hh"
239 #include "G4MuMultipleScattering.hh"
240 #include "G4hMultipleScattering.hh"
241 
242 // gamma. Use Livermore models
243 #include "G4PhotoElectricEffect.hh"
244 #include "G4LivermorePhotoElectricModel.hh"
245 #include "G4ComptonScattering.hh"
246 #include "G4LivermoreComptonModel.hh"
247 #include "G4GammaConversion.hh"
248 #include "G4LivermoreGammaConversionModel.hh"
249 #include "G4RayleighScattering.hh"
250 #include "G4LivermoreRayleighModel.hh"
251 
252 
253 // e-
254 #include "G4eMultipleScattering.hh"
255 #include "G4UniversalFluctuation.hh"
256 #include "G4UrbanMscModel.hh"
257 
258 #include "G4eIonisation.hh"
259 #include "G4LivermoreIonisationModel.hh"
260 
261 #include "G4eBremsstrahlung.hh"
262 #include "G4LivermoreBremsstrahlungModel.hh"
263 
264 // e+
265 #include "G4eplusAnnihilation.hh"
266 
267 
268 // alpha and GenericIon and deuterons, triton, He3:
269 #include "G4ionIonisation.hh"
270 #include "G4hIonisation.hh"
271 #include "G4hBremsstrahlung.hh"
272 //
273 #include "G4IonParametrisedLossModel.hh"
274 #include "G4NuclearStopping.hh"
275 #include "G4EnergyLossTables.hh"
276 
277 //muon:
278 #include "G4MuIonisation.hh"
279 #include "G4MuBremsstrahlung.hh"
280 #include "G4MuPairProduction.hh"
281 #include "G4MuonMinusCapture.hh"
282 
283 //OTHERS:
284 //#include "G4hIonisation.hh" // standard hadron ionisation
285 
286  void LBE::ConstructEM() {
287 
288  // models & processes:
289  // Use Livermore models up to 20 MeV, and standard
290  // models for higher energy
291  G4double LivermoreHighEnergyLimit = 20*CLHEP::MeV;
292  //
293  auto myParticleIterator=G4ParticleTable::GetParticleTable()->GetIterator();
294  myParticleIterator->reset();
295  while( (*(myParticleIterator))() ){
296  G4ParticleDefinition* particle = myParticleIterator->value();
297  G4ProcessManager* pmanager = particle->GetProcessManager();
298  G4String particleName = particle->GetParticleName();
299  G4String particleType = particle->GetParticleType();
300  G4double charge = particle->GetPDGCharge();
301 
302  if (particleName == "gamma")
303  {
304  G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
305  G4LivermorePhotoElectricModel* theLivermorePhotoElectricModel =
306  new G4LivermorePhotoElectricModel();
307  theLivermorePhotoElectricModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
308  thePhotoElectricEffect->AddEmModel(0, theLivermorePhotoElectricModel);
309  pmanager->AddDiscreteProcess(thePhotoElectricEffect);
310 
311  G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
312  G4LivermoreComptonModel* theLivermoreComptonModel =
313  new G4LivermoreComptonModel();
314  theLivermoreComptonModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
315  theComptonScattering->AddEmModel(0, theLivermoreComptonModel);
316  pmanager->AddDiscreteProcess(theComptonScattering);
317 
318  G4GammaConversion* theGammaConversion = new G4GammaConversion();
319  G4LivermoreGammaConversionModel* theLivermoreGammaConversionModel =
320  new G4LivermoreGammaConversionModel();
321  theLivermoreGammaConversionModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
322  theGammaConversion->AddEmModel(0, theLivermoreGammaConversionModel);
323  pmanager->AddDiscreteProcess(theGammaConversion);
324 
325  G4RayleighScattering* theRayleigh = new G4RayleighScattering();
326  G4LivermoreRayleighModel* theRayleighModel = new G4LivermoreRayleighModel();
327  theRayleighModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
328  theRayleigh->AddEmModel(0, theRayleighModel);
329  pmanager->AddDiscreteProcess(theRayleigh);
330 
331  }
332  else if (particleName == "e-")
333  {
334  //electron
335  // process ordering: AddProcess(name, at rest, along step, post step)
336  // -1 = not implemented, then ordering
337  G4eMultipleScattering* msc = new G4eMultipleScattering();
338  //msc->AddEmModel(0, new G4UrbanMscModel());
339  msc->SetStepLimitType(fUseDistanceToBoundary);
340  pmanager->AddProcess(msc, -1, 1, 1);
341 
342  // Ionisation
343  G4eIonisation* eIoni = new G4eIonisation();
344  G4LivermoreIonisationModel* theIoniLivermore = new
345  G4LivermoreIonisationModel();
346  theIoniLivermore->SetHighEnergyLimit(1*CLHEP::MeV);
347  eIoni->AddEmModel(0, theIoniLivermore, new G4UniversalFluctuation() );
348  eIoni->SetStepFunction(0.2, 100*CLHEP::um); //
349  pmanager->AddProcess(eIoni, -1, 2, 2);
350 
351  // Bremsstrahlung
352  G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
353  G4LivermoreBremsstrahlungModel* theBremLivermore = new
354  G4LivermoreBremsstrahlungModel();
355  theBremLivermore->SetHighEnergyLimit(LivermoreHighEnergyLimit);
356  eBrem->AddEmModel(0, theBremLivermore);
357  pmanager->AddProcess(eBrem, -1,-3, 3);
358  }
359  else if (particleName == "e+")
360  {
361  //positron
362  G4eMultipleScattering* msc = new G4eMultipleScattering();
363  //msc->AddEmModel(0, new G4UrbanMscModel());
364  msc->SetStepLimitType(fUseDistanceToBoundary);
365  pmanager->AddProcess(msc, -1, 1, 1);
366  G4eIonisation* eIoni = new G4eIonisation();
367  eIoni->SetStepFunction(0.2, 100*CLHEP::um);
368  pmanager->AddProcess(eIoni, -1, 2, 2);
369  pmanager->AddProcess(new G4eBremsstrahlung, -1,-3, 3);
370  pmanager->AddProcess(new G4eplusAnnihilation,0,-1, 4);
371  }
372  else if( particleName == "mu+" ||
373  particleName == "mu-" )
374  {
375  //muon
376  G4MuMultipleScattering* aMultipleScattering = new G4MuMultipleScattering();
377  pmanager->AddProcess(aMultipleScattering, -1, 1, 1);
378  pmanager->AddProcess(new G4MuIonisation(), -1, 2, 2);
379  pmanager->AddProcess(new G4MuBremsstrahlung(), -1,-1, 3);
380  pmanager->AddProcess(new G4MuPairProduction(), -1,-1, 4);
381  if( particleName == "mu-" )
382  pmanager->AddProcess(new G4MuonMinusCapture(), 0,-1,-1);
383  }
384  else if (particleName == "GenericIon")
385  {
386  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
387  G4ionIonisation* ionIoni = new G4ionIonisation();
388  ionIoni->SetEmModel(new G4IonParametrisedLossModel());
389  ionIoni->SetStepFunction(0.1, 10*CLHEP::um);
390  pmanager->AddProcess(ionIoni, -1, 2, 2);
391  pmanager->AddProcess(new G4NuclearStopping(), -1, 3,-1);
392  }
393  else if (particleName == "alpha" || particleName == "He3")
394  {
395  //MSC, ion-Ionisation, Nuclear Stopping
396  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
397 
398  G4ionIonisation* ionIoni = new G4ionIonisation();
399  ionIoni->SetStepFunction(0.1, 20*CLHEP::um);
400  pmanager->AddProcess(ionIoni, -1, 2, 2);
401  pmanager->AddProcess(new G4NuclearStopping(), -1, 3,-1);
402  }
403  else if (particleName == "proton" ||
404  particleName == "deuteron" ||
405  particleName == "triton" ||
406  particleName == "pi+" ||
407  particleName == "pi-" ||
408  particleName == "kaon+" ||
409  particleName == "kaon-")
410  {
411  //MSC, h-ionisation, bremsstrahlung
412  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
413  G4hIonisation* hIoni = new G4hIonisation();
414  hIoni->SetStepFunction(0.2, 50*CLHEP::um);
415  pmanager->AddProcess(hIoni, -1, 2, 2);
416  pmanager->AddProcess(new G4hBremsstrahlung, -1,-3, 3);
417  }
418  else if ((!particle->IsShortLived()) &&
419  (charge != 0.0) &&
420  (particle->GetParticleName() != "chargedgeantino"))
421  {
422  //all others charged particles except geantino
423  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
424  pmanager->AddProcess(new G4hIonisation, -1, 2, 2);
425  }
426 
427  }
428 }
429 
430 
431 // Optical Processes ////////////////////////////////////////////////////////
432 #include "G4Scintillation.hh"
433 #include "G4OpAbsorption.hh"
434 //#include "G4OpRayleigh.hh"
435 #include "G4OpBoundaryProcess.hh"
436 
437  void LBE::ConstructOp()
438 {
439  // default scintillation process
440  //Coverity report: check that the process is actually used, if not must delete
441  G4bool theScintProcessDefNeverUsed = true;
442  G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
443  // theScintProcessDef->DumpPhysicsTable();
444  theScintProcessDef->SetTrackSecondariesFirst(true);
445  theScintProcessDef->SetScintillationYieldFactor(1.0); //
446  theScintProcessDef->SetScintillationExcitationRatio(0.0); //
447  theScintProcessDef->SetVerboseLevel(OpVerbLevel);
448 
449  // scintillation process for alpha:
450  G4bool theScintProcessAlphaNeverUsed = true;
451  G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
452  // theScintProcessNuc->DumpPhysicsTable();
453  theScintProcessAlpha->SetTrackSecondariesFirst(true);
454  theScintProcessAlpha->SetScintillationYieldFactor(1.1);
455  theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
456  theScintProcessAlpha->SetVerboseLevel(OpVerbLevel);
457 
458  // scintillation process for heavy nuclei
459  G4bool theScintProcessNucNeverUsed = true;
460  G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
461  // theScintProcessNuc->DumpPhysicsTable();
462  theScintProcessNuc->SetTrackSecondariesFirst(true);
463  theScintProcessNuc->SetScintillationYieldFactor(0.2);
464  theScintProcessNuc->SetScintillationExcitationRatio(1.0);
465  theScintProcessNuc->SetVerboseLevel(OpVerbLevel);
466 
467  // optical processes
468  G4bool theAbsorptionProcessNeverUsed = true;
469  G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
470  // G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh();
471  G4bool theBoundaryProcessNeverUsed = true;
472  G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
473  // theAbsorptionProcess->DumpPhysicsTable();
474  // theRayleighScatteringProcess->DumpPhysicsTable();
475  theAbsorptionProcess->SetVerboseLevel(OpVerbLevel);
476  // theRayleighScatteringProcess->SetVerboseLevel(OpVerbLevel);
477  theBoundaryProcess->SetVerboseLevel(OpVerbLevel);
478 
479  auto myParticleIterator=G4ParticleTable::GetParticleTable()->GetIterator();
480  myParticleIterator->reset();
481  while( (*(myParticleIterator))() )
482  {
483  G4ParticleDefinition* particle = myParticleIterator->value();
484  G4ProcessManager* pmanager = particle->GetProcessManager();
485  G4String particleName = particle->GetParticleName();
486  if (theScintProcessDef->IsApplicable(*particle)) {
487  // if(particle->GetPDGMass() > 5.0*CLHEP::GeV)
488  if(particle->GetParticleName() == "GenericIon") {
489  pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete
490  pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest);
491  pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep);
492  theScintProcessNucNeverUsed = false;
493  }
494  else if(particle->GetParticleName() == "alpha") {
495  pmanager->AddProcess(theScintProcessAlpha);
496  pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest);
497  pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep);
498  theScintProcessAlphaNeverUsed = false;
499  }
500  else {
501  pmanager->AddProcess(theScintProcessDef);
502  pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest);
503  pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep);
504  theScintProcessDefNeverUsed = false;
505  }
506  }
507 
508  if (particleName == "opticalphoton") {
509  pmanager->AddDiscreteProcess(theAbsorptionProcess);
510  theAbsorptionProcessNeverUsed = false;
511  // pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
512  theBoundaryProcessNeverUsed = false;
513  pmanager->AddDiscreteProcess(theBoundaryProcess);
514  }
515  }
516  if ( theScintProcessDefNeverUsed ) delete theScintProcessDef;
517  if ( theScintProcessAlphaNeverUsed ) delete theScintProcessAlpha;
518  if ( theScintProcessNucNeverUsed ) delete theScintProcessNuc;
519  if ( theBoundaryProcessNeverUsed ) delete theBoundaryProcess;
520  if ( theAbsorptionProcessNeverUsed ) delete theAbsorptionProcess;
521 }
522 
523 
524 // Hadronic processes ////////////////////////////////////////////////////////
525 
526 // Elastic processes:
527 #include "G4HadronElasticProcess.hh"
528 #include "G4HadronCaptureProcess.hh"
529 #include "G4HadronElastic.hh"
530 #include "G4ChipsElasticModel.hh"
531 #include "G4ElasticHadrNucleusHE.hh"
532 #include "G4AntiNuclElastic.hh"
533 #include "G4BGGPionElasticXS.hh"
534 #include "G4CrossSectionDataSetRegistry.hh"
535 #include "G4ChipsProtonElasticXS.hh"
536 #include "G4ChipsNeutronElasticXS.hh"
537 #include "G4ComponentAntiNuclNuclearXS.hh"
538 #include "G4ChipsKaonMinusElasticXS.hh"
539 #include "G4ChipsKaonPlusElasticXS.hh"
540 #include "G4ChipsKaonZeroElasticXS.hh"
541 #include "G4BGGNucleonElasticXS.hh"
542 #include "G4CrossSectionElastic.hh"
543 
544 // Inelastic processes:
545 #include "G4PionPlusInelasticProcess.hh"
546 #include "G4PionMinusInelasticProcess.hh"
547 #include "G4KaonPlusInelasticProcess.hh"
548 #include "G4KaonZeroSInelasticProcess.hh"
549 #include "G4KaonZeroLInelasticProcess.hh"
550 #include "G4KaonMinusInelasticProcess.hh"
551 #include "G4ProtonInelasticProcess.hh"
552 #include "G4AntiProtonInelasticProcess.hh"
553 #include "G4NeutronInelasticProcess.hh"
554 #include "G4AntiNeutronInelasticProcess.hh"
555 #include "G4DeuteronInelasticProcess.hh"
556 #include "G4TritonInelasticProcess.hh"
557 #include "G4AlphaInelasticProcess.hh"
558 
559 // FTFP + BERT model
560 #include "G4TheoFSGenerator.hh"
561 #include "G4ExcitationHandler.hh"
562 #include "G4PreCompoundModel.hh"
563 #include "G4GeneratorPrecompoundInterface.hh"
564 #include "G4FTFModel.hh"
565 #include "G4LundStringFragmentation.hh"
566 #include "G4ExcitedStringDecay.hh"
567 #include "G4CascadeInterface.hh"
568 #include "G4CrossSectionInelastic.hh"
569 #include "G4BGGPionInelasticXS.hh"
570 #include "G4ChipsKaonMinusInelasticXS.hh"
571 #include "G4ChipsKaonPlusInelasticXS.hh"
572 #include "G4ChipsKaonZeroInelasticXS.hh"
573 #include "G4CrossSectionDataSetRegistry.hh"
574 #include "G4BGGNucleonInelasticXS.hh"
575 #include "G4ComponentAntiNuclNuclearXS.hh"
576 #include "G4ComponentGGNuclNuclXsc.hh"
577 
578 // Neutron high-precision models: <20 MeV
579 #include "G4ParticleHPElastic.hh"
580 #include "G4ParticleHPElasticData.hh"
581 #include "G4ParticleHPCapture.hh"
582 #include "G4ParticleHPCaptureData.hh"
583 #include "G4ParticleHPInelastic.hh"
584 #include "G4ParticleHPInelasticData.hh"
585 #include "G4NeutronCaptureXS.hh"
586 #include "G4NeutronRadCapture.hh"
587 
588 // Binary light ion cascade for alpha, deuteron and triton
589 #include "G4BinaryLightIonReaction.hh"
590 
591 // ConstructHad()
592 // Makes discrete physics processes for the hadrons, at present limited
593 // to those particles with GHEISHA interactions (INTRC > 0).
594 // The processes are: Elastic scattering and Inelastic scattering.
595 // F.W.Jones 09-JUL-1998
596  void LBE::ConstructHad()
597 {
598  // Elastic scattering
599  G4HadronElastic* elastic_lhep0 = new G4HadronElastic();
600  G4ChipsElasticModel* elastic_chip = new G4ChipsElasticModel();
601  G4ElasticHadrNucleusHE* elastic_he = new G4ElasticHadrNucleusHE();
602 
603  const G4double elastic_elimitAntiNuc = 100.0*CLHEP::MeV;
604  G4AntiNuclElastic* elastic_anuc = new G4AntiNuclElastic();
605  elastic_anuc->SetMinEnergy( elastic_elimitAntiNuc );
606  G4CrossSectionElastic* elastic_anucxs = new G4CrossSectionElastic( elastic_anuc->GetComponentCrossSection() );
607  G4HadronElastic* elastic_lhep2 = new G4HadronElastic();
608  elastic_lhep2->SetMaxEnergy( elastic_elimitAntiNuc );
609 
610  // Inelastic scattering
611  const G4double theFTFMin0 = 0.0*CLHEP::GeV;
612  const G4double theFTFMin1 = 4.0*CLHEP::GeV;
613  const G4double theFTFMax = G4HadronicParameters::Instance()->GetMaxEnergy();
614  const G4double theBERTMin0 = 0.0*CLHEP::GeV;
615  const G4double theBERTMin1 = 19.0*CLHEP::MeV;
616  const G4double theBERTMax = 5.0*CLHEP::GeV;
617  const G4double theHPMin = 0.0*CLHEP::GeV;
618  const G4double theHPMax = 20.0*CLHEP::MeV;
619  const G4double theIonBCMin = 0.0*CLHEP::GeV;
620  const G4double theIonBCMax = 5.0*CLHEP::GeV;
621 
622 
623  G4FTFModel * theStringModel = new G4FTFModel;
624  G4ExcitedStringDecay * theStringDecay = new G4ExcitedStringDecay( new G4LundStringFragmentation );
625  theStringModel->SetFragmentationModel( theStringDecay );
626  G4PreCompoundModel * thePreEquilib = new G4PreCompoundModel( new G4ExcitationHandler );
627  G4GeneratorPrecompoundInterface * theCascade = new G4GeneratorPrecompoundInterface( thePreEquilib );
628 
629  G4TheoFSGenerator * theFTFModel0 = new G4TheoFSGenerator( "FTFP" );
630  theFTFModel0->SetHighEnergyGenerator( theStringModel );
631  theFTFModel0->SetTransport( theCascade );
632  theFTFModel0->SetMinEnergy( theFTFMin0 );
633  theFTFModel0->SetMaxEnergy( theFTFMax );
634 
635  G4TheoFSGenerator * theFTFModel1 = new G4TheoFSGenerator( "FTFP" );
636  theFTFModel1->SetHighEnergyGenerator( theStringModel );
637  theFTFModel1->SetTransport( theCascade );
638  theFTFModel1->SetMinEnergy( theFTFMin1 );
639  theFTFModel1->SetMaxEnergy( theFTFMax );
640 
641  G4CascadeInterface * theBERTModel0 = new G4CascadeInterface;
642  theBERTModel0->SetMinEnergy( theBERTMin0 );
643  theBERTModel0->SetMaxEnergy( theBERTMax );
644 
645  G4CascadeInterface * theBERTModel1 = new G4CascadeInterface;
646  theBERTModel1->SetMinEnergy( theBERTMin1 );
647  theBERTModel1->SetMaxEnergy( theBERTMax );
648 
649  // Binary Cascade
650  G4BinaryLightIonReaction * theIonBC = new G4BinaryLightIonReaction( thePreEquilib );
651  theIonBC->SetMinEnergy( theIonBCMin );
652  theIonBC->SetMaxEnergy( theIonBCMax );
653 
654  G4VCrossSectionDataSet * theAntiNucleonData = new G4CrossSectionInelastic( new G4ComponentAntiNuclNuclearXS );
655  G4ComponentGGNuclNuclXsc * ggNuclNuclXsec = new G4ComponentGGNuclNuclXsc();
656  G4VCrossSectionDataSet * theGGNuclNuclData = new G4CrossSectionInelastic(ggNuclNuclXsec);
657 
658  auto myParticleIterator=G4ParticleTable::GetParticleTable()->GetIterator();
659  myParticleIterator->reset();
660  while ((*(myParticleIterator))())
661  {
662  G4ParticleDefinition* particle = myParticleIterator->value();
663  G4ProcessManager* pmanager = particle->GetProcessManager();
664  G4String particleName = particle->GetParticleName();
665 
666  if (particleName == "pi+")
667  {
668  // Elastic scattering
669  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
670  theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
671  theElasticProcess->RegisterMe( elastic_he );
672  pmanager->AddDiscreteProcess( theElasticProcess );
673  // Inelastic scattering
674  G4PionPlusInelasticProcess* theInelasticProcess = new G4PionPlusInelasticProcess("inelastic");
675  theInelasticProcess->AddDataSet( new G4BGGPionInelasticXS( G4PionPlus::Definition() ) );
676  theInelasticProcess->RegisterMe( theFTFModel1 );
677  theInelasticProcess->RegisterMe( theBERTModel0 );
678  pmanager->AddDiscreteProcess( theInelasticProcess );
679  }
680 
681  else if (particleName == "pi-")
682  {
683  // Elastic scattering
684  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
685  theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
686  theElasticProcess->RegisterMe( elastic_he );
687  pmanager->AddDiscreteProcess( theElasticProcess );
688  // Inelastic scattering
689  G4PionMinusInelasticProcess* theInelasticProcess = new G4PionMinusInelasticProcess("inelastic");
690  theInelasticProcess->AddDataSet( new G4BGGPionInelasticXS( G4PionMinus::Definition() ) );
691  theInelasticProcess->RegisterMe( theFTFModel1 );
692  theInelasticProcess->RegisterMe( theBERTModel0 );
693  pmanager->AddDiscreteProcess( theInelasticProcess );
694  }
695 
696  else if (particleName == "kaon+")
697  {
698  // Elastic scattering
699  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
700  theElasticProcess->RegisterMe( elastic_lhep0 );
701  theElasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonPlusElasticXS::Default_Name()));
702  pmanager->AddDiscreteProcess( theElasticProcess );
703  // Inelastic scattering
704  G4KaonPlusInelasticProcess* theInelasticProcess = new G4KaonPlusInelasticProcess("inelastic");
705  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonPlusInelasticXS::Default_Name()));
706  theInelasticProcess->RegisterMe( theFTFModel1 );
707  theInelasticProcess->RegisterMe( theBERTModel0 );
708  pmanager->AddDiscreteProcess( theInelasticProcess );
709  }
710 
711  else if (particleName == "kaon0S")
712  {
713  // Elastic scattering
714  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
715  theElasticProcess->RegisterMe( elastic_lhep0 );
716  theElasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroElasticXS::Default_Name()));
717  pmanager->AddDiscreteProcess( theElasticProcess );
718  // Inelastic scattering
719  G4KaonZeroSInelasticProcess* theInelasticProcess = new G4KaonZeroSInelasticProcess("inelastic");
720  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroInelasticXS::Default_Name()));
721  theInelasticProcess->RegisterMe( theFTFModel1 );
722  theInelasticProcess->RegisterMe( theBERTModel0 );
723  pmanager->AddDiscreteProcess( theInelasticProcess );
724  }
725 
726  else if (particleName == "kaon0L")
727  {
728  // Elastic scattering
729  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
730  theElasticProcess->RegisterMe( elastic_lhep0 );
731  theElasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroElasticXS::Default_Name()));
732  pmanager->AddDiscreteProcess( theElasticProcess );
733  // Inelastic scattering
734  G4KaonZeroLInelasticProcess* theInelasticProcess = new G4KaonZeroLInelasticProcess("inelastic");
735  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroInelasticXS::Default_Name()));
736  theInelasticProcess->RegisterMe( theFTFModel1 );
737  theInelasticProcess->RegisterMe( theBERTModel0 );
738  pmanager->AddDiscreteProcess( theInelasticProcess );
739  }
740 
741  else if (particleName == "kaon-")
742  {
743  // Elastic scattering
744  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
745  theElasticProcess->RegisterMe( elastic_lhep0 );
746  theElasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonMinusElasticXS::Default_Name()));
747  pmanager->AddDiscreteProcess( theElasticProcess );
748  // Inelastic scattering
749  G4KaonMinusInelasticProcess* theInelasticProcess = new G4KaonMinusInelasticProcess("inelastic");
750  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonMinusInelasticXS::Default_Name()));
751  theInelasticProcess->RegisterMe( theFTFModel1 );
752  theInelasticProcess->RegisterMe( theBERTModel0 );
753  pmanager->AddDiscreteProcess( theInelasticProcess );
754  }
755 
756  else if (particleName == "proton")
757  {
758  // Elastic scattering
759  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
760  theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name()));
761  theElasticProcess->AddDataSet( new G4BGGNucleonElasticXS( G4Proton::Proton() ) );
762  theElasticProcess->RegisterMe( elastic_chip );
763  pmanager->AddDiscreteProcess( theElasticProcess );
764  // Inelastic scattering
765  G4ProtonInelasticProcess* theInelasticProcess = new G4ProtonInelasticProcess("inelastic");
766  theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Proton::Proton() ) );
767  theInelasticProcess->RegisterMe( theFTFModel1 );
768  theInelasticProcess->RegisterMe( theBERTModel0 );
769  pmanager->AddDiscreteProcess( theInelasticProcess );
770  }
771 
772  else if (particleName == "anti_proton")
773  {
774  // Elastic scattering
775  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
776  theElasticProcess->AddDataSet( elastic_anucxs );
777  theElasticProcess->RegisterMe( elastic_lhep2 );
778  theElasticProcess->RegisterMe( elastic_anuc );
779  pmanager->AddDiscreteProcess( theElasticProcess );
780  // Inelastic scattering
781  G4AntiProtonInelasticProcess* theInelasticProcess = new G4AntiProtonInelasticProcess("inelastic");
782  theInelasticProcess->AddDataSet( theAntiNucleonData );
783  theInelasticProcess->RegisterMe( theFTFModel0 );
784  pmanager->AddDiscreteProcess( theInelasticProcess );
785  }
786 
787  else if (particleName == "neutron") {
788  // elastic scattering
789  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
790  theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name()));
791  G4HadronElastic* elastic_neutronChipsModel = new G4ChipsElasticModel();
792  elastic_neutronChipsModel->SetMinEnergy( 19.0*CLHEP::MeV );
793  theElasticProcess->RegisterMe( elastic_neutronChipsModel );
794  G4ParticleHPElastic * theElasticNeutronHP = new G4ParticleHPElastic;
795  theElasticNeutronHP->SetMinEnergy( theHPMin );
796  theElasticNeutronHP->SetMaxEnergy( theHPMax );
797  theElasticProcess->RegisterMe( theElasticNeutronHP );
798  theElasticProcess->AddDataSet( new G4ParticleHPElasticData );
799  pmanager->AddDiscreteProcess( theElasticProcess );
800  // inelastic scattering
801  G4NeutronInelasticProcess* theInelasticProcess = new G4NeutronInelasticProcess("inelastic");
802  theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Neutron::Neutron() ) );
803  theInelasticProcess->RegisterMe( theFTFModel1 );
804  theInelasticProcess->RegisterMe( theBERTModel1 );
805  G4ParticleHPInelastic * theNeutronInelasticHPModel = new G4ParticleHPInelastic;
806  theNeutronInelasticHPModel->SetMinEnergy( theHPMin );
807  theNeutronInelasticHPModel->SetMaxEnergy( theHPMax );
808  theInelasticProcess->RegisterMe( theNeutronInelasticHPModel );
809  theInelasticProcess->AddDataSet( new G4ParticleHPInelasticData );
810  pmanager->AddDiscreteProcess(theInelasticProcess);
811  // capture
812  G4HadronCaptureProcess* theCaptureProcess = new G4HadronCaptureProcess;
813  G4ParticleHPCapture * theNeutronCaptureHPModel = new G4ParticleHPCapture;
814  theNeutronCaptureHPModel->SetMinEnergy( theHPMin );
815  theNeutronCaptureHPModel->SetMaxEnergy( theHPMax );
816  G4NeutronRadCapture* theNeutronRadCapture = new G4NeutronRadCapture();
817  theNeutronRadCapture->SetMinEnergy(theHPMax*0.99);
818  theCaptureProcess->RegisterMe( theNeutronCaptureHPModel );
819  theCaptureProcess->RegisterMe( theNeutronRadCapture);
820  theCaptureProcess->AddDataSet( new G4ParticleHPCaptureData );
821  theCaptureProcess->AddDataSet((G4NeutronCaptureXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4NeutronCaptureXS::Default_Name()));
822  pmanager->AddDiscreteProcess(theCaptureProcess);
823  }
824  else if (particleName == "anti_neutron")
825  {
826  // Elastic scattering
827  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
828  theElasticProcess->AddDataSet( elastic_anucxs );
829  theElasticProcess->RegisterMe( elastic_lhep2 );
830  theElasticProcess->RegisterMe( elastic_anuc );
831  pmanager->AddDiscreteProcess( theElasticProcess );
832  // Inelastic scattering
833  G4AntiNeutronInelasticProcess* theInelasticProcess = new G4AntiNeutronInelasticProcess("inelastic");
834  theInelasticProcess->AddDataSet( theAntiNucleonData );
835  theInelasticProcess->RegisterMe( theFTFModel0 );
836  pmanager->AddDiscreteProcess( theInelasticProcess );
837  }
838 
839  else if (particleName == "deuteron")
840  {
841  // Elastic scattering
842  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
843  theElasticProcess->RegisterMe( elastic_lhep0 );
844  theElasticProcess->AddDataSet( theGGNuclNuclData );
845  pmanager->AddDiscreteProcess( theElasticProcess );
846  // Inelastic scattering
847  G4DeuteronInelasticProcess* theInelasticProcess = new G4DeuteronInelasticProcess("inelastic");
848  theInelasticProcess->AddDataSet( theGGNuclNuclData );
849  theInelasticProcess->RegisterMe( theFTFModel1 );
850  theInelasticProcess->RegisterMe( theIonBC );
851  pmanager->AddDiscreteProcess( theInelasticProcess );
852  }
853 
854  else if (particleName == "triton")
855  {
856  // Elastic scattering
857  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
858  theElasticProcess->RegisterMe( elastic_lhep0 );
859  theElasticProcess->AddDataSet( theGGNuclNuclData );
860  pmanager->AddDiscreteProcess( theElasticProcess );
861  // Inelastic scattering
862  G4TritonInelasticProcess* theInelasticProcess = new G4TritonInelasticProcess("inelastic");
863  theInelasticProcess->AddDataSet( theGGNuclNuclData );
864  theInelasticProcess->RegisterMe( theFTFModel1 );
865  theInelasticProcess->RegisterMe( theIonBC );
866  pmanager->AddDiscreteProcess( theInelasticProcess );
867  }
868 
869  else if (particleName == "alpha")
870  {
871  // Elastic scattering
872  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
873  theElasticProcess->RegisterMe( elastic_lhep0 );
874  theElasticProcess->AddDataSet( theGGNuclNuclData );
875  pmanager->AddDiscreteProcess( theElasticProcess );
876  // Inelastic scattering
877  G4AlphaInelasticProcess* theInelasticProcess = new G4AlphaInelasticProcess("inelastic");
878  theInelasticProcess->AddDataSet( theGGNuclNuclData );
879  theInelasticProcess->RegisterMe( theFTFModel1 );
880  theInelasticProcess->RegisterMe( theIonBC );
881  pmanager->AddDiscreteProcess( theInelasticProcess );
882  }
883  } // while ((*(myParticleIterator))())
884 
885  // Add stopping processes with builder
886  stoppingPhysics->ConstructProcess();
887 }
888 
889 
890 // Decays ///////////////////////////////////////////////////////////////////
891 #include "G4Decay.hh"
892 #include "G4RadioactiveDecayBase.hh"
893 #include "G4IonTable.hh"
894 #include "G4Ions.hh"
895 
896 #include "G4LossTableManager.hh"
897 #include "G4UAtomicDeexcitation.hh"
898 #include "G4NuclearLevelData.hh"
899 #include "G4NuclideTable.hh"
900 
901  void LBE::ConstructGeneral() {
902 
903  // Add Decay Process
904  G4Decay* theDecayProcess = new G4Decay();
905  G4bool theDecayProcessNeverUsed = true; //Check if theDecayProcess will be used
906  auto myParticleIterator=G4ParticleTable::GetParticleTable()->GetIterator();
907  myParticleIterator->reset();
908  while( (*(myParticleIterator))() )
909  {
910  G4ParticleDefinition* particle = myParticleIterator->value();
911  G4ProcessManager* pmanager = particle->GetProcessManager();
912 
913  if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
914  {
915  theDecayProcessNeverUsed = false;
916  pmanager ->AddProcess(theDecayProcess);
917  // set ordering for PostStepDoIt and AtRestDoIt
918  pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
919  pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
920  }
921  }
922 
923  // Declare radioactive decay to the GenericIon in the IonTable.
924  const G4IonTable *theIonTable =
925  G4ParticleTable::GetParticleTable()->GetIonTable();
926  G4RadioactiveDecayBase* theRadioactiveDecay = new G4RadioactiveDecayBase();
927 
928  //Fix for activation of RadioactiveDecay, based on G4RadioactiveDecayPhysics
929  G4EmParameters* param = G4EmParameters::Instance();
930  param->SetAugerCascade(true);
931  param->AddPhysics("world","G4RadioactiveDecay");
932 
933  G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
934  deex->SetStoreAllLevels(true);
935  deex->SetMaxLifeTime(G4NuclideTable::GetInstance()->GetThresholdOfHalfLife()
936  /std::log(2.));
937 
938  G4LossTableManager* man = G4LossTableManager::Instance();
939  G4VAtomDeexcitation* ad = man->AtomDeexcitation();
940  if(!ad) {
941  ad = new G4UAtomicDeexcitation();
942  man->SetAtomDeexcitation(ad);
943  ad->InitialiseAtomicDeexcitation();
944  }
945 
946  for (G4int i=0; i<theIonTable->Entries(); i++)
947  {
948  G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
949  G4String particleType = theIonTable->GetParticle(i)->GetParticleType();
950 
951  if (particleName == "GenericIon")
952  {
953  G4ProcessManager* pmanager =
954  theIonTable->GetParticle(i)->GetProcessManager();
955  pmanager->SetVerboseLevel(VerboseLevel);
956  pmanager ->AddProcess(theRadioactiveDecay);
957  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
958  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
959  }
960  }
961  //If we actually never used the process, delete it
962  //From Coverity report
963  if ( theDecayProcessNeverUsed ) delete theDecayProcess;
964 }
965 
966 // Cuts /////////////////////////////////////////////////////////////////////
967 void LBE::SetCuts()
968 {
969 
970  if (verboseLevel >1)
971  G4cout << "LBE::SetCuts:";
972 
973  if (verboseLevel>0){
974  G4cout << "LBE::SetCuts:";
975  G4cout << "CutLength : "
976  << G4BestUnit(defaultCutValue,"Length") << G4endl;
977  }
978 
979  //special for low energy physics
980  G4double lowlimit=250*CLHEP::eV;
981  G4ProductionCutsTable * aPCTable = G4ProductionCutsTable::GetProductionCutsTable();
982  aPCTable->SetEnergyRange(lowlimit,100*CLHEP::GeV);
983 
984  // set cut values for gamma at first and for e- second and next for e+,
985  // because some processes for e+/e- need cut values for gamma
986  SetCutValue(cutForGamma, "gamma");
987  SetCutValue(cutForElectron, "e-");
988  SetCutValue(cutForPositron, "e+");
989 
990  // SetCutValue(cutForProton, "proton");
991  // SetCutValue(cutForProton, "anti_proton");
992  // SetCutValue(cutForAlpha, "alpha");
993  // SetCutValue(cutForGenericIon, "GenericIon");
994 
995  // SetCutValueForOthers(defaultCutValue);
996 
997  if (verboseLevel>0) DumpCutValuesTable();
998 }
999