diff --git a/Analysis/plotHisto.cc b/Analysis/plotHisto.cc
new file mode 100644
index 0000000000000000000000000000000000000000..36e38074f8702ce3ef6bacc388e24404be19c6d4
--- /dev/null
+++ b/Analysis/plotHisto.cc
@@ -0,0 +1,37 @@
+#include "TCanvas.h"
+#include "TH2.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TStyle.h"
+#include <iostream>
+
+using namespace std;
+
+int main(int argc, char *argv[]){
+ TFile *f = new TFile( Form("%s.root",argv[1]) );
+ if(f->IsZombie()){
+ cout << Form("%s.root does not exist... exiting",argv[1]) << endl;
+ return 0;
+ }
+
+ TTree *t = (TTree*)f->Get("lightGuide");
+ double x,y;
+ t->SetBranchAddress("X",&x);
+ t->SetBranchAddress("Y",&y);
+
+ TCanvas *c = new TCanvas("Data","Data",1280,720);
+ TH2D *h = new TH2D("eff","Relative Efficiency",40,-82,82,40,-42.875,42.875);
+
+ int nEvents = t->GetEntries();
+ for(int ev = 0; ev < nEvents; ev++){
+ t->GetEntry(ev);
+ h->Fill(y,x);
+ }
+
+ gStyle->SetOptStat(0);
+ h->Scale(1.0/h->GetMaximum());
+ h->Draw("COLZ");
+ c->Print( Form("%s.png",argv[1]) );
+ delete f;
+ return 1;
+}
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b6e982e6be1c610a57f8cfb6d93b8725318b31be..80139f99e85c4bf937de40c9d5a5f7f750d5c8e5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -46,6 +46,9 @@ file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
add_executable(lightGuide lightGuide.cc ${sources} ${headers})
target_link_libraries(lightGuide ${Geant4_LIBRARIES} ${ROOT_LIBRARIES} cadmesh )
+add_executable(plotHisto ${DEPENDENCIES} Analysis/plotHisto.cc)
+TARGET_LINK_LIBRARIES(plotHisto ${ROOT_LIBRARIES})
+
#----------------------------------------------------------------------------
# Copy all scripts to the build directory, i.e. the directory in which we
# build lightGuide. This is so that we can run the executable directly because it
@@ -68,6 +71,8 @@ endforeach()
# Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
#
install(TARGETS lightGuide DESTINATION bin)
-install(FILES run1.mac DESTINATION bin )
-install(FILES vis.mac DESTINATION bin )
-install(DIRECTORY models/ DESTINATION models FILES_MATCHING PATTERN "*.*")
+install(FILES run1.mac DESTINATION . )
+install(FILES run2.mac DESTINATION . )
+install(FILES runAndAnalyze.sh DESTINATION . )
+install(FILES vis.mac DESTINATION . )
+install(DIRECTORY models/ DESTINATION models)
diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh
index 01cd7271e5cf12377a11737502514aa10bf35a62..79c4193a9d94c15239f1e1319ce89e08f2f06b20 100644
--- a/include/DetectorConstruction.hh
+++ b/include/DetectorConstruction.hh
@@ -55,10 +55,12 @@ class DetectorConstruction : public G4VUserDetectorConstruction
void SetCADFiletype (std::string type){filetype = type;}
void SetGDMLoutName (std::string name){GDMLoutput = name;}
+ void SetSurfaceSigmaAlpha (G4double v);
void SetSurfaceFinish (const G4OpticalSurfaceFinish finish);
void SetSurfaceType (const G4SurfaceType type);
void SetSurfaceModel (const G4OpticalSurfaceModel model);
- void SetSurfaceSigmaAlpha (G4double v);
+ void AddSurfaceMPV (const char* c, G4MaterialPropertyVector* mpv);
+
void SetRotationX (G4int arg);
void SetRotationY (G4int arg);
void SetRotationZ (G4int arg);
diff --git a/include/Materials.hh b/include/Materials.hh
index 25170bd21e05d2aec0f66f241cb0406226ce7904..23dfc245e8410f741a1fc8cbacef85de4253d7ff 100644
--- a/include/Materials.hh
+++ b/include/Materials.hh
@@ -25,6 +25,8 @@ class Materials
/*! \brief Define indexes of interest.*/
void DefineOpticalProperties(void);
+ std::vector <G4MaterialPropertiesTable*> GetMPTArray(){return MPT_Array;}
+
//Elements
/* \brief Hydrogen */
G4Element* H;
diff --git a/run1.mac b/run1.mac
index e5e896d86e6ca7af4d0e0f507fb373d2f35d226a..710e96cd8038e10d82d15ae33981d5230afdc4a9 100644
--- a/run1.mac
+++ b/run1.mac
@@ -44,8 +44,8 @@
/gps/pos/type Plane
/gps/pos/shape Rectangle
/gps/pos/centre 0. 0. 0. mm
-/gps/pos/halfx 41.875 mm
-/gps/pos/halfy 80 mm
+/gps/pos/halfx 44.875 mm
+/gps/pos/halfy 82 mm
#
# the incident surface is in the x-y plane
diff --git a/run2.mac b/run2.mac
new file mode 100644
index 0000000000000000000000000000000000000000..bca9631f81cdb8706db2c0d2922d4f075c89075e
--- /dev/null
+++ b/run2.mac
@@ -0,0 +1,67 @@
+# Macro file for JZCaPA beam test 2018
+#
+# To be run preferably in batch, without graphics:
+#
+# 31 TeV Pb ion
+
+# Initialize kernel
+/run/initialize
+
+
+/control/verbose 0
+/run/verbose 0
+/run/printProgress 0
+/tracking/verbose 0
+/run/particle/verbose 0
+/hits/verbose 0
+
+/gps/particle opticalphoton
+
+#Set reflective surface properties
+/lightGuide/surfaceModel unified
+/lightGuide/surfaceType dielectric_metal
+/lightGuide/surfaceFinish ground
+#/lightGuide/surfaceSigmaAlpha .01
+#/lightGuide/surfaceProperty SPECULARLOBECONSTANT 0.000002 1.0 0.000008 1.0
+#/lightGuide/surfaceProperty SPECULARSPIKECONSTANT 0.000002 0.99 0.000008 0.99
+#/lightGuide/surfaceProperty BACKSCATTERCONSTANT 0.000002 .05 0.000008 .05
+#/lightGuide/surfaceProperty REFLECTIVITY 0.000002 .99 0.000008 .99
+
+# the beam energy is in gaussian profile
+/gps/ene/type Gauss
+/gps/ene/mono 3.4 eV
+/gps/ene/sigma 1 eV
+/gps/polarization 0 0 0
+
+# General particle source, beam
+#/gps/pos/type Beam
+#/gps/pos/shape Circle
+#/gps/pos/centre 0. 0. 0. mm
+#/gps/pos/radius 0.75 mm
+#/gps/pos/sigma_r 0.002 mm
+
+# General particle source, rectangular plane
+/gps/pos/type Plane
+/gps/pos/shape Rectangle
+/gps/pos/centre 0. 0. 0. mm
+/gps/pos/halfx 44.875 mm
+/gps/pos/halfy 82 mm
+
+#
+# the incident surface is in the x-y plane
+/gps/pos/rot1 1 0 0
+/gps/pos/rot2 0 1 0
+#
+#
+# the beam is travelling along the z-axis
+/gps/ang/rot1 0 1 0
+/gps/ang/rot2 1 0 0
+/gps/ang/type beam1d
+# the beam angular dispersion
+# N.A. is 0.22 so acceptance cone is ~25.4 degrees
+/gps/ang/sigma_r 25.40 deg
+# 1.0milliradians = 0.057 deg
+#
+
+#if using these commands in interactive mode to visualize, use this to allow more events to be displayed (do not use in batch mode)
+/vis/ogl/set/displayListLimit 1000000
diff --git a/runAndAnalyze.sh b/runAndAnalyze.sh
new file mode 100644
index 0000000000000000000000000000000000000000..c8649e257942f2b6e5a2102e0d8e1fa92b3a243a
--- /dev/null
+++ b/runAndAnalyze.sh
@@ -0,0 +1,16 @@
+#/bin/bash!
+
+cp run2.mac run3.mac
+
+echo "/lightGuide/surfaceSigmaAlpha 0.${1}" >> run3.mac
+echo "/lightGuide/surfaceProperty SPECULARLOBECONSTANT 0.000002 .${2} 0.000008 0.${2}" >> run3.mac
+echo "/lightGuide/surfaceProperty SPECULARSPIKECONSTANT 0.000002 0.${3} 0.000008 0.${3}" >> run3.mac
+echo '/lightGuide/surfaceProperty BACKSCATTERCONSTANT 0.000002 .05 0.000008 .05' >> run3.mac
+echo '/lightGuide/surfaceProperty REFLECTIVITY 0.000002 .99 0.000008 .99' >> run3.mac
+echo '/run/beamOn 1000000 ' >> run3.mac
+
+#./lightGuide -c models/Winston_cone.stl stl -m run3.mac -o Winston_SA0-$1
+./lightGuide -m run3.mac -o SA0-$1
+./plotHisto SA0-$1
+
+rm run3.mac
diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc
index 6452bbd6436b6302754172aa19584c4d1b3cbaec..b0cfca98465edb8722f5d436ccbeaf3737f2decb 100644
--- a/src/DetectorConstruction.cc
+++ b/src/DetectorConstruction.cc
@@ -314,6 +314,19 @@ void DetectorConstruction::SetSurfaceType(const G4SurfaceType type){
G4RunManager::GetRunManager()->GeometryHasBeenModified();
}
+/*
+ *
+ */
+void DetectorConstruction::AddSurfaceMPV(const char* c,
+ G4MaterialPropertyVector* mpv) {
+ mpv->SetSpline(true);
+ materials->GetMPTArray().at(1)->AddProperty(c, mpv);
+ materials->AlSurface->SetMaterialPropertiesTable(materials->GetMPTArray().at(1));
+ G4cout << "The MPT for the surface is now: " << G4endl;
+ materials->GetMPTArray().at(1)->DumpTable();
+ G4cout << "............." << G4endl;
+}
+
/*
*
*/
diff --git a/src/DetectorMessenger.cc b/src/DetectorMessenger.cc
index 2295729aa8fbaa47dd917562e12ecc9b398485e4..ea2372d14d2288b3ae23588fbf20a5a5ed163062 100644
--- a/src/DetectorMessenger.cc
+++ b/src/DetectorMessenger.cc
@@ -298,7 +298,29 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
FatalException,ed);
}
}
-
+ // SURFACE PROPERTY
+ else if (command == fSurfaceMatPropVectorCmd) {
+ // Convert string to physics vector
+ // string format is property name, then pairs of energy, value
+ // space delimited
+ G4MaterialPropertyVector* mpv = new G4MaterialPropertyVector();
+ G4cout << newValue << G4endl;
+ std::istringstream instring(newValue);
+ G4String prop;
+ instring >> prop;
+ while (instring) {
+ G4String tmp;
+ instring >> tmp;
+ if (tmp == "") { break; }
+ G4double en = G4UIcommand::ConvertToDouble(tmp);
+ instring >> tmp;
+ G4double val;
+ val = G4UIcommand::ConvertToDouble(tmp);
+ mpv->InsertValues(en, val);
+ }
+ const char* c = prop.c_str();
+ fDetector->AddSurfaceMPV(c, mpv);
+ }
// TYPE
else if (command == fSurfaceTypeCmd) {
if (newValue == "dielectric_metal") {
diff --git a/src/EventAction.cc b/src/EventAction.cc
index d26995eee63b63b587b8d2bbf14cb3ece8fcefff..08090fd01261fec4a792344519ebeee304867744 100644
--- a/src/EventAction.cc
+++ b/src/EventAction.cc
@@ -81,6 +81,7 @@ void EventAction::EndOfEventAction(const G4Event* event){
//analysisManager->FillNtupleDColumn(4, aHit->getEnergy() );
analysisManager->FillNtupleDColumn(4, aHit->getTime() );
analysisManager->AddNtupleRow();
+ analysisManager->FillH2(1,aHit->getPos().y(), aHit->getPos().x() );
}
}
}
diff --git a/src/RunAction.cc b/src/RunAction.cc
index 18b9a8333154d2c458dd183d1c31c980038c9f02..5063fa452047d40282a75b16c01c7409e6f43f92 100644
--- a/src/RunAction.cc
+++ b/src/RunAction.cc
@@ -36,6 +36,7 @@
#include "RunAction.hh"
#include "G4VAnalysisManager.hh"
+#include "G4SystemOfUnits.hh"
#include "lgAnalysis.hh"
@@ -65,9 +66,12 @@ RunAction::RunAction(G4String fileName)
analysisManager->CreateNtupleDColumn("Y");
analysisManager->CreateNtupleDColumn("hitX");
analysisManager->CreateNtupleDColumn("hitY");
- //analysisManager->CreateNtupleDColumn("energy");
analysisManager->CreateNtupleDColumn("time");
analysisManager->FinishNtuple();
+ analysisManager->SetFirstHistoId(1);
+ analysisManager->CreateH2("1","Relative Efficiency",
+ 150, -82.0*mm, 82.0*mm,
+ 150, -44.875*mm, 44.875*mm);
}
/*