CMS 3D CMS Logo

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