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

Last change on this file since 1319 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
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: geant4-09-03-cand-01 $
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.