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 //
26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 //
34 
35 #include "Analysis.hh"
36 
37 #include "SteppingAction.hh"
38 #include "RunAction.hh"
39 #include "DetectorConstruction.hh"
40 #include "PrimaryGeneratorAction.hh"
41 
42 #include "G4SystemOfUnits.hh"
43 #include "G4SteppingManager.hh"
44 #include "G4VTouchable.hh"
45 #include "G4VPhysicalVolume.hh"
46 #include "CommandLineParser.hh"
47 
48 #include "G4Electron.hh"
49 #include "G4Proton.hh"
50 #include "G4Gamma.hh"
51 #include "G4Alpha.hh"
52 #include "G4DNAGenericIonsManager.hh"
53 #include "G4EventManager.hh"
54 #include "G4Event.hh"
55 
56 using namespace G4DNAPARSER;
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 SteppingAction::SteppingAction() : G4UserSteppingAction()
61 {}
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
65 SteppingAction::~SteppingAction()
66 {}
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
70 void SteppingAction::UserSteppingAction(const G4Step* step)
71 {
72  G4double flagParticle=-1;
73  G4double flagProcess=-1;
74  G4double x,y,z,xp,yp,zp;
75 
76  // Particle identification
77 
78  // The following method avoids the usage of string comparison
79 
80  G4ParticleDefinition * partDef =
81  step->GetTrack()->GetDynamicParticle()->GetDefinition();
82 
83  if (partDef == G4Gamma::GammaDefinition())
84  flagParticle = 0;
85 
86  if (partDef == G4Electron::ElectronDefinition())
87  flagParticle = 1;
88 
89  if (partDef == G4Proton::ProtonDefinition())
90  flagParticle = 2;
91 
92  if (partDef == G4Alpha::AlphaDefinition())
93  flagParticle = 4;
94 
95  G4DNAGenericIonsManager *instance;
96  instance = G4DNAGenericIonsManager::Instance();
97 
98  // Usage example
99  /*
100  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
101  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
102  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
103  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
104  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
105  */
106 
107  if (partDef == instance->GetIon("hydrogen"))
108  flagParticle = 3;
109 
110  if (partDef == instance->GetIon("alpha+"))
111  flagParticle = 5;
112 
113  if (partDef == instance->GetIon("helium"))
114  flagParticle = 6;
115 
116  // Alternative method (based on string comparison)
117 
118  /*
119  const G4String& particleName = step->GetTrack()->GetDynamicParticle()->
120  GetDefinition()->GetParticleName();
121 
122  if (particleName == "gamma") flagParticle = 0;
123  else if (particleName == "e-") flagParticle = 1;
124  else if (particleName == "proton") flagParticle = 2;
125  else if (particleName == "hydrogen") flagParticle = 3;
126  else if (particleName == "alpha") flagParticle = 4;
127  else if (particleName == "alpha+") flagParticle = 5;
128  else if (particleName == "helium") flagParticle = 6;
129  */
130 
131  // Process identification
132 
133  // Process sub-types are listed in G4PhysicsListHelper.cc
134  // or in Geant4-DNA process class implementation files (*.cc)
135 
136  G4StepPoint * preStep = step->GetPreStepPoint();
137  G4StepPoint * postStep = step->GetPostStepPoint();
138  G4int procID = postStep->GetProcessDefinedStep()->GetProcessSubType();
139 
140  const G4String& processName = postStep->
141  GetProcessDefinedStep()->GetProcessName();
142 
143  if (processName=="eCapture") flagProcess =1;
144  // (no subType and procID exists at the moment for this process)
145 
146  else if (flagParticle == 1)
147 
148  {
149  if (procID==58) flagProcess =10;
150  else if (procID==51) flagProcess =11;
151  else if (procID==52) flagProcess =12;
152  else if (procID==53) flagProcess =13;
153  else if (procID==55) flagProcess =14;
154  else if (procID==54) flagProcess =15;
155  else if (procID==10) flagProcess =110;
156  else if (procID==1) flagProcess =120;
157  else if (procID==2) flagProcess =130;
158  }
159 
160  else if (flagParticle == 2)
161 
162  {
163  if (procID==51) flagProcess =21;
164  else if (procID==52) flagProcess =22;
165  else if (procID==53) flagProcess =23;
166  else if (procID==56) flagProcess =24;
167  else if (procID==10) flagProcess =210;
168  else if (procID==1) flagProcess =220;
169  else if (procID==2) flagProcess =230;
170  else if (procID==8) flagProcess =240;
171  }
172 
173  else if (flagParticle == 3)
174 
175  {
176  if (procID==51) flagProcess =31;
177  else if (procID==52) flagProcess =32;
178  else if (procID==53) flagProcess =33;
179  else if (procID==57) flagProcess =35;
180  }
181 
182  else if (flagParticle == 4)
183 
184  {
185  if (procID==51) flagProcess =41;
186  else if (procID==52) flagProcess =42;
187  else if (procID==53) flagProcess =43;
188  else if (procID==56) flagProcess =44;
189  else if (procID==10) flagProcess =410;
190  else if (procID==1) flagProcess =420;
191  else if (procID==2) flagProcess =430;
192  else if (procID==8) flagProcess =440;
193  }
194 
195  else if (flagParticle == 5)
196 
197  {
198  if (procID==51) flagProcess =51;
199  else if (procID==52) flagProcess =52;
200  else if (procID==53) flagProcess =53;
201  else if (procID==56) flagProcess =54;
202  else if (procID==57) flagProcess =55;
203  else if (procID==10) flagProcess =510;
204  else if (procID==1) flagProcess =520;
205  else if (procID==2) flagProcess =530;
206  else if (procID==8) flagProcess =540;
207  }
208 
209  else if (flagParticle == 6)
210 
211  {
212  if (procID==51) flagProcess =61;
213  else if (procID==52) flagProcess =62;
214  else if (procID==53) flagProcess =63;
215  else if (procID==57) flagProcess =65;
216  }
217 
218  else if (processName=="GenericIon_G4DNAIonisation") flagProcess =73;
219  else if (processName=="msc") flagProcess =710;
220  else if (processName=="CoulombScat") flagProcess =720;
221  else if (processName=="ionIoni") flagProcess =730;
222  else if (processName=="nuclearStopping") flagProcess =740;
223  // (for all GenericIons)
224 
225  // Alternatively, using process names
226 
227  /*
228  if (processName=="e-_G4DNAElectronSolvation") flagProcess =10;
229  else if (processName=="e-_G4DNAElastic") flagProcess =11;
230  else if (processName=="e-_G4DNAExcitation") flagProcess =12;
231  else if (processName=="e-_G4DNAIonisation") flagProcess =13;
232  else if (processName=="e-_G4DNAAttachment") flagProcess =14;
233  else if (processName=="e-_G4DNAVibExcitation") flagProcess =15;
234  else if (processName=="eCapture") flagProcess =16;
235  else if (processName=="msc") flagProcess =17;
236  else if (processName=="eIoni") flagProcess =130;
237 
238  else if (processName=="proton_G4DNAElastic") flagProcess =21;
239  else if (processName=="proton_G4DNAExcitation") flagProcess =22;
240  else if (processName=="proton_G4DNAIonisation") flagProcess =23;
241  else if (processName=="proton_G4DNAChargeDecrease") flagProcess =24;
242  else if (processName=="hIoni") flagProcess =230;
243 
244  else if (processName=="hydrogen_G4DNAElastic") flagProcess =31;
245  else if (processName=="hydrogen_G4DNAExcitation") flagProcess =32;
246  else if (processName=="hydrogen_G4DNAIonisation") flagProcess =33;
247  else if (processName=="hydrogen_G4DNAChargeIncrease") flagProcess =35;
248 
249  else if (processName=="alpha_G4DNAElastic") flagProcess =41;
250  else if (processName=="alpha_G4DNAExcitation") flagProcess =42;
251  else if (processName=="alpha_G4DNAIonisation") flagProcess =43;
252  else if (processName=="alpha_G4DNAChargeDecrease") flagProcess =44;
253 
254  else if (processName=="alpha+_G4DNAElastic") flagProcess =51;
255  else if (processName=="alpha+_G4DNAExcitation") flagProcess =52;
256  else if (processName=="alpha+_G4DNAIonisation") flagProcess =53;
257  else if (processName=="alpha+_G4DNAChargeDecrease") flagProcess =54;
258  else if (processName=="alpha+_G4DNAChargeIncrease") flagProcess =55;
259 
260  else if (processName=="helium_G4DNAElastic") flagProcess =61;
261  else if (processName=="helium_G4DNAExcitation") flagProcess =62;
262  else if (processName=="helium_G4DNAIonisation") flagProcess =63;
263  else if (processName=="helium_G4DNAChargeIncrease") flagProcess =65;
264 
265  else if (processName=="GenericIon_G4DNAIonisation") flagProcess =73;
266 */
267 
268 
269  if (processName!="Transportation")
270  {
271  x=preStep->GetPosition().x()/nanometer;
272  y=preStep->GetPosition().y()/nanometer;
273  z=preStep->GetPosition().z()/nanometer;
274 
275  xp=postStep->GetPosition().x()/nanometer;
276  yp=postStep->GetPosition().y()/nanometer;
277  zp=postStep->GetPosition().z()/nanometer;
278 
279  // get analysis manager
280  CommandLineParser* parser = CommandLineParser::GetParser();
281  Command* command(0);
282  if((command = parser->GetCommandIfActive("-out"))==0) return;
283 
284  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
285 
286  // fill ntuple
287  analysisManager->FillNtupleDColumn(0, flagParticle);
288  analysisManager->FillNtupleDColumn(1, flagProcess);
289  analysisManager->FillNtupleDColumn(2, x);
290  analysisManager->FillNtupleDColumn(3, y);
291  analysisManager->FillNtupleDColumn(4, z);
292  analysisManager->FillNtupleDColumn(5,
293  step->GetTotalEnergyDeposit()/eV);
294 
295  analysisManager->FillNtupleDColumn(6,
296  std::sqrt((x-xp)*(x-xp)+
297  (y-yp)*(y-yp)+(z-zp)*(z-zp)));
298 
299  analysisManager->FillNtupleDColumn(7,
300  (preStep->
301  GetKineticEnergy() -
302  postStep->
303  GetKineticEnergy())/eV );
304 
305  analysisManager->FillNtupleDColumn(8, preStep->
306  GetKineticEnergy()/eV);
307 
308  analysisManager->FillNtupleDColumn(9,
309  preStep->GetMomentumDirection()
310  *postStep->GetMomentumDirection() );
311 
312  analysisManager->FillNtupleIColumn(10,
313  G4EventManager::GetEventManager()->
314  GetConstCurrentEvent()->GetEventID());
315 
316  analysisManager->FillNtupleIColumn(11, step->GetTrack()->GetTrackID());
317 
318  analysisManager->FillNtupleIColumn(12, step->GetTrack()->GetParentID());
319 
320  analysisManager->FillNtupleIColumn(13,
321  step->GetTrack()->GetCurrentStepNumber());
322 
323  analysisManager->AddNtupleRow();
324  }
325 }