From a85dfc7bc61c5469d0e73986f3aad26546188c42 Mon Sep 17 00:00:00 2001 From: srdejong Date: Tue, 15 Aug 2017 09:39:15 -0700 Subject: [PATCH] Added flux through spheres centered on the source for calculation of diffusion length --- README.1st | 5 ++++- include/Commands.h | 2 ++ include/EventAction.hh | 11 +++++++++++ include/SteppingAction.hh | 2 ++ src/EventAction.cc | 14 ++++++++++++-- src/RunAction.cc | 8 +++++++- src/SteppingAction.cc | 15 +++++++++++++++ 7 files changed, 53 insertions(+), 4 deletions(-) diff --git a/README.1st b/README.1st index 07d26a1..e681139 100644 --- a/README.1st +++ b/README.1st @@ -35,7 +35,8 @@ Command line parameters are optional. If none specified, PileRoomSim opens in vi Seeds for random number generator. -all - If this argument is used, all events are saved to the output ntuple. If not, only events where a neutron is captured are saved. + If this argument is used, all events are saved to the output ntuple. If not, only events where a neutron is captured are saved. + If you want to calculate the diffusion length of neutrons in the cube, this parameter must be used. -verbose Prints all the GEANT4 debugging information. @@ -47,3 +48,5 @@ Note: The environment variable XMLLOCATION should be set to the location of Grap In xml/, there is a python script called TubeCreator.py. This script can be used to create an xml file containing a tube description, or append a new description to an existing xml file. Default values (read from HE3TUBE-default.xml) are used for any parameter not specified. + + \ No newline at end of file diff --git a/include/Commands.h b/include/Commands.h index 132e425..b5b1c92 100644 --- a/include/Commands.h +++ b/include/Commands.h @@ -134,6 +134,8 @@ void help(){ G4cout<<" -all If this parameter is used, all events are saved to the \n"; G4cout<<" output ntuple. If not, only events containing a neutron\n"; G4cout<<" hit in in a helium-3 tube are saved.\n"; + G4cout<<" If you want to calculate the diffusion length of neutrons\n"; + G4cout<<" in the cube, this parameter must be used.\n"; G4cout<<" -novis If no other parameters specified, PileRoomSum runs in\n"; G4cout<<" interactive mode with no visualization set.\n"; G4cout<<" -verbose Prints GEANT4 information\n"; diff --git a/include/EventAction.hh b/include/EventAction.hh index 32956f6..9dbd1bd 100644 --- a/include/EventAction.hh +++ b/include/EventAction.hh @@ -36,6 +36,7 @@ public: inline void leftGraphite(double eKin); inline void leftWall(); inline void leftObject(int n); + inline void crossedSphere(int n); void setdataNtuple(int n) {dataNtuple=n;} void setGeoNtuple(int n) {geoNtuple=n;} @@ -44,6 +45,9 @@ public: std::vector& getEDEPvec() { return edepInHe3;} std::vector& getPIDvec() { return PIDinHe3;} std::vector& getNeutronHits() {return neutronHitVec;} + + std::vector& getDiffusionRadii() {return diffusionRadii;} + std::vector& getCrossedSphere() {return nCrossedRadii;} int getnNeutrons() {return nNeutrons;} int getnGraphite() {return GraphiteTotal;} @@ -74,6 +78,8 @@ private: std::vector isNeutronHitVec; std::vector edepInHe3; std::vector PIDinHe3; + std::vector diffusionRadii; + std::vector nCrossedRadii; double TotalEnergyDeposit; double tubeX; double tubeY; @@ -133,6 +139,11 @@ inline void EventAction::leftObject(int n){ leftObj[n] = 1; } +inline void EventAction::crossedSphere(int n){ + nCrossedRadii[n] += 1/(4*3.14159*pow(diffusionRadii[n],2)); +} + + #endif diff --git a/include/SteppingAction.hh b/include/SteppingAction.hh index 79f36c5..a80449d 100644 --- a/include/SteppingAction.hh +++ b/include/SteppingAction.hh @@ -23,6 +23,8 @@ private: int nch; int nobj; double cubeSize; + + std::vector radii; const DetectorConstruction* fDetConstruction; //SteppingActionMessenger* fSteppingActionMessenger; diff --git a/src/EventAction.cc b/src/EventAction.cc index e980964..138dac8 100644 --- a/src/EventAction.cc +++ b/src/EventAction.cc @@ -30,7 +30,12 @@ EventAction::EventAction(const DetectorConstruction* detectorConstruction, bool timeStdev(0) { //fillGeoNtuple(); - + //diffusionRadii.resize(10); + for(int i=0; i<100; i++){ + diffusionRadii.push_back(300.+700.*(double)i/100.); + } + + } @@ -64,6 +69,11 @@ void EventAction::BeginOfEventAction( const G4Event* eve) leftGrap=false; //hasn't left graphite leftConcrete=false; //hasn't left concrete ePostGraphite=-1; + + + nCrossedRadii.clear(); + nCrossedRadii.resize(diffusionRadii.size()); + leftObj.clear(); //hasn't left objects leftObj.resize(nobj); @@ -88,7 +98,7 @@ void EventAction::BeginOfEventAction( const G4Event* eve) }else if(eve->GetEventID()%10000==0){ cout<GetEventID()%1000==0){ cout<CreateNtupleDColumn(dataNtuple,"he3TubeXPos"); analysisManager->CreateNtupleDColumn(dataNtuple,"he3TubeYPos"); analysisManager->CreateNtupleDColumn(dataNtuple,"he3TubeZPos"); + @@ -78,6 +79,8 @@ RunAction::RunAction(DetectorConstruction* detConstruction) analysisManager->CreateNtupleIColumn(dataNtuple,"PIDinHe3", eventAction->getPIDvec()); analysisManager->CreateNtupleIColumn(dataNtuple,"neutronHits", eventAction->getNeutronHits()); + analysisManager->CreateNtupleDColumn(dataNtuple, "diffusionRadius", eventAction->getDiffusionRadii()); + analysisManager->CreateNtupleDColumn(dataNtuple, "diffusionFlux", eventAction->getCrossedSphere()); //no more branches added after this @@ -132,6 +135,7 @@ void RunAction::EndOfRunAction(const G4Run* run) //save ntuple G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); + analysisManager->Write(); analysisManager->CloseFile(); } @@ -277,9 +281,11 @@ void RunAction::FillGeoNtuple(G4AnalysisManager* analysisManager, const G4Run* a //print an end of run report void RunAction::report(const G4Run* run){ -const EventAction* constEventAction = static_cast(G4RunManager::GetRunManager()->GetUserEventAction()); + + const EventAction* constEventAction = static_cast(G4RunManager::GetRunManager()->GetUserEventAction()); EventAction* eventAction = const_cast(constEventAction); + int nevents = run->GetNumberOfEvent(); double meanTime, stdev; eventAction->getTime(nevents, meanTime, stdev); diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc index 1dc8777..29165ff 100644 --- a/src/SteppingAction.cc +++ b/src/SteppingAction.cc @@ -32,6 +32,8 @@ SteppingAction::SteppingAction(const DetectorConstruction* detectorConstruction, nch = fDetConstruction->GetNChannels(); nobj = fDetConstruction->GetNMiscObjects(); cubeSize = fDetConstruction->GetCubeSize(); + radii = fEventAction->getDiffusionRadii(); + } @@ -57,6 +59,19 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) } } + if(PDG==2112){ + double preR = pre_posn.r(); + double postR = post_posn.r(); + + for(int i=0; i<(int)radii.size(); i++){ + if((preRradii[i])||(preR>radii[i]&&postRGetKineticEnergy()/eV<1) fEventAction->crossedSphere(i); + } + } + } + + + //if the particle leaves the concrete walls of the room if(post_volume == fDetConstruction->GetphysiRoom()){