| 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 electromagnetic/TestEm5/src/Tracking 26 /// \file electromagnetic/TestEm5/src/TrackingAction.cc
27 /// \brief Implementation of the TrackingActio 27 /// \brief Implementation of the TrackingAction 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 "TrackingAction.hh" 34 #include "TrackingAction.hh"
35 35
36 #include "DetectorConstruction.hh" 36 #include "DetectorConstruction.hh"
>> 37 #include "Run.hh"
37 #include "EventAction.hh" 38 #include "EventAction.hh"
38 #include "HistoManager.hh" 39 #include "HistoManager.hh"
39 #include "Run.hh" <<
40 40
41 #include "G4PhysicalConstants.hh" <<
42 #include "G4RunManager.hh" 41 #include "G4RunManager.hh"
43 #include "G4SystemOfUnits.hh" <<
44 #include "G4Track.hh" 42 #include "G4Track.hh"
>> 43 #include "G4PhysicalConstants.hh"
>> 44 #include "G4SystemOfUnits.hh"
45 45
46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 47
48 TrackingAction::TrackingAction(DetectorConstru << 48 TrackingAction::TrackingAction(DetectorConstruction* DET, EventAction* EA)
49 : fDetector(det), fEventAction(event) << 49 :G4UserTrackingAction(),fDetector(DET), fEventAction(EA)
50 {} << 50 {
51 << 51 fXstartAbs = fXendAbs = fPrimaryCharge = fDirX = 0.0;
>> 52 }
>> 53
52 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 55
54 void TrackingAction::PreUserTrackingAction(con << 56 void TrackingAction::PreUserTrackingAction(const G4Track* aTrack )
55 { 57 {
56 // few initialisations 58 // few initialisations
57 // 59 //
58 if (aTrack->GetTrackID() == 1) { 60 if (aTrack->GetTrackID() == 1) {
59 fXstartAbs = fDetector->GetxstartAbs(); 61 fXstartAbs = fDetector->GetxstartAbs();
60 fXendAbs = fDetector->GetxendAbs(); << 62 fXendAbs = fDetector->GetxendAbs();
61 fPrimaryCharge = aTrack->GetDefinition()-> 63 fPrimaryCharge = aTrack->GetDefinition()->GetPDGCharge();
62 fDirX = aTrack->GetMomentumDirection().x() 64 fDirX = aTrack->GetMomentumDirection().x();
63 } 65 }
64 } 66 }
65 67
66 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 69
68 void TrackingAction::PostUserTrackingAction(co 70 void TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
69 { 71 {
70 G4AnalysisManager* analysisManager = G4Analy 72 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
71 << 73
72 Run* run = static_cast<Run*>(G4RunManager::G << 74 Run* run = static_cast<Run*>(
>> 75 G4RunManager::GetRunManager()->GetNonConstCurrentRun());
73 76
74 G4ThreeVector position = aTrack->GetPosition 77 G4ThreeVector position = aTrack->GetPosition();
75 G4ThreeVector vertex = aTrack->GetVertexPosi << 78 G4ThreeVector vertex = aTrack->GetVertexPosition();
76 G4double charge = aTrack->GetDefinition()->G << 79 G4double charge = aTrack->GetDefinition()->GetPDGCharge();
77 G4double energy = aTrack->GetKineticEnergy() << 80 G4double energy = aTrack->GetKineticEnergy();
78 << 81
79 G4bool transmit = << 82 G4bool transmit = (((fDirX >= 0.0 && position.x() >= fXendAbs) ||
80 (((fDirX >= 0.0 && position.x() >= fXendAb << 83 (fDirX < 0.0 && position.x() <= fXstartAbs))
81 && energy > 0.0); << 84 && energy > 0.0);
82 G4bool reflect = << 85 G4bool reflect = (((fDirX >= 0.0 && position.x() <= fXstartAbs) ||
83 (((fDirX >= 0.0 && position.x() <= fXstart << 86 (fDirX < 0.0 && position.x() <= fXendAbs))
84 && energy > 0.0); << 87 && energy > 0.0);
85 G4bool notabsor = (transmit || reflect); 88 G4bool notabsor = (transmit || reflect);
86 89
87 // transmitted + reflected particles counter << 90 //transmitted + reflected particles counter
88 // 91 //
89 G4int flag = 0; 92 G4int flag = 0;
90 if (charge == fPrimaryCharge) flag = 1; << 93 if (charge == fPrimaryCharge) flag = 1;
91 if (aTrack->GetTrackID() == 1) flag = 2; 94 if (aTrack->GetTrackID() == 1) flag = 2;
92 if (transmit) fEventAction->SetTransmitFlag( 95 if (transmit) fEventAction->SetTransmitFlag(flag);
93 if (reflect) fEventAction->SetReflectFlag(fl << 96 if (reflect) fEventAction->SetReflectFlag(flag);
94 << 97
95 // 98 //
96 // histograms << 99 //histograms
97 // 100 //
98 G4bool charged = (charge != 0.); 101 G4bool charged = (charge != 0.);
99 G4bool neutral = !charged; 102 G4bool neutral = !charged;
100 103
101 // energy spectrum at exit << 104 //energy spectrum at exit
102 // zero energy charged particle are not take << 105 //zero energy charged particle are not taken into account
103 // << 106 //
104 G4int id = 0; << 107 G4int id = 0;
105 if (transmit && charged) << 108 if (transmit && charged) id = 10;
106 id = 10; << 109 else if (transmit && neutral) id = 20;
107 else if (transmit && neutral) << 110 else if (reflect && charged) id = 30;
108 id = 20; << 111 else if (reflect && neutral) id = 40;
109 else if (reflect && charged) <<
110 id = 30; <<
111 else if (reflect && neutral) <<
112 id = 40; <<
113 112
114 if (id > 0) { << 113 if (id>0) { analysisManager->FillH1(id, energy); }
115 analysisManager->FillH1(id, energy); <<
116 } <<
117 114
118 // energy leakage << 115 //energy leakage
119 // 116 //
120 if (notabsor) { 117 if (notabsor) {
121 G4int trackID = aTrack->GetTrackID(); 118 G4int trackID = aTrack->GetTrackID();
122 G4int index = 0; << 119 G4int index = 0; if (trackID > 1) index = 1; //primary=0, secondaries=1
123 if (trackID > 1) index = 1; // primary=0, <<
124 G4double eleak = aTrack->GetKineticEnergy( 120 G4double eleak = aTrack->GetKineticEnergy();
125 if ((aTrack->GetDefinition() == G4Positron 121 if ((aTrack->GetDefinition() == G4Positron::Positron()) && (trackID > 1))
126 eleak += 2 * electron_mass_c2; << 122 eleak += 2*electron_mass_c2;
127 run->AddEnergyLeak(eleak, index); << 123 run->AddEnergyLeak(eleak,index);
128 } 124 }
129 125
130 // space angle distribution at exit : dN/dOm << 126 //space angle distribution at exit : dN/dOmega
131 // 127 //
132 G4ThreeVector direction = aTrack->GetMomentu << 128 G4ThreeVector direction = aTrack->GetMomentumDirection();
133 id = 0; << 129 id = 0;
134 if (transmit && charged) << 130 if (transmit && charged) id = 12;
135 id = 12; << 131 else if (transmit && neutral) id = 22;
136 else if (transmit && neutral) << 132 else if (reflect && charged) id = 32;
137 id = 22; << 133 else if (reflect && neutral) id = 42;
138 else if (reflect && charged) <<
139 id = 32; <<
140 else if (reflect && neutral) <<
141 id = 42; <<
142 134
143 if (id > 0 && std::abs(direction.x()) < 1.0) << 135 if (id>0) {
144 G4double theta = std::acos(direction.x()); << 136 G4double theta = std::acos(direction.x());
145 if (theta > 0.0) { 137 if (theta > 0.0) {
146 G4double dteta = analysisManager->GetH1W << 138 G4double dteta = analysisManager->GetH1Width(id);
147 G4double unit = analysisManager->GetH1Un << 139 G4double unit = analysisManager->GetH1Unit(id);
148 G4double weight = unit / (twopi * std::s << 140 G4double weight = unit/(twopi*std::sin(theta)*dteta);
149 /* 141 /*
150 G4cout << "theta, dteta, unit, weight: " << 142 G4cout << "theta, dteta, unit, weight: "
151 << theta << " " << 143 << theta << " "
152 << dteta << " " << 144 << dteta << " "
153 << unit << " " << 145 << unit << " "
154 << weight << G4endl; << 146 << weight << G4endl;
155 */ 147 */
156 analysisManager->FillH1(id, theta, weigh << 148 analysisManager->FillH1(id,theta,weight);
157 } 149 }
158 } 150 }
159 << 151
160 // energy fluence at exit : dE(MeV)/dOmega << 152 //energy fluence at exit : dE(MeV)/dOmega
161 // 153 //
162 id = 0; << 154 id = 0;
163 if (transmit && charged) << 155 if (transmit && charged) id = 11;
164 id = 11; << 156 else if (reflect && charged) id = 31;
165 else if (reflect && charged) << 157 else if (transmit && neutral) id = 21;
166 id = 31; << 158 else if (reflect && neutral) id = 41;
167 else if (transmit && neutral) <<
168 id = 21; <<
169 else if (reflect && neutral) <<
170 id = 41; <<
171 159
172 if (id > 0 && std::abs(direction.x()) < 1.0) << 160 if (id>0) {
173 G4double theta = std::acos(direction.x()); << 161 G4double theta = std::acos(direction.x());
174 if (theta > 0.0) { 162 if (theta > 0.0) {
175 G4double dteta = analysisManager->GetH1W << 163 G4double dteta = analysisManager->GetH1Width(id);
176 G4double unit = analysisManager->GetH1Un << 164 G4double unit = analysisManager->GetH1Unit(id);
177 G4double weight = unit / (twopi * std::s << 165 G4double weight = unit/(twopi*std::sin(theta)*dteta);
178 weight *= (aTrack->GetKineticEnergy() / << 166 weight *= (aTrack->GetKineticEnergy()/MeV);
179 analysisManager->FillH1(id, theta, weigh << 167 analysisManager->FillH1(id,theta,weight);
180 } << 168 }
181 } << 169 }
>> 170
>> 171 //projected angles distribution at exit
>> 172 //
>> 173 id = 0;
>> 174 if (transmit && charged) id = 13;
>> 175 else if (transmit && neutral) id = 23;
>> 176 else if (reflect && charged) id = 33;
>> 177 else if (reflect && neutral) id = 43;
182 178
183 // projected angles distribution at exit << 179 if (id>0) {
184 // <<
185 id = 0; <<
186 if (transmit && charged) <<
187 id = 13; <<
188 else if (transmit && neutral) <<
189 id = 23; <<
190 else if (reflect && charged) <<
191 id = 33; <<
192 else if (reflect && neutral) <<
193 id = 43; <<
194 <<
195 if (id > 0) { <<
196 if (direction.x() != 0.0) { 180 if (direction.x() != 0.0) {
197 G4double tet = std::atan(direction.y() / << 181 G4double tet = std::atan(direction.y()/std::fabs(direction.x()));
198 analysisManager->FillH1(id, tet); << 182 analysisManager->FillH1(id,tet);
199 if (transmit && (flag == 2)) run->AddMsc 183 if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
200 184
201 tet = std::atan(direction.z() / std::fab << 185 tet = std::atan(direction.z()/std::fabs(direction.x()));
202 analysisManager->FillH1(id, tet); << 186 analysisManager->FillH1(id,tet);
203 if (transmit && (flag == 2)) run->AddMsc 187 if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
204 } 188 }
205 } 189 }
206 190
207 // projected position and radius at exit << 191 //projected position and radius at exit
208 // 192 //
209 if (transmit && energy > 0.0) { 193 if (transmit && energy > 0.0) {
210 G4double y = position.y(), z = position.z( 194 G4double y = position.y(), z = position.z();
211 G4double r = std::sqrt(y * y + z * z); << 195 G4double r = std::sqrt(y*y + z*z);
212 analysisManager->FillH1(14, y); << 196 analysisManager->FillH1(14, y);
213 analysisManager->FillH1(14, z); << 197 analysisManager->FillH1(14, z);
214 analysisManager->FillH1(15, r); << 198 analysisManager->FillH1(15, r);
215 } 199 }
216 << 200
217 // x-vertex of charged secondaries << 201 //x-vertex of charged secondaries
218 // 202 //
219 if ((aTrack->GetParentID() == 1) && charged) << 203 if ((aTrack->GetParentID() == 1) && charged ) {
220 G4double xVertex = (aTrack->GetVertexPosit 204 G4double xVertex = (aTrack->GetVertexPosition()).x();
221 analysisManager->FillH1(6, xVertex); 205 analysisManager->FillH1(6, xVertex);
222 if (notabsor) analysisManager->FillH1(7, x << 206 if (notabsor) analysisManager->FillH1(7, xVertex);
223 } 207 }
224 } 208 }
225 209
226 //....oooOO0OOooo........oooOO0OOooo........oo 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 211
227 212