ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SteppingAction.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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "SteppingAction.hh"
34 #include "DetectorConstruction.hh"
35 #include "RunAction.hh"
36 #include "HistoManager.hh"
37 #include "G4ParticleTypes.hh"
38 
39 #include "G4RunManager.hh"
40 
41 #include <G4ThreeVector.hh>
42 #include <G4RotationMatrix.hh>
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
47  RunAction* RuAct)
48 :G4UserSteppingAction(),fDetector(det), fRunAction(RuAct)
49 { }
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
54 { }
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
59 {
60  G4StepPoint* prePoint = aStep->GetPreStepPoint();
61 
62  // if World --> return
63  if(prePoint->GetTouchableHandle()->GetVolume()==fDetector->GetWorld()) return;
64 
65  // here we enter in the absorber Box
66  // tag the event to be killed anyway after this step
67  //
69 
70  //count processes and keep only Multiple Scattering or gamma converion
71  //
72  G4StepPoint* endPoint = aStep->GetPostStepPoint();
73  G4String procName = endPoint->GetProcessDefinedStep()->GetProcessName();
74  fRunAction->CountProcesses(procName);
75 
76  if (procName == "msc" || procName == "muMsc" || procName == "stepMax") {
77 
78  //below, only multiple Scattering happens
79  //
80  G4ThreeVector position = endPoint->GetPosition();
81  G4ThreeVector direction = endPoint->GetMomentumDirection();
82 
83  G4double truePathLength = aStep->GetStepLength();
84  G4double geomPathLength = position.x() + 0.5*fDetector->GetBoxSize();
85  G4double ratio = geomPathLength/truePathLength;
86  fRunAction->SumPathLength(truePathLength,geomPathLength);
87  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
88  analysisManager->FillH1(1,truePathLength);
89  analysisManager->FillH1(2,geomPathLength);
90  analysisManager->FillH1(3,ratio);
91 
92  G4double yend = position.y(), zend = position.z();
93  G4double lateralDisplacement = std::sqrt(yend*yend + zend*zend);
94  fRunAction->SumLateralDisplacement(lateralDisplacement);
95  analysisManager->FillH1(4,lateralDisplacement);
96 
97  G4double psi = std::atan(lateralDisplacement/geomPathLength);
98  fRunAction->SumPsi(psi);
99  analysisManager->FillH1(5,psi);
100 
101  G4double xdir = direction.x(), ydir = direction.y(), zdir = direction.z();
102  G4double tetaPlane = std::atan2(ydir, xdir);
103  fRunAction->SumTetaPlane(tetaPlane);
104  analysisManager->FillH1(6,tetaPlane);
105  tetaPlane = std::atan2(zdir, xdir);
106  fRunAction->SumTetaPlane(tetaPlane);
107  analysisManager->FillH1(6,tetaPlane);
108 
109  G4double phiPos = std::atan2(zend, yend);
110  analysisManager->FillH1(7,phiPos);
111  G4double phiDir = std::atan2(zdir, ydir);
112  analysisManager->FillH1(8,phiDir);
113 
114  G4double phiCorrel = 0.;
115  if (lateralDisplacement > 0.)
116  phiCorrel = (yend*ydir + zend*zdir)/lateralDisplacement;
117  fRunAction->SumPhiCorrel(phiCorrel);
118  analysisManager->FillH1(9,phiCorrel);
119  } else if (procName == "conv" || procName == "GammaToMuPair" ) {
120 
121  // gamma conversion
122 
123  G4StepPoint* PrePoint = aStep->GetPreStepPoint();
124  G4double EGamma = PrePoint->GetTotalEnergy();
125  G4ThreeVector PGamma = PrePoint->GetMomentum();
126  G4ThreeVector PolaGamma = PrePoint->GetPolarization();
127 
128  G4double Eplus=-1;
129  G4ThreeVector Pplus, Pminus, Precoil;
130 
131  const G4TrackVector* secondary = fpSteppingManager->GetSecondary();
132 
133  const size_t Nsecondaries = (*secondary).size();
134 
135  //No conversion , E < threshold
136  if (Nsecondaries == 0) return;
137 
138  for (size_t lp=0; lp< std::min(Nsecondaries,size_t(2) ); lp++) {
139  if (((*secondary)[lp]->GetDefinition()==G4Electron::Definition())
140  || ((*secondary)[lp]->GetDefinition()==G4MuonMinus::Definition()) )
141  {
142  Pminus = (*secondary)[lp]->GetMomentum();
143  }
144  if (((*secondary)[lp]->GetDefinition()==G4Positron::Definition())
145  || ((*secondary)[lp]->GetDefinition()==G4MuonPlus::Definition()) )
146  {
147  Eplus = (*secondary)[lp]->GetTotalEnergy();
148  Pplus = (*secondary)[lp]->GetMomentum();
149  }
150  }
151 
152  if ( Nsecondaries >= 3 ) {
153  Precoil = (*secondary)[2]->GetMomentum();
154  }
155 
156  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
157 
158  // Fill Histograms
159 
160  G4ThreeVector z = PGamma.unit(); // gamma direction
161  G4ThreeVector x(1.,0.,0.);
162 
163  // pola perpendicular to direction
164 
165  if ( PolaGamma.mag() != 0.0 ) {
166  x = PolaGamma.unit();
167  } else { // Pola = 0 case
168  x = z.orthogonal().unit();
169  }
170 
171  G4ThreeVector y = z;
172  y = y.cross(x);
173 
174  G4RotationMatrix GtoW(x,y,z); // from gamma ref. sys. to World
175  G4RotationMatrix WtoG = inverseOf(GtoW); // from World to gamma ref. sys.
176 
177 
178  G4double angleE = Pplus.angle(Pminus) * EGamma;
179  analysisManager->FillH1(10,angleE);
180 
181  if ( Nsecondaries >= 3 ) {
182  // recoil returned
183  analysisManager->FillH1(11,std::log10(Precoil.mag()));
184  analysisManager->FillH1(12,Precoil.transform(WtoG).phi());
185  }
186  G4double phiPlus = Pplus.transform(WtoG).phi();
187  G4double phiMinus = Pminus.transform(WtoG).phi();
188  analysisManager->FillH1(13,phiPlus);
189  analysisManager->FillH1(14,std::cos(phiPlus + phiMinus) * -2.0);
190  analysisManager->FillH1(15,Eplus/EGamma);
191 
192  G4double phiPola = PolaGamma.transform(WtoG).phi();
193  analysisManager->FillH1(16, phiPola);
194  }
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......