Import Geant4 11.0.0 source tree

examples/advanced/doiPET/CMakeLists.txt

@@ -1,93 +1,103 @@
-#----------------------------------------------------------------------------
-# Setup the project
-#
-cmake_minimum_required(VERSION 3.12...3.20)
-project(doiPET)
-
-#----------------------------------------------------------------------------
-# Find Geant4 package, activating all available UI and Vis drivers by default
-# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
-# to build a batch mode only executable
-#
-option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
-if(WITH_GEANT4_UIVIS)
-  find_package(Geant4 REQUIRED ui_all vis_all)
-else()
-  find_package(Geant4 REQUIRED)
-endif()
-
-# Setup of ROOT ANALYSIS : optional.
-option(WITH_ANALYSIS_USE "Build example with analysis objects" OFF)
-if(WITH_ANALYSIS_USE)
-add_definitions(-DUSEROOT)
-else()
-endif()
-#----------------------------------------------------------------------------
-# Setup Geant4 include directories and compile definitions
-# Setup include directory for this project
-#
-include(${Geant4_USE_FILE})
-include_directories(${PROJECT_SOURCE_DIR}/include)
-
-#----------------------------------------------------------------------------
-# Find ROOT (required package)
-#
-find_package(ROOT REQUIRED)
-if(WITH_ANALYSIS_USE)
-	EXECUTE_PROCESS(COMMAND root-config --cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS OUTPUT_STRIP_TRAILING_WHITESPACE)
-	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${ROOT_CXX_FLAGS}")
-endif(WITH_ANALYSIS_USE)
-#----------------------------------------------------------------------------
-# Locate sources and headers for this project
-# NB: headers are included so they will show up in IDEs
-#
-#file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc)
-#file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
-
-include_directories(${PROJECT_SOURCE_DIR}/include
-                    ${PROJECT_SOURCE_DIR}/shared/include
-                    ${Geant4_INCLUDE_DIR}
-                    ${ROOT_INCLUDE_DIR})
-file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc
-                  ${PROJECT_SOURCE_DIR}/shared/src/*.cc)
-file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh
-                  ${PROJECT_SOURCE_DIR}/shared/include/*.hh)
-
-#----------------------------------------------------------------------------
-# Add the executable, and link it to the Geant4 libraries
-#
-add_executable(doiPET doiPET.cc ${sources} ${headers})
-
-if(WITH_ANALYSIS_USE)
-	EXECUTE_PROCESS(COMMAND root-config --libs OUTPUT_VARIABLE ROOT_LD_FLAGS OUTPUT_STRIP_TRAILING_WHITESPACE)
-	set(CMAKE_EXE_LINKER_FLAGS ${ROOT_LD_FLAGS})
-endif(WITH_ANALYSIS_USE)
-
-target_link_libraries(doiPET ${Geant4_LIBRARIES} ${ROOT_LIBRARIES})
-
-#----------------------------------------------------------------------------
-# Copy all scripts to the build directory. This is so that we can run the executable directly because it
-# relies on these scripts being in the current working directory.
-#
-set(DOIPET_SCRIPTS
-  init_vis.mac
-  vis.mac
-  run_imageQualityPhantom_wholeBody.mac
-  run_imageQualityPhantom_smallAnimal.mac
-  run_NECR.mac
-  run_spatialResolution.mac
-  run_sensitivity.mac
-  run_normalization.mac
-  run.mac
-  inputParameter.txt
-  )
-
-foreach(_script ${DOIPET_SCRIPTS})
-  configure_file(
-    ${PROJECT_SOURCE_DIR}/${_script}
-    ${PROJECT_BINARY_DIR}/${_script}
-    COPYONLY
-    )
-endforeach()
-
-install(TARGETS doiPET DESTINATION bin )
+
+#----------------------------------------------------------------------------
+# Setup the project
+#
+cmake_minimum_required(VERSION 3.16...3.21)
+project(doiPET)
+
+#----------------------------------------------------------------------------
+# Find Geant4 package, activating all available UI and Vis drivers by default
+# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
+# to build a batch mode only executable
+#
+option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
+if(WITH_GEANT4_UIVIS)
+  find_package(Geant4 REQUIRED ui_all vis_all)
+else()
+  find_package(Geant4 REQUIRED)
+endif()
+
+# Setup of ROOT ANALYSIS : optional. 
+
+#option(WITH_ANALYSIS_USE "Build example with analysis objects" OFF)
+#if(WITH_ANALYSIS_USE)
+#add_definitions(-DUSEROOT)
+#else()
+#endif()
+
+option(WITH_ANALYSIS_USE "Build example with analysis objects" OFF)
+if(WITH_ANALYSIS_USE)
+add_definitions(-DANALYSIS_USE)
+else()
+endif()
+
+#----------------------------------------------------------------------------
+# Setup Geant4 include directories and compile definitions
+# Setup include directory for this project
+#
+include(${Geant4_USE_FILE})
+include_directories(${PROJECT_SOURCE_DIR}/include)
+
+#----------------------------------------------------------------------------
+# Find ROOT (required package)
+#
+#find_package(ROOT REQUIRED)
+#if(WITH_ANALYSIS_USE)
+#	EXECUTE_PROCESS(COMMAND root-config --cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS OUTPUT_STRIP_TRAILING_WHITESPACE)
+#	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${ROOT_CXX_FLAGS}")
+#endif(WITH_ANALYSIS_USE) 
+#----------------------------------------------------------------------------
+# Locate sources and headers for this project
+# NB: headers are included so they will show up in IDEs
+#
+file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc)
+file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
+
+#include_directories(${PROJECT_SOURCE_DIR}/include 
+#                    ${PROJECT_SOURCE_DIR}/shared/include 
+#                    ${Geant4_INCLUDE_DIR}
+#                    ${ROOT_INCLUDE_DIR})
+#file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc
+#                  ${PROJECT_SOURCE_DIR}/shared/src/*.cc)
+#file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh
+#                  ${PROJECT_SOURCE_DIR}/shared/include/*.hh)
+
+#----------------------------------------------------------------------------
+# Add the executable, and link it to the Geant4 libraries
+#
+add_executable(doiPET doiPET.cc ${sources} ${headers})
+
+#if(WITH_ANALYSIS_USE)
+#	EXECUTE_PROCESS(COMMAND root-config --libs OUTPUT_VARIABLE ROOT_LD_FLAGS OUTPUT_STRIP_TRAILING_WHITESPACE)
+#	set(CMAKE_EXE_LINKER_FLAGS ${ROOT_LD_FLAGS})
+#endif(WITH_ANALYSIS_USE)
+
+target_link_libraries(doiPET ${Geant4_LIBRARIES})
+
+#----------------------------------------------------------------------------
+# Copy all scripts to the build directory. This is so that we can run the executable directly because it
+# relies on these scripts being in the current working directory.
+#
+set(DOIPET_SCRIPTS
+  init_vis.mac
+  vis.mac
+  run_imageQualityPhantom_wholeBody.mac
+  run_imageQualityPhantom_smallAnimal.mac
+  run_NECR.mac
+  run_spatialResolution.mac
+  run_sensitivity.mac 
+  run_normalization.mac 
+  run.mac
+  inputParameter.txt
+  analysis.cpp
+  )
+
+foreach(_script ${DOIPET_SCRIPTS})
+  configure_file(
+    ${PROJECT_SOURCE_DIR}/${_script}
+    ${PROJECT_BINARY_DIR}/${_script}
+    COPYONLY
+    )
+endforeach()
+
+install(TARGETS doiPET DESTINATION bin )

examples/advanced/doiPET/History

@@ -6,6 +6,26 @@
                       Category History file
                       ---------------------
 
+14.11.2021 - G. Cosmo;doiPET-V10-07-03
+           - Fixed compilation warnings for unused private variables in
+             doiPETDetectorConstruction and doiPETEventAction. Removed tabs.
+
+16.10.2021 - S. Guatelli;doiPET-V10-07-02
+           - Migration to the G4RunManagerFactory
+           - Analysis now done with G4Analysis (no dependency towards external ROOT libraries)
+           - Link to reference data and analysis scripts provided to 
+             let users compare simulatilesion results to reference data
+           - Some controls to prepare DOI-look-up table is added. 
+             If DOI look-up table is provided as a file, then it will use it. 
+             If look-up table is not provided, the program will prepare the look-up table
+              based on the information provided by the user (number of crystals, etc).
+           - In performing Anger Logic calculation, the interaction position is mapped to the center
+             of the crystal before 2D position error is added to model photomultiplier tube (PMT) error.
+           - An analysis C++ program is provided to calculate the sensitivity of the scanner for a line source.
+
+07.10.2021 - I. Hrivnacova (doiPET-V10-07-01)
+             Removed unused include of g4root.hh
+
 24.05.2021 - B. Morgan (doiPET-V10-07-00)
              Bump required CMake version range to 3.12...3.20, matching core Geant4
 

examples/advanced/doiPET/README

@@ -1,112 +0,0 @@
-===========================================================
----------------Geant4 doiPET example---------------------
-===========================================================
- Author list to be updated, with names of co-authors and contributors from National Institute of Radiological Sciences (NIRS)
-
- Abdella M. Ahmed (1, 2), Andrew Chacon (1, 2), Harley Rutherford (1, 2),
- Hideaki Tashima (3), Go Akamatsu (3), Akram Mohammadi (3), Eiji Yoshida (3), Taiga Yamaya (3)
- Susanna Guatelli (2), and Mitra Safavi-Naeini (1, 2)
-
-(1) Australian Nuclear Science and Technology Organisation, Australia
-(2) University of Wollongong, Australia
-(3) National Institute of Radiological Sciences, Japan
-
-================================================================================================
-
-This example is for simulating depth-of-interaction enabled positron emission tomography (PET) scanner
-
-1-GEOMETRY 
-The detector construction has two main parts: constructing the PET system and the phantoms.
-
-The PET system is constructed from depth-of-interaction (DOI)detectors blocks. Each detector consisted of 16 x 16 x 4 crystal array constructed from GSO scintillation material. Material are defined in the DefineMaterials() using Geant4 NIST database. The geometrical specifications are given (and can be changed) in the GlobalParameters.hh file. 
-
-The scanner has 4 ring detectors. The detectors are covered with Aluminum material. Gaps between crystal elements, as well as adjacent rings are introduced.
-
-Various types of NEMA NU phantoms has been provided and are defined in the ConstructPhantom() method. To precisely create the image quality phantom, the G4UnionSolid from the Constructive Solid Geometry (CSG) has been used. The type, position and size of the phantoms can be changed using the macro file when necessary. A macro file is provided for each type of phantom imaging. For example, to run the simulation with image quality phantom, the run_imageQualityPhantom_wholeBody.mac should be used.
-	
- 2- PHYSICS LIST
-
-The physics list contains standard electromagnetic processes and the RadioactiveDecay module for GenericIon. It is defined in the PhysicsList class as a Geant4 modular physics list with registered physics builders provided in Geant4:
-   - G4DecayPhysics - defines all particles and their decay processes
-   - G4RadioactiveDecayPhysics - defines radioactiveDecay for GenericIon
-   - G4EmStandardPhysics_option3 - defines EM standard processes
-    
- 3- ACTION INITALIZATION
-
-   The ActionInitialization class instantiates and registers to Geant4 kernel all user action classes by invoking the ActionInitialization::Build().
-
- 4- PRIMARY GENERATOR
-
-The default particle beam is F-18 ion at rest defined in the GPS (General particle Source). The GPS is used for all types of activity distribution. Various macro files are provided with the name appended on it for specific simulation. The following macro files are provided:
- 
-run_imageQualityPhantom_wholeBody.mac
-run_imageQualityPhantom_smallAnimal.mac
-run_NECR.mac
-run_sensitivity.mac
-run_spatialResolution.mac
-run_normalization.mac (This one is not given in the NEMA NU manual but it is an important part of image reconstruction) 
-   
-5-EVENT ACTION
-
-At the end of each event, the information is extracted by calling FindInteractingCrystal() function and associative container (multimap and set methods) and the containers are cleared by calling the Clear() function.
-
-6-RUN ACTION
-
-The RunAction class is used mainly to conform all the events are processed. At the beginning of the run, a file is created and opened to write the output file and the file is closed at the end of the run and other instances are deleted.
-
-   		
- 7- STEPPING ACTION
-
-The SteppingAction class is the one which is used to track the steps. In the stepping action, interaction information of the photon with the crystal and the phantoms are extracted. The interaction information (such as energy deposition, blockID, crystalID, etc) is passed into the Analysis.cc class, which outputs the result into an ASCII file. 
-
-Generation of the source (F-18 ion) is confined in the physical volume by killing the event in the SteppingAction class when it is out of the physical volume.
-
-8-ANALYSIS
-
-The interaction information from SteppingAction is passed to the Analysis class. In this class, multiple mapping of the interaction information is stored in the multimap (associative container) by taking the blockID as key value. The information obtained (such as the energy deposition, etc) is then blurred before writing it into the output file. 
-Blurring parameters are given and can be changed in the inputparameters.txt file. The inputparameter.txt file include some characteristics of the scanner such as: 
-
-Crystal dependent energy resolution: this can be set between a minimum and a maximum value. The crystals are not assigned a constant energy resolution.
-
-b)	Dead time: the default dead time for each block is 256 ns. An option for axially arranged (multiplexed) detectors can also be set (the default value is 0 ns).
-
-c)	Detection efficiency (Quantum efficiency): This parameter represents the effect of the transfer efficiency of the crystal and of the quantum efficiency of the photo-detector. To model the efficiency of the system, a coefficient (between 0 and 1) can be set. 
-
-d)	Reflector insertion: In addition to the 2D plane (tangential and axial dirction), the interaction points along the depth direction are identified by DOI technique. One method is to control the optical photons by means of reflectors. Since it takes a very long time to simulate optical photons, the response due to reflector is modeled without real reflector insertion. This is done by reading reflector patterns for each layer and then shift the response based on the existence of the reflector between crystals. The reflector patter is given in the inputparameter.txt file.
-
-(Virtual) photomultiplier tube (PMT)
-
-Four PMTs are assumed to be placed at each corners of the crystal block. These PMTs are linear in that the signals are calculated based on their distances from the interaction position.
-
-Anger Logic calculation
-
-In PET analysis, the lines of response (LOR) needed for the image reconstruction is drawn between two opposing crystals. In simulation, these crystals are identified based on the highest energy deposition. In reality, however, errors are introduced in identifying the crystals which includes the Anger Logic calculation. Without additional blurring of the crystal, simulation results will always have a better spatial resolution than experimental measurements. To include the blurring in identifying the crystals, we included the Anger Logic calculation method. The steps in identifying the crystals are as follows:
- 
-Crystal ID in 3D (tangential, axial and DOI directions) are identified with the highest energy deposition and this information along with energy deposition and position interaction is passed into the AngerLogic(c) function.
-Signals (energy deposition) on the PMTs are obtained based on the lateral distance of the PMTs from the interaction position.
-Anger Logic position calculation is performed based on the signals on each PMT to obtain position interaction in 2D. The response of the PMTs are then shifted based on whether reflector exists between crystals or not. Note that, each layer has different reflector pattern.
-The new crystal ID is then calculated by comparing the look-up-table which is provided in the code based on a publication from NIRS. 
- 
-***** ROOT ANALYSIS
-To have ROOT Analysis, Be in the build director
-/Path/doiPET/build/ and type:
-cmake -DWITH_ANALYSIS_USE=ON -DGeant4_DIR=.... ../
-
-
-***** How to run a  simulation:
-
-Be in the build director
-/Path/doiPET/build/ cmake ../
-/Path/doiPET/build/ make
-/Path/doiPET/build/ ./doiPET run.mac
-
-Simulation output:
-
-ASCII and ROOT files are created depending on the type of the output format. The following information of the event is written in the output file:
-
-EventID, BlockID, tangentialCrystalID, AxialCrystalID, DOI_ID, time, and Energy deposition in the crystal is written to the file as a list-mode format. 
-
-The user can choose to make the output either in singles or coincidence mode in the inputParameter.txt file as follows: 
-
-
-=================== end ====================

examples/advanced/doiPET/analysis.cpp

@@ -0,0 +1,144 @@
+//The following code is for anylizing axial sensitivity from coincidence list-mode data
+//It takes source position of the event (of those which are detected by the scanner) and analises axial sensitivty.
+//by Abdella M. Ahmed, 2020
+
+#define _USE_MATH_DEFINES
+#include <iostream>
+#include <cfloat>
+#include <cmath>
+#include <fstream>
+#include <string>
+#include <limits>
+#include <stdio.h>
+#include <random>
+#include <string>
+
+#define UseROOT //ASCII or UseROOT
+
+#ifdef UseROOT
+//for root
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TNtuple.h"
+#include "TROOT.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TTree.h"
+#include "TLeaf.h"
+#include "TSystem.h"
+#endif
+
+
+
+
+#define AxialLength 216 //(mm), Axial length of the scanner
+
+//Change this number based of the number of primary particles simulated in the "run.mac" file
+#define NumberOfPrimaries 3000000 //Number of particles simulated
+
+using namespace std;
+int main(){
+
+	cout<<"============Axial Sensitivity Analysis =========================="<<endl;
+
+	int eventID0, blockID0, crystalID_axial0, crystalID_tangential0, DOI_ID0; 
+	double timeStamp0, totalEdep0;
+    int eventID1, blockID1, crystalID_axial1, crystalID_tangential1, DOI_ID1;
+	double timeStamp1, totalEdep1;
+	double  spositionX, spositionY, spositionZ; //source position
+
+	double z_offset = 0.0;//Axial offset position where the plane is located
+	double planeWidth = 3;// (mm)
+	int planeNumber;
+
+	float total_sensitivity = 0.0;
+
+	int numberOfSlices = int (AxialLength/planeWidth);
+	cout<<"Number of axial planes (slices) are: " <<numberOfSlices<<endl;
+	double Counts_per_plane[numberOfSlices];
+
+	ofstream OutFile("axial_sensitivity.csv");
+	string filename = "resultCoincidence.data";
+	string filepath = "";//provide the file path if it is stored in a different location.
+	ifstream InFile(filepath+filename);
+
+	if (!InFile.is_open())
+	{
+		cout << "Unable to open input file to read .... " << endl;
+		return 1;
+	}
+	if (!OutFile.is_open())
+	{
+		cout << "Unable to open out file to write ....";
+		return 1;
+	}
+
+	OutFile << "PlaneNmber" << "," << "Z(mm)" << "," << "Sensitivity(%)" << endl;
+
+	for (int i_plane = 0; i_plane < numberOfSlices; i_plane++){
+		Counts_per_plane[i_plane] = 0;
+	}
+#ifdef ASCII
+    cout<<"\nASCII coincidence list-mode data is being analised..."<<endl;
+	while (InFile >> eventID0 >> blockID0 >> crystalID_axial0 >> crystalID_tangential0 >> DOI_ID0 >> timeStamp0 >> totalEdep0 >> 
+		             eventID1 >> blockID1 >> crystalID_axial1 >> crystalID_tangential1 >> DOI_ID1 >> timeStamp1 >> totalEdep1 >> 
+					 spositionX >> spositionY >> spositionZ){
+		
+		planeNumber = int(spositionZ/planeWidth + numberOfSlices/2 - 0.5 + 0.5);
+		Counts_per_plane[planeNumber]++;
+		total_sensitivity++;
+	}
+#endif
+	
+	
+#ifdef UseROOT
+	cout<<"\nROOT coincidence list-mode data is being analised..."<<endl;
+	TFile *f = new TFile("resultCoincidence.root","READ");
+	TTree *Singles = (TTree*)gDirectory->Get("Coincidence");
+
+	Singles->SetBranchAddress("eventID0",&eventID0);
+	Singles->SetBranchAddress("blockID0",&blockID0);
+	Singles->SetBranchAddress("crystalID_axial0",&crystalID_axial0);
+	Singles->SetBranchAddress("crystalID_tangential0",&crystalID_tangential0);
+	Singles->SetBranchAddress("DOI_ID0",&DOI_ID0);
+	Singles->SetBranchAddress("timeStamp0",&timeStamp0);
+	Singles->SetBranchAddress("totalEdep0",&totalEdep0);
+
+	Singles->SetBranchAddress("eventID1",&eventID1);
+	Singles->SetBranchAddress("blockID1",&blockID1);
+	Singles->SetBranchAddress("crystalID_axial1",&crystalID_axial1);
+	Singles->SetBranchAddress("crystalID_tangential1",&crystalID_tangential1);
+	Singles->SetBranchAddress("DOI_ID1",&DOI_ID1);
+	Singles->SetBranchAddress("timeStamp1",&timeStamp1);
+	Singles->SetBranchAddress("totalEdep1",&totalEdep1);
+
+	Singles->SetBranchAddress("spositionX",&spositionX);
+	Singles->SetBranchAddress("spositionY",&spositionY);
+	Singles->SetBranchAddress("spositionZ",&spositionZ);
+
+	//
+	int nentries = 0;
+	nentries = Singles->GetEntries();
+	for(int entry=0; entry<nentries; entry++){
+		Singles->GetEntry(entry);      
+		planeNumber = int(spositionZ/planeWidth + numberOfSlices/2 - 0.5 + 0.5);
+		Counts_per_plane[planeNumber]++;
+		total_sensitivity++;
+	}
+#endif 
+
+	
+	//Save it into CSV file
+	for (int i_plane = 0; i_plane < numberOfSlices; i_plane++){
+		//Axial mid-position of the plane or slice
+		z_offset = (i_plane - numberOfSlices/2 + 0.5)*planeWidth;
+
+		OutFile << i_plane << "," << z_offset << "," << (Counts_per_plane[i_plane]/NumberOfPrimaries)*100 << endl;
+	}
+	cout<<"\nSensitivity evaluation has completed."<<endl;
+	cout<<"\nTotal Sensitivity (Number of coinsidence per total number of pair of photons): "<<(total_sensitivity/NumberOfPrimaries)*100 << "%"<<endl;
+	InFile.close();
+	OutFile.close();
+
+	return 0;
+}

examples/advanced/doiPET/doiPET.cc

@@ -36,52 +36,43 @@
 // (2) University of Wollongong, Australia
 // (3) National Institute of Radiological Sciences, Japan
 
-
-
-//#include "doiPETGlobalParameters.hh"
 #include "doiPETDetectorConstruction.hh"
 #include "doiPETPhysicsList.hh"
 #include "doiPETAnalysis.hh"
 #include "doiPETActionInitialization.hh"
-
 #include "Randomize.hh"
 #include "G4UImanager.hh"
 #include "G4VisExecutive.hh"
 #include "G4UIExecutive.hh"
 #include "G4SystemOfUnits.hh"
-
-#ifdef G4MULTITHREADED
-#include "G4MTRunManager.hh"
-#else
-#include "G4RunManager.hh"
-#endif
-
+#include "G4RunManagerFactory.hh"
+//
 //////////////////////////////////////////////////////////////////////////////
 
 int main(int argc,char** argv)
 {
+  auto* runManager = G4RunManagerFactory::CreateRunManager();
+  G4int nThreads = 4;
+  runManager->SetNumberOfThreads(nThreads);
+   
+	G4UIExecutive* ui = 0;
+	if ( argc == 1 ) {
+    ui = new G4UIExecutive(argc, argv);
+	}
+	
 	// Choose the Random engine
+	
 	G4Random::setTheEngine(new CLHEP::RanecuEngine);
 
-#ifdef G4MULTITHREADED
-
-	G4MTRunManager* runManager = new G4MTRunManager;
-	runManager->SetNumberOfThreads(4); // Is equal to 2 by default
-#else
-
-	G4RunManager* runManager = new G4RunManager;
-#endif
-
-	runManager->SetUserInitialization(new doiPETDetectorConstruction);
-
-	runManager->SetUserInitialization(new doiPETPhysicsList);
-
-	// Set user action initialization
-	runManager->SetUserInitialization(new doiPETActionInitialization());
-
 	//Initialize analysis
 	doiPETAnalysis* ptrAnalysis = doiPETAnalysis::GetInstance();
 
+	runManager->SetUserInitialization(new doiPETDetectorConstruction());
+
+	runManager->SetUserInitialization(new doiPETPhysicsList());
+
+	// Set user action initialization
+	runManager->SetUserInitialization(new doiPETActionInitialization(ptrAnalysis));
 
 	G4double act  = 1000000 * becquerel;//Activity is set via run.mac file
 	ptrAnalysis->SetActivity(act);
@@ -102,30 +93,23 @@ int main(int argc,char** argv)
 	ptrAnalysis->ReadReflectorPattern();
 
 	// Get the pointer to the User Interface manager
-	G4UImanager* UI = G4UImanager::GetUIpointer();
+  G4UImanager* UImanager = G4UImanager::GetUIpointer();
 
-
-
-	// Process macro or start UI session
-	if (argc!=1)   // batch mode  
-	{
-		G4String command = "/control/execute ";
-		G4String fileName = argv[1];
-		UI->ApplyCommand(command+fileName);
-	}
-
-	else           //define visualization and UI terminal for interactive mode
-	{ 
-		G4VisManager* visManager = new G4VisExecutive;
-		visManager->Initialize();
-
-		G4UIExecutive * ui = new G4UIExecutive(argc,argv);      
-		ui->SessionStart();
-		delete ui;
-
-		//
-		delete visManager;  
-	}
+  if ( ! ui ) { 
+    // batch mode
+    G4String command = "/control/execute ";
+    G4String fileName = argv[1];
+    UImanager->ApplyCommand(command+fileName);
+  }
+  else { 
+    // interactive mode
+	G4VisManager* visManager = new G4VisExecutive;
+    visManager->Initialize();
+    UImanager->ApplyCommand("/control/execute init_vis.mac");
+    ui->SessionStart();
+    delete ui;
+	delete visManager;
+  }
 
 	//close the file
 	ptrAnalysis->Close();

examples/advanced/doiPET/include/doiPETActionInitialization.hh

@@ -42,12 +42,13 @@
 #define doiPETActionInitialization_h 1
 
 #include "G4VUserActionInitialization.hh"
-
+#include "doiPETAnalysis.hh"
 //
+class doiPETAnalysis;
 class doiPETActionInitialization : public G4VUserActionInitialization
 {
 public:
-	doiPETActionInitialization();
+	doiPETActionInitialization(doiPETAnalysis* );
 	virtual ~doiPETActionInitialization();
 
 	virtual void BuildForMaster() const;
@@ -56,6 +57,7 @@ public:
 	//virtual G4VSteppingVerbose* InitializeSteppingVerbose() const;
 
 private:
+	doiPETAnalysis* analysis;
 	// doiPETDetectorConstruction* fDetector;
 };
 //

examples/advanced/doiPET/include/doiPETAnalysis.hh

@@ -54,23 +54,10 @@
 #include <iterator>
 #include <vector>
 #include <algorithm> 
+#include "G4AnalysisManager.hh"
 
-#ifdef USEROOT
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wshadow"
-#include "TRandom3.h"
-#include "TFile.h"
-#include "TNtuple.h"
-#include "g4root.hh"
-#include "TROOT.h"
-#include "TFile.h"
-#include "TH1.h"
-#include "TH2.h"
-#include "TTree.h"
-#include "TLeaf.h"
-#include "TSystem.h"
-#pragma GCC diagnostic pop
-#endif
+// Define the total number of columns in the ntuple
+const G4int MaxNtCol = 17;
 
 class doiPETAnalysisMessenger;
 
@@ -96,17 +83,19 @@ public:
 
 	void GetParentParticleName(G4String);
 	void GetSizeOfDetector (G4double, G4double, G4double);
-	void GetScatterIndexInPhantom(G4double);
+	void SetScatterIndexInPhantom(G4int);
 
 	void SetSourcePosition(G4ThreeVector);//
 	void SetEventID(G4int);
 
 	void BlurringParameters();
 	void GetTimeOfAnnihilation(G4double);
+	
 
 	void PMTPosition();
-	void AngerLogic(G4int, G4int, G4int, G4double, G4double, G4double, G4double);
+	void AngerLogic(G4double, G4double, G4double, G4double, G4double, G4bool);
 	void ReadReflectorPattern();
+	void PrepareDOILookUpTable(G4String);
 
 	void SetActivity(G4double);
 	void SetIsotopeHalfLife(G4double);
@@ -117,14 +106,28 @@ public:
 	G4double QuantumEffifciency(G4double, G4int, G4int);
 	void ReadOut(G4int, G4int, G4double, G4double, G4ThreeVector, G4double);
 
+	//G4ROOT
+	void book(); // booking the ROOT file
+
+	
+	void FillListModeEvent(); //Single or Coinsidence 
+	void finish();
+	// Close the ROOT file with all the results stored in nutples 
+
+
 private:
 	static doiPETAnalysis* instance;
 	doiPETAnalysisMessenger* fAnalysisMessenger;
 	//std::multimap< G4int, InteractionInformation* > mapBlockInteraction;
 	std::set<G4int> setBlockInteraction;
 
-	G4double upperThreshold, lowerThreshold;	
+	G4double upperThreshold, lowerThreshold;
+	G4double triggerEnergy;
 
+	//G4ROOT
+	G4bool factoryOn; 
+ 	G4int  fNtColId[MaxNtCol];
+  //
 
 
 	//G4ThreeVector sourcePosition;
@@ -166,7 +169,7 @@ private:
 	G4double shiftCoeff;
 
 	G4double PositionAngerZ, PositionAngerY;
-
+	
 	//reflector pattern
 	std::vector<G4int> ireflectorLayer1_Tangential;
 	std::vector<G4int> ireflectorLayer1_Axial;
@@ -191,6 +194,11 @@ private:
 	//interaction position with respect to the crystal axis
 	G4ThreeVector interactionPos;
 
+	//interaction position
+	G4double intPosX;
+	G4double intPosY;
+	G4double intPosZ;
+
 
 	G4double interactionTime;
 
@@ -216,10 +224,16 @@ private:
 	G4int DOI_ID0,					DOI_ID1;
 	G4double timeStamp0,			timeStamp1;
 	G4double totalEdep0,			totalEdep1;
-
+	G4double sposX, sposY, sposZ;
 	//choice for the user
 	G4bool getSinglesData;
 	G4bool getCoincidenceData;
+	
+	//
+	G4bool ApplyAngerLogic;
+
+	G4double PMTblurring_tan;
+	G4double PMTblurring_axial;
 
 	G4String outputData;
 	G4int numberOfHit;
@@ -248,6 +262,7 @@ private:
 
 	//G4bool variableResolution;
 	G4bool fixedResolution;
+	G4bool isDOIlookUpTablePrepared;
 
 	G4double energyResolution_fixed;
 	std::vector<std::vector<G4double>> energyResolution_cryDependent;
@@ -282,14 +297,16 @@ private:
 	G4String asciiFileName;
 	G4String rootFileName;
 
-#ifdef USEROOT
-	TTree* tSingles;
-	TTree* tCoincidence;
-	//TH1F*hb;
-#endif
+// #ifdef USEROOT
+//     TFile* file;
+// 	TTree* tSingles;
+// 	TTree* tCoincidence;
+// 	//TH1F*hb;
+// #endif
 
 	//input file to read reflector pattern
 	std::ifstream ifs;
 };
 
 #endif
+

examples/advanced/doiPET/include/doiPETDetectorConstruction.hh

@@ -60,239 +60,241 @@ class doiPETDetectorConstructionMessenger;
 class doiPETDetectorConstruction : public G4VUserDetectorConstruction
 {
 public:
-	doiPETDetectorConstruction();
-	virtual ~doiPETDetectorConstruction();
+  doiPETDetectorConstruction();
+  virtual ~doiPETDetectorConstruction();
 
 public:
-	virtual G4VPhysicalVolume* Construct();
+  virtual G4VPhysicalVolume* Construct();
 
-	void ConstructPhantom(G4LogicalVolume*);
-	void ChangePhantom(G4String choice);
-	void SetPhantomPosition(G4ThreeVector);
-	void SetPhantomRadius(G4double);
-	void SetPhantomLength(G4double);
+  void ConstructPhantom(G4LogicalVolume*);
+  void ChangePhantom(G4String choice);
+  void SetPhantomPosition(G4ThreeVector);
+  void SetPhantomRadius(G4double);
+  void SetPhantomLength(G4double);
+  void SetNumberOfSleeves(G4int);
 
 
 private:
-	void DefineMaterials();
-	doiPETDetectorConstructionMessenger* fDetectorMessenger;
-	doiPETAnalysis* pAnalysis;
+  void DefineMaterials();
+  doiPETDetectorConstructionMessenger* fDetectorMessenger;
+  doiPETAnalysis* pAnalysis;
 
-	G4LogicalVolume* phantom_logicalV;
-	G4VPhysicalVolume* phantom_physicalV;
-	// G4LogicalVolume* gelatin_logicalV;
-	// G4VPhysicalVolume* gelatin_physicalV;
+  G4LogicalVolume* phantom_logicalV;
+  G4VPhysicalVolume* phantom_physicalV;
+  //G4LogicalVolume* gelatin_logicalV;
+  //G4VPhysicalVolume* gelatin_physicalV;
 
-	//wolrd
-	G4LogicalVolume* world_logicalV;
-	G4VPhysicalVolume* world_physicalV;
+  //world
+  G4LogicalVolume* world_logicalV;
+  G4VPhysicalVolume* world_physicalV;
 
-	//detector block
-	G4LogicalVolume* blockDetector_logicalV;
-	G4VPhysicalVolume* blockDetector_physicalV;
+  //detector block
+  G4LogicalVolume* blockDetector_logicalV;
+  G4VPhysicalVolume* blockDetector_physicalV;
 
-	//air volume to fill the detector block
-	G4LogicalVolume* airBox_logicalV;
-	G4VPhysicalVolume* airBox_physicalV;
+  //air volume to fill the detector block
+  G4LogicalVolume* airBox_logicalV;
+  G4VPhysicalVolume* airBox_physicalV;
 
-	//crystals
-	G4LogicalVolume* crystal_logicalV;
-	G4VPhysicalVolume* crystal_physicalV;
+  //crystals
+  G4LogicalVolume* crystal_logicalV;
+  G4VPhysicalVolume* crystal_physicalV;
 
 
-	//water
-	G4LogicalVolume* water_logicalV;
-	G4VPhysicalVolume* water_physicalV;
+  //water
+  G4LogicalVolume* water_logicalV;
+  G4VPhysicalVolume* water_physicalV;
 
-	//lung
-	G4LogicalVolume* lung_logicalV;
-	G4VPhysicalVolume* lung_physicalV;
+  //lung
+  G4LogicalVolume* lung_logicalV;
+  G4VPhysicalVolume* lung_physicalV;
 
-	//test line phantom
-	G4LogicalVolume* test_logicalV;
-	G4VPhysicalVolume* test_physicalV;
+  //test line phantom
+  G4LogicalVolume* test_logicalV;
+  G4VPhysicalVolume* test_physicalV;
 
-	G4LogicalVolume* lung_logicalV_PMMA;
-	G4VPhysicalVolume* lung_physicalVPMMA;
+  G4LogicalVolume* lung_logicalV_PMMA;
+  G4VPhysicalVolume* lung_physicalVPMMA;
 
-	//cold regions
-	// G4LogicalVolume* coldRegion_logicalV;
-	// G4VPhysicalVolume* coldRegion_physicalV;
+  //cold regions
+  //G4LogicalVolume* coldRegion_logicalV;
+  //G4VPhysicalVolume* coldRegion_physicalV;
 
-	//
-	//Surrounding PMMA for hot sphere
-	G4LogicalVolume* hotSpherePMMA_logicalV;
-	G4VPhysicalVolume* hotSpherePMMA_physicalV;
+  //
+  //Surrounding PMMA for hot sphere
+  G4LogicalVolume* hotSpherePMMA_logicalV;
+  G4VPhysicalVolume* hotSpherePMMA_physicalV;
 
-	//hot water phantom (activity is distributed)
-	G4LogicalVolume* hotSphereWater_logicalV;
-	G4VPhysicalVolume* hotSphereWater_physicalV;
+  //hot water phantom (activity is distributed)
+  G4LogicalVolume* hotSphereWater_logicalV;
+  G4VPhysicalVolume* hotSphereWater_physicalV;
 
-	//surrounding PMMA cold sphere
-	G4LogicalVolume* coldSpherePMMA_logicalV;
-	G4VPhysicalVolume* coldSpherePMMA_physicalV;
+  //surrounding PMMA cold sphere
+  G4LogicalVolume* coldSpherePMMA_logicalV;
+  G4VPhysicalVolume* coldSpherePMMA_physicalV;
 
-	//cold Water phantom in the cold PMMA sphere 
-	G4LogicalVolume* coldSphereWater_logicalV;
-	G4VPhysicalVolume* coldSphereWater_physicalV;
+  //cold Water phantom in the cold PMMA sphere 
+  G4LogicalVolume* coldSphereWater_logicalV;
+  G4VPhysicalVolume* coldSphereWater_physicalV;
 
-	//fillable polyethylene phantom for sensitivity
-	G4LogicalVolume* phantomPE_logicalV;
-	G4VPhysicalVolume* phantomPE_physicalV;
+  //fillable polyethylene phantom for sensitivity
+  G4LogicalVolume* phantomPE_logicalV;
+  G4VPhysicalVolume* phantomPE_physicalV;
 
-	//Image quality phantom for small animal NEMA NU-4
-	G4LogicalVolume* waterPhantom_logicalV;
-	G4VPhysicalVolume* WaterPhantom_physicalV;
+  //Image quality phantom for small animal NEMA NU-4
+  G4LogicalVolume* waterPhantom_logicalV;
+  G4VPhysicalVolume* WaterPhantom_physicalV;
 
-	G4LogicalVolume* rod_phantom_logicalV;
-	G4VPhysicalVolume* rod_phantom_physicalV;
+  G4LogicalVolume* rod_phantom_logicalV;
+  G4VPhysicalVolume* rod_phantom_physicalV;
 
-	G4LogicalVolume* chamberPMMA_logicalV;
-	G4VPhysicalVolume* chamberPMMA_physicalV;
+  G4LogicalVolume* chamberPMMA_logicalV;
+  G4VPhysicalVolume* chamberPMMA_physicalV;
 
-	//
-	G4LogicalVolume* chamberWater_logicalV;
-	G4VPhysicalVolume* chamberWater_physicalV;
+  //
+  G4LogicalVolume* chamberWater_logicalV;
+  G4VPhysicalVolume* chamberWater_physicalV;
 
-	//
-	G4LogicalVolume* chamberAir_logicalV;
-	G4VPhysicalVolume* chamberAir_physicalV;
+  //
+  G4LogicalVolume* chamberAir_logicalV;
+  G4VPhysicalVolume* chamberAir_physicalV;
 
 
-	//Dimension of the sphere
-	G4double spherePositionX, spherePositionY; // spherePositionZ;
-	G4double sphereDiameter;
-	G4double distanceFromCenter;
-	G4int numberOfSpheres;
-	G4double sphereWallThickness;
-	G4double zOffsetSpherePhantom;
+  //Dimension of the sphere
+  G4double spherePositionX, spherePositionY; //, spherePositionZ;
+  G4double sphereDiameter;
+  G4double distanceFromCenter;
+  G4int numberOfSpheres;
+  G4double sphereWallThickness;
+  G4double zOffsetSpherePhantom;
 
-	G4String PhantomType;
+  G4String PhantomType;
 
-	//materials
-	G4Material* air;
-	G4Material* pmma;
-	G4Material* water;
-	G4Material* polyethylene;
-	G4Material* polyethylene_NEMA;
-	// G4Material* inflatedLung;
-	G4Material* polystyrene;
-	G4Material* Aluminum;
+  //materials
+  G4Material* air;
+  G4Material* pmma;
+  G4Material* water;
+  G4Material* polyethylene;
+  G4Material* polyethylene_NEMA;
+  //G4Material* inflatedLung;
+  G4Material* polystyrene;
+  G4Material* Aluminum;
 
-	//elements for GSO
-	G4Element*  O;
-	G4Element* Si;
-	G4Element* Gd;
-	G4Material* GSO;
+  //elements for GSO
+  G4Element*  O;
+  G4Element* Si;
+  G4Element* Gd;
+  G4Material* GSO;
 
-	G4Material* crystalMaterial;
-	//G4Material* phantomMaterial;
+  G4Material* crystalMaterial;
+  //G4Material* phantomMaterial;
 
-	G4bool  fCheckOverlaps;
-	G4bool isotopes;
+  G4bool  fCheckOverlaps;
+  G4bool isotopes;
 
 
-	//size of world
-	G4double worldSizeX;
-	G4double worldSizeY;
-	G4double worldSizeZ;
+  //size of world
+  G4double worldSizeX;
+  G4double worldSizeY;
+  G4double worldSizeZ;
 
 
-	//The following is moved to doiPETGlobalParameters.hh
-	//G4int numberOfCrystal_DOI;
-	//G4int numberOfCrystal_tangential;
-	//G4int numberOfCrystal_axial;
+  //The following is moved to doiPETGlobalParameters.hh
+  //G4int numberOfCrystal_DOI;
+  //G4int numberOfCrystal_tangential;
+  //G4int numberOfCrystal_axial;
 
-	////
-	//G4double sizeOfCrystal_DOI; 
-	//G4double sizeOfCrystal_tangential;
-	//G4double sizeOfCrystal_axial;
+  ////
+  //G4double sizeOfCrystal_DOI; 
+  //G4double sizeOfCrystal_tangential;
+  //G4double sizeOfCrystal_axial;
 
-	////
-	//G4double crystalGap_DOI;
-	//G4double crystalGap_tangential;
-	//G4double crystalGap_axial;
+  ////
+  //G4double crystalGap_DOI;
+  //G4double crystalGap_tangential;
+  //G4double crystalGap_axial;
 
-	G4double sizeOfAirBox_DOI;
-	G4double sizeOfAirBox_axial;
-	G4double sizeOfAirBox_tangential;
+  G4double sizeOfAirBox_DOI;
+  G4double sizeOfAirBox_axial;
+  G4double sizeOfAirBox_tangential;
 
 
-	G4double sizeOfBlockDetector_DOI;
-	G4double sizeOfBlockDetector_axial;
-	G4double sizeOfBlockDetector_tangential;
+  G4double sizeOfBlockDetector_DOI;
+  G4double sizeOfBlockDetector_axial;
+  G4double sizeOfBlockDetector_tangential;
 
-	//G4double AluminumCoverThickness;
+  //G4double AluminumCoverThickness;
 
 
-	//G4int numberOfPETDetector;
-	//G4int numberOfRings;
+  //G4int numberOfPETDetector;
+  //G4int numberOfRings;
 
 
-	//G4double scannerRadius;
-	G4double thetaDetector; //The azimuthal angle for arranging the detector in the PET ring 
-	//G4double ringGap;
-	G4int blockIndex;
-	// G4int AlCase_Index;
-	G4int crystalIndex;
+  //G4double scannerRadius;
+  G4double thetaDetector; //The azimuthal angle for arranging the detector in the PET ring 
+  //G4double ringGap;
+  G4int blockIndex;
+  //G4int AlCase_Index;
+  G4int crystalIndex;
 
-	//detector position
-	G4double detectorPositionX;
-	G4double detectorPositionY;
-	G4double detectorPositionZ;
+  //detector position
+  G4double detectorPositionX;
+  G4double detectorPositionY;
+  G4double detectorPositionZ;
 
-	//crystal position
-	G4double crystalPositionX; 
-	G4double crystalPositionY;
-	G4double crystalPositionZ;
+  //crystal position
+  G4double crystalPositionX; 
+  G4double crystalPositionY;
+  G4double crystalPositionZ;
 
 
 
-	G4ThreeVector phantomPosition;
+  G4ThreeVector phantomPosition;
 
-	//
-	G4double phantomRadius;
-	G4double phantomLength;
+  //
+  G4double phantomRadius;
+  G4double phantomLength;
 
-	//Phantom dimension for rectangular box (placed for therapy study)
-	// G4double phantomSizeX, phantomSizeY, phantomSizeZ;
+  //Phantom dimension for rectangular box (placed for therapy study)
+  //G4double phantomSizeX, phantomSizeY, phantomSizeZ;
 
-	//the following is to make the body phantom 
-	G4double yOffsetBodyPhantom;
-	G4double zOffsetBodyPhantom;
-	G4double lengthOfBodyPhantom; //Interior length ( = 180m mm) + wallthickness (= 3mm)
-	G4double radiusOfBodyPhantom;
-	G4double wallThicknessOfBodyPhantom;
-	G4double radiusOfLungPhantom;
+  //the following is to make the body phantom 
+  G4double yOffsetBodyPhantom;
+  G4double zOffsetBodyPhantom;
+  G4double lengthOfBodyPhantom; //Interior length ( = 180m mm) + wallthickness (= 3mm)
+  G4double radiusOfBodyPhantom;
+  G4double wallThicknessOfBodyPhantom;
+  G4double radiusOfLungPhantom;
 
-	//Test phantom defnition. The phantom has the same as that that of NECR phantom except
-	G4double hieghtOfTestPhantom; 
-	G4double diameterOfTestPhantom;
+  //Test phantom defnition. The phantom has the same as that that of NECR phantom except
+  G4double hieghtOfTestPhantom; 
+  G4double diameterOfTestPhantom;
 
-	//To the cylindrical phantom to make the body phantom
-	G4double radiusOfSmallcyl;
-	G4double boxWidth;
-	G4double boxHeight;
+  //To the cylindrical phantom to make the body phantom
+  G4double radiusOfSmallcyl;
+  G4double boxWidth;
+  G4double boxHeight;
 
 
-	//Image quality phantom for small animals
-	G4double waterPhantomRadius;
-	G4double waterPhantomLength;
+  //Image quality phantom for small animals
+  G4double waterPhantomRadius;
+  G4double waterPhantomLength;
 
-	G4double rodPhantomLength;
-	G4double rodDiameter;
-	G4int numberOfRods;
+  G4double rodPhantomLength;
+  G4double rodDiameter;
+  G4int numberOfRods;
 
-	//Declare position for the rod phantoms
-	G4double rodPositionX, rodPositionY, rodPositionZ;
-
-	//Declare position for cold region chanmbers
-	G4double chamberPositionX, chamberPositionY, chamberPositionZ;
-	G4double chamberPhantomLength;
-	G4double chamberDiameter;
-	G4double wallThicknessOfChamber;
+  //Declare position for the rod phantoms
+  G4double rodPositionX, rodPositionY, rodPositionZ;
 
+  //Declare position for cold region chanmbers
+  G4double chamberPositionX, chamberPositionY, chamberPositionZ;
+  G4double chamberPhantomLength;
+  G4double chamberDiameter;
+  G4double wallThicknessOfChamber;
 
+  //Declare the number of concentric (Al) sleeves
+  G4int numOfSleeves;
 };
 
 /////////////////////////////////////////////////////////////////////////////////////////

examples/advanced/doiPET/include/doiPETDetectorConstructionMessenger.hh

@@ -75,6 +75,7 @@ private:
 	G4UIcmdWith3VectorAndUnit *changeThePhantomPositionCmd; 
 	G4UIcmdWithADoubleAndUnit* changePhantomRadiusCmd;
 	G4UIcmdWithADoubleAndUnit* changePhantomLengthCmd;
+	G4UIcmdWithAnInteger* changeNumOfSleevesCmd;
 };
 
 #endif

examples/advanced/doiPET/include/doiPETEventAction.hh

@@ -51,15 +51,15 @@
 class doiPETEventAction : public G4UserEventAction
 {
 public:
-	doiPETEventAction();
-	virtual ~doiPETEventAction();
+  doiPETEventAction();
+  virtual ~doiPETEventAction();
 
-	virtual void  BeginOfEventAction(const G4Event*);
-	virtual void  EndOfEventAction(const G4Event*);
+  virtual void  BeginOfEventAction(const G4Event*);
+  virtual void  EndOfEventAction(const G4Event*);
 
 private:
-	// G4int eventID;
-	G4int printModulo;
+  //G4int eventID;
+  G4int printModulo;
 };
 
 #endif

examples/advanced/doiPET/include/doiPETRun.hh

@@ -26,6 +26,7 @@
 /// \file electromagnetic/TestEm11/include/doiPETRun.hh
 /// \brief Definition of the doiPETRun class
 //
+// $Id: doiPETRun.hh 71375 2013-06-14 07:39:33Z maire $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -66,24 +67,28 @@ public:
 	void SetAnnihilationTime(G4double);
 	void SetEventID(G4int);
 	virtual void Merge(const G4Run*);
-	// void EndOfRun();     
+	G4ThreeVector CenterOfMassInteractionPos(const std::vector<G4int>&, const std::vector<G4double>&, G4double, const std::vector<G4ThreeVector>&);
+	//void EndOfRun();     
 
 private:
 
 	std::multimap< G4int, InteractionInformation* > mapBlockInteraction;
 	std::set<G4int> setBlockInteraction;
 	std::ofstream ofs;
+
+	std::vector<G4int> crystalID_vec;
+	std::vector<G4ThreeVector> posInter_vec;
+	std::vector<G4double>edepInCry_vec;
+
+	//center of mass interaction of position
+	G4ThreeVector posInterInCrystal;
+
 	G4int eventID;
 	doiPETRunAction* fRunAction;
 	doiPETAnalysis* fAnalysis;
 	G4double totalEdep;
 	G4int blockID, crystalID;
 
-	//
-	G4double activityNow;
-	G4double InitialActivity;
-	G4double halfLife;
-
 	G4double totalTime;
 	G4double prev_totalTime;
 	G4double timeInterval;

examples/advanced/doiPET/include/doiPETRunAction.hh

@@ -43,6 +43,7 @@
 
 #include "G4UserRunAction.hh"
 #include "globals.hh"
+#include "doiPETAnalysis.hh"
 
 class doiPETDetectorConstruction;
 class doiPETRun;
@@ -54,8 +55,13 @@ class doiPETPrimaryGeneratorAction;
 class doiPETRunAction : public G4UserRunAction
 {
 public:
+#ifdef ANALYSIS_USE
+    doiPETRunAction(doiPETAnalysis* analysis);
+#else
+    doiPETRunAction();
+#endif
 	//doiPETRunAction();
-	doiPETRunAction(/*doiPETDetectorConstruction*, doiPETPrimaryGeneratorAction**/);
+	//doiPETRunAction(/*doiPETDetectorConstruction*, doiPETPrimaryGeneratorAction**/);
 	virtual ~doiPETRunAction();
 	virtual G4Run* GenerateRun();  
 
@@ -63,8 +69,10 @@ public:
 	virtual void   EndOfRunAction(const G4Run*);
 
 private:
-	doiPETRun*                       fRun;    
-	//HistoManager*              fHistoManager;
+	doiPETRun* fRun;    
+#ifdef ANALYSIS_USE
+	doiPETAnalysis* analysisMan;
+#endif
 
 };
 #endif

examples/advanced/doiPET/inputParameter.txt

@@ -2,7 +2,7 @@
 #For theoretical look-up table (equally spaced position response), the out most reflector value is set to be 0
 
 #Reflector pattern for layer 1
-reflectorLayer1_Tangential:	1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
+reflectorLayer1_Tangential:	1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 
 reflectorLayer1_Axial:	1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
 #
 #Reflector pattern for layer 2
@@ -32,15 +32,22 @@ crystalResolutionMax: 0.17
 #Reference energy in keV
 crystalEnergyRef: 511 keV
 
-crystalQuantumEfficiency: 0.905
+crystalQuantumEfficiency: 1.0
 
 #The number of pixels in 2D position histogram of the response
-numberOfPixel_2D_Pixel: 32 32
 
 #Energy window in keV
 lowerThreshold: 400 keV
-upperThreshold: 600 keV 
+upperThreshold: 600 keV
+triggerEnergy:  0.05 keV # equals to 50 eV
+ 
+#Apply Anger logic calculation
+ApplyAngerLogic: true
+
+#2D position blurring (PMT position response blurring) (tangential(y) axial(z))
+PMTblurring: 1.5 1.5
+
+#Choose the type of output: singlesOutput or coincidenceOutput
+TypeOfOutput: singlesOutput
 
-#Choose the type of output: "singlesOutput" or "coincidenceOutput"
-TypeOfOutput: coincidenceOutput
 

examples/advanced/doiPET/run.mac

@@ -1,21 +1,26 @@
 #==========================- Pre-initialisation parameters -=========================
 
 #choose the type of phantom
-/changePhantom/setPhantom Phantom_spatialResolution 
+/changePhantom/setPhantom Phantom_sensitivity 
 
 
 #set the activity in the phantom
-/change/Activity 3840000 becquerel #3840000 becquerel
+/change/Activity 500000 becquerel #3840000 becquerel
 
 #Set the halflife of the isotope in units of second
 /change/HalfLife 6586.2 s #Halflife of F-18
 
-#Change the size of the phantom
-/changePhantom/Radius 2 mm 
-/changePhantom/Length 4 mm
 
-#Change the position of phantom. !!!Remember to change the source position as well.
-/changePhantom/position 10 0 0 mm
+#change the number of concentric sleeves
+# The diameter of the sleeves are changed based on the NEMA NU2 guideline
+#sleeve(minDiameter, maxDiameter) = sleeve1(3.9, 6.4)mm, sleeve2(7.0, 9.5)mm, sleeve3(10.2, 12.7)mm, sleeve4(13.4, 15.9)mm, sleeve5(16.6, 19.1)mm
+/changePhantom/NumberOfSleeves 1 #can be changed up to 5
+
+#Change the size of the phantom
+/changePhantom/Length 216 mm
+
+#Remember to change the source position along with the phantom position 
+/changePhantom/position 0 0 0 mm
 
 /run/initialize
 /run/verbose 2
@@ -25,12 +30,14 @@
 #/gps/particle e+
 /gps/energy 0
 /gps/pos/type Volume
-/gps/pos/shape Cylinder 
-/gps/pos/centre 10 0 0 mm
-/gps/pos/radius 0.5 mm
-/gps/pos/halfz  0.5 mm
+/gps/pos/shape Cylinder
+/gps/pos/centre 0 0 0 mm 
+/gps/pos/radius 0.5 mm 
+/gps/pos/halfz  108 mm
 /gps/ang/type iso
 
 #/random/setSeeds  1 2
 
-/run/beamOn 100000 #36461089
+/run/beamOn 1000000
+
+

examples/advanced/doiPET/run_imageQualityPhantom_wholeBody.mac

@@ -68,9 +68,9 @@
 /gps/pos/type Volume
 /gps/pos/shape Cylinder
 
-#The body phantom is shifted by 68 mm, 
-/gps/pos/centre 0 0 -23 mm 
-/gps/pos/radius 187 mm #increase the diameter
+#The body phantom is shifted as discribed in the NEMA standard
+/gps/pos/centre 0 -35 -23 mm 
+/gps/pos/radius 147 mm
 /gps/pos/halfz  88.5 mm
 /gps/ang/type iso
 
@@ -93,7 +93,7 @@
 
 #/random/setSeeds  1 2 
 
-/run/beamOn 836431791 #15 minuts scan
+/run/beamOn 10000
 
 
 

examples/advanced/doiPET/run_sensitivity.mac

@@ -3,15 +3,20 @@
 #choose the type of phantom
 /changePhantom/setPhantom Phantom_sensitivity 
 
+
 #set the activity in the phantom
-/change/Activity 3840000 becquerel #3840000 becquerel
+/change/Activity 500000 becquerel #3840000 becquerel
 
 #Set the halflife of the isotope in units of second
 /change/HalfLife 6586.2 s #Halflife of F-18
 
 
+#change the number of concentric sleeves
+# The diameter of the sleeves are changed based on the NEMA NU2 guideline
+#sleeve(minDiameter, maxDiameter) = sleeve1(3.9, 6.4)mm, sleeve2(7.0, 9.5)mm, sleeve3(10.2, 12.7)mm, sleeve4(13.4, 15.9)mm, sleeve5(16.6, 19.1)mm
+/changePhantom/NumberOfSleeves 5
+
 #Change the size of the phantom
-/changePhantom/Radius 4.75 mm # select 3.2, 4.75,6.35,7.95,9.55 mm
 /changePhantom/Length 705 mm
 
 #Remember to change the source position along with the phantom position 
@@ -20,9 +25,9 @@
 /run/initialize
 /run/verbose 2
 
-#/gps/particle ion
-#/gps/ion 9 18 0 0 # F-18
-/gps/particle e+
+/gps/particle ion
+/gps/ion 9 18 0 0 # F-18
+#/gps/particle e+
 /gps/energy 0
 /gps/pos/type Volume
 /gps/pos/shape Cylinder
@@ -34,3 +39,5 @@
 #/random/setSeeds  1 2
 
 /run/beamOn 1000000
+
+

examples/advanced/doiPET/run_spatialResolution.mac

@@ -31,4 +31,4 @@
 
 #/random/setSeeds  1 2
 
-/run/beamOn 15000000 #36461089
+/run/beamOn 100000 #36461089

examples/advanced/doiPET/src/doiPETActionInitialization.cc

@@ -49,8 +49,10 @@
 doiPETActionInitialization::doiPETActionInitialization()
 : G4VUserActionInitialization()
 {}*/
-doiPETActionInitialization::doiPETActionInitialization() : G4VUserActionInitialization()
-{}
+doiPETActionInitialization::doiPETActionInitialization(doiPETAnalysis* analysisMan) : G4VUserActionInitialization()
+{
+	analysis = analysisMan;
+}
 
 doiPETActionInitialization::~doiPETActionInitialization()
 {}
@@ -59,7 +61,11 @@ doiPETActionInitialization::~doiPETActionInitialization()
 //
 void doiPETActionInitialization::BuildForMaster() const
 {
-	doiPETRunAction* runAction = new doiPETRunAction();
+#ifdef ANALYSIS_USE
+	 doiPETRunAction* runAction = new doiPETRunAction(analysis);
+#else
+ 	doiPETRunAction* runAction = new doiPETRunAction();
+#endif
 	SetUserAction(runAction);
 }
 
@@ -69,7 +75,11 @@ void doiPETActionInitialization::Build() const
 	doiPETPrimaryGeneratorAction* primary = new doiPETPrimaryGeneratorAction();
 	SetUserAction(primary);
 
-	doiPETRunAction* runAction = new doiPETRunAction();
+#ifdef ANALYSIS_USE
+	 doiPETRunAction* runAction = new doiPETRunAction(analysis);
+#else
+ 	doiPETRunAction* runAction = new doiPETRunAction();
+#endif
 	SetUserAction(runAction);
 
 	SetUserAction(new doiPETEventAction());  

examples/advanced/doiPET/src/doiPETAnalysis.cc

@@ -57,20 +57,25 @@
 #include "Randomize.hh"
 #include "G4SPSRandomGenerator.hh"
 #include "doiPETAnalysisMessenger.hh"
-
+#include "G4UnitsTable.hh"
+#include "globals.hh"
 
 doiPETAnalysis* doiPETAnalysis::instance=0;
 
 /////////// Constructor /////////////////////////////////////////////
 doiPETAnalysis::doiPETAnalysis()
 {
+	factoryOn = false; //G4ROOT
+	// Initialization ntuple
+  	for (G4int k=0; k<MaxNtCol; k++) fNtColId[k] = 0;
+
 	fAnalysisMessenger = new doiPETAnalysisMessenger(this);
 
 	//Set energy window
 	lowerThreshold = 400*keV;
 	upperThreshold = 600*keV;
-
-
+	triggerEnergy = 50*eV;
+	
 	//give default initial activity. Activity strength is changed in the .mac file
 	InitialActivity = 1000000*becquerel;
 
@@ -92,8 +97,8 @@ doiPETAnalysis::doiPETAnalysis()
 	scatterIndex = 0;
 
 	//
-	numberOfPixel_tan = 32;
-	numberOfPixel_axial = 32;
+	numberOfPixel_tan = 2*numberOfCrystal_tangential; //32;
+	numberOfPixel_axial = 2*numberOfCrystal_axial; //32;
 
 	//Default value for deadtime.
 	block_DeadTime = 256*ns;
@@ -109,7 +114,7 @@ doiPETAnalysis::doiPETAnalysis()
 	crystalEnergyRef = 511 * keV;//Energy of reference in which the energy resolution of the crystal is computed
 
 	//The quantum efficiency models the probability for the event to be detected by the photo-detector.
-	//The quantum efficiency can be set inputParameter.txt file
+	//The quantum efficiency can be set in the inputParameter.txt file
 	crystalQuantumEfficiency = 1;//100% 
 	//
 
@@ -127,10 +132,18 @@ doiPETAnalysis::doiPETAnalysis()
 	//Initialize type of output. The default output is single events
 	getSinglesData  = false; //default value 
 	getCoincidenceData = false;
+
+	ApplyAngerLogic = true;
+	isDOIlookUpTablePrepared = false;
+
 	numberOfHit = 0;
 
-	//This value is based on the assumption that the shift due to the reflector is half distance from the interaction position to the air gap.
+	//This value is based on the assumption that the shift of the response due to the reflector is half distance from the interaction position to the air gap.
 	shiftCoeff = 0.5;
+
+	//
+	PMTblurring_tan = 0.0;
+	PMTblurring_axial = 0.0;
 }
 ////////// Destructor ///////////////////////////////////////////////
 doiPETAnalysis::~doiPETAnalysis()
@@ -153,9 +166,8 @@ void doiPETAnalysis::Delete()
 
 //If there is energy deposition in the phantom by the photon, the scatter index is 1, otherwise it is 0
 //Use this for checking
-void doiPETAnalysis::GetScatterIndexInPhantom(G4double edepInPhantom){
-	if(edepInPhantom>0)scatterIndex = 1;
-	else scatterIndex = 0;
+void doiPETAnalysis::SetScatterIndexInPhantom(G4int sci){
+	scatterIndex = sci;
 }
 
 //Get the source position if the process is annihilation.
@@ -174,6 +186,7 @@ void doiPETAnalysis::SetEventID(G4int evID){
 }
 
 //
+
 void doiPETAnalysis::GetSizeOfDetector(G4double detSizeDoi, G4double detSizeTan, G4double detSizeAxial){
 	sizeOfDetector_DOI = detSizeDoi;
 	sizeOfDetector_axial = detSizeTan;
@@ -244,6 +257,8 @@ void doiPETAnalysis::ReadOut(G4int blkID, G4int cryID, G4double interTime, G4dou
 	//time of the event when detected (timerTag)
 	timeStamp = totalTime + time_tof;
 
+	//triggerEnergy is the energy deposited in the detector below which the detector is insensitive to any interaction.
+	if(totalEdep<triggerEnergy)return;
 
 	//************************************** Apply dead-time ********************************************//
 	//Apply paralizable dead-time in the block beofore events are rejected by the energy window
@@ -291,15 +306,15 @@ void doiPETAnalysis::ReadOut(G4int blkID, G4int cryID, G4double interTime, G4dou
 		crystalID_axial = crystalID_2D/numberOfCrystal_axial;
 		crystalID_tangential = crystalID_2D%numberOfCrystal_tangential;
 
-		//Calculate local position of the crystal with respect to the detector. Only the lateral distances (tangential (y) and axial (z) are needed.)
-		//G4double posCrystalX = (DOI_ID-((G4double)numberOfCrystal_DOI)/2 + 0.5)*(sizeOfCrystal_DOI + crystalGap_DOI) + interactionPos.x();
-		G4double posCrystalY = (crystalID_tangential-((G4double)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential) + interactionPos.y();
-		G4double posCrystalZ = (crystalID_axial-((G4double)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial) + interactionPos.z();
+		intPosX = interactionPos.x();
+		intPosY = interactionPos.y();
+		intPosZ = interactionPos.z();
 
 
-		//shiftCoeff = 0.5 is used. This value is based on the assumption that the shift due to the reflector is half distance from the interaction position to the air gap. 
-		AngerLogic(DOI_ID, crystalID_tangential, crystalID_axial, posCrystalY, posCrystalZ, totalEdep, shiftCoeff);//
-		
+		if(ApplyAngerLogic){
+			//shiftCoeff = 0.5 is used. This value is based on the assumption that the shift due to the reflector is half distance from the interaction position to the air gap. 
+			AngerLogic(intPosX, intPosY, intPosZ, totalEdep, shiftCoeff, isDOIlookUpTablePrepared);//
+		}
 
 		//Single event output. Coincidence events can then be made using the single events.
 		if(getSinglesData) WriteOutput();
@@ -345,9 +360,11 @@ void doiPETAnalysis::Open(G4String fileName)
 {
 	if(getSinglesData){
 		asciiFileName = fileName + "Singles.data";
+		rootFileName = fileName+"Singles.root";
 	}
 	if(getCoincidenceData){
 		asciiFileName = fileName + "Coincidence.data";
+		rootFileName = fileName+"Coincidence.root";
 	}
 
 	ofs.open(asciiFileName.c_str());
@@ -355,51 +372,14 @@ void doiPETAnalysis::Open(G4String fileName)
 		G4cerr<<"=== \n File opening Error to write the output ===="<<G4endl;
 		exit(0);
 	}
-	//
-#ifdef USEROOT
-	if(getSinglesData){ 
-		rootFileName = fileName+"Singles.root";
-		tSingles = new TTree("tSingles","SinglesTree");
-		tSingles->Branch("eventID",&eventID,"eventID/I");
-		tSingles->Branch("blockID",&blockID,"blockID/I");
-		tSingles->Branch("crystalID_axial",&crystalID_axial,"crystalID_axial/I");
-		tSingles->Branch("crystalID_tangential",&crystalID_axial,"crystalID_tangential/I");
-		tSingles->Branch("DOI_ID",&DOI_ID0,"DOI_ID/I");
-		tSingles->Branch("timeStamp",&timeStamp,"timeStamp/D");
-		tSingles->Branch("totalEdep",&totalEdep,"totalEdep/D");
-	}
 
-	if(getCoincidenceData){
-		rootFileName = fileName+"Coincidence.root";
-		tCoincidence = new TTree("tCoincidence","CoincidenceTree");
-		//First Single
-		tCoincidence->Branch("eventID0",&eventID0,"eventID0/I");
-		tCoincidence->Branch("blockID0",&blockID0,"blockID0/I");
-		tCoincidence->Branch("crystalID_axial0",&crystalID_axial0,"crystalID_axial0/I");
-		tCoincidence->Branch("crystalID_tangential0",&crystalID_axial0,"crystalID_tangential0/I");
-		tCoincidence->Branch("DOI_ID0",&DOI_ID0,"DOI_ID0/I");
-		tCoincidence->Branch("timeStamp0",&timeStamp0,"timeStamp0/D");
-		tCoincidence->Branch("totalEdep0",&totalEdep0,"totalEdep0/D");
-
-		//Second Single
-		tCoincidence->Branch("eventID1",&eventID1,"eventID1/I");
-		tCoincidence->Branch("blockID1",&blockID1,"blockID1/I");
-		tCoincidence->Branch("crystalID_axial1",&crystalID_axial1,"crystalID_axial1/I");
-		tCoincidence->Branch("crystalID_tangential1",&crystalID_axial1,"crystalID_tangential1/I");
-		tCoincidence->Branch("DOI_ID1",&DOI_ID1,"DOI_ID1/I");
-		tCoincidence->Branch("timeStamp1",&timeStamp1,"timeStamp1/D");
-		tCoincidence->Branch("totalEdep1",&totalEdep1,"totalEdep1/D");
-	}
-#endif
-	//
 }
 
 void doiPETAnalysis::WriteOutput(){
 	if(getSinglesData){
-		ofs<<eventID<<" "<<blockID<<" "<<crystalID_axial<<" "<<crystalID_tangential<<" "<<DOI_ID<<" "<<std::setprecision(17)<<timeStamp/s<<" "<<std::setprecision(7)<<totalEdep/keV<<G4endl;
-#ifdef USEROOT
-		tSingles->Fill();
-#endif
+		ofs<<eventID<<" "<<blockID<<" "<<std::setprecision(17)<<timeStamp/s<<" "<<std::setprecision(7)<<totalEdep/keV<<" "<<intPosX<<" "<<intPosY<<" "<<intPosZ<<" "<<spositionX<<" "<<spositionY<<" "<<spositionZ<<G4endl;
+		//ofs<<eventID<<" "<<blockID<<" "<<crystalID_axial<<" "<<crystalID_tangential<<" "<<DOI_ID<<" "<<std::setprecision(17)<<timeStamp/s<<" "<<std::setprecision(7)<<totalEdep/keV<<G4endl;
+
 	}
 	if(getCoincidenceData){
 		//2 singles will qualify to be in coincidence within the energy window.
@@ -428,12 +408,13 @@ void doiPETAnalysis::WriteOutput(){
 		}
 
 		ofs<<eventID0<<" "<<blockID0<<" "<<crystalID_axial0<<" "<<crystalID_tangential0<<" "<<DOI_ID0<<" "<<std::setprecision(17)<<timeStamp0/s<<" "<<std::setprecision(7)<<totalEdep0/keV<<" "
-			<<eventID1<<" "<<blockID1<<" "<<crystalID_axial1<<" "<<crystalID_tangential1<<" "<<DOI_ID1<<" "<<std::setprecision(17)<<timeStamp1/s<<" "<<std::setprecision(7)<<totalEdep1/keV<<G4endl;
+			<<eventID1<<" "<<blockID1<<" "<<crystalID_axial1<<" "<<crystalID_tangential1<<" "<<DOI_ID1<<" "<<std::setprecision(17)<<timeStamp1/s<<" "<<std::setprecision(7)<<totalEdep1/keV<<" "
+			<<spositionX<<" "<<spositionY<<" "<<spositionZ<<G4endl;
 
-#ifdef USEROOT
-		tCoincidence->Fill();
-#endif
 	}
+#ifdef ANALYSIS_USE
+	FillListModeEvent();
+#endif
 
 }
 
@@ -444,33 +425,13 @@ void doiPETAnalysis::Close()
 	//close ascii file
 	ofs.close();
 
-	//
-#ifdef USEROOT
-	TFile f(rootFileName.c_str(),"RECREATE");
-	if(getSinglesData){
-		tSingles->Write();
-		delete tSingles;
-	}
-	if(getCoincidenceData){
-		tCoincidence->Write();
-		delete tCoincidence;
-	}
-
-	f.Close();
-
-
-#endif
-
 }
 
-//Place the photomultiplier tube (PMT) at each corner of the detector. 
-//The positions of the PMT is with respect to the axis of the detector block
-//All the PMTs are placed at the same doi (x) position 
-//(at +sizeOfDetector_DOI/2 which is at the top of the detector). 
+//Place the photomultiplier tube (PMT) at each corner of the detector. The positions of the PMT is with respect to the axis of the detector block
+//All the PMTs are placed at the same doi (x) position (at +sizeOfDetector_DOI/2 which is at the top of the detector). 
 
 //The PMT is placed at each corner of the crystal block and is assumed to be an ideal PMT.
-//The signal (energy deposition) of each PMT depends on  the distance of the respective 
-// PMT from the interaction point
+//The signal (energy deposition) of each PMT depends on  the distance of the respective PMT from the interaction point
 void doiPETAnalysis::PMTPosition(){
 
 	sizeOfDetector_DOI = (numberOfCrystal_DOI * sizeOfCrystal_DOI) + (numberOfCrystal_DOI - 1)*crystalGap_DOI;
@@ -568,8 +529,8 @@ void doiPETAnalysis::BlurringParameters(){
 				else {
 					fixedResolution = false;
 					//Store into a file if needed.
-					//std::string fname = "crystalDependentResolution.txt";
-					//std::ofstream outFname(fname.c_str());
+					std::string fname = "crystalDependentResolution.txt";
+					std::ofstream outFname(fname.c_str());
 
 					G4cout<<" \n Crystal dependent resolution is used. preparing look-up table .... "<<G4endl;
 					energyResolution_cryDependent.resize(numberOfBlocks_total,std::vector<G4double>(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0));	
@@ -577,11 +538,11 @@ void doiPETAnalysis::BlurringParameters(){
 						for(G4int i_cry = 0; i_cry < numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI; i_cry++){
 							energyResolution_cryDependent[i_blk][i_cry] = crystalResolutionMin + (crystalResolutionMax - crystalResolutionMin)*G4UniformRand();
 							//store into a file
-							//outFname<<i_blk<<" "<<i_cry<<" "<<energyResolution_cryDependent[i_blk][i_cry]<<G4endl;
+							outFname<<i_blk<<" "<<i_cry<<" "<<energyResolution_cryDependent[i_blk][i_cry]<<G4endl;
 						}
 					}
 					G4cout<<"Done. \n"<<G4endl;
-					//outFname.close();
+					outFname.close();
 				}
 
 			}
@@ -629,12 +590,45 @@ void doiPETAnalysis::BlurringParameters(){
 
 			}
 
-			if( (std::string::size_type)inputLine.find("numberOfPixel_2D_Pixel:")!=std::string::npos){
+			//
+			if( (std::string::size_type)inputLine.find("triggerEnergy:")!=std::string::npos){
 				std::istringstream tmpStream(inputLine);
 				tmpStream >> value[0] >> value[1] >> value[2];
-				numberOfPixel_axial = atof(value[1].c_str());
-				numberOfPixel_tan = atof(value[2].c_str());
-				G4cout<<"Number of pixels for a 2D position histogram of the response: "<<numberOfPixel_tan<<" x "<< numberOfPixel_axial <<G4endl;
+				triggerEnergy = atof(value[1].c_str());
+				if(value[2] != "keV"){
+					G4cerr<<" The unit of Trigger energy threshold is not in keV, Make it in 'keV' "<<G4endl; 
+					exit(0);
+				}
+				triggerEnergy = triggerEnergy*keV;
+				G4cout<<"Trigger energy threshold: "<<triggerEnergy/keV<<" keV."<<G4endl;
+
+			}
+
+			//Option to apply AngerLogic
+			if( (std::string::size_type)inputLine.find("ApplyAngerLogic:")!=std::string::npos){
+				std::istringstream tmpStream(inputLine);
+				tmpStream >> value[0] >> value[1];
+				if(value[1]=="true"){
+					ApplyAngerLogic = true;
+					G4cout<<"Angler Logic calculation is applied. "<<G4endl;
+				}
+				else if(value[1]=="false") {
+					ApplyAngerLogic = false;
+					G4cout<<"Angler Logic calculation is NOT applied. "<<G4endl;
+				}
+				else {
+					ApplyAngerLogic = true;
+					G4cout<<"Angler Logic calculation is applied (by defualt). "<<G4endl;
+				}
+			}
+
+			//PMT position calculation blurring at FWHM
+			if( (std::string::size_type)inputLine.find("PMTblurring:")!=std::string::npos){
+				std::istringstream tmpStream(inputLine);
+				tmpStream >> value[0] >> value[1]>>value[2];
+				PMTblurring_tan = atof(value[1].c_str());
+				PMTblurring_axial = atof(value[2].c_str());
+				G4cout<<"PMTblurring position response blurring at FWHM (tan x axial) "<<PMTblurring_tan<<" x " <<PMTblurring_axial<<" mm2"<<G4endl;
 			}
 
 			//
@@ -657,8 +651,7 @@ void doiPETAnalysis::BlurringParameters(){
 	ifs.close();
 }
 
-//The following function reads the reflector pattern for each layer. 
-//Each layer has different patterns along the tangetial and axial positions.
+//The following function reads the reflector pattern for each layer. Each layer has different patterns along the tangetial and axial positions.
 //For defualt reflector pattern, see https://link.springer.com/article/10.1007/s12194-013-0231-4
 //The patter of the reflectors can be changed in the inputParameter.txt file
 //The pattern is given as 0 and 1. If there is reflector the value is 1 and if there is no reflector, the value is 0.
@@ -698,8 +691,9 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer1_Tangential.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
+			
 
 			//Reflector patter for Layer1 in the axial direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer1_Axial:")!=std::string::npos){
@@ -711,8 +705,8 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer1_Axial.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 
 			//Reflector patter for Layer2 in the tangential direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer2_Tangential:")!=std::string::npos){
@@ -724,8 +718,8 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer2_Tangential.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 
 			//Reflector patter for Layer2 in the axial direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer2_Axial:")!=std::string::npos){
@@ -737,8 +731,8 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer2_Axial.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 
 			//Reflector patter for Layer3 in the tangential direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer3_Tangential:")!=std::string::npos){
@@ -750,8 +744,8 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer3_Tangential.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 
 			//Reflector patter for Layer3 in the axial direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer3_Axial:")!=std::string::npos){
@@ -763,8 +757,8 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer3_Axial.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 
 			//Reflector patter for Layer4 in the tangential direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer4_Tangential:")!=std::string::npos){
@@ -776,8 +770,8 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer4_Tangential.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 
 			//Reflector patter for Layer4 in the axial direction
 			if( (std::string::size_type)inputLine.find("reflectorLayer4_Axial:")!=std::string::npos){
@@ -789,34 +783,108 @@ void doiPETAnalysis::ReadReflectorPattern(){
 						ireflectorLayer4_Axial.push_back(tmp_value);
 					}
 				}
+				stringReflectorValue.clear();
 			}
-			stringReflectorValue.clear();
 		}//#
 	}//while(eof)
 
-
-	//prepare Look up table for crystal identification. 
-	G4cout<<"DOI look-up table is being prepared. "<<G4endl;
-	std::string outputFileName = "check_2Dposition.txt";// excuted only once
-	std::ofstream outFile(outputFileName.c_str());
-
-	G4double crystalPositionY;
-	G4double crystalPositionZ;
+	//
+	//for debug
+	G4cout<<"\n========= Reflector Pattern ==========="<<G4endl;
+	G4cout<<"Layer 1"<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer1_Tangential.size();i++){
+		G4cout<<ireflectorLayer1_Tangential[i]<<" ";
+	}G4cout<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer1_Axial.size();i++){
+		G4cout<<ireflectorLayer1_Axial[i]<<" ";
+	}G4cout<<G4endl;
+	G4cout<<"Layer 2"<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer2_Tangential.size();i++){
+		G4cout<<ireflectorLayer2_Tangential[i]<<" ";
+	}G4cout<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer2_Axial.size();i++){
+		G4cout<<ireflectorLayer2_Axial[i]<<" ";
+	}G4cout<<G4endl;
+	G4cout<<"Layer 3"<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer3_Tangential.size();i++){
+		G4cout<<ireflectorLayer3_Tangential[i]<<" ";
+	}G4cout<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer3_Axial.size();i++){
+		G4cout<<ireflectorLayer3_Axial[i]<<" ";
+	}G4cout<<G4endl;
+	G4cout<<"Layer 4"<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer4_Tangential.size();i++){
+		G4cout<<ireflectorLayer4_Tangential[i]<<" ";
+	}G4cout<<G4endl;
+	for(unsigned int i = 0; i<ireflectorLayer4_Axial.size();i++){
+		G4cout<<ireflectorLayer4_Axial[i]<<" ";
+	}G4cout<<G4endl;
+	G4cout<<"========= Reflector Pattern Ended ===========\n"<<G4endl;
+	
+	//Read DOI look-up-table. This look-up-table is prepared based on the assumption that the interaction is occured at the center of the crystal.
+	G4int index_doi = 0, doiID;
 	doi_table.resize(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0);
 
-	for(G4int i_DOI = 0; i_DOI<numberOfCrystal_DOI; i_DOI++){
-		//crystalPositionX=(i_DOI-((float)numberOfCrystal_DOI)/2 + 0.5)*(sizeOfCrystal_DOI + crystalGap_DOI); //Becuase only lateral distances are used
-		for(G4int i_axial=0; i_axial< numberOfCrystal_axial;i_axial++){
-			crystalPositionZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial);
-			for(G4int i_tan=0; i_tan<numberOfCrystal_tangential;i_tan++){
-				crystalPositionY=(i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential);
-				AngerLogic(i_DOI, i_tan, i_axial, crystalPositionY, crystalPositionZ, 1, 0.5);
-				outFile<<i_DOI<<" "<<i_axial<<" "<<i_tan<<" "<<crystalID_in2D_posHist<<" "<<PositionAngerZ<<" "<<PositionAngerY<<G4endl;
-				doi_table[crystalID_in2D_posHist]=i_DOI;
+	std::string LUT_FileName = "look_up_table_DOI.txt";
+	std::ifstream ifs_doiLUT;
+	std::ofstream ofs_doiLUT;
+	ifs_doiLUT.open(LUT_FileName.c_str());
+	if(ifs_doiLUT.is_open()){
+		G4cout<<" DOI Look-up table found and used: File name: "<<LUT_FileName<<G4endl;
+		//Read from file
+		while(ifs_doiLUT>>index_doi>>doiID && index_doi < int(doi_table.size())){
+			doi_table[index_doi] = doiID;
+		}
+		if(index_doi==int(doi_table.size())){
+			G4cout<<"!!!Warning: The DOI table index is greater than the total number of crystals."<<G4endl;
+			PrepareDOILookUpTable(LUT_FileName);
+		}
+		isDOIlookUpTablePrepared = true; //
+	}
+	else
+	{
+		PrepareDOILookUpTable(LUT_FileName);
+	}
+	//Write into a file.
+	ofs_doiLUT.open(LUT_FileName.c_str());
+	if(!ofs_doiLUT.is_open()){
+		G4cerr<<"Unable to open file to write doi_LUT"<<G4endl;
+		exit(0);
+	}
+	for(G4int i=0;i<int(doi_table.size()); i++){
+		ofs_doiLUT<<i<<"\t"<<doi_table[i]<<G4endl;
+	}
+	ifs_doiLUT.close();
+
+
+	ifs.close();
+}
+void doiPETAnalysis::PrepareDOILookUpTable(G4String){
+	isDOIlookUpTablePrepared = false; 
+		G4cout<<"Preparing DOI look-up table... "<<G4endl;
+		std::string outputFileName = "_check_2Dposition.txt";// excuted only once
+		std::ofstream outFile(outputFileName.c_str());
+
+		G4double crystalPositionX;
+		G4double crystalPositionY;
+		G4double crystalPositionZ;
+		//doi_table.resize(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0);
+
+		for(G4int i_DOI = 0; i_DOI<numberOfCrystal_DOI; i_DOI++){
+			crystalPositionX=(i_DOI-((float)numberOfCrystal_DOI)/2 + 0.5)*(sizeOfCrystal_DOI + crystalGap_DOI); //Becuase only lateral distances are used
+			for(G4int i_axial=0; i_axial< numberOfCrystal_axial;i_axial++){
+				crystalPositionZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial);
+				for(G4int i_tan=0; i_tan<numberOfCrystal_tangential;i_tan++){
+					crystalPositionY=(i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential);
+					AngerLogic(crystalPositionX, crystalPositionY, crystalPositionZ, 1, 0.5, isDOIlookUpTablePrepared);
+					outFile<<PositionAngerZ<<" "<<PositionAngerY<<G4endl;
+					doi_table[crystalID_in2D_posHist]=i_DOI;
+
+				}
 			}
 		}
-	}
-	ifs.close();
+		G4cout<<"done."<<G4endl;
+		isDOIlookUpTablePrepared = true;
 }
 
 
@@ -824,8 +892,22 @@ void doiPETAnalysis::ReadReflectorPattern(){
 //The reflectors shifts the response by some distance so that the response can be projected into 2D position histogram. 
 //From this 2D position histogram, the new crystal ID (in 3D along the tangential (y), axial (z) and DOI (x)) (after Anger Logic method is applied) can be obtained.
 //If the crystal ID after Anger method apllied is out of the give number of crystals (in 3D), then an error message is displayed and the event will be rejected.
-void doiPETAnalysis::AngerLogic(G4int i_doi, G4int i_tan, G4int i_axial, G4double posCrystalY, G4double posCrystalZ, G4double Edep, G4double shiftDis)
-{
+void doiPETAnalysis::AngerLogic(G4double posCrystalX, G4double posCrystalY, G4double posCrystalZ, G4double Edep, G4double shiftDis, G4bool isDOI_LUT)
+{ 
+	G4double crystalPitch_DOI = sizeOfCrystal_DOI + crystalGap_DOI;
+	G4double crystalPitch_tan = sizeOfCrystal_tangential + crystalGap_tangential;
+	G4double crystalPitch_axial = sizeOfCrystal_axial + crystalGap_axial;
+
+	//The crystal ID are calculated based on the center of mass 
+	G4int i_doi = posCrystalX/crystalPitch_DOI + (float)numberOfCrystal_DOI*0.5;
+	G4int i_tan = posCrystalY/crystalPitch_tan + (float)numberOfCrystal_tangential*0.5;
+	G4int i_axial = posCrystalZ/crystalPitch_axial + (float)numberOfCrystal_axial*0.5;
+
+	//position of interaction is shifted  the centre of the crystal. This is to use DOI-look up tables as the real scanner
+	posCrystalX = (i_doi-((float)numberOfCrystal_DOI)/2 + 0.5)*crystalPitch_DOI;
+	posCrystalY = (i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*crystalPitch_tan;
+	posCrystalZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*crystalPitch_axial;
+
 	//1z and 2z are at the same z distance; 3z and 4z are at the same z distance
 	//The signal (the energy deposition) is devided into the four PMTs depending on their lateral distances (in the axial and tangential directions) from the interaction position
 
@@ -887,14 +969,6 @@ void doiPETAnalysis::AngerLogic(G4int i_doi, G4int i_tan, G4int i_axial, G4doubl
 	PositionAngerZ = (signalZplus - signalZminus)/(signalZplus + signalZminus)*distz; 
 	PositionAngerY = (signalYplus - signalYminus)/(signalYplus + signalYminus)*disty;
 
-	//Find local position in the interacting cystal to estimate the shift due to reflector
-	//double localPosX = posCrystalX - (i_doi-((G4double)numberOfCrystal_DOI)/2 + 0.5)*(sizeOfCrystal_DOI + crystalGap_DOI);
-	G4double localPosY = posCrystalY - (G4double)(i_tan-((G4double)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential);
-	G4double localPosZ = posCrystalZ - (G4double)(i_axial-((G4double)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial);
-
-
-	G4double crystalPitch_tan = sizeOfCrystal_tangential + crystalGap_tangential;
-	G4double crystalPitch_axial = sizeOfCrystal_axial + crystalGap_axial;
 
 	//For detectors with reflector insertion (light sharing), the response is shifted depending on the reflector patter.
 	//Here, it is assumed that the shift of the response is equal to half of the distance from the interaction position to the airgap in the lateral (transversal direction)
@@ -905,45 +979,50 @@ void doiPETAnalysis::AngerLogic(G4int i_doi, G4int i_tan, G4int i_axial, G4doubl
 	//Response shift for 1st Layer
 	if(i_doi == 0){
 		//If reflector is only in one (left) side of the crystal, then response shifts to the right side (away from the reflector)
-		if(ireflectorLayer1_Tangential[i_tan] == 1 && ireflectorLayer1_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2 - localPosY)*shiftDis;
+		if(ireflectorLayer1_Tangential[i_tan] == 1 && ireflectorLayer1_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
 
 		//If reflector is only in one (right) side of the crystal, then response shifts to the left side (away from the reflector)
-		if(ireflectorLayer1_Tangential[i_tan] == 0 && ireflectorLayer1_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2 - localPosY)*shiftDis;
+		if(ireflectorLayer1_Tangential[i_tan] == 0 && ireflectorLayer1_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
 
-		if(ireflectorLayer1_Axial[i_axial] == 1 && ireflectorLayer1_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2 - localPosZ)*shiftDis;
-		if(ireflectorLayer1_Axial[i_axial] == 0 && ireflectorLayer1_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2 - localPosZ)*shiftDis;
+		if(ireflectorLayer1_Axial[i_axial] == 1 && ireflectorLayer1_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
+		if(ireflectorLayer1_Axial[i_axial] == 0 && ireflectorLayer1_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
 	}
 	if(i_doi == 1){ //Response shift for 2nd Layer
-		if(ireflectorLayer2_Tangential[i_tan] == 1 && ireflectorLayer2_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2 - localPosY)*shiftDis;
-		if(ireflectorLayer2_Tangential[i_tan] == 0 && ireflectorLayer2_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2 - localPosY)*shiftDis;
+		if(ireflectorLayer2_Tangential[i_tan] == 1 && ireflectorLayer2_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
+		if(ireflectorLayer2_Tangential[i_tan] == 0 && ireflectorLayer2_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
 
-		if(ireflectorLayer2_Axial[i_axial] == 1 && ireflectorLayer2_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2 - localPosZ)*shiftDis;
-		if(ireflectorLayer2_Axial[i_axial] == 0 && ireflectorLayer2_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2 - localPosZ)*shiftDis;
+		if(ireflectorLayer2_Axial[i_axial] == 1 && ireflectorLayer2_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
+		if(ireflectorLayer2_Axial[i_axial] == 0 && ireflectorLayer2_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
 	}
 	if(i_doi == 2){ //Response shift for 3rd Layer
-		if(ireflectorLayer3_Tangential[i_tan] == 1 && ireflectorLayer3_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2 - localPosY)*shiftDis;
-		if(ireflectorLayer3_Tangential[i_tan] == 0 && ireflectorLayer3_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2 - localPosY)*shiftDis;
+		if(ireflectorLayer3_Tangential[i_tan] == 1 && ireflectorLayer3_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
+		if(ireflectorLayer3_Tangential[i_tan] == 0 && ireflectorLayer3_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
 
-		if(ireflectorLayer3_Axial[i_axial] == 1 && ireflectorLayer3_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2 - localPosZ)*shiftDis;
-		if(ireflectorLayer3_Axial[i_axial] == 0 && ireflectorLayer3_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2 - localPosZ)*shiftDis;
+		if(ireflectorLayer3_Axial[i_axial] == 1 && ireflectorLayer3_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
+		if(ireflectorLayer3_Axial[i_axial] == 0 && ireflectorLayer3_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
 	}
 	if(i_doi == 3){ //Response shift for 4th Layer
-		if(ireflectorLayer4_Tangential[i_tan] == 1 && ireflectorLayer4_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2 - localPosY)*shiftDis;
-		if(ireflectorLayer4_Tangential[i_tan] == 0 && ireflectorLayer4_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2 - localPosY)*shiftDis;
+		if(ireflectorLayer4_Tangential[i_tan] == 1 && ireflectorLayer4_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
+		if(ireflectorLayer4_Tangential[i_tan] == 0 && ireflectorLayer4_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
 
-		if(ireflectorLayer4_Axial[i_axial] == 1 && ireflectorLayer4_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2 - localPosZ)*shiftDis;
-		if(ireflectorLayer4_Axial[i_axial] == 0 && ireflectorLayer4_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2 - localPosZ)*shiftDis;
+		if(ireflectorLayer4_Axial[i_axial] == 1 && ireflectorLayer4_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
+		if(ireflectorLayer4_Axial[i_axial] == 0 && ireflectorLayer4_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
+	}
+   
+   //Blur the 2D position (obtained by ANger Logic method) to include uncertainity of the PMT position response.
+	if(isDOI_LUT){
+		PositionAngerZ = G4RandGauss::shoot(PositionAngerZ,PMTblurring_axial/2.35);
+		PositionAngerY = G4RandGauss::shoot(PositionAngerY,PMTblurring_tan/2.35);
 	}
-
 	//The main purpose of shifting the response is to be able to project the response of all the crytal elements into a 2D position histogram so that we can identify the DOI layer 
 	//by comparing with a look-up-table which is prepared based on the reflector insertion. 
 
 	//The crystal ID in 2D position histogram along the axial (z) direction. It can have values of: 0, 1, .. , 31, in 32x32 pixel position histogram
-	crystalID_in2D_posHist_axial = (G4int)((G4double)PositionAngerZ/(G4double)(crystalPitch_axial*0.5) + (G4double)(numberOfPixel_axial - 1.0)*0.5 + 0.5);//0.5 is added for round off
+	crystalID_in2D_posHist_axial = (G4int)(PositionAngerZ/(crystalPitch_axial*0.5) + (G4double)numberOfPixel_axial*0.5);//Note! crystalPitch_axial*0.5 is the pitch for the 32x32 2D pixel space, and 0.5 is added for round off
 
 	//The crystal ID in 2D position histogram along the tangential (y) direction. It can have values of: 0, 1, .. , 31, in 32x32 pixel position histogram
-	crystalID_in2D_posHist_tan =   (G4int)((G4double)PositionAngerY/(G4double)(crystalPitch_tan*0.5) + (G4double)(numberOfPixel_tan - 1.0)*0.5 + 0.5);//y_ID
-
+	crystalID_in2D_posHist_tan =   (G4int)(PositionAngerY/(crystalPitch_tan*0.5) + (G4double)numberOfPixel_tan * 0.5);
+	
 	//continuous crystal ID in the 2D position histogram. It will be from 0 to 1023 (in the case of 16x16x4 crystal array). 
 	crystalID_in2D_posHist = crystalID_in2D_posHist_axial + crystalID_in2D_posHist_tan * numberOfPixel_tan;//32;
 
@@ -969,8 +1048,146 @@ void doiPETAnalysis::AngerLogic(G4int i_doi, G4int i_tan, G4int i_axial, G4doubl
 	CrystalIDAfterAngerLogic(crystalIDNew_tan,crystalIDNew_axial,crystalIDNew_DOI);	
 }
 
+/////
 void doiPETAnalysis::CrystalIDAfterAngerLogic(G4int i_tan, G4int i_axial, G4int i_doi){
 	crystalID_tangential = i_tan;
 	crystalID_axial = i_axial;
 	DOI_ID = i_doi;
 }
+
+void doiPETAnalysis::book() 
+{ 
+  auto manager = G4AnalysisManager::Instance();
+  
+  //manager->SetVerboseLevel(2);
+ 
+  G4bool fileOpen = manager->OpenFile(rootFileName);
+  if (!fileOpen) {
+    G4cout << "\n---> HistoManager::book(): cannot open " 
+           << rootFileName << G4endl;
+    return;
+  }
+  // Create directories  
+  //manager->SetNtupleDirectoryName("ListModeData");
+
+  manager->SetFirstNtupleId(1);
+
+  if(getSinglesData){
+	manager -> CreateNtuple("Singles", "Singles");
+	fNtColId[0] = manager -> CreateNtupleIColumn("eventID");
+	fNtColId[1] = manager -> CreateNtupleIColumn("blockID");
+	//fNtColId[2] = manager -> CreateNtupleDColumn("crystalID_axial");
+	//fNtColId[3] = manager -> CreateNtupleDColumn("crystalID_tangential");
+	//fNtColId[4] = manager -> CreateNtupleDColumn("DOI_ID");
+	fNtColId[2] = manager -> CreateNtupleDColumn("timeStamp");
+	fNtColId[3] = manager -> CreateNtupleDColumn("totalEdep");
+
+	//Interaction position of the photon with the detector
+	fNtColId[4] = manager -> CreateNtupleDColumn("intPosX");
+	fNtColId[5] = manager -> CreateNtupleDColumn("intPosY");
+	fNtColId[6] = manager -> CreateNtupleDColumn("intPosZ");
+
+	////source position (annihilation position)
+	fNtColId[7] = manager -> CreateNtupleDColumn("spositionX");
+	fNtColId[8] = manager -> CreateNtupleDColumn("spositionY");
+	fNtColId[9] = manager -> CreateNtupleDColumn("spositionZ");
+
+	manager -> FinishNtuple();	  
+  }
+  if(getCoincidenceData){
+	  manager -> CreateNtuple("Coincidence", "Coincidence");
+	fNtColId[0] = manager -> CreateNtupleIColumn("eventID0");
+	fNtColId[1] = manager -> CreateNtupleIColumn("blockID0");
+	fNtColId[2] = manager -> CreateNtupleIColumn("crystalID_axial0");
+	fNtColId[3] = manager -> CreateNtupleIColumn("crystalID_tangential0");
+	fNtColId[4] = manager -> CreateNtupleIColumn("DOI_ID0");
+	fNtColId[5] = manager -> CreateNtupleDColumn("timeStamp0");
+	fNtColId[6] = manager -> CreateNtupleDColumn("totalEdep0");
+
+	fNtColId[7] = manager -> CreateNtupleIColumn("eventID1");
+	fNtColId[8] = manager -> CreateNtupleIColumn("blockID1");
+	fNtColId[9] = manager -> CreateNtupleIColumn("crystalID_axial1");
+	fNtColId[10] = manager -> CreateNtupleIColumn("crystalID_tangential1");
+	fNtColId[11] = manager -> CreateNtupleIColumn("DOI_ID1");
+	fNtColId[12] = manager -> CreateNtupleDColumn("timeStamp1");
+	fNtColId[13] = manager -> CreateNtupleDColumn("totalEdep1");
+
+	//source position
+	fNtColId[14] = manager -> CreateNtupleDColumn("spositionX");
+	fNtColId[15] = manager -> CreateNtupleDColumn("spositionY");
+	fNtColId[16] = manager -> CreateNtupleDColumn("spositionZ");
+
+	manager -> FinishNtuple();
+
+  }
+  
+  
+  factoryOn = true;    
+}
+void doiPETAnalysis::FillListModeEvent()
+{
+
+  auto manager = G4AnalysisManager::Instance();
+  if(getSinglesData){
+	 	manager -> FillNtupleIColumn(1, fNtColId[0], G4int(eventID));
+		manager -> FillNtupleIColumn(1, fNtColId[1], G4int(blockID));
+		//manager -> FillNtupleDColumn(1, fNtColId[2], crystalID_axial);
+		//manager -> FillNtupleDColumn(1, fNtColId[3], crystalID_tangential);
+		//manager -> FillNtupleDColumn(1, fNtColId[4], DOI_ID);
+		manager -> FillNtupleDColumn(1, fNtColId[2], timeStamp/s);// in second
+		manager -> FillNtupleDColumn(1, fNtColId[3], totalEdep/keV); //in keV
+		
+		//Interaction position of the photon in the detector
+		manager -> FillNtupleDColumn(1, fNtColId[4], intPosX); //mm
+		manager -> FillNtupleDColumn(1, fNtColId[5], intPosY); //mm
+		manager -> FillNtupleDColumn(1, fNtColId[6], intPosZ); //mm
+
+		//
+		//Add source position
+		manager -> FillNtupleDColumn(1, fNtColId[7], spositionX);
+		manager -> FillNtupleDColumn(1, fNtColId[8], spositionY);
+		manager -> FillNtupleDColumn(1, fNtColId[9], spositionZ);
+
+		manager -> AddNtupleRow(1);
+  }
+
+  if(getCoincidenceData){
+	  //First Single
+		manager -> FillNtupleIColumn(1, fNtColId[0], eventID0);
+		manager -> FillNtupleIColumn(1, fNtColId[1], blockID0);
+		manager -> FillNtupleIColumn(1, fNtColId[2], crystalID_axial0);
+		manager -> FillNtupleIColumn(1, fNtColId[3], crystalID_tangential0);
+		manager -> FillNtupleIColumn(1, fNtColId[4], DOI_ID0);
+		manager -> FillNtupleDColumn(1, fNtColId[5], timeStamp0/s);
+		manager -> FillNtupleDColumn(1, fNtColId[6], totalEdep0/keV);
+	
+	//Second Single
+		manager -> FillNtupleIColumn(1, fNtColId[7], eventID1);
+		manager -> FillNtupleIColumn(1, fNtColId[8], blockID1);
+		manager -> FillNtupleIColumn(1, fNtColId[9], crystalID_axial1);
+		manager -> FillNtupleIColumn(1, fNtColId[10], crystalID_tangential1);
+		manager -> FillNtupleIColumn(1, fNtColId[11], DOI_ID1);
+		manager -> FillNtupleDColumn(1, fNtColId[12], timeStamp1/s);
+		manager -> FillNtupleDColumn(1, fNtColId[13], totalEdep1/keV);
+	
+	//Add source position
+		manager -> FillNtupleDColumn(1, fNtColId[14], spositionX);
+		manager -> FillNtupleDColumn(1, fNtColId[15], spositionY);
+		manager -> FillNtupleDColumn(1, fNtColId[16], spositionZ);
+
+		manager -> AddNtupleRow(1);
+  }
+    
+}
+void doiPETAnalysis::finish() 
+{   
+ if (factoryOn) 
+   {
+    auto manager = G4AnalysisManager::Instance();    
+    manager -> Write();
+    manager -> CloseFile();  
+      
+   // delete G4AnalysisManager::Instance();
+    factoryOn = false;
+   }
+}

examples/advanced/doiPET/src/doiPETDetectorConstructionMessenger.cc

@@ -86,6 +86,12 @@ doiPETDetectorConstructionMessenger::doiPETDetectorConstructionMessenger(doiPETD
 	changePhantomLengthCmd->SetRange("length>0");
 	changePhantomLengthCmd->AvailableForStates(G4State_PreInit);
 
+	//Change the number of sleeves to be placed
+	changeNumOfSleevesCmd = new G4UIcmdWithAnInteger("/changePhantom/NumberOfSleeves",this);
+	changeNumOfSleevesCmd->SetParameterName("choice",false);
+	changeNumOfSleevesCmd->AvailableForStates(G4State_PreInit);
+
+
 }
 
 ////////// Destructor //////////////////////////////////////////////
@@ -97,6 +103,7 @@ doiPETDetectorConstructionMessenger::~doiPETDetectorConstructionMessenger()
 	delete changeThePhantomPositionCmd;
 	delete changePhantomRadiusCmd;
 	delete changePhantomLengthCmd;
+	delete changeNumOfSleevesCmd;
 
 
 }
@@ -118,6 +125,10 @@ void doiPETDetectorConstructionMessenger::SetNewValue(G4UIcommand* command,G4Str
 	{
 		Detector->SetPhantomLength(changePhantomLengthCmd->GetNewDoubleValue(newValue));
 	}
+	else if (command == changeNumOfSleevesCmd)
+	{
+		Detector->SetNumberOfSleeves(changeNumOfSleevesCmd->GetNewIntValue(newValue));
+	}
 }
 
 

examples/advanced/doiPET/src/doiPETRun.cc

@@ -26,6 +26,7 @@
 /// \file electromagnetic/TestEm11/src/doiPETRun.cc
 /// \brief Implementation of the doiPETRun class
 //
+// $Id: doiPETRun.cc 71376 2013-06-14 07:44:50Z maire $
 // 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -69,8 +70,6 @@ void doiPETRun::GetIntractionInfomation(InteractionInformation*step)
 
 void doiPETRun::FindInteractingCrystal()
 {
-
-	//G4cout<<"totalTime: "<<totalTime/s<<" "<<InitialActivity/becquerel<<G4endl;
 	std::multimap< G4int, InteractionInformation* >::iterator mitr;
 	std::set<G4int>::iterator sitr;
 	G4int cystalIDTemp;
@@ -83,21 +82,46 @@ void doiPETRun::FindInteractingCrystal()
 		for(G4int i=0; i<numberofInteractions; i++, mitr++){
 			edep = (*mitr).second->GetEdep();
 			cystalIDTemp = (*mitr).second->GetCrystalNo();
-			//edep_AfterCrystalBlurring = doiPETAnalysis::GetInstance()->QuantumEffifciency(edep);
 			edep_AfterCrystalBlurring = doiPETAnalysis::GetInstance()->QuantumEffifciency(edep, blockID, cystalIDTemp);
 			totalEdep += edep_AfterCrystalBlurring;
+
+			crystalID_vec.push_back(cystalIDTemp);
+			posInter_vec.push_back(((*mitr).second->GetInteractionPositionInCrystal()));
+			edepInCry_vec.push_back(edep_AfterCrystalBlurring);
+
 			//This is  to identify the crystal with the highest energy deposition. 
 			if(edepMax<edep_AfterCrystalBlurring){
 				edepMax=edep_AfterCrystalBlurring;
 				crystalID = cystalIDTemp;
 				interactionTime = (*mitr).second->GetGlobalTime();
-				interactionPos = (*mitr).second->GetInteractionPositionInCrystal();
+
+				//position of interaction based on highest edep. See below for energy weighted position of interaction
+				//interactionPos = (*mitr).second->GetInteractionPositionInCrystal();
 			}
 		}
+		if(totalEdep==0)continue;
+		
+		//The interaction position is based on energy weighted center of mass calculation
+		interactionPos = CenterOfMassInteractionPos(crystalID_vec, edepInCry_vec, totalEdep, posInter_vec);
 		doiPETAnalysis::GetInstance()->ReadOut(blockID,crystalID,interactionTime,timeOfAnnihil,interactionPos,totalEdep);
+
+		crystalID_vec.clear();
+		posInter_vec.clear();
+		edepInCry_vec.clear();
 	}
+	
 }
 
+//Apply center of mass calculation in determining the interaction position
+G4ThreeVector doiPETRun::CenterOfMassInteractionPos(const std::vector<G4int> &cryID, const std::vector<G4double> &energyDep, G4double edepTotal, const std::vector<G4ThreeVector> &pos){
+	posInterInCrystal = G4ThreeVector();
+	for(std::size_t i=0; i<cryID.size();i++){
+		posInterInCrystal += pos[i] * energyDep[i];
+	}
+	return (posInterInCrystal)/edepTotal;
+}
+
+
 void doiPETRun::SetAnnihilationTime(G4double t){
 	timeOfAnnihil = t;
 }

examples/advanced/doiPET/src/doiPETRunAction.cc

@@ -57,8 +57,15 @@
 #include <fstream>
 
 //
-doiPETRunAction::doiPETRunAction() : G4UserRunAction()
-{}
+#ifdef ANALYSIS_USE
+doiPETRunAction::doiPETRunAction(doiPETAnalysis* analysis) : G4UserRunAction()
+{
+	analysisMan = analysis;
+}
+#else
+doiPETRunAction::doiPETRunAction()
+{ }
+#endif
 //
 
 
@@ -76,10 +83,14 @@ G4Run* doiPETRunAction::GenerateRun()
 
 void doiPETRunAction::BeginOfRunAction(const G4Run* aRun)
 { 
+	//pass file name without fname extention
 	//inform the runManager to save random number seed
 	G4RunManager::GetRunManager()->SetRandomNumberStore(false);
 	G4cout << "### Begin of Run " << aRun->GetRunID() << " start." << G4endl;
-
+#ifdef ANALYSIS_USE
+  // Create ROOT file, histograms and ntuple
+  analysisMan -> book();
+#endif 
 	if (isMaster) {
 		G4Random::showEngineStatus();
 	}	
@@ -91,6 +102,10 @@ void doiPETRunAction::EndOfRunAction(const G4Run*)
 	G4int nofEvents = fRun->GetNumberOfEvent();
 	if (nofEvents == 0) return;
 	G4cout<<G4endl<<"### End of Run ("<<nofEvents<<" events)"<<G4endl;
+#ifdef ANALYSIS_USE
+// Close the output ROOT file with the results
+   analysisMan -> finish(); 
+#endif
 }
 
 

examples/advanced/doiPET/src/doiPETSteppingAction.cc

@@ -75,15 +75,13 @@ void doiPETSteppingAction::UserSteppingAction(const G4Step* aStep)
 	//G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID();	
 	G4ThreeVector pos = track->GetPosition();
 
-	//G4StepPoint* point1 = aStep->GetPreStepPoint();
-	//G4ThreeVector posIntCrystal = point1->GetPosition();
 
 	//G4StepPoint* p2 = aStep->GetPostStepPoint();
 	G4StepPoint* p1 = aStep->GetPreStepPoint();
 	G4ThreeVector coord1 = p1->GetPosition();
 
-	//The following is to get the local position of the with respect to the crystal volume
-	const G4AffineTransform transformation =  p1->GetTouchable()-> GetHistory()->GetTopTransform();
+	//The following is to get the local position of the interaction with respect to the crystal block volume
+	const G4AffineTransform transformation =  p1->GetTouchable()-> GetHistory()->GetTransform(2);
 	G4ThreeVector localPosition = transformation.TransformPoint(coord1);
 
 	G4int blockID;
@@ -107,9 +105,11 @@ void doiPETSteppingAction::UserSteppingAction(const G4Step* aStep)
 	//doiPETAnalysis::GetInstance()->SetEventID (eventID);
 
 	//Get scatter information in the phantom by the annihilation photon before detected by the detector. Note that the scatter index is initialized to 0. 
-	//If there is scatter, the index is 1, and if not it is 0.
-	if(edep>0 && (volumeName == "phantom_physicalV")) scatterIndex = 1;
-	doiPETAnalysis::GetInstance()->GetScatterIndexInPhantom(scatterIndex);
+	//If there energy deposition in the phantom by the photon, then there is scatter, the index is 1, and if not it is 0.
+	if(edep > 0 && particleName == "gamma"  && volumeName == "phantom_physicalV"){
+		scatterIndex = 1;
+		doiPETAnalysis::GetInstance()->SetScatterIndexInPhantom(scatterIndex);
+	}
 
 	///////////////////    Retrive (Extract) information in the crystal ///////////////////////////
 	if(edep>0. && volumeName=="Crystal_physicalV"){
@@ -145,7 +145,6 @@ void doiPETSteppingAction::UserSteppingAction(const G4Step* aStep)
 
 		//get the global time of the interaction with the crystal
 		ExtractIntInfo->SetGlobalTime( track->GetGlobalTime()); 
-		//G4cout<<"Step: "<<eventID<<" "<<blockID<<" "<<crystalID<<" "<<edep<<G4endl;
 
 		//pass all the obtained information to the doiPETAnalysis class
 		//doiPETAnalysis::GetInstance()->GetIntractionInfomation(ExtractIntInfo);