| 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 // $Id: G4VAdjointReverseReaction.cc 87443 2014-12-04 12:26:31Z gunter $
>> 27 //
27 #include "G4VAdjointReverseReaction.hh" 28 #include "G4VAdjointReverseReaction.hh"
28 << 29 #include "G4SystemOfUnits.hh"
29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointCSManager.hh"
30 #include "G4ParticleChange.hh" << 31 #include "G4AdjointCSMatrix.hh"
>> 32 #include "G4AdjointInterpolator.hh"
>> 33 #include "G4AdjointCSMatrix.hh"
31 #include "G4VEmAdjointModel.hh" 34 #include "G4VEmAdjointModel.hh"
>> 35 #include "G4ElementTable.hh"
>> 36 #include "G4Element.hh"
>> 37 #include "G4Material.hh"
>> 38 #include "G4MaterialCutsCouple.hh"
>> 39 #include "G4AdjointCSManager.hh"
>> 40 #include "G4ParticleChange.hh"
>> 41 #include "G4AdjointElectron.hh"
32 42
33 G4VAdjointReverseReaction::G4VAdjointReverseRe <<
34 <<
35 : G4VDiscreteProcess(process_name) <<
36 { <<
37 fCSManager = G4AdjointCSManager::GetA <<
38 fIsScatProjToProj = whichScatCase; <<
39 fParticleChange = new G4ParticleChange(); <<
40 } <<
41 43
>> 44 G4VAdjointReverseReaction::
>> 45 G4VAdjointReverseReaction(G4String process_name, G4bool whichScatCase):
>> 46 G4VDiscreteProcess(process_name)
>> 47 {theAdjointCSManager = G4AdjointCSManager::GetAdjointCSManager();
>> 48 IsScatProjToProjCase=whichScatCase;
>> 49 fParticleChange=new G4ParticleChange();
>> 50 IsFwdCSUsed=false;
>> 51 IsIntegralModeUsed=false;
>> 52 lastCS=0.;
>> 53 trackid = nstep = 0;
>> 54 }
42 ////////////////////////////////////////////// 55 //////////////////////////////////////////////////////////////////////////////
43 G4VAdjointReverseReaction::~G4VAdjointReverseR << 56 //
44 { << 57 G4VAdjointReverseReaction::
45 if(fParticleChange) << 58 ~G4VAdjointReverseReaction()
46 delete fParticleChange; << 59 { if (fParticleChange) delete fParticleChange;
>> 60 }
>> 61 //////////////////////////////////////////////////////////////////////////////
>> 62 //
>> 63 void G4VAdjointReverseReaction::PreparePhysicsTable(const G4ParticleDefinition&)
>> 64 {;
47 } 65 }
48 <<
49 ////////////////////////////////////////////// 66 //////////////////////////////////////////////////////////////////////////////
>> 67 //
50 void G4VAdjointReverseReaction::BuildPhysicsTa 68 void G4VAdjointReverseReaction::BuildPhysicsTable(const G4ParticleDefinition&)
51 { 69 {
52 fCSManager->BuildCrossSectionMatrices(); // <<
53 fCSManager->BuildTotalSigmaTables(); <<
54 } <<
55 70
>> 71 theAdjointCSManager->BuildCrossSectionMatrices(); //do not worry it will be done just once
>> 72 theAdjointCSManager->BuildTotalSigmaTables();
>> 73
>> 74 }
56 ////////////////////////////////////////////// 75 //////////////////////////////////////////////////////////////////////////////
57 G4VParticleChange* G4VAdjointReverseReaction:: << 76 //
58 << 77 G4VParticleChange* G4VAdjointReverseReaction::PostStepDoIt(const G4Track& track, const G4Step& )
59 { << 78 {
>> 79
60 fParticleChange->Initialize(track); 80 fParticleChange->Initialize(track);
61 fAdjointModel->SampleSecondaries(track, fIsS << 81
62 << 82 /* if (IsFwdCSUsed && IsIntegralModeUsed){ //INtegral mode still unstable
>> 83 G4double Tkin = step.GetPostStepPoint()->GetKineticEnergy();
>> 84 G4double fwdCS = theAdjointCSManager->GetTotalForwardCS(track.GetDefinition(), Tkin, track.GetMaterialCutsCouple());
>> 85 //G4cout<<"lastCS "<<lastCS<<G4endl;
>> 86 if (fwdCS<lastCS*G4UniformRand()) { // the reaction does not take place, same integral method as the one used for forward ionisation in G4
>> 87 ClearNumberOfInteractionLengthLeft();
>> 88 return fParticleChange;
>> 89 }
>> 90
>> 91 }
>> 92 */
>> 93
>> 94 theAdjointEMModel->SampleSecondaries(track,
>> 95 IsScatProjToProjCase,
>> 96 fParticleChange);
>> 97
63 ClearNumberOfInteractionLengthLeft(); 98 ClearNumberOfInteractionLengthLeft();
64 return fParticleChange; 99 return fParticleChange;
>> 100
>> 101
>> 102
65 } 103 }
66 <<
67 ////////////////////////////////////////////// 104 //////////////////////////////////////////////////////////////////////////////
>> 105 //
68 G4double G4VAdjointReverseReaction::GetMeanFre 106 G4double G4VAdjointReverseReaction::GetMeanFreePath(const G4Track& track,
69 << 107 G4double ,
70 << 108 G4ForceCondition* condition)
71 { << 109 { *condition = NotForced;
72 *condition = NotForced; <<
73 G4double preStepKinEnergy = track.GetKinetic 110 G4double preStepKinEnergy = track.GetKineticEnergy();
74 111
75 if(track.GetTrackID() != fTrackId) << 112 if(track.GetTrackID() != trackid) {
76 { << 113 trackid = track.GetTrackID();
77 fTrackId = track.GetTrackID(); << 114 nstep = 0;
78 } 115 }
79 G4double sigma = fAdjointModel->AdjointCross << 116 ++nstep;
80 track.GetMaterialCutsCouple(), preStepKinE <<
81 117
82 G4double corr = fCSManager->GetCrossSectionC << 118
83 track.GetDefinition(), preStepKinEnergy, t << 119
84 fIsFwdCSUsed); << 120 /*G4double Sigma =
85 << 121 theAdjointEMModel->AdjointCrossSection(track.GetMaterialCutsCouple(),preStepKinEnergy,IsScatProjToProjCase);*/
86 if(std::fabs(corr) > 100.) << 122
87 { << 123 G4double Sigma =
88 sigma = 0.0; << 124 theAdjointEMModel->GetAdjointCrossSection(track.GetMaterialCutsCouple(),preStepKinEnergy,IsScatProjToProjCase);
89 } << 125
90 else << 126 //G4double sig = Sigma;
91 { << 127
92 sigma *= corr; << 128 G4double fwd_TotCS;
>> 129 G4double corr = theAdjointCSManager->GetCrossSectionCorrection(track.GetDefinition(),preStepKinEnergy,track.GetMaterialCutsCouple(),IsFwdCSUsed, fwd_TotCS);
>> 130
>> 131 if(std::fabs(corr) > 100.) { Sigma = 0.0; }
>> 132 else { Sigma *= corr; }
>> 133
>> 134 //G4cout<<fwd_TotCS<<G4endl;
>> 135 /*if (IsFwdCSUsed && IsIntegralModeUsed){ //take the maximum cross section only for charged particle
>> 136 G4double e_sigma_max, sigma_max;
>> 137 theAdjointCSManager->GetMaxFwdTotalCS(track.GetDefinition(),
>> 138 track.GetMaterialCutsCouple(), e_sigma_max, sigma_max);
>> 139 if (e_sigma_max > preStepKinEnergy){
>> 140 Sigma*=sigma_max/fwd_TotCS;
>> 141 }
93 } 142 }
>> 143 */
94 144
95 G4double mean_free_path = 1.e60; << 145 G4double mean_free_path = 1.e60 *mm;
96 if(sigma > 0.) << 146 if (Sigma>0) mean_free_path = 1./Sigma;
97 mean_free_path = 1. / sigma; << 147 lastCS=Sigma;
>> 148 /*
>> 149 if(nstep > 100) {
>> 150
>> 151 G4cout << "#* " << track.GetDefinition()->GetParticleName()
>> 152 << " " << GetProcessName()
>> 153 << " Nstep " << nstep
>> 154 << " E(MeV)= " << preStepKinEnergy << " Sig0= " << sig
>> 155 << " sig1= " << Sigma << " mfp= " << mean_free_path << G4endl;
>> 156
>> 157 }
>> 158 if (nstep > 20000) {
>> 159 exit(1);
>> 160 }
>> 161 */
>> 162 /*G4cout<<"Sigma "<<Sigma<<G4endl;
>> 163 G4cout<<"mean_free_path [mm] "<<mean_free_path/mm<<G4endl;
>> 164 */
>> 165
98 166
99 return mean_free_path; 167 return mean_free_path;
100 } << 168 }
101 169