Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 /// \file SteppingAction.cc 27 /// \brief Implementation of the SteppingActio 28 // 29 // 30 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oo 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::fkN 52 SteppingAction::fkArrayScoringVolumeNames = 53 54 const std::array<G4String, SteppingAction::fkN 55 SteppingAction::fkArrayKinematicRegionNames 56 57 const std::array<G4String, SteppingAction::fkN 58 SteppingAction::fkArrayParticleTypeNames = { 59 60 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 64 G4int SteppingAction::GetIndex(const G4int iSc 65 const G4int iPa 66 { 67 G4int index = -1; 68 if (iScoringVolume >= 0 && iScoringVolume < 69 && iKinematicRegion < fkNumberKinematicR 70 && iParticleType < fkNumberParticleTypes 71 { 72 index = iScoringVolume * fkNumberKinematic 73 + iKinematicRegion * fkNumberParti 74 } 75 if (index < 0 || index >= fkNumberCombinatio 76 G4cerr << "SteppingAction::GetIndex : WRON 77 index = 0; 78 } 79 return index; 80 } 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 83 84 SteppingAction::SteppingAction() : G4UserStepp 85 { 86 Initialize(); 87 } 88 89 //....oooOO0OOooo........oooOO0OOooo........oo 90 91 void SteppingAction::Initialize() 92 { 93 // Initialization needed at the beginning of 94 fPrimaryParticleId = 0; 95 fPrimaryParticleEnergy = 0.0; 96 fPrimaryParticleDirection = G4ThreeVector(0. 97 fAbsorberMaterialName = ""; 98 fActiveMaterialName = ""; 99 fIsFirstStepOfTheEvent = true; 100 fIsFirstStepInAbsorberLayer = true; 101 fIsFirstStepInActiveLayer = true; 102 fIsFirstStepInScoringUpDown = true; 103 fIsFirstStepInScoringSide = true; 104 fCubicVolumeScoringUpDown = 1.0; 105 fCubicVolumeScoringSide = 1.0; 106 for (G4int i = 0; i < fkNumberCombinations; 107 fArraySumStepLengths[i] = 0.0; 108 } 109 /* 110 for ( G4int i = 0; i < fkNumberCombinations; 111 G4cout << " fkNumberCombinations=" << fkNumb 112 for ( G4int i = 0; i < fkNumberScoringVolume 113 for ( G4int j = 0; j < fkNumberKinematicRe 114 for ( G4int k = 0; k < fkNumberParticleT 115 G4int index = GetIndex( i, j, k ); 116 G4cout << "(i, j, k)=(" << i << ", " < 117 if ( fArraySumStepLengths[ index ] < 1 118 else 119 G4cout << G4endl; 120 } 121 } 122 } 123 for ( G4int i = 0; i < fkNumberCombinations; 124 if ( fArraySumStepLengths[i] > 999.0 ) G4c 125 } 126 */ 127 } 128 129 //....oooOO0OOooo........oooOO0OOooo........oo 130 131 void SteppingAction::UserSteppingAction(const 132 { 133 // Get information on the primary particle 134 if (fIsFirstStepOfTheEvent) { 135 if (theStep->GetTrack()->GetParentID() == 136 fPrimaryParticleId = theStep->GetTrack() 137 fPrimaryParticleEnergy = theStep->GetPre 138 fPrimaryParticleDirection = theStep->Get 139 if (fRunPtr) { 140 fRunPtr->SetPrimaryParticleId(fPrimary 141 fRunPtr->SetPrimaryParticleEnergy(fPri 142 fRunPtr->SetPrimaryParticleDirection(f 143 } 144 fIsFirstStepOfTheEvent = false; 145 } 146 } 147 // Get information on the materials of the c 148 if (fIsFirstStepInAbsorberLayer 149 && theStep->GetPreStepPoint()->GetPhysic 150 { 151 fAbsorberMaterialName = theStep->GetPreSte 152 if (fRunPtr) fRunPtr->SetAbsorberMaterialN 153 fIsFirstStepInAbsorberLayer = false; 154 } 155 if (fIsFirstStepInActiveLayer 156 && theStep->GetPreStepPoint()->GetPhysic 157 { 158 fActiveMaterialName = theStep->GetPreStepP 159 if (fRunPtr) fRunPtr->SetActiveMaterialNam 160 fIsFirstStepInActiveLayer = false; 161 } 162 // Get information on step lengths in the sc 163 G4int iScoringVolume = -1; 164 if (theStep->GetPreStepPoint()->GetPhysicalV 165 iScoringVolume = 0; 166 if (fIsFirstStepInScoringUpDown) { 167 fCubicVolumeScoringUpDown = 168 theStep->GetTrack()->GetVolume()->GetL 169 if (fRunPtr) fRunPtr->SetCubicVolumeScor 170 fIsFirstStepInScoringUpDown = false; 171 } 172 } 173 else if (theStep->GetPreStepPoint()->GetPhys 174 iScoringVolume = 1; 175 if (fIsFirstStepInScoringSide) { 176 fCubicVolumeScoringSide = 177 theStep->GetTrack()->GetVolume()->GetL 178 if (fRunPtr) fRunPtr->SetCubicVolumeScor 179 fIsFirstStepInScoringSide = false; 180 } 181 } 182 else if (theStep->GetPreStepPoint()->GetPhys 183 iScoringVolume = 2; 184 if (fIsFirstStepInScoringUpDown) { 185 fCubicVolumeScoringUpDown = 186 theStep->GetTrack()->GetVolume()->GetL 187 if (fRunPtr) fRunPtr->SetCubicVolumeScor 188 fIsFirstStepInScoringUpDown = false; 189 } 190 } 191 if (iScoringVolume >= 0) { 192 // In the case of the upstream scoring vol 193 // is opposite with respect to the primary 194 // for avoiding to account the incoming, p 195 if (iScoringVolume == 2 196 && fPrimaryParticleDirection.dot(theSt 197 return; 198 G4double stepLength = theStep->GetTrack()- 199 G4int absPdg = theStep->GetTrack()->GetDef 200 ? 0 201 : std::abs(theStep->GetTr 202 /* 203 G4cout << std::setprecision(6) 204 << theStep->GetTrack()->GetDefiniti 205 << " Ekin[MeV]=" << theStep->GetPr 206 << " (rho,z)[mm]=(" << theStep->Ge 207 << "," << theStep->GetTrack()->GetP 208 << " " << theStep->GetTrack()->Get 209 << " " << theStep->GetTrack()->Get 210 << " L[mm]=" << stepLength << " " 211 << ( fPrimaryParticleDirection.dot( 212 theStep->GetPreStepPoint()-> 213 ? "forward" : "backward" ) 214 << G4endl; 215 */ 216 // Three kinematical regions: [0] : any v 217 G4int iKinematicRegion = theStep->GetPreSt 218 G4int iParticleType = -1; 219 if (absPdg == 11) 220 iParticleType = 1; // electron (and pos 221 else if (absPdg == 22) 222 iParticleType = 2; // gamma 223 else if (absPdg == 13) 224 iParticleType = 3; // muons (mu- and mu 225 else if (absPdg == 12 || absPdg == 14 || a 226 iParticleType = 4; // neutrinos 227 // (and anti-neutrinos), all flavors 228 else if (absPdg == 111 || absPdg == 211) 229 iParticleType = 5; // (charged) pions 230 else if (absPdg == 2112) 231 iParticleType = 6; // neutron (and anti 232 else if (absPdg == 2212) 233 iParticleType = 7; // proton (and anti 234 else if (G4IonTable::IsIon(theStep->GetTra 235 G4IonTable::IsAntiIon(theStep->Ge 236 iParticleType = 8; 237 else if (absPdg < 1000) 238 iParticleType = 9; // other mesons (e.g 239 // in most cases, bu 240 else if (absPdg > 1000) 241 iParticleType = 10; // other baryons (e 242 // etc.) 243 // Consider the specific case : scoring vo 244 G4int index = GetIndex(iScoringVolume, iKi 245 fArraySumStepLengths[index] += stepLength; 246 // Consider the "all" particle case, with 247 index = GetIndex(iScoringVolume, iKinemati 248 fArraySumStepLengths[index] += stepLength; 249 // Consider the "any" kinematic region cas 250 index = GetIndex(iScoringVolume, 0, iParti 251 fArraySumStepLengths[index] += stepLength; 252 // Consider the "any" kinematic region and 253 index = GetIndex(iScoringVolume, 0, 0); 254 fArraySumStepLengths[index] += stepLength; 255 if (fRunPtr) fRunPtr->SetSteppingArray(fAr 256 } 257 } 258 259 //....oooOO0OOooo........oooOO0OOooo........oo 260