Added flux through spheres centered on the source for calculation of …

…diffusion length

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.
 
+
+

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";

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<double>& getEDEPvec() { return edepInHe3;}
   std::vector<int>& getPIDvec() { return PIDinHe3;}
   std::vector<int>& getNeutronHits() {return neutronHitVec;}
+
+  std::vector<double>& getDiffusionRadii() {return diffusionRadii;}
+  std::vector<double>& getCrossedSphere() {return nCrossedRadii;}
 
   int getnNeutrons() {return nNeutrons;}
   int getnGraphite() {return GraphiteTotal;}
@@ -74,6 +78,8 @@ private:
   std::vector<bool> isNeutronHitVec;
   std::vector<double> edepInHe3;
   std::vector<int> PIDinHe3;
+  std::vector<double> diffusionRadii;
+  std::vector<double> 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
 
 

include/SteppingAction.hh

@@ -23,6 +23,8 @@ private:
   int nch;
   int nobj;
   double cubeSize;
+
+  std::vector<double> radii;
 
   const DetectorConstruction* fDetConstruction;
   //SteppingActionMessenger*     fSteppingActionMessenger;

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<<bold2<<blue2;
     G4cout<<"Event number: "<<nevent<<G4endl;
-    cout<<nevent<<noFormat2;
+    cout<<noFormat2;
   }else if (eve->GetEventID()%1000==0){
     cout<<bold2;
     G4cout<<"Event number: "<<nevent<<G4endl;

src/RunAction.cc

@@ -61,6 +61,7 @@ RunAction::RunAction(DetectorConstruction* detConstruction)
   analysisManager->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<const EventAction*>(G4RunManager::GetRunManager()->GetUserEventAction());
+
+  const EventAction* constEventAction = static_cast<const EventAction*>(G4RunManager::GetRunManager()->GetUserEventAction());
   EventAction* eventAction = const_cast<EventAction*>(constEventAction);  
 
+
   int nevents = run->GetNumberOfEvent();
   double meanTime, stdev;
   eventAction->getTime(nevents, meanTime, stdev);

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((preR<radii[i]&&postR>radii[i])||(preR>radii[i]&&postR<radii[i])){
+	if(a1Track->GetKineticEnergy()/eV<1) fEventAction->crossedSphere(i);
+      }
+    }
+  }
+
+
+
 
   //if the particle leaves the concrete walls of the room
   if(post_volume ==  fDetConstruction->GetphysiRoom()){