Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:31:22

0001 
0002 #include "SimG4Core/Notification/interface/CMSSteppingVerbose.h"
0003 #include "G4Event.hh"
0004 #include "G4EventManager.hh"
0005 #include "G4Track.hh"
0006 #include "G4TrackStatus.hh"
0007 #include "G4Step.hh"
0008 #include "G4SteppingManager.hh"
0009 #include "G4SteppingVerbose.hh"
0010 #include "G4ParticleDefinition.hh"
0011 #include "G4VProcess.hh"
0012 #include "G4SystemOfUnits.hh"
0013 
0014 CMSSteppingVerbose::CMSSteppingVerbose(
0015     G4int verb, G4double ekin, std::vector<G4int>& evtNum, std::vector<G4int>& primV, std::vector<G4int>& trNum)
0016     : m_PrintEvent(false),
0017       m_PrintTrack(false),
0018       m_smInitialized(false),
0019       m_verbose(verb),
0020       m_EventNumbers(evtNum),
0021       m_PrimaryVertex(primV),
0022       m_TrackNumbers(trNum),
0023       m_EkinThreshold(ekin) {
0024   m_nEvents = m_EventNumbers.size();
0025   m_nVertex = m_PrimaryVertex.size();
0026   m_nTracks = m_TrackNumbers.size();
0027   m_g4SteppingVerbose = new G4SteppingVerbose();
0028   G4VSteppingVerbose::SetInstance(m_g4SteppingVerbose);
0029   m_g4SteppingVerbose->SetSilent(1);
0030 }
0031 
0032 CMSSteppingVerbose::~CMSSteppingVerbose() {}
0033 
0034 void CMSSteppingVerbose::BeginOfEvent(const G4Event* evt) {
0035   m_PrintEvent = false;
0036   if (0 >= m_verbose) {
0037     return;
0038   }
0039   if (m_nEvents == 0) {
0040     m_PrintEvent = true;
0041   } else {
0042     for (G4int i = 0; i < m_nEvents; ++i) {
0043       // check event number
0044       if (evt->GetEventID() == m_EventNumbers[i]) {
0045         // check number of vertex
0046         if (m_nVertex == m_nEvents && evt->GetNumberOfPrimaryVertex() != m_PrimaryVertex[i]) {
0047           continue;
0048         }
0049         m_PrintEvent = true;
0050         break;
0051       }
0052     }
0053   }
0054   if (!m_PrintEvent) {
0055     return;
0056   }
0057   G4cout << "========== Event #" << evt->GetEventID() << "   " << evt->GetNumberOfPrimaryVertex()
0058          << " primary vertexes ======" << G4endl;
0059   G4cout << G4endl;
0060 }
0061 
0062 void CMSSteppingVerbose::TrackStarted(const G4Track* track, bool isKilled) {
0063   m_PrintTrack = false;
0064   if (!m_PrintEvent) {
0065     return;
0066   }
0067 
0068   if (!m_smInitialized) {
0069     G4SteppingManager* stepman = G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
0070     m_g4SteppingVerbose->SetManager(stepman);
0071     stepman->SetVerboseLevel(m_verbose);
0072     m_smInitialized = true;
0073   }
0074 
0075   if (m_nTracks == 0) {
0076     if (track->GetKineticEnergy() >= m_EkinThreshold) {
0077       m_PrintTrack = true;
0078     }
0079 
0080   } else {
0081     for (G4int i = 0; i < m_nTracks; ++i) {
0082       if (track->GetTrackID() == m_TrackNumbers[i]) {
0083         m_PrintTrack = true;
0084         break;
0085       }
0086     }
0087   }
0088   if (!m_PrintTrack) {
0089     return;
0090   }
0091 
0092   G4cout << "*********************************************************************************************************"
0093          << G4endl;
0094   const G4ParticleDefinition* pd = track->GetDefinition();
0095   G4cout << "* G4Track Information:   Particle = ";
0096   if (pd) {
0097     G4cout << pd->GetParticleName();
0098   }
0099   G4cout << ",   Track ID = " << track->GetTrackID() << ",   Parent ID = " << track->GetParentID() << G4endl;
0100   G4cout << "*********************************************************************************************************"
0101          << G4endl;
0102 
0103   G4cout << std::setw(5) << "Step#"
0104          << " " << std::setw(8) << "X(cm)"
0105          << " " << std::setw(8) << "Y(cm)"
0106          << " " << std::setw(8) << "Z(cm)"
0107          << " " << std::setw(9) << "KinE(GeV)"
0108          << " " << std::setw(8) << "dE(MeV)"
0109          << " " << std::setw(8) << "Step(mm)"
0110          << " " << std::setw(9) << "TrackL(cm)"
0111          << " " << std::setw(30) << "PhysVolume"
0112          << " " << std::setw(8) << "ProcName" << G4endl;
0113 
0114   G4int prec = G4cout.precision(4);
0115 
0116   G4cout << std::setw(5) << track->GetCurrentStepNumber() << " " << std::setw(8) << track->GetPosition().x() / CLHEP::cm
0117          << " " << std::setw(8) << track->GetPosition().y() / CLHEP::cm << " " << std::setw(8)
0118          << track->GetPosition().z() / CLHEP::cm << " " << std::setw(9) << track->GetKineticEnergy() / CLHEP::GeV << " "
0119          << std::setw(8) << "         " << std::setw(8) << "         " << std::setw(9) << "          ";
0120   if (track->GetVolume() != nullptr) {
0121     G4cout << std::setw(30) << track->GetVolume()->GetName() << " ";
0122   }
0123   if (isKilled) {
0124     G4cout << "isKilled";
0125   }
0126   G4cout << G4endl;
0127   G4cout.precision(prec);
0128 }
0129 
0130 void CMSSteppingVerbose::TrackEnded(const G4Track* track) const {
0131   if (!m_PrintTrack || 1 == m_verbose) {
0132     return;
0133   }
0134 }
0135 
0136 void CMSSteppingVerbose::StackFilled(const G4Track* track, bool isKilled) const {
0137   if (2 >= m_verbose || !m_PrintTrack || track->GetKineticEnergy() < m_EkinThreshold) {
0138     return;
0139   }
0140   G4int prec = G4cout.precision(4);
0141 
0142   G4cout << std::setw(10) << track->GetTrackID() << " " << std::setw(8) << track->GetPosition().x() / CLHEP::cm << " "
0143          << std::setw(8) << track->GetPosition().y() / CLHEP::cm << " " << std::setw(8)
0144          << track->GetPosition().z() / CLHEP::cm << " " << std::setw(9) << track->GetKineticEnergy() / CLHEP::GeV
0145          << " ";
0146   if (track->GetVolume() != nullptr) {
0147     G4cout << std::setw(24) << track->GetVolume()->GetName() << " ";
0148   }
0149   if (isKilled) {
0150     G4cout << "isKilled";
0151   }
0152   G4cout << G4endl;
0153   G4cout.precision(prec);
0154 }
0155 
0156 void CMSSteppingVerbose::NextStep(const G4Step* step, const G4SteppingManager* sManager, bool isKilled) const {
0157   if (!m_PrintTrack) {
0158     return;
0159   }
0160 
0161   G4int prec;
0162   const G4Track* track = step->GetTrack();
0163   const G4StepPoint* preStep = step->GetPreStepPoint();
0164   const G4StepPoint* postStep = step->GetPostStepPoint();
0165 
0166   if (3 <= m_verbose) {
0167     m_g4SteppingVerbose->SetSilent(0);
0168     m_g4SteppingVerbose->DPSLStarted();
0169     m_g4SteppingVerbose->DPSLAlongStep();
0170     m_g4SteppingVerbose->DPSLPostStep();
0171     if (4 <= m_verbose) {
0172       m_g4SteppingVerbose->AlongStepDoItAllDone();
0173       m_g4SteppingVerbose->PostStepDoItAllDone();
0174     }
0175     m_g4SteppingVerbose->SetSilent(1);
0176 
0177     prec = G4cout.precision(16);
0178 
0179     G4cout << G4endl;
0180     G4cout << "    ++G4Step Information " << G4endl;
0181     G4cout << "      Address of G4Track    : " << track << G4endl;
0182     G4cout << "      Step Length (mm)      : " << track->GetStepLength() << G4endl;
0183     G4cout << "      Energy Deposit (MeV)  : " << step->GetTotalEnergyDeposit() << G4endl;
0184 
0185     G4cout << "   -------------------------------------------------------"
0186            << "-------------------------------" << G4endl;
0187     G4cout << "  StepPoint Information  " << std::setw(30) << "PreStep" << std::setw(30) << "PostStep" << G4endl;
0188     G4cout << "   -------------------------------------------------------"
0189            << "-------------------------------" << G4endl;
0190     G4cout << "      Position - x (cm)   : " << std::setw(30) << preStep->GetPosition().x() / CLHEP::cm << std::setw(30)
0191            << postStep->GetPosition().x() / CLHEP::cm << G4endl;
0192     G4cout << "      Position - y (cm)   : " << std::setw(30) << preStep->GetPosition().y() / CLHEP::cm << std::setw(30)
0193            << postStep->GetPosition().y() / CLHEP::cm << G4endl;
0194     G4cout << "      Position - z (cm)   : " << std::setw(30) << preStep->GetPosition().z() / CLHEP::cm << std::setw(30)
0195            << postStep->GetPosition().z() / CLHEP::cm << G4endl;
0196     G4cout << "      Global Time (ns)    : " << std::setw(30) << preStep->GetGlobalTime() / CLHEP::ns << std::setw(30)
0197            << postStep->GetGlobalTime() / CLHEP::ns << G4endl;
0198     G4cout << "      Local Time (ns)     : " << std::setw(30) << preStep->GetLocalTime() / CLHEP::ns << std::setw(30)
0199            << postStep->GetLocalTime() / CLHEP::ns << G4endl;
0200     G4cout << "      Proper Time (ns)    : " << std::setw(30) << preStep->GetProperTime() / CLHEP::ns << std::setw(30)
0201            << postStep->GetProperTime() / CLHEP::ns << G4endl;
0202     G4cout << "      Momentum Direct - x : " << std::setw(30) << preStep->GetMomentumDirection().x() << std::setw(30)
0203            << postStep->GetMomentumDirection().x() << G4endl;
0204     G4cout << "      Momentum Direct - y : " << std::setw(30) << preStep->GetMomentumDirection().y() << std::setw(30)
0205            << postStep->GetMomentumDirection().y() << G4endl;
0206     G4cout << "      Momentum Direct - z : " << std::setw(30) << preStep->GetMomentumDirection().z() << std::setw(30)
0207            << postStep->GetMomentumDirection().z() << G4endl;
0208     G4cout << "      Momentum - x (GeV/c): " << std::setw(30) << preStep->GetMomentum().x() / CLHEP::GeV
0209            << std::setw(30) << postStep->GetMomentum().x() / CLHEP::GeV << G4endl;
0210     G4cout << "      Momentum - y (GeV/c): " << std::setw(30) << preStep->GetMomentum().y() / CLHEP::GeV
0211            << std::setw(30) << postStep->GetMomentum().y() / CLHEP::GeV << G4endl;
0212     G4cout << "      Momentum - z (GeV/c): " << std::setw(30) << preStep->GetMomentum().z() / CLHEP::GeV
0213            << std::setw(30) << postStep->GetMomentum().z() / CLHEP::GeV << G4endl;
0214     G4cout << "      Total Energy (GeV)  : " << std::setw(30) << preStep->GetTotalEnergy() / CLHEP::GeV << std::setw(30)
0215            << postStep->GetTotalEnergy() / CLHEP::GeV << G4endl;
0216     G4cout << "      Kinetic Energy (GeV): " << std::setw(30) << preStep->GetKineticEnergy() / CLHEP::GeV
0217            << std::setw(30) << postStep->GetKineticEnergy() / CLHEP::GeV << G4endl;
0218     G4cout << "      Velocity (mm/ns)    : " << std::setw(30) << preStep->GetVelocity() << std::setw(30)
0219            << postStep->GetVelocity() << G4endl;
0220     G4cout << "      Volume Name         : " << std::setw(30) << preStep->GetPhysicalVolume()->GetName();
0221     G4String volName = (postStep->GetPhysicalVolume()) ? postStep->GetPhysicalVolume()->GetName() : "OutOfWorld";
0222 
0223     G4cout << std::setw(30) << volName << G4endl;
0224     G4cout << "      Safety (mm)         : " << std::setw(30) << preStep->GetSafety() << std::setw(30)
0225            << postStep->GetSafety() << G4endl;
0226     G4cout << "      Polarization - x    : " << std::setw(30) << preStep->GetPolarization().x() << std::setw(30)
0227            << postStep->GetPolarization().x() << G4endl;
0228     G4cout << "      Polarization - y    : " << std::setw(30) << preStep->GetPolarization().y() << std::setw(30)
0229            << postStep->GetPolarization().y() << G4endl;
0230     G4cout << "      Polarization - Z    : " << std::setw(30) << preStep->GetPolarization().z() << std::setw(30)
0231            << postStep->GetPolarization().z() << G4endl;
0232     G4cout << "      Weight              : " << std::setw(30) << preStep->GetWeight() << std::setw(30)
0233            << postStep->GetWeight() << G4endl;
0234     G4cout << "      Step Status         : ";
0235     G4StepStatus tStepStatus = preStep->GetStepStatus();
0236     if (tStepStatus == fGeomBoundary) {
0237       G4cout << std::setw(30) << "Geom Limit";
0238     } else if (tStepStatus == fAlongStepDoItProc) {
0239       G4cout << std::setw(30) << "AlongStep Proc.";
0240     } else if (tStepStatus == fPostStepDoItProc) {
0241       G4cout << std::setw(30) << "PostStep Proc";
0242     } else if (tStepStatus == fAtRestDoItProc) {
0243       G4cout << std::setw(30) << "AtRest Proc";
0244     } else if (tStepStatus == fUndefined) {
0245       G4cout << std::setw(30) << "Undefined";
0246     }
0247 
0248     tStepStatus = postStep->GetStepStatus();
0249     if (tStepStatus == fGeomBoundary) {
0250       G4cout << std::setw(30) << "Geom Limit";
0251     } else if (tStepStatus == fAlongStepDoItProc) {
0252       G4cout << std::setw(30) << "AlongStep Proc.";
0253     } else if (tStepStatus == fPostStepDoItProc) {
0254       G4cout << std::setw(30) << "PostStep Proc";
0255     } else if (tStepStatus == fAtRestDoItProc) {
0256       G4cout << std::setw(30) << "AtRest Proc";
0257     } else if (tStepStatus == fUndefined) {
0258       G4cout << std::setw(30) << "Undefined";
0259     }
0260 
0261     G4cout << G4endl;
0262     G4cout << "      Process defined Step: ";
0263     if (preStep->GetProcessDefinedStep() == nullptr) {
0264       G4cout << std::setw(30) << "Undefined";
0265     } else {
0266       G4cout << std::setw(30) << preStep->GetProcessDefinedStep()->GetProcessName();
0267     }
0268     if (postStep->GetProcessDefinedStep() == nullptr) {
0269       G4cout << std::setw(30) << "Undefined";
0270     } else {
0271       G4cout << std::setw(30) << postStep->GetProcessDefinedStep()->GetProcessName();
0272     }
0273     G4cout.precision(prec);
0274 
0275     G4cout << G4endl;
0276     G4cout << "   -------------------------------------------------------"
0277            << "-------------------------------" << G4endl;
0278   }
0279 
0280   // geometry information
0281   if (4 <= m_verbose) {
0282     const G4VTouchable* tch1 = preStep->GetTouchable();
0283     const G4VTouchable* tch2 = postStep->GetTouchable();
0284 
0285     G4double x = postStep->GetPosition().x();
0286     G4double y = postStep->GetPosition().y();
0287     G4cout << "Touchable depth pre= " << tch1->GetHistoryDepth() << " post= " << tch2->GetHistoryDepth()
0288            << " trans1= " << tch1->GetTranslation(tch1->GetHistoryDepth())
0289            << " trans2= " << tch2->GetTranslation(tch2->GetHistoryDepth()) << " r= " << std::sqrt(x * x + y * y)
0290            << G4endl;
0291     const G4VPhysicalVolume* pv1 = preStep->GetPhysicalVolume();
0292     const G4VPhysicalVolume* pv2 = postStep->GetPhysicalVolume();
0293     const G4RotationMatrix* rotm = pv1->GetFrameRotation();
0294     G4cout << "PreStepVolume: " << pv1->GetName() << G4endl;
0295     G4cout << "       Translation: " << pv1->GetObjectTranslation() << G4endl;
0296     if (nullptr != rotm)
0297       G4cout << "       Rotation:    " << *rotm << G4endl;
0298     const G4VSolid* sv1 = pv1->GetLogicalVolume()->GetSolid();
0299     sv1->StreamInfo(G4cout);
0300     G4cout << G4endl;
0301     if (pv2 && pv2 != pv1) {
0302       G4cout << "PostStepVolume: " << pv2->GetName() << G4endl;
0303       G4cout << "       Translation: " << pv2->GetObjectTranslation() << G4endl;
0304       rotm = pv2->GetFrameRotation();
0305       if (nullptr != rotm)
0306         G4cout << "       Rotation:    " << *rotm << G4endl;
0307       const G4VSolid* sv2 = pv2->GetLogicalVolume()->GetSolid();
0308       sv2->StreamInfo(G4cout);
0309     }
0310     G4cout << G4endl;
0311 
0312     if (5 <= m_verbose) {
0313       G4cout << "##### Geometry Depth Analysis" << G4endl;
0314       for (G4int k = 1; k < tch1->GetHistoryDepth(); ++k) {
0315         const G4VPhysicalVolume* pv = tch1->GetVolume(k);
0316         if (pv) {
0317           const G4VSolid* sol = pv->GetLogicalVolume()->GetSolid();
0318           G4cout << " Depth # " << k << "  PhysVolume " << pv->GetName() << G4endl;
0319           G4cout << "       Translation: " << pv->GetObjectTranslation() << G4endl;
0320           const G4RotationMatrix* rotm = pv->GetFrameRotation();
0321           if (nullptr != rotm)
0322             G4cout << "       Rotation:    " << *rotm << G4endl;
0323           sol->StreamInfo(G4cout);
0324         }
0325       }
0326     }
0327     G4cout << G4endl;
0328   }
0329 
0330   // verbose == 1
0331   prec = G4cout.precision(4);
0332   G4cout << std::setw(5) << track->GetCurrentStepNumber() << " " << std::setw(8) << track->GetPosition().x() / CLHEP::cm
0333          << " " << std::setw(8) << track->GetPosition().y() / CLHEP::cm << " " << std::setw(8)
0334          << track->GetPosition().z() / CLHEP::cm << " " << std::setw(9) << track->GetKineticEnergy() / CLHEP::GeV
0335          << " ";
0336   G4cout.precision(prec);
0337   prec = G4cout.precision(3);
0338   G4cout << std::setw(8) << step->GetTotalEnergyDeposit() / CLHEP::MeV << " " << std::setw(8)
0339          << step->GetStepLength() / CLHEP::mm << " " << std::setw(9) << track->GetTrackLength() / CLHEP::cm << " ";
0340 
0341   G4bool endTracking = false;
0342   if (track->GetNextVolume() != nullptr) {
0343     G4cout << std::setw(30) << track->GetVolume()->GetName() << " ";
0344   } else {
0345     G4cout << std::setw(30) << "OutOfWorld"
0346            << " ";
0347     endTracking = true;
0348   }
0349   if (isKilled) {
0350     G4cout << "isKilled";
0351     endTracking = true;
0352   } else if (postStep->GetProcessDefinedStep() != nullptr) {
0353     G4cout << postStep->GetProcessDefinedStep()->GetProcessName();
0354   }
0355   G4cout << G4endl;
0356   G4cout.precision(prec);
0357 
0358   if (1 >= m_verbose) {
0359     return;
0360   }
0361   // verbose > 1
0362   if (!endTracking && fStopAndKill != track->GetTrackStatus()) {
0363     return;
0364   }
0365 
0366   prec = G4cout.precision(4);
0367   const G4TrackVector* tv = step->GetSecondary();
0368   G4int nt = tv->size();
0369   if (nt > 0) {
0370     G4cout << "    ++List of " << nt << " secondaries generated; "
0371            << " Ekin > " << m_EkinThreshold << " MeV are shown:" << G4endl;
0372   }
0373   for (G4int i = 0; i < nt; ++i) {
0374     if ((*tv)[i]->GetKineticEnergy() < m_EkinThreshold) {
0375       continue;
0376     }
0377     G4cout << "      (" << std::setw(9) << (*tv)[i]->GetPosition().x() / CLHEP::cm << " " << std::setw(9)
0378            << (*tv)[i]->GetPosition().y() / CLHEP::cm << " " << std::setw(9) << (*tv)[i]->GetPosition().z() / CLHEP::cm
0379            << ") cm, " << std::setw(9) << (*tv)[i]->GetKineticEnergy() / CLHEP::GeV << " GeV, " << std::setw(9)
0380            << (*tv)[i]->GetGlobalTime() / CLHEP::ns << " ns, " << std::setw(18)
0381            << (*tv)[i]->GetDefinition()->GetParticleName() << G4endl;
0382   }
0383   G4cout.precision(prec);
0384 }