| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 << 26 //
27 #include "G4VAdjointReverseReaction.hh" 27 #include "G4VAdjointReverseReaction.hh"
28 << 28 #include "G4SystemOfUnits.hh"
29 #include "G4AdjointCSManager.hh" 29 #include "G4AdjointCSManager.hh"
30 #include "G4ParticleChange.hh" << 30 #include "G4AdjointCSMatrix.hh"
>> 31 #include "G4AdjointInterpolator.hh"
>> 32 #include "G4AdjointCSMatrix.hh"
31 #include "G4VEmAdjointModel.hh" 33 #include "G4VEmAdjointModel.hh"
>> 34 #include "G4ElementTable.hh"
>> 35 #include "G4Element.hh"
>> 36 #include "G4Material.hh"
>> 37 #include "G4MaterialCutsCouple.hh"
>> 38 #include "G4AdjointCSManager.hh"
>> 39 #include "G4ParticleChange.hh"
>> 40 #include "G4AdjointElectron.hh"
32 41
33 G4VAdjointReverseReaction::G4VAdjointReverseRe <<
34 <<
35 : G4VDiscreteProcess(process_name) <<
36 { <<
37 fCSManager = G4AdjointCSManager::GetA <<
38 fIsScatProjToProj = whichScatCase; <<
39 fParticleChange = new G4ParticleChange(); <<
40 } <<
41 42
>> 43 G4VAdjointReverseReaction::
>> 44 G4VAdjointReverseReaction(G4String process_name, G4bool whichScatCase):
>> 45 G4VDiscreteProcess(process_name)
>> 46 {theAdjointCSManager = G4AdjointCSManager::GetAdjointCSManager();
>> 47 IsScatProjToProjCase=whichScatCase;
>> 48 fParticleChange=new G4ParticleChange();
>> 49 IsFwdCSUsed=false;
>> 50 IsIntegralModeUsed=false;
>> 51 lastCS=0.;
>> 52 trackid = nstep = 0;
>> 53 }
42 ////////////////////////////////////////////// 54 //////////////////////////////////////////////////////////////////////////////
43 G4VAdjointReverseReaction::~G4VAdjointReverseR << 55 //
44 { << 56 G4VAdjointReverseReaction::
45 if(fParticleChange) << 57 ~G4VAdjointReverseReaction()
46 delete fParticleChange; << 58 { if (fParticleChange) delete fParticleChange;
>> 59 }
>> 60 //////////////////////////////////////////////////////////////////////////////
>> 61 //
>> 62 void G4VAdjointReverseReaction::PreparePhysicsTable(const G4ParticleDefinition&)
>> 63 {;
47 } 64 }
48 <<
49 ////////////////////////////////////////////// 65 //////////////////////////////////////////////////////////////////////////////
>> 66 //
50 void G4VAdjointReverseReaction::BuildPhysicsTa 67 void G4VAdjointReverseReaction::BuildPhysicsTable(const G4ParticleDefinition&)
51 { 68 {
52 fCSManager->BuildCrossSectionMatrices(); // <<
53 fCSManager->BuildTotalSigmaTables(); <<
54 } <<
55 69
>> 70 theAdjointCSManager->BuildCrossSectionMatrices(); //do not worry it will be done just once
>> 71 theAdjointCSManager->BuildTotalSigmaTables();
>> 72
>> 73 }
56 ////////////////////////////////////////////// 74 //////////////////////////////////////////////////////////////////////////////
57 G4VParticleChange* G4VAdjointReverseReaction:: << 75 //
58 << 76 G4VParticleChange* G4VAdjointReverseReaction::PostStepDoIt(const G4Track& track, const G4Step& )
59 { << 77 {
>> 78
60 fParticleChange->Initialize(track); 79 fParticleChange->Initialize(track);
61 fAdjointModel->SampleSecondaries(track, fIsS << 80
62 << 81 /* if (IsFwdCSUsed && IsIntegralModeUsed){ //INtegral mode still unstable
>> 82 G4double Tkin = step.GetPostStepPoint()->GetKineticEnergy();
>> 83 G4double fwdCS = theAdjointCSManager->GetTotalForwardCS(track.GetDefinition(), Tkin, track.GetMaterialCutsCouple());
>> 84 //G4cout<<"lastCS "<<lastCS<<G4endl;
>> 85 if (fwdCS<lastCS*G4UniformRand()) { // the reaction does not take place, same integral method as the one used for forward ionisation in G4
>> 86 ClearNumberOfInteractionLengthLeft();
>> 87 return fParticleChange;
>> 88 }
>> 89
>> 90 }
>> 91 */
>> 92
>> 93 theAdjointEMModel->SampleSecondaries(track,
>> 94 IsScatProjToProjCase,
>> 95 fParticleChange);
>> 96
63 ClearNumberOfInteractionLengthLeft(); 97 ClearNumberOfInteractionLengthLeft();
64 return fParticleChange; 98 return fParticleChange;
>> 99
>> 100
>> 101
65 } 102 }
66 <<
67 ////////////////////////////////////////////// 103 //////////////////////////////////////////////////////////////////////////////
>> 104 //
68 G4double G4VAdjointReverseReaction::GetMeanFre 105 G4double G4VAdjointReverseReaction::GetMeanFreePath(const G4Track& track,
69 << 106 G4double ,
70 << 107 G4ForceCondition* condition)
71 { << 108 { *condition = NotForced;
72 *condition = NotForced; <<
73 G4double preStepKinEnergy = track.GetKinetic 109 G4double preStepKinEnergy = track.GetKineticEnergy();
74 110
75 if(track.GetTrackID() != fTrackId) << 111 if(track.GetTrackID() != trackid) {
76 { << 112 trackid = track.GetTrackID();
77 fTrackId = track.GetTrackID(); << 113 nstep = 0;
78 } 114 }
79 G4double sigma = fAdjointModel->AdjointCross << 115 ++nstep;
80 track.GetMaterialCutsCouple(), preStepKinE <<
81 116
82 G4double corr = fCSManager->GetCrossSectionC << 117
83 track.GetDefinition(), preStepKinEnergy, t << 118
84 fIsFwdCSUsed); << 119 /*G4double Sigma =
85 << 120 theAdjointEMModel->AdjointCrossSection(track.GetMaterialCutsCouple(),preStepKinEnergy,IsScatProjToProjCase);*/
86 if(std::fabs(corr) > 100.) << 121
87 { << 122 G4double Sigma =
88 sigma = 0.0; << 123 theAdjointEMModel->GetAdjointCrossSection(track.GetMaterialCutsCouple(),preStepKinEnergy,IsScatProjToProjCase);
89 } << 124
90 else << 125 //G4double sig = Sigma;
91 { << 126
92 sigma *= corr; << 127 G4double fwd_TotCS;
>> 128 G4double corr = theAdjointCSManager->GetCrossSectionCorrection(track.GetDefinition(),preStepKinEnergy,track.GetMaterialCutsCouple(),IsFwdCSUsed, fwd_TotCS);
>> 129
>> 130 if(std::fabs(corr) > 100.) { Sigma = 0.0; }
>> 131 else { Sigma *= corr; }
>> 132
>> 133 //G4cout<<fwd_TotCS<<G4endl;
>> 134 /*if (IsFwdCSUsed && IsIntegralModeUsed){ //take the maximum cross section only for charged particle
>> 135 G4double e_sigma_max, sigma_max;
>> 136 theAdjointCSManager->GetMaxFwdTotalCS(track.GetDefinition(),
>> 137 track.GetMaterialCutsCouple(), e_sigma_max, sigma_max);
>> 138 if (e_sigma_max > preStepKinEnergy){
>> 139 Sigma*=sigma_max/fwd_TotCS;
>> 140 }
93 } 141 }
>> 142 */
94 143
95 G4double mean_free_path = 1.e60; << 144 G4double mean_free_path = 1.e60 *mm;
96 if(sigma > 0.) << 145 if (Sigma>0) mean_free_path = 1./Sigma;
97 mean_free_path = 1. / sigma; << 146 lastCS=Sigma;
>> 147 /*
>> 148 if(nstep > 100) {
>> 149
>> 150 G4cout << "#* " << track.GetDefinition()->GetParticleName()
>> 151 << " " << GetProcessName()
>> 152 << " Nstep " << nstep
>> 153 << " E(MeV)= " << preStepKinEnergy << " Sig0= " << sig
>> 154 << " sig1= " << Sigma << " mfp= " << mean_free_path << G4endl;
>> 155
>> 156 }
>> 157 if (nstep > 20000) {
>> 158 exit(1);
>> 159 }
>> 160 */
>> 161 /*G4cout<<"Sigma "<<Sigma<<G4endl;
>> 162 G4cout<<"mean_free_path [mm] "<<mean_free_path/mm<<G4endl;
>> 163 */
>> 164
98 165
99 return mean_free_path; 166 return mean_free_path;
100 } << 167 }
101 168