Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/ParticleFluence/ConcentricSpheres/src/SteppingAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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::fkNumberScoringShells>
 52   SteppingAction::fkArrayScoringShellNames = {"tracker", "emCalo", "hadCalo"};
 53 
 54 const std::array<G4String, SteppingAction::fkNumberKinematicRegions>
 55   SteppingAction::fkArrayKinematicRegionNames = {"", "below 20 MeV", "above 20 MeV"};
 56 
 57 const std::array<G4String, SteppingAction::fkNumberScoringPositions>
 58   SteppingAction::fkArrayScoringPositionNames = {"forward", "backward"};
 59 
 60 const std::array<G4String, SteppingAction::fkNumberParticleTypes>
 61   SteppingAction::fkArrayParticleTypeNames = {"all",      "electron",   "gamma",      "muon",
 62                                               "neutrino", "pion",       "neutron",    "proton",
 63                                               "ion",      "otherMeson", "otherBaryon"};
 64 
 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 66 
 67 G4int SteppingAction::GetIndex(const G4int iScoringShell, const G4int iKinematicRegion,
 68                                const G4int iScoringPosition, const G4int iParticleType)
 69 {
 70   G4int index = -1;
 71   if (iScoringShell >= 0 && iScoringShell < fkNumberScoringShells && iKinematicRegion >= 0
 72       && iKinematicRegion < fkNumberKinematicRegions && iScoringPosition >= 0
 73       && iScoringPosition < fkNumberScoringPositions && iParticleType >= 0
 74       && iParticleType < fkNumberParticleTypes)
 75   {
 76     index =
 77       iScoringShell * fkNumberKinematicRegions * fkNumberScoringPositions * fkNumberParticleTypes
 78       + iKinematicRegion * fkNumberScoringPositions * fkNumberParticleTypes
 79       + iScoringPosition * fkNumberParticleTypes + iParticleType;
 80   }
 81   if (index < 0 || index >= fkNumberCombinations) {
 82     G4cerr << "SteppingAction::GetIndex : WRONG index=" << index << "  set it to 0 !" << G4endl;
 83     index = 0;
 84   }
 85   return index;
 86 }
 87 
 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 89 
 90 SteppingAction::SteppingAction() : G4UserSteppingAction()
 91 {
 92   Initialize();
 93 }
 94 
 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 96 
 97 void SteppingAction::Initialize()
 98 {
 99   // Initialization needed at the beginning of each Run
100   fPrimaryParticleId = 0;
101   fPrimaryParticleEnergy = 0.0;
102   fPrimaryParticleDirection = G4ThreeVector(0.0, 0.0, 1.0);
103   fTrackerMaterialName = fEmCaloMaterialName = fHadCaloMaterialName = "";
104   fIsFirstStepOfTheEvent = true;
105   fIsFirstStepInTracker = fIsFirstStepInEmCalo = fIsFirstStepInHadCalo = true;
106   fIsFirstStepInScoringTrackerShell = fIsFirstStepInScoringEmCaloShell =
107     fIsFirstStepInScoringHadCaloShell = true;
108   fCubicVolumeScoringTrackerShell = fCubicVolumeScoringEmCaloShell =
109     fCubicVolumeScoringHadCaloShell = 1.0;
110   for (G4int i = 0; i < fkNumberCombinations; ++i) {
111     fArraySumStepLengths[i] = 0.0;
112   }
113   /*
114   for ( G4int i = 0; i < fkNumberCombinations; ++i ) fArraySumStepLengths[i] = 999.9;
115   G4cout << " fkNumberCombinations=" << fkNumberCombinations << G4endl;
116   for ( G4int i = 0; i < fkNumberScoringShells; ++i ) {
117     for ( G4int j = 0; j < fkNumberKinematicRegions; ++j ) {
118       for ( G4int k = 0; k < fkNumberScoringPositions; ++k ) {
119         for ( G4int ll = 0; ll < fkNumberParticleTypes; ++ll ) {
120           G4int index = GetIndex( i, j, k, ll );
121           G4cout << "(i, j, k, ll)=(" << i << ", " << j << ", " << k << ", "
122                  << ll << ")  ->" << index;
123           if ( fArraySumStepLengths[ index ] < 1.0 ) G4cout << " <=== REPEATED!";
124           else                                       fArraySumStepLengths[ index ] = 0.0;
125           G4cout << G4endl;
126         }
127       }
128     }
129   }
130   for ( G4int i = 0; i < fkNumberCombinations; ++i ) {
131     if ( fArraySumStepLengths[i] > 999.0 ) G4cout << " i=" << i << " NOT COVERED !" << G4endl;
132   }
133   */
134 }
135 
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 
138 void SteppingAction::UserSteppingAction(const G4Step* theStep)
139 {
140   // Get information on the primary particle
141   if (fIsFirstStepOfTheEvent) {
142     if (theStep->GetTrack()->GetParentID() == 0) {
143       fPrimaryParticleId = theStep->GetTrack()->GetDefinition()->GetPDGEncoding();
144       fPrimaryParticleEnergy = theStep->GetPreStepPoint()->GetKineticEnergy();
145       fPrimaryParticleDirection = theStep->GetPreStepPoint()->GetMomentumDirection();
146       if (fRunPtr) {
147         fRunPtr->SetPrimaryParticleId(fPrimaryParticleId);
148         fRunPtr->SetPrimaryParticleEnergy(fPrimaryParticleEnergy);
149         fRunPtr->SetPrimaryParticleDirection(fPrimaryParticleDirection);
150       }
151       fIsFirstStepOfTheEvent = false;
152     }
153   }
154   // Get information on the materials
155   if (fIsFirstStepInTracker
156       && theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiTrackerShell")
157   {
158     fTrackerMaterialName = theStep->GetPreStepPoint()->GetMaterial()->GetName();
159     if (fRunPtr) fRunPtr->SetTrackerMaterialName(fTrackerMaterialName);
160     fIsFirstStepInTracker = false;
161   }
162   if (fIsFirstStepInEmCalo
163       && theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiEmCaloShell")
164   {
165     fEmCaloMaterialName = theStep->GetPreStepPoint()->GetMaterial()->GetName();
166     if (fRunPtr) fRunPtr->SetEmCaloMaterialName(fEmCaloMaterialName);
167     fIsFirstStepInEmCalo = false;
168   }
169   if (fIsFirstStepInHadCalo
170       && theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiHadCaloShell")
171   {
172     fHadCaloMaterialName = theStep->GetPreStepPoint()->GetMaterial()->GetName();
173     if (fRunPtr) fRunPtr->SetHadCaloMaterialName(fHadCaloMaterialName);
174     fIsFirstStepInHadCalo = false;
175   }
176   // Get information on step lengths in the scoring shells
177   G4int iScoringShell = -1;
178   if (theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiScoringTrackerShell") {
179     iScoringShell = 0;
180     if (fIsFirstStepInScoringTrackerShell) {
181       fCubicVolumeScoringTrackerShell =
182         theStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetSolid()->GetCubicVolume();
183       if (fRunPtr) fRunPtr->SetCubicVolumeScoringTrackerShell(fCubicVolumeScoringTrackerShell);
184       fIsFirstStepInScoringTrackerShell = false;
185     }
186   }
187   else if (theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiScoringEmCaloShell")
188   {
189     iScoringShell = 1;
190     if (fIsFirstStepInScoringEmCaloShell) {
191       fCubicVolumeScoringEmCaloShell =
192         theStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetSolid()->GetCubicVolume();
193       if (fRunPtr) fRunPtr->SetCubicVolumeScoringEmCaloShell(fCubicVolumeScoringEmCaloShell);
194       fIsFirstStepInScoringEmCaloShell = false;
195     }
196   }
197   else if (theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "physiScoringHadCaloShell")
198   {
199     iScoringShell = 2;
200     if (fIsFirstStepInScoringHadCaloShell) {
201       fCubicVolumeScoringHadCaloShell =
202         theStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetSolid()->GetCubicVolume();
203       if (fRunPtr) fRunPtr->SetCubicVolumeScoringHadCaloShell(fCubicVolumeScoringHadCaloShell);
204       fIsFirstStepInScoringHadCaloShell = false;
205     }
206   }
207   if (iScoringShell >= 0) {
208     G4double stepLength = theStep->GetTrack()->GetStepLength() * theStep->GetTrack()->GetWeight();
209     G4int absPdg = theStep->GetTrack()->GetDefinition() == nullptr
210                      ? 0
211                      : std::abs(theStep->GetTrack()->GetDefinition()->GetPDGEncoding());
212     /*
213     G4cout << theStep->GetTrack()->GetDefinition()->GetParticleName() << "  absPdg=" << absPdg
214            << "  Ekin[MeV]=" << theStep->GetPreStepPoint()->GetKineticEnergy()
215            << "  r[mm]=" << theStep->GetTrack()->GetPosition().mag()
216            << "  z[mm]=" << theStep->GetTrack()->GetPosition().z()
217            << "  " << theStep->GetTrack()->GetVolume()->GetName()
218            << "  " << theStep->GetTrack()->GetMaterial()->GetName()
219            << "  L[mm]=" << stepLength << "  "
220            << ( fPrimaryParticleDirection.dot( theStep->GetTrack()->GetPosition().unit() ) > 0.0
221                 ? "forward" : "backward" ) << G4endl;
222     */
223     // Three kinematical regions:  [0] : any value ;  [1] : below 20 MeV ;  [2] : above 20 MeV
224     G4int iKinematicRegion = theStep->GetPreStepPoint()->GetKineticEnergy() < 20.0 ? 1 : 2;
225     // Two scoring positions:  [0] : forward hemisphere ;  [1] : backward hemisphere
226     // (with respect to the primary particle initial direction)
227     G4int iScoringPosition =
228       fPrimaryParticleDirection.dot(theStep->GetTrack()->GetPosition().unit()) > 0.0 ? 0 : 1;
229     G4int iParticleType = -1;
230     if (absPdg == 11)
231       iParticleType = 1;  // electron (and positron)
232     else if (absPdg == 22)
233       iParticleType = 2;  // gamma
234     else if (absPdg == 13)
235       iParticleType = 3;  // muons (mu- and mu+)
236     else if (absPdg == 12 || absPdg == 14 || absPdg == 16)
237       iParticleType = 4;  // neutrinos
238     // (and anti-neutrinos), all flavors
239     else if (absPdg == 111 || absPdg == 211)
240       iParticleType = 5;  // (charged) pions
241     else if (absPdg == 2112)
242       iParticleType = 6;  // neutron (and anti-neutron)
243     else if (absPdg == 2212)
244       iParticleType = 7;  // proton  (and anti-proton)
245     else if (G4IonTable::IsIon(theStep->GetTrack()->GetDefinition()) ||  // ions (and anti-ions)
246              G4IonTable::IsAntiIon(theStep->GetTrack()->GetDefinition()))
247       iParticleType = 8;
248     else if (absPdg < 1000)
249       iParticleType = 9;  // other mesons (e.g. kaons) (Note: this works
250                           // in most cases, but not always!)
251     else if (absPdg > 1000)
252       iParticleType = 10;  // other baryons (e.g. hyperons, anti-hyperons,
253                            // etc.)
254     // Consider the specific case : scoring shell, kinematic region, scoring position, and
255     // particle type
256     G4int index = GetIndex(iScoringShell, iKinematicRegion, iScoringPosition, iParticleType);
257     fArraySumStepLengths[index] += stepLength;
258     // Consider the "all" particle case, with the same scoring shell, kinematic region and
259     // scoring position
260     index = GetIndex(iScoringShell, iKinematicRegion, iScoringPosition, 0);
261     fArraySumStepLengths[index] += stepLength;
262     // Consider the "any" kinematic region case, with the same scoring shell, scoring position
263     // and particle type
264     index = GetIndex(iScoringShell, 0, iScoringPosition, iParticleType);
265     fArraySumStepLengths[index] += stepLength;
266     // Consider the "any" kinematic region and "all" particle, with the same scoring shell and
267     // scoring position
268     index = GetIndex(iScoringShell, 0, iScoringPosition, 0);
269     fArraySumStepLengths[index] += stepLength;
270     if (fRunPtr) fRunPtr->SetSteppingArray(fArraySumStepLengths);
271   }
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275