Skip to content
Snippets Groups Projects
Commit 878512fb authored by Chad Lantz's avatar Chad Lantz
Browse files

Simplified the handling of hits and output vectors

parent ca58bf0e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
include/EventAction.hh
+ 0
− 10
View file @ 878512fb
......@@ -29,8 +29,6 @@
#include "globals.hh"
#include "G4UserEventAction.hh"
#include <vector>
class EventAction : public G4UserEventAction{
public:
......@@ -40,14 +38,6 @@ class EventAction : public G4UserEventAction{
virtual void BeginOfEventAction( const G4Event* event );
virtual void EndOfEventAction( const G4Event* event );
inline std::vector< std::vector<double>* > GetVectors( ){return fPtrVec;}
private:
G4int hitsCollID;
G4int fEventNo;
std::vector<double> x,z,xHit,zHit,time,theta,energy;
std::vector< std::vector<double>* > fPtrVec;
};
#endif
include/PMTSD.hh
+ 3
− 0
View file @ 878512fb
......@@ -48,6 +48,9 @@ public:
private:
int HCID;
HitsCollection* hitCollection;
std::vector< std::vector<double>* > fPtrVec;
bool VISUALIZE;
};
#endif
include/RunAction.hh
+ 6
− 0
View file @ 878512fb
......@@ -37,6 +37,8 @@
#include "globals.hh"
#include "G4UserRunAction.hh"
#include <vector>
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class G4Timer;
......@@ -52,9 +54,13 @@ class RunAction : public G4UserRunAction
virtual void BeginOfRunAction(const G4Run* aRun);
virtual void EndOfRunAction(const G4Run* aRun);
std::vector< std::vector<double>* > GetVectors(){ return fPtrVec; }
private:
G4Timer* fTimer;
G4String m_fileName;
std::vector< std::vector<double>* > fPtrVec;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......
src/EventAction.cc
+ 10
− 40
View file @ 878512fb
......@@ -41,62 +41,32 @@
*
*/
EventAction::EventAction()
: G4UserEventAction(),
hitsCollID(0){
fPtrVec.push_back(&x);
fPtrVec.push_back(&z);
fPtrVec.push_back(&xHit);
fPtrVec.push_back(&zHit);
fPtrVec.push_back(&time);
fPtrVec.push_back(&theta);
fPtrVec.push_back(&energy);
: G4UserEventAction()
{
}
/*
*
*/
EventAction::~EventAction(){
EventAction::~EventAction()
{
}
/*
*
*/
void EventAction::BeginOfEventAction(const G4Event* event){
fEventNo = event->GetEventID();
if(fEventNo%100000 == 0) G4cout << "Begin Event " << fEventNo << G4endl;
hitsCollID = 0;
for(uint i = 0; i < fPtrVec.size(); i++){
fPtrVec.at(i)->clear();
}
if(event->GetEventID()%100000 == 0) G4cout << "Begin Event " << event->GetEventID() << G4endl;
}
/*
*
*/
void EventAction::EndOfEventAction(const G4Event* event){
hitsCollID = G4SDManager::GetSDMpointer()->GetCollectionID("MyPMT");
G4HCofThisEvent* HCE = event->GetHCofThisEvent();
if(HCE){
HitsCollection* HC = (HitsCollection*)(HCE->GetHC(hitsCollID));
if( HC->entries() ){
for(size_t i = 0; i < HC->entries(); i++ ){
PMTHit* aHit = (*HC)[i];
// fill vectors //
x.push_back ( aHit->getPos().x() );
z.push_back ( aHit->getPos().z() );
xHit.push_back ( aHit->getHit().x() );
zHit.push_back ( aHit->getHit().z() );
time.push_back ( aHit->getTime() );
energy.push_back( aHit->getEnergy() );
G4double pTheta = atan(sqrt(pow(aHit->getMomentum().x(),2) + pow(aHit->getMomentum().z(),2) )/fabs(aHit->getMomentum().y() ));
theta.push_back( pTheta );
}
// fill ntuple //
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
analysisManager->AddNtupleRow();
}
}
// fill ntuple //
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
analysisManager->FillNtupleIColumn( 0, event->GetEventID());
analysisManager->AddNtupleRow();
}
src/PMTSD.cc
+ 33
− 10
View file @ 878512fb
......@@ -26,11 +26,14 @@
#include "PMTSD.hh"
#include "PMTHit.hh"
#include "RunAction.hh"
#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
#include "G4ThreeVector.hh"
#include "G4SDManager.hh"
#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4ios.hh"
#include "G4Poisson.hh"
#include "G4ParticleTypes.hh"
......@@ -50,6 +53,7 @@ PMTSD::PMTSD(G4String sdName)
:G4VSensitiveDetector(sdName){
collectionName.insert(sdName);
HCID = -1;
VISUALIZE = false;
}
......@@ -75,6 +79,14 @@ void PMTSD::Initialize(G4HCofThisEvent* HCE){
{ HCID = G4SDManager::GetSDMpointer()->GetCollectionID( name );}
HCE->AddHitsCollection( HCID, hitCollection );
// If the UI window exists, set visualize to true so we record and draw hits
if( G4UImanager::GetUIpointer()->GetG4UIWindow() ) VISUALIZE = true;
// Get the vectors from run action
RunAction* runAction = (RunAction*)G4RunManager::GetRunManager()->GetUserRunAction();
fPtrVec = runAction->GetVectors();
}
......@@ -82,18 +94,29 @@ void PMTSD::Initialize(G4HCofThisEvent* HCE){
*
*/
G4bool PMTSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
G4Track* theTrack = aStep->GetTrack();
PMTHit* newHit = new PMTHit();
newHit->setTrackID (theTrack->GetTrackID() );
newHit->setPos (theTrack->GetVertexPosition() );
newHit->setHit (theTrack->GetPosition() );
newHit->setEnergy (aStep->GetPreStepPoint()->GetTotalEnergy());
newHit->setMomentum (aStep->GetPreStepPoint()->GetMomentum());
newHit->setTime (theTrack->GetGlobalTime());
hitCollection->insert( newHit );
G4ThreeVector origin = theTrack->GetVertexPosition();
G4ThreeVector hitPos = theTrack->GetPosition();
G4ThreeVector momentum = aStep->GetPreStepPoint()->GetMomentum();
//Display hits if the UI is on
if(VISUALIZE){
PMTHit* newHit = new PMTHit();
newHit->setPos( theTrack->GetPosition() );
hitCollection->insert ( newHit );
}
// Fill ouput vectors
// These vectors are created in RunAction
fPtrVec[0]->push_back( origin.x() ); //Origin of the photon in x
fPtrVec[1]->push_back( origin.z() ); //Origin of the photon in z
fPtrVec[2]->push_back( hitPos.x() ); //Where the photon landed on the PMT in x
fPtrVec[3]->push_back( hitPos.z() ); //Where the photon landed on the PMT in z
fPtrVec[4]->push_back( theTrack->GetGlobalTime() ); //Time of arrival of the photon at the PMT
G4double pTheta = atan(sqrt(pow(momentum.x(),2) + pow(momentum.z(),2) )/fabs(momentum.y() ));
fPtrVec[5]->push_back( pTheta ); //Incident angle at the PMT about the y-axis
fPtrVec[6]->push_back( aStep->GetPreStepPoint()->GetTotalEnergy() ); //Photon energy
return true;
}
......
src/RunAction.cc
+ 25
− 32
View file @ 878512fb
......@@ -52,37 +52,6 @@ RunAction::RunAction(G4String fileName)
G4RunManager::GetRunManager()->SetPrintProgress(1);
// Create analysis manager. The choice of analysis
// technology is done via selection of a namespace
// in B4Analysis.hh
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
// Get event action
const EventAction* constEventAction
= static_cast<const EventAction*>(G4RunManager::GetRunManager()
->GetUserEventAction());
EventAction* eventAction
= const_cast<EventAction*>(constEventAction);
std::vector< std::vector<double>* > pVec = eventAction->GetVectors( );
// Create directories
analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
// Creating ntuple
analysisManager->CreateNtuple("lightGuide", "pos and momentum");
analysisManager->CreateNtupleDColumn( "X", *pVec[0] );
analysisManager->CreateNtupleDColumn( "Z", *pVec[1] );
analysisManager->CreateNtupleDColumn( "hitX", *pVec[2] );
analysisManager->CreateNtupleDColumn( "hitZ", *pVec[3] );
analysisManager->CreateNtupleDColumn( "time", *pVec[4] );
analysisManager->CreateNtupleDColumn( "theta", *pVec[5] );
analysisManager->CreateNtupleDColumn( "energy", *pVec[6] );
analysisManager->CreateNtupleIColumn( "EventNo" );
analysisManager->FinishNtuple();
}
/*
......@@ -91,6 +60,7 @@ RunAction::RunAction(G4String fileName)
RunAction::~RunAction(){
delete fTimer;
delete G4AnalysisManager::Instance();
for(auto vec : fPtrVec ) delete vec;
}
/*
......@@ -100,12 +70,35 @@ void RunAction::BeginOfRunAction(const G4Run* aRun){
G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
fTimer->Start();
// Get analysis manager
// Create analysis manager. The choice of analysis
// technology is done via selection of a namespace
// in B4Analysis.hh
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
// Open an output file
if(m_fileName == "") m_fileName = "output";
analysisManager->OpenFile(m_fileName);
// Create directories
analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
// Creating ntuple
// Populate the vector with valid pointers and assign them to branches
// These vectors are filled in PMTSD
for(int i = 0; i < 7; i ++ ) fPtrVec.push_back( new std::vector< double > );
analysisManager->CreateNtuple("lightGuide", "pos and momentum");
analysisManager->CreateNtupleIColumn( "EventNo" );
analysisManager->CreateNtupleDColumn( "X", *fPtrVec[0] );
analysisManager->CreateNtupleDColumn( "Z", *fPtrVec[1] );
analysisManager->CreateNtupleDColumn( "hitX", *fPtrVec[2] );
analysisManager->CreateNtupleDColumn( "hitZ", *fPtrVec[3] );
analysisManager->CreateNtupleDColumn( "time", *fPtrVec[4] );
analysisManager->CreateNtupleDColumn( "theta", *fPtrVec[5] );
analysisManager->CreateNtupleDColumn( "energy", *fPtrVec[6] );
analysisManager->FinishNtuple();
}
/*
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment