TrackingVerboseAction.cc

Go to the documentation of this file.

// File: TrackingVerboseAction.cc
// Creation: P.Arce  09/01
// Modifications: porting to CMSSW by M. Stavrianakou 22/03/06
// Description:

#include "SimG4Core/TrackingVerbose/interface/TrackingVerboseAction.h"

#include "SimG4Core/Notification/interface/BeginOfRun.h"
#include "SimG4Core/Notification/interface/BeginOfEvent.h"
#include "SimG4Core/Notification/interface/BeginOfTrack.h"
#include "SimG4Core/Notification/interface/EndOfTrack.h"
#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "SimG4Core/Application/interface/TrackingAction.h"

#include "G4Track.hh"
#include "G4Event.hh"
#include "G4ios.hh"
#include "G4TrackingManager.hh"
#include "G4EventManager.hh"
#include "G4VSteppingVerbose.hh"
#include "G4UnitsTable.hh"

#include<algorithm>

using namespace CLHEP;

TrackingVerboseAction::TrackingVerboseAction(edm::ParameterSet const & p) :
  theTrackingManager(nullptr), fVerbose(nullptr) {

  fLarge = int(1E10);
  fDEBUG = p.getUntrackedParameter<bool>("DEBUG",false);
  fHighEtPhotons = p.getUntrackedParameter<bool>("CheckForHighEtPhotons",false);
  fG4Verbose = p.getUntrackedParameter<bool>("G4Verbose",false);

  //----- Set which events are verbose
  fTVEventMin  = p.getUntrackedParameter<int>("EventMin",0);
  fTVEventMax  = p.getUntrackedParameter<int>("EventMax",fLarge);
  fTVEventStep = p.getUntrackedParameter<int>("EventStep",1);

  //----- Set which tracks of those events are verbose
  fTVTrackMin  = p.getUntrackedParameter<int>("TrackMin",0);
  fTVTrackMax  = p.getUntrackedParameter<int>("TrackMax",fLarge);
  fTVTrackStep = p.getUntrackedParameter<int>("TrackStep",1);

  //----- Set the verbosity level
  fVerboseLevel = p.getUntrackedParameter<int>("VerboseLevel",1);
  fPdgIds       = p.getUntrackedParameter<std::vector<int> >("PDGids");
  if (fDEBUG) {
    G4cout << "TV: fTVTrackMin " << fTVTrackMin 
       << " fTVTrackMax "    <<  fTVTrackMax 
       <<  " fTVTrackStep "  << fTVTrackStep  
       << " fTVEventMin "    << fTVEventMin 
       << " fTVEventMax "    << fTVEventMax   
       << " fTVEventStep "   << fTVEventStep 
       << " fVerboseLevel "  << fVerboseLevel 
       << " fG4Verbose "     << fG4Verbose 
       << " PDGIds     "     << fPdgIds.size() << G4endl;
    for (unsigned int ii=0; ii<fPdgIds.size(); ++ii) 
      G4cout << "TV: PDGId[" << ii << "] = " << fPdgIds[ii] << G4endl;
  }
  
  //----- Set verbosity off to start
  fTrackingVerboseON = false;
  fTkVerbThisEventON = false;
  
  G4cout << " TrackingVerbose constructed " << G4endl;
}

TrackingVerboseAction::~TrackingVerboseAction() {}

void TrackingVerboseAction::update(const BeginOfRun * run) {
  TrackingAction * ta = 
    dynamic_cast<TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
  theTrackingManager = ta->getTrackManager();
  fVerbose = G4VSteppingVerbose::GetInstance();
  if (fDEBUG)
    G4cout << " TV: Get the Tracking Manager: " << theTrackingManager
       << " and the SteppingVerbose: " << fVerbose << G4endl;
}

void TrackingVerboseAction::update(const BeginOfEvent * evt) {
  if (evt==nullptr) return;
  const G4Event * anEvent = (*evt)();
  if (anEvent==nullptr) return;

  //----------- Set /tracking/verbose for this event 
  int eventNo = anEvent->GetEventID();
  if (fDEBUG) G4cout << "TV: trackID: NEW EVENT " << eventNo << G4endl;

  fTkVerbThisEventON = false;
  //----- Check if event is in the selected range
  if (eventNo >= fTVEventMin && eventNo <= fTVEventMax) {
    if ((eventNo-fTVEventMin) % fTVEventStep == 0) fTkVerbThisEventON = true;
  }

  if (fDEBUG)
    G4cout << " TV: fTkVerbThisEventON " <<  fTkVerbThisEventON 
       << " fTrackingVerboseON " << fTrackingVerboseON 
       << " fTVEventMin " << fTVEventMin << " fTVEventMax " << fTVEventMax << G4endl;
  //----- check if verbosity has to be changed
  if ((fTkVerbThisEventON) && (!fTrackingVerboseON)) {
    if (fTVTrackMin == 0 && fTVTrackMax == fLarge && fTVTrackStep != 1) {
      setTrackingVerbose(fVerboseLevel);
      fTrackingVerboseON = true;
      if (fDEBUG) G4cout << "TV: VERBOSEet1 " << eventNo << G4endl;
    }
  } else if ((!fTkVerbThisEventON) && (fTrackingVerboseON) ) {
    setTrackingVerbose(0);
    fTrackingVerboseON = false;
    if (fDEBUG) G4cout << "TV: VERBOSEet0 " << eventNo << G4endl;
  }

}

void TrackingVerboseAction::update(const BeginOfTrack * trk) {
  const G4Track * aTrack = (*trk)();

  //----- High ET photon printout
  //---------- Set /tracking/verbose
  //----- track is verbose only if event is verbose
  double tkP = aTrack->GetMomentum().mag();
  double tkPx = aTrack->GetMomentum().x();
  double tkPy = aTrack->GetMomentum().y();
  double tkPz = aTrack->GetMomentum().z();

  double tvtx = aTrack->GetVertexPosition().x();
  double tvty = aTrack->GetVertexPosition().y();
  double tvtz = aTrack->GetVertexPosition().z();

  double g4t_phi=atan2(tkPy,tkPx);

  double drpart=sqrt(tkPx*tkPx + tkPy*tkPy);

  double mythetapart=acos(tkPz/sqrt(drpart*drpart+tkPz*tkPz));

  double g4t_eta=-log(tan(mythetapart/2.));
  G4int MytrackNo = aTrack->GetTrackID();
    
  if (fHighEtPhotons) {
    if (aTrack->GetDefinition()->GetParticleName() == "gamma" && aTrack->GetParentID() !=0) {
      if((tkPx*tkPx + tkPy*tkPy + tkPz*tkPz)>1000.0*1000.0 &&
     aTrack->GetCreatorProcess()->GetProcessName() == "LCapture") {
    G4cout << "MY NEW GAMMA " << G4endl;
    G4cout << "**********************************************************************"  << G4endl;
    G4cout << "MY NEW TRACK ID = " << MytrackNo << "("
           << aTrack->GetDefinition()->GetParticleName()
           <<")"<< " PARENT ="<< aTrack->GetParentID() << G4endl;
    G4cout << "Primary particle: " 
           << aTrack->GetDynamicParticle()->GetPrimaryParticle() << G4endl;
    G4cout << "Process type: " << aTrack->GetCreatorProcess()->GetProcessType()
           << " Process name: " 
           << aTrack->GetCreatorProcess()->GetProcessName() << G4endl;
    G4cout << "ToT E = " << aTrack->GetTotalEnergy() 
           << " KineE = " << aTrack->GetKineticEnergy()
           << " Tot P = " << tkP << " Pt = " << sqrt(tkPx*tkPx + tkPy*tkPy) 
           << " VTX=(" << tvtx << "," << tvty << "," << tvtz << ")" << G4endl;
    if (aTrack->GetKineticEnergy() > 1.*GeV 
        && aTrack->GetCreatorProcess()->GetProcessName() != "LCapture")
      G4cout << " KineE > 1 GeV !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << G4endl;
    const G4VTouchable* touchable=aTrack->GetTouchable();
    if (touchable!=nullptr && touchable->GetVolume()!=nullptr &&
        touchable->GetVolume()->GetLogicalVolume()!=nullptr) {
      G4Material* material=touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
      G4cout << "G4LCapture Gamma E(GeV) " 
         << aTrack->GetTotalEnergy()/GeV << "  "
         << material->GetName()<< " " 
         << touchable->GetVolume()->GetName() << G4endl;
      G4cout << "G4LCapture Gamma position(m): " 
         << aTrack->GetPosition()/m << G4endl;
      G4cout << "G4LCapture created Gamma direction " 
         << aTrack->GetMomentumDirection() << G4endl;
      G4cout << "G4LCapture gamma (eta,phi) = " 
         << "(" << g4t_eta << "," << g4t_phi << ")" << G4endl;
    }
    aTrack->GetUserInformation()->Print();
    G4cout << "**********************************************************************"  << G4endl;
      }
    }

    if (aTrack->GetDefinition()->GetParticleName() == "gamma") {
      const G4VProcess * proc = aTrack->GetCreatorProcess();
      double Tgamma = aTrack->GetKineticEnergy();
      std::string ProcName;
      const  std::string nullstr ("Null_prc");
      if (proc) ProcName = proc->GetProcessName();
      else      ProcName = nullstr;
      if (Tgamma > 2.5*GeV ) { //&& ProcName!="Decay" && ProcName!="eBrem")
    std::string volumeName("_Unknown_Vol_");
    std::string materialName("_Unknown_Mat_");
    G4Material * material = nullptr;
    G4VPhysicalVolume * pvolume = nullptr;
    G4LogicalVolume * lvolume = nullptr;
    const G4VTouchable * touchable = aTrack->GetTouchable();
    if (touchable) pvolume = touchable->GetVolume();
    if (pvolume) {
      volumeName = pvolume->GetName();
      lvolume = pvolume->GetLogicalVolume();
    }
    if (lvolume) material = lvolume->GetMaterial();
    if (material) materialName = material->GetName();
    G4cout << "**** ALL photons > 2.5 GeV ****" << G4endl;
    G4cout << ProcName << "**** ALL photons: gamma E(GeV) "
           << aTrack->GetTotalEnergy()/GeV << "  "
           <<  materialName << " " << volumeName << G4endl;
    G4cout << ProcName << "**** ALL photons: gamma position(m): " 
           << aTrack->GetPosition()/m << G4endl;
    G4cout << ProcName << "**** ALL photons: gamma direction " 
           << aTrack->GetMomentumDirection() << G4endl;
    G4cout << "**********************************************************************"  << G4endl;
      }
    }                                               
  }
    
  //---------- Set /tracking/verbose
  //----- track is verbose only if event is verbose
  if (fTkVerbThisEventON) {
    bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);

    //----- Set the /tracking/verbose for this track 
    if ((trackingVerboseThisTrack) && (!fTrackingVerboseON) ) {
      setTrackingVerbose(fVerboseLevel);
      fTrackingVerboseON = true;
      if (fDEBUG) G4cout << "TV: VERBOSEtt1 " << aTrack->GetTrackID()
             << G4endl;
      printTrackInfo(aTrack);
    } else if ((!trackingVerboseThisTrack) && ( fTrackingVerboseON )) {
      setTrackingVerbose(0);
      fTrackingVerboseON = false;
      if (fDEBUG) G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID()
             << G4endl;
    }
  }
}

void TrackingVerboseAction::update(const EndOfTrack * trk) {
  const G4Track * aTrack = (*trk)();
  if (fTkVerbThisEventON) {
    bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
    if ((trackingVerboseThisTrack) && (fTrackingVerboseON ) &&
    (fTVTrackMax < fLarge || fTVTrackStep != 1)) {
      setTrackingVerbose(0);
      fTrackingVerboseON = false;
      if (fDEBUG) G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID()
             << G4endl;
    }
  }
}

void TrackingVerboseAction::update(const G4Step* fStep) {

  if ((fG4Verbose) && (fTrackingVerboseON)) {
    G4Track* fTrack = fStep->GetTrack();
    G4cout << std::setw( 5) << fTrack->GetCurrentStepNumber() << " "
       << std::setw( 8) << G4BestUnit(fTrack->GetPosition().x() , "Length") << " "
       << std::setw( 8) << G4BestUnit(fTrack->GetPosition().y() , "Length") << " "
       << std::setw( 8) << G4BestUnit(fTrack->GetPosition().z() , "Length") << " "
       << std::setw( 9) << G4BestUnit(fTrack->GetKineticEnergy() , "Energy") << " "
       << std::setw( 8) << G4BestUnit(fStep->GetTotalEnergyDeposit(), "Energy") << " "
       << std::setw( 8) << G4BestUnit(fStep->GetStepLength() , "Length") << " "
       << std::setw( 9) << G4BestUnit(fTrack->GetTrackLength() , "Length") << " "
       << std::setw( 9) << G4BestUnit(fTrack->GetGlobalTime(), "Time") << " ";

    // Put cut comment here
    if( fTrack->GetNextVolume() != nullptr ) {
      G4cout << std::setw(11) << fTrack->GetNextVolume()->GetName() << " ";
    } else {
      G4cout << std::setw(11) << "OutOfWorld" << " ";
    }
    if(fStep->GetPostStepPoint()->GetProcessDefinedStep() != nullptr){
      G4cout << fStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
    } else {
      G4cout << "User Limit";
    }
    G4cout << G4endl;
  }
}

void TrackingVerboseAction::setTrackingVerbose(int verblev) {
  if (fDEBUG) G4cout << " setting verbose level " << verblev << G4endl;
  if (theTrackingManager!=nullptr) theTrackingManager->SetVerboseLevel(verblev);
}
 
bool TrackingVerboseAction::checkTrackingVerbose(const G4Track* aTrack) {
  int trackNo = aTrack->GetTrackID();    
  bool trackingVerboseThisTrack = false;
  //----- Check if track is in the selected range
  if (trackNo >= fTVTrackMin && trackNo <= fTVTrackMax) {
    if ((trackNo-fTVTrackMin) % fTVTrackStep == 0) trackingVerboseThisTrack = true;
  }
  if (trackingVerboseThisTrack && (!fPdgIds.empty())) {
    int pdgId = aTrack->GetDefinition()->GetPDGEncoding();
    if (std::count(fPdgIds.begin(),fPdgIds.end(),pdgId) == 0) trackingVerboseThisTrack = false;
  }
  return trackingVerboseThisTrack;
}

void TrackingVerboseAction::printTrackInfo(const G4Track* aTrack) {
  G4cout << G4endl
     << "*******************************************************"
     << "**************************************************" << G4endl
     << "* G4Track Information: "
     << "  Particle = " << aTrack->GetDefinition()->GetParticleName()
     << ","
     << "   Track ID = " << aTrack->GetTrackID()
     << ","
     << "   Parent ID = " << aTrack->GetParentID() << G4endl
     << "*******************************************************"
     << "**************************************************"
     << G4endl << G4endl;
  if (fVerbose) fVerbose->TrackingStarted();
}