ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FCALSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FCALSteppingAction.cc
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 //
27 //
28 
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
31 
32 #include <iostream>
33 
34 #include "FCALSteppingAction.hh"
35 #include "G4SteppingManager.hh"
36 
37 #include "globals.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4Track.hh"
40 #include "G4DynamicParticle.hh"
41 #include "G4Material.hh"
42 
43 #include "G4LogicalVolume.hh"
44 #include "G4VPhysicalVolume.hh"
45 #include "G4VTouchable.hh"
46 #include "G4TouchableHistory.hh"
47 
48 #include "G4Event.hh"
49 
50 #include "G4ThreeVector.hh"
51 
52 #include "G4ios.hh"
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 
56 FCALSteppingAction::FCALSteppingAction():IDold(-1),IDout(-1)
57 {;}
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60 
61 FCALSteppingAction::~FCALSteppingAction()
62 {;}
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
66 void 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 == G4ThreeVector() &&
150  PrimaryDirection != G4ThreeVector() && SecondaryDirection != G4ThreeVector())
151  {
152  G4ThreeVector Temp = Distance.cross(PrimaryDirection);
153  VectorProduct = Temp.cross(PrimaryDirection);
154  };
155 
156 VectorProductMagnitude = VectorProduct.mag();
157  if(VectorProductMagnitude == 0.)
158  {
159  VectorProductNorm = G4ThreeVector();
160  } else {
161  VectorProductNorm = (1./VectorProduct.mag()) * VectorProduct ;
162  };
163  DistOfClosestApproach = Distance * VectorProductNorm ;
164 
165  if(std::abs(DistOfClosestApproach) < DCACut)
166  {
167  NSecondaries++;
168  Secondaries[0][0] = NSecondaries;
169  Secondaries[NSecondaries][1] = aParticle->GetPDGEncoding();
170  Secondaries[NSecondaries][2] = (aTrack->GetVertexPosition()).x();
171  Secondaries[NSecondaries][3] = (aTrack->GetVertexPosition()).y();
172  Secondaries[NSecondaries][4] = (aTrack->GetVertexPosition()).z();
173  Secondaries[NSecondaries][5] =(aDynamicParticle->GetMomentum()).x();
174  Secondaries[NSecondaries][6] = (aDynamicParticle->GetMomentum()).y();
175  Secondaries[NSecondaries][7] = (aDynamicParticle->GetMomentum()).z();
176  Secondaries[NSecondaries][8] = aDynamicParticle->GetTotalMomentum();
177  Secondaries[NSecondaries][9] = aDynamicParticle->GetTotalEnergy();
178  Secondaries[NSecondaries][10] =aDynamicParticle->GetKineticEnergy();
179  };
180  };
181  };
182 
183 
184  // Get the World leaving particle
185  if(aTrack->GetNextVolume() == 0) {
186  if(IDnow != IDout) {
187  IDout = IDnow;
188 
189  NTracks++;
190 
191  OutOfWorldTracksData[0][0] = NTracks;
192 
193  OutOfWorldTracksData[NTracks][1] = aParticle->GetPDGEncoding();
194 
195  OutOfWorldTracksData[NTracks][2] = (aTrack->GetVertexPosition()).x();
196  OutOfWorldTracksData[NTracks][3] = (aTrack->GetVertexPosition()).y();
197  OutOfWorldTracksData[NTracks][4] = (aTrack->GetVertexPosition()).z();
198 
199  OutOfWorldTracksData[NTracks][5] = (aDynamicParticle->GetMomentum()).x();
200  OutOfWorldTracksData[NTracks][6] = (aDynamicParticle->GetMomentum()).y();
201  OutOfWorldTracksData[NTracks][7] = (aDynamicParticle->GetMomentum()).z();
202 
203  OutOfWorldTracksData[NTracks][8] = aDynamicParticle->GetTotalMomentum();
204 
205  OutOfWorldTracksData[NTracks][9] = aDynamicParticle->GetTotalEnergy();
206 
207  OutOfWorldTracksData[NTracks][10] = aDynamicParticle->GetKineticEnergy();
208  };
209  };
210 
211 
212 }
213 
214 void FCALSteppingAction::initialize(G4int Nev) {
215  EventNo = Nev;
216  NTracks = 0;
217  NSecondaries = 0;
218  EdepFCALEm = EdepFCALHad = 0.;
219 
220  for(G4int i=0; i<6000; i++)
221  {
222  for(G4int j=0; j<11; j++)
223  {
224  OutOfWorldTracksData[i][j] = 0.;
225  Secondaries[i][j] = 0.;
226  }
227  };
228 }
229 
230 G4double FCALSteppingAction::GetOutOfWorldTracks(G4int i, G4int j){
231  return OutOfWorldTracksData[i][j];
232 }
233 
234 G4double FCALSteppingAction::GetSecondaries(G4int i, G4int j){
235  return Secondaries[i][j];
236 }
237 
238 G4double FCALSteppingAction::GetEdepFCAL(G4String FCAL) {
239  if(strcmp(FCAL,"FCALEm") == 0) {
240  return EdepFCALEm;
241  } else {
242  if(strcmp(FCAL,"FCALHad") == 0) {
243  return EdepFCALHad;}
244  }
245  return 0.0;
246 }
247 
248 
249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 
251 
252