ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SteppingAction.cc
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 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "SteppingAction.hh"
34 #include "Run.hh"
35 #include "HistoManager.hh"
36 
37 #include "G4ParticleTypes.hh"
38 #include "G4RunManager.hh"
39 #include "G4HadronicProcess.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 
45 { }
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 
50 { }
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 
55 {
56  Run* run
58 
59  // count processes
60  //
61  const G4StepPoint* endPoint = aStep->GetPostStepPoint();
62  G4VProcess* process =
63  const_cast<G4VProcess*>(endPoint->GetProcessDefinedStep());
64  run->CountProcesses(process);
65 
66  // check that an real interaction occured (eg. not a transportation)
67  G4StepStatus stepStatus = endPoint->GetStepStatus();
68  G4bool transmit = (stepStatus==fGeomBoundary || stepStatus==fWorldBoundary);
69  if (transmit) return;
70 
71  //real processes : sum track length
72  //
73  G4double stepLength = aStep->GetStepLength();
74  run->SumTrack(stepLength);
75 
76  //energy-momentum balance initialisation
77  //
78  const G4StepPoint* prePoint = aStep->GetPreStepPoint();
79  G4double Q = - prePoint->GetKineticEnergy();
80  G4ThreeVector Pbalance = - prePoint->GetMomentum();
81 
82  //initialisation of the nuclear channel identification
83  //
85  G4String partName = particle->GetParticleName();
86  G4String nuclearChannel = partName;
87  G4HadronicProcess* hproc = dynamic_cast<G4HadronicProcess*>(process);
88  const G4Isotope* target = NULL;
89  if (hproc) target = hproc->GetTargetIsotope();
90  G4String targetName = "XXXX";
91  if (target) targetName = target->GetName();
92  nuclearChannel += " + " + targetName + " --> ";
93  if (targetName == "XXXX") run->SetTargetXXX(true);
94 
95  //scattered primary particle (if any)
96  //
98  G4int ih = 1;
99  if (aStep->GetTrack()->GetTrackStatus() == fAlive) {
100  G4double energy = endPoint->GetKineticEnergy();
101  analysis->FillH1(ih,energy);
102  //
103  G4ThreeVector momentum = endPoint->GetMomentum();
104  Q += energy;
105  Pbalance += momentum;
106  //
107  nuclearChannel += partName + " + ";
108  }
109 
110  //secondaries
111  //
112  const std::vector<const G4Track*>* secondary
113  = aStep->GetSecondaryInCurrentStep();
114  for (size_t lp=0; lp<(*secondary).size(); lp++) {
115  particle = (*secondary)[lp]->GetDefinition();
116  G4String name = particle->GetParticleName();
117  G4String type = particle->GetParticleType();
118  G4double energy = (*secondary)[lp]->GetKineticEnergy();
119  run->ParticleCount(name,energy);
120  //energy spectrum
121  ih = 0;
122  if (particle == G4Gamma::Gamma()) ih = 2;
123  else if (particle == G4Neutron::Neutron()) ih = 3;
124  else if (particle == G4Proton::Proton()) ih = 4;
125  else if (particle == G4Deuteron::Deuteron()) ih = 5;
126  else if (particle == G4Alpha::Alpha()) ih = 6;
127  else if (type == "nucleus") ih = 7;
128  else if (type == "meson") ih = 8;
129  else if (type == "baryon") ih = 9;
130  if (ih > 0) analysis->FillH1(ih,energy);
131  //atomic mass
132  if (type == "nucleus") {
133  G4int A = particle->GetAtomicMass();
134  analysis->FillH1(12, A);
135  }
136  //energy-momentum balance
137  G4ThreeVector momentum = (*secondary)[lp]->GetMomentum();
138  Q += energy;
139  Pbalance += momentum;
140  //count e- from internal conversion together with gamma
141  if (particle == G4Electron::Electron()) particle = G4Gamma::Gamma();
142  //particle flag
144  }
145 
146  //energy-momentum balance
147  G4double Pbal = Pbalance.mag();
148  run->Balance(Pbal);
149  ih = 10;
150  analysis->FillH1(ih,Q);
151  ih = 11;
152  analysis->FillH1(ih,Pbal);
153 
154  // nuclear channel
155  const G4int kMax = 16;
156  const G4String conver[] = {"0","","2 ","3 ","4 ","5 ","6 ","7 ","8 ","9 ",
157  "10 ","11 ","12 ","13 ","14 ","15 ","16 "};
158  std::map<G4ParticleDefinition*,G4int>::iterator ip;
159  for (ip = fParticleFlag.begin(); ip != fParticleFlag.end(); ip++) {
160  particle = ip->first;
161  G4String name = particle->GetParticleName();
162  G4int nb = ip->second;
163  if (nb > kMax) nb = kMax;
164  G4String Nb = conver[nb];
165  if (particle == G4Gamma::Gamma()) {
166  run->CountGamma(nb);
167  Nb = "N ";
168  name = "gamma or e-";
169  }
170  if (ip != fParticleFlag.begin()) nuclearChannel += " + ";
171  nuclearChannel += Nb + name;
172  }
173 
175  run->CountNuclearChannel(nuclearChannel, Q);
176 
177  fParticleFlag.clear();
178 
179  // kill event after first interaction
180  //
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185