CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackingVerboseAction.cc
Go to the documentation of this file.
1 // File: TrackingVerboseAction.cc
3 // Creation: P.Arce 09/01
4 // Modifications: porting to CMSSW by M. Stavrianakou 22/03/06
5 // Description:
7 
8 #include "SimG4Core/TrackingVerbose/interface/TrackingVerboseAction.h"
9 
10 #include "SimG4Core/Notification/interface/BeginOfRun.h"
11 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
12 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
13 #include "SimG4Core/Notification/interface/EndOfTrack.h"
14 #include "SimG4Core/Notification/interface/TrackInformation.h"
15 #include "SimG4Core/Application/interface/TrackingAction.h"
16 
17 #include "G4Track.hh"
18 #include "G4Event.hh"
19 #include "G4ios.hh"
20 #include "G4TrackingManager.hh"
21 #include "G4EventManager.hh"
22 #include "G4VSteppingVerbose.hh"
23 #include "G4UnitsTable.hh"
24 
25 #include<algorithm>
26 
27 TrackingVerboseAction::TrackingVerboseAction(edm::ParameterSet const & p) :
28  theTrackingManager(0), fVerbose(0) {
29 
30  fLarge = int(1E10);
31  fDEBUG = p.getUntrackedParameter<bool>("DEBUG",false);
32  fHighEtPhotons = p.getUntrackedParameter<bool>("CheckForHighEtPhotons",false);
33  fG4Verbose = p.getUntrackedParameter<bool>("G4Verbose",false);
34 
35  //----- Set which events are verbose
36  fTVEventMin = p.getUntrackedParameter<int>("EventMin",0);
37  fTVEventMax = p.getUntrackedParameter<int>("EventMax",fLarge);
38  fTVEventStep = p.getUntrackedParameter<int>("EventStep",1);
39 
40  //----- Set which tracks of those events are verbose
41  fTVTrackMin = p.getUntrackedParameter<int>("TrackMin",0);
42  fTVTrackMax = p.getUntrackedParameter<int>("TrackMax",fLarge);
43  fTVTrackStep = p.getUntrackedParameter<int>("TrackStep",1);
44 
45  //----- Set the verbosity level
46  fVerboseLevel = p.getUntrackedParameter<int>("VerboseLevel",1);
47  fPdgIds = p.getUntrackedParameter<std::vector<int> >("PDGids");
48  if (fDEBUG) {
49  G4cout << "TV: fTVTrackMin " << fTVTrackMin
50  << " fTVTrackMax " << fTVTrackMax
51  << " fTVTrackStep " << fTVTrackStep
52  << " fTVEventMin " << fTVEventMin
53  << " fTVEventMax " << fTVEventMax
54  << " fTVEventStep " << fTVEventStep
55  << " fVerboseLevel " << fVerboseLevel
56  << " fG4Verbose " << fG4Verbose
57  << " PDGIds " << fPdgIds.size() << G4endl;
58  for (unsigned int ii=0; ii<fPdgIds.size(); ++ii)
59  G4cout << "TV: PDGId[" << ii << "] = " << fPdgIds[ii] << G4endl;
60  }
61 
62  //----- Set verbosity off to start
63  fTrackingVerboseON = false;
64  fTkVerbThisEventON = false;
65 
66  G4cout << " TrackingVerbose constructed " << G4endl;
67 }
68 
69 TrackingVerboseAction::~TrackingVerboseAction() {}
70 
71 void TrackingVerboseAction::update(const BeginOfRun * run) {
72  TrackingAction * ta =
73  dynamic_cast<TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
74  theTrackingManager = ta->getTrackManager();
75  fVerbose = G4VSteppingVerbose::GetInstance();
76  if (fDEBUG)
77  G4cout << " TV: Get the Tracking Manager: " << theTrackingManager
78  << " and the SteppingVerbose: " << fVerbose << G4endl;
79 }
80 
81 void TrackingVerboseAction::update(const BeginOfEvent * evt) {
82  if (evt==0) return;
83  const G4Event * anEvent = (*evt)();
84  if (anEvent==0) return;
85 
86  //----------- Set /tracking/verbose for this event
87  int eventNo = anEvent->GetEventID();
88  if (fDEBUG) G4cout << "TV: trackID: NEW EVENT " << eventNo << G4endl;
89 
90  fTkVerbThisEventON = false;
91  //----- Check if event is in the selected range
92  if (eventNo >= fTVEventMin && eventNo <= fTVEventMax) {
93  if ((eventNo-fTVEventMin) % fTVEventStep == 0) fTkVerbThisEventON = true;
94  }
95 
96  if (fDEBUG)
97  G4cout << " TV: fTkVerbThisEventON " << fTkVerbThisEventON
98  << " fTrackingVerboseON " << fTrackingVerboseON
99  << " fTVEventMin " << fTVEventMin << " fTVEventMax " << fTVEventMax << G4endl;
100  //----- check if verbosity has to be changed
101  if ((fTkVerbThisEventON) && (!fTrackingVerboseON)) {
102  if (fTVTrackMin == 0 && fTVTrackMax == fLarge && fTVTrackStep != 1) {
103  setTrackingVerbose(fVerboseLevel);
104  fTrackingVerboseON = true;
105  if (fDEBUG) G4cout << "TV: VERBOSEet1 " << eventNo << G4endl;
106  }
107  } else if ((!fTkVerbThisEventON) && (fTrackingVerboseON) ) {
108  setTrackingVerbose(0);
109  fTrackingVerboseON = false;
110  if (fDEBUG) G4cout << "TV: VERBOSEet0 " << eventNo << G4endl;
111  }
112 
113 }
114 
115 void TrackingVerboseAction::update(const BeginOfTrack * trk) {
116  const G4Track * aTrack = (*trk)();
117 
118  //----- High ET photon printout
119  //---------- Set /tracking/verbose
120  //----- track is verbose only if event is verbose
121  double tkP = aTrack->GetMomentum().mag();
122  double tkPx = aTrack->GetMomentum().x();
123  double tkPy = aTrack->GetMomentum().y();
124  double tkPz = aTrack->GetMomentum().z();
125 
126  double tvtx = aTrack->GetVertexPosition().x();
127  double tvty = aTrack->GetVertexPosition().y();
128  double tvtz = aTrack->GetVertexPosition().z();
129 
130  double g4t_phi=atan2(tkPy,tkPx);
131 
132  double drpart=sqrt(tkPx*tkPx + tkPy*tkPy);
133 
134  double mythetapart=acos(tkPz/sqrt(drpart*drpart+tkPz*tkPz));
135 
136  double g4t_eta=-log(tan(mythetapart/2.));
137  G4int MytrackNo = aTrack->GetTrackID();
138 
139  if (fHighEtPhotons) {
140  if (aTrack->GetDefinition()->GetParticleName() == "gamma" && aTrack->GetParentID() !=0) {
141  if((tkPx*tkPx + tkPy*tkPy + tkPz*tkPz)>1000.0*1000.0 &&
142  aTrack->GetCreatorProcess()->GetProcessName() == "LCapture") {
143  G4cout << "MY NEW GAMMA " << G4endl;
144  G4cout << "**********************************************************************" << G4endl;
145  G4cout << "MY NEW TRACK ID = " << MytrackNo << "("
146  << aTrack->GetDefinition()->GetParticleName()
147  <<")"<< " PARENT ="<< aTrack->GetParentID() << G4endl;
148  G4cout << "Primary particle: "
149  << aTrack->GetDynamicParticle()->GetPrimaryParticle() << G4endl;
150  G4cout << "Process type: " << aTrack->GetCreatorProcess()->GetProcessType()
151  << " Process name: "
152  << aTrack->GetCreatorProcess()->GetProcessName() << G4endl;
153  G4cout << "ToT E = " << aTrack->GetTotalEnergy()
154  << " KineE = " << aTrack->GetKineticEnergy()
155  << " Tot P = " << tkP << " Pt = " << sqrt(tkPx*tkPx + tkPy*tkPy)
156  << " VTX=(" << tvtx << "," << tvty << "," << tvtz << ")" << G4endl;
157  if (aTrack->GetKineticEnergy() > 1.*GeV
158  && aTrack->GetCreatorProcess()->GetProcessName() != "LCapture")
159  G4cout << " KineE > 1 GeV !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << G4endl;
160  const G4VTouchable* touchable=aTrack->GetTouchable();
161  if (touchable!=0 && touchable->GetVolume()!=0 &&
162  touchable->GetVolume()->GetLogicalVolume()!=0) {
163  G4Material* material=touchable->GetVolume()->GetLogicalVolume()->GetMaterial();
164  G4cout << "G4LCapture Gamma E(GeV) "
165  << aTrack->GetTotalEnergy()/GeV << " "
166  << material->GetName()<< " "
167  << touchable->GetVolume()->GetName() << G4endl;
168  G4cout << "G4LCapture Gamma position(m): "
169  << aTrack->GetPosition()/m << G4endl;
170  G4cout << "G4LCapture created Gamma direction "
171  << aTrack->GetMomentumDirection() << G4endl;
172  G4cout << "G4LCapture gamma (eta,phi) = "
173  << "(" << g4t_eta << "," << g4t_phi << ")" << G4endl;
174  }
175  aTrack->GetUserInformation()->Print();
176  G4cout << "**********************************************************************" << G4endl;
177  }
178  }
179 
180  if (aTrack->GetDefinition()->GetParticleName() == "gamma") {
181  const G4VProcess * proc = aTrack->GetCreatorProcess();
182  double Tgamma = aTrack->GetKineticEnergy();
183  std::string ProcName;
184  const std::string nullstr ("Null_prc");
185  if (proc) ProcName = proc->GetProcessName();
186  else ProcName = nullstr;
187  if (Tgamma > 2.5*GeV ) { //&& ProcName!="Decay" && ProcName!="eBrem")
188  std::string volumeName("_Unknown_Vol_");
189  std::string materialName("_Unknown_Mat_");
190  G4Material * material = 0;
191  G4VPhysicalVolume * pvolume = 0;
192  G4LogicalVolume * lvolume = 0;
193  const G4VTouchable * touchable = aTrack->GetTouchable();
194  if (touchable) pvolume = touchable->GetVolume();
195  if (pvolume) {
196  volumeName = pvolume->GetName();
197  lvolume = pvolume->GetLogicalVolume();
198  }
199  if (lvolume) material = lvolume->GetMaterial();
200  if (material) materialName = material->GetName();
201  G4cout << "**** ALL photons > 2.5 GeV ****" << G4endl;
202  G4cout << ProcName << "**** ALL photons: gamma E(GeV) "
203  << aTrack->GetTotalEnergy()/GeV << " "
204  << materialName << " " << volumeName << G4endl;
205  G4cout << ProcName << "**** ALL photons: gamma position(m): "
206  << aTrack->GetPosition()/m << G4endl;
207  G4cout << ProcName << "**** ALL photons: gamma direction "
208  << aTrack->GetMomentumDirection() << G4endl;
209  G4cout << "**********************************************************************" << G4endl;
210  }
211  }
212  }
213 
214  //---------- Set /tracking/verbose
215  //----- track is verbose only if event is verbose
216  if (fTkVerbThisEventON) {
217  bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
218 
219  //----- Set the /tracking/verbose for this track
220  if ((trackingVerboseThisTrack) && (!fTrackingVerboseON) ) {
221  setTrackingVerbose(fVerboseLevel);
222  fTrackingVerboseON = true;
223  if (fDEBUG) G4cout << "TV: VERBOSEtt1 " << aTrack->GetTrackID()
224  << G4endl;
225  printTrackInfo(aTrack);
226  } else if ((!trackingVerboseThisTrack) && ( fTrackingVerboseON )) {
227  setTrackingVerbose(0);
228  fTrackingVerboseON = false;
229  if (fDEBUG) G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID()
230  << G4endl;
231  }
232  }
233 }
234 
235 void TrackingVerboseAction::update(const EndOfTrack * trk) {
236  const G4Track * aTrack = (*trk)();
237  if (fTkVerbThisEventON) {
238  bool trackingVerboseThisTrack = checkTrackingVerbose(aTrack);
239  if ((trackingVerboseThisTrack) && (fTrackingVerboseON ) &&
240  (fTVTrackMax < fLarge || fTVTrackStep != 1)) {
241  setTrackingVerbose(0);
242  fTrackingVerboseON = false;
243  if (fDEBUG) G4cout << "TV: VERBOSEtt0 " << aTrack->GetTrackID()
244  << G4endl;
245  }
246  }
247 }
248 
249 void TrackingVerboseAction::update(const G4Step* fStep) {
250 
251  if ((fG4Verbose) && (fTrackingVerboseON)) {
252  G4Track* fTrack = fStep->GetTrack();
253  G4cout << std::setw( 5) << fTrack->GetCurrentStepNumber() << " "
254  << std::setw( 8) << G4BestUnit(fTrack->GetPosition().x() , "Length") << " "
255  << std::setw( 8) << G4BestUnit(fTrack->GetPosition().y() , "Length") << " "
256  << std::setw( 8) << G4BestUnit(fTrack->GetPosition().z() , "Length") << " "
257  << std::setw( 9) << G4BestUnit(fTrack->GetKineticEnergy() , "Energy") << " "
258  << std::setw( 8) << G4BestUnit(fStep->GetTotalEnergyDeposit(), "Energy") << " "
259  << std::setw( 8) << G4BestUnit(fStep->GetStepLength() , "Length") << " "
260  << std::setw( 9) << G4BestUnit(fTrack->GetTrackLength() , "Length") << " "
261  << std::setw( 9) << G4BestUnit(fTrack->GetGlobalTime(), "Time") << " ";
262 
263  // Put cut comment here
264  if( fTrack->GetNextVolume() != 0 ) {
265  G4cout << std::setw(11) << fTrack->GetNextVolume()->GetName() << " ";
266  } else {
267  G4cout << std::setw(11) << "OutOfWorld" << " ";
268  }
269  if(fStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL){
270  G4cout << fStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
271  } else {
272  G4cout << "User Limit";
273  }
274  G4cout << G4endl;
275  }
276 }
277 
278 void TrackingVerboseAction::setTrackingVerbose(int verblev) {
279  if (fDEBUG) G4cout << " setting verbose level " << verblev << G4endl;
280  if (theTrackingManager!=0) theTrackingManager->SetVerboseLevel(verblev);
281 }
282 
283 bool TrackingVerboseAction::checkTrackingVerbose(const G4Track* aTrack) {
284  int trackNo = aTrack->GetTrackID();
285  bool trackingVerboseThisTrack = false;
286  //----- Check if track is in the selected range
287  if (trackNo >= fTVTrackMin && trackNo <= fTVTrackMax) {
288  if ((trackNo-fTVTrackMin) % fTVTrackStep == 0) trackingVerboseThisTrack = true;
289  }
290  if (trackingVerboseThisTrack && (fPdgIds.size()>0)) {
291  int pdgId = aTrack->GetDefinition()->GetPDGEncoding();
292  if (std::count(fPdgIds.begin(),fPdgIds.end(),pdgId) == 0) trackingVerboseThisTrack = false;
293  }
294  return trackingVerboseThisTrack;
295 }
296 
297 void TrackingVerboseAction::printTrackInfo(const G4Track* aTrack) {
298  G4cout << G4endl
299  << "*******************************************************"
300  << "**************************************************" << G4endl
301  << "* G4Track Information: "
302  << " Particle = " << aTrack->GetDefinition()->GetParticleName()
303  << ","
304  << " Track ID = " << aTrack->GetTrackID()
305  << ","
306  << " Parent ID = " << aTrack->GetParentID() << G4endl
307  << "*******************************************************"
308  << "**************************************************"
309  << G4endl << G4endl;
310  if (fVerbose) fVerbose->TrackingStarted();
311 }
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
proc
TrainProcessor *const proc
Definition: MVATrainer.cc:101
TrackingVerboseAction::fPdgIds
std::vector< int > fPdgIds
Definition: TrackingVerboseAction.h:70
TrackingVerboseAction::checkTrackingVerbose
bool checkTrackingVerbose(const G4Track *)
Definition: TrackingVerboseAction.cc:283
NULL
#define NULL
Definition: scimark2.h:8
TrackingVerboseAction::printTrackInfo
void printTrackInfo(const G4Track *)
Definition: TrackingVerboseAction.cc:297
TrackingVerboseAction::theTrackingManager
G4TrackingManager * theTrackingManager
Definition: TrackingVerboseAction.h:71
TrackingAction::getTrackManager
G4TrackingManager * getTrackManager()
Definition: TrackingAction.cc:118
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
TrackingVerboseAction::fVerbose
G4VSteppingVerbose * fVerbose
Definition: TrackingVerboseAction.h:72
funct::log
Log< T >::type log(const T &t)
Definition: Log.h:22
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
TrackingVerboseAction::update
void update(const BeginOfRun *)
This routine will be called when the appropriate signal arrives.
Definition: TrackingVerboseAction.cc:71
benchmark_cfg.pdgId
tuple pdgId
Definition: benchmark_cfg.py:18
TrackingVerboseAction::setTrackingVerbose
void setTrackingVerbose(int verblev)
Definition: TrackingVerboseAction.cc:278
TrackingVerboseAction::TrackingVerboseAction
TrackingVerboseAction(edm::ParameterSet const &p)
Definition: TrackingVerboseAction.cc:27
Generated for CMSSW Reference Manual by doxygen 1.8.5