Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:31 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: TrackingVerboseAction.cc
0003 // Creation: P.Arce  09/01
0004 // Modifications: porting to CMSSW by M. Stavrianakou 22/03/06
0005 // Description:
0006 ///////////////////////////////////////////////////////////////////////////////
0007 
0008 #include "SimG4Core/TrackingVerbose/interface/TrackingVerboseAction.h"
0009 
0010 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0011 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0012 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0013 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0014 #include "SimG4Core/Notification/interface/TrackInformation.h"
0015 #include "SimG4Core/Application/interface/TrackingAction.h"
0016 
0017 #include "G4Track.hh"
0018 #include "G4Event.hh"
0019 #include "G4ios.hh"
0020 #include "G4TrackingManager.hh"
0021 #include "G4EventManager.hh"
0022 #include "G4VSteppingVerbose.hh"
0023 #include "G4UnitsTable.hh"
0024 
0025 #include <algorithm>
0026 
0027 using namespace CLHEP;
0028 
0029 TrackingVerboseAction::TrackingVerboseAction(edm::ParameterSet const& p)
0030     : theTrackingManager(nullptr), fVerbose(nullptr) {
0031   fLarge = int(1E10);
0032   fDEBUG = p.getUntrackedParameter<bool>("DEBUG", false);
0033   fHighEtPhotons = p.getUntrackedParameter<bool>("CheckForHighEtPhotons", false);
0034   fG4Verbose = p.getUntrackedParameter<bool>("G4Verbose", false);
0035 
0036   //----- Set which events are verbose
0037   fTVEventMin = p.getUntrackedParameter<int>("EventMin", 0);
0038   fTVEventMax = p.getUntrackedParameter<int>("EventMax", fLarge);
0039   fTVEventStep = p.getUntrackedParameter<int>("EventStep", 1);
0040 
0041   //----- Set which tracks of those events are verbose
0042   fTVTrackMin = p.getUntrackedParameter<int>("TrackMin", 0);
0043   fTVTrackMax = p.getUntrackedParameter<int>("TrackMax", fLarge);
0044   fTVTrackStep = p.getUntrackedParameter<int>("TrackStep", 1);
0045 
0046   //----- Set the verbosity level
0047   fVerboseLevel = p.getUntrackedParameter<int>("VerboseLevel", 1);
0048   fPdgIds = p.getUntrackedParameter<std::vector<int> >("PDGids");
0049   if (fDEBUG) {
0050     G4cout << "TV: fTVTrackMin " << fTVTrackMin << " fTVTrackMax " << fTVTrackMax << " fTVTrackStep " << fTVTrackStep
0051            << " fTVEventMin " << fTVEventMin << " fTVEventMax " << fTVEventMax << " fTVEventStep " << fTVEventStep
0052            << " fVerboseLevel " << fVerboseLevel << " fG4Verbose " << fG4Verbose << " PDGIds     " << fPdgIds.size()
0053            << G4endl;
0054     for (unsigned int ii = 0; ii < fPdgIds.size(); ++ii)
0055       G4cout << "TV: PDGId[" << ii << "] = " << fPdgIds[ii] << G4endl;
0056   }
0057 
0058   //----- Set verbosity off to start
0059   fTrackingVerboseON = false;
0060   fTkVerbThisEventON = false;
0061 
0062   G4cout << " TrackingVerbose constructed " << G4endl;
0063 }
0064 
0065 TrackingVerboseAction::~TrackingVerboseAction() {}
0066 
0067 void TrackingVerboseAction::update(const BeginOfRun* run) {
0068   /*
0069   TrackingAction * ta = 
0070     dynamic_cast<TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
0071   theTrackingManager = ta->getTrackManager();
0072   */
0073   fVerbose = G4VSteppingVerbose::GetInstance();
0074   if (fDEBUG)
0075     G4cout << " TV: Get the Tracking Manager: " << theTrackingManager << " and the SteppingVerbose: " << fVerbose
0076            << G4endl;
0077 }
0078 
0079 void TrackingVerboseAction::update(const BeginOfEvent* evt) {
0080   if (evt == nullptr)
0081     return;
0082   const G4Event* anEvent = (*evt)();
0083   if (anEvent == nullptr)
0084     return;
0085 
0086   //----------- Set /tracking/verbose for this event
0087   int eventNo = anEvent->GetEventID();
0088   if (fDEBUG)
0089     G4cout << "TV: trackID: NEW EVENT " << eventNo << G4endl;
0090 
0091   fTkVerbThisEventON = false;
0092   //----- Check if event is in the selected range
0093   if (eventNo >= fTVEventMin && eventNo <= fTVEventMax) {
0094     if ((eventNo - fTVEventMin) % fTVEventStep == 0)
0095       fTkVerbThisEventON = true;
0096   }
0097 
0098   if (fDEBUG)
0099     G4cout << " TV: fTkVerbThisEventON " << fTkVerbThisEventON << " fTrackingVerboseON " << fTrackingVerboseON
0100            << " fTVEventMin " << fTVEventMin << " fTVEventMax " << fTVEventMax << G4endl;
0101   //----- check if verbosity has to be changed
0102   if ((fTkVerbThisEventON) && (!fTrackingVerboseON)) {
0103     if (fTVTrackMin == 0 && fTVTrackMax == fLarge && fTVTrackStep != 1) {
0104       setTrackingVerbose(fVerboseLevel);
0105       fTrackingVerboseON = true;
0106       if (fDEBUG)
0107         G4cout << "TV: VERBOSEet1 " << eventNo << G4endl;
0108     }
0109   } else if ((!fTkVerbThisEventON) && (fTrackingVerboseON)) {
0110     setTrackingVerbose(0);
0111     fTrackingVerboseON = false;
0112     if (fDEBUG)
0113       G4cout << "TV: VERBOSEet0 " << eventNo << G4endl;
0114   }
0115 }
0116 
0117 void TrackingVerboseAction::update(const BeginOfTrack* trk) {
0118   const G4Track* aTrack = (*trk)();
0119 
0120   //----- High ET photon printout
0121   //---------- Set /tracking/verbose
0122   //----- track is verbose only if event is verbose
0123   double tkP = aTrack->GetMomentum().mag();
0124   double tkPx = aTrack->GetMomentum().x();
0125   double tkPy = aTrack->GetMomentum().y();
0126   double tkPz = aTrack->GetMomentum().z();
0127 
0128   double tvtx = aTrack->GetVertexPosition().x();
0129   double tvty = aTrack->GetVertexPosition().y();
0130   double tvtz = aTrack->GetVertexPosition().z();
0131 
0132   double g4t_phi = atan2(tkPy, tkPx);
0133 
0134   double drpart = sqrt(tkPx * tkPx + tkPy * tkPy);
0135 
0136   double mythetapart = acos(tkPz / sqrt(drpart * drpart + tkPz * tkPz));
0137 
0138   double g4t_eta = -log(tan(mythetapart / 2.));
0139   G4int MytrackNo = aTrack->GetTrackID();
0140 
0141   if (fHighEtPhotons) {
0142     if (aTrack->GetDefinition()->GetParticleName() == "gamma" && aTrack->GetParentID() != 0) {
0143       if ((tkPx * tkPx + tkPy * tkPy + tkPz * tkPz) > 1000.0 * 1000.0 &&
0144           aTrack->GetCreatorProcess()->GetProcessName() == "LCapture") {
0145         G4cout << "MY NEW GAMMA " << G4endl;
0146         G4cout << "**********************************************************************" << G4endl;
0147         G4cout << "MY NEW TRACK ID = " << MytrackNo << "(" << aTrack->GetDefinition()->GetParticleName() << ")"
0148                << " PARENT =" << aTrack->GetParentID() << G4endl;
0149         G4cout << "Primary particle: " << aTrack->GetDynamicParticle()->GetPrimaryParticle() << G4endl;
0150         G4cout << "Process type: " << aTrack->GetCreatorProcess()->GetProcessType()
0151                << " Process name: " << aTrack->GetCreatorProcess()->GetProcessName() << G4endl;
0152         G4cout << "ToT E = " << aTrack->GetTotalEnergy() << " KineE = " << aTrack->GetKineticEnergy()
0153                << " Tot P = " << tkP << " Pt = " << sqrt(tkPx * tkPx + tkPy * tkPy) << " VTX=(" << tvtx << "," << tvty
0154                << "," << tvtz << ")" << G4endl;
0155         if (aTrack->GetKineticEnergy() > 1. * GeV && aTrack->GetCreatorProcess()->GetProcessName() != "LCapture")
0156           G4cout << " KineE > 1 GeV !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << G4endl;
0157         const G4VTouchable* touchable = aTrack->GetTouchable();
0158         if (touchable != nullptr && touchable->GetVolume() != nullptr &&
0159             touchable->GetVolume()->GetLogicalVolume() != nullptr) {
0160           G4Material* material = touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
0161           G4cout << "G4LCapture Gamma E(GeV) " << aTrack->GetTotalEnergy() / GeV << "  " << material->GetName() << " "
0162                  << touchable->GetVolume()->GetName() << G4endl;
0163           G4cout << "G4LCapture Gamma position(m): " << aTrack->GetPosition() / m << G4endl;
0164           G4cout << "G4LCapture created Gamma direction " << aTrack->GetMomentumDirection() << G4endl;
0165           G4cout << "G4LCapture gamma (eta,phi) = "
0166                  << "(" << g4t_eta << "," << g4t_phi << ")" << G4endl;
0167         }
0168         aTrack->GetUserInformation()->Print();
0169         G4cout << "**********************************************************************" << G4endl;
0170       }
0171     }
0172 
0173     if (aTrack->GetDefinition()->GetParticleName() == "gamma") {
0174       const G4VProcess* proc = aTrack->GetCreatorProcess();
0175       double Tgamma = aTrack->GetKineticEnergy();
0176       std::string ProcName;
0177       const std::string nullstr("Null_prc");
0178       if (proc)
0179         ProcName = proc->GetProcessName();
0180       else
0181         ProcName = nullstr;
0182       if (Tgamma > 2.5 * GeV) {  //&& ProcName!="Decay" && ProcName!="eBrem")
0183         std::string volumeName("_Unknown_Vol_");
0184         std::string materialName("_Unknown_Mat_");
0185         G4Material* material = nullptr;
0186         G4VPhysicalVolume* pvolume = nullptr;
0187         G4LogicalVolume* lvolume = nullptr;
0188         const G4VTouchable* touchable = aTrack->GetTouchable();
0189         if (touchable)
0190           pvolume = touchable->GetVolume();
0191         if (pvolume) {
0192           volumeName = pvolume->GetName();
0193           lvolume = pvolume->GetLogicalVolume();
0194         }
0195         if (lvolume)
0196           material = lvolume->GetMaterial();
0197         if (material)
0198           materialName = material->GetName();
0199         G4cout << "**** ALL photons > 2.5 GeV ****" << G4endl;
0200         G4cout << ProcName << "**** ALL photons: gamma E(GeV) " << aTrack->GetTotalEnergy() / GeV << "  "
0201                << materialName << " " << volumeName << G4endl;
0202         G4cout << ProcName << "**** ALL photons: gamma position(m): " << aTrack->GetPosition() / m << G4endl;
0203         G4cout << ProcName << "**** ALL photons: gamma direction " << aTrack->GetMomentumDirection() << G4endl;
0204         G4cout << "**********************************************************************" << G4endl;
0205       }
0206     }
0207   }
0208 
0209   //---------- Set /tracking/verbose
0210   //----- track is verbose only if event is verbose
0211   if (fTkVerbThisEventON) {
0212     bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
0213 
0214     //----- Set the /tracking/verbose for this track
0215     if ((trackingVerboseThisTrack) && (!fTrackingVerboseON)) {
0216       setTrackingVerbose(fVerboseLevel);
0217       fTrackingVerboseON = true;
0218       if (fDEBUG)
0219         G4cout << "TV: VERBOSEtt1 " << aTrack->GetTrackID() << G4endl;
0220       printTrackInfo(aTrack);
0221     } else if ((!trackingVerboseThisTrack) && (fTrackingVerboseON)) {
0222       setTrackingVerbose(0);
0223       fTrackingVerboseON = false;
0224       if (fDEBUG)
0225         G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID() << G4endl;
0226     }
0227   }
0228 }
0229 
0230 void TrackingVerboseAction::update(const EndOfTrack* trk) {
0231   const G4Track* aTrack = (*trk)();
0232   if (fTkVerbThisEventON) {
0233     bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
0234     if ((trackingVerboseThisTrack) && (fTrackingVerboseON) && (fTVTrackMax < fLarge || fTVTrackStep != 1)) {
0235       setTrackingVerbose(0);
0236       fTrackingVerboseON = false;
0237       if (fDEBUG)
0238         G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID() << G4endl;
0239     }
0240   }
0241 }
0242 
0243 void TrackingVerboseAction::update(const G4Step* fStep) {
0244   if ((fG4Verbose) && (fTrackingVerboseON)) {
0245     G4Track* fTrack = fStep->GetTrack();
0246     G4cout << std::setw(5) << fTrack->GetCurrentStepNumber() << " " << std::setw(8)
0247            << G4BestUnit(fTrack->GetPosition().x(), "Length") << " " << std::setw(8)
0248            << G4BestUnit(fTrack->GetPosition().y(), "Length") << " " << std::setw(8)
0249            << G4BestUnit(fTrack->GetPosition().z(), "Length") << " " << std::setw(9)
0250            << G4BestUnit(fTrack->GetKineticEnergy(), "Energy") << " " << std::setw(8)
0251            << G4BestUnit(fStep->GetTotalEnergyDeposit(), "Energy") << " " << std::setw(8)
0252            << G4BestUnit(fStep->GetStepLength(), "Length") << " " << std::setw(9)
0253            << G4BestUnit(fTrack->GetTrackLength(), "Length") << " " << std::setw(9)
0254            << G4BestUnit(fTrack->GetGlobalTime(), "Time") << " ";
0255 
0256     // Put cut comment here
0257     if (fTrack->GetNextVolume() != nullptr) {
0258       G4cout << std::setw(11) << fTrack->GetNextVolume()->GetName() << " ";
0259     } else {
0260       G4cout << std::setw(11) << "OutOfWorld"
0261              << " ";
0262     }
0263     if (fStep->GetPostStepPoint()->GetProcessDefinedStep() != nullptr) {
0264       G4cout << fStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
0265     } else {
0266       G4cout << "User Limit";
0267     }
0268     G4cout << G4endl;
0269   }
0270 }
0271 
0272 void TrackingVerboseAction::setTrackingVerbose(int verblev) {
0273   if (fDEBUG)
0274     G4cout << " setting verbose level " << verblev << G4endl;
0275   if (theTrackingManager != nullptr)
0276     theTrackingManager->SetVerboseLevel(verblev);
0277 }
0278 
0279 bool TrackingVerboseAction::checkTrackingVerbose(const G4Track* aTrack) {
0280   int trackNo = aTrack->GetTrackID();
0281   bool trackingVerboseThisTrack = false;
0282   //----- Check if track is in the selected range
0283   if (trackNo >= fTVTrackMin && trackNo <= fTVTrackMax) {
0284     if ((trackNo - fTVTrackMin) % fTVTrackStep == 0)
0285       trackingVerboseThisTrack = true;
0286   }
0287   if (trackingVerboseThisTrack && (!fPdgIds.empty())) {
0288     int pdgId = aTrack->GetDefinition()->GetPDGEncoding();
0289     if (std::count(fPdgIds.begin(), fPdgIds.end(), pdgId) == 0)
0290       trackingVerboseThisTrack = false;
0291   }
0292   return trackingVerboseThisTrack;
0293 }
0294 
0295 void TrackingVerboseAction::printTrackInfo(const G4Track* aTrack) {
0296   G4cout << G4endl << "*******************************************************"
0297          << "**************************************************" << G4endl << "* G4Track Information: "
0298          << "  Particle = " << aTrack->GetDefinition()->GetParticleName() << ","
0299          << "   Track ID = " << aTrack->GetTrackID() << ","
0300          << "   Parent ID = " << aTrack->GetParentID() << G4endl
0301          << "*******************************************************"
0302          << "**************************************************" << G4endl << G4endl;
0303   if (fVerbose)
0304     fVerbose->TrackingStarted();
0305 }
[ 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!