| 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 hadronic/Hadr03/src/SteppingAction.cc
27 /// \brief Implementation of the SteppingActio 27 /// \brief Implementation of the SteppingAction class
28 // 28 //
29 // << 29 // $Id$
>> 30 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 33
33 #include "SteppingAction.hh" 34 #include "SteppingAction.hh"
34 << 35 #include "PrimaryGeneratorAction.hh"
>> 36 #include "RunAction.hh"
35 #include "HistoManager.hh" 37 #include "HistoManager.hh"
36 #include "Run.hh" <<
37 38
38 #include "G4HadronicProcess.hh" <<
39 #include "G4ParticleTypes.hh" 39 #include "G4ParticleTypes.hh"
40 #include "G4RunManager.hh" 40 #include "G4RunManager.hh"
>> 41 #include "G4HadronicProcess.hh"
>> 42
>> 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 44
>> 45 SteppingAction::SteppingAction(PrimaryGeneratorAction* prim, RunAction* RuAct)
>> 46 :fPrimary(prim),fRunAction(RuAct)
>> 47 { }
41 48
42 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 50
44 // SteppingAction::SteppingAction() << 51 SteppingAction::~SteppingAction()
45 //: G4UserSteppingAction() << 52 { }
46 //{ } <<
47 53
48 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 55
50 void SteppingAction::UserSteppingAction(const 56 void SteppingAction::UserSteppingAction(const G4Step* aStep)
51 { 57 {
52 // check trackID and stepNumber <<
53 G4int trackID = aStep->GetTrack()->GetTrackI <<
54 G4int stepNb = aStep->GetTrack()->GetCurrent <<
55 if (trackID * stepNb != 1) return; <<
56 // ok, we are at first interaction of the pr <<
57 <<
58 Run* run = static_cast<Run*>(G4RunManager::G <<
59 <<
60 // count processes 58 // count processes
61 // << 59 //
62 const G4StepPoint* endPoint = aStep->GetPost 60 const G4StepPoint* endPoint = aStep->GetPostStepPoint();
63 G4VProcess* process = const_cast<G4VProcess* << 61 const G4VProcess* process = endPoint->GetProcessDefinedStep();
64 run->CountProcesses(process); << 62 fRunAction->CountProcesses(process);
65 << 63
66 // check that an real interaction occured (e 64 // check that an real interaction occured (eg. not a transportation)
67 G4StepStatus stepStatus = endPoint->GetStepS 65 G4StepStatus stepStatus = endPoint->GetStepStatus();
68 G4bool transmit = (stepStatus == fGeomBounda << 66 G4bool transmit = (stepStatus==fGeomBoundary || stepStatus==fWorldBoundary);
69 if (transmit) return; 67 if (transmit) return;
70 << 68
71 // real processes : sum track length << 69 //real processes : sum track length
72 // 70 //
73 G4double stepLength = aStep->GetStepLength() 71 G4double stepLength = aStep->GetStepLength();
74 run->SumTrack(stepLength); << 72 fRunAction->SumTrack(stepLength);
75 << 73
76 // energy-momentum balance initialisation << 74 //energy-momentum balance initialisation
77 // 75 //
78 const G4StepPoint* prePoint = aStep->GetPreS 76 const G4StepPoint* prePoint = aStep->GetPreStepPoint();
79 G4double Q = -prePoint->GetKineticEnergy(); << 77 G4double Q = - prePoint->GetKineticEnergy();
80 G4ThreeVector Pbalance = -prePoint->GetMomen << 78 G4ThreeVector Pbalance = - prePoint->GetMomentum();
81 << 79
82 // initialisation of the nuclear channel ide << 80 //initialisation of the nuclear channel identification
83 // 81 //
84 G4ParticleDefinition* particle = aStep->GetT 82 G4ParticleDefinition* particle = aStep->GetTrack()->GetDefinition();
85 G4String partName = particle->GetParticleNam 83 G4String partName = particle->GetParticleName();
86 G4String nuclearChannel = partName; 84 G4String nuclearChannel = partName;
87 G4HadronicProcess* hproc = dynamic_cast<G4Ha << 85 G4HadronicProcess* hproc = (G4HadronicProcess*) process;
88 const G4Isotope* target = NULL; << 86 const G4Isotope* target = hproc->GetTargetIsotope();
89 if (hproc) target = hproc->GetTargetIsotope( << 87 G4String targetName = "XXXX";
90 G4String targetName = "XXXX"; <<
91 if (target) targetName = target->GetName(); 88 if (target) targetName = target->GetName();
92 nuclearChannel += " + " + targetName + " --> 89 nuclearChannel += " + " + targetName + " --> ";
93 if (targetName == "XXXX") run->SetTargetXXX( << 90
94 << 91 //scattered primary particle (if any)
95 // scattered primary particle (if any) <<
96 // 92 //
97 G4AnalysisManager* analysis = G4AnalysisMana 93 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
98 G4int ih = 1; 94 G4int ih = 1;
99 if (aStep->GetTrack()->GetTrackStatus() == f 95 if (aStep->GetTrack()->GetTrackStatus() == fAlive) {
100 G4double energy = endPoint->GetKineticEner << 96 G4double energy = endPoint->GetKineticEnergy();
101 analysis->FillH1(ih, energy); << 97 analysis->FillH1(ih,energy);
102 // 98 //
103 G4ThreeVector momentum = endPoint->GetMome 99 G4ThreeVector momentum = endPoint->GetMomentum();
104 Q += energy; << 100 Q += energy;
105 Pbalance += momentum; 101 Pbalance += momentum;
106 // 102 //
107 nuclearChannel += partName + " + "; 103 nuclearChannel += partName + " + ";
108 } << 104 }
109 << 105
110 // secondaries << 106 //secondaries
111 // << 107 //
112 const std::vector<const G4Track*>* secondary << 108 const G4TrackVector* secondary = fpSteppingManager->GetSecondary();
113 for (size_t lp = 0; lp < (*secondary).size() << 109 for (size_t lp=0; lp<(*secondary).size(); lp++) {
114 particle = (*secondary)[lp]->GetDefinition << 110 particle = (*secondary)[lp]->GetDefinition();
115 G4String name = particle->GetParticleName( << 111 G4String name = particle->GetParticleName();
116 G4String type = particle->GetParticleType( << 112 G4String type = particle->GetParticleType();
>> 113 G4double charge = particle->GetPDGCharge();
117 G4double energy = (*secondary)[lp]->GetKin 114 G4double energy = (*secondary)[lp]->GetKineticEnergy();
118 run->ParticleCount(name, energy); << 115 fRunAction->ParticleCount(name,energy);
119 // energy spectrum << 116 //energy spectrum
120 ih = 0; << 117 if (charge > 3.) ih = 2;
121 if (particle == G4Gamma::Gamma()) << 118 else if (particle == G4Gamma::Gamma()) ih = 3;
122 ih = 2; << 119 else if (particle == G4Neutron::Neutron()) ih = 4;
123 else if (particle == G4Electron::Electron( << 120 else if (particle == G4Proton::Proton()) ih = 5;
124 ih = 3; << 121 else if (particle == G4Deuteron::Deuteron()) ih = 6;
125 else if (particle == G4Neutron::Neutron()) << 122 else if (particle == G4Alpha::Alpha()) ih = 7;
126 ih = 4; << 123 else if (type == "nucleus") ih = 8;
127 else if (particle == G4Proton::Proton()) << 124 else if (type == "meson") ih = 9;
128 ih = 5; << 125 else if (type == "baryon") ih = 10;
129 else if (particle == G4Deuteron::Deuteron( << 126 analysis->FillH1(ih,energy);
130 ih = 6; << 127 //energy-momentum balance
131 else if (particle == G4Alpha::Alpha()) <<
132 ih = 7; <<
133 else if (type == "nucleus") <<
134 ih = 8; <<
135 else if (type == "meson") <<
136 ih = 9; <<
137 else if (type == "baryon") <<
138 ih = 10; <<
139 if (ih > 0) analysis->FillH1(ih, energy); <<
140 // atomic mass <<
141 if (type == "nucleus") { <<
142 G4int A = particle->GetAtomicMass(); <<
143 analysis->FillH1(13, A); <<
144 } <<
145 // energy-momentum balance <<
146 G4ThreeVector momentum = (*secondary)[lp]- 128 G4ThreeVector momentum = (*secondary)[lp]->GetMomentum();
147 Q += energy; << 129 Q += energy;
148 Pbalance += momentum; 130 Pbalance += momentum;
149 // count e- from internal conversion toget << 131 //particle flag
150 if (particle == G4Electron::Electron()) pa <<
151 // particle flag <<
152 fParticleFlag[particle]++; 132 fParticleFlag[particle]++;
153 } 133 }
154 << 134
155 // energy-momentum balance << 135 //energy-momentum balance
156 G4double Pbal = Pbalance.mag(); 136 G4double Pbal = Pbalance.mag();
157 run->Balance(Pbal); << 137 fRunAction->Balance(Pbal);
158 ih = 11; 138 ih = 11;
159 analysis->FillH1(ih, Q); << 139 analysis->FillH1(ih,Q);
160 ih = 12; 140 ih = 12;
161 analysis->FillH1(ih, Pbal); << 141 analysis->FillH1(ih,Pbal);
162 << 142
163 // nuclear channel 143 // nuclear channel
164 const G4int kMax = 16; << 144 const G4int kMax = 16;
165 const G4String conver[] = {"0", "", "2 " << 145 const G4String conver[] = {"0","","2 ","3 ","4 ","5 ","6 ","7 ","8 ","9 ",
166 "9 ", "10 ", "11 << 146 "10 ","11 ","12 ","13 ","14 ","15 ","16 "};
167 std::map<G4ParticleDefinition*, G4int>::iter << 147 std::map<G4ParticleDefinition*,G4int>::iterator ip;
168 for (ip = fParticleFlag.begin(); ip != fPart 148 for (ip = fParticleFlag.begin(); ip != fParticleFlag.end(); ip++) {
169 particle = ip->first; 149 particle = ip->first;
170 G4String name = particle->GetParticleName( << 150 G4String name = particle->GetParticleName();
171 G4int nb = ip->second; 151 G4int nb = ip->second;
172 if (nb > kMax) nb = kMax; << 152 if (nb > kMax) nb = kMax;
173 G4String Nb = conver[nb]; << 153 G4String Nb = conver[nb];
174 if (particle == G4Gamma::Gamma()) { 154 if (particle == G4Gamma::Gamma()) {
175 run->CountGamma(nb); << 155 fRunAction->CountGamma(nb);
176 Nb = "N "; << 156 Nb = "N ";
177 name = "gamma or e-"; << 157 }
178 } <<
179 if (ip != fParticleFlag.begin()) nuclearCh 158 if (ip != fParticleFlag.begin()) nuclearChannel += " + ";
180 nuclearChannel += Nb + name; 159 nuclearChannel += Nb + name;
181 } 160 }
182 << 161
183 /// G4cout << "\n nuclear channel: " << nucl << 162 ///G4cout << "\n nuclear channel: " << nuclearChannel << G4endl;
184 run->CountNuclearChannel(nuclearChannel, Q); << 163 fRunAction->CountNuclearChannel(nuclearChannel, Q);
185 << 164
186 fParticleFlag.clear(); 165 fParticleFlag.clear();
187 << 166
188 // kill event after first interaction 167 // kill event after first interaction
189 // 168 //
190 G4RunManager::GetRunManager()->AbortEvent(); << 169 G4RunManager::GetRunManager()->AbortEvent();
191 } 170 }
192 171
193 //....oooOO0OOooo........oooOO0OOooo........oo 172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 173
>> 174
194 175