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