diff --git a/CMakeLists.txt b/CMakeLists.txt index 4b6b130..fcfd4a1 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,17 +6,21 @@ cmake_minimum_required(VERSION 2.6 FATAL_ERROR) #---------------------------------------------------------------------------- # Adding a variable to use root through out the G4 project. # -set(mytarget detectionSystems) # equivalent to "name := myexample" in G4 GNUmakefiles -set(useROOT true ) # use true or false (or comment to set to false) -set(myROOTclassDir dataRootClass ) # directory of the root classes -set(myROOTclass true ) # comment if none (please see http://root.cern.ch/phpBB3//viewtopic.php?t=6172) -set(myROOTS3class TS3Data.cpp ) # comment if none -set(myROOTFragmentclass TFragment.cpp ) # comment if none -set(myROOTTigFragmentclass TTigFragment.cpp) # comment if none -set(myROOTSpiceclass TSpiceData.cpp ) # comment if none -set(myROOTGriffclass TGriffinData.cpp) # comment if none -#set(myROOTNewclass TNewData.cpp) # comment if none +set(mytarget detectionSystems) # equivalent to "name := myexample" in G4 GNUmakefiles +set(useROOT true ) # use true or false (or comment to set to false) +set(myROOTclassDir dataRootClass ) # directory of the root classes +set(myROOTclass true ) # comment if none (please see http://root.cern.ch/phpBB3//viewtopic.php?t=6172) +set(myROOTS3class TS3Data.cpp ) # comment if none +set(myROOTSpiceclass TSpiceData.cpp ) # comment if none +set(myROOTPacesclass TPacesData.cpp ) # comment if none +set(myROOTNewclass TNewData.cpp) # comment if none +set(myROOTSceptarclass TSceptarData.cpp) # comment if none +set(myROOTGriffclass TGriffinData.cpp ) # comment if none +set(myROOTHistoryclass THistoryData.cpp ) # comment if none +set(myROOTFragmentclass TFragment.cpp ) # comment if none +set(myROOTTigFragmentclass TTigFragment.cpp) # comment if none +#set(myROOTNewclass TNewData.cpp) # comment if none # http://www.cmake.org/cmake/help/cmake_tutorial.html # http://www.cmake.org/cmake/help/cmake2.6docs.html @@ -57,7 +61,7 @@ include_directories(${PROJECT_SOURCE_DIR}/include # if(useROOT) if(myROOTclass) - + if(myROOTFragmentclass) message(" --- compiling Fragment --- ") EXECUTE_PROCESS(COMMAND rootcint -f dictFragment.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTFragmentclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/Fragment_linkdef.h ) @@ -67,7 +71,7 @@ if(useROOT) message(" --- compiling TigFragment --- ") EXECUTE_PROCESS(COMMAND rootcint -f dictTigFragment.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTTigFragmentclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/TigFragment_linkdef.h ) endif(myROOTTigFragmentclass) - + if(myROOTS3class) message(" --- compiling S3 --- ") EXECUTE_PROCESS(COMMAND rootcint -f dictS3.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTS3class} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/S3_linkdef.h ) @@ -78,11 +82,33 @@ if(useROOT) EXECUTE_PROCESS(COMMAND rootcint -f dictSpice.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTSpiceclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/Spice_linkdef.h ) endif(myROOTSpiceclass) + if(myROOTPacesclass) + message(" --- compiling Paces --- ") + EXECUTE_PROCESS(COMMAND rootcint -f dictPaces.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTPacesclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/Paces_linkdef.h ) + endif(myROOTPacesclass) + + if(myROOTNewclass) + message(" --- compiling New --- ") + EXECUTE_PROCESS(COMMAND rootcint -f dictNew.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTNewclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/New_linkdef.h) + endif(myROOTNewclass) + + if(myROOTSceptarclass) + message(" --- compiling Sceptar --- ") + EXECUTE_PROCESS(COMMAND rootcint -f dictSceptar.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTSceptarclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/Sceptar_linkdef.h) + endif(myROOTSceptarclass) + if(myROOTGriffclass) message(" --- compiling Griffin --- ") EXECUTE_PROCESS(COMMAND rootcint -f dictGriffin.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTGriffclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/Griffin_linkdef.h ) endif(myROOTGriffclass) + if(myROOTHistoryclass) + message(" --- compiling History --- ") + EXECUTE_PROCESS(COMMAND rootcint -f dictHistory.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTHistoryclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/History_linkdef.h ) + endif(myROOTHistoryclass) + + + #if(myROOTNewclass) #message(" --- compiling New --- ") #EXECUTE_PROCESS(COMMAND rootcint -f dictNew.cxx -c ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/${myROOTNewclass} ${PROJECT_SOURCE_DIR}/${myROOTclassDir}/New_linkdef.h) @@ -175,16 +201,36 @@ if(useROOT) set(myROOTDict ${myROOTDict} dictSpice.cxx) endif(myROOTSpiceclass) + if(myROOTPacesclass) + message(" --- Linking : dictPaces.cxx --- ") + set(myROOTDict ${myROOTDict} dictPaces.cxx) + endif(myROOTPacesclass) + if(myROOTGriffclass) message(" --- Linking : dictGriffin.cxx --- ") set(myROOTDict ${myROOTDict} dictGriffin.cxx) endif(myROOTGriffclass) - + + if(myROOTHistoryclass) + message(" --- Linking : dictHistory.cxx --- ") + set(myROOTDict ${myROOTDict} dictHistory.cxx) + endif(myROOTHistoryclass) + + if(myROOTNewclass) + message(" --- Linking : dictNew.cxx --- ") + set(myROOTDict ${myROOTDict} dictNew.cxx) + endif(myROOTNewclass) + + if(myROOTSceptarclass) + message(" --- Linking : dictSceptar.cxx --- ") + set(myROOTDict ${myROOTDict} dictSceptar.cxx) + endif(myROOTSceptarclass) + #if(myROOTNewclass) # message(" --- Linking : dictNew.cxx --- ") #set(myROOTDict ${myROOTDict} dictNew.cxx) #endif(myROOTNewclass) - + # This will replace the semicolons in a cmake list ";" to spaces " " foreach(arg ${myROOTDict}) set(myROOTDictSPACED "${myROOTDictSPACED} ${arg}") diff --git a/dataRootClass/History_linkdef.h b/dataRootClass/History_linkdef.h new file mode 100644 index 0000000..b17738f --- /dev/null +++ b/dataRootClass/History_linkdef.h @@ -0,0 +1,3 @@ +#ifdef __CINT__ +#pragma link C++ class THistoryData ; +#endif diff --git a/dataRootClass/Makefile b/dataRootClass/Makefile index e79381a..49497a4 100644 --- a/dataRootClass/Makefile +++ b/dataRootClass/Makefile @@ -31,10 +31,9 @@ GLIBS = $(ROOTGLIBS) $(SYSLIBS) # Specify Shared libraries #------------------------- -#SHARELIB = libDetector.so libNewData.so libNew2Data.so - -SHARELIB = libSpiceData.so libS3Data.so libGriffinData.so\ - libFragment.so libTigFragment.so +#SHARELIB = libDetector.so libNewData.so libNew2Data.so (example) +SHARELIB = libSpiceData.so libS3Data.so libPacesData.so libGriffinData.so libNewData.so\ + libHistoryData.so libFragment.so libTigFragment.so all: $(SHARELIB) @@ -62,6 +61,35 @@ libGriffinData.so: TGriffinData.o TGriffinDataDict.o TGriffinDataDict.cpp: TGriffinData.h rootcint -f $@ -c $^ +# P A C E S +libPacesData.so: TPacesData.o TPacesDataDict.o + $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ + +TPacesDataDict.cpp: TPacesData.h + rootcint -f $@ -c $^ + +# N E W +libNewData.so: TNewData.o TNewDataDict.o + $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ + +TNewDataDict.cpp: TNewData.h + rootcint -f $@ -c $^ + +# S C E P T A R +libSceptarData.so: TSceptarData.o TSceptarDataDict.o + $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ + +TSceptarDataDict.cpp: TSceptarData.h + rootcint -f $@ -c $^ + +# H I S T O R Y +libHistoryData.so: THistoryData.o THistoryDataDict.o + $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ + +THistoryDataDict.cpp: THistoryData.h + rootcint -f $@ -c $^ + + # F R A G M E N T (G R S I S P O O N) libFragment.so: TFragment.o TFragmentDict.o $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ @@ -76,7 +104,6 @@ libTigFragment.so: TTigFragment.o TTigFragmentDict.o TTigFragmentDict.cpp: TTigFragment.h rootcint -f $@ -c $^ - #libNewData.so: TNewData.o TNewDataDict.o # $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ # @@ -98,8 +125,13 @@ distclean: # Manage dependancies #------------------------- +TSpiceData.o: TSpiceData.cpp TSpiceData.h +TS3Data.o: TS3Data.cpp TS3Data.h +TPacesData.o: TPacesData.cpp TPacesData.h +TNewData.o: TNewData.cpp TNewData.h +TSceptarData.o: TSceptarData.cpp TSceptarData.h +TGriffinData.o: TGriffinData.cpp TGriffinData.h +THistoryData.o: THistoryData.cpp THistoryData.h TFragment.o: TFragment.cpp TFragment.h TTigFragment.o: TTigFragment.cpp TTigFragment.h -TSpiceData.o: TSpiceData.cpp TSpiceData.h -TS3Data.o: TS3Data.cpp TS3Data.h -TGriffinData.o: TGriffinData.cpp TGriffinData.h + diff --git a/dataRootClass/New_lindef.h b/dataRootClass/New_lindef.h deleted file mode 100644 index 5ffef14..0000000 --- a/dataRootClass/New_lindef.h +++ /dev/null @@ -1,3 +0,0 @@ -#ifdef __CINT__ -#pragma link C++ class TNewData ; -#endif diff --git a/dataRootClass/Paces_linkdef.h b/dataRootClass/Paces_linkdef.h new file mode 100644 index 0000000..c1de5b3 --- /dev/null +++ b/dataRootClass/Paces_linkdef.h @@ -0,0 +1,3 @@ +#ifdef __CINT__ +#pragma link C++ class TPacesData ; +#endif diff --git a/dataRootClass/README.md b/dataRootClass/README.md index 167907b..082b7e1 100644 --- a/dataRootClass/README.md +++ b/dataRootClass/README.md @@ -32,7 +32,7 @@ To add another "Setter" class for another detector : - Replace "New" by the name of your detector (the dictionary file must be the same) -- Finally, in ./src/RootManager.hh, add the header of the new class and use it +- Finally, in ./include/RootManager.hh, add the header of the new class and use it e.g. diff --git a/dataRootClass/Sceptar_linkdef.h b/dataRootClass/Sceptar_linkdef.h new file mode 100644 index 0000000..bdbb46b --- /dev/null +++ b/dataRootClass/Sceptar_linkdef.h @@ -0,0 +1,3 @@ +#ifdef __CINT__ +#pragma link C++ class TSceptarData ; +#endif diff --git a/dataRootClass/TFragment.cpp b/dataRootClass/TFragment.cpp index 897eea0..ba2c594 100644 --- a/dataRootClass/TFragment.cpp +++ b/dataRootClass/TFragment.cpp @@ -1,6 +1,6 @@ #include"TFragment.h" -ClassImp(TFragment); +ClassImp(TFragment) TFragment::TFragment() { diff --git a/dataRootClass/TGriffinData.cpp b/dataRootClass/TGriffinData.cpp index d442b09..ff09ff0 100644 --- a/dataRootClass/TGriffinData.cpp +++ b/dataRootClass/TGriffinData.cpp @@ -24,13 +24,13 @@ ClassImp(TGriffinData) TGriffinData::TGriffinData() { // Default constructor - Clear(); + ClearVariables(); } TGriffinData::~TGriffinData() {} -void TGriffinData::Clear() +void TGriffinData::ClearVariables() { fDetNbr.clear(); fPrimaryEnergy.clear(); diff --git a/dataRootClass/TGriffinData.h b/dataRootClass/TGriffinData.h index 2066d8b..fd4ddfb 100644 --- a/dataRootClass/TGriffinData.h +++ b/dataRootClass/TGriffinData.h @@ -40,14 +40,14 @@ class TGriffinData : public TObject { TGriffinData(); virtual ~TGriffinData(); - void Clear(); + void ClearVariables(); void Dump() const; ///////////////////// GETTERS //////////////////////// Int_t GetGriffinEDetectorNbr(Int_t i) {return fDetNbr.at(i);} - Int_t GetGriffinEMult(Int_t i) {return fDetNbr.size();} + Int_t GetGriffinEMult(Int_t) {return fDetNbr.size();} TVector3 GetPositionFirstHit(Int_t i) {return fPositionFirstHit.at(i);} Int_t GetEventNumber(void) {return fEventNumber;} diff --git a/dataRootClass/THistoryData.cpp b/dataRootClass/THistoryData.cpp new file mode 100644 index 0000000..c0c925f --- /dev/null +++ b/dataRootClass/THistoryData.cpp @@ -0,0 +1,134 @@ +/***************************************************************************** + * Original Author: Mhd Moukaddam contact address: moukaddam@triumf.ca * + *---------------------------------------------------------------------------* + * Decription: This class stores selected information of the G4 simulation * + * to help in analysing the final simulated physics spectra. * + * Information such as, the number of generation in a cascade, * + * This class derives from TObject (ROOT) and its aim is to be * + * stored in the output TTree of the G4 simulation * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include +using namespace std; + +#include "THistoryData.h" + + +ClassImp(THistoryData) + +THistoryData::THistoryData() { + ClearVariables(); +} + +THistoryData::~THistoryData() {} + + +void THistoryData::ClearVariables(){ + + fHistoryPrimaryID.clear(); + fHistoryPrimaryPdg.clear(); + fHistoryPrimaryEnergy.clear(); + + fHistoryPrimaryPositionVertex.clear(); + fHistoryPrimaryMomentumVertex.clear(); + fHistoryPrimaryPosition1stImpact.clear(); + fHistoryPrimaryMomentum1stImpact.clear(); + fHistoryPrimaryVolume1stImpact.clear(); + + fHistoryParentID.clear(); + + fHistoryCurrentCreatorProcess.clear(); + fHistoryCurrentID.clear(); + fHistoryCurrentPdg.clear(); + fHistoryCurrentEnergy.clear(); + fHistoryCurrentTime.clear(); + fHistoryCurrentVolume1stImpact.clear(); + + fHistoryCurrentPositionVertex.clear(); + fHistoryCurrentPosition1stImpact.clear(); + fHistoryCurrentPositionDetector.clear(); + fHistoryCurrentPositionDeath.clear(); + + fHistoryCurrentMomentumVertex.clear(); + fHistoryCurrentMomentum1stImpact.clear(); + fHistoryCurrentMomentumDetector.clear(); + fHistoryCurrentMomentumDeath.clear(); + +} + + +void THistoryData::Dump() const +{ + cout << "XXXXXXXXXXXXXXXXXXXXXXXX HISTORY XXXXXXXXXXXXXXXXXXXXXXXX" << endl; + + cout << "History_Mult = " < + +//ROOT +#include "TObject.h" +#include "TVector3.h" + +class THistoryData : public TObject { + + private: +vector fHistoryPrimaryID; +vector fHistoryPrimaryPdg; // the type of the particle generated in the source e.g. pdg = 22 (gamma) in 0 --> gamma (primary) --> e- --> gamma --> e- +vector fHistoryPrimaryEnergy; + +vector fHistoryPrimaryPositionVertex; +vector fHistoryPrimaryMomentumVertex; +vector fHistoryPrimaryPosition1stImpact; +vector fHistoryPrimaryMomentum1stImpact; +vector fHistoryPrimaryVolume1stImpact; + +vector fHistoryParentID; + +vector fHistoryCurrentID; +vector fHistoryCurrentPdg; //the type of the last particle produced e.g. pdg = 11 (e-) in 0 --> gamma (primary) --> e- --> gamma --> e- +vector fHistoryCurrentEnergy; +vector fHistoryCurrentTime; // time when the detector is hit +vector fHistoryCurrentCreatorProcess; +vector fHistoryCurrentVolume1stImpact; + +vector fHistoryCurrentPositionVertex; +vector fHistoryCurrentPosition1stImpact; +vector fHistoryCurrentPositionDetector; +vector fHistoryCurrentPositionDeath; + +vector fHistoryCurrentMomentumVertex; +vector fHistoryCurrentMomentum1stImpact; +vector fHistoryCurrentMomentumDetector; +vector fHistoryCurrentMomentumDeath; + + + public: + THistoryData(); + virtual ~THistoryData(); + + void ClearVariables(); + void Dump() const; + + + + ///////////////////// GETTERS //////////////////////// +inline Int_t GetHistoryMult(void) { return fHistoryPrimaryEnergy.size(); } +inline Double_t GetHistoryPrimaryID(int i) { return fHistoryPrimaryID.at(i); } +inline Int_t GetHistoryPrimaryPdg(int i) { return fHistoryPrimaryPdg.at(i); } +inline Double_t GetHistoryPrimaryEnergy(int i) { return fHistoryPrimaryEnergy.at(i); } + +inline TVector3 GetHistoryPrimaryPositionVertex(int i) { return fHistoryPrimaryPositionVertex.at(i); } +inline TVector3 GetHistoryPrimaryMomentumVertex(int i) { return fHistoryPrimaryMomentumVertex.at(i); } +inline TVector3 GetHistoryPrimaryMomentum1stImpact(int i) { return fHistoryPrimaryMomentum1stImpact.at(i); } +inline TVector3 GetHistoryPrimaryPosition1stImpact(int i) { return fHistoryPrimaryPosition1stImpact.at(i); } +inline TString GetHistoryPrimaryVolume1stImpact(int i) { return fHistoryPrimaryVolume1stImpact.at(i); } + +inline Int_t GetHistoryParentID(int i) { return fHistoryParentID.at(i); } + +inline Double_t GetHistoryCurrentID(int i) { return fHistoryCurrentID.at(i); } +inline Int_t GetHistoryCurrentPdg(int i) { return fHistoryCurrentPdg.at(i); } +inline Double_t GetHistoryCurrentEnergy(int i) { return fHistoryCurrentEnergy.at(i); } +inline Double_t GetHistoryCurrentTime(int i) { return fHistoryCurrentTime.at(i); } +inline TString GetHistoryCurrentCreatorProcess(int i) { return fHistoryCurrentCreatorProcess.at(i); } + +inline TVector3 GetHistoryCurrentPositionVertex(int i) { return fHistoryCurrentPositionVertex.at(i); } +inline TVector3 GetHistoryCurrentPosition1stImpact(int i) { return fHistoryCurrentPosition1stImpact.at(i); } +inline TVector3 GetHistoryCurrentPositionDetector(int i) { return fHistoryCurrentPositionDetector.at(i); } +inline TVector3 GetHistoryCurrentPositionDeath(int i) { return fHistoryCurrentPositionDeath.at(i); } + +inline TVector3 GetHistoryCurrentMomentumVertex(int i) { return fHistoryCurrentMomentumVertex.at(i); } +inline TVector3 GetHistoryCurrentMomentum1stImpact(int i) { return fHistoryCurrentMomentum1stImpact.at(i); } +inline TVector3 GetHistoryCurrentMomentumDetector(int i) { return fHistoryCurrentMomentumDetector.at(i); } +inline TVector3 GetHistoryCurrentMomentumDeath(int i) { return fHistoryCurrentMomentumDeath.at(i); } + + + ///////////////////// SETTERS //////////////////////// +inline void SetHistoryPrimaryID(Double_t i) { fHistoryPrimaryID.push_back(i); } +inline void SetHistoryPrimaryPdg(Int_t i) { fHistoryPrimaryPdg.push_back(i); } +inline void SetHistoryPrimaryEnergy(Double_t i) { fHistoryPrimaryEnergy.push_back(i); } + +inline void SetHistoryPrimaryPositionVertex(Double_t x, Double_t y, Double_t z) { fHistoryPrimaryPositionVertex.push_back(TVector3(x,y,z)); } +inline void SetHistoryPrimaryMomentumVertex(Double_t x, Double_t y, Double_t z) { fHistoryPrimaryMomentumVertex.push_back(TVector3(x,y,z)); } +inline void SetHistoryPrimaryMomentum1stImpact(Double_t x, Double_t y, Double_t z) { fHistoryPrimaryMomentum1stImpact.push_back(TVector3(x,y,z)); } +inline void SetHistoryPrimaryPosition1stImpact(Double_t x, Double_t y, Double_t z) { fHistoryPrimaryPosition1stImpact.push_back(TVector3(x,y,z)); } +inline void SetHistoryPrimaryVolume1stImpact(TString i) { fHistoryPrimaryVolume1stImpact.push_back(i); } + +inline void SetHistoryParentID(Int_t i) { fHistoryParentID.push_back(i); } + +inline void SetHistoryCurrentCreatorProcess(TString i) { fHistoryCurrentCreatorProcess.push_back(i); } +inline void SetHistoryCurrentID(Double_t i) { fHistoryCurrentID.push_back(i); } +inline void SetHistoryCurrentPdg(Int_t i) { fHistoryCurrentPdg.push_back(i); } +inline void SetHistoryCurrentEnergy(Double_t i) { fHistoryCurrentEnergy.push_back(i); } +inline void SetHistoryCurrentTime(Double_t i) { fHistoryCurrentTime.push_back(i); } // Time whe the detector is hit +inline void SetHistoryCurrentVolume1stImpact(TString i) { fHistoryCurrentVolume1stImpact.push_back(i); } + + +inline void SetHistoryCurrentPositionVertex(Double_t x, Double_t y, Double_t z) { fHistoryCurrentPositionVertex.push_back(TVector3(x,y,z)); } +inline void SetHistoryCurrentPosition1stImpact(Double_t x, Double_t y, Double_t z) { fHistoryCurrentPosition1stImpact.push_back(TVector3(x,y,z)); } +inline void SetHistoryCurrentPositionDetector(Double_t x, Double_t y, Double_t z) { fHistoryCurrentPositionDetector.push_back(TVector3(x,y,z)); } +inline void SetHistoryCurrentPositionDeath(Double_t x, Double_t y, Double_t z) { fHistoryCurrentPositionDeath.push_back(TVector3(x,y,z)); } + +inline void SetHistoryCurrentMomentumVertex(Double_t x, Double_t y, Double_t z) { fHistoryCurrentMomentumVertex.push_back(TVector3(x,y,z)); } +inline void SetHistoryCurrentMomentum1stImpact(Double_t x, Double_t y, Double_t z) { fHistoryCurrentMomentum1stImpact.push_back(TVector3(x,y,z)); } +inline void SetHistoryCurrentMomentumDetector(Double_t x, Double_t y, Double_t z) { fHistoryCurrentMomentumDetector.push_back(TVector3(x,y,z)); } +inline void SetHistoryCurrentMomentumDeath(Double_t x, Double_t y, Double_t z) { fHistoryCurrentMomentumDeath.push_back(TVector3(x,y,z)); } + + ClassDef(THistoryData,1) // HistoryData structure +}; + +#endif diff --git a/dataRootClass/TNewData.cpp b/dataRootClass/TNewData.cpp new file mode 100644 index 0000000..1633353 --- /dev/null +++ b/dataRootClass/TNewData.cpp @@ -0,0 +1,69 @@ + + +#include +using namespace std; + +#include "TNewData.h" + + +ClassImp(TNewData) + +TNewData::TNewData() +{ + ClearVariables(); +} + +TNewData::~TNewData() +{ +} + +void TNewData::ClearVariables() +{ + + fNew_DetNbr.clear(); + fNew_Energy.clear(); + + fPrimaryTheta.clear(); + fPrimaryPhi.clear(); + fPrimaryEnergy.clear(); + fPrimaryPdg.clear(); + fPdg.clear(); + fPositionFirstHit.clear(); + fEventNumber = -1 ; + +} + +void TNewData::Dump() const +{ + + cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event : " << fEventNumber << "XXXXXXXXXXXXXXXXXXXXXXXX" << endl; + + // DSSD + // (Th,E) + cout << "New_MultE = " << fNew_DetNbr.size() << endl; + for (UShort_t i = 0; i < fNew_DetNbr.size(); i++) + cout << "DetNbr: " << fNew_DetNbr[i] << "Energy: " << fNew_Energy[i] << endl; + + // First hit + cout << "Position of First Hit Mult = " << fPositionFirstHit.size() << endl; + for (UShort_t i = 0; i < fPositionFirstHit.size(); i++) + cout << " X - " << fPositionFirstHit.at(i).X() << " Y - " < + +//ROOT +#include "TObject.h" +#include "TVector3.h" + +class TNewData : public TObject +{ + +private: + + vector fNew_DetNbr; + vector fNew_Energy; + + vector fPositionFirstHit; + + vector fPrimaryTheta; + vector fPrimaryPhi; + vector fPrimaryEnergy; + vector fPrimaryPdg; + vector fPdg; + + Int_t fEventNumber; + +public: + + TNewData(); + virtual ~TNewData(); + + void ClearVariables(); + void Dump() const; + + ///////////////////// GETTERS //////////////////////// + Int_t GetNewMult() {return fNew_DetNbr.size();} + Int_t GetNewDetectorNbr(Int_t i) {return fNew_DetNbr.at(i);} + Double_t GetNewEnergy(Int_t i) {return fNew_Energy.at(i);} + + + TVector3 GetPositionFirstHit(Int_t i) {return fPositionFirstHit.at(i);} + Int_t GetEventNumber(void) {return fEventNumber;} + + Double_t GetPrimaryTheta(Int_t i) {return fPrimaryTheta.at(i);} + Double_t GetPrimaryPhi(Int_t i) {return fPrimaryPhi.at(i);} + Double_t GetPrimaryEnergy(Int_t i) {return fPrimaryEnergy.at(i);} + Int_t GetPrimaryPdg(Int_t i) {return fPrimaryPdg.at(i);} + Int_t GetPdg(Int_t i) {return fPdg.at(i);} + + ///////////////////// SETTERS //////////////////////// + void SetNewDetectorNbr(Int_t det) {fNew_DetNbr.push_back(det);} + void SetNewEnergy(Double_t E) {fNew_Energy.push_back(E);} + + void SetPrimaryTheta(Double_t theta) {fPrimaryTheta.push_back(theta);} + void SetPrimaryPhi(Double_t phi) {fPrimaryPhi.push_back(phi);} + void SetPrimaryEnergy(Double_t energy) {fPrimaryEnergy.push_back(energy);} + void SetPrimaryPdg(Int_t pdg) {fPrimaryPdg.push_back(pdg);} + void SetPdg(Int_t pdg) {fPdg.push_back(pdg);} + + void SetPositionFirstHit(TVector3 position) {fPositionFirstHit.push_back(position);} + void SetEventNumber(Int_t i) {fEventNumber = i;} + + ClassDef(TNewData,1) + +}; + +#endif + + + + + + diff --git a/dataRootClass/TPacesData.cpp b/dataRootClass/TPacesData.cpp new file mode 100644 index 0000000..6b8419a --- /dev/null +++ b/dataRootClass/TPacesData.cpp @@ -0,0 +1,79 @@ +/***************************************************************************** + * Original Author : Mhd Moukaddam contact address: moukaddam@triumf.ca * + *---------------------------------------------------------------------------* + * Decription: This class stores the results of the G4 simulation for the * + * PACESdetector. And was adapted from S1 detector Class. * + * This class derives from TObject (ROOT) and its aim is to be * + * stored in the output TTree of the G4 simulation * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include +using namespace std; + +#include "TPacesData.h" + + +ClassImp(TPacesData) + +TPacesData::TPacesData() { + ClearVariables(); +} + +TPacesData::~TPacesData() { +} + + +void TPacesData::ClearVariables() { + + fPaces_DetNbr.clear(); + fPaces_Energy.clear(); + + fPrimaryTheta.clear(); + fPrimaryPhi.clear(); + fPrimaryEnergy.clear(); + fPrimaryPdg.clear(); + fPdg.clear(); + fPositionFirstHit.clear(); + fEventNumber = -1 ; +} + + +void TPacesData::Dump() const +{ + cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event : " << fEventNumber << " XXXXXXXXXXXXXXXXXXXXXXXX" << endl; + + // DSSD + // (Th,E) + cout << "Paces_MultE = " << fPaces_DetNbr.size() << endl; + for (UShort_t i = 0; i < fPaces_DetNbr.size(); i++) + cout << "DetNbr: " << fPaces_DetNbr[i] << endl; + for (UShort_t i = 0; i < fPaces_Energy.size(); i++) + cout << "Energy: " << fPaces_Energy[i] << endl; + + // First hit + cout << "Position of First Hit Mult = " << fPositionFirstHit.size() << endl; + for (UShort_t i = 0; i < fPositionFirstHit.size(); i++) + cout << " X - " << fPositionFirstHit.at(i).X() << " Y - " < + +//ROOT +#include "TObject.h" +#include "TVector3.h" + +class TPacesData : public TObject { + private: + // SiLi + vector fPaces_DetNbr; + vector fPaces_Energy; + + vector fPositionFirstHit; + + vector fPrimaryTheta; + vector fPrimaryPhi; + vector fPrimaryEnergy; + vector fPrimaryPdg; + vector fPdg; + + Int_t fEventNumber; + + public: + TPacesData(); + virtual ~TPacesData(); + + void ClearVariables(); + void Dump() const; + + ///////////////////// GETTERS //////////////////////// + Int_t GetPacesMult() {return fPaces_DetNbr.size();} + Int_t GetPacesDetectorNbr(Int_t i) {return fPaces_DetNbr.at(i);} + Double_t GetPacesEnergy(Int_t i) {return fPaces_Energy.at(i);} + + + TVector3 GetPositionFirstHit(Int_t i) {return fPositionFirstHit.at(i);} + Int_t GetEventNumber(void) {return fEventNumber;} + + Double_t GetPrimaryTheta(Int_t i) {return fPrimaryTheta.at(i);} + Double_t GetPrimaryPhi(Int_t i) {return fPrimaryPhi.at(i);} + Double_t GetPrimaryEnergy(Int_t i) {return fPrimaryEnergy.at(i);} + Int_t GetPrimaryPdg(Int_t i) {return fPrimaryPdg.at(i);} + Int_t GetPdg(Int_t i) {return fPdg.at(i);} + + ///////////////////// SETTERS //////////////////////// + void SetPacesDetectorNbr(Int_t det) {fPaces_DetNbr.push_back(det);} + void SetPacesEnergy(Double_t E) {fPaces_Energy.push_back(E);} + + void SetPositionFirstHit(TVector3 position) {fPositionFirstHit.push_back(position);} + void SetEventNumber(Int_t i) {fEventNumber = i;} + + void SetPrimaryTheta(Double_t theta) {fPrimaryTheta.push_back(theta);} + void SetPrimaryPhi(Double_t phi) {fPrimaryPhi.push_back(phi);} + void SetPrimaryEnergy(Double_t energy) {fPrimaryEnergy.push_back(energy);} + void SetPrimaryPdg(Int_t pdg) {fPrimaryPdg.push_back(pdg);} + void SetPdg(Int_t pdg) {fPdg.push_back(pdg);} + + ClassDef(TPacesData,1) // PacesData structure +}; + +#endif diff --git a/dataRootClass/TS3Data.cpp b/dataRootClass/TS3Data.cpp index a776d22..7262b17 100644 --- a/dataRootClass/TS3Data.cpp +++ b/dataRootClass/TS3Data.cpp @@ -22,14 +22,14 @@ using namespace std; ClassImp(TS3Data) TS3Data::TS3Data() { - Clear(); + ClearVariables(); } TS3Data::~TS3Data() { } -void TS3Data::Clear() { +void TS3Data::ClearVariables() { // DSSD // (Th,E) fS3_Theta_DetNbr.clear(); diff --git a/dataRootClass/TS3Data.h b/dataRootClass/TS3Data.h index ddf5934..c832b7a 100644 --- a/dataRootClass/TS3Data.h +++ b/dataRootClass/TS3Data.h @@ -51,7 +51,7 @@ class TS3Data : public TObject { TS3Data(); virtual ~TS3Data(); - void Clear(); + void ClearVariables(); void Dump() const; diff --git a/dataRootClass/TSceptarData.cpp b/dataRootClass/TSceptarData.cpp new file mode 100644 index 0000000..670d914 --- /dev/null +++ b/dataRootClass/TSceptarData.cpp @@ -0,0 +1,67 @@ + + +#include +using namespace std; + +#include "TSceptarData.h" + + +ClassImp(TSceptarData) + +TSceptarData::TSceptarData() +{ + ClearVariables(); +} + +TSceptarData::~TSceptarData() +{ +} + +void TSceptarData::ClearVariables() +{ + fSceptar_DetNbr.clear(); + fSceptar_Energy.clear(); + + fPrimaryTheta.clear(); + fPrimaryPhi.clear(); + fPrimaryEnergy.clear(); + fPrimaryPdg.clear(); + fPdg.clear(); + fPositionFirstHit.clear(); + fEventNumber = -1 ; +} + +void TSceptarData::Dump() const +{ + + cout << "XXXXXXXXXXXXXXXXXXXXXXXX Sceptar Event : " << fEventNumber << "XXXXXXXXXXXXXXXXXXXXXXXX" << endl; + + // DSSD + // (Th,E) + cout << "Sceptar_MultE = " << fSceptar_DetNbr.size() << endl; + for (UShort_t i = 0; i < fSceptar_DetNbr.size(); i++) + cout << "DetNbr: " << fSceptar_DetNbr[i] << "Energy: " << fSceptar_Energy[i] << endl; + + // First hit + cout << "Position of First Hit Mult = " << fPositionFirstHit.size() << endl; + for (UShort_t i = 0; i < fPositionFirstHit.size(); i++) + cout << " X - " << fPositionFirstHit.at(i).X() << " Y - " < + +//ROOT +#include "TObject.h" +#include "TVector3.h" + +class TSceptarData : public TObject +{ + +private: + + vector fSceptar_DetNbr; + vector fSceptar_Energy; + + vector fPositionFirstHit; + + vector fPrimaryTheta; + vector fPrimaryPhi; + vector fPrimaryEnergy; + vector fPrimaryPdg; + vector fPdg; + + Int_t fEventNumber; + +public: + + TSceptarData(); + virtual ~TSceptarData(); + + void ClearVariables(); + void Dump() const; + + ///////////////////// GETTERS //////////////////////// + Int_t GetSceptarMult() {return fSceptar_DetNbr.size();} + Int_t GetSceptarDetectorNbr(Int_t i) {return fSceptar_DetNbr.at(i);} + Double_t GetSceptarEnergy(Int_t i) {return fSceptar_Energy.at(i);} + + + TVector3 GetPositionFirstHit(Int_t i) {return fPositionFirstHit.at(i);} + Int_t GetEventNumber(void) {return fEventNumber;} + + Double_t GetPrimaryTheta(Int_t i) {return fPrimaryTheta.at(i);} + Double_t GetPrimaryPhi(Int_t i) {return fPrimaryPhi.at(i);} + Double_t GetPrimaryEnergy(Int_t i) {return fPrimaryEnergy.at(i);} + Int_t GetPrimaryPdg(Int_t i) {return fPrimaryPdg.at(i);} + Int_t GetPdg(Int_t i) {return fPdg.at(i);} + + ///////////////////// SETTERS //////////////////////// + void SetSceptarDetectorNbr(Int_t det) {fSceptar_DetNbr.push_back(det);} + void SetSceptarEnergy(Double_t E) {fSceptar_Energy.push_back(E);} + + void SetPrimaryTheta(Double_t theta) {fPrimaryTheta.push_back(theta);} + void SetPrimaryPhi(Double_t phi) {fPrimaryPhi.push_back(phi);} + void SetPrimaryEnergy(Double_t energy) {fPrimaryEnergy.push_back(energy);} + void SetPrimaryPdg(Int_t pdg) {fPrimaryPdg.push_back(pdg);} + void SetPdg(Int_t pdg) {fPdg.push_back(pdg);} + + void SetPositionFirstHit(TVector3 position) {fPositionFirstHit.push_back(position);} + void SetEventNumber(Int_t i) {fEventNumber = i;} + + ClassDef(TSceptarData,1) + +}; + +#endif + + + + + + diff --git a/dataRootClass/TSpiceData.cpp b/dataRootClass/TSpiceData.cpp index 7516b2d..f41107c 100644 --- a/dataRootClass/TSpiceData.cpp +++ b/dataRootClass/TSpiceData.cpp @@ -24,13 +24,13 @@ ClassImp(TSpiceData) TSpiceData::TSpiceData() { // Default constructor - Clear(); + ClearVariables(); } TSpiceData::~TSpiceData() {} -void TSpiceData::Clear() +void TSpiceData::ClearVariables() { // DSSD // (Th,E) diff --git a/dataRootClass/TSpiceData.h b/dataRootClass/TSpiceData.h index c0d96aa..cc0509c 100644 --- a/dataRootClass/TSpiceData.h +++ b/dataRootClass/TSpiceData.h @@ -53,7 +53,7 @@ class TSpiceData : public TObject { TSpiceData(); virtual ~TSpiceData(); - void Clear(); + void ClearVariables(); void Dump() const; diff --git a/dataRootClass/TTigFragment.cpp b/dataRootClass/TTigFragment.cpp index 985c702..cea6dab 100644 --- a/dataRootClass/TTigFragment.cpp +++ b/dataRootClass/TTigFragment.cpp @@ -1,6 +1,6 @@ #include "TTigFragment.h" -ClassImp(TTigFragment); +ClassImp(TTigFragment) TTigFragment::TTigFragment() { diff --git a/include/ApparatusSpiceTargetChamber.hh b/include/ApparatusSpiceTargetChamber.hh index e078e73..64fee1f 100755 --- a/include/ApparatusSpiceTargetChamber.hh +++ b/include/ApparatusSpiceTargetChamber.hh @@ -96,11 +96,9 @@ private: G4LogicalVolume* mclamp_chamber_log; G4LogicalVolume* mclamp_shield_log; G4LogicalVolume* electro_box_log; - G4LogicalVolume* s3_mount_log; G4LogicalVolume* shield_cover_log; - G4LogicalVolume* detector_mount_log; - G4LogicalVolume* annular_clamp_log; - + G4LogicalVolume* cold_finger_log; + private: //////////////////////////////////////////// // Physical Volumes used in ApparatusSpiceTargetChamber @@ -130,11 +128,9 @@ private: G4VPhysicalVolume* mclamp_chamber_phys; G4VPhysicalVolume* mclamp_shield_phys; G4VPhysicalVolume* electro_box_phys; - G4VPhysicalVolume* s3_mount_phys; G4VPhysicalVolume* shield_cover_phys; - G4VPhysicalVolume* detector_mount_phys; - G4VPhysicalVolume* annular_clamp_phys; - + G4VPhysicalVolume* cold_finger_phys; + private: //////////////////////////////////////////// // Properties used in ApparatusSpiceTargetChamber @@ -166,14 +162,11 @@ private: G4String gear_plate_material; G4String target_mount_plate_material; G4String electro_box_material; - G4String s3_mount_material; G4String small_bolt_material; G4String large_bolt_material; G4String shield_cover_material; G4String magnet_cover_material; - G4String detector_mount_material; - G4String annular_clamp_material; - + G4String cold_finger_material; //------------------------- // Dimensions: @@ -324,45 +317,22 @@ private: G4double electrobox_lip_length; G4double electrobox_lip_inner_radius; G4double electrobox_z_offset; - + // ------------------------- - // Dimensions of Si-CD Mount + // Dimensions of ColdFinger // ------------------------- - G4double s3_mount_length; - G4double s3_mount_thickness; - G4double s3_active_radius; - G4double s3_mount_chamfer; - G4double s3_mount_centre_to_chamfer; - G4double s3_mount_angular_offset; - G4double s3_mount_z_offset; - + G4double coldfinger_thickness; + G4double coldfinger_length; + G4double coldfinger_width; + G4double coldfinger_z_offset; + G4double coldfinger_hole_radius; + // ------------------------------ // Dimensions of Plastic Coatings // ------------------------------ G4double magnet_coating_thickness; G4double shield_coating_thickness; - // ---------------------------- - // Dimensions of Detector Mount - // ---------------------------- - G4double detector_mount_length; - G4double detector_mount_width; - G4double detector_mount_thickness; - G4double detector_mount_inner_radius; - G4double detector_mount_lip_radius; - G4double detector_mount_lip_thickness; - G4double detector_mount_angular_offset; - G4double detector_to_target_distance; - G4double detector_thickness; - - // --------------------------- - // Dimensions of Annular Clamp - // --------------------------- - G4double annular_clamp_thickness; - G4double annular_clamp_length; - G4double annular_clamp_width; - G4double annular_clamp_plane_offset; - // --------------------------- // individual offsets for visualisation // --------------------------- @@ -382,16 +352,13 @@ private: G4int targetWheelCopyNumber; G4int biasPlateCopyNumber; G4int gearCopyNumber; - G4int targetMountCopyNumber; G4int gearStickCopyNumber; G4int electroBoxCopyNumber; - G4int s3MountCopyNumber; G4int photonShieldClampBoltCopyNumber; G4int shieldCoveringCopyNumber; G4int magnetCoveringCopyNumber; - G4int detectorMountCopyNumber; - G4int annularClampCopyNumber; - + G4int coldFingerCopyNumber; + private: ////////////////////////////////////////////////////// // internal methods and functions in ApparatusSpiceTargetChamber::Build() @@ -405,8 +372,8 @@ private: void BuildTargetWheelGears(); void BuildTargetWheelGearPlates(); void BuildGearStick(); - void BuildTargetMountPlate(); void BuildBiasPlate(); + void BuildTargetMountPlate(); void BuildPhotonShield(); void BuildPhotonShieldClamps(); void BuildPhotonShieldClampBolts(); @@ -415,10 +382,9 @@ private: void BuildMagnetClampChamber(); void BuildMagnetClampPhotonShield(); void BuildElectroBox(); - void BuildS3Mount(); void BuildShieldCovering(); - void BuildDetectorMount(); - void BuildAnnularClamps(); + void BuildColdFinger(); + void PlaceTargetChamberFrontRing(); void PlaceTargetChamberSphere(); void PlaceTargetChamberCylinderDownstream(); @@ -436,11 +402,9 @@ private: void PlaceMagnetClampChamber(G4int); void PlaceMagnetClampPhotonShield(G4int); void PlaceElectroBox(); - void PlaceS3Mount(); void PlaceShieldCovering(); - void PlaceDetectorMount(); - void PlaceAnnularClamps(); - + void PlaceColdFinger(); + // functions G4RotationMatrix* RotateMagnets(G4int); G4ThreeVector TranslateMagnets(G4int,G4double,G4double); diff --git a/include/DetectionSystemNew.hh b/include/DetectionSystemNew.hh new file mode 100644 index 0000000..0d5cc43 --- /dev/null +++ b/include/DetectionSystemNew.hh @@ -0,0 +1,194 @@ + +#ifndef DetectionSystemNew_h +#define DetectionSystemNew_h 1 + +#include "G4VUserDetectorConstruction.hh" +#include "globals.hh" + +class G4LogicalVolume; +class G4VPhysicalVolume; +class G4Material; +class DetectorMessenger; +//class MagneticField; + + + +// define custom colours for visualisation +#define CHAMBER_COL 0.5,0.5,0.0 +#define LAYER1_COL 0.7,0.0,0.8 +#define LAYER2_COL 0.7,0.2,0.8 +#define LAYER3_COL 0.7,0.4,0.8 +#define NDFEB_COL 0.7,0.1,0.1 +#define MAGCOV_COL 0.3,0.3,0.3 + + + +class DetectionSystemNew +{ + + + public: + DetectionSystemNew(); + ~DetectionSystemNew(); + + + public: + void Build(G4LogicalVolume*); + + + private: + G4LogicalVolume* expHallLog; + + + private: + + //////////////////////////////////////////////// + // Logical Volumes used in NewDetectionSystem // + //////////////////////////////////////////////// + + G4LogicalVolume* final_chamber_log; + G4LogicalVolume* magnet_log; + G4LogicalVolume* magnet_cover_log; + G4LogicalVolume* photon_shield_layer_one_log; + G4LogicalVolume* photon_shield_layer_two_log; + G4LogicalVolume* photon_shield_layer_three_log; + + private: + + ///////////////////////////////////////////////// + // Physical Volumes used in NewDetectionSystem // + ///////////////////////////////////////////////// + + G4VPhysicalVolume* final_chamber_phys; + G4VPhysicalVolume* magnet_phys; + G4VPhysicalVolume* magnet_cover_phys; + G4VPhysicalVolume* photon_shield_layer_one_phys; + G4VPhysicalVolume* photon_shield_layer_two_phys; + G4VPhysicalVolume* photon_shield_layer_three_phys; + + private: + + /////////////////////////////////////////// + // Properties used in NewDetectionSystem // + /////////////////////////////////////////// + + //////////////////////// + // Magnet properties: // + //////////////////////// + + + G4int NUMBER_OF_MAGNETS; + G4int detector_alignment; + G4int no_magnet_layer; + + + //////////////// + // Materials: // + //////////////// + + G4String magnet_material; + G4String magnet_cover_material; + G4String chamber_material; + G4String photon_shield_layer_one_material; + G4String photon_shield_layer_two_material; + G4String photon_shield_layer_three_material; + + ///////////////// + // Dimensions: // + ///////////////// + + ////////////// + // Chamber: // + ////////////// + + G4double chamber_hemisphere_inner_radius; + G4double chamber_hemisphere_outer_radius; + G4double chamber_cylinder_inner_radius; + G4double chamber_cylinder_outer_radius; + G4double chamber_cylinder_height; + G4double chamber_cut_out_inner_radius; + G4double chamber_cut_out_outer_radius; + G4double chamber_cut_out_height; + + //////////////////// + // Photon Shield: // + //////////////////// + + G4double photon_shield_layer_one_thickness; + G4double photon_shield_layer_two_thickness; + G4double photon_shield_layer_three_thickness; + + G4double photon_shield_target_distance ; + G4double detector_full_width ; + G4double detector_target_distance; + + G4double photon_shield_beam_cut_outer_radius; + G4double photon_shield_beam_cut_inner_radius; + G4double photon_shield_beam_cut_height; + + G4double photon_shield_rotation; + G4double photon_shield_cut_angle; + + + ////////////// + // Magnets: // + ////////////// + + G4double plate_one_thickness; + G4double plate_one_length; + G4double plate_one_height; + G4double plate_one_one_lower_height; + G4double plate_one_edge_x; + G4double plate_two_thickness; + G4double plate_two_length; + G4double plate_two_height; + + G4double plate_two_one_thickness; + G4double plate_two_one_length; + G4double plate_two_one_height; + G4double plate_two_one_edge_x; + G4double plate_two_two_thickness; + G4double plate_two_two_length; + G4double plate_two_two_height; + + G4double plate_three_one_thickness; + G4double plate_three_one_length; + G4double plate_three_one_height; + G4double plate_three_one_edge_x; + G4double plate_three_two_thickness; + G4double plate_three_two_length; + G4double plate_three_two_height; + + G4double distance_from_target; + G4double top_cutting_box_angle; + G4double bottom_cutting_box_angle; + G4double top_cutting_box_lowering_height; + G4double bottom_cutting_box_raising_height; + G4double magnet_cover_thickness; + + + private: + + /////////////////////////////////////////////////////////////////// + // internal methods and functions in NewDetectionSystem::Build() // + /////////////////////////////////////////////////////////////////// + + // methods + void BuildChamber(); + void BuildPhotonShield(); + void BuildMagnet(); + void BuildMagnetCover(); + + void PlaceChamber(); + void PlacePhotonShield(); + void PlaceMagnet(G4int); + void PlaceMagnetCover(G4int); + + // functions + G4RotationMatrix* RotateMagnets(G4int); + G4ThreeVector TranslateMagnets(G4int,G4double,G4double); + + +}; + +#endif diff --git a/include/DetectionSystemS3.hh b/include/DetectionSystemS3.hh index 54379d0..cb83bb7 100755 --- a/include/DetectionSystemS3.hh +++ b/include/DetectionSystemS3.hh @@ -39,6 +39,8 @@ #include "globals.hh" #define AL_COL 0.5,0.5,0.5 +#define PEEK_COL 0.5, 0.5, 0.0 + class DetectionSystemS3 { @@ -52,26 +54,24 @@ public: private: G4AssemblyVolume* assembly; G4AssemblyVolume* assemblyS3Ring[24]; + G4VPhysicalVolume* s3_mount_phys; -public: - G4int Build(); - G4int PlaceDetector(G4LogicalVolume* exp_hall_log, G4ThreeVector move, - G4int ringNumber, G4int nRadSeg, G4int detectorNumber); - G4int PlaceGuardRing(G4LogicalVolume* exp_hall_log, G4ThreeVector move); - + private: G4ThreeVector GetDirectionXYZ(G4double theta, G4double phi); G4LogicalVolume* S3InnerGuardRing_log; G4LogicalVolume* S3OuterGuardRing_log; G4LogicalVolume* siDetS3Ring_log[24]; - + G4LogicalVolume* s3_mount_log; + //--------------------------------------------------------// // SPICE physical properties // OBS: crystal properties are public, others are private //--------------------------------------------------------// private: G4String wafer_material; + G4String s3_mount_material; //-----------------------------// // parameters for the annular // @@ -90,7 +90,21 @@ public: private: G4double S3DetGuardRingInnerDiameter; G4double S3DetGuardRingOuterDiameter; - + + // ------------------------- + // Dimensions of Si-CD Mount + // ------------------------- + G4double s3_mount_length; + G4double s3_mount_thickness; + G4double s3_active_radius; + G4double s3_mount_chamfer; + G4double s3_mount_centre_to_chamfer; + + // ------------------------- + // Copy number + // ------------------------- + G4int s3MountCopyNumber; + //------------------------------------------------// // internal methods in Build() //------------------------------------------------// @@ -98,8 +112,17 @@ private: G4int BuildSiliconWafer(G4int ringID); G4int BuildInnerGuardRing(); G4int BuildOuterGuardRing(); - - G4Tubs* BuildCrystal(G4int myRingID); + void BuildS3Mount(); + + + G4Tubs* BuildCrystal(G4int myRingID); + + public: + G4int Build(); + G4int PlaceDetector( G4LogicalVolume* exp_hall_log, G4ThreeVector move, G4double angle_offset, G4int ringNumber, G4int nRadSeg, G4int detectorNumber); + G4int PlaceGuardRing(G4LogicalVolume* exp_hall_log, G4ThreeVector move); + void PlaceS3Mount(G4LogicalVolume* exp_hall_log, G4ThreeVector move, G4double angle_offset); + }; #endif diff --git a/include/DetectionSystemSpice.hh b/include/DetectionSystemSpice.hh index 23b5a45..1bdba43 100755 --- a/include/DetectionSystemSpice.hh +++ b/include/DetectionSystemSpice.hh @@ -39,6 +39,7 @@ #include "globals.hh" #define AL_COL 0.5,0.5,0.5 +#define PEEK_COL 0.5, 0.5, 0.0 class DetectionSystemSpice { @@ -59,23 +60,57 @@ public: G4int PlaceDetector(G4LogicalVolume* exp_hall_log, G4ThreeVector move, G4int ringNumber, G4int nRadSeg, G4int detectorNumber); G4int PlaceGuardRing(G4LogicalVolume* exp_hall_log, G4ThreeVector move); - + void PlaceDetectorMount(G4LogicalVolume* exp_hall_log, G4ThreeVector move); + void PlaceAnnularClamps(G4LogicalVolume* exp_hall_log, G4ThreeVector move); private: G4ThreeVector GetDirectionXYZ(G4double theta, G4double phi); - + + G4LogicalVolume* detector_mount_log; + G4LogicalVolume* annular_clamp_log; G4LogicalVolume* siInnerGuardRing_log; G4LogicalVolume* siOuterGuardRing_log; G4LogicalVolume* siDetSpiceRing_log[10]; - + G4VPhysicalVolume* detector_mount_phys; + G4VPhysicalVolume* annular_clamp_phys; + //--------------------------------------------------------// // SPICE physical properties // OBS: crystal properties are public, others are private //--------------------------------------------------------// private: G4String wafer_material; - + G4String detector_mount_material; + G4String annular_clamp_material; + + // ---------------------------- + // Dimensions of Detector Mount + // ---------------------------- + G4double detector_mount_length; + G4double detector_mount_width; + G4double detector_mount_thickness; + G4double detector_mount_inner_radius; + G4double detector_mount_lip_radius; + G4double detector_mount_lip_thickness; + G4double detector_mount_angular_offset; + G4double detector_to_target_distance; + G4double detector_thickness; + + // --------------------------- + // Dimensions of Annular Clamp + // --------------------------- + G4double annular_clamp_thickness; + G4double annular_clamp_length; + G4double annular_clamp_width; + G4double annular_clamp_plane_offset; + + //----------------------------- + // copy numbers + //----------------------------- + G4int detectorMountCopyNumber; + G4int annularClampCopyNumber; + //-----------------------------// // parameters for the annular // // planar detector crystal // @@ -101,9 +136,11 @@ private: G4int BuildSiliconWafer(G4int ringID); G4int BuildInnerGuardRing(); G4int BuildOuterGuardRing(); - + void BuildDetectorMount(); + void BuildAnnularClamps(); + G4Tubs* BuildCrystal(G4int myRingID); - + }; #endif diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh index 6511804..da795bf 100755 --- a/include/DetectorConstruction.hh +++ b/include/DetectorConstruction.hh @@ -57,6 +57,7 @@ class DetectionSystemS3; class DetectionSystemPaces; class DetectionSystemSodiumIodide; class DetectionSystemLanthanumBromide; +class DetectionSystemNew; class DetectionSystemBox; @@ -77,7 +78,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction void SetWorldDimensions( G4ThreeVector ); void SetWorldVis( G4bool ); void SetWorldMagneticField( G4ThreeVector ); - void SetTabMagneticField(G4String); + void SetTabMagneticField(G4String, G4double, G4double); void SetGenericTargetMaterial( G4String ); void SetGenericTargetDimensions( G4ThreeVector ); @@ -135,10 +136,13 @@ class DetectorConstruction : public G4VUserDetectorConstruction void AddDetectionSystemSceptar(G4int ndet); void AddDetectionSystemPaces(G4int ndet); + + void AddDetectionSystemNew(G4int ndet); + void AddNewSquareDetector(G4int ndet); void SetSpiceResolutionVariables(G4double intercept, G4double gain); void AddDetectionSystemSpice(G4int nRings); - void AddDetectionSystemS3(G4int nRings); + void AddDetectionSystemS3(G4int nRings, G4double posX, G4double posY, G4int posZ, G4double AngleOffset); void UseTIGRESSPositions( G4bool input ) {useTigressPositions = input;}; private: diff --git a/include/DetectorMessenger.hh b/include/DetectorMessenger.hh index b045990..de0b29b 100755 --- a/include/DetectorMessenger.hh +++ b/include/DetectorMessenger.hh @@ -110,12 +110,15 @@ class DetectorMessenger: public G4UImessenger G4UIcmdWithAnInteger* AddDetectionSystemGriffinForwardDetectorCmd; G4UIcmdWithAnInteger* AddDetectionSystemGriffinBackCmd; G4UIcmdWithAnInteger* AddDetectionSystemGriffinBackDetectorCmd; - //G4UIcmdWith3Vector* AddDetectionSystemGriffinPositionConfigCmd; - G4UIcommand* SetSpiceResolutionVariablesCmd; + //G4UIcmdWith3Vector* AddDetectionSystemGriffinPositionConfigCmd; + G4UIcommand* SetSpiceResolutionVariablesCmd; G4UIcmdWithAnInteger* AddDetectionSystemSpiceCmd; - G4UIcmdWithAnInteger* AddDetectionSystemS3Cmd; + G4UIcommand* AddDetectionSystemS3Cmd; G4UIcmdWithAnInteger* AddDetectionSystemPacesCmd; + G4UIcmdWithAnInteger* AddDetectionSystemNewCmd; + G4UIcmdWithAnInteger* AddNewSquareDetectorCmd; G4UIcmdWithAnInteger* AddDetectionSystemGriffinHevimetCmd ; + G4UIcmdWithAnInteger* AddDetectionSystemGriffinCustomDetectorCmd ; G4UIcmdWithAnInteger* AddDetectionSystemGriffinCustomCmd ; diff --git a/include/EventAction.hh b/include/EventAction.hh index 691ef3f..ea88eb7 100755 --- a/include/EventAction.hh +++ b/include/EventAction.hh @@ -31,6 +31,10 @@ #include "globals.hh" #include "HistoManager.hh" +#include "TrackInformation.hh" + +//c++ +#include class RunAction; class HistoManager; @@ -44,7 +48,9 @@ static const int NUMSTEPVARS = 16; class EventAction : public G4UserEventAction { + public: + EventAction(RunAction*, HistoManager*); virtual ~EventAction(); @@ -52,39 +58,12 @@ public: void EndOfEventAction(const G4Event*); G4int GetEventNumber(){return evtNb;}; - + void AddStepTracker(G4double eventNumber, G4double stepNumber, G4String volume, G4double cryNumber, G4double detNumber, G4double depEnergy, G4double posx, G4double posy, G4double posz, G4double time, G4double originDirectionX, G4double originDirectionY, G4double originDirectionZ, - G4double originEnergy, G4int originPdg, G4int originID, G4int trackID) { - - if(histoManager->GetStepTrackerBool()) { - stepTracker[0][stepIndex] = eventNumber; - stepTracker[1][stepIndex] = stepNumber; - stepTracker[2][stepIndex] = cryNumber; - stepTracker[3][stepIndex] = detNumber; - stepTracker[4][stepIndex] = depEnergy; - stepTracker[5][stepIndex] = posx; - stepTracker[6][stepIndex] = posy; - stepTracker[7][stepIndex] = posz; - stepTracker[8][stepIndex] = time; - stepTracker[9][stepIndex] = originDirectionX; - stepTracker[10][stepIndex] = originDirectionY; - stepTracker[11][stepIndex] = originDirectionZ; - stepTracker[12][stepIndex] = originEnergy; - stepTracker[13][stepIndex] = originPdg; - stepTracker[14][stepIndex] = originID; - stepTracker[15][stepIndex] = trackID; - - stepVolume[stepIndex] = volume ; - stepIndex++; - if(stepIndex == MAXSTEPS) { - G4cout << "\n ----> error 13423549 \n" << G4endl; - exit(1); - } - }; - }; + G4double originEnergy, G4int originPdg, G4int originID, G4int trackID) ; // particle types void AddParticleType(G4int index) {particleTypes[index] += 1;}; @@ -113,7 +92,9 @@ public: void AddSpiceCrystDet(G4double de, G4double dl, G4int det) {SpiceCrystEnergyDet[det] += de; SpiceCrystTrackDet[det] += dl;} ; void AddS3CrystDet(G4double de, G4double dl, G4int det) {SpiceCrystEnergyDet[det] += de; SpiceCrystTrackDet[det] += dl;} ; void AddPacesCrystDet(G4double de, G4double dl, G4int det) {PacesCrystEnergyDet[det] += de; PacesCrystTrackDet[det] += dl;} ; + void AddNewCrystDet(G4double de, G4double dl, G4int det) {NewCrystEnergyDet[det] += de; NewCrystTrackDet[det] += dl;} ; + void SetPrimaryInfo(TrackInformation* info); private: @@ -130,6 +111,7 @@ private: void FillS3Cryst(); void FillPacesCryst() ; void Fill8piCryst() ; + void FillNewCryst(); RunAction* runAct; HistoManager* histoManager; @@ -142,6 +124,7 @@ private: G4double stepTracker[NUMSTEPVARS][MAXSTEPS]; G4String stepVolume[MAXSTEPS]; // volume at each step G4int stepIndex; + vector PrimaryInfo; // Particle types in simulation G4int particleTypes[NUMPARTICLETYPES]; @@ -195,6 +178,9 @@ private: G4double PacesCrystEnergyDet[MAXNUMDET] ; G4double PacesCrystTrackDet[MAXNUMDET] ; + + G4double NewCrystEnergyDet[MAXNUMDET] ; + G4double NewCrystTrackDet[MAXNUMDET] ; }; diff --git a/include/HistoManager.hh b/include/HistoManager.hh index c4d1b45..3ff3c21 100755 --- a/include/HistoManager.hh +++ b/include/HistoManager.hh @@ -456,7 +456,30 @@ enum HISTONAME paces_crystal_edep_det16, paces_crystal_edep_det17, paces_crystal_edep_det18, - paces_crystal_edep_det19 + paces_crystal_edep_det19, + new_crystal_edep, +// This is probably the wrong number of detectors for new. + new_crystal_edep_sum, + new_crystal_edep_det0, + new_crystal_edep_det1, + new_crystal_edep_det2, + new_crystal_edep_det3, + new_crystal_edep_det4, + new_crystal_edep_det5, + new_crystal_edep_det6, + new_crystal_edep_det7, + new_crystal_edep_det8, + new_crystal_edep_det9, + new_crystal_edep_det10, + new_crystal_edep_det11, + new_crystal_edep_det12, + new_crystal_edep_det13, + new_crystal_edep_det14, + new_crystal_edep_det15, + new_crystal_edep_det16, + new_crystal_edep_det17, + new_crystal_edep_det18, + new_crystal_edep_det19 }; diff --git a/include/NewSquareDetector.hh b/include/NewSquareDetector.hh new file mode 100644 index 0000000..7ff9b35 --- /dev/null +++ b/include/NewSquareDetector.hh @@ -0,0 +1,94 @@ +#ifndef NewSquareDetector_h +#define NewSquareDetector_h 1 + +#include "G4VUserDetectorConstruction.hh" +#include "globals.hh" + +#define SIL_COL 0.1,0.9,0.9 +#define GUARD_COL 0.1,0.1,0.9 + +class NewSquareDetector +{ + + public: + NewSquareDetector(); + ~NewSquareDetector(); + + ////////////////////////////////// + // logical and physical volumes // + ////////////////////////////////// + + private: + G4AssemblyVolume* assembly; + G4AssemblyVolume* assemblySquareDet[5]; + G4AssemblyVolume* assemblyGuardRing[5]; + + + public: + G4int Build(); + G4int PlaceDetector(G4LogicalVolume* exp_hall_log, G4int detector); + G4int PlaceGuardRing(G4LogicalVolume* exp_hall_log, G4int detector); + + + private: + G4LogicalVolume* expHallLog; + + private: + G4LogicalVolume* square_guard_ring_log[5]; + G4LogicalVolume* SquareDetect_log[5]; + + + //////////////////////////////////////////////////////////// + // OBS: crystal properties are public, others are private // + //////////////////////////////////////////////////////////// + + private: + G4String detector_material; + + //////////////////////////////////////////////// + // parameters for the planar detector crystal // + //////////////////////////////////////////////// + + public: + G4double square_segment_element_length; + G4double square_segment_element_width; + G4double square_segment_thickness; + G4double no_x_segments; + G4double no_y_segments; + G4double no_detectors; + G4double SquareDetectorDistance; + G4double BeamPipeXDistanceGR; + G4double BeamPipeYDistanceGR; + G4double BeamPipeXDistanceD; + G4double BeamPipeYDistanceD; + G4int detector_alignment; + G4double SquareDetectorRotationXY; + G4double SquareDetectorRotationZ; + G4double SquareDetectorCorrection; + + G4ThreeVector TranslateDetectors(G4int,G4double,G4double); + + + /////////////////////////////////// + // parameters for the guard ring // + /////////////////////////////////// + + private: + G4double square_guard_ring_depth; + G4double square_guard_ring_thickness; + + ///////////////////////////////// + // internal methods in Build() // + ///////////////////////////////// + + private: + G4int BuildDetectorFace1(G4int detectorID); + G4int BuildGuardRing1(G4int detectorID); + G4Box* BuildSegment1(); + G4int BuildDetectorFace2(G4int detectorID); + G4int BuildGuardRing2(G4int detectorID); + G4Box* BuildSegment2(); + +}; + +#endif diff --git a/include/PrimaryGeneratorMessenger.hh b/include/PrimaryGeneratorMessenger.hh index c0ef1e8..37ca277 100755 --- a/include/PrimaryGeneratorMessenger.hh +++ b/include/PrimaryGeneratorMessenger.hh @@ -70,7 +70,7 @@ class PrimaryGeneratorMessenger: public G4UImessenger G4UIcmdWithADoubleAndUnit* radiusCmd; G4UIcmdWithAString* particleCmd; G4UIcommand* ionCmd; - G4UIcommand* rangeCmd; + G4UIcommand* rangeCmd; G4UIcmdWithAString* betaPlusEmissionCmd; G4UIcmdWithAString* betaMinusEmissionCmd; G4UIcmdWithAString* radioactiveBetaDecayCmd; diff --git a/include/RawG4Event.hh b/include/RawG4Event.hh index 0f7daa0..d1eeadb 100644 --- a/include/RawG4Event.hh +++ b/include/RawG4Event.hh @@ -3,6 +3,7 @@ //Root #include "TVector3.h" +#include "TString.h" //c++ #include "vector" @@ -14,6 +15,7 @@ using namespace std; //g4 #include "globals.hh" + class RawG4Event { @@ -75,6 +77,8 @@ void FillVectors(int pdg, // particle pdg double PrimEnergy,//original(primary) energy double Mx, double My, double Mz); // primary particle momentum vector + + /************************* G E T T E R S ************************/ @@ -87,6 +91,7 @@ TVector3 GetHCPrimaryMomentum(int i); Double_t GetHCPrimaryTheta(int i); Double_t GetHCPrimaryPhi(int i); + // after treatment Double_t GetFullEnergy(void) ; // full energy in pad TVector3 GetFirstHitPosition(void) ; // first hit position @@ -110,6 +115,7 @@ Double_t GetPrimaryEnergyForID(int ID) ; // for a specific type of particle PDG Double_t GetDetectedEnergyForID(int ID) ; // full energy in pad for a specific type of particle PDG Double_t GetAngleOfEmissionForID(int ID) ; // for a specific type of particle PDG // from the moment + /* // Per particle type (PDG) diff --git a/include/RootManager.hh b/include/RootManager.hh index 87a6232..3fbe5c8 100644 --- a/include/RootManager.hh +++ b/include/RootManager.hh @@ -9,24 +9,32 @@ #include #include #include +#include using namespace std ; +// CLHEP +#include "CLHEP/Random/RandGauss.h" + +//G4 +#include "TrackInformation.hh" + //ROOT #include "TFile.h" #include "TH1.h" #include "TTree.h" -// CLHEP -#include "CLHEP/Random/RandGauss.h" - //User #include "RawG4Event.hh" - #include "../dataRootClass/TTigFragment.h" #include "../dataRootClass/TSpiceData.h" #include "../dataRootClass/TS3Data.h" #include "../dataRootClass/TGriffinData.h" +#include "../dataRootClass/TPacesData.h" +#include "../dataRootClass/THistoryData.h" +#include "../dataRootClass/TNewData.h" +#include "../dataRootClass/TSceptarData.h" + class RootManager { @@ -51,13 +59,15 @@ class RootManager { map fGeantEvent; //Writing Class for detectors goes here - TTigFragment* fFragment; - - TSpiceData* fSpiceData; - TS3Data* fS3Data; - TGriffinData* fGriffinData; - - + TSpiceData* fSpiceData; + TS3Data* fS3Data; + TPacesData* fPacesData; + TGriffinData* fGriffinData; + THistoryData* fHistoryData; + TNewData* fNewData; + TSceptarData* fSceptarData; + TTigFragment* fFragment; + public: // fill the histograms void FillHist(double); @@ -74,21 +84,27 @@ class RootManager { double,//original (primary) particle energy double, double, double);// primary particle momentum vector + // clear all vectors + void ClearVariables(void); + //Build the mnemonic used in TRIUMF string BuildMnemonic(string volume, int detector, int crystal); - - //Set event number - void SetEventNumber(int) ; // Set the branches on the tree void SetTree(); //SortEvent - void SortEvent(); + void SortEvent(int eventNb); //Set the data in Spice writing Class - void SetSpiceEvent(string mnemonic, int ring, int seg); - void SetS3Event(string mnemonic, int ring, int seg); + void SetHistory( vector info ); // History of the LAST particles in a cascade of events (Whether a part or all of this cascade ended in the detector or ) + void SetSpiceEvent(int eventNb, string mnemonic, int Ring, int Seg); + void SetS3Event(int eventNb, string mnemonic, int Ring, int Seg); + void SetPacesEvent(int eventNb, string mnemonic, int Ring, int Seg); + void SetNewEvent(int eventNb, string mnemonic, int detector, int Seg); + void SetSceptarEvent(int eventNb, string mnemonic, int detector, int Seg); + + void SetGriffinEvent(int key); void SetFragmentEvent(string key); diff --git a/include/SteppingAction.hh b/include/SteppingAction.hh index eca6393..087e827 100755 --- a/include/SteppingAction.hh +++ b/include/SteppingAction.hh @@ -60,11 +60,20 @@ private: void SetDetAndCryNumberForDeadLayerSpecificGriffinCrystal(G4String); void SetDetNumberForGenericDetector( G4String ); + //Paces + void SetDetAndCryNumberForPacesDetector( G4String ); + //Spice void SetDetAndCryNumberForSpiceDetector( G4String ) ; //S3 void SetDetAndCryNumberForS3Detector( G4String ) ; + + //New + void SetDetAndCryNumberForNewDetector( G4String ) ; + + //Sceptar + void SetDetAndCryNumberForSceptarDetector( G4String ) ; G4int FindTrueGriffinDetector(G4int); diff --git a/include/TabulatedMagneticField.hh b/include/TabulatedMagneticField.hh index 58b249f..a8c2dca 100644 --- a/include/TabulatedMagneticField.hh +++ b/include/TabulatedMagneticField.hh @@ -39,10 +39,11 @@ class TabulatedMagneticField // The physical extent of the defined region double dx, dy, dz; double fZoffset; + double fZrotation; bool invertX, invertY, invertZ; public: - TabulatedMagneticField(const char* filename, G4double zOffset ); + TabulatedMagneticField(const char* filename, G4double zOffset, G4double zRotation ); void GetFieldValue( const double Point[4], double *Bfield ) const; }; diff --git a/include/TrackInformation.hh b/include/TrackInformation.hh index 22a54f1..e9ddd13 100644 --- a/include/TrackInformation.hh +++ b/include/TrackInformation.hh @@ -5,8 +5,10 @@ #include "G4ThreeVector.hh" #include "G4ParticleDefinition.hh" #include "G4Track.hh" +#include "G4Step.hh" #include "G4Allocator.hh" #include "G4VUserTrackInformation.hh" +#include "G4VProcess.hh" //c++ #include @@ -27,47 +29,142 @@ class TrackInformation : public G4VUserTrackInformation {return (this==&right);} void Print() const; - + void Clear() ; + void PartialClear() ; + private: + // Original particle at source + G4int originalParentID; G4int originalTrackID; G4int originalPdg; G4ThreeVector originalPosition; G4ThreeVector originalMomentum; G4double originalEnergy; G4double originalTime; + // At first impact + G4String originalImpactVolume ; + G4ThreeVector originalImpactPosition; + G4ThreeVector originalImpactMomentum; + +// secondaries + G4String currentProcess; + G4int currentTrackID; + G4int currentparentTrackID; + G4int currentPdg; + G4ThreeVector currentPositionAtVertex; + G4ThreeVector currentPositionAtDetector; + G4ThreeVector currentPositionAtDeath; // death or out of world + G4ThreeVector currentMomentumAtDeath; // death or out of world + + //At first impact + G4String currentImpactVolume ; + G4ThreeVector currentImpactPosition; + G4ThreeVector currentImpactMomentum; + G4ThreeVector currentMomentumAtVertex; // contains the angle at the emission + G4ThreeVector currentMomentumAtDetector; // contains angle at the impact on the detector - vector AncestorsBirthVolume; // MHD 02 May 2013 - vector AncestorsDeathVolume; // MHD 02 May 2013 - vector AncestorsPdg; // MHD 02 May 2013 - + G4double currentEnergyAtVertex; + G4double currentTimeAtVertex; + + vector SecondariesProcess; + vector SecondariesBirthVolume; + vector SecondariesDeathVolume; + vector SecondariesPdg; + + G4bool Tagged; + public: + //Getters - inline G4int GetOriginalTrackID() const {return originalTrackID;} - inline G4int GetOriginalPdg() const {return originalPdg;} - inline G4ThreeVector GetOriginalPosition() const {return originalPosition;} - inline G4ThreeVector GetOriginalMomentum() const {return originalMomentum;} - inline G4double GetOriginalEnergy() const {return originalEnergy;} - inline G4double GetOriginalTime() const {return originalTime;} - inline vector GetAncestorsBirthVolume() const {return AncestorsBirthVolume;} - inline vector GetAncestorsDeathVolume() const {return AncestorsDeathVolume;} - inline vector GetAncestorsPdg() const {return AncestorsPdg;} - + //Original particle at source + inline G4bool GetTagged(void) const {return Tagged;} + inline G4int GetOriginalParentID(void) const {return originalParentID;} + inline G4int GetOriginalTrackID(void) const {return originalTrackID;} + inline G4int GetOriginalPdg(void) const {return originalPdg;} + inline G4ThreeVector GetOriginalPosition(void) const {return originalPosition;} + inline G4ThreeVector GetOriginalMomentum(void) const {return originalMomentum;} + inline G4double GetOriginalEnergy(void) const {return originalEnergy;} + inline G4double GetOriginalTime(void) const {return originalTime;} + + inline G4String GetOriginalImpactVolume(void) { return originalImpactVolume; } + inline G4ThreeVector GetOriginalImpactMomentum(void) { return originalImpactMomentum; } + inline G4ThreeVector GetOriginalImpactPosition(void) { return originalImpactPosition; } + +// Secondaries + inline G4String GetCurrentProcess(void) const {return currentProcess;} + inline G4int GetCurrentParentID(void) const {return currentparentTrackID;} + inline G4int GetCurrentTrackID(void) const {return currentTrackID;} + inline G4int GetCurrentPdg(void) const {return currentPdg;} + inline G4ThreeVector GetCurrentPositionAtVertex(void) const {return currentPositionAtVertex;} + inline G4ThreeVector GetCurrentMomentumAtVertex(void) const {return currentMomentumAtVertex;} + inline G4ThreeVector GetCurrentPositionAtDetector(void) const {return currentPositionAtDetector;} + inline G4ThreeVector GetCurrentMomentumAtDetector(void) const {return currentMomentumAtDetector;} + inline G4ThreeVector GetCurrentPositionAtDeath(void) const {return currentPositionAtDeath;} + inline G4ThreeVector GetCurrentMomentumAtDeath(void) const {return currentMomentumAtDeath;} + inline G4double GetCurrentEnergyAtVertex(void) const {return currentEnergyAtVertex;} + inline G4double GetCurrentTimeAtVertex(void) const {return currentTimeAtVertex;} + // all the vector + inline vector GetSecondariesProcess(void) const {return SecondariesProcess;} + inline vector GetSecondariesBirthVolume(void) const {return SecondariesBirthVolume;} + inline vector GetSecondariesDeathVolume(void) const {return SecondariesDeathVolume;} + inline vector GetSecondariesPdg(void) const {return SecondariesPdg;} + // 1 element + inline G4String GetSecondariesProcessAt(unsigned i) const {return SecondariesProcess.at(i);} + inline G4String GetSecondariesBirthVolumeAt(unsigned i) const {return SecondariesBirthVolume.at(i);} + inline G4String GetSecondariesDeathVolumeAt(unsigned i) const {return SecondariesDeathVolume.at(i);} + inline G4int GetSecondariesPdgAt(unsigned i) const {return SecondariesPdg.at(i);} + // size of vectors + inline G4int GetSecondariesProcessSize(void) const {return SecondariesProcess.size();} + inline G4int GetSecondariesBirthVolumeSize(void) const {return SecondariesBirthVolume.size();} + inline G4int GetSecondariesDeathVolumeSize(void) const {return SecondariesDeathVolume.size();} + inline G4int GetSecondariesPdgSize(void) const {return SecondariesPdg.size();} + + inline G4String GetCurrentImpactVolume(void) { return currentImpactVolume; } + inline G4ThreeVector GetCurrentImpactMomentum(void) { return currentImpactMomentum; } + inline G4ThreeVector GetCurrentImpactPosition(void) { return currentImpactPosition; } + //Setters - void SetAncestorsBirthVolumeElement(G4String ans_birth_vol) { AncestorsBirthVolume.push_back(ans_birth_vol);} - void SetAncestorsDeathVolumeElement(G4String ans_death_vol) { AncestorsDeathVolume.push_back(ans_death_vol);} - void SetAncestorsPdgElement(G4int ans_pdg) { AncestorsPdg.push_back(ans_pdg);} + inline void SetTagged( G4bool tag) { Tagged = tag ;} + + // original + inline void SetOriginalImpactVolume(G4String vol ) { originalImpactVolume = vol ; } + inline void SetOriginalImpactMomentum(G4ThreeVector momentum) { originalImpactMomentum = momentum;} + inline void SetOriginalImpactPosition(G4ThreeVector position) { originalImpactPosition = position;} + // secondaries + inline void SetCurrentProcess(G4String process) { currentProcess = process ;} + inline void SetCurrentParentID(int id) { currentparentTrackID = id ;} + inline void SetCurrentTrackID(int id) { currentTrackID = id ;} + inline void SetCurrentPdg(int pdg) { currentPdg = pdg ;} + inline void SetCurrentPositionAtVertex(G4ThreeVector position) { currentPositionAtVertex = position;} + inline void SetCurrentMomentumAtVertex(G4ThreeVector momentum) { currentMomentumAtVertex = momentum;} + inline void SetCurrentPositionAtDetector(G4ThreeVector position) { currentPositionAtDetector = position ;} + inline void SetCurrentMomentumAtDetector(G4ThreeVector momentum) { currentMomentumAtDetector = momentum;} + inline void SetCurrentPositionAtDeath(G4ThreeVector position) { currentPositionAtDeath = position ;} + inline void SetCurrentMomentumAtDeath(G4ThreeVector position) { currentMomentumAtDeath = position ;} + inline void SetCurrentEnergyAtVertex(G4double energy ) {currentEnergyAtVertex = energy ;} + inline void SetCurrentTimeAtVertex(G4double time ) {currentTimeAtVertex = time ;} + + inline void SetCurrentImpactVolume(G4String vol ) { currentImpactVolume = vol ; } + inline void SetCurrentImpactMomentum(G4ThreeVector momentum) { currentImpactMomentum = momentum;} + inline void SetCurrentImpactPosition(G4ThreeVector position) { currentImpactPosition = position;} + + inline void SetSecondariesProcessElement(G4String ans_birth_vol) { SecondariesProcess.push_back(ans_birth_vol);} + inline void SetSecondariesBirthVolumeElement(G4String ans_birth_vol) { SecondariesBirthVolume.push_back(ans_birth_vol);} + inline void SetSecondariesDeathVolumeElement(G4String ans_death_vol) { SecondariesDeathVolume.push_back(ans_death_vol);} + inline void SetSecondariesPdgElement(G4int ans_pdg) { SecondariesPdg.push_back(ans_pdg);} + }; extern G4Allocator aTrackInformationAllocator; -inline void* TrackInformation::operator new(size_t) -{ void* aTrackInfo; - aTrackInfo = (void*)aTrackInformationAllocator.MallocSingle(); - return aTrackInfo; +inline void* TrackInformation::operator new(size_t){ + void* aTrackInfo; + aTrackInfo = (void*)aTrackInformationAllocator.MallocSingle(); + return aTrackInfo; } -inline void TrackInformation::operator delete(void *aTrackInfo) -{ aTrackInformationAllocator.FreeSingle((TrackInformation*)aTrackInfo);} +inline void TrackInformation::operator delete(void *aTrackInfo){ + aTrackInformationAllocator.FreeSingle((TrackInformation*)aTrackInfo);} #endif diff --git a/include/TrackingAction.hh b/include/TrackingAction.hh index c7bbd32..253e794 100644 --- a/include/TrackingAction.hh +++ b/include/TrackingAction.hh @@ -29,13 +29,6 @@ class TrackingAction : public G4UserTrackingAction I got this version from another Tip by Makoto Asai http://geant4.slac.stanford.edu/Tips/event/5.html -#ifndef T01TrackingAction_h -#define T01TrackingAction_h 1 - -class G4Track; -#include "G4UserTrackingAction.hh" -#include "globals.hh" - class T01TrackingAction : public G4UserTrackingAction { public: diff --git a/include/nonUniformMagneticField.hh b/include/nonUniformMagneticField.hh index c53dac3..b02cc49 100755 --- a/include/nonUniformMagneticField.hh +++ b/include/nonUniformMagneticField.hh @@ -14,7 +14,7 @@ class nonUniformMagneticField public: - nonUniformMagneticField(const char* fieldName, double zOffset); // field read from file + nonUniformMagneticField(const char* fieldName, double zOffset, double zRotation); // field read from file ~nonUniformMagneticField() ; void SetStepperType( G4int i) { fStepperType = i ; } diff --git a/src/ApparatusSpiceTargetChamber.cc b/src/ApparatusSpiceTargetChamber.cc index 0f5fdc8..8568b6b 100755 --- a/src/ApparatusSpiceTargetChamber.cc +++ b/src/ApparatusSpiceTargetChamber.cc @@ -51,7 +51,7 @@ ApparatusSpiceTargetChamber::ApparatusSpiceTargetChamber() this->photon_shield_layer_one_material = "WTa"; this->photon_shield_layer_two_material = "Tin"; this->photon_shield_layer_three_material = "Copper"; - this->downstream_cone_material = "Aluminum"; + this->downstream_cone_material = "Aluminum"; this->ps_clamp_material = "Aluminum"; this->target_wheel_material = "Aluminum"; this->target_wheel_gear_material_a = "Delrin"; @@ -61,14 +61,12 @@ ApparatusSpiceTargetChamber::ApparatusSpiceTargetChamber() this->bias_plate_material = "Aluminum"; this->gear_stick_material = "Peek"; this->electro_box_material = "Aluminum"; - this->s3_mount_material = "Peek"; //? this->small_bolt_material = "Brass"; this->large_bolt_material = "Titanium"; this->shield_cover_material = "Kapton"; this->magnet_cover_material = "Peek"; - this->detector_mount_material = "Aluminum"; //? - this->annular_clamp_material = "Peek"; // ? - + this->cold_finger_material = "Copper"; + //----------------------------- // Dimensions of Target Chamber //----------------------------- @@ -224,45 +222,13 @@ ApparatusSpiceTargetChamber::ApparatusSpiceTargetChamber() this->electrobox_lip_length = 21.5*mm; this->electrobox_lip_inner_radius = 102*mm; this->electrobox_z_offset = -90*mm; - - // ------------------------- - // Dimensions of Si-CD Mount - // ------------------------- - this->s3_mount_length = 120*mm; - this->s3_mount_thickness = 2.3*mm; - this->s3_active_radius = 35*mm; - this->s3_mount_chamfer = 28.284*mm; - this->s3_mount_centre_to_chamfer = 70.711*mm; - this->s3_mount_angular_offset = 20*deg; - this->s3_mount_z_offset = 21*mm; - + // ------------------------------ // Dimensions of Plastic Coatings // ------------------------------ this->magnet_coating_thickness = 1.25*mm; this->shield_coating_thickness = 0.2*mm; - - // ---------------------------- - // Dimensions of Detector Mount - // ---------------------------- - this->detector_mount_length = 148*mm; - this->detector_mount_width = 130*mm; - this->detector_mount_thickness = 8*mm; - this->detector_mount_inner_radius = 52.8*mm; - this->detector_mount_lip_radius = 48.3*mm; - this->detector_mount_lip_thickness = 1*mm; - this->detector_mount_angular_offset = 0*deg; - this->detector_to_target_distance = 115*mm; - this->detector_thickness = 6*mm; - - // --------------------------- - // Dimensions of Annular Clamp - // --------------------------- - this->annular_clamp_thickness = 2*mm; - this->annular_clamp_length = 19.2*mm; - this->annular_clamp_width = 20*mm; - this->annular_clamp_plane_offset = 47.8*mm; //from beam line to first edge - + // ------------------------------------- // individual offsets for visualisation // ------------------------------------- @@ -271,6 +237,15 @@ ApparatusSpiceTargetChamber::ApparatusSpiceTargetChamber() this->middleRingOffset = 0*mm; this->backDetAndAlBoxOffset = 0*mm; + // ------------------------- + // Dimensions of ColdFinger + // ------------------------- + this->coldfinger_thickness = 5*mm ; // CHECK + this->coldfinger_length = 150*mm ; // CHECK + this->coldfinger_width = 200*mm ; // CHECK + this->coldfinger_z_offset = -122*mm ; // CHECK + this->coldfinger_hole_radius = 10*mm ; // CHECK + } // end ApparatusSpiceTargetChamber ////////////////////////////////////////////////// @@ -305,11 +280,10 @@ ApparatusSpiceTargetChamber::~ApparatusSpiceTargetChamber() delete mclamp_chamber_log; delete mclamp_shield_log; delete electro_box_log; - delete s3_mount_log; delete shield_cover_log; delete magnet_cover_log; - delete detector_mount_log; - delete annular_clamp_log; + delete cold_finger_log; + // physical volumes in ApparatusSpiceTargetChamber delete target_chamber_front_ring_phys; delete target_chamber_front_cone_phys; @@ -335,11 +309,9 @@ ApparatusSpiceTargetChamber::~ApparatusSpiceTargetChamber() delete mclamp_chamber_phys; delete mclamp_shield_phys; delete electro_box_phys; - delete s3_mount_phys; delete shield_cover_phys; delete magnet_cover_phys; - delete detector_mount_phys; - delete annular_clamp_phys; + delete cold_finger_phys; } // end ~ApparatusSpiceTargetChamber @@ -367,12 +339,9 @@ void ApparatusSpiceTargetChamber::Build(G4LogicalVolume* exp_hall_log) BuildMagnetClampChamber(); BuildMagnetClampPhotonShield(); BuildElectroBox(); - BuildS3Mount(); BuildShieldCovering(); BuildMagnetCovering(); - BuildDetectorMount(); - BuildAnnularClamps(); - + BuildColdFinger(); PlaceTargetChamberFrontRing(); PlaceTargetChamberSphere(); @@ -387,11 +356,9 @@ void ApparatusSpiceTargetChamber::Build(G4LogicalVolume* exp_hall_log) //PlaceShieldCovering(); PlacePhotonShieldClamps(); PlaceElectroBox(); - PlaceS3Mount(); - PlaceDetectorMount(); - PlaceAnnularClamps(); - for(G4int copyID=0; copyIDNUMBER_OF_MAGNETS; copyID++) - { + PlaceColdFinger(); + + for(G4int copyID=0; copyIDNUMBER_OF_MAGNETS; copyID++) { PlaceCollectorMagnet(copyID); PlaceMagnetClampChamber(copyID); PlaceMagnetClampPhotonShield(copyID); @@ -852,7 +819,7 @@ void ApparatusSpiceTargetChamber::BuildShieldCovering() vis_att->SetVisibility(true); // ** Build Photon Shield Solid - G4double inner_radius = this->photon_shield_inner_radius; + //G4double inner_radius = this->photon_shield_inner_radius; G4double front_outer_radius = this->photon_shield_front_radius; G4double back_outer_radius = this->photon_shield_back_radius; G4double half_length = this->photon_shield_length/2.; @@ -1325,128 +1292,31 @@ void ApparatusSpiceTargetChamber::BuildElectroBox(){ } // end:BuildElectroBox() -void ApparatusSpiceTargetChamber::BuildS3Mount() { - - // ** Visualisation - G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(PEEK_COL)); - vis_att->SetVisibility(true); - - // ** Dimensions - G4double half_length = this->s3_mount_length/2.; - G4double half_thickness = this->s3_mount_thickness/2.; - // Chamfer - G4double chamfer_half_length = this->s3_mount_chamfer/2.; - // Active Area Cut - G4double active_radius = this->s3_active_radius; - - // ** Shapes - G4Box* s3mount_box = new G4Box("s3_mount_box", half_length, half_length, half_thickness); - G4Box* chamfer_cut = new G4Box("chamber_cut", chamfer_half_length, chamfer_half_length, chamfer_half_length); - G4Tubs* active_cut = new G4Tubs("active_cut", 0, active_radius, chamfer_half_length, 0, 360*deg); - - G4double plane_offset = (this->s3_mount_centre_to_chamfer + chamfer_half_length) / sqrt(2.); - G4ThreeVector trans(plane_offset, plane_offset, 0); - G4RotationMatrix* rotate = new G4RotationMatrix(45*deg, 0, 0); - G4SubtractionSolid* s3_mount0 = new G4SubtractionSolid("s3_mount0", s3mount_box, chamfer_cut, rotate, trans); - trans.setX(-plane_offset); - G4SubtractionSolid* s3_mount1 = new G4SubtractionSolid("s3_mount1", s3_mount0, chamfer_cut, rotate, trans); - trans.setY(-plane_offset); - G4SubtractionSolid* s3_mount2 = new G4SubtractionSolid("s3_mount2", s3_mount1, chamfer_cut, rotate, trans); - trans.setX(plane_offset); - G4SubtractionSolid* s3_mount3 = new G4SubtractionSolid("s3_mount3", s3_mount2, chamfer_cut, rotate, trans); - G4SubtractionSolid* s3_mount = new G4SubtractionSolid("s3_mount", s3_mount3, active_cut); - - // ** Logical - G4Material* s3_mount_material = G4Material::GetMaterial(this->s3_mount_material); - s3_mount_log = new G4LogicalVolume(s3_mount, s3_mount_material, "s3_mount_log", 0, 0, 0); - s3_mount_log->SetVisAttributes(vis_att); - -} // end::BuildS3Mount() -void ApparatusSpiceTargetChamber::BuildDetectorMount() { - - // ** Visualisation - G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(AL_COL)); - vis_att->SetVisibility(true); - - // ** Dimensions - // Box - G4double box_half_width = this->detector_mount_width/2.; - G4double box_half_length = this->detector_mount_length/2.; - G4double box_half_thickness = this->detector_mount_thickness/2.; - // Inner Radius - G4double lip_radius = this->detector_mount_lip_radius; - G4double lip_half_thickness = this->detector_mount_lip_thickness/2.; - G4double box_cut_radius = this->detector_mount_inner_radius; - // Annular Clamp - G4double clamp_half_thickness = this->annular_clamp_thickness; - G4double clamp_half_width = this->annular_clamp_width/2.; - G4double clamp_half_length = this->annular_clamp_length/2.; - - // ** Shapes - G4Box* mount_box = new G4Box("mount_box", box_half_width, box_half_length, box_half_thickness); - G4Tubs* inner_radius_cut = new G4Tubs("inner_radius_cut", 0, lip_radius, 2*box_half_thickness, 0, 360*deg); - G4Tubs* lip_cut = new G4Tubs("lip_cut", 0, box_cut_radius, box_half_thickness, 0, 360*deg); - G4Box* annular_clamp = new G4Box("annular_clamp", clamp_half_width, clamp_half_length, clamp_half_thickness); - - G4SubtractionSolid* detector_mount_pre = new G4SubtractionSolid("detector_mount_pre", mount_box, inner_radius_cut); - G4ThreeVector trans(0, 0, this->detector_mount_lip_thickness); - G4SubtractionSolid* detector_mount = new G4SubtractionSolid("detector_mount", detector_mount_pre, lip_cut, 0, trans); - - G4double plane_offset = (this->annular_clamp_plane_offset + clamp_half_length) / sqrt(2.); - G4double z_offset = box_half_thickness; - G4ThreeVector move(plane_offset, plane_offset, z_offset); - G4RotationMatrix* rotate = new G4RotationMatrix(45*deg, 0, 0); - G4SubtractionSolid* detector_mount2 = new G4SubtractionSolid("detector_mount2", detector_mount, annular_clamp, rotate, move); - move.setX(-plane_offset); - rotate->rotateZ(90*deg); - G4SubtractionSolid* detector_mount3 = new G4SubtractionSolid("detector_mount3", detector_mount2, annular_clamp, rotate, move); - move.setY(-plane_offset); - rotate->rotateZ(90*deg); - G4SubtractionSolid* detector_mount4 = new G4SubtractionSolid("detector_mount4", detector_mount3, annular_clamp, rotate, move); - move.setX(plane_offset); - rotate->rotateZ(90*deg); - G4SubtractionSolid* detector_mount5 = new G4SubtractionSolid("detector_mount5", detector_mount4, annular_clamp, rotate, move); - - // ** Logical - G4Material* detector_mount_material = G4Material::GetMaterial(this->detector_mount_material); - detector_mount_log = new G4LogicalVolume(detector_mount5, detector_mount_material, "detector_mount_log", 0, 0, 0); - detector_mount_log->SetVisAttributes(vis_att); - -} // end::BuildDetectorMount() - -void ApparatusSpiceTargetChamber::BuildAnnularClamps() { +void ApparatusSpiceTargetChamber::BuildColdFinger(){ // ** Visualisation - G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(PEEK_COL)); + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(CU_COL)); vis_att->SetVisibility(true); // ** Dimensions - G4double clamp_half_length = this->annular_clamp_length/2.; - G4double clamp_half_width = this->annular_clamp_width/2.; - G4double clamp_half_thickness = this->annular_clamp_thickness/2.; - // Distance - G4double beam_clamp_distance = this->annular_clamp_plane_offset + clamp_half_length; - + G4double coldfinger_half_thickness = this->coldfinger_thickness/2.; + G4double coldfinger_half_length = this->coldfinger_length/2.; + G4double coldfinger_half_width = this->coldfinger_width/2.; + G4double inner_hole_radius = this->coldfinger_hole_radius; + // ** Shapes - G4Box* annular_clamp = new G4Box("annular_clamp", clamp_half_width, clamp_half_length, clamp_half_thickness); - - G4ThreeVector move(2*beam_clamp_distance, 0, 0); - G4UnionSolid* double_clamps = new G4UnionSolid("double_clamps", annular_clamp, annular_clamp, 0, move); - - G4Box* annular_clamp2 = new G4Box("annular_clamp2", clamp_half_length, clamp_half_width, clamp_half_thickness); - G4ThreeVector trans(0, 2*beam_clamp_distance, 0); - G4UnionSolid* double_clamps2 = new G4UnionSolid("double_clamps2", annular_clamp2, annular_clamp2, 0, trans); - - G4ThreeVector trans2(beam_clamp_distance, -beam_clamp_distance, 0); - G4UnionSolid* four_clamps = new G4UnionSolid("four_clamps", double_clamps, double_clamps2, 0, trans2); + G4Box* plate_box = new G4Box("plate_box", coldfinger_half_length, coldfinger_half_width, coldfinger_half_thickness); + G4Tubs* inner_hole = new G4Tubs("inner_hole", 0, inner_hole_radius, 2*coldfinger_half_thickness, 0, 360*deg); + G4SubtractionSolid* plate_sub_hole = new G4SubtractionSolid("plate-hole", plate_box , inner_hole, 0, G4ThreeVector(0,0,0)); // ** Logical - G4Material* annular_clamp_material = G4Material::GetMaterial(this->annular_clamp_material); - annular_clamp_log = new G4LogicalVolume(four_clamps, annular_clamp_material, "annular_clamp_log", 0, 0, 0); - annular_clamp_log->SetVisAttributes(vis_att); + G4Material* cold_finger_material = G4Material::GetMaterial(this->cold_finger_material); + cold_finger_log = new G4LogicalVolume(plate_sub_hole, cold_finger_material, " cold_finger_log", 0, 0, 0); + cold_finger_log->SetVisAttributes(vis_att); -} // end::BuildAnnularClamps() +} // end:BuildColdFinger() + void ApparatusSpiceTargetChamber::PlaceTargetChamberFrontRing() { @@ -1610,17 +1480,17 @@ void ApparatusSpiceTargetChamber::PlacePhotonShield() + this->photon_shield_layer_three_thickness/2. + this->photon_shield_layer_two_thickness/2. + this->photon_shield_layer_one_thickness/2.); - // photon_shield_layer_two_phys = new G4PVPlacement(rotate,move,photon_shield_layer_two_log, - // "photon_shield_layer_two", expHallLog, - // false,0); + photon_shield_layer_two_phys = new G4PVPlacement(rotate,move,photon_shield_layer_two_log, + "photon_shield_layer_two", expHallLog, + false,0); move.set(0,0, this->photon_shield_back_face_pos + this->photon_shield_layer_three_thickness/2. + this->photon_shield_layer_two_thickness/2. + this->photon_shield_layer_one_thickness/2.); - // photon_shield_layer_three_phys = new G4PVPlacement(rotate,move,photon_shield_layer_three_log, - // "photon_shield_layer_three", expHallLog, - // false,0); + photon_shield_layer_three_phys = new G4PVPlacement(rotate,move,photon_shield_layer_three_log, + "photon_shield_layer_three", expHallLog, + false,0); }// end:PlacePhotonShield() @@ -1781,54 +1651,18 @@ void ApparatusSpiceTargetChamber::PlaceElectroBox() } // end:PlaceElectroBox() -void ApparatusSpiceTargetChamber::PlaceS3Mount() +void ApparatusSpiceTargetChamber::PlaceColdFinger() { - G4double z_offset = this->s3_mount_z_offset + this->s3_mount_thickness/2. + this->frontDomeOffset; - G4double angular_offset = this->s3_mount_angular_offset; + G4double z_offset = coldfinger_z_offset - coldfinger_thickness/2. ; G4ThreeVector move(0, 0, z_offset); - G4RotationMatrix* rotate = new G4RotationMatrix(angular_offset, 0, 0); - - s3_mount_phys = new G4PVPlacement(rotate, move, s3_mount_log, - "s3_mount", expHallLog, - false, 0); -} // end::PlaceS3Mount() - -void ApparatusSpiceTargetChamber::PlaceDetectorMount() -{ - - G4double detector_mount_gap = this->detector_mount_thickness - - this->detector_mount_lip_thickness - this->detector_thickness; - G4double z_offset = -this->detector_to_target_distance - - this->detector_mount_thickness/2. + detector_mount_gap + this->backDetAndAlBoxOffset; - - G4RotationMatrix* rotate = new G4RotationMatrix(this->detector_mount_angular_offset, 0, 0); - G4ThreeVector move(0, 0, z_offset); - detector_mount_phys = new G4PVPlacement(rotate, move, detector_mount_log, - "detector_mount", expHallLog, - false, 0); - -} // end::PlaceDetectorMount() - -void ApparatusSpiceTargetChamber::PlaceAnnularClamps() { - - G4double z_offset = -this->detector_to_target_distance + this->backDetAndAlBoxOffset; - G4double x_offset = (this->annular_clamp_plane_offset - + this->annular_clamp_length/2.) - * cos(this->detector_mount_angular_offset + 45*deg); - G4double y_offset = (this->annular_clamp_plane_offset - + this->annular_clamp_length/2.) - * sin(this->detector_mount_angular_offset + 45*deg); - - G4RotationMatrix* rotate = new G4RotationMatrix(this->detector_mount_angular_offset + 45*deg, 0, 0); - G4ThreeVector move(-x_offset, -y_offset, z_offset); - annular_clamp_phys = new G4PVPlacement(rotate, move, annular_clamp_log, - "annular_clamp", expHallLog, - false,0); + cold_finger_phys = new G4PVPlacement(0, move, cold_finger_log, + "cold_finger", expHallLog, + false, 0); -} // end::PlaceAnnularClamps() +} // end:PlaceColdFinger() //////////////////////////////////////////////////////////////////// diff --git a/src/DetectionSystemNew.cc b/src/DetectionSystemNew.cc new file mode 100644 index 0000000..ff5f22e --- /dev/null +++ b/src/DetectionSystemNew.cc @@ -0,0 +1,779 @@ +#include "DetectorConstruction.hh" + +#include "DetectorMessenger.hh" + +#include "G4Material.hh" + +#include "G4Tubs.hh" +#include "G4Box.hh" +#include "G4Sphere.hh" + +#include "G4Cons.hh" +#include "G4Trd.hh" +#include "G4LogicalVolume.hh" +#include "G4PVPlacement.hh" +#include "G4SubtractionSolid.hh" +#include "G4UnionSolid.hh" + +#include "G4GeometryManager.hh" +#include "G4PhysicalVolumeStore.hh" +#include "G4LogicalVolumeStore.hh" +#include "G4SolidStore.hh" +#include "G4AssemblyVolume.hh" +#include "G4VSolid.hh" +#include "G4Trap.hh" + +#include "G4VisAttributes.hh" +#include "G4Colour.hh" + +#include "DetectionSystemNew.hh" + +#include + + +//////////////////////////////////////////////// +// building and placing the different parts // +// of the chamber // +// including chamber itself, magnets and // +// mounts for magnets and detectors // +// not including detectors // +//////////////////////////////////////////////// + + +DetectionSystemNew::DetectionSystemNew() +{ + + this->detector_alignment = 1; // 1 for parallel to beam axis, 2 for perpindicular to beam axis + this->NUMBER_OF_MAGNETS = 8; + this->no_magnet_layer = 3; + + /////////////////////////////////////////////////////////////////////////// + // Assigning materials used to their relevant parts inside the detector // + /////////////////////////////////////////////////////////////////////////// + + this->chamber_material = "Delrin"; +// this->cold_finger_material = "Copper"; + this->magnet_material = "NdFeB"; + this->photon_shield_layer_one_material = "WTa" ; //MAY NEED TO CHANGE + this->photon_shield_layer_two_material = "Tin" ; //MAY NEED TO CHANGE + this->photon_shield_layer_three_material = "Aluminum" ; //MAY NEED TO CHANGE + // this->bias_plate_material = "Aluminum"; +// this->shield_cover_material = "Delrin"; +// this->magnet_cover_material = "Delrin"; +// this-> + + + + + //////////////////////////////////////////////////////// + // Defining the dimensions of the detector components // + //////////////////////////////////////////////////////// + + ////////////////////////////////// + // Chamber Hemisphere Component // + ////////////////////////////////// + + this->chamber_hemisphere_inner_radius = 100.00*mm; + this->chamber_hemisphere_outer_radius = 110.00*mm; + + ////////////////////////////////////////// + // Chamber Cylinder/Beam Pipe Component // + ////////////////////////////////////////// + + this->chamber_cylinder_inner_radius = 0*mm; + this->chamber_cylinder_outer_radius = 40*mm; + this->chamber_cylinder_height = 100*mm; + + ////////////////////////////// + // Chamber Cut Out Cylinder // + ////////////////////////////// + + this->chamber_cut_out_inner_radius = 0*mm; + this->chamber_cut_out_outer_radius = 35*mm; + this->chamber_cut_out_height = 300*mm; + + + ///////////// + // Magnets // + ///////////// + + this->plate_one_thickness = 3.4*mm; + this->plate_one_length = 82*mm; //SPICE = 75*mm; + this->plate_one_height = 60*mm; // z-component + + this->plate_two_thickness = 5*mm; + this->plate_two_length = 32*mm; //SPICE = 55*mm; + this->plate_two_height = 60*mm; + + this->top_cutting_box_angle = 65*deg; + this->bottom_cutting_box_angle = 65*deg; + this->top_cutting_box_lowering_height = 55*mm; + this->bottom_cutting_box_raising_height = 25*mm; + this->distance_from_target = 4*mm; // in z-direction + + + /////////////////// + // Magnet Covers // + /////////////////// + +// this->magnet_cover_thickness =2*mm; + + + ////////////////////// + // Photon Shielding // + ////////////////////// + + this->photon_shield_layer_one_thickness = (25)*mm; + this->photon_shield_layer_two_thickness = (4)*mm; + this->photon_shield_layer_three_thickness = (1)*mm; + this->photon_shield_target_distance = (7)*mm; + this->detector_target_distance = (75)*mm; + + if(this->detector_alignment == 1){ + this->detector_full_width = (62)*mm; + this->photon_shield_rotation = 22.5*deg ;} + else{ + this->detector_full_width = (86)*mm; + this->photon_shield_rotation = -22.5*deg ;} + + /////////////////////////////////// + // Photon Shielding Beam Cut Out // + /////////////////////////////////// + + this->photon_shield_beam_cut_outer_radius = (5)*mm; + this->photon_shield_beam_cut_inner_radius = (0)*mm; + this->photon_shield_beam_cut_height = (100)*mm; + + + //////////////////////// + // Cold Finger Piping // + //////////////////////// + +// this->cold_finger_inner_radius = (15)*mm; +// this->cold_finger_thickness = (3)*mm; +// this->cold_finger_height = (100)*mm; + + /////////////////////////////////////////////////////// + // Defining the positions of the detector components // + /////////////////////////////////////////////////////// + +// this->chamber_hemisphere_dist = (0.)*mm; +// this->chamber_cylinder_dist = (-75.)*mm; +// this->cold_finger_dist = (-75)*mm; + + +} //end NewDetectionSystem + +DetectionSystemNew::~DetectionSystemNew() +{ + + + /////////////////////////////////////////////////////////// + // Deleting the logical space of the detector components // + /////////////////////////////////////////////////////////// + + +// delete chamber_log; + delete final_chamber_log; +// delete cold_finger_log; + delete magnet_log; + delete photon_shield_layer_one_log; + delete photon_shield_layer_two_log; + delete photon_shield_layer_three_log; +// delete bias_plate_log; +// delete shield_cover_log; +// delete magnet_cover_log; + + //////////////////////////////////////////////////////////// + // Deleting the physical space of the detector components // + //////////////////////////////////////////////////////////// + +// delete chamber_phys; + delete final_chamber_phys; +// delete cold_finger_phys; + delete magnet_phys; + delete photon_shield_layer_one_phys; + delete photon_shield_layer_two_phys; + delete photon_shield_layer_three_phys; +// delete bias_plate_phys; +// delete shield_cover_phys; +// delete magnet_cover_phys; + +} //end ~NewDetectionSystem + + //////////////////////////////////////////////////////// + // Need to now build and place components of detector // + //////////////////////////////////////////////////////// + +void DetectionSystemNew::Build(G4LogicalVolume* exp_hall_log) +{ + this->expHallLog = exp_hall_log; + /////////////////////////////////////////////////// + // Building functions for components of detector // + /////////////////////////////////////////////////// + + BuildChamber(); +// BuildColdFinger(); + BuildMagnet(); + BuildPhotonShield(); +/// BuildBiasPlate(); +// BuildShieldCover(); +// BuildMagnetCover(); + + ////////////////////////////////////////////////// + // Placing functions for components of detector // + ////////////////////////////////////////////////// + + PlaceChamber(); + // PlaceColdFinger(); + PlacePhotonShield(); +// PlaceBiasPlate(); +// PlaceShieldCover(); + + for(G4int copyID=0; copyIDNUMBER_OF_MAGNETS; copyID++) + { + PlaceMagnet(copyID); +// PlaceMagnetCover(copyID); + } + +} //End Build + + +///////////////////////////////// +// Building components methods // +///////////////////////////////// + +//////////////////// +// Target Chamber // +//////////////////// + + + +void DetectionSystemNew::BuildChamber() +{ + //////////////////////////////// + // Visualisation of Component // + //////////////////////////////// + + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(CHAMBER_COL)); + vis_att->SetVisibility(true); + + ///////////////////////////// + // Dimensions of Component // + ///////////////////////////// + + + //Main hemisphere shell + G4double hemisphere_inner_radius = this->chamber_hemisphere_inner_radius; + G4double hemisphere_outer_radius = this->chamber_hemisphere_outer_radius; + + //Main Cylinder + G4double cylinder_inner_radius = this->chamber_cylinder_inner_radius; + G4double cylinder_outer_radius = this->chamber_cylinder_outer_radius; + G4double cylinder_height = this->chamber_cylinder_height/2; + + //Cut Out Cylinder + G4double cut_out_inner_radius = this->chamber_cut_out_inner_radius; + G4double cut_out_outer_radius = this->chamber_cut_out_outer_radius; + G4double cut_out_height = this->chamber_cut_out_height; + + //Primitive Volumes + G4Sphere* chamber_hemisphere = new G4Sphere("chamber_hemisphere", hemisphere_inner_radius, hemisphere_outer_radius, 0*pi, 2*pi, 90*deg,90*deg); + G4Tubs* chamber_cylinder = new G4Tubs("chamber_cylinder", cylinder_inner_radius, cylinder_outer_radius, cylinder_height, 0*deg, 360*deg); + G4Tubs* chamber_cut_out = new G4Tubs("chamber_cut_out", cut_out_inner_radius, cut_out_outer_radius, cut_out_height, 0*deg, 360*deg); + + //Putting Chamber components together + G4ThreeVector zTrans(0,0,-((cylinder_height*2)+51*mm)); + G4UnionSolid* chamber_solid = new G4UnionSolid("chamber_solid", chamber_hemisphere, chamber_cylinder,0,zTrans); + G4SubtractionSolid* final_chamber = new G4SubtractionSolid("final_chamber", chamber_solid, chamber_cut_out); + + //Logical Volume + G4Material* chamber_material = G4Material::GetMaterial(this->chamber_material); + final_chamber_log = new G4LogicalVolume(final_chamber,chamber_material,"final_chamber_log",0,0,0); + final_chamber_log->SetVisAttributes(vis_att); + + + +} //end of additions + + + + +/////////////////////// +// A Singular Magnet // +/////////////////////// + + + +void DetectionSystemNew::BuildMagnet() +{ + + //Visualisation + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(NDFEB_COL)); + vis_att->SetVisibility(true); + + + //Building the middle plate primary volume, middle plate made with cylinder with same radius as chamber and then cut to size + G4Tubs* plate_one_pre_cut = new G4Tubs("plate_one_pre_cut", 0., this->chamber_hemisphere_inner_radius, this->plate_one_thickness/2, 0., 90*deg); + G4Box* plate_one_cut_box = new G4Box("plate_one_cut_box", this->chamber_hemisphere_inner_radius, this->chamber_hemisphere_inner_radius, this->plate_one_thickness +1); + + //Setting the cut translations to make middle plate the correct size + G4double translate_cut_x_height_top = this->chamber_hemisphere_inner_radius -1 ; + G4double translate_cut_z_height_top = -(this->chamber_hemisphere_inner_radius) + this->distance_from_target; + G4ThreeVector translate_cut_height_top(translate_cut_x_height_top, translate_cut_z_height_top, 0); + + G4double translate_cut_x_height_bottom = this->chamber_hemisphere_inner_radius -1 ; + G4double translate_cut_z_height_bottom = this->chamber_hemisphere_inner_radius + this->distance_from_target + this->plate_one_height; + G4ThreeVector translate_cut_height_bottom(translate_cut_x_height_bottom, translate_cut_z_height_bottom, 0); + + G4double translate_cut_x_length = -this->chamber_hemisphere_inner_radius + (this->chamber_hemisphere_inner_radius - this->plate_one_length) ; + G4double translate_cut_z_length = this->chamber_hemisphere_inner_radius - 1 ; + G4ThreeVector translate_cut_length(translate_cut_x_length, translate_cut_z_length, 0); + + //Making the cuts from the magnets + G4SubtractionSolid* plate_one_temp1 = new G4SubtractionSolid("plate_one_temp1", plate_one_pre_cut, plate_one_cut_box, 0 , translate_cut_height_top); + G4SubtractionSolid* plate_one_temp2 = new G4SubtractionSolid("plate_one_temp2", plate_one_temp1, plate_one_cut_box, 0 , translate_cut_height_bottom); + G4SubtractionSolid* plate_one_temp3 = new G4SubtractionSolid("plate_one_temp3", plate_one_temp2, plate_one_cut_box, 0 , translate_cut_length); + + //Build angled cut box + G4Box* plate_one_angle_cut_box = new G4Box("plate_one_angle_cut_box", this->plate_one_length/2., this->plate_one_height/2., this->plate_one_thickness/2. + 1.); + + // ( Top Cut box is rotated and moved, exact movement calculated below ) + G4double angle_to_corner = atan(this->plate_one_length/this->plate_one_height); + G4double hypotenuse_to_corner = sqrt(pow(this->plate_one_height/2., 2) + pow(this->plate_one_length/2., 2)); + G4double angle_difference, translate_cut_x, translate_cut_z; + + if (angle_to_corner < this->top_cutting_box_angle) + { + angle_difference = this->top_cutting_box_angle - angle_to_corner; + translate_cut_x = -(hypotenuse_to_corner * sin(angle_difference)) +(this->chamber_hemisphere_inner_radius - this->plate_one_length) ; + translate_cut_z = (hypotenuse_to_corner * cos(angle_difference))- this->distance_from_target - this->top_cutting_box_lowering_height; + + } + else if (angle_to_corner > this->top_cutting_box_angle) + { + angle_difference = angle_to_corner - this->top_cutting_box_angle; + translate_cut_x = (this->chamber_hemisphere_inner_radius - this->plate_one_length) + (hypotenuse_to_corner * sin(angle_difference)); + translate_cut_z = (hypotenuse_to_corner * cos(angle_difference)) - this->distance_from_target - this->top_cutting_box_lowering_height; + + } + + G4ThreeVector translate_cut_box(translate_cut_x, -translate_cut_z, 0); + G4RotationMatrix* rotate_cut_box = new G4RotationMatrix; + rotate_cut_box->rotateZ(this->top_cutting_box_angle); + G4SubtractionSolid* plate_one_temp4 = new G4SubtractionSolid("plate_one_temp4", plate_one_temp3, plate_one_angle_cut_box, + rotate_cut_box, translate_cut_box); + + + + // ( Bottom Cut box is rotated and moved, exact movement calculated below ) + G4double angle_to_corner_bottom = atan(this->plate_one_length/this->plate_one_height); + G4double hypotenuse_to_corner_bottom = sqrt(pow(this->plate_one_height/2., 2) + pow(this->plate_one_length/2., 2)); + G4double angle_difference_bottom, translate_cut_x_bottom, translate_cut_z_bottom; + + if (angle_to_corner < this->bottom_cutting_box_angle) + { + angle_difference_bottom = this->bottom_cutting_box_angle - angle_to_corner_bottom; + translate_cut_x_bottom = -(hypotenuse_to_corner_bottom * sin(angle_difference_bottom)) +(this->chamber_hemisphere_inner_radius - this->plate_one_length) ; + translate_cut_z_bottom = -(hypotenuse_to_corner_bottom * cos(angle_difference_bottom))- this->distance_from_target - this->plate_one_height + this->bottom_cutting_box_raising_height; + + } + else if (angle_to_corner > this->bottom_cutting_box_angle) + { + angle_difference_bottom = angle_to_corner_bottom - this->bottom_cutting_box_angle; + translate_cut_x_bottom = (this->chamber_hemisphere_inner_radius - this->plate_one_length) + (hypotenuse_to_corner_bottom * sin(angle_difference_bottom)); + translate_cut_z_bottom = -(hypotenuse_to_corner_bottom * cos(angle_difference_bottom)) - this->distance_from_target - this->plate_one_height + this->bottom_cutting_box_raising_height; + + } + + G4ThreeVector translate_cut_box_bottom(translate_cut_x_bottom, -translate_cut_z_bottom, 0); + G4RotationMatrix* rotate_cut_box_bottom = new G4RotationMatrix; + rotate_cut_box_bottom->rotateZ(360*deg-this->bottom_cutting_box_angle); + G4SubtractionSolid* plate_one = new G4SubtractionSolid("plate_one", plate_one_temp4, plate_one_angle_cut_box, + rotate_cut_box_bottom, translate_cut_box_bottom); + + + + //Building the paddles primary volume, paddles made with cylinder with same radius as chamber and then cut to size + G4double paddle_thickness = (this->plate_two_thickness + this->plate_one_thickness/2.); + + G4Tubs* plate_two_pre_cut = new G4Tubs("plate_two_pre_cut", 0., this->chamber_hemisphere_inner_radius, paddle_thickness, 0., 90*deg); + G4Box* plate_two_cut_box = new G4Box("plate_two_cut_box", this->chamber_hemisphere_inner_radius, this->chamber_hemisphere_inner_radius, paddle_thickness +1); + + //Setting the cut translations to make middle plate the correct size + G4double translate_cut_x_height_top2 = this->chamber_hemisphere_inner_radius -1 ; + G4double translate_cut_z_height_top2 = -(this->chamber_hemisphere_inner_radius) + this->distance_from_target; + G4ThreeVector translate_cut_height_top2(translate_cut_x_height_top2, translate_cut_z_height_top2, 0); + + G4double translate_cut_x_height_bottom2 = this->chamber_hemisphere_inner_radius -1 ; + G4double translate_cut_z_height_bottom2 = this->chamber_hemisphere_inner_radius + this->distance_from_target + this->plate_two_height; + G4ThreeVector translate_cut_height_bottom2(translate_cut_x_height_bottom2, translate_cut_z_height_bottom2, 0); + + G4double translate_cut_x_length2 = -this->chamber_hemisphere_inner_radius + (this->chamber_hemisphere_inner_radius - this->plate_two_length) ; + G4double translate_cut_z_length2 = this->chamber_hemisphere_inner_radius - 1 ; + G4ThreeVector translate_cut_length2(translate_cut_x_length2, translate_cut_z_length2, 0); + + //Making the cuts from the magnets + G4SubtractionSolid* plate_two_temp1 = new G4SubtractionSolid("plate_two_temp1", plate_two_pre_cut, plate_two_cut_box, 0 , translate_cut_height_top2); + G4SubtractionSolid* plate_two_temp2 = new G4SubtractionSolid("plate_two_temp2", plate_two_temp1, plate_two_cut_box, 0 , translate_cut_height_bottom2); + G4SubtractionSolid* plate_two = new G4SubtractionSolid("plate_two", plate_two_temp2, plate_two_cut_box, 0 , translate_cut_length2); + + //Combine the plates to make the magnet + G4UnionSolid* fullmagnet = new G4UnionSolid("fullmagnet", plate_one, plate_two); + + + //Logical Volume + G4Material* magnet_material = G4Material::GetMaterial(this->magnet_material); + magnet_log = new G4LogicalVolume(fullmagnet, magnet_material, "magnet_log", 0, 0, 0); + magnet_log->SetVisAttributes(vis_att); + + +} // end:BuildCollectorMagnet() + + + +/////////////////// +// Photon Shield // +/////////////////// + +void DetectionSystemNew::BuildPhotonShield() +{ + + //Visualisation + G4VisAttributes* vis_att1 = new G4VisAttributes(G4Colour(LAYER1_COL)); + vis_att1->SetVisibility(true); + G4VisAttributes* vis_att2 = new G4VisAttributes(G4Colour(LAYER2_COL)); + vis_att2->SetVisibility(true); + G4VisAttributes* vis_att3 = new G4VisAttributes(G4Colour(LAYER3_COL)); + vis_att3->SetVisibility(true); + + //Bring in some dimensions pre set at start of code + G4double detector_full_width = this->detector_full_width/2; + G4double target_distance = this->photon_shield_target_distance; + G4double detector_target_distance = this->detector_target_distance; + + G4double layer_one_thickness = this->photon_shield_layer_one_thickness; + G4double layer_two_thickness = this->photon_shield_layer_two_thickness; + G4double layer_three_thickness = this->photon_shield_layer_three_thickness; + + //Cut Out Dimensions + G4double cut_inner_radius = this->photon_shield_beam_cut_inner_radius; + G4double cut_outer_radius = this->photon_shield_beam_cut_outer_radius; + G4double cut_height = this->photon_shield_beam_cut_height; + + + //Defining the angle made with the z plane from target edge to detector edge + G4double phi = atan((detector_full_width - cut_outer_radius/2)/detector_target_distance); + + //Defining the First layer x-y plane dimensions using Phi (Uses trig to add on extra length due to small right angle triangle in trapezium) + G4double layer_one_face_xy_dimension = (cut_outer_radius/2 + (tan(phi)*target_distance)); + G4double layer_one_base_xy_dimension = (layer_one_face_xy_dimension + (tan(phi)*layer_one_thickness)); + + //Defining the Second layer x-y plane dimensions using Phi (Uses trig to add on extra length due to small right angle triangle in trapezium) + G4double layer_two_face_xy_dimension = layer_one_base_xy_dimension; + G4double layer_two_base_xy_dimension = (layer_two_face_xy_dimension + (tan(phi)*layer_two_thickness)); + + //Defining the Third layer x-y plane dimensions using Phi (Uses trig to add on extra length due to small right angle triangle in trapezium) + G4double layer_three_face_xy_dimension = layer_two_base_xy_dimension; + G4double layer_three_base_xy_dimension = (layer_three_face_xy_dimension + (tan(phi)*layer_three_thickness)); + + + //Primitive Volumes + G4Trd* photon_shield_layer_one_pre = new G4Trd("photon_shield_layer_one_pre", layer_one_base_xy_dimension, + layer_one_face_xy_dimension, layer_one_base_xy_dimension, layer_one_face_xy_dimension, (layer_one_thickness/2)); + G4Trd* photon_shield_layer_two_pre = new G4Trd("photon_shield_layer_two_pre", layer_two_base_xy_dimension, + layer_two_face_xy_dimension, layer_two_base_xy_dimension, layer_two_face_xy_dimension, (layer_two_thickness/2)); + G4Trd* photon_shield_layer_three_pre = new G4Trd("photon_shield_layer_three_pre", layer_three_base_xy_dimension, + layer_three_face_xy_dimension, layer_three_base_xy_dimension, layer_three_face_xy_dimension, (layer_three_thickness/2)); + + G4Tubs* photon_shield_cut_out = new G4Tubs("photon_shield_cut_out", cut_inner_radius, cut_outer_radius, cut_height, 0*deg, 360*deg); + + //Putting Photon Shield components together + G4SubtractionSolid* photon_shield_layer_one = new G4SubtractionSolid("photon_shield_layer_one", photon_shield_layer_one_pre, photon_shield_cut_out); + G4SubtractionSolid* photon_shield_layer_two = new G4SubtractionSolid("photon_shield_layer_two", photon_shield_layer_two_pre, photon_shield_cut_out); + G4SubtractionSolid* photon_shield_layer_three = new G4SubtractionSolid("photon_shield_layer_three", photon_shield_layer_three_pre, photon_shield_cut_out); + + //Logical Volumes + G4Material* photon_shield_layer_one_material = G4Material::GetMaterial(this->photon_shield_layer_one_material); + photon_shield_layer_one_log = new G4LogicalVolume(photon_shield_layer_one,photon_shield_layer_one_material,"photon_shield_layer_one_log",0,0,0); + photon_shield_layer_one_log->SetVisAttributes(vis_att1); + + G4Material* photon_shield_layer_two_material = G4Material::GetMaterial(this->photon_shield_layer_two_material); + photon_shield_layer_two_log = new G4LogicalVolume(photon_shield_layer_two,photon_shield_layer_two_material,"photon_shield_layer_two_log",0,0,0); + photon_shield_layer_two_log->SetVisAttributes(vis_att2); + + G4Material* photon_shield_layer_three_material = G4Material::GetMaterial(this->photon_shield_layer_three_material); + photon_shield_layer_three_log = new G4LogicalVolume(photon_shield_layer_three,photon_shield_layer_three_material,"photon_shield_layer_three_log",0,0,0); + photon_shield_layer_three_log->SetVisAttributes(vis_att3); + +} + + + +///////////////// +// Cold Finger // +///////////////// + +//void DetectionSystemNew::BuildColdFinger() +//{ + + //////////////////////////////// + // Visualisation of Component // + //////////////////////////////// + +// G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(COLD_FINGER_COL)); +// vis_att->SetVisibility(true); + + + ///////////////////////////// + // Dimensions of Component // + ///////////////////////////// + + //Main Cylinder +// G4double inner_radius = this->cold_finger_inner_radius; +// G4double outer_radius = this->(cold_finger_inner_radius + this->cold_finger_thickness); +// G4double height = this->cold_finger_height/2; + + //Primitive Volumes +// G4Tubs* cold_finger = new G4Tubs("cold_finger", inner_radius, outer_radius, height, 0*deg, 360*deg); + + //Logical Volume +// G4Material* cold_finger_material = G4Material::GetMaterial(this->cold_finger_material); +// cold_finger_log = new G4LogicalVolume(cold_finger,cold_finger_material,"cold_finger_log",0,0,0); +// cold_finger_log->SetVisAttributes(vis_att); + +//}//end of building cold finger + + +//////////////////////////////// +// Placing components methods // +//////////////////////////////// + +///////////// +// Chamber // +///////////// + + +void DetectionSystemNew::PlaceChamber() +{ + + G4ThreeVector move(0,0,0); + final_chamber_phys = new G4PVPlacement(0,move,final_chamber_log, "final_chamber",expHallLog, false,0); + +}// end:PlaceChamber + +/////////////////// +// Photon Shield // +/////////////////// + +void DetectionSystemNew::PlacePhotonShield() +{ + + G4double z_position_layer_one = -(this->photon_shield_target_distance+this->photon_shield_layer_one_thickness/2) ; + G4ThreeVector move1(0,0,z_position_layer_one); + G4RotationMatrix* rotate1= new G4RotationMatrix; + if(this->detector_alignment == 1){ + rotate1->rotateZ(this->photon_shield_rotation);} + else{ + rotate1->rotateZ(this->photon_shield_rotation);} + photon_shield_layer_one_phys = new G4PVPlacement(rotate1,move1,photon_shield_layer_one_log, "photon_shield_layer_one",expHallLog, false,0); + + G4double z_position_layer_two = -(this->photon_shield_target_distance+this->photon_shield_layer_one_thickness+(this->photon_shield_layer_two_thickness/2)); + G4ThreeVector move2(0,0,z_position_layer_two); + G4RotationMatrix* rotate2= new G4RotationMatrix; + if(this->detector_alignment == 1){ + rotate2->rotateZ(this->photon_shield_rotation);} + else{ + rotate2->rotateZ(this->photon_shield_rotation);} + photon_shield_layer_two_phys = new G4PVPlacement(rotate2,move2,photon_shield_layer_two_log, "photon_shield_layer_two",expHallLog, false,0); + + G4double z_position_layer_three = -(this->photon_shield_target_distance+this->photon_shield_layer_one_thickness+this->photon_shield_layer_two_thickness+(this->photon_shield_layer_three_thickness/2)); + G4ThreeVector move3(0,0,z_position_layer_three); + G4RotationMatrix* rotate3= new G4RotationMatrix; + if(this->detector_alignment == 1){ + rotate3->rotateZ(this->photon_shield_rotation);} + else{ + rotate3->rotateZ(this->photon_shield_rotation);} + photon_shield_layer_three_phys = new G4PVPlacement(rotate3,move3,photon_shield_layer_three_log, "photon_shield_layer_three",expHallLog, false,0); + + + +}// end:PlacePhotonShield + +///////////// +// Magnets // +///////////// + +void DetectionSystemNew::PlaceMagnet(G4int copyID) +{ + + // ** Position Co-ordinates + G4double magnetPosX = 0*mm; + G4double magnetPosZ = 0*mm; + + G4RotationMatrix* rotation; + rotation = RotateMagnets(copyID); + G4double radial_position = magnetPosX; + + G4ThreeVector move = TranslateMagnets(copyID, radial_position, magnetPosZ); + + // ** Physical Volume + magnet_phys = new G4PVPlacement(rotation, move, magnet_log, + "collector_magnet",expHallLog,false,0); + +} // end:PlaceMagnet() + + +///////////// +// Magnets // +///////////// +/* +void DetectionSystemNew::PlaceMagnetCover(G4int copyID) +{ + + // ** Position Co-ordinates + G4double magnetPosX = this->plate_one_edge_x + this->plate_one_length/2.; - this->magnet_cover_thickness; + G4double magnetPosZ =(this->distance_from_target + this->plate_one_height/2.); + + G4RotationMatrix* rotation; + rotation = RotateMagnets(copyID); + G4double radial_position = magnetPosX; + + G4ThreeVector move = TranslateMagnets(copyID, radial_position, magnetPosZ); + + // ** Physical Volume + magnet_cover_phys = new G4PVPlacement(rotation, move, magnet_cover_log, + "collector_magnet_cover",expHallLog, false,0); + +} // end:PlaceMagnetCover() +*/ + + + + +//////////////////////////////////////////////////////////////// +// Translations and Rotations of Magnets (same for all parts) // +//////////////////////////////////////////////////////////////// + +G4RotationMatrix* DetectionSystemNew::RotateMagnets(G4int copyID) +{ + + G4RotationMatrix* rotate = new G4RotationMatrix; + if(this->NUMBER_OF_MAGNETS == 8) +{ + rotate->rotateZ(-(copyID)*45.*deg); + rotate->rotateX(90.*deg); + /*if(copyID==0){ + rotate->rotateZ(350.*deg);} + else if(copyID==1){ + rotate->rotateZ(280*deg);} + else if(copyID==2){ + rotate->rotateZ(260.*deg);} + else if(copyID==3){ + rotate->rotateZ(190.*deg);} + else if(copyID==4){ + rotate->rotateZ(170.*deg);} + else if(copyID==5){ + rotate->rotateZ(100*deg);} + else if(copyID==6){ + rotate->rotateZ(80.*deg);} + else if(copyID==7){ + rotate->rotateZ(10.*deg);} +if(copyID==0){ + rotate->rotateZ(345.*deg);} + else if(copyID==1){ + rotate->rotateZ(270*deg);} + else if(copyID==2){ + rotate->rotateZ(195.*deg);} + else if(copyID==3){ + rotate->rotateZ(165.*deg);} + else if(copyID==4){ + rotate->rotateZ(90.*deg);} + else if(copyID==5){ + rotate->rotateZ(15*deg);}*/ +} + else { + rotate->rotateZ(-(copyID+0.5)*90.*deg); + rotate->rotateX(90.*deg); + /*if(copyID==0){ + rotate->rotateZ(330.*deg);} + else if(copyID==1){ + rotate->rotateZ(210*deg);} + else if(copyID==2){ + rotate->rotateZ(-210.*deg);} + else if(copyID==3){ + rotate->rotateZ(30.*deg);}*/ + +} + + return rotate; +} + +G4ThreeVector DetectionSystemNew::TranslateMagnets(G4int copyID, G4double radial_position, G4double z_position) +{ + G4double x_position(0); + G4double y_position(0); + if(this->NUMBER_OF_MAGNETS ==8){ + x_position = radial_position*cos((copyID)*45.*deg); + y_position = radial_position*sin((copyID)*45.*deg); + /* if(copyID==0){ + x_position = radial_position*cos(10.*deg); + y_position = radial_position*sin(10.*deg);} + else if(copyID==1){ + x_position = radial_position*cos(80.*deg); + y_position = radial_position*sin(80.*deg);} + else if(copyID==2){ + x_position = radial_position*cos(100.*deg); + y_position = radial_position*sin(100.*deg);} + else if(copyID==3){ + x_position = radial_position*cos(170.*deg); + y_position = radial_position*sin(170.*deg);} + else if(copyID==4){ + x_position = radial_position*cos(190.*deg); + y_position = radial_position*sin(190.*deg);} + else if(copyID==5){ + x_position = radial_position*cos(260.*deg); + y_position = radial_position*sin(260.*deg);} + else if(copyID==6){ + x_position = radial_position*cos(280.*deg); + y_position = radial_position*sin(280.*deg);} + else if(copyID==7){ + x_position = radial_position*cos(-10.*deg); + y_position = radial_position*sin(-10.*deg);} +if(copyID==0){ + x_position = radial_position*cos(15.*deg); + y_position = radial_position*sin(15.*deg);} + else if(copyID==1){ + x_position = radial_position*cos(90.*deg); + y_position = radial_position*sin(90.*deg);} + else if(copyID==2){ + x_position = radial_position*cos(165.*deg); + y_position = radial_position*sin(165.*deg);} + else if(copyID==3){ + x_position = radial_position*cos(195.*deg); + y_position = radial_position*sin(195.*deg);} + else if(copyID==4){ + x_position = radial_position*cos(270.*deg); + y_position = radial_position*sin(270.*deg);} + else if(copyID==5){ + x_position = radial_position*cos(345.*deg); + y_position = radial_position*sin(345.*deg);}*/ + } + else{ + x_position = radial_position*cos((copyID+0.5)*90.*deg); + y_position = radial_position*sin((copyID+0.5)*90.*deg); + /*if(copyID==0){ + x_position = radial_position*cos(30.*deg); + y_position = radial_position*sin(30.*deg);} + else if(copyID==1){ + x_position = radial_position*cos(150.*deg); + y_position = radial_position*sin(150.*deg);} + else if(copyID==2){ + x_position = radial_position*cos(210.*deg); + y_position = radial_position*sin(210.*deg);} + else if(copyID==3){ + x_position = radial_position*cos(-30.*deg); + y_position = radial_position*sin(-30.*deg);}*/ + } + return G4ThreeVector(x_position, y_position, z_position); +} diff --git a/src/DetectionSystemS3.cc b/src/DetectionSystemS3.cc index 6094dd8..4171685 100755 --- a/src/DetectionSystemS3.cc +++ b/src/DetectionSystemS3.cc @@ -38,7 +38,7 @@ DetectionSystemS3::DetectionSystemS3() : //-----------------------------// this->S3DetCrystalOuterDiameter = 70.*mm; this->S3DetCrystalInnerDiameter = 22.*mm; - this->S3DetCrystalThickness = .15*mm; + this->S3DetCrystalThickness = 1.0*mm; // .15*mm this->S3DetRadialSegments = 24.; this->S3DetPhiSegments = 32.; @@ -46,15 +46,27 @@ DetectionSystemS3::DetectionSystemS3() : // parameters for guard ring // //-----------------------------// this->S3DetGuardRingOuterDiameter = 76*mm; - this->S3DetGuardRingInnerDiameter = 20*mm; + this->S3DetGuardRingInnerDiameter = 20*mm; + + // ------------------------- + // Dimensions of Si-CD Mount + // ------------------------- + this->s3_mount_material = "Peek"; //? + this->s3_mount_length = 120*mm; + this->s3_mount_thickness = 2.3*mm; + this->s3_active_radius = 35*mm; + this->s3_mount_chamfer = 28.284*mm; + this->s3_mount_centre_to_chamfer = 70.711*mm; + } DetectionSystemS3::~DetectionSystemS3() { delete [] siDetS3Ring_log; delete S3InnerGuardRing_log; - delete S3OuterGuardRing_log; - + delete S3OuterGuardRing_log; + delete s3_mount_log; + delete s3_mount_phys; } //---------------------------------------------------------// @@ -71,23 +83,25 @@ G4int DetectionSystemS3::Build() BuildSiliconWafer(ringID); // Build Silicon Ring } // end for(int ringID) - // Build Guard Rings + // Build Guard Rings and Mounting PCB BuildInnerGuardRing(); BuildOuterGuardRing(); - + BuildS3Mount(); + return 1; } //---------------------------------------------------------// // "place" function called in DetectorMessenger // // if detector is added // -//---------------------------------------------------------// -G4int DetectionSystemS3::PlaceDetector(G4LogicalVolume* exp_hall_log, G4ThreeVector move, G4int ringNumber, G4int Seg, G4int detectorNumber) -{ - G4RotationMatrix* rotate = new G4RotationMatrix; +//---------------------------------------------------------//logicWorld, pos, rotate , ring, Seg, detID +G4int DetectionSystemS3::PlaceDetector(G4LogicalVolume* exp_hall_log, G4ThreeVector move, G4double angle_offset, G4int ringNumber, G4int Seg, G4int detectorNumber) +{ G4int NumberSeg = (G4int)this->S3DetPhiSegments; - G4double angle = (360./NumberSeg)*(Seg-0.5)*deg; - rotate->rotateZ(angle); + G4RotationMatrix* rotate = new G4RotationMatrix; + G4double angle = ( (360./NumberSeg)*(Seg-0.5) + angle_offset )*deg ; + rotate->rotateZ(angle); + assemblyS3Ring[ringNumber]->MakeImprint(exp_hall_log, move, rotate, detectorNumber); return 1; @@ -187,7 +201,7 @@ G4int DetectionSystemS3::BuildOuterGuardRing() if( !material ) { G4cout << " ----> Material " << this->wafer_material << " not found, cannot build the outer guard ring of the S3 detector! " << G4endl; return 0; - } + } // Set visualization attributes G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(0.0,1.0,1.0)); @@ -218,6 +232,22 @@ G4int DetectionSystemS3::BuildOuterGuardRing() } + void DetectionSystemS3::PlaceS3Mount(G4LogicalVolume* exp_hall_log, G4ThreeVector move, G4double angular_offset) +{ + + G4RotationMatrix* rotate = new G4RotationMatrix(angular_offset, 0, 0); + G4double z_offset = this->S3DetCrystalThickness/2.; //this->s3_mount_thickness/2. + this->S3DetCrystalThickness/2. /*+ this->frontDomeOffset (this offset was used historically for visu checking as far as I understood)*/; + G4ThreeVector move_offset(0, 0, -z_offset); + move = move + move_offset ; + + s3_mount_phys = new G4PVPlacement(rotate, move, s3_mount_log, + "s3_mount", exp_hall_log, false, 0); + +} // end::PlaceS3Mount() + + + + /////////////////////////////////////////////////////// // Build one segment of S3 // the geometry depends on the distance from the center @@ -237,4 +267,43 @@ G4Tubs* DetectionSystemS3::BuildCrystal(G4int RingID) return crystal_block; }//end ::BuildCrystal + +void DetectionSystemS3::BuildS3Mount() { + + // ** Visualisation + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(PEEK_COL)); + vis_att->SetVisibility(true); + + // ** Dimensions + G4double half_length = this->s3_mount_length/2.; + G4double half_thickness = this->s3_mount_thickness/2.; + // Chamfer + G4double chamfer_half_length = this->s3_mount_chamfer/2.; + // Active Area Cut + G4double active_radius = this->s3_active_radius; + + // ** Shapes + G4Box* s3mount_box = new G4Box("s3_mount_box", half_length, half_length, half_thickness); + G4Box* chamfer_cut = new G4Box("chamber_cut", chamfer_half_length, chamfer_half_length, chamfer_half_length); + G4Tubs* active_cut = new G4Tubs("active_cut", 0, active_radius, chamfer_half_length, 0, 360*deg); + + G4double plane_offset = (this->s3_mount_centre_to_chamfer + chamfer_half_length) / sqrt(2.); + G4ThreeVector trans(plane_offset, plane_offset, 0); + G4RotationMatrix* rotate = new G4RotationMatrix(45*deg, 0, 0); + G4SubtractionSolid* s3_mount0 = new G4SubtractionSolid("s3_mount0", s3mount_box, chamfer_cut, rotate, trans); + trans.setX(-plane_offset); + G4SubtractionSolid* s3_mount1 = new G4SubtractionSolid("s3_mount1", s3_mount0, chamfer_cut, rotate, trans); + trans.setY(-plane_offset); + G4SubtractionSolid* s3_mount2 = new G4SubtractionSolid("s3_mount2", s3_mount1, chamfer_cut, rotate, trans); + trans.setX(plane_offset); + G4SubtractionSolid* s3_mount3 = new G4SubtractionSolid("s3_mount3", s3_mount2, chamfer_cut, rotate, trans); + G4SubtractionSolid* s3_mount = new G4SubtractionSolid("s3_mount", s3_mount3, active_cut); + + // ** Logical + G4Material* s3_mount_material = G4Material::GetMaterial(this->s3_mount_material); + s3_mount_log = new G4LogicalVolume(s3_mount, s3_mount_material, "s3_mount_log", 0, 0, 0); + s3_mount_log->SetVisAttributes(vis_att); + +} // end::BuildS3Mount() + diff --git a/src/DetectionSystemSpice.cc b/src/DetectionSystemSpice.cc index 6e511e4..33686bb 100755 --- a/src/DetectionSystemSpice.cc +++ b/src/DetectionSystemSpice.cc @@ -9,6 +9,7 @@ #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4SubtractionSolid.hh" +#include "G4UnionSolid.hh" #include "G4GeometryManager.hh" #include "G4PhysicalVolumeStore.hh" @@ -29,9 +30,34 @@ DetectionSystemSpice::DetectionSystemSpice() : ///////////////////////////////////////////////////////////////////// // SPICE Physical Properties ///////////////////////////////////////////////////////////////////// - - this->wafer_material = "Silicon"; - + //-----------------------------// + // Materials // + //-----------------------------// + this->wafer_material = "Silicon"; + this->detector_mount_material = "Aluminum"; // CHECK ? + this->annular_clamp_material = "Peek"; // CHECK ? + + // ---------------------------- + // Dimensions of Detector Mount + // ---------------------------- + this->detector_mount_length = 148*mm; + this->detector_mount_width = 130*mm; + this->detector_mount_thickness = 8*mm; + this->detector_mount_inner_radius = 52.8*mm; + this->detector_mount_lip_radius = 48.3*mm; + this->detector_mount_lip_thickness = 1*mm; + this->detector_mount_angular_offset = 0*deg; + this->detector_to_target_distance = 115*mm; + this->detector_thickness = 6*mm; + + // --------------------------- + // Dimensions of Annular Clamp + // --------------------------- + this->annular_clamp_thickness = 2*mm; + this->annular_clamp_length = 19.2*mm; + this->annular_clamp_width = 20*mm; + this->annular_clamp_plane_offset = 47.8*mm; //from beam line to first edge + //-----------------------------// // parameters for the annular // // planar detector crystal // @@ -51,12 +77,17 @@ DetectionSystemSpice::DetectionSystemSpice() : } DetectionSystemSpice::~DetectionSystemSpice() -{ // LogicalVolumes in ConstructSPICEDetectionSystem - delete [] siDetSpiceRing_log; - - delete siInnerGuardRing_log; - delete siOuterGuardRing_log; - +{ + // LogicalVolumes in ConstructSPICEDetectionSystem + delete detector_mount_log; + delete annular_clamp_log; + delete [] siDetSpiceRing_log; + delete siInnerGuardRing_log; + delete siOuterGuardRing_log; + + //Physical volumes + delete detector_mount_phys; + delete annular_clamp_phys; } //---------------------------------------------------------// @@ -80,7 +111,9 @@ G4int DetectionSystemSpice::Build() BuildInnerGuardRing(); BuildOuterGuardRing(); - + BuildDetectorMount(); + BuildAnnularClamps(); + return 1; } // end Build @@ -99,7 +132,7 @@ G4int DetectionSystemSpice::PlaceDetector(G4LogicalVolume* exp_hall_log, G4Three rotate->rotateZ(-210*deg-angle); // the axis are inverted, this operation will correct for it [MHD : 03 April 2014] assemblySiRing[ringNumber]->MakeImprint(exp_hall_log, move, rotate, SegmentNumber); - + return 1; } @@ -111,6 +144,41 @@ G4int DetectionSystemSpice::PlaceGuardRing(G4LogicalVolume* exp_hall_log, G4Thre return 1; } + +void DetectionSystemSpice::PlaceDetectorMount(G4LogicalVolume* exp_hall_log, G4ThreeVector move) +{ + + G4double detector_mount_gap = this->detector_mount_thickness + - this->detector_mount_lip_thickness - this->detector_thickness; + + G4double z_offset = - this->detector_mount_thickness/2. + detector_mount_gap ; + G4ThreeVector offset(0, 0, z_offset); + + move = move + offset ; + G4RotationMatrix* rotate = new G4RotationMatrix(this->detector_mount_angular_offset, 0, 0); + detector_mount_phys = new G4PVPlacement(rotate, move, detector_mount_log, + "detector_mount", exp_hall_log, + false, 0); + +} // end::PlaceDetectorMount() + +void DetectionSystemSpice::PlaceAnnularClamps(G4LogicalVolume* exp_hall_log, G4ThreeVector move) { + + G4double z_offset = this->annular_clamp_thickness/2. ; + G4double x_offset = (this->annular_clamp_plane_offset + + this->annular_clamp_length/2.) + * cos(this->detector_mount_angular_offset + 45*deg); + G4double y_offset = (this->annular_clamp_plane_offset + + this->annular_clamp_length/2.) + * sin(this->detector_mount_angular_offset + 45*deg); + G4ThreeVector offset(-x_offset, -y_offset, z_offset); + + move = move + offset ; + G4RotationMatrix* rotate = new G4RotationMatrix(this->detector_mount_angular_offset + 45*deg, 0, 0); + annular_clamp_phys = new G4PVPlacement(rotate, move, annular_clamp_log,"annular_clamp", exp_hall_log, + false,0); + +} // end::PlaceAnnularClamps() //---------------------------------------------------------// // build functions for different parts // @@ -155,7 +223,6 @@ G4int DetectionSystemSpice::BuildSiliconWafer(G4int RingID) // RingID = { 0, 9 this->assemblySiRing[RingID]->AddPlacedVolume(siDetSpiceRing_log[RingID], move, rotate); - return 1; } @@ -229,7 +296,96 @@ G4int DetectionSystemSpice::BuildOuterGuardRing() this->assembly->AddPlacedVolume(siOuterGuardRing_log, move, rotate); return 1; -} +} + + +void DetectionSystemSpice::BuildDetectorMount() { + + // ** Visualisation + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(AL_COL)); + vis_att->SetVisibility(true); + + // ** Dimensions + // Box + G4double box_half_width = this->detector_mount_width/2.; + G4double box_half_length = this->detector_mount_length/2.; + G4double box_half_thickness = this->detector_mount_thickness/2.; + // Inner Radius + G4double lip_radius = this->detector_mount_lip_radius; + //G4double lip_half_thickness = this->detector_mount_lip_thickness/2.; + G4double box_cut_radius = this->detector_mount_inner_radius; + // Annular Clamp + G4double clamp_half_thickness = this->annular_clamp_thickness; + G4double clamp_half_width = this->annular_clamp_width/2.; + G4double clamp_half_length = this->annular_clamp_length/2.; + + // ** Shapes + G4Box* mount_box = new G4Box("mount_box", box_half_width, box_half_length, box_half_thickness); + G4Tubs* inner_radius_cut = new G4Tubs("inner_radius_cut", 0, lip_radius, 2*box_half_thickness, 0, 360*deg); + G4Tubs* lip_cut = new G4Tubs("lip_cut", 0, box_cut_radius, box_half_thickness, 0, 360*deg); + G4Box* annular_clamp = new G4Box("annular_clamp", clamp_half_width, clamp_half_length, clamp_half_thickness); + + G4SubtractionSolid* detector_mount_pre = new G4SubtractionSolid("detector_mount_pre", mount_box, inner_radius_cut); + G4ThreeVector trans(0, 0, this->detector_mount_lip_thickness); + G4SubtractionSolid* detector_mount = new G4SubtractionSolid("detector_mount", detector_mount_pre, lip_cut, 0, trans); + + G4double plane_offset = (this->annular_clamp_plane_offset + clamp_half_length) / sqrt(2.); + G4double z_offset = box_half_thickness; + G4ThreeVector move(plane_offset, plane_offset, z_offset); + G4RotationMatrix* rotate = new G4RotationMatrix(45*deg, 0, 0); + G4SubtractionSolid* detector_mount2 = new G4SubtractionSolid("detector_mount2", detector_mount, annular_clamp, rotate, move); + move.setX(-plane_offset); + rotate->rotateZ(90*deg); + G4SubtractionSolid* detector_mount3 = new G4SubtractionSolid("detector_mount3", detector_mount2, annular_clamp, rotate, move); + move.setY(-plane_offset); + rotate->rotateZ(90*deg); + G4SubtractionSolid* detector_mount4 = new G4SubtractionSolid("detector_mount4", detector_mount3, annular_clamp, rotate, move); + move.setX(plane_offset); + rotate->rotateZ(90*deg); + G4SubtractionSolid* detector_mount5 = new G4SubtractionSolid("detector_mount5", detector_mount4, annular_clamp, rotate, move); + + // ** Logical + G4Material* detector_mount_material = G4Material::GetMaterial(this->detector_mount_material); + detector_mount_log = new G4LogicalVolume(detector_mount5, detector_mount_material, "detector_mount_log", 0, 0, 0); + detector_mount_log->SetVisAttributes(vis_att); + + //this->assembly->AddPlacedVolume(detector_mount_log, move, rotate); CHECK!! + +} // end::BuildDetectorMount() + +void DetectionSystemSpice::BuildAnnularClamps() { + + // ** Visualisation + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(PEEK_COL)); + vis_att->SetVisibility(true); + + // ** Dimensions + G4double clamp_half_length = this->annular_clamp_length/2.; + G4double clamp_half_width = this->annular_clamp_width/2.; + G4double clamp_half_thickness = this->annular_clamp_thickness/2.; + // Distance + G4double beam_clamp_distance = this->annular_clamp_plane_offset + clamp_half_length; + + // ** Shapes + G4Box* annular_clamp = new G4Box("annular_clamp", clamp_half_width, clamp_half_length, clamp_half_thickness); + + G4ThreeVector move(2*beam_clamp_distance, 0, 0); + G4UnionSolid* double_clamps = new G4UnionSolid("double_clamps", annular_clamp, annular_clamp, 0, move); + + G4Box* annular_clamp2 = new G4Box("annular_clamp2", clamp_half_length, clamp_half_width, clamp_half_thickness); + G4ThreeVector trans(0, 2*beam_clamp_distance, 0); + G4UnionSolid* double_clamps2 = new G4UnionSolid("double_clamps2", annular_clamp2, annular_clamp2, 0, trans); + + G4ThreeVector trans2(beam_clamp_distance, -beam_clamp_distance, 0); + G4UnionSolid* four_clamps = new G4UnionSolid("four_clamps", double_clamps, double_clamps2, 0, trans2); + + // ** Logical + G4Material* annular_clamp_material = G4Material::GetMaterial(this->annular_clamp_material); + annular_clamp_log = new G4LogicalVolume(four_clamps, annular_clamp_material, "annular_clamp_log", 0, 0, 0); + annular_clamp_log->SetVisAttributes(vis_att); + +} // end::BuildAnnularClamps() + /////////////////////////////////////////////////////// // Build one segment of Spice, diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc index f0b2792..9b8de43 100755 --- a/src/DetectorConstruction.cc +++ b/src/DetectorConstruction.cc @@ -70,6 +70,8 @@ #include "DetectionSystemSpice.hh" #include "DetectionSystemS3.hh" #include "DetectionSystemPaces.hh" +#include "DetectionSystemNew.hh" +#include "NewSquareDetector.hh" #include "DetectionSystemSodiumIodide.hh" #include "DetectionSystemLanthanumBromide.hh" #include "ApparatusGenericTarget.hh" @@ -242,9 +244,9 @@ void DetectorConstruction::SetWorldMagneticField( G4ThreeVector vec ) //expHallMagField->SetFieldValue(G4ThreeVector(vec.x(),vec.y(),vec.z())); } -void DetectorConstruction::SetTabMagneticField(G4String PathAndTableName) +void DetectorConstruction::SetTabMagneticField(G4String PathAndTableName, G4double z_offset, G4double z_rotation) { - nonUniformMagneticField* tabulatedField = new nonUniformMagneticField(PathAndTableName,0); + nonUniformMagneticField* tabulatedField = new nonUniformMagneticField(PathAndTableName,z_offset,z_rotation); } void DetectorConstruction::UpdateGeometry() @@ -740,7 +742,11 @@ void DetectorConstruction::AddDetectionSystemSpice(G4int nRings) G4int segmentID=0; G4double annularDetectorDistance = 115*mm /*+ 150*mm*/; G4ThreeVector pos(0,0,-annularDetectorDistance); + + pSpice->PlaceDetectorMount(logicWorld,pos); + pSpice->PlaceAnnularClamps(logicWorld,pos); pSpice->PlaceGuardRing(logicWorld, pos); + for(int ring = 0; ringBuild(); - G4int NumberSeg = 32; // Segments in Phi - G4int detID=0; - G4double S3DetectorDistance = 21*mm /*+ 150*mm*/; - G4ThreeVector pos(0,0,S3DetectorDistance); - pS3->PlaceGuardRing(logicWorld, pos); + G4int detID=0; + G4int NumberSeg = 32; // Segments in Phi + + G4ThreeVector pos(posX*mm,posY*mm,posZ*mm); + pS3->PlaceS3Mount(logicWorld, pos, AngleOffset); + pS3->PlaceGuardRing(logicWorld, pos); + for(int ring = 0; ringPlaceDetector(logicWorld, pos, ring, Seg, detID); + pS3->PlaceDetector(logicWorld, pos, AngleOffset , ring, Seg, detID); detID++; } // end for(int Seg=0; SegPlaceDetector( logicWorld, ndet ) ; } + +void DetectorConstruction::AddDetectionSystemNew(G4int ndet) +{ + //Create Target Chamber + DetectionSystemNew* pNew = new DetectionSystemNew(); + pNew->Build( logicWorld ); +// pNew->Place( logicWorld) ; + +} + +void DetectorConstruction::AddNewSquareDetector(G4int nDet) +{ + NewSquareDetector* pNewSquare = new NewSquareDetector(); + pNewSquare->Build(); + + + for(int detector = 1; detector<(nDet+1); detector++) + { + pNewSquare->PlaceGuardRing(logicWorld, detector); + + + pNewSquare->PlaceDetector(logicWorld, detector); + + } +} //end NewSquareDetector + + diff --git a/src/DetectorMessenger.cc b/src/DetectorMessenger.cc index 0633551..63f1420 100755 --- a/src/DetectorMessenger.cc +++ b/src/DetectorMessenger.cc @@ -46,6 +46,7 @@ #include "G4UIcommand.hh" #include "G4UIparameter.hh" #include "G4UIcmdWithADouble.hh" +#include "G4String.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -268,22 +269,40 @@ DetectorMessenger::DetectorMessenger(DetectorConstruction* Det) SetSpiceResolutionVariablesCmd = new G4UIcommand("/DetSys/Spice/setResolution",this); SetSpiceResolutionVariablesCmd->SetGuidance("Set Resolution of SPICE Detection System"); - G4UIparameter *parameter1,*parameter2, *parameter3; - G4bool omitable; - parameter1 = new G4UIparameter ("inter", 'd', omitable = false); + G4UIparameter *parameter1,*parameter2;//, *parameter3; + parameter1 = new G4UIparameter ("inter", 'd', false); SetSpiceResolutionVariablesCmd->SetParameter(parameter1); - parameter2 = new G4UIparameter ("gain", 'd', omitable = false); + parameter2 = new G4UIparameter ("gain", 'd', false); SetSpiceResolutionVariablesCmd->SetParameter(parameter2); SetSpiceResolutionVariablesCmd->AvailableForStates(G4State_PreInit,G4State_Idle); - AddDetectionSystemS3Cmd = new G4UIcmdWithAnInteger("/DetSys/det/addS3",this); - AddDetectionSystemS3Cmd->SetGuidance("Add Detection System S3"); + AddDetectionSystemS3Cmd = new G4UIcommand("/DetSys/det/addS3",this); + AddDetectionSystemS3Cmd->SetGuidance("Add Detection System S3 [Nb of rings, xyz position of the center in mm]"); + G4UIparameter *S3parameter1,*S3parameter2, *S3parameter3, *S3parameter4, *S3parameter5; + S3parameter1 = new G4UIparameter ("rings", 'i', true); + AddDetectionSystemS3Cmd->SetParameter(S3parameter1); + S3parameter2 = new G4UIparameter ("posX", 'd', true); + AddDetectionSystemS3Cmd->SetParameter(S3parameter2); + S3parameter3 = new G4UIparameter ("posY", 'd', true); + AddDetectionSystemS3Cmd->SetParameter(S3parameter3); + S3parameter4 = new G4UIparameter ("posZ", 'd', true); + AddDetectionSystemS3Cmd->SetParameter(S3parameter4); + S3parameter5 = new G4UIparameter ("AngleOffset", 'd', true); + AddDetectionSystemS3Cmd->SetParameter(S3parameter5); AddDetectionSystemS3Cmd->AvailableForStates(G4State_PreInit,G4State_Idle); AddDetectionSystemPacesCmd = new G4UIcmdWithAnInteger("/DetSys/det/addPaces",this); AddDetectionSystemPacesCmd->SetGuidance("Add Detection System Paces"); AddDetectionSystemPacesCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + AddDetectionSystemNewCmd = new G4UIcmdWithAnInteger("/DetSys/det/addNew",this); + AddDetectionSystemNewCmd->SetGuidance("Add Detection System New"); + AddDetectionSystemNewCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + AddNewSquareDetectorCmd = new G4UIcmdWithAnInteger("/DetSys/det/addSquare",this); + AddNewSquareDetectorCmd->SetGuidance("Add New Square Detector"); + AddNewSquareDetectorCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + UseTIGRESSPositionsCmd = new G4UIcmdWithABool("/DetSys/det/UseTIGRESSPositions",this); UseTIGRESSPositionsCmd->SetGuidance("Use TIGRESS detector positions rather than GRIFFIN"); UseTIGRESSPositionsCmd->AvailableForStates(G4State_PreInit,G4State_Idle); @@ -334,6 +353,8 @@ DetectorMessenger::~DetectorMessenger() delete AddDetectionSystemSpiceCmd; delete AddDetectionSystemS3Cmd; delete AddDetectionSystemPacesCmd; + delete AddDetectionSystemNewCmd; + delete AddNewSquareDetectorCmd; delete AddDetectionSystemGriffinForwardCmd; delete AddDetectionSystemGriffinForwardDetectorCmd; @@ -371,7 +392,12 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue) Detector->SetWorldMagneticField(WorldMagneticFieldCmd->GetNew3VectorValue(newValue)); } if( command == WorldTabMagneticFieldCmd ) { - Detector->SetTabMagneticField(newValue); + G4String path; + G4double z_offset, zrotation_offset; + const char* s = newValue; + std::istringstream is ((char*)s); + is>>path>>z_offset>>zrotation_offset; + Detector->SetTabMagneticField(path, z_offset, zrotation_offset ); // z in mm, angle in degree } if( command == UpdateCmd ) { Detector->UpdateGeometry(); @@ -501,11 +527,22 @@ void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue) Detector->SetSpiceResolutionVariables(inter,gain); } if( command == AddDetectionSystemS3Cmd ) { - Detector->AddDetectionSystemS3(AddDetectionSystemS3Cmd->GetNewIntValue(newValue)); + G4int rings; + G4double x,y,z, angle_offset; + const char* s = newValue; + std::istringstream is ((char*)s); + is>>rings>>x>>y>>z>>angle_offset; + Detector->AddDetectionSystemS3(rings,x,y,z,angle_offset); // values in mm } if( command == AddDetectionSystemPacesCmd ) { Detector->AddDetectionSystemPaces(AddDetectionSystemPacesCmd->GetNewIntValue(newValue)); } + if( command == AddDetectionSystemNewCmd ) { + Detector->AddDetectionSystemNew(AddDetectionSystemNewCmd->GetNewIntValue(newValue)); + } + if( command == AddNewSquareDetectorCmd ) { + Detector->AddNewSquareDetector(AddNewSquareDetectorCmd->GetNewIntValue(newValue)); + } if( command == UseTIGRESSPositionsCmd ) { Detector->UseTIGRESSPositions(UseTIGRESSPositionsCmd->GetNewBoolValue(newValue)); } diff --git a/src/EventAction.cc b/src/EventAction.cc index f494f1c..db55274 100755 --- a/src/EventAction.cc +++ b/src/EventAction.cc @@ -73,77 +73,75 @@ EventAction::~EventAction() //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void EventAction::BeginOfEventAction(const G4Event* evt) -{ - evtNb = evt->GetEventID(); - if (evtNb%printModulo == 0) -// G4cout << "\n---> Begin of event: " << evtNb << G4endl; - printf( " ---> Ev.# %5d\r", evtNb); - G4cout.flush(); - - RootManager::instance()->SetEventNumber(evtNb); - ClearVariables(); +void EventAction::BeginOfEventAction(const G4Event* evt) { + + //G4cout << " ------------------- \n\n NEW EVENT \n\n -------------------" << evtNb << G4endl ; + evtNb = evt->GetEventID(); + if (evtNb%printModulo == 0) printf( " ---> Ev.# %5d\r", evtNb); + G4cout.flush(); + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void EventAction::EndOfEventAction(const G4Event*) -{ - G4double eventNumber = 0 ; - G4double stepNumber = 0; - G4int cryNumber = 0; - G4int detNumber = 0; - G4double depEnergy = 0; - G4double posx = 0; - G4double posy = 0; - G4double posz = 0; - G4double time = 0; - G4double originDirectionX = 0; - G4double originDirectionY = 0; - G4double originDirectionZ = 0; - G4double originEnergy = 0; - G4int originPdg = 0; - G4int originID = 0; - G4int trackID = 0; - - G4String volume = ""; +{ + + G4double eventNumber = 0 ; + G4double stepNumber = 0; + G4int cryNumber = 0; + G4int detNumber = 0; + G4double depEnergy = 0; + G4double posx = 0; + G4double posy = 0; + G4double posz = 0; + G4double time = 0; + G4double originDirectionX = 0; + G4double originDirectionY = 0; + G4double originDirectionZ = 0; + G4double originEnergy = 0; + G4int originPdg = 0; + G4int originID = 0; + G4int trackID = 0; + G4String volume = ""; for(G4int i = 0; i < MAXSTEPS; i++) { - if(stepTracker[1][i] != 0 && histoManager->GetStepTrackerBool()) { - - eventNumber = stepTracker[0][i] ; - stepNumber = stepTracker[1][i]; - cryNumber = stepTracker[2][i]; - detNumber = stepTracker[3][i]; - depEnergy = stepTracker[4][i]; - posx = stepTracker[5][i]/mm; - posy = stepTracker[6][i]/mm; - posz = stepTracker[7][i]/mm; - time = stepTracker[8][i]/second; - // primary particle info - originDirectionX = stepTracker[9][i]; - originDirectionY = stepTracker[10][i]; - originDirectionZ = stepTracker[11][i]; - originEnergy = stepTracker[12][i]; - originPdg = stepTracker[13][i]; - originID = stepTracker[14][i]; - trackID = stepTracker[15][i]; - - volume = stepVolume[i]; - - histoManager->FillNtuple(eventNumber, stepNumber, cryNumber, detNumber, depEnergy, posx, posy, posz, time ); - RootManager::instance()->FillG4Hit(volume, detNumber, cryNumber, trackID /* this should be particle pdg */, depEnergy, posx, posy, posz, originID, originPdg, originEnergy, originDirectionX, originDirectionY, originDirectionZ); - - //RootManager::instance()->FillHist(1000/keV); //optional - } - + eventNumber = stepTracker[0][i] ; + stepNumber = stepTracker[1][i]; + cryNumber = stepTracker[2][i]; + detNumber = stepTracker[3][i]; + depEnergy = stepTracker[4][i]; + posx = stepTracker[5][i]/mm; + posy = stepTracker[6][i]/mm; + posz = stepTracker[7][i]/mm; + time = stepTracker[8][i]/second; + // primary particle info + originDirectionX = stepTracker[9][i]; + originDirectionY = stepTracker[10][i]; + originDirectionZ = stepTracker[11][i]; + originEnergy = stepTracker[12][i]; + originPdg = stepTracker[13][i]; + originID = stepTracker[14][i]; + trackID = stepTracker[15][i]; + volume = stepVolume[i]; + + histoManager->FillNtuple(eventNumber, stepNumber, cryNumber, detNumber, depEnergy, posx, posy, posz, time ); + //RootManager::instance()->FillHist(1000/keV); //optional + RootManager::instance()->FillG4Hit(volume, detNumber, cryNumber, trackID , // trackID <- this should be particle pdg + depEnergy, posx, posy, posz, + originID, originPdg, originEnergy, + originDirectionX, originDirectionY, originDirectionZ); + + } } - if (depEnergy>0.0) { // if condition satisfied Sort the HitCollection and make a physical event - RootManager::instance()->SortEvent(); + if (1/*depEnergy>0.0*/) { // if condition satisfied Sort the HitCollection and make a physical event + RootManager::instance()->ClearVariables(); + RootManager::instance()->SetHistory( PrimaryInfo ); + RootManager::instance()->SortEvent(evtNb); } - + FillParticleType() ; FillGridEkin() ; FillGriffinCryst() ; @@ -157,14 +155,13 @@ void EventAction::EndOfEventAction(const G4Event*) FillSpiceCryst() ; FillPacesCryst() ; Fill8piCryst() ; - + FillNewCryst() ; //G4int i =0; //histoManager->FillNtuple(stepTracker[0][i], stepTracker[1][i], stepTracker[2][i], stepTracker[3][i], stepTracker[4][i]/keV, stepTracker[5][i]/mm, stepTracker[6][i]/mm, stepTracker[7][i]/mm, stepTracker[8][i]/second ); - //histoManager->FillNtuple(stepTracker[0][i], stepTracker[1][i], stepTracker[2][i], stepTracker[3][i], stepTracker[4][i]/keV, stepTracker[5][i]/mm, stepTracker[6][i]/mm, stepTracker[7][i]/mm, stepTracker[8][i]/second ); + ClearVariables(); - ClearVariables() ; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -175,12 +172,15 @@ void EventAction::EndOfEventAction(const G4Event*) // //histoManager->FillHisto2D(1, x, y); //} -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oo/*oOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void EventAction::ClearVariables() { + + PrimaryInfo.clear(); + if(histoManager->GetStepTrackerBool()) { stepIndex = 0; for (G4int i = 0 ; i < MAXSTEPS; i++) { @@ -215,6 +215,8 @@ void EventAction::ClearVariables() SpiceCrystTrackDet[i] = 0 ; PacesCrystEnergyDet[i] = 0 ; PacesCrystTrackDet[i] = 0 ; + NewCrystTrackDet[i] =0 ; + NewCrystEnergyDet[i] =0 ; } @@ -241,6 +243,143 @@ void EventAction::ClearVariables() } + void EventAction::SetPrimaryInfo(TrackInformation * info ){ + + /* + bool family = false ; + bool daughter = false ; + bool copy = false ; + + unsigned iSize = 0; + unsigned nSize = PrimaryInfo.size() ; + */ + + //cout << " =============== Candidate info " << endl ; + //cout << " Adress " << info << endl ; + //info->Print(); + //cout << " ================================ " << endl ; + + + /* cout << " ++++++++++++++++++++ content At the begining : " << nSize << endl ; + for ( unsigned iSize = 0 ; iSize < nSize ; iSize++) { + cout << " Adress " << PrimaryInfo.at(iSize) << endl ; + PrimaryInfo.at(iSize)->Print(); + } + cin.get(); + + cout << " original parent current " << endl ; + cout << " new " << info->GetOriginalTrackID() << " " ; + cout << " " << info->GetParentTrackID() << " " ; + cout << " " << info->GetCurrentTrackID() << endl ; + */ + + /* + for ( iSize = 0 ; iSize < nSize ; iSize++) { + + + // cout << " old " << PrimaryInfo.at(iSize)->GetOriginalTrackID() << " " ; + // cout << " " << PrimaryInfo.at(iSize)->GetParentTrackID() << " " ; + // cout << " " << PrimaryInfo.at(iSize)->GetCurrentTrackID() << endl ; + + if (PrimaryInfo.at(iSize)->GetOriginalTrackID()==info->GetOriginalTrackID() ) { // same family + //cout << "Same family" << endl ; + family = true ; + } + //else { cout << "Diff. family" << endl ;} + + if (family && info->GetParentTrackID() == PrimaryInfo.at(iSize)->GetCurrentTrackID()) { // daughter contains all the history up to the family tree + //cout << "Daughter -> Update info " << endl ; + PrimaryInfo.at(iSize) = new TrackInformation(info) ; // replace this info with the updated one + //PrimaryInfo.at(iSize)->Print(); + daughter = true ; + copy = false ; + break; + } + //else { cout << "Not Daughter" << endl ;} + + if (family && !daughter && info->GetCurrentTrackID() == PrimaryInfo.at(iSize)->GetCurrentTrackID() ) { // if it's the same information skip! + //cout << "Same Information" << endl ; + copy = true ; + break; + } + //else { cout << " Not Same Information" << endl ;} + + } + */ + + //if (nSize>0 && iSize == nSize) cout << " at end of loop iSize : " << iSize << " nSize : "<< nSize<< endl ; + + //if its a new track add it to the vector + //if (!daughter && !copy) { + //cout << ">>>>>>>>>>>>>>>>>>>>> Add this info to the vector" << endl ; + PrimaryInfo.push_back( new TrackInformation(info)); + // } + + /* + cout << " ++++++++++++++++++++ content size : " << PrimaryInfo.size() << endl ; + + for ( unsigned iSize = 0 ; iSize < PrimaryInfo.size() ; iSize++) { + PrimaryInfo.at(iSize)->Print(); + } + cin.get(); + + */ + + } + +void EventAction::AddStepTracker(G4double eventNumber, G4double stepNumber, G4String volume, + G4double cryNumber, G4double detNumber, + G4double depEnergy, G4double posx, G4double posy, G4double posz, G4double time, + G4double originDirectionX, G4double originDirectionY, G4double originDirectionZ, + G4double originEnergy, G4int originPdg, G4int originID, G4int trackID) { + + /*cout << "eventNumber " << eventNumber << endl + << "step Number " << stepNumber << endl ; + << "volume " << volume << endl + << "cryNumber " << cryNumber << endl + << "detNumber " << detNumber << endl + << "depEnergy " << depEnergy << endl + << "posx " << posx << endl + << "posy " << posy << endl + << "posz " << posz << endl + << "time " << time << endl + << "originDirectionX " << originDirectionX << endl + << "originDirectionY " << originDirectionY << endl + << "originDirectionZ " << originDirectionZ << endl + << "originEnergy " << originEnergy << endl + << "originPdg " << originPdg << endl + << "originID " << originID << endl + << "trackID " << trackID << endl ; + */ + if(histoManager->GetStepTrackerBool()) { + stepTracker[0][stepIndex] = eventNumber; + stepTracker[1][stepIndex] = stepNumber; + stepTracker[2][stepIndex] = cryNumber; + stepTracker[3][stepIndex] = detNumber; + stepTracker[4][stepIndex] = depEnergy; + stepTracker[5][stepIndex] = posx; + stepTracker[6][stepIndex] = posy; + stepTracker[7][stepIndex] = posz; + stepTracker[8][stepIndex] = time; + stepTracker[9][stepIndex] = originDirectionX; + stepTracker[10][stepIndex] = originDirectionY; + stepTracker[11][stepIndex] = originDirectionZ; + stepTracker[12][stepIndex] = originEnergy; + stepTracker[13][stepIndex] = originPdg; + stepTracker[14][stepIndex] = originID; + stepTracker[15][stepIndex] = trackID; + + stepVolume[stepIndex] = volume ; + + stepIndex++; + if(stepIndex == MAXSTEPS) { + G4cout << "\n ----> error 13423549 \n" << G4endl; + exit(1); + } + }; + } + + void EventAction::FillParticleType() { G4int numParticleTypes = 0; @@ -468,19 +607,10 @@ void EventAction::FillPacesCryst() } } -//void AddStepTracker(G4int eventNumber, G4int stepNumber, G4int cryNumber, G4int detNumber, G4double depEnergy, G4double posx, G4double posy, G4double posz, G4double time) -//{ -// stepTracker[0][stepIndex] = (G4double)(eventNumber); -// stepTracker[1][stepIndex] = (G4double)(stepNumber); -// stepTracker[2][stepIndex] = (G4double)(cryNumber); -// stepTracker[3][stepIndex] = (G4double)(detNumber); -// stepTracker[4][stepIndex] = depEnergy; -// stepTracker[5][stepIndex] = posx; -// stepTracker[6][stepIndex] = posy; -// stepTracker[7][stepIndex] = posz; -// stepTracker[8][stepIndex] = time; -// stepIndex++; - +void EventAction::FillNewCryst() +{ + +} //} diff --git a/src/HistoManager.cc b/src/HistoManager.cc index 3e3e0b5..09a165f 100755 --- a/src/HistoManager.cc +++ b/src/HistoManager.cc @@ -333,6 +333,20 @@ void HistoManager::book() filename = "paces_crystal_edep_det" + detString; MakeHisto(analysisManager, filename, title, xmin, xmax, nbins); } + + // new detector + filename = "new_crystal_edep"; + MakeHisto(analysisManager, filename, title, xmin, xmax, nbins); + + filename = "new_crystal_edep_sum"; + MakeHisto(analysisManager, filename, title, xmin, xmax, nbins); + + for (G4int i=0; i < MAXNUMDET; i++) { + detString = G4intToG4String(i); + + filename = "new_crystal_edep_det" + detString; + MakeHisto(analysisManager, filename, title, xmin, xmax, nbins); + } } /////////////////////////////////////////////////////////////////// diff --git a/src/MapEvent.cc b/src/MapEvent.cc index 4d4b4e3..9d21426 100644 --- a/src/MapEvent.cc +++ b/src/MapEvent.cc @@ -14,13 +14,11 @@ MapEvent::~MapEvent(void) void MapEvent::FillVectors( int pdg, double Energy, // depositid energy double Px, double Py, double Pz, // position vector - int ID, // original(primary) TrackID int PrimPdg,// primary particle pdg PDG encoding double PrimEnergy,//original(primary) energy - double Mx, double My, double Mz) // primary particle momentum vector - - { + double Mx, double My, double Mz) { // primary particle momentum vector + fHCPdg.push_back(pdg); fHCEnergy.push_back(Energy) ; fHCPosition.push_back(TVector3(Px,Py,Pz)) ; diff --git a/src/NewSquareDetector.cc b/src/NewSquareDetector.cc new file mode 100644 index 0000000..a5a6094 --- /dev/null +++ b/src/NewSquareDetector.cc @@ -0,0 +1,561 @@ +#include "DetectorConstruction.hh" +#include "DetectorMessenger.hh" + +#include "G4Material.hh" + +#include "G4Tubs.hh" +#include "G4Box.hh" +#include "G4Sphere.hh" +#include "G4LogicalVolume.hh" +#include "G4PVPlacement.hh" +#include "G4SubtractionSolid.hh" + +#include "G4GeometryManager.hh" +#include "G4PhysicalVolumeStore.hh" +#include "G4LogicalVolumeStore.hh" +#include "G4SolidStore.hh" +#include "G4AssemblyVolume.hh" + +#include "G4VisAttributes.hh" +#include "G4Colour.hh" + +#include "NewSquareDetector.hh" + + +NewSquareDetector::NewSquareDetector() + + +{ + + this->detector_alignment = 1; // 1 for parallel to beam axis, 2 for perpindicular to beam axis + + + /////////////////////////////// + // Square Detector Materials // + /////////////////////////////// + + this->detector_material = "Silicon"; + + + if(this->detector_alignment == 1){ + + //////////////////////////// + // Dimensions for Crystal // + //////////////////////////// + + this->square_segment_element_length = 25*mm ; //50*mm for perpindicular to beam axis detector size, 25mm for parallel to beam axis detector + this->square_segment_element_width = 55*mm ; // 36*mm for perpindicular to beam axis detector size , 50mm for parallel to beam axis detector + this->square_segment_thickness = 5*mm ; //thickness in X-Y plane 5mm for perpindicular to beam axis and 1mm possible for parallel to beam axis + + ////////////////////////////////////////////////// + // How far back the detector is from the target // + ////////////////////////////////////////////////// + + this->SquareDetectorDistance = - 93*mm; // -75*mm for perpindicular to beam axis detector size, -93mm for parallel to beam axis detector + this->SquareDetectorRotationXY = 22.5*deg; + this->SquareDetectorRotationZ = 22.5*deg; + + this->SquareDetectorCorrection = ((this->square_segment_element_width*sin(this->SquareDetectorRotationZ))+(this->square_segment_thickness*sin(this->SquareDetectorRotationZ)) )*mm; + + } + + else{ + + //////////////////////////// + // Dimensions for Crystal // + //////////////////////////// + + this->square_segment_element_length = 50*mm ; //50*mm for perpindicular to beam axis detector size, 25mm for parallel to beam axis detector + this->square_segment_element_width = 36*mm ; // 36*mm for perpindicular to beam axis detector size , 50mm for parallel to beam axis detector + this->square_segment_thickness = 5*mm ; //thickness in X-Y plane 5mm for perpindicular to beam axis and 1mm possible for parallel to beam axis + + ////////////////////////////////////////////////// + // How far back the detector is from the target // + ////////////////////////////////////////////////// + + this->SquareDetectorDistance = - 75*mm; // -75*mm for perpindicular to beam axis detector size, -93mm for parallel to beam axis detector + this->SquareDetectorRotationXY = 22.5*deg; + + } + + + + ////////////////////////////////////// + // No. of Segments in x-y dimension // + ////////////////////////////////////// + + this->no_x_segments = 4. ; + this->no_y_segments = 3. ; + this->no_detectors = 4. ; + + ////////////////////////////// + // Dimensions of Guard Ring // + ////////////////////////////// + + this->square_guard_ring_depth = 6*mm ; + this->square_guard_ring_thickness = this->square_segment_thickness ; + + + + ///////////////////////////////////////////////////// + // Relative position in X-Y Plane around Beam Axis // + ///////////////////////////////////////////////////// + + //Guard Ring of Detector + this->BeamPipeXDistanceGR = -7*mm; + this->BeamPipeYDistanceGR = -7*mm; + + // Active Area of Detector + this->BeamPipeXDistanceD = (BeamPipeXDistanceGR + square_guard_ring_depth)*mm; + this->BeamPipeYDistanceD = (BeamPipeYDistanceGR - square_guard_ring_depth)*mm; + +} + + +//////////////////////////// +// Delete Logical Volumes // +//////////////////////////// + +NewSquareDetector::~NewSquareDetector() +{ + + delete [] SquareDetect_log; + delete [] square_guard_ring_log; + +} + + +///////////////////////////////////////// +// Main Build Function in Construction // +///////////////////////////////////////// + +G4int NewSquareDetector::Build() +{ + + //Build assembly volumes + this->assembly = new G4AssemblyVolume(); + + if(this->detector_alignment == 1){ + for(int detectorID=1; detectorID<(this->no_detectors+1); detectorID++) + { + + this->assemblySquareDet[detectorID] = new G4AssemblyVolume(); + this->assemblyGuardRing[detectorID] = new G4AssemblyVolume(); + //Build logical segment of square detector + BuildDetectorFace1(detectorID); + BuildGuardRing1(detectorID); + } + + } + + else{ + for(int detectorID=1; (this->no_detectors+1); detectorID++) + { + + this->assemblySquareDet[detectorID] = new G4AssemblyVolume(); + this->assemblyGuardRing[detectorID] = new G4AssemblyVolume(); + //Build logical segment of square detector + BuildDetectorFace2(detectorID); + BuildGuardRing2(detectorID); + } + + } + + + return 1; +} // end Build + + +///////////////////////////////////////////////////////////// +// Building a Singular Segment of a Singular Detector Face // +///////////////////////////////////////////////////////////// + +G4Box* NewSquareDetector::BuildSegment1() +{ + + //Define the geometry of a single segment + G4double segment_element_length = ((this->square_segment_element_length - (this->square_guard_ring_depth*2))/no_x_segments)/2; + G4double segment_element_width = ((this->square_segment_element_width - (this->square_guard_ring_depth*2))/no_y_segments)/2; + G4double segment_thickness = this->square_segment_thickness/2; + + //Creating the logical volume of a single segment + G4Box* crystal_segment = new G4Box("cystal_segment", segment_element_length, segment_thickness, segment_element_width); + + return crystal_segment; + +}//end Build Segment + + +G4Box* NewSquareDetector::BuildSegment2() +{ + + //Define the geometry of a single segment + G4double segment_element_length = ((this->square_segment_element_length - (this->square_guard_ring_depth*2))/no_x_segments)/2; + G4double segment_element_width = ((this->square_segment_element_width - (this->square_guard_ring_depth*2))/no_y_segments)/2; + G4double segment_thickness = this->square_segment_thickness/2; + + + //Creating the logical volume of a single segment + G4Box* crystal_segment = new G4Box("cystal_segment", segment_element_length, segment_element_width, segment_thickness); + + return crystal_segment; + +}//end Build Segment + +//////////////////////////////////////////////////////////////// +// Building a Singular Guard Ring of a Singular Detector Face // +//////////////////////////////////////////////////////////////// + +G4int NewSquareDetector::BuildGuardRing1(G4int detectorID) +{ + + // Set visualization attributes + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(GUARD_COL)); + vis_att->SetVisibility(true); + + //Define the material + G4Material* material = G4Material::GetMaterial(this->detector_material); + + //Define Outer Edge Geometry of the Guard Ring + G4double guard_ring_outer_length = this->square_segment_element_length/2 ; + G4double guard_ring_outer_width = this->square_segment_element_width /2; + G4double guard_ring_outer_thickness = this->square_segment_thickness /2; + + //Define Inner Edge Geometry of the Guard Ring + G4double guard_ring_inner_length = (this->square_segment_element_length - (this->square_guard_ring_depth *2))/2 ; + G4double guard_ring_inner_width = (this->square_segment_element_width - (this->square_guard_ring_depth *2))/2 ; + G4double guard_ring_inner_thickness = (this->square_guard_ring_thickness)/2 ; + + //Create the two Sections of a Single Guard Ring + G4Box* outer_guard_ring = new G4Box("outer_guard_ring", guard_ring_outer_length, guard_ring_outer_thickness, guard_ring_outer_width); + G4Box* inner_guard_ring = new G4Box("inner_guard_ring", guard_ring_inner_length, guard_ring_inner_thickness, guard_ring_inner_width); + + //Defining the Combination of the two Sections forming a singular Guard Ring + G4ThreeVector Trans = G4ThreeVector(0,0,0); + + G4RotationMatrix* rot = new G4RotationMatrix(); + rot->rotateZ(0*deg); + G4SubtractionSolid* square_guard_ring = new G4SubtractionSolid("squard_guard_ring", outer_guard_ring, inner_guard_ring,rot, Trans); + + square_guard_ring_log[detectorID] = new G4LogicalVolume(square_guard_ring, material, "square_guard_ring",0,0,0); + square_guard_ring_log[detectorID]->SetVisAttributes(vis_att); + + // Define rotation and movement of a singular Guard Ring before creating logical volume + G4ThreeVector xdirection = G4ThreeVector(1,0,0); + G4double x_position = (this->square_guard_ring_thickness/2)*mm +5*mm ;// - this->square_guard_ring_depth+2*mm; + G4ThreeVector xmove = x_position * xdirection; + + G4ThreeVector zdirection = G4ThreeVector(0,0,1); + G4double z_position = (this->square_segment_element_width/2) -2*mm ; + G4ThreeVector zmove = z_position * zdirection; + + G4ThreeVector ydirection = G4ThreeVector(0,1,0); + G4double y_position = (this->square_segment_element_length/2) +1*mm ; + G4ThreeVector ymove = y_position * ydirection; + + + G4ThreeVector move = (xmove + ymove + zmove); + G4RotationMatrix* rotate = new G4RotationMatrix; + rotate->rotateX(0*deg); + + + //Creating the logical volume of singular guard ring + + this->assemblyGuardRing[detectorID]->AddPlacedVolume(square_guard_ring_log[detectorID], move, rotate); + + return 1; + +}//end Build Guard Ring 1 + + +G4int NewSquareDetector::BuildGuardRing2(G4int detectorID) +{ + + // Set visualization attributes + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(GUARD_COL)); + vis_att->SetVisibility(true); + + //Define the material + G4Material* material = G4Material::GetMaterial(this->detector_material); + + //Define Outer Edge Geometry of the Guard Ring + G4double guard_ring_outer_length = this->square_segment_element_length/2 ; + G4double guard_ring_outer_width = this->square_segment_element_width /2; + G4double guard_ring_outer_thickness = this->square_segment_thickness /2; + + //Define Inner Edge Geometry of the Guard Ring + G4double guard_ring_inner_length = (this->square_segment_element_length - (this->square_guard_ring_depth *2))/2 ; + G4double guard_ring_inner_width = (this->square_segment_element_width - (this->square_guard_ring_depth *2))/2 ; + G4double guard_ring_inner_thickness = (this->square_guard_ring_thickness)/2 ; + + //Create the two Sections of a Single Guard Ring + G4Box* outer_guard_ring = new G4Box("outer_guard_ring", guard_ring_outer_length, guard_ring_outer_width, guard_ring_outer_thickness); + G4Box* inner_guard_ring = new G4Box("inner_guard_ring", guard_ring_inner_length, guard_ring_inner_width, guard_ring_inner_thickness); + + //Defining the Combination of the two Sections forming a singular Guard Ring + G4ThreeVector Trans = G4ThreeVector(0,0,0); + + G4RotationMatrix* rot = new G4RotationMatrix(); + rot->rotateZ(0*deg); + G4SubtractionSolid* square_guard_ring = new G4SubtractionSolid("squard_guard_ring", outer_guard_ring, inner_guard_ring,rot, Trans); + + square_guard_ring_log[detectorID] = new G4LogicalVolume(square_guard_ring, material, "square_guard_ring",0,0,0); + square_guard_ring_log[detectorID]->SetVisAttributes(vis_att); + + // Define rotation and movement of a singular Guard Ring before creating logical volume + G4double BeamPipeXDistance = this->BeamPipeXDistanceGR; + G4double BeamPipeYDistance = this->BeamPipeYDistanceGR; + + G4ThreeVector xdirection = G4ThreeVector(1,0,0); + G4double x_position = (this->square_segment_element_length/2) + BeamPipeXDistanceGR ; + G4ThreeVector xmove = x_position * xdirection; + + G4ThreeVector ydirection = G4ThreeVector(0,1,0); + G4double y_position = -(this->square_segment_element_width/2.)+BeamPipeYDistanceGR; + G4ThreeVector ymove = y_position * ydirection; + + G4ThreeVector zdirection = G4ThreeVector(0,0,1); + G4double z_position = -(this->square_segment_thickness/2.); + G4ThreeVector zmove = z_position * zdirection; + + + G4ThreeVector move = (xmove + ymove + zmove); + G4RotationMatrix* rotate = new G4RotationMatrix; + rotate->rotateX(0*deg); + + //cout << "moveGR" << move << endl; + + //Creating the logical volume of singular guard ring + + this->assemblyGuardRing[detectorID]->AddPlacedVolume(square_guard_ring_log[detectorID], move, rotate); + + return 1; + +}//end Build Guard Ring 2 + +/////////////////////////////////////// +// Building a Singular Detector Face // +/////////////////////////////////////// + +G4int NewSquareDetector::BuildDetectorFace1(G4int detectorID) +{ + + //Define the material to be used + G4Material* material = G4Material::GetMaterial(this->detector_material); + + //Set Visualisation Attributes + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(SIL_COL)); + vis_att->SetVisibility(true); + + // Define rotation and movement of the first segment away from the origin + G4ThreeVector xdirection = G4ThreeVector(1,0,0); + G4double x_position = -((this->square_segment_element_length/2) - this->square_guard_ring_depth - (((this->square_segment_element_length - (this->square_guard_ring_depth *2))/this->no_x_segments)/2) - this->square_guard_ring_thickness/2) *mm +5*mm; + + G4ThreeVector xalign = x_position * xdirection; + + G4ThreeVector zdirection = G4ThreeVector(0,0,1); + G4double z_position =(((this->square_segment_element_width )/this->no_y_segments)/2)+(this->square_guard_ring_depth/2)+(this->square_guard_ring_depth/6)*mm-2*mm; + G4ThreeVector zalign = z_position * zdirection; + + G4ThreeVector ydirection = G4ThreeVector(0,1,0); + G4double y_position = (this->square_segment_element_length/2)*mm +1*mm ; + G4ThreeVector yalign = y_position * ydirection; + + + G4ThreeVector align = xalign + yalign + zalign; + G4RotationMatrix* rotate = new G4RotationMatrix(); + rotate->rotateZ(0*deg); + + //Build the logical volume of a single segment + + G4Box* SquareDetect = BuildSegment1(); + + // construct logical volume + + G4String detectName = "SquareDetect_"; + G4String detectorNo = G4UIcommand::ConvertToString(detectorID); + detectName += detectorNo; + detectName += "_Log"; + + SquareDetect_log[detectorID] = new G4LogicalVolume(SquareDetect, material, detectName, 0, 0, 0); + SquareDetect_log[detectorID]->SetVisAttributes(vis_att); + + //Using for loops to define the array of segments for one detector crystal + + for(int xseg = 0; xsegno_x_segments; xseg++) + { + for(int zseg=0; zsegno_y_segments; zseg++) + { + + G4ThreeVector xdirection = G4ThreeVector(1,0,0); + G4double xposition = ((this->square_segment_element_length - (this->square_guard_ring_depth *2))/this->no_x_segments)*mm; + G4ThreeVector xmove = xposition * xdirection ; + + G4ThreeVector zdirection = G4ThreeVector(0,0,1); + G4double zposition = ((this->square_segment_element_width - (this->square_guard_ring_depth *2))/this->no_y_segments)*mm; + G4ThreeVector zmove = zposition * zdirection; + + G4ThreeVector allmove = ((xmove*xseg) + (zmove*zseg)); + + G4ThreeVector move = align + allmove; + + //Create the logical volume of the array + this->assemblySquareDet[detectorID]->AddPlacedVolume(SquareDetect_log[detectorID], move, rotate); + + } + } + + return 1; + +} //end build + +G4int NewSquareDetector::BuildDetectorFace2(G4int detectorID) +{ + + //Define the material to be used + G4Material* material = G4Material::GetMaterial(this->detector_material); + + //Set Visualisation Attributes + G4VisAttributes* vis_att = new G4VisAttributes(G4Colour(SIL_COL)); + vis_att->SetVisibility(true); + + // Define rotation and movement of the first segment away from the origin + G4double BeamPipeXDistanceD = this->BeamPipeXDistanceD; + G4double BeamPipeYDistanceD = this->BeamPipeYDistanceD; + + G4ThreeVector xdirection = G4ThreeVector(1,0,0); + G4double x_position = (((this->square_segment_element_length - (this->square_guard_ring_depth*2))/this->no_x_segments)/2) + BeamPipeXDistanceD; + G4ThreeVector xrealign = x_position * xdirection; + + G4ThreeVector ydirection = G4ThreeVector(0,1,0); + G4double y_position =(-(((this->square_segment_element_width- (this->square_guard_ring_depth*2))/this->no_y_segments)/2 )) + BeamPipeYDistanceD; + G4ThreeVector yrealign = y_position * ydirection; + + G4ThreeVector zdirection = G4ThreeVector(0,0,1); + G4double z_position = -(this->square_segment_thickness/2.); + G4ThreeVector zrealign = z_position * zdirection; + + G4ThreeVector realign = (xrealign + yrealign + zrealign); + G4RotationMatrix* rotate = new G4RotationMatrix(); + rotate->rotateZ(0*deg); + + //Build the logical volume of a single segment + + G4Box* SquareDetect = BuildSegment2(); + + // construct logical volume + + G4String detectName = "SquareDetect_"; + G4String detectorNo = G4UIcommand::ConvertToString(detectorID); + detectName += detectorNo; + detectName += "_Log"; + + SquareDetect_log[detectorID] = new G4LogicalVolume(SquareDetect, material, detectName, 0, 0, 0); + SquareDetect_log[detectorID]->SetVisAttributes(vis_att); + + //Using for loops to define the array of segments for one detector crystal + + for(int xseg = 0; xsegno_x_segments; xseg++) + { + for(int yseg=0; ysegno_y_segments; yseg++) + { + + G4ThreeVector xdirection = G4ThreeVector(1,0,0); + G4double xposition = ((this->square_segment_element_length - (this->square_guard_ring_depth *2))/this->no_x_segments)*mm; + G4ThreeVector xmove = xposition * xdirection ; + + G4ThreeVector ydirection = G4ThreeVector(0,-1,0); + G4double yposition = ((this->square_segment_element_width - (this->square_guard_ring_depth *2))/this->no_y_segments)*mm; + G4ThreeVector ymove = yposition * ydirection; + + G4ThreeVector allmove = ((xmove*xseg) + (ymove*yseg)); + + G4ThreeVector move = realign + allmove; + + //Create the logical volume of the array + this->assemblySquareDet[detectorID]->AddPlacedVolume(SquareDetect_log[detectorID], move, rotate); + + } + } + + return 1; + +} //end build2 + +////////////////////////////////////////////////// +// "place" function called in DetectorMessenger // +////////////////////////////////////////////////// + +G4int NewSquareDetector::PlaceDetector(G4LogicalVolume* exp_hall_log, G4int detector) +{ + + //Define movement of Logical volume in relation to target/source (z-axis) + G4double SquareDetectorTargetDistance = this->SquareDetectorDistance; + G4double SquareDetectorCorrection = this->SquareDetectorCorrection; + G4ThreeVector move = TranslateDetectors(detector,SquareDetectorCorrection, SquareDetectorTargetDistance); + + //Define the rotation to create the detector plates + G4RotationMatrix* rotate = new G4RotationMatrix; + + if(this->detector_alignment == 1){ + + rotate->rotateX(this->SquareDetectorRotationZ); + + rotate->rotateZ((((360.*deg/this->no_detectors) * detector)-this->SquareDetectorRotationXY));} + + + else{ + rotate->rotateZ((((360.*deg/this->no_detectors) * detector)+this->SquareDetectorRotationXY));} + + // Build physical volume of the detectors + assemblySquareDet[detector]->MakeImprint(exp_hall_log, move, rotate, detector); + + return 1; +} //end place + +G4int NewSquareDetector::PlaceGuardRing(G4LogicalVolume* exp_hall_log , G4int detector) +{ + + //Define movement of Logical volume in relation to target/source (z-axis) + G4double SquareDetectorTargetDistance = this->SquareDetectorDistance; + G4double SquareDetectorCorrection = this->SquareDetectorCorrection; + G4ThreeVector move = TranslateDetectors(detector,SquareDetectorCorrection, SquareDetectorTargetDistance); + + + //Define the rotation to create the detector guard rings + G4RotationMatrix* rotate = new G4RotationMatrix; + + if(this->detector_alignment == 1){ + + + rotate->rotateX(this->SquareDetectorRotationZ); + rotate->rotateZ((((360.*deg/this->no_detectors) * detector)-this->SquareDetectorRotationXY)); + } + else{ + rotate->rotateZ((((360.*deg/this->no_detectors) * detector)+this->SquareDetectorRotationXY));} + + // Build physical volume of the guard rings + assemblyGuardRing[detector]->MakeImprint(exp_hall_log, move, rotate, detector); + + return 1; +} //end place + + + + + +G4ThreeVector NewSquareDetector::TranslateDetectors(G4int copyID, G4double movement, G4double z_position) +{ + G4double x_position(0); + G4double y_position(0); + x_position = movement*sin(-((copyID)*90.*deg)); + y_position = movement*cos((copyID)*90.*deg); + return G4ThreeVector(x_position, y_position, z_position); +} + + + + + + + diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc index 12f598e..779de6d 100755 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -148,14 +148,14 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) G4double monteCarloRand = UniformRand48(); // energy distribution - for(G4int i = 0 ; i < energyDist.size() ; i++ ) + for(unsigned i = 0 ; i < energyDist.size() ; i++ ) { if(monteCarloRand <= monteCarlo[i]) { selectedEnergy = energyDist[i]*1.0*keV; break; } - } // end for(G4int i + } // end for(unsigned i energy = selectedEnergy; @@ -204,7 +204,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) myGammaAndICParticle->GenerateParticles(this->outputTimeInSeconds); // generate gamma and IC partile decay //myGammaAndICParticle->PrintAllGeneratedParticles(); numberOfParticles = myGammaAndICParticle->GetNumberOfGeneratedParticles(); - for( G4int i = 0 ; i < numberOfParticles ; i++ ) // get the number of particles in decay and loop over them + for( unsigned i = 0 ; i < numberOfParticles ; i++ ) // get the number of particles in decay and loop over them { if(emit_gamma_ic_flag) { @@ -525,12 +525,12 @@ void PrimaryGeneratorAction::SetBetaPlusEmission( G4String file ) exit (1); } - for(G4int i = 0 ; i < weightDist.size() ; i++ ) + for(unsigned i = 0 ; i < weightDist.size() ; i++ ) { sumWeight = sumWeight + weightDist[i]; } - for(G4int i = 0 ; i < weightDist.size() ; i++ ) + for(unsigned i = 0 ; i < weightDist.size() ; i++ ) { percentage = weightDist[i]/sumWeight; if(i == 0) @@ -564,12 +564,12 @@ void PrimaryGeneratorAction::SetBetaMinusEmission( G4String file ) exit (1); } - for(G4int i = 0 ; i < weightDist.size() ; i++ ) + for(unsigned i = 0 ; i < weightDist.size() ; i++ ) { sumWeight = sumWeight + weightDist[i]; } - for(G4int i = 0 ; i < weightDist.size() ; i++ ) + for(unsigned i = 0 ; i < weightDist.size() ; i++ ) { percentage = weightDist[i]/sumWeight; if(i == 0) @@ -609,13 +609,17 @@ void PrimaryGeneratorAction::EmissionForRadioactiveSourceDecay( G4Event* myEvent // RandFlat::shoot(m,n) excludes m and n itself! betaDecayRandomiser = RandFlat::shoot(0.,100.); levelProbSum=0; + + + while( (levelProbSum += levelScheme[i][4]) < betaDecayRandomiser ) { i++; + // make sure that i is running only through table elements if( i==nLevels ) { - G4cout << "Total beta decay probability found = " <> nLevels >> nTransPerLevel >> nParam >> bindingEnergyK >> bindingEnergyL1 >> bindingEnergyL2 >> bindingEnergyL3; @@ -745,8 +749,9 @@ void PrimaryGeneratorAction::LevelSchemeReader( const char* filename ) G4cout << " [binding energies "<< bindingEnergyK<< " " << bindingEnergyL1 << " " <> KXRayEnergy[k]; KXRayOrigin[k]=shellFrom[i]; + fKXRayEnergy.push_back(KXRayEnergy[k]); fKXRayOrigin.push_back(shellFrom[i]); @@ -773,7 +779,8 @@ void PrimaryGeneratorAction::LevelSchemeReader( const char* filename ) { file >> L1XRayEnergy[l1]; L1XRayOrigin[l1]=shellFrom[i]; - + + fL1XRayEnergy.push_back(L1XRayEnergy[l1]); fL1XRayOrigin.push_back(shellFrom[i]); @@ -783,7 +790,8 @@ void PrimaryGeneratorAction::LevelSchemeReader( const char* filename ) { file >> L2XRayEnergy[l2]; L2XRayOrigin[l2]=shellFrom[i]; - + + fL2XRayEnergy.push_back(L2XRayEnergy[l2]); fL2XRayOrigin.push_back(shellFrom[i]); @@ -793,7 +801,8 @@ void PrimaryGeneratorAction::LevelSchemeReader( const char* filename ) { file >> L3XRayEnergy[l3]; L3XRayOrigin[l3]=shellFrom[i]; - + + fL3XRayEnergy.push_back(L3XRayEnergy[l3]); fL3XRayOrigin.push_back(shellFrom[i]); @@ -801,7 +810,7 @@ void PrimaryGeneratorAction::LevelSchemeReader( const char* filename ) } } - getline(file,buffer); G4cout <<"skipping this line :" <> dumDouble; fAugerIntensity.push_back(dumDouble) ; } + { file >> dumDouble; + + fAugerIntensity.push_back(dumDouble) ; + } //skip line - getline(file,buffer); G4cout <<"skipping this line :" <> dumDouble; fAugerEnergy.push_back(dumDouble) ; } + { file >> dumDouble; + + fAugerEnergy.push_back(dumDouble) ; + } // set up storage space for table for magnetic field @@ -854,10 +869,10 @@ void PrimaryGeneratorAction::LevelSchemeReader( const char* filename ) G4cout << "Level scheme storage set." << endl; // skip another 3 lines - getline(file,buffer); G4cout <<"skipping this line :" <= xOrARandom) { - G4int i = 0; + unsigned i = 0; addedIntensity = fKXRayIntensity.at(0); while( addedIntensity < xOrARandom ) { @@ -1041,7 +1056,7 @@ void PrimaryGeneratorAction::EmissionForVacantShell(int shell, G4Event* myEvent) } else if( (totalXRayIntensity+totalAugerIntensity) >= xOrARandom) { - G4int i = 0; + unsigned i = 0; addedIntensity = totalXRayIntensity+fAugerIntensity.at(0); while( addedIntensity < xOrARandom ) { @@ -1071,14 +1086,14 @@ void PrimaryGeneratorAction::EmissionForVacantShell(int shell, G4Event* myEvent) { // sum up the X-ray intensities and Auger intensities totalXRayIntensity = 0.; - for(G4int i=0;i xRayRandom) { - G4int i = 0; + unsigned i = 0; addedIntensity = fL1XRayIntensity.at(0); while( addedIntensity < xRayRandom ) { @@ -1106,14 +1121,14 @@ void PrimaryGeneratorAction::EmissionForVacantShell(int shell, G4Event* myEvent) { // sum up the X-ray intensities and Auger intensities totalXRayIntensity = 0.; - for(G4int i=0;i xRayRandom) { - G4int i = 0; + unsigned i = 0; addedIntensity = fL2XRayIntensity.at(0); while( addedIntensity < xRayRandom ) { @@ -1139,14 +1154,14 @@ void PrimaryGeneratorAction::EmissionForVacantShell(int shell, G4Event* myEvent) { // sum up the X-ray intensities and Auger intensities totalXRayIntensity = 0.; - for(G4int i=0;i xRayRandom) { - G4int i = 0; + unsigned i = 0; addedIntensity = fL3XRayIntensity.at(0); while( addedIntensity < xRayRandom ) { diff --git a/src/RawG4Event.cc b/src/RawG4Event.cc index bd73a19..82ea2b4 100644 --- a/src/RawG4Event.cc +++ b/src/RawG4Event.cc @@ -7,8 +7,7 @@ RawG4Event::~RawG4Event(void){ } -void RawG4Event::FillVectors( - int pdg, +void RawG4Event::FillVectors( int pdg, double Energy, // depositid energy int detector, // detector int crystal, // crystal @@ -25,7 +24,6 @@ void RawG4Event::FillVectors( fHCDetector.push_back(detector) ; fHCCrystal.push_back(crystal) ; - fHCPrimaryID.push_back(ID) ; fHCPrimaryPdg.push_back(PrimPdg) ; fHCPrimaryEnergy.push_back(PrimEnergy) ; diff --git a/src/RootManager.cc b/src/RootManager.cc index beddbe7..5175159 100644 --- a/src/RootManager.cc +++ b/src/RootManager.cc @@ -19,26 +19,41 @@ RootManager::RootManager() printf("RootManager has been created.\n"); fflush(stdout); - // open root file int t = (int) time(0); // get time now fRootfile = new TFile(Form("output_at_%d_seconds.root",t ),"RECREATE"); + if(!fRootfile) { cout << "\nCould not open file " <SetAutoSave(100000); - /* + /* At this stage you can define what branches are written in the tree */ + //fOutputTree->Branch("S3Branch","TS3Data",&fS3Data); + //fOutputTree->Branch("SpiceBranch","TSpiceData",&fSpiceData); + //---------------- //fOutputTree->Branch("TTigFragment","TTigFragment",&fFragment, 1000, 99); //---------------- - fOutputTree->Branch("SpiceBranch","TSpiceData",&fSpiceData); - fOutputTree->Branch("S3Branch","TS3Data",&fS3Data); + //fOutputTree->Branch("PacesBranch","TPacesData",&fPacesData); + //---------------- + fOutputTree->Branch("NewBranch", "TNewData", &fNewData); + //---------------- + //fOutputTree->Branch("SceptarBranch","TSceptarData",&fSceptarData); //---------------- // fOutputTree->Branch("GriffinBranch","TGriffinData",&fGriffinData); //---------------- + fOutputTree->Branch("HistoryBranch","THistoryData",&fHistoryData); + //---------------- + /* Other detector branches goes here */ + + } -} - -void RootManager::FillHist(double temp) { - float tempf = (float)temp; - fHist->Fill(tempf); -} - - +void RootManager::FillHist(double temp) { + float tempf = (float)temp; + fHist->Fill(tempf); + } + void RootManager::FillG4Hit(string volume, // Word representing the key used to identify the detector, and to build mnemonics -int detector, -int crystal, -int pdg, // integer representing the type of particle -double Energy, // depositid energy in a step -double Px, double Py, double Pz,// position vector -int Id, // original(primary) ID -int PrimPdg, // primary particle definition PDG encoding -double PrimEnergy, // original(primary) energy -double Mx, double My, double Mz) { // primary particle momentum vector + int detector, + int crystal, + int pdg, // integer representing the type of particle + double Energy, // depositid energy in a step + double Px, double Py, double Pz,// position vector + int Id, // original(primary) ID + int PrimPdg, // primary particle definition PDG encoding + double PrimEnergy, // original(primary) energy + double Mx, double My, double Mz) { // primary particle momentum vector /* Some other functions */ - //Get the menmonic string mnemonic = BuildMnemonic(volume,detector,crystal); //Fill the event, calculate the full deposited energy etc... @@ -108,77 +130,127 @@ double Mx, double My, double Mz) { // primary particle momentum vector -void RootManager::SortEvent(void) { -// Sort Data from the map second element by getters and set them - std::map::iterator it; - for (std::map::iterator it=fGeantEvent.begin(); it!=fGeantEvent.end(); ++it) { - - string system (it->first, 0, 3); // first three letters of the mnemonic (it->first) defines the system +void RootManager::SortEvent(int eventNb) { - //Fragment - //if (1) SetFragmentEvent(it->first); // take all the event in the fragment tree - - //Spice - if (system=="SPI") SetSpiceEvent(it->first, it->second.GetDetector(), it->second.GetCrystal()); - if (system=="SPE") SetS3Event(it->first, it->second.GetDetector(), it->second.GetCrystal()); + // clear the Fragment + //fFragment->ClearVariables(); - //Griffin - //Need to put Griffin key here + // clear the SpiceData object + //fSpiceData->ClearVariables(); + //fS3Data->ClearVariables(); - //other - //if (system==XYZ) SetOtherDetectorEvent(key); - } + // clear the PacesData object + //fPacesData->ClearVariables(); - -// fill the tree by SpiceData - fOutputTree->Fill(); - -// clear the map - fGeantEvent.clear(); + //clear the New Object + //fNewData->ClearVariables(); -// clear the Fragment - //fFragment->Clear(); - -// clear the SpiceData object - fSpiceData->Clear(); - fS3Data->Clear(); + //clear the Sceptar Object + //fSceptarData->ClearVariables(); -// clear the GriffinData object -// fGriffinData->Clear(); -} + // clear the GriffinData object + // fGriffinData->ClearVariables(); -string RootManager::BuildMnemonic(string volume, int detector, int crystal) { + // Sort Data from the map second element by getters and set them + std::map::iterator it; + for (std::map::iterator it=fGeantEvent.begin(); it!=fGeantEvent.end(); ++it) { + + string system (it->first, 0, 3); // first three letters of the mnemonic (it->first) defines the system + // cout << " RootManager::SortEvent -- system " << system << endl ; + //Fragment + //if (1) SetFragmentEvent(it->first); // take all the event in the fragment tree + + //Spice + //if (system=="SPI") SetSpiceEvent(eventNb, it->first, it->second.GetDetector(), it->second.GetCrystal()); + //if (system=="SPE") SetS3Event(eventNb, it->first, it->second.GetDetector(), it->second.GetCrystal()); + + //Paces + //if (system=="PAC") { + // SetPacesEvent(eventNb, it->first, it->second.GetDetector(), it->second.GetCrystal()); + // fPacesData->Dump(); + // } -std::string system (volume, 0, 2); -std::string sub_system (volume, 2, 1); -std::string number = "dummy" ; + //NEW + if (system=="NEW") SetNewEvent(eventNb, it->first, it->second.GetDetector(), it->second.GetCrystal()); -//Build the mnemonic for spice -if (system == "SP") - if (sub_system == "I") { + //Sceptar + //if (system=="SEP") SetSceptarEvent(eventNb, it->first, it->second.GetDetector(), it->second.GetCrystal()); + + //Griffin + //Need to put Griffin key here + + //New Detector + //if (system==XYZ) SetOtherDetectorEvent(key); + } - detector = 9 - detector ; // spice is upstream the target, according to the mnemonics outer ring is '0' inner ring is 9. + +// fill the tree + fOutputTree->Fill(); - ostringstream convert; // stream used for the conversion - convert << std::setw(3) << std::setfill('0') << (detector * 12 + crystal); // set the width to 3 (xyz), fill the blanks by zeros - number = convert.str(); // set 'Result' to the contents of the stream +// ClearVariables the map + fGeantEvent.clear(); - return (volume + "00XN" + number) ; - } - else if (sub_system == "E") { - - //some operations - return volume+"00XN"+ number ; - } +} +string RootManager::BuildMnemonic(string volume, int detector, int crystal) { -//Griffin + std::string system (volume, 0, 2); + std::string sub_system (volume, 2, 1); + std::string number = "00" ; + + //Build the mnemonic for spice + //--------------------- + if (system == "SP"){ // SPICE SiLi + if (sub_system == "I") { + detector = 9 - detector ; // (detector here = ring) spice is upstream the target, according to the mnemonics outer ring is '0' inner ring is 9. + ostringstream convert; // stream used for the conversion + convert << std::setw(3) << std::setfill('0') << (detector * 12 + crystal); // set the width to 3 (xyz), fill the blanks by zeros + number = convert.str(); // set 'Result' to the contents of the stream + return ( system + sub_system + "00XN" + number) ; + } + else if (sub_system == "E") { // SPICE S3 + return system + sub_system +"00"+"XN"+ number ; + } + } + //--------------------- PACES + else if (system == "PA" && sub_system == "C") { + // (detector here = ring), Paces has only one "ring", => detector = 0 + ostringstream convert; // stream used for the conversion + convert << std::setw(2) << std::setfill('0') << (crystal); // set the width to 2, i.e. (xy), fill the blanks by zeros = { 01, ... , 05 } + number = convert.str(); // set 'Result' to the contents of the stream + //cout << " RootManager::BuildMnemonic : " << system + sub_system << number <<"XN00" << " | " << system + sub_system + number+"XN"+"00" << endl ; + return system + sub_system + number + "XN" + "00"; + } + //--------------------- New PACES + else if (system == "NE" && sub_system == "W") { + detector = 4 - detector ; // New is upstream the target + ostringstream convert; // stream used for the conversion + convert << std::setw(3) << std::setfill('0') << (detector * 12 + crystal); // set the width to 3 (xyz), fill the blanks by zeros + number = convert.str(); // set 'Result' to the contents of the stream + return (volume + "00XN" + number) ; + } + //--------------------- Sceptar + else if (system == "SE"){ + if (sub_system == "P") { + ostringstream convert; // stream used for the conversion + convert << std::setw(3) << std::setfill('0') << (crystal ); // set the width to 3 (xyz), fill the blanks by zeros + number = convert.str(); // set 'Result' to the contents of the stream + return (volume + "00XN" + number) ; + } + } + //--------------------- + //Griffin -//else + //--------------------- + //Other detectors + + //------------------ + //Default + return volume; } @@ -186,22 +258,20 @@ if (system == "SP") void RootManager::SetFragmentEvent(string mnemonic) { - // treat - fGeantEvent.at(mnemonic).SortPrimary(); - - // get the Energy - double energy = fGeantEvent.at(mnemonic).GetFullEnergy(); - - // Set the data into the fragment , NB : TTigFragment Class has the all the memebers "public", no setters are used, instead we could assigne the values immediatly. - fFragment->ChannelName = mnemonic ; - fFragment->ChargeCal = energy ; -} + // treat + fGeantEvent.at(mnemonic).SortPrimary(); + // get the Energy + double energy = fGeantEvent.at(mnemonic).GetFullEnergy(); + // Set the data into the fragment , NB : TTigFragment Class has the all the memebers "public", no setters are used, instead we could assigne the values immediatly. + fFragment->ChannelName = mnemonic ; + fFragment->ChargeCal = energy ; + } -void RootManager::SetSpiceEvent(string mnemonic, int Ring, int Seg) -{ +void RootManager::SetSpiceEvent(int eventNb, string mnemonic, int Ring, int Seg) { // treat fGeantEvent.at(mnemonic).SortPrimary(); + fSpiceData->SetEventNumber(eventNb) ; // get primary // Pdg int mult = fGeantEvent.at(mnemonic).GetPrimaryPdgMult(); // inside this particular pad @@ -231,24 +301,25 @@ void RootManager::SetSpiceEvent(string mnemonic, int Ring, int Seg) // (Th,E) fSpiceData->SetSpiceThetaEDetectorNbr(1) ; fSpiceData->SetSpiceThetaEStripNbr(Ring) ; - fSpiceData->SetSpiceThetaEEnergy(energy) ; + fSpiceData->SetSpiceThetaEEnergy(energy/keV) ; fSpiceData->SetSpiceThetaEResEnergy(applied_resolution) ; // (Ph,E) fSpiceData->SetSpicePhiEDetectorNbr(1) ; fSpiceData->SetSpicePhiEStripNbr(Seg) ; - fSpiceData->SetSpicePhiEEnergy(energy) ; + fSpiceData->SetSpicePhiEEnergy(energy/keV) ; fSpiceData->SetSpicePhiEResEnergy(applied_resolution) ; fSpiceData->SetPositionFirstHit(pos) ; } -void RootManager::SetS3Event(string mnemonic, int Ring, int Seg) { +void RootManager::SetS3Event(int eventNb, string mnemonic, int Ring, int Seg) { // treat fGeantEvent.at(mnemonic).SortPrimary(); + fS3Data->SetEventNumber(eventNb) ; // get primary // Pdg int mult = fGeantEvent.at(mnemonic).GetPrimaryPdgMult(); // inside this particular pad @@ -257,53 +328,297 @@ void RootManager::SetS3Event(string mnemonic, int Ring, int Seg) { // Energy mult = fGeantEvent.at(mnemonic).GetPrimaryEnergyMult(); // this should be the same as above - for (int i = 0 ; iSetPrimaryEnergy( fGeantEvent.at(mnemonic).GetPrimaryEnergy(i) ) ; } // Momentum mult = fGeantEvent.at(mnemonic).GetPrimaryThetaMult(); // this should be the same as above - for (int i = 0 ; iSetPrimaryTheta(fGeantEvent.at(mnemonic).GetPrimaryTheta(i) ) ; fS3Data->SetPrimaryPhi( fGeantEvent.at(mnemonic).GetPrimaryPhi(i) ) ; } // get the energy double energy = fGeantEvent.at(mnemonic).GetFullEnergy(); + double stDev = (SpiceResolution[1]*keV) * energy + (SpiceResolution[0]*keV); // the result is in MeV + double applied_resolution = CLHEP::RandGauss::shoot(energy, stDev); TVector3 pos = fGeantEvent.at(mnemonic).GetSecondHitPosition() ; // fill the S3Data object // (Th,E) fS3Data->SetS3ThetaEDetectorNbr(1) ; - fS3Data->SetS3ThetaEStripNbr(Ring) ; - fS3Data->SetS3ThetaEEnergy(energy) ; - + fS3Data->SetS3ThetaEStripNbr(Ring) ; + fS3Data->SetS3ThetaEEnergy(energy/keV) ; + // (Ph,E) fS3Data->SetS3PhiEDetectorNbr(1) ; fS3Data->SetS3PhiEStripNbr(Seg) ; - fS3Data->SetS3PhiEEnergy(energy) ; + fS3Data->SetS3PhiEEnergy(applied_resolution/keV) ; fS3Data->SetPositionFirstHit(pos) ; } -void RootManager::SetGriffinEvent(int Crystal) + +void RootManager::SetPacesEvent(int eventNb, string mnemonic, int Ring, int Seg) { - //Stuff goes here - cout << " " << endl; + + // treat + fGeantEvent.at(mnemonic).SortPrimary(); + + //cout << " mnemonic " << mnemonic << endl; + //cout << " Seg " << Seg << endl; + //cout << " Ring " << Ring << endl; + + //fGeantEvent.at(mnemonic).ShowVectorsContent(); + //G4cin.get(); + + fPacesData->SetEventNumber(eventNb) ; + // get primary + // Pdg + int mult = fGeantEvent.at(mnemonic).GetPrimaryPdgMult(); // inside this particular pad + for (int i = 0 ; iSetPrimaryPdg( fGeantEvent.at(mnemonic).GetPrimaryPdg(i) ) ; + + // Energy + mult = fGeantEvent.at(mnemonic).GetPrimaryEnergyMult(); // this should be the same as above + for (int i = 0 ; iSetPrimaryEnergy( fGeantEvent.at(mnemonic).GetPrimaryEnergy(i) ) ; + } + + // Momentum + mult = fGeantEvent.at(mnemonic).GetPrimaryThetaMult(); // this should be the same as above + for (int i = 0 ; iSetPrimaryTheta( fGeantEvent.at(mnemonic).GetPrimaryTheta(i) ) ; + fPacesData->SetPrimaryPhi( fGeantEvent.at(mnemonic).GetPrimaryPhi(i) ) ; + } + + // get the energy + double energy = fGeantEvent.at(mnemonic).GetFullEnergy(); + TVector3 pos = fGeantEvent.at(mnemonic).GetFirstHitPosition() ; + + // fill the PacesData object + Ring = Ring ; // irrelevant, in PACES we only have one "ring", should take off this parameter at some point. + fPacesData->SetPacesDetectorNbr(Seg) ; + fPacesData->SetPacesEnergy(energy/keV) ; + + //pos.Dump(); + //getchar(); + fPacesData->SetPositionFirstHit(pos) ; + +} + + +void RootManager::SetNewEvent(int eventNb, string mnemonic, int detector, int Seg) +{ + + //treat + fGeantEvent.at(mnemonic).SortPrimary(); + + fNewData->SetEventNumber(eventNb) ; + //get primary + //Pdg + int mult = fGeantEvent.at(mnemonic).GetPrimaryPdgMult(); + for( int i = 0 ; iSetPrimaryPdg(fGeantEvent.at(mnemonic).GetPrimaryPdg(i) ) ; + } + + //Energy + mult = fGeantEvent.at(mnemonic).GetPrimaryEnergyMult(); + for(int i = 0 ; iSetPrimaryEnergy(fGeantEvent.at(mnemonic).GetPrimaryEnergy(i) ); + } + + //Momentum + mult = fGeantEvent.at(mnemonic).GetPrimaryThetaMult(); + for(int i = 0 ; iSetPrimaryTheta(fGeantEvent.at(mnemonic).GetPrimaryTheta(i) ) ; + fNewData->SetPrimaryPhi(fGeantEvent.at(mnemonic).GetPrimaryPhi(i) ); + } + + //get the energy + double energy = fGeantEvent.at(mnemonic).GetFullEnergy(); + TVector3 pos = fGeantEvent.at(mnemonic).GetFirstHitPosition(); + + //fill the NewData Object + fNewData->SetNewDetectorNbr(detector); + fNewData->SetNewEnergy(energy); + fNewData->SetPositionFirstHit(pos); + + //fNewData->Dump(); } -// set event number -void RootManager::SetEventNumber(int i ) +void RootManager::SetSceptarEvent(int eventNb, string mnemonic, int detector, int Seg) { -fSpiceData->SetEventNumber(i) ; -fS3Data->SetEventNumber(i) ; -//fGriffinData->SetEventNumber(i); -} + + //treat + fGeantEvent.at(mnemonic).SortPrimary(); + + fSceptarData->SetEventNumber(eventNb) ; + //get primary + //Pdg + int mult = fGeantEvent.at(mnemonic).GetPrimaryPdgMult(); + for( int i = 0 ; iSetPrimaryPdg(fGeantEvent.at(mnemonic).GetPrimaryPdg(i) ) ; + } + + //Energy + mult = fGeantEvent.at(mnemonic).GetPrimaryEnergyMult(); + for(int i = 0 ; iSetPrimaryEnergy(fGeantEvent.at(mnemonic).GetPrimaryEnergy(i) ); + } + + //Momentum + mult = fGeantEvent.at(mnemonic).GetPrimaryThetaMult(); + for(int i = 0 ; iSetPrimaryTheta(fGeantEvent.at(mnemonic).GetPrimaryTheta(i) ) ; + fSceptarData->SetPrimaryPhi(fGeantEvent.at(mnemonic).GetPrimaryPhi(i) ); + } + + //get the energy + double energy = fGeantEvent.at(mnemonic).GetFullEnergy(); + TVector3 pos = fGeantEvent.at(mnemonic).GetFirstHitPosition(); + + //fill the SceptarData Object + fSceptarData->SetSceptarDetectorNbr(detector); + fSceptarData->SetSceptarEnergy(energy); + fSceptarData->SetPositionFirstHit(pos); + + //fSceptarData->Dump(); + +} + + void RootManager::ClearVariables(void){ + + //printf("RootManager : ClearVariables .\n"); + + // clear the Fragment + //fFragment->ClearVariables(); + + // clear the HistoryData object + fHistoryData->ClearVariables(); + + // ClearVariables the SpiceData object + //fSpiceData->ClearVariables(); + //fS3Data->ClearVariables(); + + // ClearVariables the PacesData object + //fPacesData->ClearVariables(); + + // clear the NewPaces + fNewData->ClearVariables(); + + // clear Sceptar + //fSceptarData->ClearVariables(); + + // ClearVariables the GriffinData object + //fGriffinData->ClearVariables(); + } + +void RootManager::SetHistory( vector info ){ + + for (unsigned iInfo = 0 ; iInfo < info.size() ; iInfo ++ ) { + + //info.at(iInfo)->Print(); + double x, y, z ;// dummy variabales + + //primary information + fHistoryData->SetHistoryPrimaryID(info.at(iInfo)->GetOriginalTrackID()) ; + fHistoryData->SetHistoryPrimaryPdg(info.at(iInfo)->GetOriginalPdg()) ; + fHistoryData->SetHistoryPrimaryEnergy(info.at(iInfo)->GetOriginalEnergy()) ; + + x = info.at(iInfo)->GetOriginalPosition().getX(); + y = info.at(iInfo)->GetOriginalPosition().getY(); + z = info.at(iInfo)->GetOriginalPosition().getZ(); + fHistoryData->SetHistoryPrimaryPositionVertex(x,y,z) ; + + x = info.at(iInfo)->GetOriginalMomentum().getX(); + y = info.at(iInfo)->GetOriginalMomentum().getY(); + z = info.at(iInfo)->GetOriginalMomentum().getZ(); + fHistoryData->SetHistoryPrimaryMomentumVertex(x,y,z) ; + + // 1s timpact of the primary + x = info.at(iInfo)->GetOriginalImpactMomentum().getX(); + y = info.at(iInfo)->GetOriginalImpactMomentum().getY(); + z = info.at(iInfo)->GetOriginalImpactMomentum().getZ(); + fHistoryData->SetHistoryPrimaryMomentum1stImpact(x,y,z) ; + + x = info.at(iInfo)->GetOriginalImpactPosition().getX(); + y = info.at(iInfo)->GetOriginalImpactPosition().getY(); + z = info.at(iInfo)->GetOriginalImpactPosition().getZ(); + fHistoryData->SetHistoryPrimaryPosition1stImpact(x,y,z) ; + fHistoryData->SetHistoryPrimaryVolume1stImpact(info.at(iInfo)->GetOriginalImpactVolume()) ; + + // Parent of the current track hitting the target + fHistoryData->SetHistoryParentID(info.at(iInfo)->GetCurrentParentID()) ; + + // the current track hitting the target information + fHistoryData->SetHistoryCurrentID(info.at(iInfo)->GetCurrentTrackID()) ; + fHistoryData->SetHistoryCurrentPdg(info.at(iInfo)->GetCurrentPdg()) ; + fHistoryData->SetHistoryCurrentEnergy(info.at(iInfo)->GetCurrentEnergyAtVertex()) ; + fHistoryData->SetHistoryCurrentCreatorProcess(info.at(iInfo)->GetCurrentProcess()) ; + fHistoryData->SetHistoryCurrentTime(info.at(iInfo)->GetCurrentTimeAtVertex()) ; + fHistoryData->SetHistoryCurrentVolume1stImpact(info.at(iInfo)->GetCurrentImpactVolume()) ; + + + x = info.at(iInfo)->GetCurrentPositionAtVertex().getX(); + y = info.at(iInfo)->GetCurrentPositionAtVertex().getY(); + z = info.at(iInfo)->GetCurrentPositionAtVertex().getZ(); + fHistoryData->SetHistoryCurrentPositionVertex(x,y,z) ; + + x = info.at(iInfo)->GetCurrentImpactPosition().getX(); + y = info.at(iInfo)->GetCurrentImpactPosition().getY(); + z = info.at(iInfo)->GetCurrentImpactPosition().getZ(); + fHistoryData->SetHistoryCurrentPosition1stImpact(x,y,z) ; + + x = info.at(iInfo)->GetCurrentPositionAtDetector().getX(); + y = info.at(iInfo)->GetCurrentPositionAtDetector().getY(); + z = info.at(iInfo)->GetCurrentPositionAtDetector().getZ(); + fHistoryData->SetHistoryCurrentPositionDetector(x,y,z) ; + + x = info.at(iInfo)->GetCurrentPositionAtDeath().getX(); + y = info.at(iInfo)->GetCurrentPositionAtDeath().getY(); + z = info.at(iInfo)->GetCurrentPositionAtDeath().getZ(); + fHistoryData->SetHistoryCurrentPositionDeath(x,y,z) ; + + + + x = info.at(iInfo)->GetCurrentMomentumAtVertex().getX(); + y = info.at(iInfo)->GetCurrentMomentumAtVertex().getY(); + z = info.at(iInfo)->GetCurrentMomentumAtVertex().getZ(); + fHistoryData->SetHistoryCurrentMomentumVertex(x,y,z) ; + + x = info.at(iInfo)->GetCurrentImpactMomentum().getX(); + y = info.at(iInfo)->GetCurrentImpactMomentum().getY(); + z = info.at(iInfo)->GetCurrentImpactMomentum().getZ(); + fHistoryData->SetHistoryCurrentMomentum1stImpact(x,y,z) ; + + x = info.at(iInfo)->GetCurrentMomentumAtDetector().getX(); + y = info.at(iInfo)->GetCurrentMomentumAtDetector().getY(); + z = info.at(iInfo)->GetCurrentMomentumAtDetector().getZ(); + fHistoryData->SetHistoryCurrentMomentumDetector(x,y,z) ; + + x = info.at(iInfo)->GetCurrentMomentumAtDeath().getX(); + y = info.at(iInfo)->GetCurrentMomentumAtDeath().getY(); + z = info.at(iInfo)->GetCurrentMomentumAtDeath().getZ(); + fHistoryData->SetHistoryCurrentMomentumDeath(x,y,z) ; + + } + + //cin.get(); +} + +void RootManager::SetGriffinEvent(int Crystal) +{ + Crystal = Crystal; + //Stuff goes here + //cout << " " << endl; + +} + void RootManager::Close() @@ -315,3 +630,5 @@ void RootManager::Close() + + diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc index a67a4db..505ae10 100755 --- a/src/SteppingAction.cc +++ b/src/SteppingAction.cc @@ -81,30 +81,23 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) { G4int particleType = 0; - G4int volNameOver9; - G4int evntNb; - + G4int evntNb = 0 ; det = 0; cry = 0; - stepNumber++; - // Get volume of the current step G4VPhysicalVolume* volume = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(); G4String volname = volume->GetName(); - - // Get the process of the current step (PreStep or Poststep??) - //G4String process = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); - + // collect energy and track length step by step // As it's called more than once, get the Track and assign to variable G4double edep = aStep->GetTotalEnergyDeposit(); G4double ekin = aStep->GetPreStepPoint()->GetKineticEnergy(); G4Track* theTrack = aStep->GetTrack(); - G4double stepl = 0.; - if (theTrack->GetDefinition()->GetPDGCharge() != 0.) - stepl = aStep->GetStepLength(); + G4double stepl = 0.; + if (theTrack->GetDefinition()->GetPDGCharge() != 0.) + stepl = aStep->GetStepLength(); // Track particle type in EVERY step //G4cout << "Particle name = " << aStep->GetTrack()->GetParticleDefinition()->GetParticleName() << G4endl; @@ -117,45 +110,129 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) eventaction->AddParticleType(particleType); evntNb = eventaction->GetEventNumber(); - //G4cout << "Found Edep = " << edep/keV << " keV in " << volname << G4endl; - // example volname - //volname = av_1_impr_6_sodium_iodide_crystal_block_log_pv_0 - - // Get initial momentum direction & energy of particle - G4int trackID = theTrack->GetTrackID(); - G4int parentID = theTrack->GetParentID(); - - // The vertex corresponds to where the particle has been created anywhere in the chamber - //G4double initialDirectionX = theTrack->GetVertexMomentumDirection().getX(); - //G4double initialDirectionY = theTrack->GetVertexMomentumDirection().getY(); - //G4double initialDirectionZ = theTrack->GetVertexMomentumDirection().getZ(); - //G4double initialEnergy = theTrack->GetVertexKineticEnergy(); - // if (parentID == 0) initialEnergy = theTrack->GetVertexKineticEnergy(); - - + // Get the track (extra) information , mainly about the primary particle, but could be used for other stuff TrackInformation* info = (TrackInformation*)(aStep->GetTrack()->GetUserInformation()); - //info->Print(); // for inspection - // The Origin corresponds to the information about the primary particle (exclusively from the source) + //info->SetTagged(false); + + +// ------------- add information to the track info -------------------- + +// After the very first step + if (aStep->GetTrack()->GetCurrentStepNumber()== 1 ){ + + if ( aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()==0) { // if out of world + + if (aStep->GetTrack()->GetCreatorProcess() == 0 ) { // if original particle (this information will be past by to all the descendants) + info->SetTagged(true); + info->SetOriginalImpactVolume("OOW");// Out Of World + info->SetOriginalImpactPosition(aStep->GetTrack()->GetPosition()) ; + info->SetOriginalImpactMomentum(aStep->GetTrack()->GetMomentum()) ; + } + else { // if other particle (this information will change with every new track) + info->SetCurrentImpactVolume("OOW"); + info->SetCurrentImpactPosition(aStep->GetTrack()->GetPosition()) ; + info->SetCurrentImpactMomentum(aStep->GetTrack()->GetMomentum()) ; + //cin.get() ; + } + } + else { // if NOT out of world + // if original particle (this information will be past by to all the descendants) + if (aStep->GetTrack()->GetCreatorProcess() == 0 ) { + info->SetOriginalImpactVolume(aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName()); + info->SetOriginalImpactPosition(aStep->GetTrack()->GetPosition()) ; + info->SetOriginalImpactMomentum(aStep->GetTrack()->GetMomentum()) ; + } + else {// if other particle (this information will change with every new track) + info->SetCurrentImpactVolume(aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName()); + info->SetCurrentImpactPosition(aStep->GetTrack()->GetPosition()) ; + info->SetCurrentImpactMomentum(aStep->GetTrack()->GetMomentum()) ; + } + } + + } + + +//Tag this particle for any event + info->SetTagged(true); + +// Tag this track if the particle will hit the detector +// Get volume before and after for each step + G4String after = "VolAfter" ; + G4String before = "VolBefore" ; + if ( aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()!=0) { + before = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() ; + after = aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName() ; + } + else { + before = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() ; + after = "OOW" ; + } + +// If the detector is touched get position and momentum (and time?) + + // Build the selection + G4bool selection = 0 ; + selection += (before.contains("SquareDetect") || after.contains("SquareDetect")) && !info->GetTagged(); + selection += (before.contains("siDetS3") || after.contains("siDetS3")) && !info->GetTagged() ; + + + if( selection ){ + // Get position and momentum + info->SetCurrentPositionAtDetector(aStep->GetTrack()->GetPosition()) ; + info->SetCurrentMomentumAtDetector(aStep->GetTrack()->GetMomentum()) ; + //info->SetTagged(true); + } + +// this is the last step, in case the particle is tagged => add some new info and keep the track + if (aStep->GetTrack()->GetTrackStatus() == 2 && info->GetTagged() ) { + + // Get position + info->SetCurrentPositionAtDeath(aStep->GetTrack()->GetPosition()) ; + info->SetCurrentMomentumAtDeath(aStep->GetTrack()->GetMomentum()) ; + // Get Pdg + info->SetCurrentPdg(aStep->GetTrack()->GetDefinition()->GetPDGEncoding()); + // Get Energy and time of the track at vertex + info->SetCurrentEnergyAtVertex(aStep->GetTrack()->GetVertexKineticEnergy()); + info->SetCurrentTimeAtVertex(aStep->GetTrack()->GetGlobalTime()); // in ns + // Get Vertex of the track + info->SetCurrentPositionAtVertex(aStep->GetTrack()->GetVertexPosition()); + // Get momentum at the vertex + info->SetCurrentMomentumAtVertex(aStep->GetTrack()->GetVertexMomentumDirection()); + // Get the creator process + if(aStep->GetTrack()->GetCreatorProcess()!=0) + info->SetCurrentProcess(aStep->GetTrack()->GetCreatorProcess()->GetProcessName()); + else + info->SetCurrentProcess("source") ; + + // Keep the track in the record + //cout << " <<<<<< Setting info >>>>>> " << endl ; + eventaction->SetPrimaryInfo( info ) ; + //cout << " <<<<<< <<<<<< >>>>>> >>>>>> " << endl ; + + //clear the parts unrelated to the original particle information, and untag for the new track + info->PartialClear(); // info->SetTagged(false); included + } + + G4int OriginID = info->GetOriginalTrackID() ; G4int OriginPdg = info->GetOriginalPdg() ; - G4double OriginEnergy = info->GetOriginalEnergy() ; // Kinetic Energy + G4double OriginEnergy = info->GetOriginalEnergy() ; // Kinetic Energy G4ThreeVector OriginMoment = info->GetOriginalMomentum() ; - + G4StepPoint* point1 = aStep->GetPreStepPoint(); G4StepPoint* point2 = aStep->GetPostStepPoint(); G4ThreeVector pos1 = point1->GetPosition(); G4ThreeVector pos2 = point2->GetPosition(); - G4double time1 = point1->GetGlobalTime(); + //G4double time1 = point1->GetGlobalTime(); G4double time2 = point2->GetGlobalTime(); - - + G4int trackID = theTrack->GetTrackID(); + size_t found; G4String search; - G4int searchLength; // Grid Cell found = volname.find("gridcell"); @@ -241,21 +318,42 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) // Sceptar found = volname.find("sceptar_square_scintillator_log"); if (edep != 0 && found!=G4String::npos) { - SetDetNumberForGenericDetector(volname); - eventaction->AddSceptarSquareCrystDet(edep,stepl,det-1); + SetDetAndCryNumberForSceptarDetector(volname); + eventaction->AddSceptarSquareCrystDet(edep,stepl,det); + eventaction->AddStepTracker(evntNb, stepNumber, "SEP", cry, det, edep, pos2.x(), pos2.y(), pos2.z(), time2, OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), OriginEnergy, OriginPdg, OriginID, trackID); } found = volname.find("sceptar_angled_scintillator_log"); if (edep != 0 && found!=G4String::npos) { - SetDetNumberForGenericDetector(volname); - eventaction->AddSceptarAngledCrystDet(edep,stepl,det-1); + SetDetAndCryNumberForSceptarDetector(volname); + eventaction->AddSceptarAngledCrystDet(edep,stepl,det); + eventaction->AddStepTracker(evntNb, stepNumber, "SEP", cry, det, edep, pos2.x(), pos2.y(), pos2.z(), time2, OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), OriginEnergy, OriginPdg, OriginID, trackID); } // Paces found = volname.find("paces_silicon_block_log"); if (edep != 0 && found!=G4String::npos) { - SetDetNumberForGenericDetector(volname); - eventaction->AddPacesCrystDet(edep,stepl,det-1); + SetDetAndCryNumberForPacesDetector(volname); + eventaction->AddPacesCrystDet(edep,stepl,det); + eventaction->AddStepTracker(evntNb, stepNumber, "PAC", cry, det, edep, + pos2.x(), pos2.y(), pos2.z(), time2, + OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), + OriginEnergy, OriginPdg, OriginID, trackID); + } + + // New Paces + found = volname.find("SquareDetect"); + if(edep !=0 && found!=G4String::npos){ + SetDetAndCryNumberForNewDetector(volname); + eventaction->AddNewCrystDet(edep,stepl,det); + eventaction->AddStepTracker(evntNb, stepNumber, "NEW", cry, det, edep, + pos2.x(), pos2.y(), pos2.z(), time2, + OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), + OriginEnergy, OriginPdg, OriginID, trackID); + +//cout << edep << endl ; +//cin.get() ; + } //SPICE @@ -263,7 +361,10 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) if (edep != 0 && found!=G4String::npos) { SetDetAndCryNumberForSpiceDetector(volname); eventaction->AddSpiceCrystDet(edep,stepl,det); - eventaction->AddStepTracker(evntNb, stepNumber, "SPI", cry, det, edep, pos2.x(), pos2.y(), pos2.z(), time2, OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), OriginEnergy, OriginPdg, OriginID, trackID); + eventaction->AddStepTracker(evntNb, stepNumber, "SPI", cry, det, edep, + pos2.x(), pos2.y(), pos2.z(), time2, + OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), + OriginEnergy, OriginPdg, OriginID, trackID); } //S3 of SPICE @@ -271,9 +372,16 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) if (edep != 0 && found!=G4String::npos) { SetDetAndCryNumberForS3Detector(volname); eventaction->AddSpiceCrystDet(edep,stepl,det); - eventaction->AddStepTracker(evntNb, stepNumber, "SPE", cry, det, edep, pos2.x(), pos2.y(), pos2.z(), time2, OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), OriginEnergy, OriginPdg, OriginID, trackID); + eventaction->AddStepTracker(evntNb, stepNumber, "SPE", cry, det, edep, + pos2.x(), pos2.y(), pos2.z(), time2, + OriginMoment.getX(), OriginMoment.getY(), OriginMoment.getZ(), + OriginEnergy, OriginPdg, OriginID, trackID); } - + + + stepNumber++; + + } void SteppingAction::SetDetAndCryNumberForGriffinComponent(G4String volname) @@ -349,7 +457,6 @@ void SteppingAction::SetDetAndCryNumberForSpiceDetector(G4String volname) dummy = volname.substr (UnderScoreIndex[4]+1,UnderScoreIndex[5]-UnderScoreIndex[4]-1); det = atoi(dummy.c_str()); // ring - //G4cout << " (Stepping action) in " << volname << " segment = " << cry << " ring = " << det << G4endl; //G4cout << " in " << volname << " segment = " << cry << " ring = " << det << G4endl; //G4cin.get(); @@ -377,6 +484,84 @@ void SteppingAction::SetDetAndCryNumberForS3Detector(G4String volname) } +void SteppingAction::SetDetAndCryNumberForPacesDetector(G4String volname) +{ + // the volume name contains five underscrores : av_xxx_impr_SegmentID_siDetSpiceRing_RingID_etc... + G4String dummy=""; + size_t UnderScoreIndex[6]; + size_t old = -1 ; + for (int i = 0 ; i < 6 ; i++ ){ + UnderScoreIndex[i] = volname.find_first_of("_",old+1); + old = UnderScoreIndex[i] ; + } + + dummy = volname.substr (UnderScoreIndex[2]+1,UnderScoreIndex[3]-UnderScoreIndex[2]-1); // select the substring between the underscores + cry = atoi(dummy.c_str()) ; // in paces, start counting from 1 + + dummy = volname.substr (UnderScoreIndex[4]+1,UnderScoreIndex[5]-UnderScoreIndex[4]-1); + det = atoi(dummy.c_str()); // ring + + //G4cout << " (Stepping action) in " << volname << " segment = " << cry << " ring = " << det << G4endl; + //G4cin.get(); +} + +void SteppingAction::SetDetAndCryNumberForNewDetector(G4String volname) +{ + //the volume name contains eight underscores + G4String dummy=""; + size_t UnderScoreIndex[8]; + size_t old = -1; + for (int i=0; i<8 ; i++){ + UnderScoreIndex[i] = volname.find_first_of("_",old+1); + old=UnderScoreIndex[i]; + } + + +// cout << volname << endl ; +//cin.get() ; + + dummy = volname.substr (UnderScoreIndex[7]+1, UnderScoreIndex[8]-UnderScoreIndex[7]-1); + cry = atoi(dummy.c_str()) ; + + dummy = volname.substr (UnderScoreIndex[4]+1, UnderScoreIndex[5]-UnderScoreIndex[4]-1); + det = atoi(dummy.c_str()) ; + + +// cout << cry << " " << det << endl ; +//cin.get() ; + + +} + +void SteppingAction::SetDetAndCryNumberForSceptarDetector(G4String volname) +{ + //the volume name contains eight underscores + G4String dummy=""; + size_t UnderScoreIndex[9]; + size_t old = -1; + for (int i=0; i<9 ; i++){ + UnderScoreIndex[i] = volname.find_first_of("_",old+1); + old=UnderScoreIndex[i]; + } + + +// cout << volname << endl ; +//cin.get() ; + + dummy = volname.substr (UnderScoreIndex[2]+1, UnderScoreIndex[3]-UnderScoreIndex[2]-1); + cry = atoi(dummy.c_str()) ; + + dummy = volname.substr (UnderScoreIndex[0]+1, UnderScoreIndex[1]-UnderScoreIndex[0]-1); + det = atoi(dummy.c_str()) ; + + + //cout << cry << " " << det << endl ; +//cin.get() ; + + +} + + void SteppingAction::SetDetNumberForGenericDetector(G4String volname) { const char *cstr = volname.c_str(); diff --git a/src/TabulatedMagneticField.cc b/src/TabulatedMagneticField.cc index 5aa4857..7918f8b 100644 --- a/src/TabulatedMagneticField.cc +++ b/src/TabulatedMagneticField.cc @@ -1,9 +1,10 @@ #include "TabulatedMagneticField.hh" #include "G4SystemOfUnits.hh" #include "G4PhysicalConstants.hh" +#include "G4ThreeVector.hh" -TabulatedMagneticField::TabulatedMagneticField( const char* filename, G4double zOffset ) - :fZoffset(zOffset),invertX(false),invertY(false),invertZ(false) +TabulatedMagneticField::TabulatedMagneticField( const char* filename, G4double zOffset, G4double zRotation ) + :fZoffset(zOffset),fZrotation(zRotation),invertX(false),invertY(false),invertZ(false) { double lenUnit= mm; @@ -22,6 +23,12 @@ TabulatedMagneticField::TabulatedMagneticField( const char* filename, G4double z char buffer[256]; file.getline(buffer,256); + if (!file.is_open()) { + G4cout <<"\n\ncannot open file : " << filename <<"\n\n" ; + G4cin.get(); + } + + // Read table dimensions file >> nx >> ny >> nz; // Note dodgy order @@ -75,6 +82,7 @@ TabulatedMagneticField::TabulatedMagneticField( const char* filename, G4double z } file.close(); + maxx = xval * lenUnit; maxy = yval * lenUnit; maxz = zval * lenUnit; @@ -114,12 +122,28 @@ void TabulatedMagneticField::GetFieldValue(const double point[4], double x = point[0]; double y = point[1]; double z = point[2] + fZoffset; + + //Rotation treatment Mhd : 25 Mar 2015 + /* + x,y,z ---Rot (angle) ---> x',y',z' + | + v + Bx,By,Bz <---Rot (-angle)--- Bx',By',Bz' + */ + + // Rotate the position here : Mhd : 25 Mar 2015 + G4ThreeVector R( x, y, z); + R.rotateZ(-fZrotation*deg); // rotation made in the opposite direction of the lens + x = R.getX(); + y = R.getY(); + z = R.getZ(); + // Check that the point is within the defined region if ( x>=minx && x<=maxx && y>=miny && y<=maxy && z>=minz && z<=maxz ) { - + // Position of given point within region, normalized to the range // [0,1] double xfraction = (x - minx) / dx; @@ -186,11 +210,19 @@ void TabulatedMagneticField::GetFieldValue(const double point[4], zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal + zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) + zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ; + + // Rotate the BField here : Mhd : 25 Mar 2015 + G4ThreeVector B(Bfield[0], Bfield[1], Bfield[2]); + B.rotateZ(fZrotation*deg); // rotation made in the same direction of the lens + Bfield[0] = B.getX(); + Bfield[1] = B.getY(); + Bfield[2] = B.getZ(); } else { Bfield[0] = 0.0; Bfield[1] = 0.0; Bfield[2] = 0.0; } + } diff --git a/src/TrackInformation.cc b/src/TrackInformation.cc index 2ea6f3c..ec47396 100644 --- a/src/TrackInformation.cc +++ b/src/TrackInformation.cc @@ -6,49 +6,112 @@ G4Allocator aTrackInformationAllocator; TrackInformation::TrackInformation() : G4VUserTrackInformation() { - originalTrackID = 0; - originalPdg = 0; - originalPosition = G4ThreeVector(0.,0.,0.); - originalMomentum = G4ThreeVector(0.,0.,0.); - originalEnergy = 0.; - originalTime = 0.; + +Clear() ; + +} + + +void TrackInformation::Clear() +{ + Tagged = false; + + originalTrackID = -1; + originalPdg = -1; + originalPosition = G4ThreeVector(-1.,-1.,-1.); + originalMomentum = G4ThreeVector(-1.,-1.,-1.); + originalEnergy = -1; + originalTime = -1; + + originalImpactVolume = "None"; + originalImpactPosition = G4ThreeVector(-1.,-1.,-1.); + originalImpactMomentum = G4ThreeVector(-1.,-1.,-1.); + + currentProcess = "None"; + currentTrackID = -1; + currentparentTrackID = -1 ; + currentPdg = -1 ; + + currentPositionAtVertex= G4ThreeVector(-1.,-1.,-1.); + currentPositionAtDetector= G4ThreeVector(-1.,-1.,-1.); + currentPositionAtDeath= G4ThreeVector(-1.,-1.,-1.); // death or out of world + + currentMomentumAtVertex= G4ThreeVector(-1.,-1.,-1.); // contains the angle at the emission + currentMomentumAtDetector= G4ThreeVector(-1.,-1.,-1.); // contains angle at the impact on the detector + currentMomentumAtDeath= G4ThreeVector(-1.,-1.,-1.); // if != 0, then the particle has escaped + + currentEnergyAtVertex = -1; + currentTimeAtVertex = -1; + + currentImpactVolume = "None"; + currentImpactPosition = G4ThreeVector(-1.,-1.,-1.); + currentImpactMomentum = G4ThreeVector(-1.,-1.,-1.); + + SecondariesProcess.clear(); + SecondariesBirthVolume.clear(); + SecondariesDeathVolume.clear(); + SecondariesPdg.clear(); + +} + +void TrackInformation::PartialClear() +{ + Tagged = false; + + currentProcess = -1; + currentTrackID = -1; + currentparentTrackID = -1 ; + currentPdg = -1 ; + + currentPositionAtVertex= G4ThreeVector(-1.,-1.,-1.); + currentPositionAtDetector= G4ThreeVector(-1.,-1.,-1.); + currentPositionAtDeath= G4ThreeVector(-1.,-1.,-1.); // death or out of world + + currentMomentumAtVertex= G4ThreeVector(-1.,-1.,-1.); // contains the angle at the emission + currentMomentumAtDetector= G4ThreeVector(-1.,-1.,-1.); // contains angle at the impact on the detector + currentMomentumAtDeath= G4ThreeVector(-1.,-1.,-1.); // death or out of world - AncestorsBirthVolume.clear(); - AncestorsDeathVolume.clear(); - AncestorsPdg.clear(); + currentEnergyAtVertex = -1; + currentTimeAtVertex = -1; + + SecondariesProcess.clear(); + SecondariesBirthVolume.clear(); + SecondariesDeathVolume.clear(); + SecondariesPdg.clear(); + } + TrackInformation::TrackInformation(const G4Track* aTrack) : G4VUserTrackInformation() { //G4cout<< " this is a track with a parent ID = " <GetParentID()<GetParentID(); originalTrackID = aTrack->GetTrackID(); originalPdg = aTrack->GetDefinition()->GetPDGEncoding(); - originalPosition = aTrack->GetPosition(); - originalMomentum = aTrack->GetMomentum(); - originalEnergy = aTrack->GetKineticEnergy(); + originalPosition = aTrack->GetVertexPosition(); + originalMomentum = aTrack->GetVertexMomentumDirection(); + originalEnergy = aTrack->GetVertexKineticEnergy(); originalTime = aTrack->GetGlobalTime(); - // MHD : 2 May 2013 ; - // Not sure if I need to add something here.. + currentTrackID = aTrack->GetTrackID(); + currentparentTrackID = aTrack->GetParentID(); + + // Process name at birth + if (aTrack->GetCreatorProcess() == 0){ + SecondariesProcess.push_back("Source") ; + } + else { + SecondariesProcess.push_back(aTrack->GetCreatorProcess()->GetProcessName()) ; + } // this is only used once when the Track informations are retrived from the primary track - AncestorsBirthVolume.push_back(aTrack->GetNextVolume()->GetName()) ; - // Mhd : Redundant information for the first element but i'll keep it for the sake of coherence' - AncestorsPdg.push_back(aTrack->GetDefinition()->GetPDGEncoding()); + SecondariesBirthVolume.push_back(aTrack->GetNextVolume()->GetName()) ; + // Mhd : Redundant information for the first element but i'll keep it for the sake of coherence + SecondariesPdg.push_back(aTrack->GetDefinition()->GetPDGEncoding()); // We cant do this step since the particle isn't dead yet - // AncestorsBirthVolume.push_back(aTrack->GetNextVolume()->GetName()) ; - - /* - // MHD 25 april 2013 - // detect change of material and make sure the information object is already - // const G4Step* GetStep() const; aTrack->GetStep() ; // Not working - if( ( aTrack->GetStep()->GetPreStepPoint()->GetMaterial() != aTrack->GetStep()->GetPostStepPoint()->GetMaterial() ) && aTrack->GetUserInformation()!=0 ) - { - G4cout<GetStep()->GetPreStepPoint()->GetMaterial() << " -> " <GetStep()->GetPreStepPoint()->GetMaterial()<GetStep()->GetPostStepPoint()->GetMaterial()->GetDensity()); - } - */ + // SecondariesDeathVolume.push_back(aTrack->GetNextVolume()->GetName()) ; } @@ -56,16 +119,42 @@ TrackInformation::TrackInformation(const G4Track* aTrack) TrackInformation::TrackInformation(const TrackInformation* aTrackInfo) : G4VUserTrackInformation() { + Tagged = aTrackInfo->Tagged; + + //original at source + originalParentID = aTrackInfo->originalParentID; originalTrackID = aTrackInfo->originalTrackID; originalPdg = aTrackInfo->originalPdg; originalPosition = aTrackInfo->originalPosition; originalMomentum = aTrackInfo->originalMomentum; originalEnergy = aTrackInfo->originalEnergy; originalTime = aTrackInfo->originalTime; + // after impact + originalImpactVolume = aTrackInfo->originalImpactVolume; + originalImpactPosition = aTrackInfo->originalImpactPosition; + originalImpactMomentum = aTrackInfo->originalImpactMomentum; + + // secondaries + currentProcess = aTrackInfo->currentProcess; + currentTrackID = aTrackInfo->currentTrackID; + currentparentTrackID = aTrackInfo->currentparentTrackID ; + currentPdg = aTrackInfo->currentPdg ; - AncestorsBirthVolume = aTrackInfo->AncestorsBirthVolume; - AncestorsDeathVolume = aTrackInfo->AncestorsDeathVolume; - AncestorsPdg = aTrackInfo->AncestorsPdg; + currentPositionAtVertex= aTrackInfo->currentPositionAtVertex; + currentPositionAtDetector= aTrackInfo->currentPositionAtDetector; + currentPositionAtDeath= aTrackInfo->currentPositionAtDeath; // death or out of world + + currentMomentumAtVertex= aTrackInfo->currentMomentumAtVertex; // contains the angle at the emission + currentMomentumAtDetector= aTrackInfo->currentMomentumAtDetector; // contains angle at the impact on the detector + currentMomentumAtDeath= aTrackInfo->currentMomentumAtDeath; // death or out of world + + currentEnergyAtVertex = aTrackInfo->currentEnergyAtVertex; + currentTimeAtVertex = aTrackInfo->currentTimeAtVertex; + + SecondariesProcess = aTrackInfo->SecondariesProcess; + SecondariesBirthVolume = aTrackInfo->SecondariesBirthVolume; + SecondariesDeathVolume = aTrackInfo->SecondariesDeathVolume; + SecondariesPdg = aTrackInfo->SecondariesPdg; } @@ -75,25 +164,68 @@ TrackInformation::~TrackInformation(){;} void TrackInformation::Print() const { - G4cout - << " at " << originalPosition - << "Original track ID " << originalTrackID << G4endl - << "Original track Pdg " << originalPdg << G4endl + + G4cout << " ============================================================= " << G4endl ; + + G4cout << " +++++ Original at source +++++ " << G4endl ; + if (Tagged) + G4cout << " This Track is tagged !! " << Tagged << G4endl ; + else + G4cout << " No tagging is set " << Tagged << G4endl ; + + G4cout << " at Source " << G4endl + << " Position " << originalPosition << G4endl + << " Momentum " << originalMomentum << G4endl + << "Original parent ID " << originalParentID << G4endl + << "Original track ID " << originalTrackID << G4endl + << "Original track Pdg " << originalPdg << G4endl << "Original track Energy " << originalEnergy << G4endl - << "Original track Time " << originalTime << G4endl ; + << "Original track Time " << originalTime << G4endl ; + + G4cout << " ----- Original after first impact ----- " << G4endl ; + G4cout << " at Volume " << originalImpactVolume << G4endl + << " Position " << originalImpactPosition << G4endl ; + G4cout << " Momentum " << originalImpactMomentum << G4endl ; + G4cout << " +++++ Secondaries +++++ " << G4endl ; + G4cout << "Current track Process " << currentProcess << G4endl + << "Current parent ID " << currentparentTrackID << G4endl + << "Current track ID " << currentTrackID << G4endl + << "Current track Pdg " << currentPdg << G4endl + << "Current track Energy " << currentEnergyAtVertex << G4endl + << "Current track time " << currentTimeAtVertex << G4endl ; + G4cout << " ----- Secondary key impact ----- " << G4endl ; + G4cout << " at Vertex " << G4endl + << " Position " << currentPositionAtVertex << G4endl + << " Momentum " << currentMomentumAtVertex << G4endl ; + + G4cout << " at 1st impact " << currentImpactVolume << G4endl + << " Position " << currentImpactPosition << G4endl ; + G4cout << " Momentum " << currentImpactMomentum << G4endl ; - for(unsigned i = 0 ; i < AncestorsBirthVolume.size() ; i++ ) + G4cout << " at Detector " << G4endl + << " Position " << currentPositionAtDetector << G4endl + << " Momentum " << currentMomentumAtDetector << G4endl ; + + G4cout << " at Death " << G4endl + << " Position " << currentPositionAtDeath << G4endl + << " Momentum " << currentMomentumAtDeath << G4endl ; + + for(unsigned i = 0 ; i < SecondariesProcess.size() ; i++ ) + G4cout + << " SecondariesProcess at " <SetUserInformation(anInfo); - } // Mhd 02 May 2013 - if(aTrack->GetParentID()!=0 /*&& aTrack->GetUserInformation()!=0*/) // I think we dont need the second condition - { + if(aTrack->GetParentID()!=0) { - TrackInformation* oldInfo = (TrackInformation*)aTrack->GetUserInformation(); - TrackInformation* newInfo = oldInfo; - G4Track* theTrack = (G4Track*)aTrack; - - //Set (append) the Pdg of the ancestor particle - newInfo->SetAncestorsPdgElement(theTrack->GetDefinition()->GetPDGEncoding()); - - //Set (append) the Birth volume of the ancestor particle - if( theTrack->GetNextVolume() != 0 ) - { - newInfo->SetAncestorsBirthVolumeElement(theTrack->GetNextVolume()->GetName()); - } -else - { - newInfo->SetAncestorsBirthVolumeElement("OutOfWorld"); - } + TrackInformation* oldInfo = (TrackInformation*)aTrack->GetUserInformation(); + TrackInformation* newInfo = oldInfo; + G4Track* theTrack = (G4Track*)aTrack; + + //Set (append) the Pdg of the ancestor particle + //newInfo->SetCurrentParentID(newInfo->GetCurrentTrackID()); + newInfo->SetCurrentParentID(theTrack->GetParentID()); + + //Set (append) the Pdg of the ancestor particle + newInfo->SetCurrentTrackID(theTrack->GetTrackID()); + + //Set (append) the Pdg of the ancestor particle + newInfo->SetSecondariesPdgElement(theTrack->GetDefinition()->GetPDGEncoding()); + + //Set (append) the process name at birth of the ancestor particle + if( theTrack->GetNextVolume() != 0 ) { + newInfo->SetSecondariesProcessElement(theTrack->GetCreatorProcess()->GetProcessName()); + } + else { + newInfo->SetSecondariesProcessElement("OOW/Dead"); // OOW = OutOfWorld + } + + //Set (append) the Birth volume of the ancestor particle + if( theTrack->GetNextVolume() != 0 ) { + newInfo->SetSecondariesBirthVolumeElement(theTrack->GetNextVolume()->GetName()); + } + else { + newInfo->SetSecondariesBirthVolumeElement("OOW/Dead"); + } - // set (append) the new user information - theTrack->SetUserInformation(newInfo); - + // set (append) the new user information + theTrack->SetUserInformation(newInfo); } @@ -63,12 +73,12 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) if( theTrack->GetNextVolume() != 0 ) { //G4cout << " Adding the next DEATH volume "<< theTrack->GetNextVolume()->GetName() <SetAncestorsDeathVolumeElement(theTrack->GetNextVolume()->GetName()); + newInfo->SetSecondariesDeathVolumeElement(theTrack->GetNextVolume()->GetName()); } else { //G4cout << std::setw(11) << "OutOfWorld" << " "<SetAncestorsDeathVolumeElement("OutOfWorld"); + newInfo->SetSecondariesDeathVolumeElement("OOW/Dead"); } // set (append) the new user information @@ -92,5 +102,6 @@ else } } } + } diff --git a/src/nonUniformMagneticField.cc b/src/nonUniformMagneticField.cc index 7c7b25f..110eabb 100755 --- a/src/nonUniformMagneticField.cc +++ b/src/nonUniformMagneticField.cc @@ -24,10 +24,10 @@ // // Constructors: -nonUniformMagneticField::nonUniformMagneticField(const char* fieldName, G4double zOffset) +nonUniformMagneticField::nonUniformMagneticField(const char* fieldName="./", G4double zOffset=0., G4double zRotation=0.) : fChordFinder(0), fStepper(0) { - fMagneticField = new TabulatedMagneticField(fieldName, zOffset); + fMagneticField = new TabulatedMagneticField(fieldName, zOffset, zRotation); GetGlobalFieldManager()->CreateChordFinder(fMagneticField); // fFieldMessenger = new F03FieldMessenger(this) ;