| 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 // $Id: FCALSteppingAction.cc,v 1.7 2006/06/29 16:03:15 gunter Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-04 $
26 // 28 //
27 // 29 //
28 30
29 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
30 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
31 33
32 #include <iostream> <<
33 <<
34 #include "FCALSteppingAction.hh" 34 #include "FCALSteppingAction.hh"
35 #include "G4SteppingManager.hh" 35 #include "G4SteppingManager.hh"
36 36
37 #include "globals.hh" <<
38 #include "G4SystemOfUnits.hh" <<
39 #include "G4Track.hh" 37 #include "G4Track.hh"
40 #include "G4DynamicParticle.hh" 38 #include "G4DynamicParticle.hh"
41 #include "G4Material.hh" 39 #include "G4Material.hh"
42 40
43 #include "G4LogicalVolume.hh" 41 #include "G4LogicalVolume.hh"
44 #include "G4VPhysicalVolume.hh" 42 #include "G4VPhysicalVolume.hh"
45 #include "G4VTouchable.hh" 43 #include "G4VTouchable.hh"
46 #include "G4TouchableHistory.hh" 44 #include "G4TouchableHistory.hh"
47 45
48 #include "G4Event.hh" 46 #include "G4Event.hh"
49 47
50 #include "G4ThreeVector.hh" 48 #include "G4ThreeVector.hh"
51 49
52 #include "G4ios.hh" 50 #include "G4ios.hh"
>> 51 #include <iostream>
>> 52 #include "globals.hh"
53 53
54 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 55
56 FCALSteppingAction::FCALSteppingAction():IDold 56 FCALSteppingAction::FCALSteppingAction():IDold(-1),IDout(-1)
57 {;} 57 {;}
58 58
59 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60 60
61 FCALSteppingAction::~FCALSteppingAction() 61 FCALSteppingAction::~FCALSteppingAction()
62 {;} 62 {;}
63 63
64 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 65
66 void FCALSteppingAction::UserSteppingAction(co 66 void FCALSteppingAction::UserSteppingAction(const G4Step* astep)
67 { 67 {
68 // Get Edep 68 // Get Edep
69 G4double Edep = astep->GetTotalEnergyDeposi 69 G4double Edep = astep->GetTotalEnergyDeposit();
70 70
71 // Get Track 71 // Get Track
72 G4Track* aTrack = astep->GetTrack(); 72 G4Track* aTrack = astep->GetTrack();
73 73
74 // Get Touchable History 74 // Get Touchable History
75 G4TouchableHistory* theTouchable = (G4Touch 75 G4TouchableHistory* theTouchable = (G4TouchableHistory*)(aTrack->GetTouchable());
76 76
77 // Energy deposit in FCAL1 and FCAL2 77 // Energy deposit in FCAL1 and FCAL2
78 if(Edep != 0.) 78 if(Edep != 0.)
79 { 79 {
80 G4VPhysicalVolume* physVol = theTouchabl 80 G4VPhysicalVolume* physVol = theTouchable->GetVolume();
81 81
82 if(strcmp(physVol->GetName(),"FCALEmModu 82 if(strcmp(physVol->GetName(),"FCALEmModulePhysical")== 0 ||
83 strcmp(physVol->GetName(),"F1LArGapPhysical 83 strcmp(physVol->GetName(),"F1LArGapPhysical") == 0)
84 { 84 {
85 EdepFCALEm = EdepFCALEm + Edep; 85 EdepFCALEm = EdepFCALEm + Edep;
86 }; 86 };
87 87
88 if( (strcmp(physVol->GetName(), "FCALHad 88 if( (strcmp(physVol->GetName(), "FCALHadModulePhysical") == 0) ||
89 (strcmp(physVol->GetName(), "CuPlateAPhysi 89 (strcmp(physVol->GetName(), "CuPlateAPhysical") == 0) ||
90 (strcmp(physVol->GetName(), "CuPlateBPhysi 90 (strcmp(physVol->GetName(), "CuPlateBPhysical") == 0) ||
91 (strcmp(physVol->GetName(), "WAbsorberPhys 91 (strcmp(physVol->GetName(), "WAbsorberPhysical") == 0) ||
92 (strcmp(physVol->GetName(), "F2RodPhysical 92 (strcmp(physVol->GetName(), "F2RodPhysical") == 0) ||
93 (strcmp(physVol->GetName(), "F2LArGapPhysi 93 (strcmp(physVol->GetName(), "F2LArGapPhysical") == 0) )
94 { 94 {
95 EdepFCALHad = EdepFCALHad + Edep; 95 EdepFCALHad = EdepFCALHad + Edep;
96 }; 96 };
97 }; 97 };
98 98
99 // Get Tracks properties 99 // Get Tracks properties
100 G4int TrackID = aTrack->GetTrackID(); 100 G4int TrackID = aTrack->GetTrackID();
101 G4int ParentID = aTrack->GetParentID(); 101 G4int ParentID = aTrack->GetParentID();
102 // Get Associated particle 102 // Get Associated particle
103 const G4DynamicParticle * aDynamicParticle = 103 const G4DynamicParticle * aDynamicParticle = aTrack->GetDynamicParticle();
104 G4ParticleDefinition * aParticle = aTrack->G 104 G4ParticleDefinition * aParticle = aTrack->GetDefinition();
105 G4String ParticleName = aParticle->GetPartic 105 G4String ParticleName = aParticle->GetParticleName();
106 106
107 IDnow = EventNo + 10000*TrackID+ 100000000*P 107 IDnow = EventNo + 10000*TrackID+ 100000000*ParentID;
108 108
109 if(IDnow != IDold) 109 if(IDnow != IDold)
110 { 110 {
111 IDold = IDnow; 111 IDold = IDnow;
112 112
113 // Get the primary particle 113 // Get the primary particle
114 if(TrackID==1 && ParentID==0 && (aTrack- 114 if(TrackID==1 && ParentID==0 && (aTrack->GetCurrentStepNumber()) == 1)
115 { 115 {
116 PrimaryVertex = aTrack->GetVertexPositi 116 PrimaryVertex = aTrack->GetVertexPosition();
117 PrimaryDirection = aTrack->GetVertexMoment 117 PrimaryDirection = aTrack->GetVertexMomentumDirection();
118 118
119 NSecondaries = 1; 119 NSecondaries = 1;
120 Secondaries[NSecondaries][1] = aParticle-> 120 Secondaries[NSecondaries][1] = aParticle->GetPDGEncoding();
121 Secondaries[NSecondaries][2] = PrimaryVert 121 Secondaries[NSecondaries][2] = PrimaryVertex.x();
122 Secondaries[NSecondaries][3] = PrimaryVert 122 Secondaries[NSecondaries][3] = PrimaryVertex.y();
123 Secondaries[NSecondaries][4] = PrimaryVert 123 Secondaries[NSecondaries][4] = PrimaryVertex.z();
124 Secondaries[NSecondaries][5] = (aDynamicPa 124 Secondaries[NSecondaries][5] = (aDynamicParticle->GetMomentum()).x();
125 Secondaries[NSecondaries][6] = (aDynamicPa 125 Secondaries[NSecondaries][6] = (aDynamicParticle->GetMomentum()).y();
126 Secondaries[NSecondaries][7] = (aDynamicPa 126 Secondaries[NSecondaries][7] = (aDynamicParticle->GetMomentum()).z();
127 Secondaries[NSecondaries][8] = aDynamicPar 127 Secondaries[NSecondaries][8] = aDynamicParticle->GetTotalMomentum();
128 Secondaries[NSecondaries][9] = aDynamicPar 128 Secondaries[NSecondaries][9] = aDynamicParticle->GetTotalEnergy();
129 Secondaries[NSecondaries][10] = aDynamicPa 129 Secondaries[NSecondaries][10] = aDynamicParticle->GetKineticEnergy();
130 130
131 G4cout << " **** Primary : " << EventNo < 131 G4cout << " **** Primary : " << EventNo << G4endl;
132 G4cout << " Vertex : " << PrimaryVertex << 132 G4cout << " Vertex : " << PrimaryVertex << G4endl;
133 } 133 }
134 134
135 135
136 // Get secondaries in air close to the p 136 // Get secondaries in air close to the primary tracks (DCA < 2.mm)
137 G4double DCACut = 2.*mm; 137 G4double DCACut = 2.*mm;
138 G4String Material = aTrack->GetMaterial( 138 G4String Material = aTrack->GetMaterial()->GetName();
139 G4ThreeVector TrackPos = aTrack->GetVert 139 G4ThreeVector TrackPos = aTrack->GetVertexPosition();
140 140
141 if(TrackID != 1 && ParentID == 1 && (str 141 if(TrackID != 1 && ParentID == 1 && (strcmp(Material,"Air")==0) && (TrackPos.z() > 135.*cm))
142 { 142 {
143 SecondaryVertex = aTrack->GetVertexPositio 143 SecondaryVertex = aTrack->GetVertexPosition();
144 SecondaryDirection = aTrack->GetVertexMome 144 SecondaryDirection = aTrack->GetVertexMomentumDirection();
145 145
146 // calculate DCA of secondries to primary 146 // calculate DCA of secondries to primary particle
147 Distance = PrimaryVertex - SecondaryVertex 147 Distance = PrimaryVertex - SecondaryVertex ;
148 VectorProduct = PrimaryDirection.cross(Sec 148 VectorProduct = PrimaryDirection.cross(SecondaryDirection);
149 if(VectorProduct == G4ThreeVector() && << 149 if(VectorProduct == 0. &&
150 PrimaryDirection != G4ThreeVector() && << 150 PrimaryDirection != 0. && SecondaryDirection != 0.)
151 { 151 {
152 G4ThreeVector Temp = Distance.cross(Pr 152 G4ThreeVector Temp = Distance.cross(PrimaryDirection);
153 VectorProduct = Temp.cross(PrimaryDire 153 VectorProduct = Temp.cross(PrimaryDirection);
154 }; 154 };
155 << 155 VectorProductMagnitude = VectorProduct.mag();
156 VectorProductMagnitude = VectorProduct.mag(); <<
157 if(VectorProductMagnitude == 0.) 156 if(VectorProductMagnitude == 0.)
158 { 157 {
159 VectorProductNorm = G4ThreeVector(); << 158 VectorProductNorm = 0.;
160 } else { 159 } else {
161 VectorProductNorm = (1./VectorProduct. 160 VectorProductNorm = (1./VectorProduct.mag()) * VectorProduct ;
162 }; 161 };
163 DistOfClosestApproach = Distance * VectorP 162 DistOfClosestApproach = Distance * VectorProductNorm ;
164 163
165 if(std::abs(DistOfClosestApproach) < DCACu 164 if(std::abs(DistOfClosestApproach) < DCACut)
166 { 165 {
167 NSecondaries++; 166 NSecondaries++;
168 Secondaries[0][0] = NSecondaries; 167 Secondaries[0][0] = NSecondaries;
169 Secondaries[NSecondaries][1] = aPartic 168 Secondaries[NSecondaries][1] = aParticle->GetPDGEncoding();
170 Secondaries[NSecondaries][2] = (aTrack 169 Secondaries[NSecondaries][2] = (aTrack->GetVertexPosition()).x();
171 Secondaries[NSecondaries][3] = (aTrack 170 Secondaries[NSecondaries][3] = (aTrack->GetVertexPosition()).y();
172 Secondaries[NSecondaries][4] = (aTrack 171 Secondaries[NSecondaries][4] = (aTrack->GetVertexPosition()).z();
173 Secondaries[NSecondaries][5] =(aDynami 172 Secondaries[NSecondaries][5] =(aDynamicParticle->GetMomentum()).x();
174 Secondaries[NSecondaries][6] = (aDynam 173 Secondaries[NSecondaries][6] = (aDynamicParticle->GetMomentum()).y();
175 Secondaries[NSecondaries][7] = (aDynam 174 Secondaries[NSecondaries][7] = (aDynamicParticle->GetMomentum()).z();
176 Secondaries[NSecondaries][8] = aDynami 175 Secondaries[NSecondaries][8] = aDynamicParticle->GetTotalMomentum();
177 Secondaries[NSecondaries][9] = aDynami 176 Secondaries[NSecondaries][9] = aDynamicParticle->GetTotalEnergy();
178 Secondaries[NSecondaries][10] =aDynami 177 Secondaries[NSecondaries][10] =aDynamicParticle->GetKineticEnergy();
179 }; 178 };
180 }; 179 };
181 }; 180 };
182 181
183 182
184 // Get the World leaving particle 183 // Get the World leaving particle
185 if(aTrack->GetNextVolume() == 0) { 184 if(aTrack->GetNextVolume() == 0) {
186 if(IDnow != IDout) { 185 if(IDnow != IDout) {
187 IDout = IDnow; 186 IDout = IDnow;
188 187
189 NTracks++; 188 NTracks++;
190 189
191 OutOfWorldTracksData[0][0] = NTracks; 190 OutOfWorldTracksData[0][0] = NTracks;
192 191
193 OutOfWorldTracksData[NTracks][1] = aPart 192 OutOfWorldTracksData[NTracks][1] = aParticle->GetPDGEncoding();
194 193
195 OutOfWorldTracksData[NTracks][2] = (aTra 194 OutOfWorldTracksData[NTracks][2] = (aTrack->GetVertexPosition()).x();
196 OutOfWorldTracksData[NTracks][3] = (aTra 195 OutOfWorldTracksData[NTracks][3] = (aTrack->GetVertexPosition()).y();
197 OutOfWorldTracksData[NTracks][4] = (aTra 196 OutOfWorldTracksData[NTracks][4] = (aTrack->GetVertexPosition()).z();
198 197
199 OutOfWorldTracksData[NTracks][5] = (aDyn 198 OutOfWorldTracksData[NTracks][5] = (aDynamicParticle->GetMomentum()).x();
200 OutOfWorldTracksData[NTracks][6] = (aDyn 199 OutOfWorldTracksData[NTracks][6] = (aDynamicParticle->GetMomentum()).y();
201 OutOfWorldTracksData[NTracks][7] = (aDyn 200 OutOfWorldTracksData[NTracks][7] = (aDynamicParticle->GetMomentum()).z();
202 201
203 OutOfWorldTracksData[NTracks][8] = aDyna 202 OutOfWorldTracksData[NTracks][8] = aDynamicParticle->GetTotalMomentum();
204 203
205 OutOfWorldTracksData[NTracks][9] = aDyna 204 OutOfWorldTracksData[NTracks][9] = aDynamicParticle->GetTotalEnergy();
206 205
207 OutOfWorldTracksData[NTracks][10] = aDyn 206 OutOfWorldTracksData[NTracks][10] = aDynamicParticle->GetKineticEnergy();
208 }; 207 };
209 }; 208 };
210 209
211 210
212 } 211 }
213 212
214 void FCALSteppingAction::initialize(G4int Nev) 213 void FCALSteppingAction::initialize(G4int Nev) {
215 EventNo = Nev; 214 EventNo = Nev;
216 NTracks = 0; 215 NTracks = 0;
217 NSecondaries = 0; 216 NSecondaries = 0;
218 EdepFCALEm = EdepFCALHad = 0.; 217 EdepFCALEm = EdepFCALHad = 0.;
219 218
220 for(G4int i=0; i<6000; i++) 219 for(G4int i=0; i<6000; i++)
221 { 220 {
222 for(G4int j=0; j<11; j++) 221 for(G4int j=0; j<11; j++)
223 { 222 {
224 OutOfWorldTracksData[i][j] = 0.; 223 OutOfWorldTracksData[i][j] = 0.;
225 Secondaries[i][j] = 0.; 224 Secondaries[i][j] = 0.;
226 } 225 }
227 }; 226 };
228 } 227 }
229 228
230 G4double FCALSteppingAction::GetOutOfWorldTrac 229 G4double FCALSteppingAction::GetOutOfWorldTracks(G4int i, G4int j){
231 return OutOfWorldTracksData[i][j]; 230 return OutOfWorldTracksData[i][j];
232 } 231 }
233 232
234 G4double FCALSteppingAction::GetSecondaries(G4 233 G4double FCALSteppingAction::GetSecondaries(G4int i, G4int j){
235 return Secondaries[i][j]; 234 return Secondaries[i][j];
236 } 235 }
237 236
238 G4double FCALSteppingAction::GetEdepFCAL(G4Str 237 G4double FCALSteppingAction::GetEdepFCAL(G4String FCAL) {
239 if(strcmp(FCAL,"FCALEm") == 0) { 238 if(strcmp(FCAL,"FCALEm") == 0) {
240 return EdepFCALEm; 239 return EdepFCALEm;
241 } else { 240 } else {
242 if(strcmp(FCAL,"FCALHad") == 0) { 241 if(strcmp(FCAL,"FCALHad") == 0) {
243 return EdepFCALHad;} 242 return EdepFCALHad;}
244 } 243 }
245 return 0.0; 244 return 0.0;
246 } 245 }
247 246
248 247
249 //....oooOO0OOooo........oooOO0OOooo........oo 248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 249
251 250
252 251
253 252