source: trunk/examples/advanced/lAr_calorimeter/src/FCALSteppingAction.cc@ 1194

Last change on this file since 1194 was 807, checked in by garnier, 17 years ago

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