diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx index 6f6bbe9e598a8052fef38f06753f6114a3416255..b7ef5c16f6dec659ed8bb107ec6796b5c1a26b46 100644 --- a/NPLib/Core/NPDetectorManager.cxx +++ b/NPLib/Core/NPDetectorManager.cxx @@ -260,6 +260,7 @@ void NPL::DetectorManager::BuildPhysicalEvent(){ static std::map<std::string,VDetector*>::iterator end= m_Detector.end(); for (it =begin; it != end; ++it) { + cout << " " <<endl; (it->second->*m_ClearEventPhysicsPtr)(); (it->second->*m_BuildPhysicalPtr)(); if(m_FillSpectra){ diff --git a/NPLib/Detectors/Mugast/CMakeLists.txt b/NPLib/Detectors/Mugast/CMakeLists.txt index 31661284318085dd85fa9d20dcd273a526653967..e11a41e88cd4b9110560a0ceb913d6f962f95fb9 100644 --- a/NPLib/Detectors/Mugast/CMakeLists.txt +++ b/NPLib/Detectors/Mugast/CMakeLists.txt @@ -4,4 +4,4 @@ add_custom_command(OUTPUT TMugastDataDict.cxx COMMAND ../../scripts/build_dict.s add_library(NPMugast SHARED TMugastData.cxx TMugastDataDict.cxx TMugastPhysics.cxx TMugastPhysicsDict.cxx TMugastSpectra.h TMugastSpectra.cxx) target_link_libraries(NPMugast ${ROOT_LIBRARIES} NPCore) #install(FILES TMugastData.h TMugastPhysics.h TMugastSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) -install(FILES MugastMap.h TMugastSpectra.h TMugastData.h TMugastPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) +install(FILES MugastMap.h MugastReverseMap.h TMugastSpectra.h TMugastData.h TMugastPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) diff --git a/NPLib/Detectors/Mugast/MugastMap.h b/NPLib/Detectors/Mugast/MugastMap.h index d590baaf19b20af5ee8bfc0e56b1b489d0b08366..a0424c457e8f29bf79c37167c09972fe7f29931c 100644 --- a/NPLib/Detectors/Mugast/MugastMap.h +++ b/NPLib/Detectors/Mugast/MugastMap.h @@ -1,6 +1,7 @@ #ifndef MUGASTMAP #define MUGASTMAP +#include "MugastReverseMap.h" namespace MUGAST_MAP{ diff --git a/NPLib/Detectors/Mugast/MugastReverseMap.h b/NPLib/Detectors/Mugast/MugastReverseMap.h new file mode 100644 index 0000000000000000000000000000000000000000..b9c4d26cd56af31ebcf1a0df87b8cc6ee55036bd --- /dev/null +++ b/NPLib/Detectors/Mugast/MugastReverseMap.h @@ -0,0 +1,792 @@ +#ifndef MUGASTREVERSEMAP +#define MUGASTREVERSEMAP +namespace MUGAST_MAP{ + + const int ReverseTrapezeX[128] ={ + 64, + 57, + 62, + 55, + 60, + 53, + 58, + 51, + 56, + 49, + 54, + 59, + 52, + 61, + 50, + 63, + 47, + 48, + 45, + 46, + 43, + 44, + 41, + 42, + 39, + 40, + 37, + 38, + 35, + 36, + 33, + 34, + 29, + 30, + 27, + 28, + 25, + 26, + 23, + 24, + 21, + 22, + 19, + 20, + 17, + 18, + 32, + 16, + 13, + 2, + 11, + 4, + 9, + 6, + 1, + 8, + 3, + 10, + 7, + 12, + 5, + 14, + 15, + 31, + 127, + 128, + 125, + 126, + 123, + 124, + 121, + 122, + 119, + 120, + 117, + 118, + 115, + 116, + 113, + 114, + 111, + 112, + 109, + 110, + 107, + 108, + 105, + 106, + 103, + 104, + 101, + 102, + 99, + 100, + 97, + 98, + 93, + 94, + 91, + 92, + 89, + 90, + 87, + 88, + 85, + 86, + 83, + 84, + 81, + 82, + 79, + 80, + 77, + 78, + 75, + 76, + 73, + 74, + 71, + 72, + 69, + 70, + 67, + 68, + 65, + 66, + 96, + 95, + }; + + const int ReverseTrapezeY[128] ={ + 94, + 97, + 92, + 100, + 90, + 102, + 88, + 104, + 86, + 106, + 84, + 108, + 82, + 110, + 80, + 112, + 95, + 114, + 78, + 116, + 76, + 118, + 74, + 120, + 72, + 122, + 70, + 124, + 68, + 126, + 65, + 127, + 73, + 98, + 71, + 99, + 69, + 101, + 67, + 103, + 66, + 128, + 79, + 125, + 77, + 123, + 75, + 121, + 81, + 119, + 83, + 117, + 85, + 115, + 87, + 113, + 89, + 111, + 91, + 109, + 93, + 107, + 96, + 105, + 32, + 38, + 31, + 36, + 27, + 35, + 24, + 33, + 21, + 39, + 29, + 37, + 22, + 41, + 34, + 42, + 19, + 40, + 30, + 43, + 20, + 46, + 28, + 48, + 17, + 45, + 25, + 50, + 18, + 9, + 26, + 44, + 23, + 47, + 13, + 52, + 11, + 54, + 15, + 49, + 16, + 56, + 14, + 58, + 12, + 60, + 10, + 62, + 8, + 63, + 6, + 64, + 4, + 61, + 1, + 59, + 2, + 57, + 3, + 55, + 5, + 53, + 7, + 51, + }; + + const int ReverseSquareX[128] ={ + 4, + 128, + 6, + 8, + 10, + 12, + 14, + 16, + 18, + 20, + 22, + 24, + 26, + 28, + 30, + 32, + 34, + 1, + 3, + 5, + 7, + 9, + 11, + 13, + 15, + 17, + 19, + 21, + 23, + 25, + 27, + 29, + 31, + 33, + 35, + 36, + 38, + 40, + 42, + 44, + 46, + 48, + 50, + 52, + 54, + 56, + 58, + 60, + 62, + 64, + 37, + 39, + 41, + 43, + 45, + 47, + 49, + 51, + 53, + 55, + 57, + 59, + 61, + 63, + 2, + 65, + 67, + 69, + 71, + 73, + 75, + 77, + 79, + 81, + 83, + 85, + 87, + 89, + 91, + 93, + 95, + 97, + 99, + 101, + 103, + 105, + 107, + 109, + 111, + 113, + 115, + 117, + 119, + 121, + 123, + 125, + 127, + 66, + 68, + 70, + 72, + 74, + 76, + 78, + 80, + 82, + 84, + 86, + 88, + 90, + 92, + 94, + 96, + 98, + 100, + 102, + 104, + 106, + 108, + 110, + 112, + 114, + 116, + 118, + 120, + 122, + 124, + 126, + }; + + const int ReverseSquareY[128] ={ + 125, + 19, + 123, + 121, + 119, + 117, + 115, + 113, + 111, + 109, + 107, + 105, + 103, + 101, + 99, + 98, + 96, + 128, + 127, + 126, + 124, + 122, + 120, + 118, + 116, + 114, + 112, + 110, + 108, + 106, + 104, + 102, + 100, + 97, + 95, + 93, + 91, + 89, + 87, + 85, + 83, + 81, + 79, + 77, + 75, + 73, + 71, + 69, + 67, + 66, + 94, + 92, + 90, + 88, + 86, + 84, + 82, + 80, + 78, + 76, + 74, + 72, + 70, + 68, + 65, + 63, + 62, + 60, + 58, + 56, + 54, + 52, + 50, + 48, + 46, + 44, + 42, + 40, + 38, + 36, + 33, + 31, + 30, + 28, + 26, + 24, + 22, + 20, + 18, + 16, + 14, + 12, + 10, + 8, + 6, + 4, + 1, + 64, + 61, + 59, + 57, + 55, + 53, + 51, + 49, + 47, + 45, + 43, + 41, + 39, + 37, + 35, + 34, + 32, + 29, + 27, + 25, + 23, + 21, + 17, + 15, + 13, + 11, + 9, + 7, + 5, + 3, + 2, + }; + + const int ReverseAnnularX[128] ={ + 64, + 62, + 63, + 60, + 61, + 58, + 59, + 56, + 57, + 54, + 55, + 52, + 53, + 50, + 51, + 48, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + }; + + const int ReverseAnnularY[128] ={ + 46, + 50, + 48, + 44, + 40, + 36, + 26, + 22, + 18, + 14, + 16, + 20, + 24, + 28, + 38, + 42, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + }; + +} +#endif \ No newline at end of file diff --git a/NPLib/Detectors/Mugast/ReverseMap.cxx b/NPLib/Detectors/Mugast/ReverseMap.cxx new file mode 100644 index 0000000000000000000000000000000000000000..74425e45808add7e0ad4c2dbcddd7a65497ec95d --- /dev/null +++ b/NPLib/Detectors/Mugast/ReverseMap.cxx @@ -0,0 +1,110 @@ +// Generate the reversed map based on on the MugastMap.h file + +#include "MugastMap.h" +#include <fstream> +#include <iostream> +void ReverseMap(){ + std::ofstream f("MugastReverseMap.h"); + ///////////////// + // File header // + ///////////////// + f<< "#ifndef MUGASTREVERSEMAP" << endl; + f<< "#define MUGASTREVERSEMAP" << endl; + f<< "namespace MUGAST_MAP{" << endl << endl; + + ////////////////////////////////////////////////////////////////////////////// + // Trapeze X + unsigned int size = 128; + // build the array + int ReverseTX[128]; + for(unsigned int i = 0 ; i < size ; i++){ + ReverseTX[MUGAST_MAP::TrapezeX[i]-1]=i+1; + } + // write array + f << "\tconst int ReverseTrapezeX[128] ={" << endl; + for(unsigned int i = 0 ; i < size ; i++){ + f <<"\t\t" << ReverseTX[i] << ","<< endl ; + } + f<<"\t\t};" << endl << endl; + + ////////////////////////////////////////////////////////////////////////////// + // Trapeze Y + // build the array + int ReverseTY[128]; + for(unsigned int i = 0 ; i < size ; i++){ + ReverseTY[MUGAST_MAP::TrapezeY[i]-1]=i+1; + } + // write array + f << "\tconst int ReverseTrapezeY[128] ={" << endl; + for(unsigned int i = 0 ; i < size ; i++){ + f <<"\t\t" << ReverseTY[i] << ","<< endl ; + } + f<<"\t\t};" << endl << endl; + + ////////////////////////////////////////////////////////////////////////////// + // Square X + // build the array + int ReverseSX[128]; + for(unsigned int i = 0 ; i < size ; i++){ + ReverseSX[MUGAST_MAP::SquareX[i]-1]=i+1; + } + // write array + f << "\tconst int ReverseSquareX[128] ={" << endl; + for(unsigned int i = 0 ; i < size ; i++){ + f <<"\t\t" << ReverseSX[i] << ","<< endl ; + } + f<<"\t\t};" << endl << endl; + + ////////////////////////////////////////////////////////////////////////////// + // Square Y + // build the array + int ReverseSY[128]; + for(unsigned int i = 0 ; i < size ; i++){ + ReverseSY[MUGAST_MAP::SquareY[i]-1]=i+1; + } + // write array + f << "\tconst int ReverseSquareY[128] ={" << endl; + for(unsigned int i = 0 ; i < size ; i++){ + f <<"\t\t" << ReverseSY[i] << ","<< endl ; + } + f<<"\t\t};" << endl << endl; + + ////////////////////////////////////////////////////////////////////////////// + // Annular X + // build the array + int ReverseAX[128]={128}; + for(unsigned int i = 0 ; i < size ; i++){ + if(MUGAST_MAP::AnnularX[i]-1<16) + ReverseAX[MUGAST_MAP::AnnularX[i]-1]=i+1; + } + // write array + f << "\tconst int ReverseAnnularX[128] ={" << endl; + for(unsigned int i = 0 ; i < size ; i++){ + f <<"\t\t" << ReverseAX[i] << ","<< endl ; + } + f<<"\t\t};" << endl << endl; + + ////////////////////////////////////////////////////////////////////////////// + // Annular Y + // build the array + int ReverseAY[128]={128}; + for(unsigned int i = 0 ; i < size ; i++){ + if(MUGAST_MAP::AnnularY[i]-1<16) + ReverseAY[MUGAST_MAP::AnnularY[i]-1]=i+1; + } + // write array + f << "\tconst int ReverseAnnularY[128] ={" << endl; + for(unsigned int i = 0 ; i < size ; i++){ + f <<"\t\t" << ReverseAY[i] << ","<< endl ; + } + f<<"\t\t};" << endl << endl; + + + + ///////////////// + // file footer // + ///////////////// + f<<"}" << endl; + f<<"#endif"; + f.close(); +} diff --git a/NPLib/Detectors/Mugast/TMugastData.cxx b/NPLib/Detectors/Mugast/TMugastData.cxx index 5ad743a5a0bd2dd925980e3911b587953b1ff133..341b8ae05536444bb5efbd84c4194535458c83ed 100644 --- a/NPLib/Detectors/Mugast/TMugastData.cxx +++ b/NPLib/Detectors/Mugast/TMugastData.cxx @@ -6,9 +6,9 @@ *****************************************************************************/ /***************************************************************************** - * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk * + * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr * * * - * Creation Date : November 2018 * + * Creation Date : Novembre 2018 * * Last update : * *---------------------------------------------------------------------------* * Decription: * diff --git a/NPLib/Detectors/Mugast/TMugastData.h b/NPLib/Detectors/Mugast/TMugastData.h index 39aea243bc3db0aba8e939460de37d68f6f567b4..00dcfa0f7030a463da4f6a2732852deb1c8ef704 100644 --- a/NPLib/Detectors/Mugast/TMugastData.h +++ b/NPLib/Detectors/Mugast/TMugastData.h @@ -8,9 +8,9 @@ *****************************************************************************/ /***************************************************************************** - * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk * + * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr * * * - * Creation Date : November 2018 * + * Creation Date : Novembre 2018 * * Last update : * *---------------------------------------------------------------------------* * Decription: * @@ -200,12 +200,12 @@ public: // (X,E) inline unsigned short GetDSSDXEMult() const {return fMG_DSSDXE_DetectorNbr.size();} inline unsigned short GetDSSDXEDetectorNbr(const int& i) const {return fMG_DSSDXE_DetectorNbr[i];} - inline unsigned short GetDSSDXEStripNbr(const int& i) const {return fMG_DSSDXE_StripNbr[i];} + inline unsigned short GetDSSDXEStripNbr(const int& i) const {return fMG_DSSDXE_StripNbr[i];} inline double GetDSSDXEEnergy(const int& i) const {return fMG_DSSDXE_Energy[i];} // (X,T) inline unsigned short GetDSSDXTMult() const {return fMG_DSSDXT_DetectorNbr.size();} inline unsigned short GetDSSDXTDetectorNbr(const int& i) const {return fMG_DSSDXT_DetectorNbr[i];} - inline unsigned short GetDSSDXTStripNbr(const int& i) const {return fMG_DSSDXT_StripNbr[i];} + inline unsigned short GetDSSDXTStripNbr(const int& i) const {return fMG_DSSDXT_StripNbr[i];} inline double GetDSSDXTTime(const int& i) const {return fMG_DSSDXT_Time[i];} // (Y,E) inline unsigned short GetDSSDYEMult() const {return fMG_DSSDYE_DetectorNbr.size();} diff --git a/NPLib/Detectors/Mugast/TMugastPhysics.cxx b/NPLib/Detectors/Mugast/TMugastPhysics.cxx index 41d6ec299f0225697d73ab6b0903241096e0f29f..147815084499e866d5ec2a52b812c9c13e2ba9ef 100644 --- a/NPLib/Detectors/Mugast/TMugastPhysics.cxx +++ b/NPLib/Detectors/Mugast/TMugastPhysics.cxx @@ -6,13 +6,13 @@ *****************************************************************************/ /***************************************************************************** - * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk * + * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr * * * - * Creation Date : febuary 2009 * - * Last update : october 2010 * + * Creation Date : novembre 2018 * + * Last update : * *---------------------------------------------------------------------------* * Decription: * - * This class hold Mugast treated data * + * This class hold Mugast treated data * * * *---------------------------------------------------------------------------* * Comment: * @@ -544,7 +544,7 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) { TVector3 C = blocks[i]->GetTVector3("X1_Y128", "mm"); TVector3 D = blocks[i]->GetTVector3("X128_Y128", "mm"); - AddTelescope(A, B, C, D); + AddTelescope(DetectorType[detectorNbr],A, B, C, D); } else if (blocks[i]->HasTokenList(annular)) { @@ -561,7 +561,7 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) { det = i+1; m_DetectorNumberIndex[detectorNbr]=det; TVector3 Center = blocks[i]->GetTVector3("Center", "mm"); - AddTelescope(Center); + AddTelescope(DetectorType[detectorNbr],Center); } @@ -583,7 +583,7 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) { double Phi = blocks[i]->GetDouble("PHI", "deg"); double R = blocks[i]->GetDouble("R", "mm"); vector<double> beta = blocks[i]->GetVectorDouble("BETA", "deg"); - AddTelescope(Theta, Phi, R, beta[0], beta[1], beta[2]); + AddTelescope(DetectorType[detectorNbr],Theta, Phi, R, beta[0], beta[1], beta[2]); } else { @@ -728,7 +728,7 @@ void TMugastPhysics::InitializeRootOutput() { } ///// Specific to MugastArray //// -void TMugastPhysics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1, +void TMugastPhysics::AddTelescope(MG_DetectorType type,TVector3 C_X1_Y1, TVector3 C_X128_Y1, TVector3 C_X1_Y128, TVector3 C_X128_Y128) { // To avoid warning C_X128_Y128 *= 1; @@ -750,14 +750,16 @@ void TMugastPhysics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1, // Geometry Parameter double Base,Height; - if(DetectorType[m_NumberOfTelescope-1]==MG_TRAPEZE){ + if(type==MG_TRAPEZE){ Base = 91.48; // mm Height = 104.688; // mm } - if(DetectorType[m_NumberOfTelescope-1]==MG_SQUARE){ + if(type==MG_SQUARE){ Base = 91.716; // mm Height = 94.916; // mm +// + // Height = 194.916; // mm } //double Face = 98; // mm double NumberOfStrip = 128; @@ -773,10 +775,13 @@ void TMugastPhysics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1, vector<vector<double>> OneTelescopeStripPositionZ; // Moving StripCenter to 1.1 corner: - //Strip_1_1 = C_X1_Y1 + (U + V) * (StripPitch / 2.); - Strip_1_1 = C_X1_Y1 + U * (StripPitchBase / 2.) + V * (StripPitchHeight / 2.); - //Strip_1_1 += U + V ; - + // Strip_1_1 = C_X1_Y1 + U * (StripPitchBase / 2.) + V * (StripPitchHeight / 2.); + // This calculation recenter the strip around the detector center. + // This account for cases where the provided corner coordinates + // does not match the detector size + TVector3 Center = 0.25*(C_X1_Y128 + C_X128_Y128 + C_X1_Y1 + C_X128_Y1); + Strip_1_1 = Center-(0.5*Base*U+0.5*Height*V) + U * (StripPitchBase / 2.) + V * (StripPitchHeight / 2.); + for (int i = 0; i < 128; ++i) { lineX.clear(); lineY.clear(); @@ -834,7 +839,7 @@ void TMugastPhysics::InitializeStandardParameter() { } //////////////////////////////////////////////////////////////////////////////// -void TMugastPhysics::AddTelescope(TVector3 C_Center) { +void TMugastPhysics::AddTelescope(MG_DetectorType type,TVector3 C_Center) { // To avoid warning m_NumberOfTelescope++; double Z = C_Center.Z(); @@ -938,7 +943,7 @@ void TMugastPhysics::AddTelescope(TVector3 C_Center) { } //////////////////////////////////////////////////////////////////////////////// -void TMugastPhysics::AddTelescope(double theta, double phi, double distance, +void TMugastPhysics::AddTelescope(MG_DetectorType type,double theta, double phi, double distance, double beta_u, double beta_v, double beta_w) { m_NumberOfTelescope++; diff --git a/NPLib/Detectors/Mugast/TMugastPhysics.h b/NPLib/Detectors/Mugast/TMugastPhysics.h index 7094a14c55f070ec807dcfc1783c73a346982071..a93f38af9d4b75736cb9bea21608f6b1b5d18ac2 100644 --- a/NPLib/Detectors/Mugast/TMugastPhysics.h +++ b/NPLib/Detectors/Mugast/TMugastPhysics.h @@ -1,26 +1,24 @@ #ifndef TMUGASTPHYSICS_H #define TMUGASTPHYSICS_H /***************************************************************************** - * Copyright (C) 2009-2016 this file is part of the NPTool Project * + * Copyright (C) 2009-2016 this file is part of the NPTool Project * * * * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * * For the list of contributors see $NPTOOL/Licence/Contributors * *****************************************************************************/ /***************************************************************************** - * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk * + * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr * * * - * Creation Date : febuary 2009 * + * Creation Date : novembre 2018 * * Last update : * *---------------------------------------------------------------------------* * Decription: * - * This class hold must2 treated data * + * This class hold Mugast treated data * * * *---------------------------------------------------------------------------* * Comment: * * * - * * - * * *****************************************************************************/ // STL #include <map> @@ -167,15 +165,15 @@ class TMugastPhysics : public TObject, public NPL::VDetector { void ReadAnalysisConfig(); // Add a Telescope using Corner Coordinate information - void AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1, TVector3 C_X1_Y128, + void AddTelescope(MG_DetectorType type,TVector3 C_X1_Y1, TVector3 C_X128_Y1, TVector3 C_X1_Y128, TVector3 C_X128_Y128); // Add a Telescope using R Theta Phi of Si center information - void AddTelescope(double theta, double phi, double distance, double beta_u, + void AddTelescope(MG_DetectorType type,double theta, double phi, double distance, double beta_u, double beta_v, double beta_w); // Special Method for Annular S1 - void AddTelescope(TVector3 C_Center); + void AddTelescope(MG_DetectorType type,TVector3 C_Center); // Give and external TMustData object to TMugastPhysics. Needed for online // analysis for example. diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc index 53bede6d67ac64bf5571cde5e8d89724e374eb01..117b02d53a3179190a62479d9bdee6f78bd3643a 100644 --- a/NPSimulation/Detectors/MUST2/MUST2Array.cc +++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc @@ -903,7 +903,7 @@ void MUST2Array::ReadSensitive(const G4Event*) { for(it=mapBack.begin();it!=mapBack.end();it++){ double energyY = RandGauss::shoot(it->second.first, ResoStrip); - double timeX = TimeOffset - RandGauss::shoot(it->second.second, ResoTimeMust); + double timeY = TimeOffset - RandGauss::shoot(it->second.second, ResoTimeMust); unsigned int strip = it->first-1000000*(it->first/1000000); unsigned int det = it->first/1000000; if (energyY > ThresholdSi) { @@ -912,7 +912,7 @@ void MUST2Array::ReadSensitive(const G4Event*) { m_Event->SetStripYE(det, strip , NPL::EnergyToADC(energyY, 0, 63, 8192, 0)); m_Event->SetStripYT(det, strip , - NPL::EnergyToADC(timeX, 0, 1000, 8192, 16384)); + NPL::EnergyToADC(timeY, 0, 1000, 8192, 16384)); } } diff --git a/NPSimulation/Detectors/Mugast/Mugast.cc b/NPSimulation/Detectors/Mugast/Mugast.cc index 0fd6ec3a3116705725f816d2bc757b745802efc1..18e370814824a179fb19efea49bcde47026abba1 100644 --- a/NPSimulation/Detectors/Mugast/Mugast.cc +++ b/NPSimulation/Detectors/Mugast/Mugast.cc @@ -41,13 +41,16 @@ // NPTool header #include "Mugast.hh" -#include "CalorimeterScorers.hh" +#include "MugastReverseMap.h" +#include "DSSDScorers.hh" #include "InteractionScorers.hh" #include "RootOutput.h" #include "MaterialManager.hh" #include "NPSDetectorFactory.hh" #include "NPOptionManager.h" #include "NPSHitsMap.hh" +#include "NPCore.h" + // CLHEP header #include "CLHEP/Random/RandGauss.h" @@ -58,20 +61,25 @@ using namespace CLHEP; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... namespace Mugast_NS{ // Resolution - const G4double SigmaTime = 0.212765957 ;// = 500ps - const G4double SigmaEnergy = 0.0149 ;// 0.0223 = 52keV of Resolution // Unit is MeV/2.35 14.861996 + const G4double SigmaTime = 0.212765957 ;// = 500ps + const G4double SigmaEnergy = 0.0149 ;// 0.0223 = 52keV of Resolution // Unit is MeV/2.35 14.861996 const G4double TimeOffset = 500 ;// 500 ns stop // Threshold - const G4double ThresholdSi = 500 * keV; + const G4double EnergyThreshold = 1 * MeV; // Geometry const G4double AluStripThickness = 0.4*micrometer ; - const G4double SiliconThickness = 300*micrometer ; - const G4double SiliconFace = 100.42*mm ; + const G4double SiliconThickness = 3000*micrometer ; + // Square + + // const G4double SquareBase = 88*mm; +// const G4double SquareHeight = 87*mm; + const G4double SquareBase = 91.716*mm; const G4double SquareHeight = 94.916*mm; + // const G4double SquareHeight = 194.916*mm; const G4double SquareLength = 1*cm; // Trapezoid const G4double TrapezoidBaseLarge = 91.48*mm; @@ -87,7 +95,9 @@ using namespace Mugast_NS; // Mugast Specific Method Mugast::Mugast(){ m_Event = new TMugastData() ; - m_MugastScorer = 0; + m_SquareScorer= 0; + m_TrapezoidScorer= 0; + m_AnnularScorer= 0; m_SquareDetector = 0; m_TrapezoidDetector = 0; m_AnnularDetector = 0; @@ -116,7 +126,7 @@ G4LogicalVolume* Mugast::BuildSquareDetector(){ if(!m_SquareDetector){ G4String Name = "MugastSquare"; - G4Box* solidSquare = new G4Box(Name, 0.5*SquareHeight, 0.5*SquareHeight, 0.5*SquareLength); + G4Box* solidSquare = new G4Box(Name, 0.5*SquareBase, 0.5*SquareHeight, 0.5*SquareLength); G4LogicalVolume* logicSquare = new G4LogicalVolume(solidSquare, MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"), Name, 0, 0, 0); @@ -126,7 +136,7 @@ G4LogicalVolume* Mugast::BuildSquareDetector(){ logicSquare->SetVisAttributes(SquareVisAtt); // Silicon detector itself - G4ThreeVector positionFirstStage = G4ThreeVector(0, 0, 0.5*cm); + G4ThreeVector positionFirstStage = G4ThreeVector(0, 0, 0); G4Box* solidFirstStage = new G4Box("solidFirstStage", 0.5*SquareBase, 0.5*SquareHeight, 0.5*SiliconThickness); G4LogicalVolume* logicFirstStage = new G4LogicalVolume(solidFirstStage, @@ -137,10 +147,10 @@ G4LogicalVolume* Mugast::BuildSquareDetector(){ new G4PVPlacement(0, positionFirstStage, logicFirstStage, Name + "_FirstStage", logicSquare, false, 0); m_SquareDetector = logicSquare; // Set First Stage sensible - // logicFirstStage->SetSensitiveDetector(m_FirstStageScorer); + logicFirstStage->SetSensitiveDetector(m_SquareScorer); ///Visualisation of FirstStage Strip - G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); // blue + G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); logicFirstStage->SetVisAttributes(FirstStageVisAtt); } @@ -150,6 +160,7 @@ G4LogicalVolume* Mugast::BuildSquareDetector(){ G4LogicalVolume* Mugast::BuildTrapezoidDetector(){ if(!m_SquareDetector){ G4String Name = "MugastTrapezoid"; + // Definition of the volume containing the sensitive detector G4Trap* solidTrapezoid = new G4Trap(Name, TrapezoidLength*0.5, 0*deg, 0*deg, @@ -165,7 +176,7 @@ G4LogicalVolume* Mugast::BuildTrapezoidDetector(){ logicTrapezoid->SetVisAttributes(TrapezoideVisAtt); // Silicon detector itself - G4ThreeVector positionFirstStage = G4ThreeVector(0, 0, 0.5*SiliconThickness); + G4ThreeVector positionFirstStage = G4ThreeVector(0, 0, 0); G4Trap* solidFirstStage = new G4Trap("solidFirstStage", 0.5*SiliconThickness, 0*deg, 0*deg, @@ -185,10 +196,10 @@ G4LogicalVolume* Mugast::BuildTrapezoidDetector(){ 0); // Set First Stage sensible - logicFirstStage->SetSensitiveDetector(m_MugastScorer); + logicFirstStage->SetSensitiveDetector(m_TrapezoidScorer); ///Visualisation of FirstStage Strip - G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); // blue + G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); logicFirstStage->SetVisAttributes(FirstStageVisAtt); m_TrapezoidDetector=logicTrapezoid; @@ -270,7 +281,7 @@ void Mugast::ConstructDetector(G4LogicalVolume* world){ u = m_X128_Y1[i] - m_X1_Y1[i]; u = u.unit(); - v = m_X1_Y128[i] - m_X1_Y1[i]; + v = (m_X1_Y128[i] + m_X128_Y128[i] - m_X1_Y1[i] - m_X128_Y1[i]); v = v.unit(); w = u.cross(v); @@ -281,7 +292,7 @@ void Mugast::ConstructDetector(G4LogicalVolume* world){ // Passage Matrix from Lab Referential to Telescope Referential rot = new G4RotationMatrix(u, v, w); // translation to place Telescope - pos = w * SquareLength * 0.5 + Center; + pos = w * SiliconThickness* 0.5 + Center; new G4PVPlacement(G4Transform3D(*rot, pos), BuildSquareDetector(), "MugastSquare", world, false, m_DetectorNumber[i]); } else if(m_Shape[i]=="Trapezoid"){ @@ -297,9 +308,9 @@ void Mugast::ConstructDetector(G4LogicalVolume* world){ u = m_X128_Y1[i] - m_X1_Y1[i]; u = u.unit(); - v = (m_X1_Y128[i] + m_X128_Y128[i] - m_X1_Y1[i] - m_X128_Y1[i] ); + v = (m_X1_Y128[i] + m_X128_Y128[i] - m_X1_Y1[i] - m_X128_Y1[i] ); v = v.unit(); - cout << u << " " << v << " " << u.dot(v) << endl; + w = u.cross(v); w = w.unit(); @@ -309,7 +320,7 @@ void Mugast::ConstructDetector(G4LogicalVolume* world){ rot = new G4RotationMatrix(u, v, w); rot->rotate(180*deg,w); // translation to place Telescope - pos = w * TrapezoidLength * 0.5 + Center; + pos = w * SiliconThickness* 0.5 + Center; new G4PVPlacement(G4Transform3D(*rot, pos), BuildTrapezoidDetector(), "MugastTrapezoid", world, false, m_DetectorNumber[i]); } @@ -331,23 +342,108 @@ void Mugast::InitializeRootOutput(){ // Read sensitive part and fill the Root tree. // Called at in the EventAction::EndOfEventAvtion void Mugast::ReadSensitive(const G4Event* ){ -/* m_Event->Clear(); + m_Event->Clear(); + + /////////// + // Square + DSSDScorers::PS_Rectangle* SquareScorer = (DSSDScorers::PS_Rectangle*) m_SquareScorer->GetPrimitive(0); + + + // Loop on the Square map + unsigned int sizeFront= SquareScorer->GetLengthMult(); + + for (unsigned int i=0 ; i<sizeFront ; i++){ + + double Energy = RandGauss::shoot(SquareScorer->GetEnergyLength(i), SigmaEnergy); + + if(Energy>EnergyThreshold){ + double Time = RandGauss::shoot(SquareScorer->GetTimeLength(i), SigmaTime); + int DetNbr = SquareScorer->GetDetectorLength(i); + int StripFront = SquareScorer->GetStripLength(i); + + m_Event->SetDSSDXE(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseSquareX[StripFront-1], + NPL::EnergyToADC(Energy, 0, 63, 8192, 16384)); + + m_Event->SetDSSDXT(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseSquareX[StripFront-1], + NPL::EnergyToADC(Time, 0, 1000, 8192, 16384)); + + } + } + + unsigned int sizeBack= SquareScorer->GetWidthMult(); + for (unsigned int i=0 ; i<sizeBack ; i++){ + double Energy = RandGauss::shoot(SquareScorer->GetEnergyWidth(i), SigmaEnergy); + + if(Energy>EnergyThreshold){ + double Time = RandGauss::shoot(SquareScorer->GetTimeWidth(i), SigmaTime); + int DetNbr = SquareScorer->GetDetectorWidth(i); + int StripBack = SquareScorer->GetStripWidth(i); + + m_Event->SetDSSDYE(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseSquareY[StripBack-1], + NPL::EnergyToADC(Energy, 0, 63, 8192, 0)); + + m_Event->SetDSSDYT(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseSquareY[StripBack-1], + NPL::EnergyToADC(Time, 0, 1000, 8192, 16384)); + } + } + // clear map for next event + SquareScorer->clear(); /////////// - // Calorimeter scorer - CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_MugastScorer->GetPrimitive(0); - - unsigned int size = Scorer->GetMult(); - for(unsigned int i = 0 ; i < size ; i++){ - vector<unsigned int> level = Scorer->GetLevel(i); - double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Mugast_NS::SigmaEnergy); - if(Energy>Mugast_NS::ThresholdSi){ - double Time = RandGauss::shoot(Scorer->GetTime(i),Mugast_NS::SigmaTime); - int DetectorNbr = level[0]; - // m_Event->SetEnergy(DetectorNbr,Energy); - // m_Event->SetTime(DetectorNbr,Time); + // Trapezoid + DSSDScorers::PS_Rectangle* TrapezoidScorer = (DSSDScorers::PS_Rectangle*) m_TrapezoidScorer->GetPrimitive(0); + + + // Loop on the Trapezoid map + sizeFront= TrapezoidScorer->GetLengthMult(); + + for (unsigned int i=0 ; i<sizeFront ; i++){ + + double Energy = RandGauss::shoot(TrapezoidScorer->GetEnergyLength(i), SigmaEnergy); + + if(Energy>EnergyThreshold){ + double Time = RandGauss::shoot(TrapezoidScorer->GetTimeLength(i), SigmaTime); + int DetNbr = TrapezoidScorer->GetDetectorLength(i); + int StripFront = TrapezoidScorer->GetStripLength(i); + + m_Event->SetDSSDXE(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseTrapezeX[StripFront-1], + NPL::EnergyToADC(Energy, 0, 63, 8192, 16384)); + + m_Event->SetDSSDXT(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseTrapezeX[StripFront-1], + NPL::EnergyToADC(Time, 0, 1000, 8192, 16384)); + } - }*/ + } + + sizeBack= TrapezoidScorer->GetWidthMult(); + for (unsigned int i=0 ; i<sizeBack ; i++){ + + double Energy = RandGauss::shoot(TrapezoidScorer->GetEnergyWidth(i), SigmaEnergy); + + if(Energy>EnergyThreshold){ + double Time = RandGauss::shoot(TrapezoidScorer->GetTimeWidth(i), SigmaTime); + int DetNbr = TrapezoidScorer->GetDetectorWidth(i); + int StripBack = 128-TrapezoidScorer->GetStripWidth(i)+1; + + m_Event->SetDSSDYE(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseTrapezeY[StripBack-1], + NPL::EnergyToADC(Energy, 0, 63, 8192, 0)); + + m_Event->SetDSSDYT(MG_DetectorType::MG_NOCHANGE,DetNbr, + MUGAST_MAP::ReverseTrapezeY[StripBack-1], + NPL::EnergyToADC(Time, 0, 1000, 8192, 16384)); + } + } + // clear map for next event + TrapezoidScorer->clear(); + + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -355,19 +451,41 @@ void Mugast::ReadSensitive(const G4Event* ){ void Mugast::InitializeScorers() { // This check is necessary in case the geometry is reloaded bool already_exist = false; - m_MugastScorer = CheckScorer("MugastScorer",already_exist) ; + m_SquareScorer= CheckScorer("ScorerSquare",already_exist) ; + m_TrapezoidScorer= CheckScorer("ScorerTrapezoid",already_exist) ; + m_AnnularScorer= CheckScorer("MugastScorerAnnular",already_exist) ; if(already_exist) return ; // Otherwise the scorer is initialised - vector<int> level; level.push_back(0); - G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",level, 0) ; - G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ; + G4VPrimitiveScorer* SquareScorer = new DSSDScorers::PS_Rectangle("MugastSquareScorer",1, + SquareBase, + SquareHeight, + 128, + 128); + + + G4VPrimitiveScorer* TrapezoidScorer = new DSSDScorers::PS_Rectangle("MugastTrapezoidScorer",1, + TrapezoidBaseLarge, + TrapezoidHeight, + 128, + 128); + + G4VPrimitiveScorer* InteractionS= new InteractionScorers::PS_Interactions("InteractionS",ms_InterCoord, 0) ; + G4VPrimitiveScorer* InteractionT= new InteractionScorers::PS_Interactions("InteractionT",ms_InterCoord, 0) ; + G4VPrimitiveScorer* InteractionA= new InteractionScorers::PS_Interactions("InteractionA",ms_InterCoord, 0) ; //and register it to the multifunctionnal detector - m_MugastScorer->RegisterPrimitive(Calorimeter); - m_MugastScorer->RegisterPrimitive(Interaction); - G4SDManager::GetSDMpointer()->AddNewDetector(m_MugastScorer) ; + m_SquareScorer->RegisterPrimitive(SquareScorer); + m_SquareScorer->RegisterPrimitive(InteractionS); + m_TrapezoidScorer->RegisterPrimitive(TrapezoidScorer); + m_TrapezoidScorer->RegisterPrimitive(InteractionT); +// m_AnnularScorer->RegisterPrimitive(AnnularScorer); +// m_AnnularScorer->RegisterPrimitive(Interaction); + + G4SDManager::GetSDMpointer()->AddNewDetector(m_SquareScorer) ; + G4SDManager::GetSDMpointer()->AddNewDetector(m_TrapezoidScorer) ; +// G4SDManager::GetSDMpointer()->AddNewDetector(m_MugastScorerAnnular) ; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/NPSimulation/Detectors/Mugast/Mugast.hh b/NPSimulation/Detectors/Mugast/Mugast.hh index 691101c404be7fcbc61a07c3136b75c0e756dfb5..2a77fcbe10c2fe44a119517e8fbb46dba170d9d7 100644 --- a/NPSimulation/Detectors/Mugast/Mugast.hh +++ b/NPSimulation/Detectors/Mugast/Mugast.hh @@ -91,7 +91,9 @@ class Mugast : public NPS::VDetector{ void InitializeScorers() ; // Associated Scorer - G4MultiFunctionalDetector* m_MugastScorer ; + G4MultiFunctionalDetector* m_SquareScorer; + G4MultiFunctionalDetector* m_TrapezoidScorer; + G4MultiFunctionalDetector* m_AnnularScorer; //////////////////////////////////////////////////// ///////////Event class to store Data//////////////// //////////////////////////////////////////////////// diff --git a/NPSimulation/Process/PhysicsList.cc b/NPSimulation/Process/PhysicsList.cc index 536feb997cb29cd2e6357f75cfbbc35cf38e7ffc..cd31f68fc05934cba0435b23d3082b26e84c0f6e 100644 --- a/NPSimulation/Process/PhysicsList.cc +++ b/NPSimulation/Process/PhysicsList.cc @@ -486,6 +486,7 @@ void PhysicsList::AddLevelData(){ levelData->AddPrivateData(Z, A, level_path); const G4LevelManager* levelManager = levelData->GetLevelManager(Z, A); G4int Nentries = levelManager->NumberOfTransitions()+1; + cout << " Found " << Nentries << " states including g.s." << endl; for(G4int j = 1; j < Nentries; j++){ // Excited states cout << " - Level " << j << " energy = " << levelManager->LevelEnergy(j) diff --git a/Projects/MUGAST_LISE/Analysis.cxx b/Projects/MUGAST_LISE/Analysis.cxx index 153cb7ccd01e9189f3b84e92860764f8dec7ffb6..c45ac477ce18c296c066b6a500972f5aded426df 100644 --- a/Projects/MUGAST_LISE/Analysis.cxx +++ b/Projects/MUGAST_LISE/Analysis.cxx @@ -55,15 +55,15 @@ void Analysis::Init() { string WindowsMaterial = m_DetectorManager->GetWindowsMaterial(); // energy losses - string light=NPL::ChangeNameToG4Standard(myReaction.GetNucleus3().GetName()); - string beam=NPL::ChangeNameToG4Standard(myReaction.GetNucleus1().GetName()); + string light=NPL::ChangeNameToG4Standard(myReaction.GetNucleus3()->GetName()); + string beam=NPL::ChangeNameToG4Standard(myReaction.GetNucleus1()->GetName()); LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 ); LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100); LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100); BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100); - if(WindowsThickness){ +/* if(WindowsThickness){ BeamWindow= new NPL::EnergyLoss(beam+"_"+WindowsMaterial+".G4table","G4Table",100); LightWindow= new NPL::EnergyLoss(light+"_"+WindowsMaterial+".G4table","G4Table",100); } @@ -72,7 +72,7 @@ void Analysis::Init() { BeamWindow= NULL; LightWindow=NULL; } - +*/ // initialize various parameters Rand = TRandom3(); DetectorNumber = 0; @@ -100,7 +100,7 @@ void Analysis::TreatEvent() { BeamImpact = TVector3(0,0,zImpact); // determine beam energy for a randomized interaction point in target // 1% FWHM randominastion (E/100)/2.35 - myReaction.SetBeamEnergy(Rand.Gaus(myInit->GetIncidentFinalKineticEnergy(),myInit->GetIncidentFinalKineticEnergy()/235)); + // myReaction.SetBeamEnergy(Rand.Gaus(myInit->GetIncidentFinalKineticEnergy(),myInit->GetIncidentFinalKineticEnergy()/235)); OriginalThetaLab = myInit->GetThetaLab_WorldFrame(0); OriginalELab = myInit->GetKineticEnergy(0); @@ -118,8 +118,8 @@ void Analysis::TreatEvent() { X = GD -> GetPositionOfInteraction().X(); Y = GD -> GetPositionOfInteraction().Y(); Z = GD -> GetPositionOfInteraction().Z(); - cout << " -- " << endl; - cout << X-Coord->GetDetectedPositionX(0) << " " << Y-Coord->GetDetectedPositionY(0) << " " << Z << endl; + //cout << " -- " << endl; + //cout << X-Coord->GetDetectedPositionX(0) << " " << Y-Coord->GetDetectedPositionY(0) << " " << Z << endl; ThetaGDSurface = HitDirection.Angle( TVector3(0,0,1) ) ; ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;