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