Commit c12a18cf authored by Chad Lantz's avatar Chad Lantz
Browse files

Fixed string checking and index problems

parent 541f66a7
......@@ -23,10 +23,11 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: RunAction.hh 93886 2015-11-03 08:28:26Z gcosmo $
//
/// \file RunAction.hh
/// \brief Definition of the RunAction class
/// \file AnalysisManager.hh
/// \brief Definition of the AnalysisManager class
/// \author Chad Lantz
/// \date 16 April 2020
#include "AnalysisManager.hh"
#include "DetectorConstruction.hh"
......@@ -51,13 +52,6 @@ AnalysisManager* AnalysisManager::getInstance(void)
AnalysisManager::AnalysisManager()
: m_FactoryOn(false)
{
// Get a pointer to DetectorConstruction
// const DetectorConstruction* constDetectorConstruction
// = static_cast<const DetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
// m_detectorConstruction
// = const_cast<DetectorConstruction*>(constDetectorConstruction);
m_analysisManager = G4AnalysisManager::Instance();
}
......@@ -96,38 +90,48 @@ void AnalysisManager::Book()
return;
}
//Get information about the detector configuration from DetectorConstruction
DetectorConstruction* detectorConstruction = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction();
int nZDCs = detectorConstruction->GetnZDCs();
int nRPDs = detectorConstruction->GetnRPDs();
CLUSTER = detectorConstruction->GetClusterFlag();
//Make vectors for the detectors we have
//Indecies are [module#][dataType][dataPoint]
m_ZDCdblVec = new std::vector< std::vector< std::vector<double> > >(nZDCs);
m_ZDCintVec = new std::vector< std::vector< std::vector< int > > >(nZDCs);
m_RPDdblVec = new std::vector< std::vector< std::vector<double> > >(nRPDs);
m_RPDintVec = new std::vector< std::vector< std::vector< int > > >(nRPDs);
//Create ZDC trees and branches
G4cout << "Number of ZDCs " << nZDCs << G4endl;
G4cout << "Number of RPDs " << nRPDs << G4endl;
FiberSD* sd;
char name[20];
for(int zdcNo = 0; zdcNo < nZDCs; zdcNo++){
sprintf(name,"ZDC%d_SD",zdcNo+1);
//Find out from the SD if optical is on for this detector
sprintf(name,"ZDC%d_SD",zdcNo+1);
sd = (FiberSD*)G4SDManager::GetSDMpointer()->FindSensitiveDetector( name );
if( sd->OpticalIsOn() ) MakeZDCOpticalTree( zdcNo, zdcNo );
else MakeZDCTree( zdcNo, zdcNo );
else MakeZDCTree( zdcNo, zdcNo );
}//end ZDC loop
//Create RPD trees and branches
for(int rpdNo = 0; rpdNo < nRPDs; rpdNo++){
sprintf(name,"RPD%d_SD",rpdNo+1);
int nTuple = nZDCs + rpdNo;
//Find out from the SD if optical is on for this detector
sprintf(name,"RPD%d_SD",rpdNo+1);
sd = (FiberSD*)G4SDManager::GetSDMpointer()->FindSensitiveDetector( name );
if( sd->OpticalIsOn() ) MakeRPDOpticalTree( nTuple, rpdNo );
else MakeRPDTree( nTuple, rpdNo );
else MakeRPDTree( nTuple, rpdNo );
}//end RPD loop
//Turn on the data factory I guess
m_FactoryOn = true;
G4cout << "\n----> Output file is open in "
......@@ -156,21 +160,26 @@ void AnalysisManager::FillNtuples(){
// fill ntuples //
for(int i = 0; i < m_analysisManager->GetNofNtuples(); i++){
G4cout << "Writing nTuple " << i << G4endl;
m_analysisManager->AddNtupleRow(i);
}
// Clear ZDC vectors //
for(uint i = 0; i < m_ZDCdblVec->size(); i++){
std::cout << "Clearing ZDC" << i << std::endl;
m_ZDCdblVec->at(i).clear();
m_ZDCintVec->at(i).clear();
for(uint j = 0; j < m_ZDCdblVec->at(i).size(); j++){
m_ZDCdblVec->at(i).at(j).clear();
}
for(uint j = 0; j < m_ZDCintVec->at(i).size(); j++){
m_ZDCintVec->at(i).at(j).clear();
}
}
// Clear RPD vectors //
for(uint i = 0; i < m_RPDdblVec->size(); i++){
std::cout << "Clearing RPD" << i << std::endl;
m_RPDdblVec->at(i).clear();
m_RPDintVec->at(i).clear();
for(uint j = 0; j < m_RPDdblVec->at(i).size(); j++){
m_RPDdblVec->at(i).at(j).clear();
}
for(uint j = 0; j < m_RPDintVec->at(i).size(); j++){
m_RPDintVec->at(i).at(j).clear();
}
}
}
......@@ -181,12 +190,10 @@ void AnalysisManager::MakeZDCTree( G4int nTupleNo, G4int zdcNo ){
sprintf(name,"ZDC%dtree",zdcNo+1);
m_analysisManager->CreateNtuple( name, "ZDC data");
if(!CLUSTER){
G4cout << "making non cluster tree" << G4endl;
//Resize the vector for the number of branches storing double vectors
m_ZDCdblVec->at(zdcNo).resize(11);
//Make double branches
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStep" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastSteptest" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosX" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosY" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosZ" );
......@@ -214,12 +221,10 @@ void AnalysisManager::MakeZDCTree( G4int nTupleNo, G4int zdcNo ){
m_analysisManager->CreateNtupleIColumn( nTupleNo, "pid", m_ZDCintVec->at(zdcNo).at(5) );
} else { //Fewer vector branches to save storage space
std::cout << "making cluster tree" << std::endl;
//Resize the vector for the number of branches storing double vectors
m_ZDCdblVec->at(zdcNo).resize(6);
//Make double branches
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStep" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStepTest" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosX" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosY" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosZ" );
......@@ -243,13 +248,11 @@ void AnalysisManager::MakeZDCOpticalTree( G4int nTupleNo, G4int zdcNo ){
char name[20];
sprintf(name,"ZDC%dtree",zdcNo+1);
m_analysisManager->CreateNtuple( name, "ZDC data");
G4cout << "making optical tree" << G4endl;
//Resize the vector for the number of branches storing double vectors
m_ZDCdblVec->at(zdcNo).resize(6);
//Make double branches
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStep" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStepTest" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosX" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosY" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosZ" );
......@@ -281,7 +284,6 @@ void AnalysisManager::MakeRPDTree( G4int nTupleNo, G4int rpdNo ){
m_RPDdblVec->at(rpdNo).resize(11);
//Make double branches
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStep" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStepTest" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosX" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosY" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosZ" );
......@@ -313,7 +315,6 @@ void AnalysisManager::MakeRPDTree( G4int nTupleNo, G4int rpdNo ){
m_RPDdblVec->at(rpdNo).resize(3);
//Make double branches
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStep" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "lastStepTest" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosX" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosY" );
m_analysisManager->CreateNtupleDColumn( nTupleNo, "gunPosZ" );
......
......@@ -23,10 +23,10 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// Author: Chad Lantz
//
// \file RunAction.cc
// \brief Implementation of the RunAction class
/// \file RunAction.cc
/// \brief Implementation of the RunAction class
/// \author Chad Lantz
/// \date 16 April 2020
#include "EventAction.hh"
#include "FiberHit.hh"
......@@ -87,7 +87,7 @@ void EventAction::EndOfEventAction(const G4Event* evt){
FiberSD* sd = (FiberSD*)G4SDManager::GetSDMpointer()->FindSensitiveDetector( name );
if( sd->OpticalIsOn() ) ProcessOpticalHitCollection( HC );
else ProcessHitCollection( HC );
else ProcessHitCollection( HC );
hitsCollID++;
}// end while < nCollections
......@@ -98,9 +98,9 @@ void EventAction::EndOfEventAction(const G4Event* evt){
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
// fill ntuples //
for(int i = 0; i < analysisManager->GetNofNtuples(); i++){
analysisManager->FillNtupleDColumn(i,2, pVert->GetX0() );
analysisManager->FillNtupleDColumn(i,3, pVert->GetY0() );
analysisManager->FillNtupleDColumn(i,4, pVert->GetZ0() );
analysisManager->FillNtupleDColumn(i,1, pVert->GetX0() );
analysisManager->FillNtupleDColumn(i,2, pVert->GetY0() );
analysisManager->FillNtupleDColumn(i,3, pVert->GetZ0() );
analysisManager->FillNtupleIColumn(i,0, fEventNo );
}
......@@ -121,24 +121,14 @@ void EventAction::ProcessHitCollection( FiberHitsCollection* HC ){
int n_hit = HC->entries();
G4cout << "nHits = " << n_hit << G4endl;
G4cout << name << " 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
G4int trackID = (*HC)[i]->getTrackID();
//Sum energy from all steps a particle takes in a single scoring volume
if (trackID == prevTrackId || radiatorNo == prevRadiatorNo || eDep != 0) {
......@@ -150,76 +140,88 @@ void EventAction::ProcessHitCollection( FiberHitsCollection* HC ){
continue;
}// end summation
if(name.compare(0,3,"RPD")){ //RPD hits
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 charge = (*HC)[i]->getCharge();
G4int rodNo = (*HC)[i]->getRodNb();
G4int modNb = (*HC)[i]->getModNb();
G4int pid = (*HC)[i]->getParticle()->GetPDGEncoding();
if(name.compare(0,3,"RPD") == 0){ //RPD hits
int rpdNo = atoi( name.substr(3,1).c_str() );
if(!CLUSTER){
//doubles
m_RPDdblVec->at(rpdNo).at(0). push_back( position.x() );
m_RPDdblVec->at(rpdNo).at(1). push_back( position.y() );
m_RPDdblVec->at(rpdNo).at(2). push_back( position.z() );
m_RPDdblVec->at(rpdNo).at(3). push_back( momentum.x() );
m_RPDdblVec->at(rpdNo).at(4). push_back( momentum.y() );
m_RPDdblVec->at(rpdNo).at(5). push_back( momentum.z() );
m_RPDdblVec->at(rpdNo).at(6). push_back( energy );
m_RPDdblVec->at(rpdNo).at(7). push_back( velocity );
m_RPDdblVec->at(rpdNo).at(8). push_back( beta );
m_RPDdblVec->at(rpdNo).at(9). push_back( eDepSum );
m_RPDdblVec->at(rpdNo).at(10).push_back( charge );
m_RPDdblVec->at(rpdNo-1).at(0). push_back( position.x() );
m_RPDdblVec->at(rpdNo-1).at(1). push_back( position.y() );
m_RPDdblVec->at(rpdNo-1).at(2). push_back( position.z() );
m_RPDdblVec->at(rpdNo-1).at(3). push_back( momentum.x() );
m_RPDdblVec->at(rpdNo-1).at(4). push_back( momentum.y() );
m_RPDdblVec->at(rpdNo-1).at(5). push_back( momentum.z() );
m_RPDdblVec->at(rpdNo-1).at(6). push_back( energy );
m_RPDdblVec->at(rpdNo-1).at(7). push_back( velocity );
m_RPDdblVec->at(rpdNo-1).at(8). push_back( beta );
m_RPDdblVec->at(rpdNo-1).at(9). push_back( eDepSum );
m_RPDdblVec->at(rpdNo-1).at(10).push_back( charge );
//ints
m_RPDintVec->at(rpdNo).at(0).push_back( modNb );
m_RPDintVec->at(rpdNo).at(1).push_back( radiatorNo );
m_RPDintVec->at(rpdNo).at(2).push_back( rodNo );
m_RPDintVec->at(rpdNo).at(3).push_back( nCherenkovsSum );
m_RPDintVec->at(rpdNo).at(4).push_back( trackID );
m_RPDintVec->at(rpdNo).at(5).push_back( pid );
m_RPDintVec->at(rpdNo-1).at(0).push_back( modNb );
m_RPDintVec->at(rpdNo-1).at(1).push_back( radiatorNo );
m_RPDintVec->at(rpdNo-1).at(2).push_back( rodNo );
m_RPDintVec->at(rpdNo-1).at(3).push_back( nCherenkovsSum );
m_RPDintVec->at(rpdNo-1).at(4).push_back( trackID );
m_RPDintVec->at(rpdNo-1).at(5).push_back( pid );
} else{
//doubles
m_RPDdblVec->at(rpdNo).at(0).push_back( position.x() );
m_RPDdblVec->at(rpdNo).at(1).push_back( position.y() );
m_RPDdblVec->at(rpdNo).at(2).push_back( position.z() );
m_RPDdblVec->at(rpdNo-1).at(0).push_back( position.x() );
m_RPDdblVec->at(rpdNo-1).at(1).push_back( position.y() );
m_RPDdblVec->at(rpdNo-1).at(2).push_back( position.z() );
//ints
m_RPDintVec->at(rpdNo).at(0).push_back( nCherenkovsSum );
m_RPDintVec->at(rpdNo).at(1).push_back( rodNo );
m_RPDintVec->at(rpdNo-1).at(0).push_back( nCherenkovsSum );
m_RPDintVec->at(rpdNo-1).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
if( name.compare(0,3,"ZDC") == 0 ){//ZDC hitsCollID, check to be sure/symmetric
int zdcNo = atoi( name.substr(3,1).c_str() );
G4cout << "Processing ZDC" << zdcNo << " hits" << G4endl;
if(!CLUSTER){
//doubles
m_ZDCdblVec->at(zdcNo).at(0). push_back( position.x() );
m_ZDCdblVec->at(zdcNo).at(1). push_back( position.y() );
m_ZDCdblVec->at(zdcNo).at(2). push_back( position.z() );
m_ZDCdblVec->at(zdcNo).at(3). push_back( momentum.x() );
m_ZDCdblVec->at(zdcNo).at(4). push_back( momentum.y() );
m_ZDCdblVec->at(zdcNo).at(5). push_back( momentum.z() );
m_ZDCdblVec->at(zdcNo).at(6). push_back( energy );
m_ZDCdblVec->at(zdcNo).at(7). push_back( velocity );
m_ZDCdblVec->at(zdcNo).at(8). push_back( beta );
m_ZDCdblVec->at(zdcNo).at(9). push_back( eDepSum );
m_ZDCdblVec->at(zdcNo).at(10).push_back( charge );
m_ZDCdblVec->at(zdcNo-1).at(0). push_back( position.x() );
m_ZDCdblVec->at(zdcNo-1).at(1). push_back( position.y() );
m_ZDCdblVec->at(zdcNo-1).at(2). push_back( position.z() );
m_ZDCdblVec->at(zdcNo-1).at(3). push_back( momentum.x() );
m_ZDCdblVec->at(zdcNo-1).at(4). push_back( momentum.y() );
m_ZDCdblVec->at(zdcNo-1).at(5). push_back( momentum.z() );
m_ZDCdblVec->at(zdcNo-1).at(6). push_back( energy );
m_ZDCdblVec->at(zdcNo-1).at(7). push_back( velocity );
m_ZDCdblVec->at(zdcNo-1).at(8). push_back( beta );
m_ZDCdblVec->at(zdcNo-1).at(9). push_back( eDepSum );
m_ZDCdblVec->at(zdcNo-1).at(10).push_back( charge );
//ints
m_ZDCintVec->at(zdcNo).at(0).push_back( modNb );
m_ZDCintVec->at(zdcNo).at(1).push_back( radiatorNo );
m_ZDCintVec->at(zdcNo).at(2).push_back( rodNo );
m_ZDCintVec->at(zdcNo).at(3).push_back( nCherenkovsSum );
m_ZDCintVec->at(zdcNo).at(4).push_back( trackID );
m_ZDCintVec->at(zdcNo).at(5).push_back( pid );
m_ZDCintVec->at(zdcNo-1).at(0).push_back( modNb );
m_ZDCintVec->at(zdcNo-1).at(1).push_back( radiatorNo );
m_ZDCintVec->at(zdcNo-1).at(2).push_back( rodNo );
m_ZDCintVec->at(zdcNo-1).at(3).push_back( nCherenkovsSum );
m_ZDCintVec->at(zdcNo-1).at(4).push_back( trackID );
m_ZDCintVec->at(zdcNo-1).at(5).push_back( pid );
} else{
//doubles
m_ZDCdblVec->at(zdcNo).at(0). push_back( position.x() );
m_ZDCdblVec->at(zdcNo).at(1). push_back( position.y() );
m_ZDCdblVec->at(zdcNo).at(2). push_back( position.z() );
m_ZDCdblVec->at(zdcNo-1).at(0). push_back( position.x() );
m_ZDCdblVec->at(zdcNo-1).at(1). push_back( position.y() );
m_ZDCdblVec->at(zdcNo-1).at(2). push_back( position.z() );
//ints
m_ZDCintVec->at(zdcNo).at(0).push_back( nCherenkovsSum );
m_ZDCintVec->at(zdcNo).at(1).push_back( radiatorNo );
m_ZDCintVec->at(zdcNo-1).at(0).push_back( nCherenkovsSum );
m_ZDCintVec->at(zdcNo-1).at(1).push_back( radiatorNo );
}// end if !CLUSTER
}// end if ZDC
}// end else (!RPD)
......@@ -264,29 +266,29 @@ void EventAction::ProcessOpticalHitCollection ( FiberHitsCollection* HC ){
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
if( trackID != prevTrackId ){
if( name.compare(0,3,"ZDC") ){//ZDC hitsCollID, check to be sure/symmetric
if( name.compare(0,3,"ZDC") == 0 ){//ZDC hitsCollID, check to be sure/symmetric
int zdcNo = atoi( name.substr(3,1).c_str() );
m_ZDCdblVec->at(zdcNo).at(0). push_back( origin.x() );
m_ZDCdblVec->at(zdcNo).at(1). push_back( origin.y() );
m_ZDCdblVec->at(zdcNo).at(2). push_back( origin.z() );
m_ZDCdblVec->at(zdcNo).at(3). push_back( momentum.x() );
m_ZDCdblVec->at(zdcNo).at(4). push_back( momentum.y() );
m_ZDCdblVec->at(zdcNo).at(5). push_back( momentum.z() );
m_ZDCdblVec->at(zdcNo-1).at(0). push_back( origin.x() );
m_ZDCdblVec->at(zdcNo-1).at(1). push_back( origin.y() );
m_ZDCdblVec->at(zdcNo-1).at(2). push_back( origin.z() );
m_ZDCdblVec->at(zdcNo-1).at(3). push_back( momentum.x() );
m_ZDCdblVec->at(zdcNo-1).at(4). push_back( momentum.y() );
m_ZDCdblVec->at(zdcNo-1).at(5). push_back( momentum.z() );
analysisManager->FillNtupleIColumn( zdcNo, 1, nCherenkovsSum );
m_ZDCintVec->at(zdcNo).at(0).push_back( rodNo );
m_ZDCintVec->at(zdcNo-1).at(0).push_back( rodNo );
}//end fill ZDC vectors
if( name.compare(0,3,"RPD") ){//RPD hitsCollID, check to be sure/symmetric
if( name.compare(0,3,"RPD") == 0 ){//RPD hitsCollID, check to be sure/symmetric
int rpdNo = atoi( name.substr(3,1).c_str() );
m_RPDdblVec->at(rpdNo).at(0). push_back( origin.x() );
m_RPDdblVec->at(rpdNo).at(1). push_back( origin.y() );
m_RPDdblVec->at(rpdNo).at(2). push_back( origin.z() );
m_RPDdblVec->at(rpdNo).at(3). push_back( momentum.x() );
m_RPDdblVec->at(rpdNo).at(4). push_back( momentum.y() );
m_RPDdblVec->at(rpdNo).at(5). push_back( momentum.z() );
m_RPDdblVec->at(rpdNo-1).at(0). push_back( origin.x() );
m_RPDdblVec->at(rpdNo-1).at(1). push_back( origin.y() );
m_RPDdblVec->at(rpdNo-1).at(2). push_back( origin.z() );
m_RPDdblVec->at(rpdNo-1).at(3). push_back( momentum.x() );
m_RPDdblVec->at(rpdNo-1).at(4). push_back( momentum.y() );
m_RPDdblVec->at(rpdNo-1).at(5). push_back( momentum.z() );
analysisManager->FillNtupleIColumn( rpdNo + m_ZDCdblVec->size(), 1, nCherenkovsSum );
m_RPDintVec->at(rpdNo).at(0). push_back( rodNo );
m_RPDintVec->at(rpdNo-1).at(0). push_back( rodNo );
}//end fill RPD vectors
}// end if trackID
}// end hit loop
......
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