Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:26 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!