Geant4 Cross Reference |
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 /// \file SteppingAction.cc 27 /// \brief Implementation of the SteppingAction class 28 // 29 // 30 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 33 34 #include "SteppingAction.hh" 35 36 #include "Run.hh" 37 38 #include "G4IonTable.hh" 39 #include "G4LossTableManager.hh" 40 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleTypes.hh" 42 #include "G4Step.hh" 43 #include "G4StepPoint.hh" 44 #include "G4SystemOfUnits.hh" 45 #include "G4TouchableHistory.hh" 46 #include "G4Track.hh" 47 #include "G4VPhysicalVolume.hh" 48 #include "G4VSolid.hh" 49 #include "G4VTouchable.hh" 50 51 const std::array<G4String, SteppingAction::fkNumberScoringVolumes> 52 SteppingAction::fkArrayScoringVolumeNames = {"downstream", "side", "upstream"}; 53 54 const std::array<G4String, SteppingAction::fkNumberKinematicRegions> 55 SteppingAction::fkArrayKinematicRegionNames = {"", "below 20 MeV", "above 20 MeV"}; 56 57 const std::array<G4String, SteppingAction::fkNumberParticleTypes> 58 SteppingAction::fkArrayParticleTypeNames = {"all", "electron", "gamma", "muon", 59 "neutrino", "pion", "neutron", "proton", 60 "ion", "otherMeson", "otherBaryon"}; 61 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 63 64 G4int SteppingAction::GetIndex(const G4int iScoringVolume, const G4int iKinematicRegion, 65 const G4int iParticleType) 66 { 67 G4int index = -1; 68 if (iScoringVolume >= 0 && iScoringVolume < fkNumberScoringVolumes && iKinematicRegion >= 0 69 && iKinematicRegion < fkNumberKinematicRegions && iParticleType >= 0 70 && iParticleType < fkNumberParticleTypes) 71 { 72 index = iScoringVolume * fkNumberKinematicRegions * fkNumberParticleTypes 73 + iKinematicRegion * fkNumberParticleTypes + iParticleType; 74 } 75 if (index < 0 || index >= fkNumberCombinations) { 76 G4cerr << "SteppingAction::GetIndex : WRONG index=" << index << " set it to 0 !" << G4endl; 77 index = 0; 78 } 79 return index; 80 } 81 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 83 84 SteppingAction::SteppingAction() : G4UserSteppingAction() 85 { 86 Initialize(); 87 } 88 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 90 91 void SteppingAction::Initialize() 92 { 93 // Initialization needed at the beginning of each Run 94 fPrimaryParticleId = 0; 95 fPrimaryParticleEnergy = 0.0; 96 fPrimaryParticleDirection = G4ThreeVector(0.0, 0.0, 1.0); 97 fTargetMaterialName = ""; 98 fIsFirstStepOfTheEvent = true; 99 fIsFirstStepInTarget = true; 100 fIsFirstStepInScoringUpDown = true; 101 fIsFirstStepInScoringSide = true; 102 fCubicVolumeScoringUpDown = 1.0; 103 fCubicVolumeScoringSide = 1.0; 104 for (G4int i = 0; i < fkNumberCombinations; ++i) { 105 fArraySumStepLengths[i] = 0.0; 106 } 107 /* 108 for ( G4int i = 0; i < fkNumberCombinations; ++i ) fArraySumStepLengths[i] = 999.9; 109 G4cout << " fkNumberCombinations=" << fkNumberCombinations << G4endl; 110 for ( G4int i = 0; i < fkNumberScoringVolumes; ++i ) { 111 for ( G4int j = 0; j < fkNumberKinematicRegions; ++j ) { 112 for ( G4int k = 0; k < fkNumberParticleTypes; ++k ) { 113 G4int index = GetIndex( i, j, k ); 114 G4cout << "(i, j, k)=(" << i << ", " << j << ", " << k << ") ->" << index; 115 if ( fArraySumStepLengths[ index ] < 1.0 ) G4cout << " <=== REPEATED!"; 116 else fArraySumStepLengths[ index ] = 0.0; 117 G4cout << G4endl; 118 } 119 } 120 } 121 for ( G4int i = 0; i < fkNumberCombinations; ++i ) { 122 if ( fArraySumStepLengths[i] > 999.0 ) G4cout << " i=" << i << " NOT COVERED !" << G4endl; 123 } 124 */ 125 } 126 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 128 129 void SteppingAction::UserSteppingAction(const G4Step* theStep) 130 { 131 // Get information on the primary particle 132 if (fIsFirstStepOfTheEvent) { 133 if (theStep->GetTrack()->GetParentID() == 0) { 134 fPrimaryParticleId = theStep->GetTrack()->GetDefinition()->GetPDGEncoding(); 135 fPrimaryParticleEnergy = theStep->GetPreStepPoint()->GetKineticEnergy(); 136 fPrimaryParticleDirection = theStep->GetPreStepPoint()->GetMomentumDirection(); 137 if (fRunPtr) { 138 fRunPtr->SetPrimaryParticleId(fPrimaryParticleId); 139 fRunPtr->SetPrimaryParticleEnergy(fPrimaryParticleEnergy); 140 fRunPtr->SetPrimaryParticleDirection(fPrimaryParticleDirection); 141 } 142 fIsFirstStepOfTheEvent = false; 143 } 144 } 145 // Get information on the target material 146 if (fIsFirstStepInTarget 147 && theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiLayer") 148 { 149 fTargetMaterialName = theStep->GetPreStepPoint()->GetMaterial()->GetName(); 150 if (fRunPtr) fRunPtr->SetTargetMaterialName(fTargetMaterialName); 151 fIsFirstStepInTarget = false; 152 } 153 // Get information on step lengths in the scoring volumes 154 G4int iScoringVolume = -1; 155 if (theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiScoringDownstream") { 156 iScoringVolume = 0; 157 if (fIsFirstStepInScoringUpDown) { 158 fCubicVolumeScoringUpDown = 159 theStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetSolid()->GetCubicVolume(); 160 if (fRunPtr) fRunPtr->SetCubicVolumeScoringUpDown(fCubicVolumeScoringUpDown); 161 fIsFirstStepInScoringUpDown = false; 162 } 163 } 164 else if (theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiScoringSide") { 165 iScoringVolume = 1; 166 if (fIsFirstStepInScoringSide) { 167 fCubicVolumeScoringSide = 168 theStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetSolid()->GetCubicVolume(); 169 if (fRunPtr) fRunPtr->SetCubicVolumeScoringSide(fCubicVolumeScoringSide); 170 fIsFirstStepInScoringSide = false; 171 } 172 } 173 else if (theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiScoringUpstream") { 174 iScoringVolume = 2; 175 if (fIsFirstStepInScoringUpDown) { 176 fCubicVolumeScoringUpDown = 177 theStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetSolid()->GetCubicVolume(); 178 if (fRunPtr) fRunPtr->SetCubicVolumeScoringUpDown(fCubicVolumeScoringUpDown); 179 fIsFirstStepInScoringUpDown = false; 180 } 181 } 182 if (iScoringVolume >= 0) { 183 // In the case of the upstream scoring volume, consider only particles whose direction 184 // is opposite with respect to the primary particle (this is needed, in particular, 185 // for avoiding to account the incoming, primary beam particle in the "upstream" fluence). 186 if (iScoringVolume == 2 187 && fPrimaryParticleDirection.dot(theStep->GetPreStepPoint()->GetMomentumDirection()) > 0.0) 188 { 189 return; 190 } 191 G4double stepLength = theStep->GetTrack()->GetStepLength() * theStep->GetTrack()->GetWeight(); 192 G4int absPdg = theStep->GetTrack()->GetDefinition() == nullptr 193 ? 0 194 : std::abs(theStep->GetTrack()->GetDefinition()->GetPDGEncoding()); 195 /* 196 G4cout << std::setprecision(6) 197 << theStep->GetTrack()->GetDefinition()->GetParticleName() << " absPdg=" << absPdg 198 << " Ekin[MeV]=" << theStep->GetPreStepPoint()->GetKineticEnergy() 199 << " (rho,z)[mm]=(" << theStep->GetTrack()->GetPosition().perp() 200 << "," << theStep->GetTrack()->GetPosition().z() << ")" 201 << " " << theStep->GetTrack()->GetVolume()->GetName() 202 << " " << theStep->GetTrack()->GetMaterial()->GetName() 203 << " L[mm]=" << stepLength << " " 204 << ( fPrimaryParticleDirection.dot( 205 theStep->GetPreStepPoint()->GetMomentumDirection() ) > 0.0 206 ? "forward" : "backward" ) 207 << G4endl; 208 */ 209 // Three kinematical regions: [0] : any value ; [1] : below 20 MeV ; [2] : above 20 MeV 210 G4int iKinematicRegion = theStep->GetPreStepPoint()->GetKineticEnergy() < 20.0 ? 1 : 2; 211 G4int iParticleType = -1; 212 if (absPdg == 11) 213 iParticleType = 1; // electron (and positron) 214 else if (absPdg == 22) 215 iParticleType = 2; // gamma 216 else if (absPdg == 13) 217 iParticleType = 3; // muons (mu- and mu+) 218 else if (absPdg == 12 || absPdg == 14 || absPdg == 16) 219 iParticleType = 4; // neutrinos 220 //(and anti-neutrinos), all flavors 221 else if (absPdg == 111 || absPdg == 211) 222 iParticleType = 5; // (charged) pions 223 else if (absPdg == 2112) 224 iParticleType = 6; // neutron (and anti-neutron) 225 else if (absPdg == 2212) 226 iParticleType = 7; // proton (and anti-proton) 227 else if (G4IonTable::IsIon(theStep->GetTrack()->GetDefinition()) || // ions (and anti-ions) 228 G4IonTable::IsAntiIon(theStep->GetTrack()->GetDefinition())) 229 iParticleType = 8; 230 else if (absPdg < 1000) 231 iParticleType = 9; // other mesons (e.g. kaons) (Note: this works 232 // in most cases, but not always!) 233 else if (absPdg > 1000) 234 iParticleType = 10; // other baryons (e.g. hyperons, anti-hyperons, 235 // etc.) 236 // Consider the specific case : scoring volume, kinematic region and particle type 237 G4int index = GetIndex(iScoringVolume, iKinematicRegion, iParticleType); 238 fArraySumStepLengths[index] += stepLength; 239 // Consider the "all" particle case, with the same scoring volume and kinematic region 240 index = GetIndex(iScoringVolume, iKinematicRegion, 0); 241 fArraySumStepLengths[index] += stepLength; 242 // Consider the "any" kinematic region case, with the same scoring volume and particle type 243 index = GetIndex(iScoringVolume, 0, iParticleType); 244 fArraySumStepLengths[index] += stepLength; 245 // Consider the "any" kinematic region and "all" particle, with the same scoring volume 246 index = GetIndex(iScoringVolume, 0, 0); 247 fArraySumStepLengths[index] += stepLength; 248 if (fRunPtr) fRunPtr->SetSteppingArray(fArraySumStepLengths); 249 } 250 } 251 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 253