Commit 1084409b authored by Chad Lantz's avatar Chad Lantz
Browse files

Created methods for handling OPTICAL and CLUSTER conditions.

parent 27b73860
......@@ -32,6 +32,7 @@
#define EventAction_h 1
#include "G4UserEventAction.hh"
#include "FiberHit.hh"
#include "globals.hh"
#include <vector>
......@@ -45,24 +46,26 @@ class EventAction : public G4UserEventAction
EventAction( );
virtual ~EventAction();
virtual void BeginOfEventAction (const G4Event* event);
virtual void EndOfEventAction (const G4Event* event);
virtual void SetnZDCs (G4int nZDCs);
virtual void SetnRPDs (G4int nRPDs);
inline void SetClusterFlag (G4bool arg){CLUSTER = arg;}
virtual void BeginOfEventAction ( const G4Event* event );
virtual void EndOfEventAction ( const G4Event* event );
virtual void ProcessHitCollection ( FiberHitsCollection* HC );
virtual void ProcessOpticalHitCollection ( FiberHitsCollection* HC );
virtual void SetnZDCs ( G4int nZDCs );
virtual void SetnRPDs ( G4int nRPDs );
inline void SetClusterFlag ( G4bool arg ){CLUSTER = arg;}
inline std::vector< std::vector< std::vector<double>* > >* GetRPDdoubleVectors( ){return fRPDdblVec;}
inline std::vector< std::vector< std::vector< int >* > >* GetRPDintVectors ( ){return fRPDintVec;}
inline std::vector< std::vector< std::vector<double>* > >* GetZDCdoubleVectors( ){return fZDCdblVec;}
inline std::vector< std::vector< std::vector< int >* > >* GetZDCintVectors ( ){return fZDCintVec;}
inline std::vector< std::vector< std::vector<double> > >* GetRPDdoubleVectors( ){return fRPDdblVec;}
inline std::vector< std::vector< std::vector< int > > >* GetRPDintVectors ( ){return fRPDintVec;}
inline std::vector< std::vector< std::vector<double> > >* GetZDCdoubleVectors( ){return fZDCdblVec;}
inline std::vector< std::vector< std::vector< int > > >* GetZDCintVectors ( ){return fZDCintVec;}
private:
// Indecies are [mod#][dataSet][dataValue]
// data set order will need to be coordinated betweeen here and Run Action
std::vector< std::vector< std::vector<double>* > > *fRPDdblVec, *fZDCdblVec;
std::vector< std::vector< std::vector< int >* > > *fRPDintVec, *fZDCintVec;
std::vector< std::vector< std::vector<double> > > *fRPDdblVec, *fZDCdblVec;
std::vector< std::vector< std::vector< int > > > *fRPDintVec, *fZDCintVec;
G4int hitsCollID, fEventNo, fnZDCs, fnRPDs;
G4double gunPosX, gunPosY, gunPosZ;
......
......@@ -31,6 +31,8 @@
#ifndef RunAction_h
#define RunAction_h 1
#include "Analysis.hh"
#include "G4UserRunAction.hh"
#include "globals.hh"
#include "G4ThreeVector.hh"
......@@ -55,10 +57,19 @@ class RunAction : public G4UserRunAction
virtual void BeginOfRunAction(const G4Run*);
virtual void EndOfRunAction(const G4Run*);
virtual void MakeZDCTree ( G4int nTupleNo, G4int zdcNo );
virtual void MakeZDCOpticalTree( G4int nTupleNo, G4int zdcNo );
virtual void MakeRPDTree ( G4int nTupleNo, G4int rpdNo );
virtual void MakeRPDOpticalTree( G4int nTupleNo, G4int rpdNo );
private:
G4String m_fileName;
G4bool CLUSTER;
G4AnalysisManager* m_analysisManager;
std::vector< std::vector< std::vector<double> > > *m_ZDCdblVec, *m_RPDdblVec;
std::vector< std::vector< std::vector< int > > > *m_ZDCintVec, *m_RPDintVec;
};
......
......@@ -30,6 +30,7 @@
#include "EventAction.hh"
#include "FiberHit.hh"
#include "FiberSD.hh"
#include "Analysis.hh"
#include "G4UnitsTable.hh"
......@@ -39,6 +40,9 @@
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4PrimaryVertex.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4SystemOfUnits.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -57,24 +61,13 @@ EventAction::~EventAction()
void EventAction::SetnZDCs(G4int nZDCs){
//Make a new vector of vectors of pointers to vectors of doubles
fZDCdblVec = new std::vector< std::vector< std::vector<double>* > >(nZDCs);
fZDCintVec = new std::vector< std::vector< std::vector< int >* > >(nZDCs);
//If cluster flag is false, we have more vectors to store
G4int nDblVecs = (!CLUSTER) ? 11 : 0;
G4int nIntVecs = (!CLUSTER) ? 7 : 2;
fZDCdblVec = new std::vector< std::vector< std::vector<double> > >(nZDCs);
fZDCintVec = new std::vector< std::vector< std::vector< int > > >(nZDCs);
// Resize to the largest number of branches we will fill from vectors
for(G4int i = 0; i < nZDCs; i++){
fZDCdblVec->at(i).resize(nDblVecs);
for(G4int j = 0; j < nDblVecs; j++){
fZDCdblVec->at(i).at(j) = new std::vector<double>;
}//end double vector loop
fZDCintVec->at(i).resize(nIntVecs);
for(G4int j = 0; j < nIntVecs; j++){
fZDCintVec->at(i).at(j) = new std::vector< int >;
}//end int vector loop
fZDCdblVec->at(i).resize(10);
fZDCintVec->at(i).resize(8);
}//end ZDC loop
}
......@@ -83,24 +76,13 @@ void EventAction::SetnZDCs(G4int nZDCs){
void EventAction::SetnRPDs(G4int nRPDs){
//Make a new vector of vectors of pointers to vectors of doubles
fRPDdblVec = new std::vector< std::vector< std::vector<double>* > >(nRPDs);
fRPDintVec = new std::vector< std::vector< std::vector< int >* > >(nRPDs);
//If cluster flag is false, we have more vectors to store
G4int nDblVecs = (!CLUSTER) ? 11 : 3;
G4int nIntVecs = (!CLUSTER) ? 7 : 2;
fRPDdblVec = new std::vector< std::vector< std::vector<double> > >(nRPDs);
fRPDintVec = new std::vector< std::vector< std::vector< int > > >(nRPDs);
for(G4int i = 0; i < nRPDs; i++){
fRPDdblVec->at(i).resize(10);
fRPDintVec->at(i).resize(8);
fRPDdblVec->at(i).resize(nDblVecs);
for(G4int j = 0; j < nDblVecs; j++){
fRPDdblVec->at(i).at(j) = new std::vector<double>;
}//end double vector loop
fRPDintVec->at(i).resize(nIntVecs);
for(G4int j = 0; j < nIntVecs; j++){
fRPDintVec->at(i).at(j) = new std::vector< int >;
}//end int vector loop
}//end RPD loop
}
......@@ -130,150 +112,216 @@ void EventAction::EndOfEventAction(const G4Event* evt){
HC = (FiberHitsCollection*)(HCE->GetHC(hitsCollID));
G4String name = HC->GetSDname();
G4int prevTrackId = 0;
G4int prevRadiatorNo = 0;
G4int TotalnCherenkovs = 0;
G4double TotaleDep = 0.0;
int n_hit = HC->entries();
for ( G4int i = 0 ; i < n_hit; i++){
G4ThreeVector position = (*HC)[i]->getPos();
G4ThreeVector momentum = (*HC)[i]->getMomentum();
G4double energy = (*HC)[i]->getEnergy();
G4double velocity = (*HC)[i]->getVelocity();
G4double beta = (*HC)[i]->getBeta();
G4double eDep = (*HC)[i]->getEdep();
G4double charge = (*HC)[i]->getCharge();
G4int radiatorNo = (*HC)[i]->getRadNb();
G4int rodNo = (*HC)[i]->getRodNb();
G4int modNb = (*HC)[i]->getModNb();
G4int trackID = (*HC)[i]->getTrackID();
G4int pid = (*HC)[i]->getParticle()->GetPDGEncoding();
G4int nCherenkovs = (*HC)[i]->getNCherenkovs(); // This is the number of cherenkovs in a single step within the SD
//Sum energy from all steps a particle takes in a single scoring volume
if (trackID == prevTrackId || radiatorNo == prevRadiatorNo || eDep != 0) {
TotalnCherenkovs += nCherenkovs;
TotaleDep += eDep;
prevTrackId = trackID;
prevRadiatorNo = radiatorNo;
continue;
}// end summation
if(name.compare(0,3,"RPD")){ //RPD hits
int rpdNo = atoi( name.substr(3,1).c_str() );
if(!CLUSTER){
//doubles
fRPDdblVec->at(rpdNo).at(0)-> push_back( position.x() );
fRPDdblVec->at(rpdNo).at(1)-> push_back( position.y() );
fRPDdblVec->at(rpdNo).at(2)-> push_back( position.z() );
fRPDdblVec->at(rpdNo).at(3)-> push_back( momentum.x() );
fRPDdblVec->at(rpdNo).at(4)-> push_back( momentum.y() );
fRPDdblVec->at(rpdNo).at(5)-> push_back( momentum.z() );
fRPDdblVec->at(rpdNo).at(6)-> push_back( energy );
fRPDdblVec->at(rpdNo).at(7)-> push_back( velocity );
fRPDdblVec->at(rpdNo).at(8)-> push_back( beta );
fRPDdblVec->at(rpdNo).at(9)-> push_back( TotaleDep );
fRPDdblVec->at(rpdNo).at(10)->push_back( charge );
//ints
fRPDintVec->at(rpdNo).at(0)->push_back( fEventNo );
fRPDintVec->at(rpdNo).at(1)->push_back( modNb );
fRPDintVec->at(rpdNo).at(2)->push_back( radiatorNo );
fRPDintVec->at(rpdNo).at(3)->push_back( rodNo );
fRPDintVec->at(rpdNo).at(4)->push_back( TotalnCherenkovs );
fRPDintVec->at(rpdNo).at(5)->push_back( trackID );
fRPDintVec->at(rpdNo).at(6)->push_back( pid );
} else{
//doubles
fRPDdblVec->at(rpdNo).at(0)->push_back( position.x() );
fRPDdblVec->at(rpdNo).at(1)->push_back( position.y() );
fRPDdblVec->at(rpdNo).at(2)->push_back( position.z() );
//ints
fRPDintVec->at(rpdNo).at(0)->push_back( rodNo );
fRPDintVec->at(rpdNo).at(1)->push_back( TotalnCherenkovs );
}// end if !CLUSTER
// end if RPD
} else{
if( name.compare(0,3,"ZDC") ){//ZDC hitsCollID, check to be sure/symmetric
int zdcNo = atoi( name.substr(3,1).c_str() );
if(!CLUSTER){
//doubles
fZDCdblVec->at(zdcNo).at(0)-> push_back( position.x() );
fZDCdblVec->at(zdcNo).at(1)-> push_back( position.y() );
fZDCdblVec->at(zdcNo).at(2)-> push_back( position.z() );
fZDCdblVec->at(zdcNo).at(3)-> push_back( momentum.x() );
fZDCdblVec->at(zdcNo).at(4)-> push_back( momentum.y() );
fZDCdblVec->at(zdcNo).at(5)-> push_back( momentum.z() );
fZDCdblVec->at(zdcNo).at(6)-> push_back( energy );
fZDCdblVec->at(zdcNo).at(7)-> push_back( velocity );
fZDCdblVec->at(zdcNo).at(8)-> push_back( beta );
fZDCdblVec->at(zdcNo).at(9)-> push_back( TotaleDep );
fZDCdblVec->at(zdcNo).at(10)->push_back( charge );
//ints
fZDCintVec->at(zdcNo).at(0)->push_back( fEventNo );
fZDCintVec->at(zdcNo).at(1)->push_back( modNb );
fZDCintVec->at(zdcNo).at(2)->push_back( radiatorNo );
fZDCintVec->at(zdcNo).at(3)->push_back( rodNo );
fZDCintVec->at(zdcNo).at(4)->push_back( TotalnCherenkovs );
fZDCintVec->at(zdcNo).at(5)->push_back( trackID );
fZDCintVec->at(zdcNo).at(6)->push_back( pid );
} else{
//ints
fZDCintVec->at(zdcNo).at(0)->push_back( radiatorNo );
fZDCintVec->at(zdcNo).at(1)->push_back( TotalnCherenkovs );
}// end if !CLUSTER
}// end if ZDC
}// end else (!RPD)
TotalnCherenkovs = 0;
TotaleDep = 0.0;
hitsCollID++;
}//end of hit loop
// fill ntuples //
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
for(int i = 0; i < analysisManager->GetNofNtuples(); i++){
analysisManager->FillNtupleDColumn(i,1,LastStepInVolume);
analysisManager->FillNtupleDColumn(i,2, pVert->GetX0() );
analysisManager->FillNtupleDColumn(i,3, pVert->GetY0() );
analysisManager->FillNtupleDColumn(i,4, pVert->GetZ0() );
analysisManager->AddNtupleRow(i);
}
// Clear RPD vectors //
for(uint i = 0; i < fRPDdblVec->size(); i++){
//Clear double vectors
for(uint j = 0; j < fRPDdblVec->at(0).size(); j++){
fRPDdblVec->at(i).at(j)->clear();
}
//Clear int vectors
for(uint j = 0; j < fRPDintVec->at(0).size(); j++){
fRPDintVec->at(i).at(j)->clear();
}
}
// Clear ZDC vectors //
for(uint i = 0; i < fZDCdblVec->size(); i++){
//Clear double vectors
for(uint j = 0; j < fZDCdblVec->at(0).size(); j++){
fZDCdblVec->at(i).at(j)->clear();
}
//Clear int vectors
for(uint j = 0; j < fZDCintVec->at(0).size(); j++){
fZDCintVec->at(i).at(j)->clear();
}
}
FiberSD* sd = (FiberSD*)G4SDManager::GetSDMpointer()->FindSensitiveDetector( name );
if( sd->OpticalIsOn() ) ProcessOpticalHitCollection( HC );
else ProcessHitCollection( HC );
hitsCollID++;
}// end while < nCollections
// fill ntuples //
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
for(int i = 0; i < analysisManager->GetNofNtuples(); i++){
analysisManager->FillNtupleDColumn(i,1,LastStepInVolume);
analysisManager->FillNtupleDColumn(i,2, pVert->GetX0() );
analysisManager->FillNtupleDColumn(i,3, pVert->GetY0() );
analysisManager->FillNtupleDColumn(i,4, pVert->GetZ0() );
analysisManager->FillNtupleIColumn(i,0, fEventNo );
analysisManager->AddNtupleRow(i);
}
// Clear ZDC vectors //
for(uint i = 0; i < fZDCdblVec->size(); i++){
fZDCdblVec->at(i).clear();
fZDCintVec->at(i).clear();
}
// Clear RPD vectors //
for(uint i = 0; i < fRPDdblVec->size(); i++){
fRPDdblVec->at(i).clear();
fRPDintVec->at(i).clear();
}
}// end if HCE
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventAction::ProcessHitCollection( FiberHitsCollection* HC ){
G4int prevTrackId = 0;
G4int prevRadiatorNo = 0;
G4int nCherenkovsSum = 0;
G4double eDepSum = 0.0;
G4String name = HC->GetSDname();
int n_hit = HC->entries();
G4cout << "nHits = " << n_hit << G4endl;
for ( G4int i = 0 ; i < n_hit; i++){
G4ThreeVector position = (*HC)[i]->getPos();
G4ThreeVector origin = (*HC)[i]->getOrigin();
G4ThreeVector momentum = (*HC)[i]->getMomentum();
G4double energy = (*HC)[i]->getEnergy();
G4double velocity = (*HC)[i]->getVelocity();
G4double beta = (*HC)[i]->getBeta();
G4double eDep = (*HC)[i]->getEdep();
G4double charge = (*HC)[i]->getCharge();
G4int radiatorNo = (*HC)[i]->getRadNb();
G4int rodNo = (*HC)[i]->getRodNb();
G4int modNb = (*HC)[i]->getModNb();
G4int trackID = (*HC)[i]->getTrackID();
G4int pid = (*HC)[i]->getParticle()->GetPDGEncoding();
G4int nCherenkovs = (*HC)[i]->getNCherenkovs(); // This is the number of cherenkovs in a single step within the SD
//Sum energy from all steps a particle takes in a single scoring volume
if (trackID == prevTrackId || radiatorNo == prevRadiatorNo || eDep != 0) {
nCherenkovsSum += nCherenkovs;
eDepSum += eDep;
prevTrackId = trackID;
prevRadiatorNo = radiatorNo;
continue;
}// end summation
if(name.compare(0,3,"RPD")){ //RPD hits
int rpdNo = atoi( name.substr(3,1).c_str() );
if(!CLUSTER){
//doubles
fRPDdblVec->at(rpdNo).at(0). push_back( position.x() );
fRPDdblVec->at(rpdNo).at(1). push_back( position.y() );
fRPDdblVec->at(rpdNo).at(2). push_back( position.z() );
fRPDdblVec->at(rpdNo).at(3). push_back( momentum.x() );
fRPDdblVec->at(rpdNo).at(4). push_back( momentum.y() );
fRPDdblVec->at(rpdNo).at(5). push_back( momentum.z() );
fRPDdblVec->at(rpdNo).at(6). push_back( energy );
fRPDdblVec->at(rpdNo).at(7). push_back( velocity );
fRPDdblVec->at(rpdNo).at(8). push_back( beta );
fRPDdblVec->at(rpdNo).at(9). push_back( eDepSum );
fRPDdblVec->at(rpdNo).at(10).push_back( charge );
//ints
fRPDintVec->at(rpdNo).at(0).push_back( modNb );
fRPDintVec->at(rpdNo).at(1).push_back( radiatorNo );
fRPDintVec->at(rpdNo).at(2).push_back( rodNo );
fRPDintVec->at(rpdNo).at(3).push_back( nCherenkovsSum );
fRPDintVec->at(rpdNo).at(4).push_back( trackID );
fRPDintVec->at(rpdNo).at(5).push_back( pid );
} else{
//doubles
fRPDdblVec->at(rpdNo).at(0).push_back( position.x() );
fRPDdblVec->at(rpdNo).at(1).push_back( position.y() );
fRPDdblVec->at(rpdNo).at(2).push_back( position.z() );
//ints
fRPDintVec->at(rpdNo).at(0).push_back( nCherenkovsSum );
fRPDintVec->at(rpdNo).at(1).push_back( rodNo );
}// end if !CLUSTER
// end if RPD
} else{
if( name.compare(0,3,"ZDC") ){//ZDC hitsCollID, check to be sure/symmetric
int zdcNo = atoi( name.substr(3,1).c_str() );
if(!CLUSTER){
//doubles
fZDCdblVec->at(zdcNo).at(0). push_back( position.x() );
fZDCdblVec->at(zdcNo).at(1). push_back( position.y() );
fZDCdblVec->at(zdcNo).at(2). push_back( position.z() );
fZDCdblVec->at(zdcNo).at(3). push_back( momentum.x() );
fZDCdblVec->at(zdcNo).at(4). push_back( momentum.y() );
fZDCdblVec->at(zdcNo).at(5). push_back( momentum.z() );
fZDCdblVec->at(zdcNo).at(6). push_back( energy );
fZDCdblVec->at(zdcNo).at(7). push_back( velocity );
fZDCdblVec->at(zdcNo).at(8). push_back( beta );
fZDCdblVec->at(zdcNo).at(9). push_back( eDepSum );
fZDCdblVec->at(zdcNo).at(10).push_back( charge );
//ints
fZDCintVec->at(zdcNo).at(0).push_back( modNb );
fZDCintVec->at(zdcNo).at(1).push_back( radiatorNo );
fZDCintVec->at(zdcNo).at(2).push_back( rodNo );
fZDCintVec->at(zdcNo).at(3).push_back( nCherenkovsSum );
fZDCintVec->at(zdcNo).at(4).push_back( trackID );
fZDCintVec->at(zdcNo).at(5).push_back( pid );
} else{
//doubles
fZDCdblVec->at(zdcNo).at(0). push_back( position.x() );
fZDCdblVec->at(zdcNo).at(1). push_back( position.y() );
fZDCdblVec->at(zdcNo).at(2). push_back( position.z() );
//ints
fZDCintVec->at(zdcNo).at(0).push_back( nCherenkovsSum );
fZDCintVec->at(zdcNo).at(1).push_back( radiatorNo );
}// end if !CLUSTER
}// end if ZDC
}// end else (!RPD)
nCherenkovsSum = 0;
eDepSum = 0.0;
}//end of hit loop
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventAction::ProcessOpticalHitCollection ( FiberHitsCollection* HC ){
G4int prevTrackId = 0;
G4int nCherenkovsSum = 0;
G4String name = HC->GetSDname();
FiberSD* sd = (FiberSD*)G4SDManager::GetSDMpointer()->FindSensitiveDetector( name );
G4double topOfVolume = sd->GetTopOfVolume();
int n_hit = HC->entries();
G4cout << name << " nHits = " << n_hit << G4endl;
for ( G4int i = 0 ; i < n_hit; i++){
G4int nCherenkovs = (*HC)[i]->getNCherenkovs();
//Grab nCherenkovs and skip anything that isn't a photon
if( (*HC)[i]->getParticle() != G4OpticalPhoton::OpticalPhotonDefinition() ){
nCherenkovsSum += nCherenkovs;
continue;
}
//Skip any photons that we've already recorded
G4ThreeVector position = (*HC)[i]->getPos();
G4int trackID = (*HC)[i]->getTrackID();
if( position.y() < topOfVolume - 0.2*mm || trackID == prevTrackId ){
continue;
}
trackID = prevTrackId;
G4ThreeVector origin = (*HC)[i]->getOrigin();
G4ThreeVector momentum = (*HC)[i]->getMomentum();
G4int rodNo = (*HC)[i]->getRodNb();
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
if( trackID != prevTrackId ){
if( name.compare(0,3,"ZDC") ){//ZDC hitsCollID, check to be sure/symmetric
int zdcNo = atoi( name.substr(3,1).c_str() );
fZDCdblVec->at(zdcNo).at(0). push_back( origin.x() );
fZDCdblVec->at(zdcNo).at(1). push_back( origin.y() );
fZDCdblVec->at(zdcNo).at(2). push_back( origin.z() );
fZDCdblVec->at(zdcNo).at(3). push_back( momentum.x() );
fZDCdblVec->at(zdcNo).at(4). push_back( momentum.y() );
fZDCdblVec->at(zdcNo).at(5). push_back( momentum.z() );
analysisManager->FillNtupleIColumn( zdcNo, 1, nCherenkovsSum );
fZDCintVec->at(zdcNo).at(0).push_back( rodNo );
}//end fill ZDC vectors
if( name.compare(0,3,"RPD") ){//RPD hitsCollID, check to be sure/symmetric
int rpdNo = atoi( name.substr(3,1).c_str() );
fRPDdblVec->at(rpdNo).at(0). push_back( origin.x() );
fRPDdblVec->at(rpdNo).at(1). push_back( origin.y() );
fRPDdblVec->at(rpdNo).at(2). push_back( origin.z() );
fRPDdblVec->at(rpdNo).at(3). push_back( momentum.x() );
fRPDdblVec->at(rpdNo).at(4). push_back( momentum.y() );
fRPDdblVec->at(rpdNo).at(5). push_back( momentum.z() );
analysisManager->FillNtupleIColumn( rpdNo + fZDCdblVec->size(), 1, nCherenkovsSum );
fRPDintVec->at(rpdNo).at(0). push_back( rodNo );
}//end fill RPD vectors
}// end if trackID
}// end hit loop
}
This diff is collapsed.
Supports Markdown
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