CMS 3D CMS Logo

SteppingAction.cc
Go to the documentation of this file.
1 
5 
6 #include "G4LogicalVolumeStore.hh"
7 #include "G4ParticleTable.hh"
8 #include "G4PhysicalVolumeStore.hh"
9 #include "G4RegionStore.hh"
10 #include "G4UnitsTable.hh"
11 #include "G4SystemOfUnits.hh"
12 
14 
15 //#define DebugLog
16 
18  const CMSSteppingVerbose* sv, bool hasW)
19  : eventAction_(e), tracker(nullptr), calo(nullptr), steppingVerbose(sv),
20  initialized(false), killBeamPipe(false),hasWatcher(hasW)
21 {
23  (p.getParameter<double>("CriticalEnergyForVacuum")*CLHEP::MeV);
24  if(0.0 < theCriticalEnergyForVacuum) { killBeamPipe = true; }
26  (p.getParameter<double>("CriticalDensity")*CLHEP::g/CLHEP::cm3);
27  maxTrackTime = p.getParameter<double>("MaxTrackTime")*CLHEP::ns;
28  maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
29  maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
30  deadRegionNames = p.getParameter<std::vector<std::string> >("DeadRegions");
31  ekinMins = p.getParameter<std::vector<double> >("EkinThresholds");
32  ekinNames = p.getParameter<std::vector<std::string> >("EkinNames");
33  ekinParticles = p.getParameter<std::vector<std::string> >("EkinParticles");
34 
35  edm::LogInfo("SimG4CoreApplication") << "SteppingAction:: KillBeamPipe = "
36  << killBeamPipe << " CriticalDensity = "
37  << theCriticalDensity*CLHEP::cm3/CLHEP::g
38  << " g/cm3;"
39  << " CriticalEnergyForVacuum = "
41  << " Mev;"
42  << " MaxTrackTime = "
43  << maxTrackTime/CLHEP::ns
44  << " ns";
45 
46  numberTimes = maxTrackTimes.size();
47  if(numberTimes > 0) {
48  for (unsigned int i=0; i<numberTimes; i++) {
49  edm::LogInfo("SimG4CoreApplication") << "SteppingAction::MaxTrackTime for "
50  << maxTimeNames[i] << " is "
51  << maxTrackTimes[i] << " ns ";
52  maxTrackTimes[i] *= ns;
53  }
54  }
55 
56  ndeadRegions = deadRegionNames.size();
57  if(ndeadRegions > 0) {
58  edm::LogInfo("SimG4CoreApplication")
59  << "SteppingAction: Number of DeadRegions where all trackes are killed "
60  << ndeadRegions;
61  for(unsigned int i=0; i<ndeadRegions; ++i) {
62  edm::LogInfo("SimG4CoreApplication")
63  << "SteppingAction: DeadRegion " << i << ". " << deadRegionNames[i];
64  }
65  }
66  numberEkins = ekinNames.size();
67  numberPart = ekinParticles.size();
68  if(0 == numberPart) { numberEkins = 0; }
69 
70  if(numberEkins > 0) {
71 
72  edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Kill following "
73  << numberPart
74  << " particles in " << numberEkins
75  << " volumes";
76  for (unsigned int i=0; i<numberPart; ++i) {
77  edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Particle " << i
78  << " " << ekinParticles[i]
79  << " Threshold = " << ekinMins[i]
80  << " MeV";
81  ekinMins[i] *= CLHEP::MeV;
82  }
83  for (unsigned int i=0; i<numberEkins; ++i) {
84  edm::LogInfo("SimG4CoreApplication") << "SteppingAction::LogVolume[" << i
85  << "] = " << ekinNames[i];
86  }
87  }
88 }
89 
91 
92 void SteppingAction::UserSteppingAction(const G4Step * aStep)
93 {
94  if (!initialized) { initialized = initPointer(); }
95 
96  // if(hasWatcher) { m_g4StepSignal(aStep); }
97  m_g4StepSignal(aStep);
98 
99  G4Track * theTrack = aStep->GetTrack();
100  TrackStatus tstat = (theTrack->GetTrackStatus() == fAlive) ? sAlive : sKilledByProcess;
101  G4StepPoint* postStep = aStep->GetPostStepPoint();
102 
103  if(0 == tstat && postStep->GetPhysicalVolume() != nullptr) {
104 
105  G4StepPoint* preStep = aStep->GetPreStepPoint();
106  const G4Region* theRegion =
107  preStep->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
108 
109  // kill in dead regions
110  if(isInsideDeadRegion(theRegion)) { tstat = sDeadRegion; }
111 
112  // kill out of time
113  if(0 == tstat && isOutOfTimeWindow(theTrack, theRegion)) { tstat = sOutOfTime; }
114 
115  // kill low-energy in volumes on demand
116  if(0 == tstat && numberEkins > 0 && isLowEnergy(aStep)) { tstat = sLowEnergy; }
117 
118  // kill low-energy in vacuum
119  G4double kinEnergy = theTrack->GetKineticEnergy();
120  if(0 == tstat && killBeamPipe && kinEnergy < theCriticalEnergyForVacuum
121  && theTrack->GetDefinition()->GetPDGCharge() != 0.0 && kinEnergy > 0.0
122  && theTrack->GetNextVolume()->GetLogicalVolume()->GetMaterial()->GetDensity()
123  <= theCriticalDensity) {
124  tstat = sLowEnergyInVacuum;
125  }
126 
127  // check transition tracker/calo
128  if(0 == tstat) {
129 
130  if(isThisVolume(preStep->GetTouchable(),tracker) &&
131  isThisVolume(postStep->GetTouchable(),calo)) {
132 
133  math::XYZVectorD pos((preStep->GetPosition()).x(),
134  (preStep->GetPosition()).y(),
135  (preStep->GetPosition()).z());
136 
137  math::XYZTLorentzVectorD mom((preStep->GetMomentum()).x(),
138  (preStep->GetMomentum()).y(),
139  (preStep->GetMomentum()).z(),
140  preStep->GetTotalEnergy());
141 
142  uint32_t id = theTrack->GetTrackID();
143 
144  std::pair<math::XYZVectorD,math::XYZTLorentzVectorD> p(pos,mom);
146  }
147  } else {
148  theTrack->SetTrackStatus(fStopAndKill);
149 #ifdef DebugLog
150  PrintKilledTrack(theTrack, tstat);
151 #endif
152  }
153  }
154  if(nullptr != steppingVerbose) {
155  steppingVerbose->NextStep(aStep, fpSteppingManager, (1 < tstat));
156  }
157 }
158 
159 bool SteppingAction::isLowEnergy(const G4Step * aStep) const
160 {
161  bool flag = false;
162  const G4StepPoint* sp = aStep->GetPostStepPoint();
163  G4LogicalVolume* lv = sp->GetPhysicalVolume()->GetLogicalVolume();
164  for (unsigned int i=0; i<numberEkins; ++i) {
165  if (lv == ekinVolumes[i]) {
166  flag = true;
167  break;
168  }
169  }
170  if (flag) {
171  double ekin = sp->GetKineticEnergy();
172  int pCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
173  for (unsigned int i=0; i<numberPart; ++i) {
174  if (pCode == ekinPDG[i]) {
175  return (ekin <= ekinMins[i]) ? true : false;
176  }
177  }
178  }
179  return false;
180 }
181 
183 {
184  const G4PhysicalVolumeStore * pvs = G4PhysicalVolumeStore::GetInstance();
185  if (pvs) {
186  std::vector<G4VPhysicalVolume*>::const_iterator pvcite;
187  for (pvcite = pvs->begin(); pvcite != pvs->end(); ++pvcite) {
188  if ((*pvcite)->GetName() == "Tracker") tracker = (*pvcite);
189  if ((*pvcite)->GetName() == "CALO") calo = (*pvcite);
190  if (tracker && calo) break;
191  }
192  if (tracker || calo) {
193  edm::LogInfo("SimG4CoreApplication") << "Pointer for Tracker " << tracker
194  << " and for Calo " << calo;
195  if (tracker) LogDebug("SimG4CoreApplication") << "Tracker vol name "
196  << tracker->GetName();
197  if (calo) LogDebug("SimG4CoreApplication") << "Calorimeter vol name "
198  << calo->GetName();
199  }
200  }
201 
202  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
203  if (numberEkins > 0) {
204  if (lvs) {
205  ekinVolumes.resize(numberEkins, 0);
206  std::vector<G4LogicalVolume*>::const_iterator lvcite;
207  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite) {
208  for (unsigned int i=0; i<numberEkins; ++i) {
209  if ((*lvcite)->GetName() == (G4String)(ekinNames[i])) {
210  ekinVolumes[i] = (*lvcite);
211  break;
212  }
213  }
214  }
215  }
216  for (unsigned int i=0; i<numberEkins; ++i) {
217  edm::LogInfo("SimG4CoreApplication") << ekinVolumes[i]->GetName()
218  <<" with pointer " << ekinVolumes[i];
219  }
220  }
221 
222  if(numberPart > 0) {
223  G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
224  G4String partName;
225  ekinPDG.resize(numberPart, 0);
226  for (unsigned int i=0; i<numberPart; ++i) {
227  ekinPDG[i] =
228  theParticleTable->FindParticle(partName=ekinParticles[i])->GetPDGEncoding();
229  edm::LogInfo("SimG4CoreApplication") << "Particle " << ekinParticles[i]
230  << " with PDG code " << ekinPDG[i]
231  << " and KE cut off "
232  << ekinMins[i]/MeV << " MeV";
233  }
234  }
235 
236  const G4RegionStore * rs = G4RegionStore::GetInstance();
237  if (numberTimes > 0) {
238  maxTimeRegions.resize(numberTimes, 0);
239  std::vector<G4Region*>::const_iterator rcite;
240  for (rcite = rs->begin(); rcite != rs->end(); ++rcite) {
241  for (unsigned int i=0; i<numberTimes; ++i) {
242  if ((*rcite)->GetName() == (G4String)(maxTimeNames[i])) {
243  maxTimeRegions[i] = (*rcite);
244  break;
245  }
246  }
247  }
248  }
249  if (ndeadRegions > 0) {
250  deadRegions.resize(ndeadRegions, 0);
251  std::vector<G4Region*>::const_iterator rcite;
252  for (rcite = rs->begin(); rcite != rs->end(); ++rcite) {
253  for (unsigned int i=0; i<ndeadRegions; ++i) {
254  if ((*rcite)->GetName() == (G4String)(deadRegionNames[i])) {
255  deadRegions[i] = (*rcite);
256  break;
257  }
258  }
259  }
260  }
261  return true;
262 }
263 
264 void SteppingAction::PrintKilledTrack(const G4Track* aTrack,
265  const TrackStatus& tst) const
266 {
267  std::string vname = "";
268  std::string rname = "";
269  std::string typ = " ";
270  switch (tst) {
271  case sKilledByProcess:
272  typ = " G4Process ";
273  break;
274  case sDeadRegion:
275  typ = " in dead region ";
276  break;
277  case sOutOfTime:
278  typ = " out of time window ";
279  break;
280  case sLowEnergy:
281  typ = " low energy limit ";
282  break;
283  case sLowEnergyInVacuum:
284  typ = " low energy limit in vacuum ";
285  break;
286  default:
287  break;
288  }
289  G4VPhysicalVolume* pv = aTrack->GetNextVolume();
290  if(pv) {
291  vname = pv->GetLogicalVolume()->GetName();
292  rname = pv->GetLogicalVolume()->GetRegion()->GetName();
293  }
294 
295  edm::LogInfo("SimG4CoreApplication")
296  << "Track #" << aTrack->GetTrackID()
297  << " " << aTrack->GetDefinition()->GetParticleName()
298  << " E(MeV)= " << aTrack->GetKineticEnergy()/MeV
299  << " T(ns)= " << aTrack->GetGlobalTime()/ns
300  << " is killed due to " << typ
301  << " inside LV: " << vname << " (" << rname
302  << ") at " << aTrack->GetPosition();
303 }
#define LogDebug(id)
virtual ~SteppingAction()
T getParameter(std::string const &) const
virtual void UserSteppingAction(const G4Step *aStep) final
const G4VPhysicalVolume * tracker
std::vector< int > ekinPDG
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
std::vector< const G4Region * > deadRegions
double theCriticalDensity
SimActivityRegistry::G4StepSignal m_g4StepSignal
void NextStep(const G4Step *, const G4SteppingManager *ptr, bool isKilled) const
unsigned int numberPart
const CMSSteppingVerbose * steppingVerbose
bool isOutOfTimeWindow(G4Track *theTrack, const G4Region *reg) const
unsigned int numberEkins
#define nullptr
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
const double MeV
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
Definition: Vector3D.h:8
TrackStatus
std::vector< std::string > deadRegionNames
const G4VPhysicalVolume * calo
def pv(vc)
Definition: MetAnalyzer.py:6
unsigned int numberTimes
std::vector< double > maxTrackTimes
std::vector< G4LogicalVolume * > ekinVolumes
double theCriticalEnergyForVacuum
EventAction * eventAction_
std::vector< double > ekinMins
std::vector< std::string > maxTimeNames
std::vector< std::string > ekinNames
void addTkCaloStateInfo(uint32_t t, const std::pair< math::XYZVectorD, math::XYZTLorentzVectorD > &p)
Definition: EventAction.cc:74
bool isLowEnergy(const G4Step *aStep) const
bool isInsideDeadRegion(const G4Region *reg) const
SteppingAction(EventAction *ea, const edm::ParameterSet &ps, const CMSSteppingVerbose *, bool hasW)
void PrintKilledTrack(const G4Track *, const TrackStatus &) const
unsigned int ndeadRegions
std::vector< std::string > ekinParticles
bool isThisVolume(const G4VTouchable *touch, const G4VPhysicalVolume *pv) const
std::vector< const G4Region * > maxTimeRegions