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

Working toward getting all information and killing optical photons in FiberSD

parent e65a1c8e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 191 additions and 298 deletions
MonteCarlo/include/EventAction.hh
+ 6
− 3
View file @ b9b9dce3
......@@ -34,13 +34,15 @@
#include "G4UserEventAction.hh"
#include "globals.hh"
#include <vector>
/// Event action class
///
class EventAction : public G4UserEventAction
{
public:
EventAction(RunAction* runAction);
EventAction( );
virtual ~EventAction();
virtual void BeginOfEventAction (const G4Event* event);
......@@ -55,11 +57,12 @@ class EventAction : public G4UserEventAction
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;
......
MonteCarlo/include/FiberHit.hh
+ 42
− 38
View file @ b9b9dce3
......@@ -58,45 +58,49 @@ public:
void Print();
public:
void setTrackID (G4int track) { trackID = track; }
void setModNb (G4int mod) { modNb = mod; }
void setRadNb (G4int rad) { radNb = rad; }
void setRodNb (G4int rod) { rodNb = rod; }
void setEdep (G4double de) { edep = de; }
void setPos (G4ThreeVector xyz) { pos = xyz; }
void setMomentum (G4ThreeVector mom) { momentum = mom;}
void setParticle (G4ParticleDefinition *part) {particle = part;}
void setEnergy (G4double e) {energy = e;}
void setCharge (G4double c) {charge = c;}
void setNCherenkovs (G4int n) {nCherenkovs = n;}
G4int getTrackID() { return trackID; }
G4int getModNb() { return modNb; }
G4int getRadNb() { return radNb; }
G4int getRodNb() { return rodNb; }
G4int getNCherenkovs() { return nCherenkovs; }
G4double getEdep() { return edep; }
G4ThreeVector getPos() { return pos; }
G4ParticleDefinition* getParticle() { return particle; }
G4double getEnergy() { return energy; }
G4ThreeVector getMomentum() { return momentum; }
G4double getCharge() {return charge;}
G4double getVelocity() {return velocity;}
G4double getBeta() {return beta;}
void setParticle (G4ParticleDefinition *part) { particle = part; }
void setPos (G4ThreeVector xyz) { pos = xyz; }
void setMomentum (G4ThreeVector mom) { momentum = mom; }
void setTrackID (G4int track) { trackID = track; }
void setModNb (G4int mod) { modNb = mod; }
void setRadNb (G4int rad) { radNb = rad; }
void setRodNb (G4int rod) { rodNb = rod; }
void setNCherenkovs (G4int n) { nCherenkovs = n; }
void setEdep (G4double de) { edep = de; }
void setEnergy (G4double e) { energy = e; }
void setCharge (G4double c) { charge = c; }
G4ParticleDefinition* getParticle (){ return particle; }
G4ThreeVector getPos (){ return pos; }
G4ThreeVector getMomentum (){ return momentum; }
G4int getTrackID (){ return trackID; }
G4int getModNb (){ return modNb; }
G4int getRadNb (){ return radNb; }
G4int getRodNb (){ return rodNb; }
G4int getNCherenkovs (){ return nCherenkovs; }
G4double getEdep (){ return edep; }
G4double getEnergy (){ return energy; }
G4double getCharge (){ return charge; }
G4double getVelocity (){ return velocity; }
G4double getBeta (){ return beta; }
private:
G4int trackID;
G4int modNb;
G4int radNb;
G4int rodNb;
G4double velocity;
G4double beta;
G4int nCherenkovs;
G4double edep;
G4ThreeVector pos;
G4ParticleDefinition *particle;
G4ThreeVector momentum;
G4double energy;
G4double charge;
G4ParticleDefinition* particle;
G4ThreeVector pos;
G4ThreeVector momentum;
G4int trackID;
G4int modNb;
G4int radNb;
G4int rodNb;
G4int nCherenkovs;
G4double velocity;
G4double beta;
G4double energy;
G4double charge;
G4double edep;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......
MonteCarlo/include/FiberSD.hh
+ 6
− 5
View file @ b9b9dce3
......@@ -32,6 +32,7 @@
#define FiberSD_h 1
#include "G4VSensitiveDetector.hh"
#include "FiberHit.hh"
class G4Step;
class G4HCofThisEvent;
......@@ -41,21 +42,21 @@ class G4HCofThisEvent;
class FiberSD : public G4VSensitiveDetector
{
public:
FiberSD(G4String, G4int);
FiberSD(G4String, G4int,G4bool);
~FiberSD();
void HistInitialize();
void Initialize(G4HCofThisEvent*);
int CalculateCherenkovs(G4Step*,int);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
void Initialize ( G4HCofThisEvent* );
G4bool ProcessHits ( G4Step*, G4TouchableHistory* );
void EndOfEvent ( G4HCofThisEvent* );
private:
int HCID;
G4double m_modCoreIndexRefraction;
FiberHitsCollection* fiberCollection;
G4int m_modNum;
G4bool OPTICAL;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......
MonteCarlo/include/SteppingAction.hh
+ 1
− 0
View file @ b9b9dce3
......@@ -48,6 +48,7 @@ class SteppingAction : public G4UserSteppingAction
private:
G4double lastStep;
G4bool OPTICAL;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......
MonteCarlo/src/EventAction.cc
+ 90
− 85
View file @ b9b9dce3
......@@ -30,6 +30,7 @@
#include "EventAction.hh"
#include "FiberHit.hh"
#include "Analysis.hh"
#include "G4UnitsTable.hh"
......@@ -41,10 +42,8 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventAction::EventAction(RunAction* runAction)
: G4UserEventAction(),
fRunAction(runAction)
{
EventAction::EventAction( )
: G4UserEventAction(){
hitsCollID = -1;
}
......@@ -62,8 +61,8 @@ void EventAction::SetnZDCs(G4int 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) ? 10 : 0;
G4int nIntVecs = (!CLUSTER) ? 8 : 2;
G4int nDblVecs = (!CLUSTER) ? 11 : 0;
G4int nIntVecs = (!CLUSTER) ? 7 : 2;
for(G4int i = 0; i < nZDCs; i++){
......@@ -74,7 +73,7 @@ void EventAction::SetnZDCs(G4int nZDCs){
fZDCintVec->at(i).resize(nIntVecs);
for(G4int j = 0; j < nIntVecs; j++){
fZDCdblVec->at(i).at(j) = new std::vector< int >;
fZDCintVec->at(i).at(j) = new std::vector< int >;
}//end int vector loop
}//end ZDC loop
}
......@@ -88,8 +87,8 @@ void EventAction::SetnRPDs(G4int 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) ? 10 : 3;
G4int nIntVecs = (!CLUSTER) ? 8 : 2;
G4int nDblVecs = (!CLUSTER) ? 11 : 3;
G4int nIntVecs = (!CLUSTER) ? 7 : 2;
for(G4int i = 0; i < nRPDs; i++){
......@@ -100,7 +99,7 @@ void EventAction::SetnRPDs(G4int nRPDs){
fRPDintVec->at(i).resize(nIntVecs);
for(G4int j = 0; j < nIntVecs; j++){
fRPDdblVec->at(i).at(j) = new std::vector< int >;
fRPDintVec->at(i).at(j) = new std::vector< int >;
}//end int vector loop
}//end RPD loop
}
......@@ -121,39 +120,42 @@ void EventAction::EndOfEventAction(const G4Event* evt){
gunPosY = pVert->GetY0();
gunPosZ = pVert->GetZ0();
// Last step in volume?
G4double LastStepInVolume = pVert->GetPosition().z();
G4cout << ">>> Event " << evt->GetEventID() << G4endl;
G4HCofThisEvent * HCE = evt->GetHCofThisEvent();
QuartzHitsCollection* HC = 0;
FiberHitsCollection* HC = 0;
G4int nCollections = HCE->GetNumberOfCollections();
int IDholder = 0;
if(HCE) {
while (hitsCollID < nCollections) {
HC = (FiberHitsCollection*)(HCE->GetHC(hitsCollID));
// Turn this into how we filter what we want to fill
G4String name = HC->GetSDname();
int n_hit = HC->entries();
IDholder = hitsCollID;
int prevTrackId = 0;
int prevRadiatorNo = 0;
G4int prevTrackId = 0;
G4int prevRadiatorNo = 0;
G4int TotalnCherenkovs = 0;
G4double TotaleDep = 0;
std::cout << Form("%s hits collection, nHits = %d" name.c_str(),n_hit) << std::endl;
for ( int i = 0 ; i < n_hit; i++){
G4double TotaleDep = 0.0;
int n_hit = HC->entries();
std::cout << printf("%s hits collection, nHits = %d", name.c_str(),n_hit) << std::endl;
for ( G4int i = 0 ; i < n_hit; i++){
G4int radiatorNo = (*HC)[i]->getRadNb();
G4int rodNo = (*HC)[i]->getRodNb();
G4double eDep = (*HC)[i]->getEdep();
G4int modNb = (*HC)[i]->getModNb();
G4int trackID = (*HC)[i]->getTrackID();
G4ThreeVector position = (*HC)[i]->getPos();
G4ThreeVector momentum = (*HC)[i]->getMomentum();
G4double energy = (*HC)[i]->getEnergy();
G4int pid = (*HC)[i]->getParticle()->GetPDGEncoding();
G4int nCherenkovs = (*HC)[i]->getNCherenkovs();
G4double charge = (*HC)[i]->getCharge();
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) {
......@@ -165,70 +167,72 @@ void EventAction::EndOfEventAction(const G4Event* evt){
continue;
}// end summation
if(name.strcmp("RPD",3)){ //RPD hits
int rpdNo = atoi(name.c_str(),3,1);
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(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 );
fRPDintVec->at(rpdNo).at(7).push_back( charge );
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() );
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 );
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.strcmp("ZDC",3) ){//ZDC hitsCollID, check to be sure/symmetric
int zdcNo = atoi(name.c_str(),3,1);
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(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 );
fZDCintVec->at(zdcNo).at(7).push_back( charge );
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 );
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)
......@@ -239,31 +243,32 @@ void EventAction::EndOfEventAction(const G4Event* evt){
// fill ntuples //
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
for(int i = 0; i < fnZDCs + fnRPDs; i++){
for(int i = 0; i < analysisManager->GetNofNtuples(); i++){
analysisManager->FillNtupleDColumn(i,1,LastStepInVolume);
analysisManager->AddNtupleRow(i);
}
// Clear RPD vectors //
for(int i = 0; i < fRPDdblVec->size(); i++){
for(uint i = 0; i < fRPDdblVec->size(); i++){
//Clear double vectors
for(int j = 0; j < fRPDdblVec->at(0).size(); j++){
fRPDdblVec->at(i).at(j).clear();
for(uint j = 0; j < fRPDdblVec->at(0).size(); j++){
fRPDdblVec->at(i).at(j)->clear();
}
//Clear int vectors
for(int j = 0; j < fRPDintVec->at(0).size(); j++){
fRPDintVec->at(i).at(j).clear();
for(uint j = 0; j < fRPDintVec->at(0).size(); j++){
fRPDintVec->at(i).at(j)->clear();
}
}
// Clear ZDC vectors //
for(int i = 0; i < fZDCdblVec->size(); i++){
for(uint i = 0; i < fZDCdblVec->size(); i++){
//Clear double vectors
for(int j = 0; j < fZDCdblVec->at(0).size(); j++){
fZDCdblVec->at(i).at(j).clear();
for(uint j = 0; j < fZDCdblVec->at(0).size(); j++){
fZDCdblVec->at(i).at(j)->clear();
}
//Clear int vectors
for(int j = 0; j < fZDCintVec->at(0).size(); j++){
fZDCintVec->at(i).at(j).clear();
for(uint j = 0; j < fZDCintVec->at(0).size(); j++){
fZDCintVec->at(i).at(j)->clear();
}
}
......
MonteCarlo/src/FiberSD.cc
+ 26
− 148
View file @ b9b9dce3
......@@ -26,7 +26,7 @@
// Author: Michael Phipps
#include "FiberSD.hh"
#include "SharedData.hh"
#include "SteppingAction.hh"
#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
......@@ -37,16 +37,14 @@
#include "G4ParticleTypes.hh"
#include "G4ParticleDefinition.hh"
#include "TString.h"
#include <string>
#include <iostream>
#include <cmath>
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
FiberSD::FiberSD(G4String sdName, SharedData* sd, G4int modNum)
:G4VSensitiveDetector(sdName), m_sd(sd), m_modNum(modNum) {
FiberSD::FiberSD(G4String sdName, G4int modNum, G4bool optical)
:G4VSensitiveDetector(sdName), m_modNum(modNum), OPTICAL(optical) {
collectionName.insert(sdName);
HCID = -1;
......@@ -71,8 +69,6 @@ void FiberSD::Initialize(G4HCofThisEvent* HCE){
std::string name = collectionName[0];
// static G4int HCID = -1;
if(HCID<0)
{ HCID = G4SDManager::GetSDMpointer()->GetCollectionID( name );}
......@@ -83,56 +79,46 @@ void FiberSD::Initialize(G4HCofThisEvent* HCE){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool FiberSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
G4double eDep = aStep->GetTotalEnergyDeposit();
TEnv* config = m_sd->GetConfig();
G4Track* theTrack = aStep->GetTrack();
if (config->GetValue("OPTICAL_ON",false) == 1){
if(theTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return false;}
G4int modNum = m_modNum;
//G4int modNStripsPerGap;
//modNStripsPerGap = 64;
char variable[256];
sprintf(variable,"mod%dCoreIndexRefraction",6);
m_modCoreIndexRefraction = config->GetValue( variable,1.46);
G4double eDep = aStep->GetTotalEnergyDeposit();
// Figure out if this is necessary
//
G4int totalRodNum = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(0);
G4String radNum_s = "7" + aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() + "0";
G4String radNum_s = "7" + aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() + "0";
G4int rodNum;
G4int radNum;
//radNum = totalRodNum / modNStripsPerGap;
//rodNum = totalRodNum % modNStripsPerGap;
rodNum = totalRodNum;
radNum = (std::stoi (radNum_s)*100)+rodNum;
/*
G4cout << G4endl << "TrackID: " << aStep->GetTrack()->GetTrackID() << G4endl;
G4cout << "RodNum = " << rodNum << G4endl;
G4cout << "physvol " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << G4endl;
G4cout << "logvol " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName() << G4endl;
*/
// ^^^^ Figure out if this is necessary ^^^^
G4double energy = aStep->GetPreStepPoint()->GetTotalEnergy();
G4ThreeVector momentum = aStep->GetPreStepPoint()->GetMomentum();
G4ParticleDefinition *particle = aStep->GetTrack()->GetDefinition();
G4double charge = aStep->GetPreStepPoint()->GetCharge();
//Get the number of Cherenkov photons created in this step
//If optical is off, stop that track
int capturedPhotons =0;
if (config->GetValue("OPTICAL_ON",false) == 0){
capturedPhotons = CalculateCherenkovs(aStep, modNum);}
const std::vector< const G4Track* >* secVec = aStep->GetSecondaryInCurrentStep();
const G4Track* track;
//for(const G4Track* track : *secVec){
for(uint i = 0; i < secVec->size(); i++){
track = secVec->at(i);
if( track->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()){
capturedPhotons++;
if(!OPTICAL){
track->SetTrackStatus( fStopAndKill );
}//end if optical
}//end if pid==0
}//end secondary track loop
FiberHit* newHit = new FiberHit();
newHit->setCharge ( charge );
newHit->setCharge (charge );
newHit->setTrackID (aStep->GetTrack()->GetTrackID() );
newHit->setModNb (modNum );
newHit->setModNb (m_modNum );
newHit->setRadNb (radNum );
newHit->setRodNb (rodNum );
newHit->setEdep (eDep );
......@@ -141,120 +127,12 @@ if (config->GetValue("OPTICAL_ON",false) == 1){
newHit->setEnergy (energy);
newHit->setMomentum (momentum);
newHit->setNCherenkovs (capturedPhotons);
// if (capturedPhotons != 0 ) std::cout << " capturedPhotons " << capturedPhotons << std::endl;
if( particle->GetPDGEncoding() == 0 ) fiberCollection->insert (newHit ); // only want to record photons
// only want to record photons if optical flag is on
if( OPTICAL && particle->GetPDGEncoding() == 0 ) fiberCollection->insert (newHit );
return true;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int FiberSD::CalculateCherenkovs(G4Step* aStep,__attribute__((unused)) int modNum) {
const G4DynamicParticle* aParticle = aStep->GetTrack()->GetDynamicParticle();
const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
// LogStream<<MSG::INFO<<"ZDC stripCharge "<< charge << endreq;
if (charge==0) { return false; }
const G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
const G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
const G4double beta = (pPreStepPoint ->GetBeta() + pPostStepPoint->GetBeta())/2.;
if (beta==0) { return false; }
TEnv* config = m_sd->GetConfig();
G4double minWavelength = config->GetValue("cherenkovMinWavelength",250) * CLHEP::nanometer ;
G4double maxWavelength = config->GetValue("cherenkovMaxWavelength",600) * CLHEP::nanometer ;
G4double inverseWLDiff = (1./minWavelength) - (1./maxWavelength);
const float step_length(aStep->GetStepLength());
G4double MeanNumberOfPhotons = 2*TMath::Pi()*(1./137.)*step_length*inverseWLDiff*(charge)*(charge)*(1.0 - 1.0/(beta*m_modCoreIndexRefraction*beta*m_modCoreIndexRefraction));
if (MeanNumberOfPhotons <= 0.0) { return false; }
// std::cout << " n photons " << MeanNumberOfPhotons << std::endl;
const G4int NumPhotons = (G4int)G4Poisson(MeanNumberOfPhotons);
if (NumPhotons <= 0) { return false; }
const G4ThreeVector pos = pPreStepPoint->GetPosition();
// float yPos = pos.y();
const G4ThreeVector p0 = aStep->GetDeltaPosition().unit();
const float BetaInverse = 1./beta;
float coreIndexRefraction;
float claddingIndexRefraction;
std::string modCladding;
bool cladding;
char name[256];
sprintf(name,"mod%dCoreIndexRefraction",6);
coreIndexRefraction = config->GetValue(name,1.46) ;
sprintf(name,"mod%dCladdingIndexRefraction",6);
claddingIndexRefraction = config->GetValue(name,1.43) ;
sprintf(name,"mod%dCladding",6);
modCladding = config->GetValue(name,"true") ;
cladding = modCladding == "true" ? true : false;
if (!cladding) claddingIndexRefraction = 1.;
float criticalAngle = asin(claddingIndexRefraction/coreIndexRefraction);
// std::cout << " mod " << modNum <<" core n " << coreIndexRefraction << " cladding n " << claddingIndexRefraction << " critical angle " << criticalAngle << std::endl;
int photonCount = 0;
for (G4int I = 0; I < NumPhotons; I++) {
G4double rand;
// float sampledEnergy;
float cosTheta, sin2Theta;
// sample an energy for Photon
rand = G4UniformRand();
// sampledEnergy = Pmin + rand * dp;
cosTheta = BetaInverse / coreIndexRefraction;
sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
// Generate random position of photon on cone surface defined by Theta
rand = G4UniformRand();
const float phi = 2.*M_PI*rand;
const float sinPhi = sin(phi);
const float cosPhi = cos(phi);
// calculate x,y, and z components of photon momentum
// (in coord system with primary particle direction aligned with the z axis)
// and Rotate momentum direction back to global reference system
const float sinTheta = sqrt(sin2Theta);
const float px = sinTheta*cosPhi;
const float py = sinTheta*sinPhi;
const float pz = cosTheta;
G4ParticleMomentum photonMomentum(px, py, pz);
photonMomentum.rotateUz(p0);
bool Transmission=0;
const float PT = sqrt(photonMomentum.getX()*photonMomentum.getX() + photonMomentum.getZ()*photonMomentum.getZ());
const float PY = photonMomentum.getY();
if (PY<=0) Transmission = 1; // If photon is travelling with -ve PY, Its not going to reach the PMT
else if (PT==0) Transmission = 1; // if photon is along Y-axis it will be completely transmitted into the PMT
else {
const float Theta = M_PI/2.0-atan(PT/PY);
if (Theta < criticalAngle) Transmission = 0;
else Transmission=1.0;
}
if (Transmission == 1.0) ++photonCount;
}
// std::cout << " totalPhotonsCaptured " << photonCount << " % " << (float) photonCount/NumPhotons << " critical angle " << criticalAngle << std::endl;
return photonCount;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void FiberSD::EndOfEvent(G4HCofThisEvent*)
......
MonteCarlo/src/RunAction.cc
+ 7
− 1
View file @ b9b9dce3
......@@ -89,6 +89,8 @@ RunAction::RunAction()
analysisManager->CreateNtuple( Form("ZDC%dtree",zdcNo), "ZDC data");
if(!CLUSTER){
//Make double branches
analysisManager->CreateNtupleDColumn( zdcNo, "lastStep" );
analysisManager->CreateNtupleDColumn( zdcNo, "lastSteptest" );
analysisManager->CreateNtupleDColumn( zdcNo, "x", *ZDCdblVec->at(zdcNo).at(0) );
analysisManager->CreateNtupleDColumn( zdcNo, "y", *ZDCdblVec->at(zdcNo).at(1) );
analysisManager->CreateNtupleDColumn( zdcNo, "z", *ZDCdblVec->at(zdcNo).at(2) );
......@@ -111,8 +113,9 @@ RunAction::RunAction()
analysisManager->CreateNtupleIColumn( zdcNo, "charge", *ZDCintVec->at(zdcNo).at(7) );
} else { // There's only two branches to save space on cluster jobs
analysisManager->CreateNtupleDColumn( zdcNo, "lastStep" );
analysisManager->CreateNtupleIColumn( zdcNo, "radNo", *ZDCintVec->at(zdcNo).at(0) );
analysisManager->CreateNtupleIColumn( zdcNo, "nCherenkovs", *ZDCintVec->at(zdcNo).at(1) );
analysisManager->CreateNtupleIColumn( zdcNo, "nCherenkovs", *ZDCintVec->at(zdcNo).at(1) );
}//end if !CLUSTER
}//end ZDC loop
......@@ -122,6 +125,8 @@ RunAction::RunAction()
int nTuple = nZDCs + rpdNo;
if(!CLUSTER){
//Make double branches
analysisManager->CreateNtupleDColumn( nTuple, "lastStep" );
analysisManager->CreateNtupleDColumn( nTuple, "lastSteptest" );
analysisManager->CreateNtupleDColumn( nTuple, "x", *RPDdblVec->at(rpdNo).at(0) );
analysisManager->CreateNtupleDColumn( nTuple, "y", *RPDdblVec->at(rpdNo).at(1) );
analysisManager->CreateNtupleDColumn( nTuple, "z", *RPDdblVec->at(rpdNo).at(2) );
......@@ -145,6 +150,7 @@ RunAction::RunAction()
} else { // There's only two branches to save space on cluster jobs
//Make double branches
analysisManager->CreateNtupleDColumn( nTuple, "lastStep" );
analysisManager->CreateNtupleDColumn( nTuple, "x", *RPDdblVec->at(rpdNo).at(0) );
analysisManager->CreateNtupleDColumn( nTuple, "y", *RPDdblVec->at(rpdNo).at(1) );
analysisManager->CreateNtupleDColumn( nTuple, "z", *RPDdblVec->at(rpdNo).at(2) );
......
MonteCarlo/src/SteppingAction.cc
+ 13
− 18
View file @ b9b9dce3
......@@ -30,21 +30,23 @@
#include "SteppingAction.hh"
#include "EventAction.hh"
#include "Analysis.hh"
#include "G4Step.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "g4root.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
SteppingAction::SteppingAction(EventAction* eventAction)
SteppingAction::SteppingAction( )
: G4UserSteppingAction(){
// Get number of OPTICAL flag from DetectorConstruction
const DetectorConstruction* constDetectorConstruction
= static_cast<const DetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
DetectorConstruction* detectorConstruction
= const_cast<DetectorConstruction*>(constDetectorConstruction);
m_sd->AddOutputToRPDTree("LastStepInVolume", &lastStep);
OPTICAL = detectorConstruction->GetOpticalFlag();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -60,18 +62,11 @@ void SteppingAction::UserSteppingAction(__attribute__((unused)) const G4Step* th
if(theTrack->GetParentID()==0 && theStep->IsLastStepInVolume()){
lastStep = theTrack->GetPosition().getZ();
}
auto analysisManager = G4AnalysisManager::Instance();
for(int i = 0; i < analysisManager->GetNofNtuples(); i++{
analysisManager->FillNtupleDColumn(i,0,lastStep);
}
}
G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
G4StepPoint* thePrePoint = theStep->GetPreStepPoint();
G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
G4ParticleDefinition* particleType = theTrack->GetDefinition();
if(particleType==G4OpticalPhoton::OpticalPhotonDefinition() && thePostPV != NULL){
if( ( strncmp("ph",thePrePV->GetName().c_str(), 2) == 0) ){
//std::cout << theTrack->GetParticleDefinition()->GetParticleName() << std::endl;
theTrack->SetTrackStatus(fStopAndKill);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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