Reference Manual

Go to the documentation of this file.
00001 
00002 // File: TrackingVerboseAction.cc
00003 // Creation: P.Arce  09/01
00004 // Modifications: porting to CMSSW by M. Stavrianakou 22/03/06
00005 // Description:
00007 
00008 #include "SimG4Core/TrackingVerbose/interface/TrackingVerboseAction.h"
00009 
00010 #include "SimG4Core/Notification/interface/BeginOfRun.h"
00011 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
00012 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
00013 #include "SimG4Core/Notification/interface/EndOfTrack.h"
00014 #include "SimG4Core/Notification/interface/TrackInformation.h"
00015 #include "SimG4Core/Application/interface/TrackingAction.h"
00016 
00017 #include "G4Track.hh"
00018 #include "G4Event.hh"
00019 #include "G4ios.hh"
00020 #include "G4TrackingManager.hh"
00021 #include "G4EventManager.hh"
00022 #include "G4VSteppingVerbose.hh"
00023 #include "G4UnitsTable.hh"
00024 
00025 #include<algorithm>
00026 
00027 TrackingVerboseAction::TrackingVerboseAction(edm::ParameterSet const & p) :
00028   theTrackingManager(0), fVerbose(0) {
00029 
00030   fLarge = int(1E10);
00031   fDEBUG = p.getUntrackedParameter<bool>("DEBUG",false);
00032   fHighEtPhotons = p.getUntrackedParameter<bool>("CheckForHighEtPhotons",false);
00033   fG4Verbose = p.getUntrackedParameter<bool>("G4Verbose",false);
00034 
00035   //----- Set which events are verbose
00036   fTVEventMin  = p.getUntrackedParameter<int>("EventMin",0);
00037   fTVEventMax  = p.getUntrackedParameter<int>("EventMax",fLarge);
00038   fTVEventStep = p.getUntrackedParameter<int>("EventStep",1);
00039 
00040   //----- Set which tracks of those events are verbose
00041   fTVTrackMin  = p.getUntrackedParameter<int>("TrackMin",0);
00042   fTVTrackMax  = p.getUntrackedParameter<int>("TrackMax",fLarge);
00043   fTVTrackStep = p.getUntrackedParameter<int>("TrackStep",1);
00044 
00045   //----- Set the verbosity level
00046   fVerboseLevel = p.getUntrackedParameter<int>("VerboseLevel",1);
00047   fPdgIds       = p.getUntrackedParameter<std::vector<int> >("PDGids");
00048   if (fDEBUG) {
00049     G4cout << "TV: fTVTrackMin " << fTVTrackMin 
00050            << " fTVTrackMax "    <<  fTVTrackMax 
00051            <<  " fTVTrackStep "  << fTVTrackStep  
00052            << " fTVEventMin "    << fTVEventMin 
00053            << " fTVEventMax "    << fTVEventMax   
00054            << " fTVEventStep "   << fTVEventStep 
00055            << " fVerboseLevel "  << fVerboseLevel 
00056            << " fG4Verbose "     << fG4Verbose 
00057            << " PDGIds     "     << fPdgIds.size() << G4endl;
00058     for (unsigned int ii=0; ii<fPdgIds.size(); ++ii) 
00059       G4cout << "TV: PDGId[" << ii << "] = " << fPdgIds[ii] << G4endl;
00060   }
00061   
00062   //----- Set verbosity off to start
00063   fTrackingVerboseON = false;
00064   fTkVerbThisEventON = false;
00065   
00066   G4cout << " TrackingVerbose constructed " << G4endl;
00067 }
00068 
00069 TrackingVerboseAction::~TrackingVerboseAction() {}
00070 
00071 void TrackingVerboseAction::update(const BeginOfRun * run) {
00072   TrackingAction * ta = 
00073     dynamic_cast<TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
00074   theTrackingManager = ta->getTrackManager();
00075   fVerbose = G4VSteppingVerbose::GetInstance();
00076   if (fDEBUG)
00077     G4cout << " TV: Get the Tracking Manager: " << theTrackingManager
00078            << " and the SteppingVerbose: " << fVerbose << G4endl;
00079 }
00080 
00081 void TrackingVerboseAction::update(const BeginOfEvent * evt) {
00082   if (evt==0) return;
00083   const G4Event * anEvent = (*evt)();
00084   if (anEvent==0) return;
00085 
00086   //----------- Set /tracking/verbose for this event 
00087   int eventNo = anEvent->GetEventID();
00088   if (fDEBUG) G4cout << "TV: trackID: NEW EVENT " << eventNo << G4endl;
00089 
00090   fTkVerbThisEventON = false;
00091   //----- Check if event is in the selected range
00092   if (eventNo >= fTVEventMin && eventNo <= fTVEventMax) {
00093     if ((eventNo-fTVEventMin) % fTVEventStep == 0) fTkVerbThisEventON = true;
00094   }
00095 
00096   if (fDEBUG)
00097     G4cout << " TV: fTkVerbThisEventON " <<  fTkVerbThisEventON 
00098            << " fTrackingVerboseON " << fTrackingVerboseON 
00099            << " fTVEventMin " << fTVEventMin << " fTVEventMax " << fTVEventMax << G4endl;
00100   //----- check if verbosity has to be changed
00101   if ((fTkVerbThisEventON) && (!fTrackingVerboseON)) {
00102     if (fTVTrackMin == 0 && fTVTrackMax == fLarge && fTVTrackStep != 1) {
00103       setTrackingVerbose(fVerboseLevel);
00104       fTrackingVerboseON = true;
00105       if (fDEBUG) G4cout << "TV: VERBOSEet1 " << eventNo << G4endl;
00106     }
00107   } else if ((!fTkVerbThisEventON) && (fTrackingVerboseON) ) {
00108     setTrackingVerbose(0);
00109     fTrackingVerboseON = false;
00110     if (fDEBUG) G4cout << "TV: VERBOSEet0 " << eventNo << G4endl;
00111   }
00112 
00113 }
00114 
00115 void TrackingVerboseAction::update(const BeginOfTrack * trk) {
00116   const G4Track * aTrack = (*trk)();
00117 
00118   //----- High ET photon printout
00119   //---------- Set /tracking/verbose
00120   //----- track is verbose only if event is verbose
00121   double tkP = aTrack->GetMomentum().mag();
00122   double tkPx = aTrack->GetMomentum().x();
00123   double tkPy = aTrack->GetMomentum().y();
00124   double tkPz = aTrack->GetMomentum().z();
00125 
00126   double tvtx = aTrack->GetVertexPosition().x();
00127   double tvty = aTrack->GetVertexPosition().y();
00128   double tvtz = aTrack->GetVertexPosition().z();
00129 
00130   double g4t_phi=atan2(tkPy,tkPx);
00131 
00132   double drpart=sqrt(tkPx*tkPx + tkPy*tkPy);
00133 
00134   double mythetapart=acos(tkPz/sqrt(drpart*drpart+tkPz*tkPz));
00135 
00136   double g4t_eta=-log(tan(mythetapart/2.));
00137   G4int MytrackNo = aTrack->GetTrackID();
00138     
00139   if (fHighEtPhotons) {
00140     if (aTrack->GetDefinition()->GetParticleName() == "gamma" && aTrack->GetParentID() !=0) {
00141       if((tkPx*tkPx + tkPy*tkPy + tkPz*tkPz)>1000.0*1000.0 &&
00142          aTrack->GetCreatorProcess()->GetProcessName() == "LCapture") {
00143         G4cout << "MY NEW GAMMA " << G4endl;
00144         G4cout << "**********************************************************************"  << G4endl;
00145         G4cout << "MY NEW TRACK ID = " << MytrackNo << "("
00146                << aTrack->GetDefinition()->GetParticleName()
00147                <<")"<< " PARENT ="<< aTrack->GetParentID() << G4endl;
00148         G4cout << "Primary particle: " 
00149                << aTrack->GetDynamicParticle()->GetPrimaryParticle() << G4endl;
00150         G4cout << "Process type: " << aTrack->GetCreatorProcess()->GetProcessType()
00151                << " Process name: " 
00152                << aTrack->GetCreatorProcess()->GetProcessName() << G4endl;
00153         G4cout << "ToT E = " << aTrack->GetTotalEnergy() 
00154                << " KineE = " << aTrack->GetKineticEnergy()
00155                << " Tot P = " << tkP << " Pt = " << sqrt(tkPx*tkPx + tkPy*tkPy) 
00156                << " VTX=(" << tvtx << "," << tvty << "," << tvtz << ")" << G4endl;
00157         if (aTrack->GetKineticEnergy() > 1.*GeV 
00158             && aTrack->GetCreatorProcess()->GetProcessName() != "LCapture")
00159           G4cout << " KineE > 1 GeV !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << G4endl;
00160         const G4VTouchable* touchable=aTrack->GetTouchable();
00161         if (touchable!=0 && touchable->GetVolume()!=0 &&
00162             touchable->GetVolume()->GetLogicalVolume()!=0) {
00163           G4Material* material=touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
00164           G4cout << "G4LCapture Gamma E(GeV) " 
00165                  << aTrack->GetTotalEnergy()/GeV << "  "
00166                  << material->GetName()<< " " 
00167                  << touchable->GetVolume()->GetName() << G4endl;
00168           G4cout << "G4LCapture Gamma position(m): " 
00169                  << aTrack->GetPosition()/m << G4endl;
00170           G4cout << "G4LCapture created Gamma direction " 
00171                  << aTrack->GetMomentumDirection() << G4endl;
00172           G4cout << "G4LCapture gamma (eta,phi) = " 
00173                  << "(" << g4t_eta << "," << g4t_phi << ")" << G4endl;
00174         }
00175         aTrack->GetUserInformation()->Print();
00176         G4cout << "**********************************************************************"  << G4endl;
00177       }
00178     }
00179 
00180     if (aTrack->GetDefinition()->GetParticleName() == "gamma") {
00181       const G4VProcess * proc = aTrack->GetCreatorProcess();
00182       double Tgamma = aTrack->GetKineticEnergy();
00183       std::string ProcName;
00184       const  std::string nullstr ("Null_prc");
00185       if (proc) ProcName = proc->GetProcessName();
00186       else      ProcName = nullstr;
00187       if (Tgamma > 2.5*GeV ) { //&& ProcName!="Decay" && ProcName!="eBrem")
00188         std::string volumeName("_Unknown_Vol_");
00189         std::string materialName("_Unknown_Mat_");
00190         G4Material * material = 0;
00191         G4VPhysicalVolume * pvolume = 0;
00192         G4LogicalVolume * lvolume = 0;
00193         const G4VTouchable * touchable = aTrack->GetTouchable();
00194         if (touchable) pvolume = touchable->GetVolume();
00195         if (pvolume) {
00196           volumeName = pvolume->GetName();
00197           lvolume = pvolume->GetLogicalVolume();
00198         }
00199         if (lvolume) material = lvolume->GetMaterial();
00200         if (material) materialName = material->GetName();
00201         G4cout << "**** ALL photons > 2.5 GeV ****" << G4endl;
00202         G4cout << ProcName << "**** ALL photons: gamma E(GeV) "
00203                << aTrack->GetTotalEnergy()/GeV << "  "
00204                <<  materialName << " " << volumeName << G4endl;
00205         G4cout << ProcName << "**** ALL photons: gamma position(m): " 
00206                << aTrack->GetPosition()/m << G4endl;
00207         G4cout << ProcName << "**** ALL photons: gamma direction " 
00208                << aTrack->GetMomentumDirection() << G4endl;
00209         G4cout << "**********************************************************************"  << G4endl;
00210       }
00211     }                                               
00212   }
00213     
00214   //---------- Set /tracking/verbose
00215   //----- track is verbose only if event is verbose
00216   if (fTkVerbThisEventON) {
00217     bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
00218 
00219     //----- Set the /tracking/verbose for this track 
00220     if ((trackingVerboseThisTrack) && (!fTrackingVerboseON) ) {
00221       setTrackingVerbose(fVerboseLevel);
00222       fTrackingVerboseON = true;
00223       if (fDEBUG) G4cout << "TV: VERBOSEtt1 " << aTrack->GetTrackID()
00224                          << G4endl;
00225       printTrackInfo(aTrack);
00226     } else if ((!trackingVerboseThisTrack) && ( fTrackingVerboseON )) {
00227       setTrackingVerbose(0);
00228       fTrackingVerboseON = false;
00229       if (fDEBUG) G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID()
00230                          << G4endl;
00231     }
00232   }
00233 }
00234 
00235 void TrackingVerboseAction::update(const EndOfTrack * trk) {
00236   const G4Track * aTrack = (*trk)();
00237   if (fTkVerbThisEventON) {
00238     bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
00239     if ((trackingVerboseThisTrack) && (fTrackingVerboseON ) &&
00240         (fTVTrackMax < fLarge || fTVTrackStep != 1)) {
00241       setTrackingVerbose(0);
00242       fTrackingVerboseON = false;
00243       if (fDEBUG) G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID()
00244                          << G4endl;
00245     }
00246   }
00247 }
00248 
00249 void TrackingVerboseAction::update(const G4Step* fStep) {
00250 
00251   if ((fG4Verbose) && (fTrackingVerboseON)) {
00252     G4Track* fTrack = fStep->GetTrack();
00253     G4cout << std::setw( 5) << fTrack->GetCurrentStepNumber() << " "
00254            << std::setw( 8) << G4BestUnit(fTrack->GetPosition().x() , "Length") << " "
00255            << std::setw( 8) << G4BestUnit(fTrack->GetPosition().y() , "Length") << " "
00256            << std::setw( 8) << G4BestUnit(fTrack->GetPosition().z() , "Length") << " "
00257            << std::setw( 9) << G4BestUnit(fTrack->GetKineticEnergy() , "Energy") << " "
00258            << std::setw( 8) << G4BestUnit(fStep->GetTotalEnergyDeposit(), "Energy") << " "
00259            << std::setw( 8) << G4BestUnit(fStep->GetStepLength() , "Length") << " "
00260            << std::setw( 9) << G4BestUnit(fTrack->GetTrackLength() , "Length") << " "
00261            << std::setw( 9) << G4BestUnit(fTrack->GetGlobalTime(), "Time") << " ";
00262 
00263     // Put cut comment here
00264     if( fTrack->GetNextVolume() != 0 ) {
00265       G4cout << std::setw(11) << fTrack->GetNextVolume()->GetName() << " ";
00266     } else {
00267       G4cout << std::setw(11) << "OutOfWorld" << " ";
00268     }
00269     if(fStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL){
00270       G4cout << fStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
00271     } else {
00272       G4cout << "User Limit";
00273     }
00274     G4cout << G4endl;
00275   }
00276 }
00277 
00278 void TrackingVerboseAction::setTrackingVerbose(int verblev) {
00279   if (fDEBUG) G4cout << " setting verbose level " << verblev << G4endl;
00280   if (theTrackingManager!=0) theTrackingManager->SetVerboseLevel(verblev);
00281 }
00282  
00283 bool TrackingVerboseAction::checkTrackingVerbose(const G4Track* aTrack) {
00284   int trackNo = aTrack->GetTrackID();    
00285   bool trackingVerboseThisTrack = false;
00286   //----- Check if track is in the selected range
00287   if (trackNo >= fTVTrackMin && trackNo <= fTVTrackMax) {
00288     if ((trackNo-fTVTrackMin) % fTVTrackStep == 0) trackingVerboseThisTrack = true;
00289   }
00290   if (trackingVerboseThisTrack && (fPdgIds.size()>0)) {
00291     int pdgId = aTrack->GetDefinition()->GetPDGEncoding();
00292     if (std::count(fPdgIds.begin(),fPdgIds.end(),pdgId) == 0) trackingVerboseThisTrack = false;
00293   }
00294   return trackingVerboseThisTrack;
00295 }
00296 
00297 void TrackingVerboseAction::printTrackInfo(const G4Track* aTrack) {
00298   G4cout << G4endl
00299          << "*******************************************************"
00300          << "**************************************************" << G4endl
00301          << "* G4Track Information: "
00302          << "  Particle = " << aTrack->GetDefinition()->GetParticleName()
00303          << ","
00304          << "   Track ID = " << aTrack->GetTrackID()
00305          << ","
00306          << "   Parent ID = " << aTrack->GetParentID() << G4endl
00307          << "*******************************************************"
00308          << "**************************************************"
00309          << G4endl << G4endl;
00310   if (fVerbose) fVerbose->TrackingStarted();
00311 }
/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4Core/TrackingVerbose/src/TrackingVerboseAction.cc
