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