ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackingAction.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 //
28 //
29 //
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "TrackingAction.hh"
35 
36 #include "DetectorConstruction.hh"
37 #include "Run.hh"
38 #include "EventAction.hh"
39 #include "HistoManager.hh"
40 
41 #include "G4RunManager.hh"
42 #include "G4Track.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
49 :G4UserTrackingAction(),fDetector(DET), fEventAction(EA)
50 {
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
57 {
58  // few initialisations
59  //
60  if (aTrack->GetTrackID() == 1) {
64  fDirX = aTrack->GetMomentumDirection().x();
65  }
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
71 {
73 
74  Run* run = static_cast<Run*>(
76 
78  G4ThreeVector vertex = aTrack->GetVertexPosition();
80  G4double energy = aTrack->GetKineticEnergy();
81 
82  G4bool transmit = (((fDirX >= 0.0 && position.x() >= fXendAbs) ||
83  (fDirX < 0.0 && position.x() <= fXstartAbs))
84  && energy > 0.0);
85  G4bool reflect = (((fDirX >= 0.0 && position.x() <= fXstartAbs) ||
86  (fDirX < 0.0 && position.x() <= fXendAbs))
87  && energy > 0.0);
88  G4bool notabsor = (transmit || reflect);
89 
90  //transmitted + reflected particles counter
91  //
92  G4int flag = 0;
93  if (charge == fPrimaryCharge) flag = 1;
94  if (aTrack->GetTrackID() == 1) flag = 2;
95  if (transmit) fEventAction->SetTransmitFlag(flag);
96  if (reflect) fEventAction->SetReflectFlag(flag);
97 
98  //
99  //histograms
100  //
101  G4bool charged = (charge != 0.);
102  G4bool neutral = !charged;
103 
104  //energy spectrum at exit
105  //zero energy charged particle are not taken into account
106  //
107  G4int id = 0;
108  if (transmit && charged) id = 10;
109  else if (transmit && neutral) id = 20;
110  else if (reflect && charged) id = 30;
111  else if (reflect && neutral) id = 40;
112 
113  if (id>0) { analysisManager->FillH1(id, energy); }
114 
115  //energy leakage
116  //
117  if (notabsor) {
118  G4int trackID = aTrack->GetTrackID();
119  G4int index = 0; if (trackID > 1) index = 1; //primary=0, secondaries=1
120  G4double eleak = aTrack->GetKineticEnergy();
121  if ((aTrack->GetDefinition() == G4Positron::Positron()) && (trackID > 1))
122  eleak += 2*electron_mass_c2;
123  run->AddEnergyLeak(eleak,index);
124  }
125 
126  //space angle distribution at exit : dN/dOmega
127  //
128  G4ThreeVector direction = aTrack->GetMomentumDirection();
129  id = 0;
130  if (transmit && charged) id = 12;
131  else if (transmit && neutral) id = 22;
132  else if (reflect && charged) id = 32;
133  else if (reflect && neutral) id = 42;
134 
135  if (id>0) {
136  G4double theta = std::acos(direction.x());
137  if (theta > 0.0) {
138  G4double dteta = analysisManager->GetH1Width(id);
139  G4double unit = analysisManager->GetH1Unit(id);
140  G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
141  /*
142  G4cout << "theta, dteta, unit, weight: "
143  << theta << " "
144  << dteta << " "
145  << unit << " "
146  << weight << G4endl;
147  */
148  analysisManager->FillH1(id,theta,weight);
149  }
150  }
151 
152  //energy fluence at exit : dE(MeV)/dOmega
153  //
154  id = 0;
155  if (transmit && charged) id = 11;
156  else if (reflect && charged) id = 31;
157  else if (transmit && neutral) id = 21;
158  else if (reflect && neutral) id = 41;
159 
160  if (id>0) {
161  G4double theta = std::acos(direction.x());
162  if (theta > 0.0) {
163  G4double dteta = analysisManager->GetH1Width(id);
164  G4double unit = analysisManager->GetH1Unit(id);
165  G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
166  weight *= (aTrack->GetKineticEnergy()/MeV);
167  analysisManager->FillH1(id,theta,weight);
168  }
169  }
170 
171  //projected angles distribution at exit
172  //
173  id = 0;
174  if (transmit && charged) id = 13;
175  else if (transmit && neutral) id = 23;
176  else if (reflect && charged) id = 33;
177  else if (reflect && neutral) id = 43;
178 
179  if (id>0) {
180  if (direction.x() != 0.0) {
181  G4double tet = std::atan(direction.y()/std::fabs(direction.x()));
182  analysisManager->FillH1(id,tet);
183  if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
184 
185  tet = std::atan(direction.z()/std::fabs(direction.x()));
186  analysisManager->FillH1(id,tet);
187  if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
188  }
189  }
190 
191  //projected position and radius at exit
192  //
193  if (transmit && energy > 0.0) {
194  G4double y = position.y(), z = position.z();
195  G4double r = std::sqrt(y*y + z*z);
196  analysisManager->FillH1(14, y);
197  analysisManager->FillH1(14, z);
198  analysisManager->FillH1(15, r);
199  }
200 
201  //x-vertex of charged secondaries
202  //
203  if ((aTrack->GetParentID() == 1) && charged && energy > 0.0) {
204  G4double xVertex = (aTrack->GetVertexPosition()).x();
205  analysisManager->FillH1(6, xVertex);
206  if (notabsor) analysisManager->FillH1(7, xVertex);
207  }
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
211