CMS 3D CMS Logo
CMSSteppingVerbose.cc
Go to the documentation of this file.
1 
2 #include "SimG4Core/Notification/interface/CMSSteppingVerbose.h"
3 #include "G4Event.hh"
4 #include "G4EventManager.hh"
5 #include "G4Track.hh"
6 #include "G4TrackStatus.hh"
7 #include "G4Step.hh"
8 #include "G4SteppingManager.hh"
9 #include "G4SteppingVerbose.hh"
10 #include "G4ParticleDefinition.hh"
11 #include "G4VProcess.hh"
12 #include "G4SystemOfUnits.hh"
13 
14 CMSSteppingVerbose::CMSSteppingVerbose(
15  G4int verb, G4double ekin, std::vector<G4int>& evtNum, std::vector<G4int>& primV, std::vector<G4int>& trNum)
16  : m_PrintEvent(false),
17  m_PrintTrack(false),
18  m_smInitialized(false),
19  m_verbose(verb),
20  m_EventNumbers(evtNum),
21  m_PrimaryVertex(primV),
22  m_TrackNumbers(trNum),
23  m_EkinThreshold(ekin) {
24  m_nEvents = m_EventNumbers.size();
25  m_nVertex = m_PrimaryVertex.size();
26  m_nTracks = m_TrackNumbers.size();
27  m_g4SteppingVerbose = new G4SteppingVerbose();
28  G4VSteppingVerbose::SetInstance(m_g4SteppingVerbose);
29  m_g4SteppingVerbose->SetSilent(1);
30 }
31 
32 void CMSSteppingVerbose::beginOfEvent(const G4Event* evt) {
33  m_PrintEvent = false;
34  if (0 >= m_verbose) {
35  return;
36  }
37  if (m_nEvents == 0) {
38  m_PrintEvent = true;
39  } else {
40  for (G4int i = 0; i < m_nEvents; ++i) {
41  // check event number
42  if (evt->GetEventID() == m_EventNumbers[i]) {
43  // check number of vertex
44  if (m_nVertex == m_nEvents && evt->GetNumberOfPrimaryVertex() != m_PrimaryVertex[i]) {
45  continue;
46  }
47  m_PrintEvent = true;
48  break;
49  }
50  }
51  }
52  if (!m_PrintEvent) {
53  return;
54  }
55  G4cout << "========== Event #" << evt->GetEventID() << " " << evt->GetNumberOfPrimaryVertex()
56  << " primary vertexes ======" << G4endl;
57  G4cout << G4endl;
58 }
59 
60 void CMSSteppingVerbose::stopEventPrint() {
61  m_PrintEvent = false;
62  m_PrintTrack = false;
63  m_verbose = 0;
64 }
65 
66 void CMSSteppingVerbose::trackStarted(const G4Track* track, bool isKilled) {
67  m_PrintTrack = false;
68  if (!m_PrintEvent) {
69  return;
70  }
71 
72  if (!m_smInitialized) {
73  G4SteppingManager* stepman = G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
74  m_g4SteppingVerbose->SetManager(stepman);
75  stepman->SetVerboseLevel(m_verbose);
76  m_smInitialized = true;
77  }
78 
79  if (m_nTracks == 0) {
80  if (track->GetKineticEnergy() >= m_EkinThreshold) {
81  m_PrintTrack = true;
82  }
83 
84  } else {
85  for (G4int i = 0; i < m_nTracks; ++i) {
86  if (track->GetTrackID() == m_TrackNumbers[i]) {
87  m_PrintTrack = true;
88  break;
89  }
90  }
91  }
92  if (!m_PrintTrack) {
93  return;
94  }
95 
96  G4cout << "*********************************************************************************************************"
97  << G4endl;
98  const G4ParticleDefinition* pd = track->GetDefinition();
99  G4cout << "* G4Track Information: Particle = ";
100  if (pd) {
101  G4cout << pd->GetParticleName();
102  }
103  G4cout << ", Track ID = " << track->GetTrackID() << ", Parent ID = " << track->GetParentID() << G4endl;
104  G4cout << "*********************************************************************************************************"
105  << G4endl;
106 
107  G4cout << std::setw(5) << "Step#"
108  << " " << std::setw(8) << "X(cm)"
109  << " " << std::setw(8) << "Y(cm)"
110  << " " << std::setw(8) << "Z(cm)"
111  << " " << std::setw(9) << "KinE(GeV)"
112  << " " << std::setw(8) << "dE(MeV)"
113  << " " << std::setw(8) << "Step(mm)"
114  << " " << std::setw(9) << "TrackL(cm)"
115  << " " << std::setw(30) << "PhysVolume"
116  << " " << std::setw(8) << "ProcName" << G4endl;
117 
118  G4int prec = G4cout.precision(4);
119 
120  G4cout << std::setw(5) << track->GetCurrentStepNumber() << " " << std::setw(8) << track->GetPosition().x() / CLHEP::cm
121  << " " << std::setw(8) << track->GetPosition().y() / CLHEP::cm << " " << std::setw(8)
122  << track->GetPosition().z() / CLHEP::cm << " " << std::setw(9) << track->GetKineticEnergy() / CLHEP::GeV << " "
123  << std::setw(8) << " " << std::setw(8) << " " << std::setw(9) << " ";
124  if (track->GetVolume() != nullptr) {
125  G4cout << std::setw(30) << track->GetVolume()->GetName() << " ";
126  }
127  if (isKilled) {
128  G4cout << "isKilled";
129  }
130  G4cout << G4endl;
131  G4cout.precision(prec);
132 }
133 
134 void CMSSteppingVerbose::stackFilled(const G4Track* track, bool isKilled) const {
135  if (2 >= m_verbose || !m_PrintTrack || track->GetKineticEnergy() < m_EkinThreshold) {
136  return;
137  }
138  G4int prec = G4cout.precision(4);
139 
140  G4cout << std::setw(10) << track->GetTrackID() << " " << std::setw(8) << track->GetPosition().x() / CLHEP::cm << " "
141  << std::setw(8) << track->GetPosition().y() / CLHEP::cm << " " << std::setw(8)
142  << track->GetPosition().z() / CLHEP::cm << " " << std::setw(9) << track->GetKineticEnergy() / CLHEP::GeV
143  << " ";
144  if (track->GetVolume() != nullptr) {
145  G4cout << std::setw(24) << track->GetVolume()->GetName() << " ";
146  }
147  if (isKilled) {
148  G4cout << "isKilled";
149  }
150  G4cout << G4endl;
151  G4cout.precision(prec);
152 }
153 
154 void CMSSteppingVerbose::nextStep(const G4Step* step, const G4SteppingManager* sManager, bool isKilled) const {
155  if (!m_PrintTrack) {
156  return;
157  }
158 
159  G4int prec;
160  const G4Track* track = step->GetTrack();
161  const G4StepPoint* preStep = step->GetPreStepPoint();
162  const G4StepPoint* postStep = step->GetPostStepPoint();
163 
164  if (3 <= m_verbose) {
165  m_g4SteppingVerbose->SetSilent(0);
166  m_g4SteppingVerbose->DPSLStarted();
167  m_g4SteppingVerbose->DPSLAlongStep();
168  m_g4SteppingVerbose->DPSLPostStep();
169  if (4 <= m_verbose) {
170  m_g4SteppingVerbose->AlongStepDoItAllDone();
171  m_g4SteppingVerbose->PostStepDoItAllDone();
172  }
173  m_g4SteppingVerbose->SetSilent(1);
174 
175  prec = G4cout.precision(16);
176 
177  G4cout << G4endl;
178  G4cout << " ++G4Step Information " << G4endl;
179  G4cout << " Address of G4Track : " << track << G4endl;
180  G4cout << " Step Length (mm) : " << track->GetStepLength() << G4endl;
181  G4cout << " Energy Deposit (MeV) : " << step->GetTotalEnergyDeposit() << G4endl;
182 
183  G4cout << " -------------------------------------------------------"
184  << "-------------------------------" << G4endl;
185  G4cout << " StepPoint Information " << std::setw(30) << "PreStep" << std::setw(30) << "PostStep" << G4endl;
186  G4cout << " -------------------------------------------------------"
187  << "-------------------------------" << G4endl;
188  G4cout << " Position - x (cm) : " << std::setw(30) << preStep->GetPosition().x() / CLHEP::cm << std::setw(30)
189  << postStep->GetPosition().x() / CLHEP::cm << G4endl;
190  G4cout << " Position - y (cm) : " << std::setw(30) << preStep->GetPosition().y() / CLHEP::cm << std::setw(30)
191  << postStep->GetPosition().y() / CLHEP::cm << G4endl;
192  G4cout << " Position - z (cm) : " << std::setw(30) << preStep->GetPosition().z() / CLHEP::cm << std::setw(30)
193  << postStep->GetPosition().z() / CLHEP::cm << G4endl;
194  G4cout << " Global Time (ns) : " << std::setw(30) << preStep->GetGlobalTime() / CLHEP::ns << std::setw(30)
195  << postStep->GetGlobalTime() / CLHEP::ns << G4endl;
196  G4cout << " Local Time (ns) : " << std::setw(30) << preStep->GetLocalTime() / CLHEP::ns << std::setw(30)
197  << postStep->GetLocalTime() / CLHEP::ns << G4endl;
198  G4cout << " Proper Time (ns) : " << std::setw(30) << preStep->GetProperTime() / CLHEP::ns << std::setw(30)
199  << postStep->GetProperTime() / CLHEP::ns << G4endl;
200  G4cout << " Momentum Direct - x : " << std::setw(30) << preStep->GetMomentumDirection().x() << std::setw(30)
201  << postStep->GetMomentumDirection().x() << G4endl;
202  G4cout << " Momentum Direct - y : " << std::setw(30) << preStep->GetMomentumDirection().y() << std::setw(30)
203  << postStep->GetMomentumDirection().y() << G4endl;
204  G4cout << " Momentum Direct - z : " << std::setw(30) << preStep->GetMomentumDirection().z() << std::setw(30)
205  << postStep->GetMomentumDirection().z() << G4endl;
206  G4cout << " Momentum - x (GeV/c): " << std::setw(30) << preStep->GetMomentum().x() / CLHEP::GeV
207  << std::setw(30) << postStep->GetMomentum().x() / CLHEP::GeV << G4endl;
208  G4cout << " Momentum - y (GeV/c): " << std::setw(30) << preStep->GetMomentum().y() / CLHEP::GeV
209  << std::setw(30) << postStep->GetMomentum().y() / CLHEP::GeV << G4endl;
210  G4cout << " Momentum - z (GeV/c): " << std::setw(30) << preStep->GetMomentum().z() / CLHEP::GeV
211  << std::setw(30) << postStep->GetMomentum().z() / CLHEP::GeV << G4endl;
212  G4cout << " Total Energy (GeV) : " << std::setw(30) << preStep->GetTotalEnergy() / CLHEP::GeV << std::setw(30)
213  << postStep->GetTotalEnergy() / CLHEP::GeV << G4endl;
214  G4cout << " Kinetic Energy (GeV): " << std::setw(30) << preStep->GetKineticEnergy() / CLHEP::GeV
215  << std::setw(30) << postStep->GetKineticEnergy() / CLHEP::GeV << G4endl;
216  G4cout << " Velocity (mm/ns) : " << std::setw(30) << preStep->GetVelocity() << std::setw(30)
217  << postStep->GetVelocity() << G4endl;
218  G4cout << " Volume Name : " << std::setw(30) << preStep->GetPhysicalVolume()->GetName();
219  G4String volName = (postStep->GetPhysicalVolume()) ? postStep->GetPhysicalVolume()->GetName() : "OutOfWorld";
220 
221  G4cout << std::setw(30) << volName << G4endl;
222  G4cout << " Safety (mm) : " << std::setw(30) << preStep->GetSafety() << std::setw(30)
223  << postStep->GetSafety() << G4endl;
224  G4cout << " Polarization - x : " << std::setw(30) << preStep->GetPolarization().x() << std::setw(30)
225  << postStep->GetPolarization().x() << G4endl;
226  G4cout << " Polarization - y : " << std::setw(30) << preStep->GetPolarization().y() << std::setw(30)
227  << postStep->GetPolarization().y() << G4endl;
228  G4cout << " Polarization - Z : " << std::setw(30) << preStep->GetPolarization().z() << std::setw(30)
229  << postStep->GetPolarization().z() << G4endl;
230  G4cout << " Weight : " << std::setw(30) << preStep->GetWeight() << std::setw(30)
231  << postStep->GetWeight() << G4endl;
232  G4cout << " Step Status : ";
233  G4StepStatus tStepStatus = preStep->GetStepStatus();
234  if (tStepStatus == fGeomBoundary) {
235  G4cout << std::setw(30) << "Geom Limit";
236  } else if (tStepStatus == fAlongStepDoItProc) {
237  G4cout << std::setw(30) << "AlongStep Proc.";
238  } else if (tStepStatus == fPostStepDoItProc) {
239  G4cout << std::setw(30) << "PostStep Proc";
240  } else if (tStepStatus == fAtRestDoItProc) {
241  G4cout << std::setw(30) << "AtRest Proc";
242  } else if (tStepStatus == fUndefined) {
243  G4cout << std::setw(30) << "Undefined";
244  }
245 
246  tStepStatus = postStep->GetStepStatus();
247  if (tStepStatus == fGeomBoundary) {
248  G4cout << std::setw(30) << "Geom Limit";
249  } else if (tStepStatus == fAlongStepDoItProc) {
250  G4cout << std::setw(30) << "AlongStep Proc.";
251  } else if (tStepStatus == fPostStepDoItProc) {
252  G4cout << std::setw(30) << "PostStep Proc";
253  } else if (tStepStatus == fAtRestDoItProc) {
254  G4cout << std::setw(30) << "AtRest Proc";
255  } else if (tStepStatus == fUndefined) {
256  G4cout << std::setw(30) << "Undefined";
257  }
258 
259  G4cout << G4endl;
260  G4cout << " Process defined Step: ";
261  if (preStep->GetProcessDefinedStep() == nullptr) {
262  G4cout << std::setw(30) << "Undefined";
263  } else {
264  G4cout << std::setw(30) << preStep->GetProcessDefinedStep()->GetProcessName();
265  }
266  if (postStep->GetProcessDefinedStep() == nullptr) {
267  G4cout << std::setw(30) << "Undefined";
268  } else {
269  G4cout << std::setw(30) << postStep->GetProcessDefinedStep()->GetProcessName();
270  }
271  G4cout.precision(prec);
272 
273  G4cout << G4endl;
274  G4cout << " -------------------------------------------------------"
275  << "-------------------------------" << G4endl;
276  }
277 
278  // geometry information
279  if (4 <= m_verbose) {
280  const G4VTouchable* tch1 = preStep->GetTouchable();
281  const G4VTouchable* tch2 = postStep->GetTouchable();
282 
283  G4double x = postStep->GetPosition().x();
284  G4double y = postStep->GetPosition().y();
285  G4cout << "Touchable depth pre= " << tch1->GetHistoryDepth() << " post= " << tch2->GetHistoryDepth()
286  << " trans1= " << tch1->GetTranslation(tch1->GetHistoryDepth())
287  << " trans2= " << tch2->GetTranslation(tch2->GetHistoryDepth()) << " r= " << std::sqrt(x * x + y * y)
288  << G4endl;
289  const G4VPhysicalVolume* pv1 = preStep->GetPhysicalVolume();
290  const G4VPhysicalVolume* pv2 = postStep->GetPhysicalVolume();
291  const G4RotationMatrix* rotm = pv1->GetFrameRotation();
292  G4cout << "PreStepVolume: " << pv1->GetName() << G4endl;
293  G4cout << " Translation: " << pv1->GetObjectTranslation() << G4endl;
294  if (nullptr != rotm)
295  G4cout << " Rotation: " << *rotm << G4endl;
296  const G4VSolid* sv1 = pv1->GetLogicalVolume()->GetSolid();
297  sv1->StreamInfo(G4cout);
298  G4cout << G4endl;
299  if (pv2 && pv2 != pv1) {
300  G4cout << "PostStepVolume: " << pv2->GetName() << G4endl;
301  G4cout << " Translation: " << pv2->GetObjectTranslation() << G4endl;
302  rotm = pv2->GetFrameRotation();
303  if (nullptr != rotm)
304  G4cout << " Rotation: " << *rotm << G4endl;
305  const G4VSolid* sv2 = pv2->GetLogicalVolume()->GetSolid();
306  sv2->StreamInfo(G4cout);
307  }
308  G4cout << G4endl;
309 
310  if (5 <= m_verbose) {
311  G4cout << "##### Geometry Depth Analysis" << G4endl;
312  for (G4int k = 1; k < tch1->GetHistoryDepth(); ++k) {
313  const G4VPhysicalVolume* pv = tch1->GetVolume(k);
314  if (pv) {
315  const G4VSolid* sol = pv->GetLogicalVolume()->GetSolid();
316  G4cout << " Depth # " << k << " PhysVolume " << pv->GetName() << G4endl;
317  G4cout << " Translation: " << pv->GetObjectTranslation() << G4endl;
318  const G4RotationMatrix* rotm = pv->GetFrameRotation();
319  if (nullptr != rotm)
320  G4cout << " Rotation: " << *rotm << G4endl;
321  sol->StreamInfo(G4cout);
322  }
323  }
324  }
325  G4cout << G4endl;
326  }
327 
328  // verbose == 1
329  prec = G4cout.precision(4);
330  G4cout << std::setw(5) << track->GetCurrentStepNumber() << " " << std::setw(8) << track->GetPosition().x() / CLHEP::cm
331  << " " << std::setw(8) << track->GetPosition().y() / CLHEP::cm << " " << std::setw(8)
332  << track->GetPosition().z() / CLHEP::cm << " " << std::setw(9) << track->GetKineticEnergy() / CLHEP::GeV
333  << " ";
334  G4cout.precision(prec);
335  prec = G4cout.precision(3);
336  G4cout << std::setw(8) << step->GetTotalEnergyDeposit() / CLHEP::MeV << " " << std::setw(8)
337  << step->GetStepLength() / CLHEP::mm << " " << std::setw(9) << track->GetTrackLength() / CLHEP::cm << " ";
338 
339  G4bool endTracking = false;
340  if (track->GetNextVolume() != nullptr) {
341  G4cout << std::setw(30) << track->GetVolume()->GetName() << " ";
342  } else {
343  G4cout << std::setw(30) << "OutOfWorld"
344  << " ";
345  endTracking = true;
346  }
347  if (isKilled) {
348  G4cout << "isKilled";
349  endTracking = true;
350  } else if (postStep->GetProcessDefinedStep() != nullptr) {
351  G4cout << postStep->GetProcessDefinedStep()->GetProcessName();
352  }
353  G4cout << G4endl;
354  G4cout.precision(prec);
355 
356  if (1 >= m_verbose) {
357  return;
358  }
359  // verbose > 1
360  if (!endTracking && fStopAndKill != track->GetTrackStatus()) {
361  return;
362  }
363 
364  prec = G4cout.precision(4);
365  const G4TrackVector* tv = step->GetSecondary();
366  G4int nt = tv->size();
367  if (nt > 0) {
368  G4cout << " ++List of " << nt << " secondaries generated; "
369  << " Ekin > " << m_EkinThreshold << " MeV are shown:" << G4endl;
370  }
371  for (G4int i = 0; i < nt; ++i) {
372  if ((*tv)[i]->GetKineticEnergy() < m_EkinThreshold) {
373  continue;
374  }
375  G4cout << " (" << std::setw(9) << (*tv)[i]->GetPosition().x() / CLHEP::cm << " " << std::setw(9)
376  << (*tv)[i]->GetPosition().y() / CLHEP::cm << " " << std::setw(9) << (*tv)[i]->GetPosition().z() / CLHEP::cm
377  << ") cm, " << std::setw(9) << (*tv)[i]->GetKineticEnergy() / CLHEP::GeV << " GeV, " << std::setw(9)
378  << (*tv)[i]->GetGlobalTime() / CLHEP::ns << " ns, " << std::setw(18)
379  << (*tv)[i]->GetDefinition()->GetParticleName() << G4endl;
380  }
381  G4cout.precision(prec);
382 }
mps_fire.i
i
Definition: mps_fire.py:429
CMSSteppingVerbose::beginOfEvent
void beginOfEvent(const G4Event *)
Definition: CMSSteppingVerbose.cc:32
CMSSteppingVerbose::m_EventNumbers
std::vector< G4int > m_EventNumbers
Definition: CMSSteppingVerbose.h:48
ecalTB2006H4_GenSimDigiReco_cfg.G4cout
G4cout
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:285
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CMSSteppingVerbose::m_TrackNumbers
std::vector< G4int > m_TrackNumbers
Definition: CMSSteppingVerbose.h:50
CMSSteppingVerbose::CMSSteppingVerbose
CMSSteppingVerbose(G4int verb, G4double ekin, std::vector< G4int > &evtNum, std::vector< G4int > &primV, std::vector< G4int > &trNum)
Definition: CMSSteppingVerbose.cc:14
CMSSteppingVerbose::stackFilled
void stackFilled(const G4Track *, bool isKilled) const
Definition: CMSSteppingVerbose.cc:134
nt
int nt
Definition: AMPTWrapper.h:42
CMSSteppingVerbose::m_g4SteppingVerbose
G4SteppingVerbose * m_g4SteppingVerbose
Definition: CMSSteppingVerbose.h:52
mkfit::Const::sol
constexpr float sol
Definition: Config.h:13
CMSSteppingVerbose::m_PrimaryVertex
std::vector< G4int > m_PrimaryVertex
Definition: CMSSteppingVerbose.h:49
CMSSteppingVerbose::nextStep
void nextStep(const G4Step *, const G4SteppingManager *ptr, bool isKilled) const
Definition: CMSSteppingVerbose.cc:154
CMSSteppingVerbose::trackStarted
void trackStarted(const G4Track *, bool isKilled)
Definition: CMSSteppingVerbose.cc:66
step
step
Definition: StallMonitor.cc:98
HLT_2023v12_cff.track
track
Definition: HLT_2023v12_cff.py:9616
Generated for CMSSW Reference Manual by doxygen 1.8.14