RockCutterAction.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Use Geant4's user "hooks" to kill particles in the rock
3 ///
4 /// \author rhatcher@fnal.gov
5 ////////////////////////////////////////////////////////////////////////
6 
7 // NOvA Includes
9 
10 #include "nug4/G4Base/PrimaryParticleInformation.h"
11 #include "nug4/G4Base/UserActionFactory.h"
12 USERACTIONREG3(g4n,RockCutterAction,g4n::RockCutterAction)
13 
14 // G4 includes
15 #include "Geant4/G4Event.hh"
16 #include "Geant4/G4Track.hh"
17 #include "Geant4/G4ThreeVector.hh"
18 #include "Geant4/G4ParticleDefinition.hh"
19 #include "Geant4/G4PrimaryParticle.hh"
20 #include "Geant4/G4DynamicParticle.hh"
21 #include "Geant4/G4VUserPrimaryParticleInformation.hh"
22 #include "Geant4/G4Step.hh"
23 #include "Geant4/G4StepPoint.hh"
24 #include "Geant4/G4VProcess.hh"
25 #include "Geant4/G4VPhysicalVolume.hh"
26 #include "Geant4/G4VTouchable.hh"
27 
28 // ROOT includes
29 #include <TLorentzVector.h>
30 #include <TString.h>
31 
32 // C/C++ includes
33 #include <algorithm>
34 #include <string>
35 
36 // ART includes
38 
39 namespace g4n {
40 
41  // Initialize static members.
42 
43  //-------------------------------------------------------------
44  // Constructor.
46  : fVerbose(0)
47  , fX0 (0)
48  , fY0 (0)
49  , fZ0 (728)
50  , fX_max (304.8)
51  , fX_min (-1620)
52  , fY_max (407.9)
53  , fY_min (-178.9)
54  , fZ_max (1800)
55  , fZ_min (-840)
56 
57  {
58 
59  }
60 
61  //-------------------------------------------------------------
62  // Destructor.
64  {
65  // Delete anything that we created with "new'.
66  }
67 
68  //-------------------------------------------------------------
70  {
71  fVerbose = pset.get< int >("Verbose",0);
72 
73  //Center of ND:
74  fX0 = pset.get< double >("X0",0);
75  fY0 = pset.get< double >("Y0",0);
76  fZ0 = pset.get< double >("Z0",728);
77 
78  //Grand Det Enclosure:
79  fX_max = pset.get< double >("X_max",304.8);
80  fY_max = pset.get< double >("Y_max",407.9);
81  fZ_max = pset.get< double >("Z_max",1800);
82  fX_min = pset.get< double >("X_min",-1620);
83  fY_min = pset.get< double >("Y_min",-178.9);
84  fZ_min = pset.get< double >("Z_min",-840);
85  if(fVerbose>0) PrintConfig("");
86  }
87 
88  //-------------------------------------------------------------
90  {
91  mf::LogInfo("RockCutterAction")
92  << "RockCutterAction::PrintConfig \n"
93  << " Verbose " << fVerbose << "\n"
94  << " fX0: " <<fX0 << " fY0: " <<fY0 << " fZ0:" << fZ0 << "\n"
95  << " fX_min: " <<fX_min << " fX_max: " << fX_max << "\n"
96  << " fY_min: " <<fY_min << " fY_max: " << fY_max << "\n"
97  << " fZ_min: " <<fZ_min << " fZ_max: " << fZ_max << "\n" ;
98 
99  }
100 
101  //-------------------------------------------------------------
102  G4ClassificationOfNewTrack
104  {
105 
106  if(fVerbose>10) std::cout << " in RockCutter Action " << std::endl;
107  //I have to do this to prevent duplicated trackIDes
108  //BUT, TIDIUS CODING, DO NOT MODIFY OTHER G4TRACK INFO!
109  G4Track *track = const_cast<G4Track*>(gtrack);
110 
111  const G4ThreeVector& g4pos = track->GetPosition(); // g4 in mm by default
112  // assuming we want positions in cm ...
113  const TVector3 startPos(g4pos.x()/CLHEP::cm,g4pos.y()/CLHEP::cm,g4pos.z()/CLHEP::cm);
114 
115  //Only do this to charged particles
116  if (fabs(track->GetParticleDefinition()->GetPDGCharge())<0.1) {
117 
118  //Define grand enclosure: the box containing the tunnel (of course partially)
119  //See if it is outside the grand detector enclosure:
120  if (startPos.X()>fX_max ||
121  startPos.X()<fX_min ||
122  startPos.Y()>fY_max ||
123  startPos.Y()<fY_min ||
124  startPos.Z()>fZ_max ||
125  startPos.Z()<fZ_min )
126  {
127 
128  //OK, surgery is needed:
129 
130  //Plan A: Fast
131 
132  const double energy = (track->GetTotalEnergy())/CLHEP::GeV;
133  const double dEdx = 4.4732e-3;
134  const G4ThreeVector& dir = track->GetMomentumDirection();
135  const TVector3 direction(dir.x(),dir.y(),dir.z());
136  if (startPos.Dot(direction)>0){
137 
138  // std::cout<<"wrong direction ID:"<<track->GetTrackID()<<std::endl;
139  track->SetTrackStatus(fStopAndKill);
140  }
141  else if (Distance_Wall(startPos.X(),startPos.Y(),startPos.Z())*dEdx>energy){
142 
143  // std::cout<<"not fast enough ID:"<<track->GetTrackID()<<std::endl;
144  track->SetTrackStatus(fStopAndKill);
145  }
146 
147  }
148  }
149 
150  return fUrgent;
151  }
152 
153 
154  double RockCutterAction::Distance_Wall(double x,double y,double z)
155  {
156  double distance = 0;
157  double dx = std::min(fabs(x-fX_max),fabs(x-fX_min));
158  double dy = std::min(fabs(y-fY_max),fabs(y-fY_min));
159  double dz = std::min(fabs(z-fZ_max),fabs(z-fZ_min));
160 
161  if (x<fX_max && x>fX_min && y<fY_max && y>fY_min)
162  distance = dz;
163 
164  else if (y<fY_max && y>-fY_min && z<fZ_max && z>fZ_min)
165  distance = dx;
166 
167  else if (z<fZ_max && z>-fZ_min && x<fX_max && x>fX_min)
168  distance = dy;
169 
170  else if (x<fX_max && x>fX_min)
171  distance = sqrt(dy*dy+dz*dz);
172 
173  else if (y<fY_max && y>fY_min)
174  distance = sqrt(dx*dx+dz*dz);
175 
176  else if (z<fZ_max && z>fZ_min)
177  distance = sqrt(dx*dx+dy*dy);
178 
179  else distance = sqrt(dx*dx+dy*dy+dz*dz);
180 
181  return distance;
182 
183  }
184 
185 
186  /*
187  //-------------------------------------------------------------
188  void RockCutterAction::BeginOfEventAction(const G4Event*)
189  {
190 
191  }
192 
193  //-------------------------------------------------------------
194  void RockCutterAction::EndOfEventAction(const G4Event*)
195  {
196 
197  }
198 
199  //-------------------------------------------------------------
200  void RockCutterAction::PreTrackingAction(const G4Track* track)
201  {
202  }
203 
204  //-------------------------------------------------------------
205  void RockCutterAction::PostTrackingAction( const G4Track* track)
206  {
207 
208  }
209 
210  //-------------------------------------------------------------
211  void RockCutterAction::SteppingAction(const G4Step* step)
212  {
213 
214  }
215 
216  //-------------------------------------------------------------
217  void RockCutterAction::StackNewStage()
218  {
219 
220  }
221 
222  //-------------------------------------------------------------
223  void RockCutterAction::StackPrepareNewEvent()
224  {
225 
226  }
227  */
228 
229 } // end namespace
TruthSlim – remove generated objects that don&#39;t contribute.
void Config(fhicl::ParameterSet const &pset)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
T sqrt(T number)
Definition: d0nt_math.hpp:156
Definition: event.h:19
double Distance_Wall(double, double, double)
unsigned distance(const T &t1, const T &t2)
static constexpr double cm
Definition: SystemOfUnits.h:99
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
void PrintConfig(std::string const &opt)
T get(std::string const &key) const
Definition: ParameterSet.h:231
double energy
Definition: plottest35.C:25
z
Definition: test.py:28
G4ClassificationOfNewTrack StackClassifyNewTrack(const G4Track *)
G4UserStackingAction interfaces.
OStream cout
Definition: OStream.cxx:6
TDirectory * dir
Definition: macro.C:5
T min(const caf::Proxy< T > &a, T b)
static constexpr double GeV
enum BeamMode string