diff --git a/Examples/Example1/ShowResults.C b/Examples/Example1/ShowResults.C index 6c2d5867b0550fafc5bc08cc53d8a0a17037d5cb..4a6f99d23035e0e6c11d8635981fcf1c4d2dcaa4 100644 --- a/Examples/Example1/ShowResults.C +++ b/Examples/Example1/ShowResults.C @@ -42,7 +42,7 @@ ETOF->Draw("same"); // Kinematical Line // c1->cd(3); //chain->Draw("ELab:ThetaLab>>hKine(500,0,45,400,0,40)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); -chain->Draw("ELab:ThetaLab","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); +chain->Draw("ELab:ThetaLab>>h(1000,0,90,1000,0,30)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); NPL::Reaction r("11Li(d,3He)10He@553"); r.SetExcitationHeavy(1.4); diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source index e7453e1420abd32340d849183716ec78d1bad10e..6cb81d2825f67c2a17088cd685af319ea22cd7a2 100644 --- a/Inputs/EventGenerator/alpha.source +++ b/Inputs/EventGenerator/alpha.source @@ -4,10 +4,10 @@ % Energy are given in MeV , Position in mm % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic - EnergyLow= 0 MeV - EnergyHigh= 850 MeV + EnergyLow= 30 MeV + EnergyHigh= 30 MeV HalfOpenAngleMin= 0 deg - HalfOpenAngleMax= 10 deg + HalfOpenAngleMax= 180 deg x0= 0 mm y0= 0 mm z0= 0 mm diff --git a/NPLib/Core/CMakeLists.txt b/NPLib/Core/CMakeLists.txt index e5ef4df3709217bd4f734fe5ddeab98312e3d072..fc38abc16db7b1326ce3af8ad9ad8859fd05ba52 100644 --- a/NPLib/Core/CMakeLists.txt +++ b/NPLib/Core/CMakeLists.txt @@ -1,5 +1,5 @@ add_custom_command(OUTPUT TAsciiFileDict.cxx COMMAND ../scripts/build_dict.sh TAsciiFile.h TAsciiFileDict.cxx TAsciiFile.rootmap libNPCore.so) add_custom_command(OUTPUT NPVDetectorDict.cxx COMMAND ../scripts/build_dict.sh NPVDetector.h NPVDetectorDict.cxx NPVDetector.rootmap libNPCore.so NPCoreLinkdef.h) -add_library(NPCore SHARED NPVAnalysis.cxx NPAnalysisFactory.cxx NPCalibrationManager.cxx NPOptionManager.cxx RootOutput.cxx RootInput.cxx TAsciiFile.cxx TAsciiFileDict.cxx NPDetectorManager.cxx NPVDetector.cxx NPVDetectorDict.cxx NPVSpectra.cxx NPDetectorFactory.cxx NPSpectraServer.cxx NPInputParser.cxx) +add_library(NPCore SHARED NPVAnalysis.cxx NPAnalysisFactory.cxx NPCalibrationManager.cxx NPOptionManager.cxx RootOutput.cxx RootInput.cxx TAsciiFile.cxx TAsciiFileDict.cxx NPDetectorManager.cxx NPVDetector.cxx NPVDetectorDict.cxx NPVSpectra.cxx NPDetectorFactory.cxx NPSpectraServer.cxx NPInputParser.cxx NPImage.cxx) target_link_libraries(NPCore ${ROOT_LIBRARIES}) -install(FILES NPVAnalysis.h NPAnalysisFactory.h NPCalibrationManager.h NPOptionManager.h RootInput.h RootOutput.h TAsciiFile.h NPDetectorManager.h NPVDetector.h NPGlobalSystemOfUnits.h NPPhysicalConstants.h NPSystemOfUnits.h NPVSpectra.h NPDetectorFactory.h NPSpectraServer.h NPInputParser.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) +install(FILES NPVAnalysis.h NPAnalysisFactory.h NPCalibrationManager.h NPOptionManager.h RootInput.h RootOutput.h TAsciiFile.h NPDetectorManager.h NPVDetector.h NPGlobalSystemOfUnits.h NPPhysicalConstants.h NPSystemOfUnits.h NPVSpectra.h NPDetectorFactory.h NPSpectraServer.h NPInputParser.h NPImage.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) diff --git a/NPLib/Core/NPImage.cxx b/NPLib/Core/NPImage.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b14777de0581133dfc6e667e49670326781b4564 --- /dev/null +++ b/NPLib/Core/NPImage.cxx @@ -0,0 +1,159 @@ +/***************************************************************************** + * 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: matta@lpccaen.in2p3.fr * + * * + * Creation Date : * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class is wrapper of root TASImage to manipulate and generate * + * detector pixel map * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +//NPTool +#include "NPImage.h" + +//STL +#include<iostream> + +//////////////////////////////////////////////////////////////////////////////// +NPL::Image::Image(){ + m_Image = NULL; + m_ARGB = NULL; +} +//////////////////////////////////////////////////////////////////////////////// +NPL::Image::Image(std::string filename,double Xscalling,double Yscalling){ + m_Image = new TASImage(filename.c_str()); + m_ARGB = m_Image->GetArgbArray(); + m_Xscaling = Xscalling; + m_Yscaling = Yscalling; +} +//////////////////////////////////////////////////////////////////////////////// +NPL::Image::~Image(){ + if(m_Image) + delete m_Image; + if(m_ARGB) + delete m_ARGB; +} +//////////////////////////////////////////////////////////////////////////////// +void NPL::Image::Open(std::string filename){ + if(m_Image) + delete m_Image; + if(m_ARGB) + delete m_ARGB; + + m_Image = new TASImage(filename.c_str()); + m_ARGB = m_Image->GetArgbArray(); +} + +//////////////////////////////////////////////////////////////////////////////// +void NPL::Image::Print(double scale){ + unsigned int yPixel = m_Image->GetHeight(); + unsigned int xPixel = m_Image->GetWidth(); + + m_Image->Scale(scale*xPixel,scale*yPixel); + + for(unsigned int i = 0 ; i < m_Image->GetWidth() ; i++){ + for(unsigned int j = 0 ; j < m_Image->GetHeight() ; j++){ + unsigned int index = i*xPixel+j; + unsigned int b = GetBlueAtPixel(i,j); + unsigned int g = GetGreenAtPixel(i,j); + unsigned int r = GetRedAtPixel(i,j); + unsigned int a = GetAlphaAtPixel(i,j); + + if(r||b||g) + std::cout << "x"<< " " ; + else + std::cout << "."<<" " ; + } + std::cout << std::endl; + } +} +//////////////////////////////////////////////////////////////////////////////// +void NPL::Image::Save(std::string filename){ + m_Image->WriteImage(filename.c_str()); +} + + +//////////////////////////////////////////////////////////////////////////////// +void NPL::Image::Draw(){ + m_Image->Draw(); +} + +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetRedAtCoordinate(double x, double y){ + unsigned int xp = x/m_Xscaling+m_Image->GetWidth()/2 ; + unsigned int yp = y/m_Yscaling+m_Image->GetHeight()/2 ; + return GetRedAtPixel(xp,yp); +} + +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetGreenAtCoordinate(double x, double y){ + unsigned int xp = x/m_Xscaling+m_Image->GetWidth()/2 ; + unsigned int yp = y/m_Yscaling+m_Image->GetHeight()/2 ; + return GetRedAtPixel(xp,yp); +} +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetBlueAtCoordinate(double x, double y){ + unsigned int xp = x/m_Xscaling+m_Image->GetWidth()/2 ; + unsigned int yp = y/m_Yscaling+m_Image->GetHeight()/2 ; + return GetRedAtPixel(xp,yp); +} +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetAlphaAtCoordinate(double x, double y){ + unsigned int xp = x/m_Xscaling+m_Image->GetWidth()/2 ; + unsigned int yp = y/m_Yscaling+m_Image->GetHeight()/2 ; + return GetRedAtPixel(xp,yp); +} +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetPixelAtCoordinate(double x, double y){ + unsigned int xp = x/m_Xscaling+m_Image->GetWidth()/2 ; + unsigned int yp = y/m_Yscaling+m_Image->GetHeight()/2 ; + return m_ARGB[PixelToIndex(xp,yp)]; +} + + +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::PixelToIndex(unsigned int x, unsigned int y){ + return x*m_Image->GetWidth()+y; +} + +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetRedAtPixel(unsigned int x, unsigned int y){ + unsigned int v = m_ARGB[PixelToIndex(x,y)]; + v= (v>>16)&0xff; + return v; +} +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetGreenAtPixel(unsigned int x, unsigned int y){ + unsigned int v = m_ARGB[PixelToIndex(x,y)]; + v= (v>>8)&0xff; + return v; + +} +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetBlueAtPixel(unsigned int x, unsigned int y){ + unsigned int v = m_ARGB[PixelToIndex(x,y)]; + v= v&0xff; + return v; + +} +//////////////////////////////////////////////////////////////////////////////// +unsigned int NPL::Image::GetAlphaAtPixel(unsigned int x, unsigned int y){ + unsigned int v = m_ARGB[PixelToIndex(x,y)]; + v= (v>>24)&0xff; + return v; + +} + diff --git a/NPLib/Core/NPImage.h b/NPLib/Core/NPImage.h new file mode 100644 index 0000000000000000000000000000000000000000..39000431e42849d1462f03b0393f5e85150c2c62 --- /dev/null +++ b/NPLib/Core/NPImage.h @@ -0,0 +1,71 @@ +/***************************************************************************** + * 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: matta@lpccaen.in2p3.fr * + * * + * Creation Date : * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class is wrapper of root TASImage to manipulate and generate * + * detector pixel map * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + + +// ROOT +#include"TASImage.h" + +//STL +#include<string> + +namespace NPL{ + + class Image{ + + public: + Image(); + Image(std::string filename,double Xscalling,double Yscalling); + ~Image(); + + public: + void Open(std::string filename); + void Print(double scale=0); // print mockup in terminal + void Draw(); // draw in a new canvas (interactive only) + void Save(std::string filename); // save the image + + public: // Get the value based on coordinate + unsigned int GetRedAtCoordinate(double x, double y); + unsigned int GetGreenAtCoordinate(double x, double y); + unsigned int GetBlueAtCoordinate(double x, double y); + unsigned int GetAlphaAtCoordinate(double x, double y); + unsigned int GetPixelAtCoordinate(double x, double y); + + + public: // Get the value based on pixel + unsigned int PixelToIndex(unsigned int x, unsigned int y); + unsigned int GetRedAtPixel(unsigned int x, unsigned int y); + unsigned int GetGreenAtPixel(unsigned int x, unsigned int y); + unsigned int GetBlueAtPixel(unsigned int x, unsigned int y); + unsigned int GetAlphaAtPixel(unsigned int x, unsigned int y); + + public: + void GenerateStrip(unsigned int Nbr, double x, double y , double interstrip , unsigned int pixel); + + private: + TASImage* m_Image; // the image itself + unsigned int* m_ARGB; // the alpha red green blue array of pixel value + double m_Xscaling; // pixel per mm + double m_Yscaling; // pixel per mm + }; + +} diff --git a/NPLib/Detectors/Chio/TChio_digPhysics.cxx b/NPLib/Detectors/Chio/TChio_digPhysics.cxx index dc23e116a4c2bdb980fa0e09718e4d8535880774..7b238a1db016e3c3f5fa17649e1ad324a5e2ce31 100644 --- a/NPLib/Detectors/Chio/TChio_digPhysics.cxx +++ b/NPLib/Detectors/Chio/TChio_digPhysics.cxx @@ -62,14 +62,11 @@ void TChio_digPhysics::ReadConfiguration(NPL::InputParser parser){ if(NPOptionManager::getInstance()->GetVerboseLevel()) cout << "//// " << blocks.size() << " detectors found " << endl; - vector<string> token = {"A","B","C","D"}; + vector<string> token = {"POS"}; for(unsigned int i = 0 ; i < blocks.size() ; i++){ if(blocks[i]->HasTokenList(token)){ - TVector3 A = blocks[i]->GetTVector3("A","mm"); - TVector3 B = blocks[i]->GetTVector3("B","mm"); - TVector3 C = blocks[i]->GetTVector3("C","mm"); - TVector3 D = blocks[i]->GetTVector3("D","mm"); + TVector3 Pos = blocks[i]->GetTVector3("POS","mm"); // AddChio(A,B,C,D); } diff --git a/NPLib/Detectors/MUST2/ressources/Mask.cxx b/NPLib/Detectors/MUST2/ressources/Mask.cxx new file mode 100644 index 0000000000000000000000000000000000000000..33c72bc1d99ea77f2c99759b6319b781f95eed92 --- /dev/null +++ b/NPLib/Detectors/MUST2/ressources/Mask.cxx @@ -0,0 +1,58 @@ +void Mask(){ + + double dimension = 100.42; + double active = 97.30; + double pitch = 97.30/128; + double width = 0.7; +// +// double dimension = 100; +// double active = 97; +// double pitch = 97/97; +// double width = 0.5; + + // mm per pixel + double scale = 0.01; + //pitch in pixel + unsigned int spitch = pitch/scale; + unsigned int swidth = width/scale; + unsigned int sinter = (spitch - swidth)/2; + cout << spitch << " " << swidth << " " << sinter << endl; + // image size + unsigned int size = dimension/scale; + cout << "Image size: " << size << "x" << size<< endl ; + double* zargb = new double[size*size]; + TASImage* zero = new TASImage("zero",zargb,size,size,0); + zero->WriteImage("mask.png"); + TASImage* mask = new TASImage("mask.png"); + delete[] zargb; + unsigned int* argb = mask->GetArgbArray(); + unsigned int* argb2 = mask->GetArgbArray(); + unsigned int index = 0; + double border1 = (dimension-active)/scale; + double border2 = (active)/scale; + + for(unsigned int px = 0 ; px < size ; px++){ + for(unsigned int py = 0 ; py < size ; py++){ + if(px%1000==0) + cout << "\r" << px << "/" << size << flush; + // Compute array index + index = px * size + py; + // Inactive sides + if(px < border1|| py < border1 || px>border2 || py>border2) + argb[index] = 0xffff0000; + else{ // strips + unsigned int coord = px-border1; + unsigned int nbr = coord/spitch; + // cout << coord << " " << nbr*spitch+sinter << " " << spitch << " " << sinter << " " << spitch << " " << sinter << endl; + if(coord<(nbr*spitch+sinter) || coord >(((nbr+1)*spitch-sinter))){// interstrip + argb[index] = 0xffff0000+(((nbr+1))<<8)+nbr; + } + else + argb[index] = 0xff000000 + nbr; + } + } + } + mask->WriteImage("mask.png"); + //mask->Draw(); + +} diff --git a/NPLib/Detectors/MUST2/ressources/maskBack.png b/NPLib/Detectors/MUST2/ressources/maskBack.png new file mode 100644 index 0000000000000000000000000000000000000000..40753f6d6e226ee99d1d9a27524e1576f56a9a5f Binary files /dev/null and b/NPLib/Detectors/MUST2/ressources/maskBack.png differ diff --git a/NPLib/Detectors/MUST2/ressources/maskFront.png b/NPLib/Detectors/MUST2/ressources/maskFront.png new file mode 100644 index 0000000000000000000000000000000000000000..6644c19a7279d326859c83ab956ac8ffc3fa3681 Binary files /dev/null and b/NPLib/Detectors/MUST2/ressources/maskFront.png differ diff --git a/NPLib/ressources/CMake/Root.cmake b/NPLib/ressources/CMake/Root.cmake index 245789a62d628452a751dd78ce3620b8b62fb545..d6b72cb0d510032846e6c2e8903c16c1cf8f3241 100644 --- a/NPLib/ressources/CMake/Root.cmake +++ b/NPLib/ressources/CMake/Root.cmake @@ -36,7 +36,7 @@ endif() ## Collect the different information ## # List of Root dependencies -set(ROOT_LIBRARIES dl Gui Core RIO Net Hist Gpad Tree Physics MathCore Thread ) +set(ROOT_LIBRARIES dl Gui Core RIO Net Hist Gpad Tree Physics MathCore Thread ASImage) # Lib directories exec_program(${NPTOOL_ROOT_CONFIG} ARGS "--libdir" OUTPUT_VARIABLE ROOT_LIBRARY_DIR) diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc index bfe806930eedede351090fc25b3820d4a70410b6..52625d30d7b981b38bc4df9d12fbca1fae13277d 100644 --- a/NPSimulation/Detectors/MUST2/MUST2Array.cc +++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc @@ -37,13 +37,13 @@ #include "G4VisAttributes.hh" #include "G4Colour.hh" #include "G4PVDivision.hh" - +#include "Randomize.hh" // NPS #include "MaterialManager.hh" #include "NPSDetectorFactory.hh" #include "ObsoleteGeneralScorers.hh" -#include "MUST2Scorers.hh" - +#include "SiliconScorers.hh" +#include "CalorimeterScorers.hh" #include "NPOptionManager.h" //ROOT @@ -161,7 +161,7 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, Name , world , false , - 0 ); + TelescopeNumber ); if (m_non_sensitive_part_visiualisation){ G4VisAttributes* FrameVisAtt = new G4VisAttributes(G4Colour(0.80, 0.80, 0.80)); @@ -176,7 +176,7 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0); - new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, 0); + new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, TelescopeNumber ); logicVacBox->SetVisAttributes(G4VisAttributes::Invisible); @@ -193,8 +193,8 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, G4Box* solidAluStrip = new G4Box("AluBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*AluStripThickness); G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0); - new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicMM, false, 0); - new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicMM, false, 0); + new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicMM, false, TelescopeNumber ); + new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicMM, false, TelescopeNumber ); logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible); @@ -203,7 +203,7 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, G4Box* solidSilicon = new G4Box("solidSilicon", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconThickness); G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0); - new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicMM, false, 0); + new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicMM, false, TelescopeNumber ); ///Set Silicon strip sensible @@ -245,148 +245,158 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, G4ThreeVector ShiftSiLiDown = G4ThreeVector(0.25 * SiliconFace + 0.5 * SiLiSpace, 0, 0) ; // SiLi : left side of SiLi detector - G4Box* solidSiLi_LT = new G4Box("SiLi_LT" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); - G4Box* solidSiLi_RT = new G4Box("SiLi_RT" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); - G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); - G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); - G4Box* solidSiLi_LB = new G4Box("SiLi_LB" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); - G4Box* solidSiLi_RB = new G4Box("SiLi_RB" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); - G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); - G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); - - G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume(solidSiLi_LT , m_MaterialSilicon , "SiLi_LT" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume(solidSiLi_RT , m_MaterialSilicon , "SiLi_RT" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume(solidSiLi_LC1 , m_MaterialSilicon , "SiLi_LC1" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume(solidSiLi_RC1 , m_MaterialSilicon , "SiLi_RC1" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume(solidSiLi_LB , m_MaterialSilicon , "SiLi_LB" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume(solidSiLi_RB , m_MaterialSilicon , "SiLi_RB" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume(solidSiLi_LC2 , m_MaterialSilicon , "SiLi_LC2" , 0 , 0 , 0); - G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume(solidSiLi_RC2 , m_MaterialSilicon , "SiLi_RC2" , 0 , 0 , 0); + G4Box* solidSiLi_LT = new G4Box("SiLi_LT", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); + G4Box* solidSiLi_RT = new G4Box("SiLi_RT", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); + G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); + G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); + G4Box* solidSiLi_LB = new G4Box("SiLi_LB", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); + G4Box* solidSiLi_RB = new G4Box("SiLi_RB", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); + G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); + G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); + + G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume(solidSiLi_LT , m_MaterialSilicon, "SiLi_LT" , 0, 0, 0); + G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume(solidSiLi_RT , m_MaterialSilicon, "SiLi_RT" , 0, 0, 0); + G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume(solidSiLi_LC1, m_MaterialSilicon, "SiLi_LC1", 0, 0, 0); + G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume(solidSiLi_RC1, m_MaterialSilicon, "SiLi_RC1", 0, 0, 0); + G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume(solidSiLi_LB , m_MaterialSilicon, "SiLi_LB" , 0, 0, 0); + G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume(solidSiLi_RB , m_MaterialSilicon, "SiLi_RB" , 0, 0, 0); + G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume(solidSiLi_LC2, m_MaterialSilicon, "SiLi_LC2", 0, 0, 0); + G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume(solidSiLi_RC2, m_MaterialSilicon, "SiLi_RC2", 0, 0, 0); G4double interSiLi = 0.5 * mm; // Top - G4ThreeVector positionSiLi_LT_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LT_up = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, + 0); positionSiLi_LT_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_RT_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RT_up = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, + 0); positionSiLi_RT_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_LC1_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LC1_up = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, + 0); positionSiLi_LC1_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_RC1_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RC1_up = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, + 0); positionSiLi_RC1_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_LB_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , - 0); + G4ThreeVector positionSiLi_LB_up = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , + 0); positionSiLi_LB_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_RB_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX , - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , - 0); + G4ThreeVector positionSiLi_RB_up = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX , + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , + 0); positionSiLi_RB_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_LC2_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LC2_up = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, + 0); positionSiLi_LC2_up += ShiftSiLiUp ; - G4ThreeVector positionSiLi_RC2_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX , - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi , - 0); + G4ThreeVector positionSiLi_RC2_up = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX , + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi , + 0); positionSiLi_RC2_up += ShiftSiLiUp ; // Down - G4ThreeVector positionSiLi_LT_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LT_down = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, + 0); positionSiLi_LT_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_RT_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RT_down = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, + 0); positionSiLi_RT_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_LC1_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LC1_down = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, + 0); positionSiLi_LC1_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_RC1_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi , - 0); + G4ThreeVector positionSiLi_RC1_down = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi , + 0); positionSiLi_RC1_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_LB_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LB_down = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, + 0); positionSiLi_LB_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_RB_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RB_down = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, + 0); positionSiLi_RB_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_LC2_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LC2_down = + G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, + 0); positionSiLi_LC2_down += ShiftSiLiDown ; - G4ThreeVector positionSiLi_RC2_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RC2_down = + G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, + 0); positionSiLi_RC2_down += ShiftSiLiDown ; + new G4PVPlacement(0 , positionSiLi_RT_down , logicSiLi_RT , Name + "_SiLi_Pad1" , logicSiLi , false ,1); + new G4PVPlacement(0 , positionSiLi_LT_down , logicSiLi_LT , Name + "_SiLi_Pad2" , logicSiLi , false ,2); + new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" , logicSiLi , false ,3); + new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" , logicSiLi , false ,4); + new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" , logicSiLi , false ,5); + new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" , logicSiLi , false ,6); + new G4PVPlacement(0 , positionSiLi_LB_down , logicSiLi_LB , Name + "_SiLi_Pad7" , logicSiLi , false ,7); + new G4PVPlacement(0 , positionSiLi_RB_down , logicSiLi_RB , Name + "_SiLi_Pad8" , logicSiLi , false ,8); + new G4PVPlacement(0 , positionSiLi_LT_up , logicSiLi_LT , Name + "_SiLi_Pad9" , logicSiLi , false ,9); + new G4PVPlacement(0 , positionSiLi_RT_up , logicSiLi_RT , Name + "_SiLi_Pad10" , logicSiLi , false ,10); + new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" , logicSiLi , false ,11); + new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" , logicSiLi , false ,12); + new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" , logicSiLi , false ,13); + new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" , logicSiLi , false ,14); + new G4PVPlacement(0 , positionSiLi_RB_up , logicSiLi_RB , Name + "_SiLi_Pad15" , logicSiLi , false ,15); + new G4PVPlacement(0 , positionSiLi_LB_up , logicSiLi_LB , Name + "_SiLi_Pad16" , logicSiLi , false ,16); - // up - new G4PVPlacement(0 , positionSiLi_LT_up , logicSiLi_LT , Name + "_SiLi_Pad9" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_RT_up , logicSiLi_RT , Name + "_SiLi_Pad10" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" , logicSiLi , false , 0); - new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" , logicSiLi , false , 0); - - new G4PVPlacement(0 , positionSiLi_LB_up , logicSiLi_LB , Name + "_SiLi_Pad16" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_RB_up , logicSiLi_RB , Name + "_SiLi_Pad15" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" , logicSiLi , false , 0); - new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" , logicSiLi , false , 0); - - - // down - new G4PVPlacement(0 , positionSiLi_LT_down , logicSiLi_LT , Name + "_SiLi_Pad2" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_RT_down , logicSiLi_RT , Name + "_SiLi_Pad1" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" , logicSiLi , false , 0) ; - - new G4PVPlacement(0 , positionSiLi_LB_down , logicSiLi_LB , Name + "_SiLi_Pad7" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_RB_down , logicSiLi_RB , Name + "_SiLi_Pad8" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" , logicSiLi , false , 0) ; - new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" , logicSiLi , false , 0) ; // Set SiLi sensible logicSiLi_LT->SetSensitiveDetector(m_SiLiScorer); @@ -563,13 +573,13 @@ void MUST2Array::ReadConfiguration(NPL::InputParser parser){ cout << "//// " << blocks.size() << " telescope found" << endl; for(unsigned int i = 0 ; i < blocks.size() ; i++){ if(NPOptionManager::getInstance()->GetVerboseLevel()) - cout << endl << "//// Must 2 Telecope " << i+1 << endl; + cout << endl << "//// Must 2 Telecope " << i+1 << endl; // Cartesian Case vector<string> cart = {"X1_Y1","X1_Y128","X128_Y1","X128_Y128","SI","SILI","CSI"}; // Spherical Case vector<string> sphe= {"R","THETA","PHI","BETA","SI","SILI","CSI"}; - if(blocks[i]->HasTokenList(cart)){ + if(blocks[i]->HasTokenList(cart)){ G4ThreeVector A = NPS::ConvertVector( blocks[i]->GetTVector3("X1_Y1","mm")); G4ThreeVector B = NPS::ConvertVector( blocks[i]->GetTVector3("X128_Y1","mm")); G4ThreeVector C = NPS::ConvertVector( blocks[i]->GetTVector3("X1_Y128","mm")); @@ -598,7 +608,7 @@ void MUST2Array::ReadConfiguration(NPL::InputParser parser){ } if(blocks[i]->GetString("VIS")=="all") - m_non_sensitive_part_visiualisation = true; + m_non_sensitive_part_visiualisation = true; } } @@ -609,13 +619,13 @@ void MUST2Array::ReadConfiguration(NPL::InputParser parser){ void MUST2Array::ConstructDetector(G4LogicalVolume* world){ G4RotationMatrix* MMrot = NULL ; G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ; G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ; - bool Si = true ; - bool SiLi = true ; - bool CsI = true ; + bool Si = true ; + bool SiLi = true ; + bool CsI = true ; G4int NumberOfTelescope = m_DefinitionType.size() ; @@ -708,320 +718,184 @@ void MUST2Array::InitializeRootOutput(){ // Read sensitive part and fill the Root tree. // Called at in the EventAction::EndOfEventAvtion void MUST2Array::ReadSensitive(const G4Event* event){ - G4String DetectorNumber ; - m_Event->Clear() ; + G4String DetectorNumber; + m_Event->Clear(); ////////////////////////////////////////////////////////////////////////////////////// //////////////////////// Used to Read Event Map of detector ////////////////////////// ////////////////////////////////////////////////////////////////////////////////////// - // Si - std::map<G4int, G4int*>::iterator DetectorNumber_itr; - std::map<G4int, G4double*>::iterator Energy_itr; - std::map<G4int, G4double*>::iterator Time_itr; - std::map<G4int, G4int*>::iterator X_itr; - std::map<G4int, G4int*>::iterator Y_itr; - std::map<G4int, G4double*>::iterator Pos_X_itr; - std::map<G4int, G4double*>::iterator Pos_Y_itr; - std::map<G4int, G4double*>::iterator Pos_Z_itr; - std::map<G4int, G4double*>::iterator Ang_Theta_itr; - std::map<G4int, G4double*>::iterator Ang_Phi_itr; - - NPS::HitsMap<G4int>* DetectorNumberHitMap; - NPS::HitsMap<G4double>* EnergyHitMap; - NPS::HitsMap<G4double>* TimeHitMap; - NPS::HitsMap<G4int>* XHitMap; - NPS::HitsMap<G4int>* YHitMap; - NPS::HitsMap<G4double>* PosXHitMap; - NPS::HitsMap<G4double>* PosYHitMap; - NPS::HitsMap<G4double>* PosZHitMap; - NPS::HitsMap<G4double>* AngThetaHitMap; - NPS::HitsMap<G4double>* AngPhiHitMap; - - // Si(Li) - std::map<G4int, G4double*>::iterator SiLiEnergy_itr; - std::map<G4int, G4int*>::iterator SiLiPadNbr_itr; - NPS::HitsMap<G4double>* SiLiEnergyHitMap; - NPS::HitsMap<G4int>* SiLiPadNbrHitMap; - - // CsI - std::map<G4int, G4double*>::iterator CsIEnergy_itr; - std::map<G4int, G4int*>::iterator CsICristalNbr_itr; - NPS::HitsMap<G4double>* CsIEnergyHitMap; - NPS::HitsMap<G4int>* CsICristalNbrHitMap; - ////////////////////////////////////////////////////////////////////////////////////// - ////////////////////////////////////////////////////////////////////////////////////// - - // Read the Scorer associate to the Silicon Strip - - //Detector Number - G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/DetectorNumber"); - DetectorNumberHitMap = (NPS::HitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)); - DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin(); - - //Energy - G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripEnergy") ; - EnergyHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ; - Energy_itr = EnergyHitMap->GetMap()->begin(); - - //Time of Flight - G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripTime"); - TimeHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)); - Time_itr = TimeHitMap->GetMap()->begin(); - - //Strip Number X - G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberX") ; - XHitMap = (NPS::HitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)); - X_itr = XHitMap->GetMap()->begin(); - - //Strip Number Y - G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberY"); - YHitMap = (NPS::HitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)); - Y_itr = YHitMap->GetMap()->begin(); - - //Interaction Coordinate X - G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordX"); - PosXHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)); - Pos_X_itr = PosXHitMap->GetMap()->begin(); - - //Interaction Coordinate Y - G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordY"); - PosYHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)); - Pos_Y_itr = PosYHitMap->GetMap()->begin(); - - //Interaction Coordinate Z - G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordZ"); - PosZHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)); - Pos_Z_itr = PosXHitMap->GetMap()->begin(); - - //Interaction Coordinate Angle Theta - G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngTheta"); - AngThetaHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)); - Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); - - //Interaction Coordinate Angle Phi - G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngPhi"); - AngPhiHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)); - Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); - - // Read the Scorer associate to the SiLi - //Energy - G4int SiLiEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiEnergy"); - SiLiEnergyHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiEnergyCollectionID)); - SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin(); - - G4int SiLiPadNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiPadNbr"); - SiLiPadNbrHitMap = (NPS::HitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(SiLiPadNbrCollectionID)); - SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin(); - - // Read the Scorer associate to the CsI crystal - //Energy - G4int CsIEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIEnergy"); - CsIEnergyHitMap = (NPS::HitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIEnergyCollectionID)); - CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin(); - - G4int CsICristalNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsICristalNbr"); - CsICristalNbrHitMap = (NPS::HitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICristalNbrCollectionID)); - CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin(); - - ///////////////////// - - G4int sizeN = DetectorNumberHitMap->entries(); - G4int sizeE = EnergyHitMap->entries(); - G4int sizeT = TimeHitMap->entries(); - G4int sizeX = XHitMap->entries(); - G4int sizeY = YHitMap->entries(); - - // Loop on Telescope Number entry - for (G4int l = 0 ; l < sizeN ; l++) { - G4int N = *(DetectorNumber_itr->second); - G4int NTrackID = DetectorNumber_itr->first - N; - - - if (N > 0) { - - m_Event->SetMMStripXEDetectorNbr(N) ; - m_Event->SetMMStripYEDetectorNbr(N) ; - m_Event->SetMMStripXTDetectorNbr(N) ; - m_Event->SetMMStripYTDetectorNbr(N) ; - - // Energy - Energy_itr = EnergyHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeE ; h++) { - G4int ETrackID = Energy_itr->first - N; - G4double E = *(Energy_itr->second); - - if (ETrackID == NTrackID) { - m_Event->SetMMStripXEEnergy(RandGauss::shoot(E, ResoStrip)); - m_Event->SetMMStripYEEnergy(RandGauss::shoot(E, ResoStrip)); - } - - Energy_itr++; - } - - - // Time - Time_itr = TimeHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeT ; h++) { - G4int TTrackID = Time_itr->first - N; - G4double T = *(Time_itr->second); - - if (TTrackID == NTrackID) { - m_Event->SetMMStripXTTime(RandGauss::shoot(T, ResoTimeMust)) ; - m_Event->SetMMStripYTTime(RandGauss::shoot(T, ResoTimeMust)) ; - } - - Time_itr++; - } - - - // X - X_itr = XHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { - G4int XTrackID = X_itr->first - N; - G4int X = *(X_itr->second); - if (XTrackID == NTrackID) { - m_Event->SetMMStripXEStripNbr(X); - m_Event->SetMMStripXTStripNbr(X); - } - - X_itr++; + // Read the Scorer associate to the Silicon Strip + SILICONSCORERS::PS_Silicon_Images* SiScorer = (SILICONSCORERS::PS_Silicon_Images*) m_StripScorer->GetPrimitive(0); + + set<int> trig; // list of telescope that got a Si trigger + vector<unsigned int> indexes = SiScorer->GetIndexes(); + unsigned int size = indexes.size(); + for(unsigned int i = 0 ;i<size;i++){ + unsigned int index = indexes[i]; + double energy = SiScorer->GetEnergy(index); + double energyX = RandGauss::shoot(energy, ResoStrip); + double energyY= RandGauss::shoot(energy, ResoStrip); + int detectorNbr = SiScorer->GetDetectorNbr(index); + double time = SiScorer->GetTime(index); + // X + if(energyX>0.1*keV){ // above threshold + // Pixel value at interaction point + unsigned int a,r,g,b; + // pixel + SiScorer->GetARGBFront(indexes[i],a,r,g,b); + if(r==0){ + trig.insert(detectorNbr); + // Energy + m_Event->SetMMStripXEDetectorNbr(detectorNbr); + m_Event->SetMMStripXEStripNbr(b+1); + m_Event->SetMMStripXEEnergy(energyX); + + // Time + m_Event->SetMMStripXTDetectorNbr(detectorNbr); + m_Event->SetMMStripXTStripNbr(b+1); + m_Event->SetMMStripXTTime(RandGauss::shoot(time, ResoTimeMust)); } - - // Y - Y_itr = YHitMap->GetMap()->begin() ; - for (G4int h = 0 ; h < sizeY ; h++) { - G4int YTrackID = Y_itr->first - N ; - G4int Y = *(Y_itr->second); - if (YTrackID == NTrackID) { - m_Event->SetMMStripYEStripNbr(Y); - m_Event->SetMMStripYTStripNbr(Y); + else{ // Interstrip X, keep maximum shared energy + double rand = G4UniformRand(); + if(rand>0.5){ + energyX = rand*energyX; + if(energyX>0.1*keV){ + trig.insert(detectorNbr); + // Energy + m_Event->SetMMStripXEDetectorNbr(detectorNbr); + m_Event->SetMMStripXEStripNbr(b+1); + m_Event->SetMMStripXEEnergy(energyX); + + // Time + m_Event->SetMMStripXTDetectorNbr(detectorNbr); + m_Event->SetMMStripXTStripNbr(b+1); + m_Event->SetMMStripXTTime(RandGauss::shoot(time, ResoTimeMust)); + } } - - Y_itr++; - } - - // Pos X - Pos_X_itr = PosXHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < PosXHitMap->entries() ; h++) { - G4int PosXTrackID = Pos_X_itr->first - N ; - G4double PosX = *(Pos_X_itr->second) ; - if (PosXTrackID == NTrackID) { - ms_InterCoord->SetDetectedPositionX(PosX) ; + else{ + energyX = (1-rand)*energyX; + if(energyX>0.1*keV){ + trig.insert(detectorNbr); + + // Energy + m_Event->SetMMStripXEDetectorNbr(detectorNbr); + m_Event->SetMMStripXEStripNbr(g+1); + m_Event->SetMMStripXEEnergy(energyX); + + // Time + m_Event->SetMMStripXTDetectorNbr(detectorNbr); + m_Event->SetMMStripXTStripNbr(g+1); + m_Event->SetMMStripXTTime(RandGauss::shoot(time, ResoTimeMust)); + } } - Pos_X_itr++; } - - // Pos Y - Pos_Y_itr = PosYHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < PosYHitMap->entries() ; h++) { - G4int PosYTrackID = Pos_Y_itr->first - N ; - G4double PosY = *(Pos_Y_itr->second) ; - if (PosYTrackID == NTrackID) { - ms_InterCoord->SetDetectedPositionY(PosY) ; - } - Pos_Y_itr++; + } + // Y + if(energyY>0.1*keV){ // above threshold + // Pixel value at interaction point + unsigned int a,r,g,b; + // pixel + SiScorer->GetARGBBack(indexes[i],a,r,g,b); + if(r==0){ + trig.insert(detectorNbr); + // Energy + m_Event->SetMMStripYEDetectorNbr(detectorNbr); + m_Event->SetMMStripYEStripNbr(b+1); + m_Event->SetMMStripYEEnergy(energyY); + + // Time + m_Event->SetMMStripYTDetectorNbr(detectorNbr); + m_Event->SetMMStripYTStripNbr(b+1); + m_Event->SetMMStripYTTime(RandGauss::shoot(time, ResoTimeMust)); } - - // Pos Z - Pos_Z_itr = PosZHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < PosZHitMap->entries() ; h++) { - G4int PosZTrackID = Pos_Z_itr->first - N ; - G4double PosZ = *(Pos_Z_itr->second) ; - if (PosZTrackID == NTrackID) { - ms_InterCoord->SetDetectedPositionZ(PosZ) ; + else{ // Interstrip Y, keep both strip with shared energy + double rand = G4UniformRand(); + double energyY1 = rand*energyY; + if(energyY1>0.1*keV){ + trig.insert(detectorNbr); + + // Energy + m_Event->SetMMStripYEDetectorNbr(detectorNbr); + m_Event->SetMMStripYEStripNbr(b+1); + m_Event->SetMMStripYEEnergy(energyY1); + + // Time + m_Event->SetMMStripYTDetectorNbr(detectorNbr); + m_Event->SetMMStripYTStripNbr(b+1); + m_Event->SetMMStripYTTime(RandGauss::shoot(time, ResoTimeMust)); } - Pos_Z_itr++; - } - // Angle Theta - Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < AngThetaHitMap->entries(); h++) { - G4int AngThetaTrackID = Ang_Theta_itr->first - N ; - G4double AngTheta = *(Ang_Theta_itr->second) ; - if (AngThetaTrackID == NTrackID) { - ms_InterCoord->SetDetectedAngleTheta(AngTheta) ; + if(energyY1>0.1*keV){ + trig.insert(detectorNbr); + double energyY2 = (1-rand)*energyY; + // Energy + m_Event->SetMMStripYEDetectorNbr(detectorNbr); + m_Event->SetMMStripYEStripNbr(g+1); + m_Event->SetMMStripYEEnergy(energyY2); + + // Time + m_Event->SetMMStripYTDetectorNbr(detectorNbr); + m_Event->SetMMStripYTStripNbr(g+1); + m_Event->SetMMStripYTTime(RandGauss::shoot(time, ResoTimeMust)); } - Ang_Theta_itr++; } + } + } - // Angle Phi - Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < AngPhiHitMap->entries() ; h++) { - G4int AngPhiTrackID = Ang_Phi_itr->first - N ; - G4double AngPhi = *(Ang_Phi_itr->second) ; - if (AngPhiTrackID == NTrackID) { - ms_InterCoord->SetDetectedAnglePhi(AngPhi) ; - } - Ang_Phi_itr++; - } - - // Si(Li) - SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ; - SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin() ; - - for (G4int h = 0 ; h < SiLiEnergyHitMap->entries() ; h++) { - G4int SiLiEnergyTrackID = SiLiEnergy_itr->first -N ; - G4double SiLiEnergy = *(SiLiEnergy_itr->second) ; - - G4int PadNbr = *(SiLiPadNbr_itr->second) ; + // Si(Li) + NPS::HitsMap<G4double*>* SiLiHitMap; + std::map<G4int, G4double**>::iterator SiLi_itr; - if (SiLiEnergyTrackID == NTrackID) { - m_Event->SetMMSiLiEEnergy(RandGauss::shoot(SiLiEnergy, ResoSiLi)) ; - m_Event->SetMMSiLiEPadNbr(PadNbr); - m_Event->SetMMSiLiTPadNbr(PadNbr); - m_Event->SetMMSiLiTTime(1); - m_Event->SetMMSiLiTDetectorNbr(N); - m_Event->SetMMSiLiEDetectorNbr(N); - } + G4int SiLiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiScorer"); + SiLiHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(SiLiCollectionID)); - SiLiEnergy_itr++ ; - SiLiPadNbr_itr++ ; + // CsI + NPS::HitsMap<G4double*>* CsIHitMap; + std::map<G4int, G4double**>::iterator CsI_itr; + + G4int CsICollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIScorer"); + CsIHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CsICollectionID)); + + + // Look for SiLi data in Trigg Telescope + std::set<int>::iterator itr; + for(itr=trig.begin();itr!=trig.end();itr++){ + for(SiLi_itr = SiLiHitMap->GetMap()->begin(); SiLi_itr!=SiLiHitMap->GetMap()->end() ; SiLi_itr++){ + G4double* Info = *(SiLi_itr->second); + if(Info[7]==*itr){//matching telescope number + m_Event->SetMMSiLiEDetectorNbr(Info[7]); + m_Event->SetMMSiLiEEnergy(RandGauss::shoot(Info[0],ResoSiLi)); + m_Event->SetMMSiLiEPadNbr(Info[8]); + m_Event->SetMMSiLiTDetectorNbr(Info[7]); + m_Event->SetMMSiLiTTime(RandGauss::shoot(Info[1],ResoTimeMust)); + m_Event->SetMMSiLiTPadNbr(Info[8]); } + } + } - // CsI - CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ; - CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin() ; - - for (G4int h = 0 ; h < CsIEnergyHitMap->entries() ; h++) { - G4int CsIEnergyTrackID = CsIEnergy_itr->first-N ; - G4double CsIEnergy = *(CsIEnergy_itr->second) ; - - G4int CristalNumber = *(CsICristalNbr_itr->second) ; - if (CsIEnergyTrackID == NTrackID) { - m_Event->SetMMCsIEEnergy(RandGauss::shoot(CsIEnergy, ResoCsI*sqrt(CsIEnergy))); - m_Event->SetMMCsIECristalNbr(CristalNumber); - m_Event->SetMMCsITCristalNbr(CristalNumber); - m_Event->SetMMCsITTime(1); - m_Event->SetMMCsITDetectorNbr(N); - m_Event->SetMMCsIEDetectorNbr(N); - } - - CsIEnergy_itr++ ; - CsICristalNbr_itr++ ; + // Look for CsI data in Trigg Telescope + for(itr=trig.begin();itr!=trig.end();itr++){ + for(CsI_itr = CsIHitMap->GetMap()->begin(); CsI_itr!=CsIHitMap->GetMap()->end() ; CsI_itr++){ + G4double* Info = *(CsI_itr->second); + + if(Info[7]==*itr){//matching telescope number + m_Event->SetMMCsIEDetectorNbr(Info[7]); + m_Event->SetMMCsIEEnergy(RandGauss::shoot(Info[0],ResoCsI)); + m_Event->SetMMCsIECristalNbr(Info[8]); + m_Event->SetMMCsITDetectorNbr(Info[7]); + m_Event->SetMMCsITTime(RandGauss::shoot(Info[1],ResoTimeMust)); + m_Event->SetMMCsITCristalNbr(Info[8]); } - } - - DetectorNumber_itr++; } + + SiScorer->clear(); // clear map for next event - DetectorNumberHitMap ->clear(); - EnergyHitMap ->clear() ; - TimeHitMap ->clear() ; - XHitMap ->clear() ; - YHitMap ->clear() ; - SiLiEnergyHitMap ->clear() ; - SiLiPadNbrHitMap ->clear() ; - CsIEnergyHitMap ->clear() ; - CsICristalNbrHitMap ->clear() ; - PosXHitMap ->clear() ; - PosYHitMap ->clear() ; - PosZHitMap ->clear() ; - AngThetaHitMap ->clear() ; - AngPhiHitMap ->clear() ; + SiLiHitMap->clear() ; + CsIHitMap->clear() ; } @@ -1037,49 +911,26 @@ void MUST2Array::InitializeScorers() { // if the scorer were created previously nothing else need to be made if(already_exist) return; - G4VPrimitiveScorer* DetNbr = new OBSOLETEGENERALSCORERS::PSDetectorNumber("DetectorNumber","MUST2Telescope", 0); - G4VPrimitiveScorer* Energy = new OBSOLETEGENERALSCORERS::PSEnergy("StripEnergy","MUST2Telescope", 0); - G4VPrimitiveScorer* TOF = new OBSOLETEGENERALSCORERS::PSTOF("StripTime","MUST2Telescope", 0); - - G4VPrimitiveScorer* StripPositionX = new PSStripNumberX("StripNumberX", 0, SiliconFace, 128); - G4VPrimitiveScorer* StripPositionY = new PSStripNumberY("StripNumberY", 0, SiliconFace, 128); - - G4VPrimitiveScorer* InteractionCoordinatesX = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","MUST2Telescope", 0); - G4VPrimitiveScorer* InteractionCoordinatesY = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","MUST2Telescope", 0); - G4VPrimitiveScorer* InteractionCoordinatesZ = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","MUST2Telescope", 0); - - G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","MUST2Telescope", 0); - G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","MUST2Telescope", 0) ; - + G4VPrimitiveScorer* SiScorer = new SILICONSCORERS::PS_Silicon_Images("SiScorer","/scratch/nptool/NPLib/Detectors/MUST2/maskFront.png","/scratch/nptool/NPLib/Detectors/MUST2/maskBack.png",0.01,0.01,0,0,0xffff0000,0); //and register it to the multifunctionnal detector - m_StripScorer->RegisterPrimitive(DetNbr); - m_StripScorer->RegisterPrimitive(Energy); - m_StripScorer->RegisterPrimitive(TOF); - m_StripScorer->RegisterPrimitive(StripPositionX); - m_StripScorer->RegisterPrimitive(StripPositionY); - m_StripScorer->RegisterPrimitive(InteractionCoordinatesX); - m_StripScorer->RegisterPrimitive(InteractionCoordinatesY); - m_StripScorer->RegisterPrimitive(InteractionCoordinatesZ); - m_StripScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta); - m_StripScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi); - + m_StripScorer->RegisterPrimitive(SiScorer); // SiLi Associate Scorer - G4VPrimitiveScorer* SiLiEnergy = new OBSOLETEGENERALSCORERS::PSEnergy("SiLiEnergy","MUST2Telescope", 0) ; - G4VPrimitiveScorer* SiLiPadNbr = new PSPadOrCristalNumber("SiLiPadNbr",0) ; - m_SiLiScorer->RegisterPrimitive(SiLiEnergy); - m_SiLiScorer->RegisterPrimitive(SiLiPadNbr); + vector<int> SiLi_nesting={3,0}; + G4VPrimitiveScorer* SiLiScorer= + new CALORIMETERSCORERS::PS_CalorimeterWithInteraction("SiLiScorer",SiLi_nesting) ; + m_SiLiScorer->RegisterPrimitive(SiLiScorer); // CsI Associate Scorer - G4VPrimitiveScorer* CsIEnergy = new OBSOLETEGENERALSCORERS::PSEnergy("CsIEnergy","MUST2Telescope", 0) ; - G4VPrimitiveScorer* CsICristalNbr = new PSPadOrCristalNumber("CsICristalNbr",0) ; - m_CsIScorer->RegisterPrimitive(CsIEnergy) ; - m_CsIScorer->RegisterPrimitive(CsICristalNbr) ; + vector<int> CsI_nesting = {2,0}; + G4VPrimitiveScorer* CsIScorer= + new CALORIMETERSCORERS::PS_CalorimeterWithInteraction("CsIScorer",CsI_nesting, 0) ; + m_CsIScorer->RegisterPrimitive(CsIScorer) ; // Add All Scorer to the Global Scorer Manager - G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer) ; - G4SDManager::GetSDMpointer()->AddNewDetector(m_SiLiScorer) ; - G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ; + G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer); + G4SDManager::GetSDMpointer()->AddNewDetector(m_SiLiScorer); + G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.hh b/NPSimulation/Detectors/MUST2/MUST2Array.hh index ff06606aead1d8659c207fdad7dc337441b541c7..68914a94e16a12eab14c52e0242c9a41fc3c3f7b 100644 --- a/NPSimulation/Detectors/MUST2/MUST2Array.hh +++ b/NPSimulation/Detectors/MUST2/MUST2Array.hh @@ -47,7 +47,7 @@ namespace MUST2 const G4double AluStripThickness = 0.4*micrometer ; const G4double SiliconThickness = 300*micrometer ; - const G4double SiliconFace = 98*mm ; + const G4double SiliconFace = 100.42*mm ; const G4double VacBoxThickness = 3*cm ; const G4double SiLiThickness = 5.1*mm; // Must be checked diff --git a/NPSimulation/Scorers/CalorimeterScorers.cc b/NPSimulation/Scorers/CalorimeterScorers.cc index c892f401e67b94e6c11a0dfb8c3fe70744ccf510..ef102b9ce931c7df2c5f245340c84b739e1ef157 100644 --- a/NPSimulation/Scorers/CalorimeterScorers.cc +++ b/NPSimulation/Scorers/CalorimeterScorers.cc @@ -126,7 +126,6 @@ G4bool PS_CalorimeterWithInteraction::ProcessHits(G4Step* aStep, G4TouchableHist Infos[5] = m_Position.theta(); Infos[6] = m_Position.phi(); - for(unsigned int i = 0 ; i < mysize ; i++){ Infos[i+7] = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_NestingLevel[i]); } diff --git a/NPSimulation/Scorers/SiliconScorers.cc b/NPSimulation/Scorers/SiliconScorers.cc index acc039beb89504717eb1c3d2ad7a04c2bd1419ea..f3e176108635ceecfdc1125ba233a170a33f964a 100644 --- a/NPSimulation/Scorers/SiliconScorers.cc +++ b/NPSimulation/Scorers/SiliconScorers.cc @@ -22,7 +22,159 @@ #include "SiliconScorers.hh" #include "G4UnitsTable.hh" using namespace SILICONSCORERS ; +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +PS_Silicon_Images::PS_Silicon_Images(G4String name, string imageFront,string imageBack,double scalingFront, double scalingBack, double centerOffsetX,double centerOffsetY,unsigned int ignoreValue, G4int depth) :G4VPrimitiveScorer(name, depth),HCID(-1){ + m_ImageFront = new NPL::Image(imageFront,scalingFront,scalingFront); + m_ImageBack = new NPL::Image(imageBack,scalingBack,scalingBack); + m_ScalingFront = scalingFront; + m_ScalingBack = scalingBack; + m_CenterOffsetX = centerOffsetX; + m_CenterOffsetY = centerOffsetY; + m_IgnoreValue = ignoreValue; + m_Level = depth; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +G4bool PS_Silicon_Images::ProcessHits(G4Step* aStep, G4TouchableHistory*){ + // contain Energy Time, DetNbr, PixelFront and PixelBack + t_Energy = aStep->GetTotalEnergyDeposit(); + t_Time = aStep->GetPreStepPoint()->GetGlobalTime(); + + t_DetectorNbr = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); + t_Position = aStep->GetPreStepPoint()->GetPosition(); + + // Interaction coordinates (used to fill the InteractionCoordinates branch) + t_X = t_Position.x(); + t_Y = t_Position.y(); + t_Z = t_Position.z(); + t_Theta = t_Position.theta(); + t_Phi = t_Position.phi(); + + t_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(t_Position); + + t_PixelFront = m_ImageFront->GetPixelAtCoordinate(t_Position.x(),t_Position.y()); + t_PixelBack = m_ImageBack->GetPixelAtCoordinate(t_Position.x(),t_Position.y()); + + unsigned int front = (m_ImageFront->GetPixelAtCoordinate(t_Position.x(),t_Position.y()))&0xff; + unsigned int back = (m_ImageBack->GetPixelAtCoordinate(t_Position.x(),t_Position.y()))&0xff; + + + // If front and back are in inactive part of the wafer, + // nothing is added to the map + if(t_PixelFront == m_IgnoreValue && t_PixelBack == m_IgnoreValue) + return FALSE; + + t_Index = aStep->GetTrack()->GetTrackID() + t_DetectorNbr * 1e3 ; + // Check if the particle has interact before, if yes, add up the energies. + map<unsigned int, double>::iterator it; + it= m_Energy.find(t_Index); + if(it!=m_Energy.end()){ + double dummy = it->second; + t_Energy+=dummy; + } + + m_Energy[t_Index] = t_Energy; + m_Time[t_Index] = t_Time; + m_DetectorNbr[t_Index] = t_DetectorNbr; + m_PixelFront[t_Index] = t_PixelFront; + m_PixelBack[t_Index] = t_PixelBack; + m_X[t_Index] = t_X; + m_Y[t_Index] = t_Y; + m_Z[t_Index] = t_Z; + m_Theta[t_Index] = t_Theta; + m_Phi[t_Index] = t_Phi; + + + return TRUE; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_Silicon_Images::Initialize(G4HCofThisEvent* HCE){ +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_Silicon_Images::EndOfEvent(G4HCofThisEvent*){ +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_Silicon_Images::clear(){ + m_Energy.clear(); + m_Time.clear(); + m_DetectorNbr.clear(); + m_PixelFront.clear(); + m_PixelBack.clear(); + m_X.clear(); + m_Y.clear(); + m_Z.clear(); + m_Theta.clear(); + m_Phi.clear(); +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +vector<unsigned int> PS_Silicon_Images::GetIndexes(){ + vector<unsigned int > indexes; + map<unsigned int , double>::iterator it; + for(it=m_Energy.begin(); it!=m_Energy.end();it++) + indexes.push_back(it->first); + + return indexes; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +unsigned int PS_Silicon_Images::GetPixelFront(unsigned int index){ + return m_PixelFront[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +unsigned int PS_Silicon_Images::GetPixelBack(unsigned int index){ + return m_PixelBack[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_Silicon_Images::GetARGBFront(unsigned int index,unsigned int& a,unsigned int& r,unsigned int& g,unsigned int& b){ + unsigned int Info =m_PixelFront[index]; + a = (Info>>24)&0xff; + r = (Info>>16)&0xff; + g = (Info>>8)&0xff; + b = (Info>>0)&0xff; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void PS_Silicon_Images::GetARGBBack(unsigned int index,unsigned int& a,unsigned int& r,unsigned int& g,unsigned int& b){ + unsigned int Info =m_PixelBack[index]; + + a = (Info>>24)&0xff; + r = (Info>>16)&0xff; + g = (Info>>8)&0xff; + b = (Info>>0)&0xff; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetEnergy(unsigned int index){ + return m_Energy[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetTime(unsigned int index){ + return m_Time[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +unsigned int PS_Silicon_Images::GetDetectorNbr(unsigned int index){ + return m_DetectorNbr[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetX(unsigned int index){ + return m_X[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetY(unsigned int index){ + return m_Y[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetZ(unsigned int index){ + return m_Z[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetTheta(unsigned int index){ + return m_Theta[index]; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +double PS_Silicon_Images::GetPhi(unsigned int index){ + return m_Phi[index]; +} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... PS_Silicon_Rectangle::PS_Silicon_Rectangle(G4String name,G4int Level, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4int depth,G4String axis) :G4VPrimitiveScorer(name, depth),HCID(-1){ @@ -192,7 +344,7 @@ G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){ phi+=2*M_PI; m_StripSectorNumber = (int) ((phi - m_PhiStart) / m_StripPitchSector ) + 1 ; m_StripQuadrantNumber = (int) ((phi - m_PhiStart) /m_StripPitchQuadrant) + 1 ; - + //Rare case where particle is close to edge of silicon plan if (m_StripRingNumber > m_NumberOfStripRing) m_StripRingNumber = m_NumberOfStripRing; if (m_StripSectorNumber > m_NumberOfStripSector) m_StripSectorNumber = m_NumberOfStripSector; diff --git a/NPSimulation/Scorers/SiliconScorers.hh b/NPSimulation/Scorers/SiliconScorers.hh index 01b4fbffca1f668610bb9369e89c675e0ed44f39..ccbc133f0ce6ec799af195853df46c76b02ade17 100644 --- a/NPSimulation/Scorers/SiliconScorers.hh +++ b/NPSimulation/Scorers/SiliconScorers.hh @@ -26,13 +26,87 @@ *****************************************************************************/ #include "G4VPrimitiveScorer.hh" #include "NPSHitsMap.hh" - +#include "NPImage.h" #include <map> using namespace std; using namespace CLHEP; namespace SILICONSCORERS { -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + class PS_Silicon_Images : public G4VPrimitiveScorer{ + + public: // with description + PS_Silicon_Images(G4String name, string imageFront,string imageBack,double scalingFront, double scalingBack, double centerOffsetX,double centerOffsetY,unsigned int ignoreValue, G4int depth=0); + ~PS_Silicon_Images(){}; + + protected: // with description + G4bool ProcessHits(G4Step*, G4TouchableHistory*); + + public: + void Initialize(G4HCofThisEvent*); + void EndOfEvent(G4HCofThisEvent*); + void clear(); + void DrawAll(){}; + void PrintAll(){}; + + private: // Geometry of the detector + NPL::Image* m_ImageFront; + NPL::Image* m_ImageBack; + double m_ScalingFront; + double m_ScalingBack; + double m_CenterOffsetX; + double m_CenterOffsetY; + unsigned int m_IgnoreValue; + + // Level at which to find the copy number linked to the detector number + G4int m_Level; + + private: // inherited from G4VPrimitiveScorer + G4int HCID; + map<unsigned int,double> m_Energy; + map<unsigned int,double> m_Time; + map<unsigned int,unsigned int> m_DetectorNbr; + map<unsigned int,unsigned int> m_PixelFront; + map<unsigned int,unsigned int> m_PixelBack; + map<unsigned int,double> m_X; + map<unsigned int,double> m_Y; + map<unsigned int,double> m_Z; + map<unsigned int,double> m_Theta; + map<unsigned int,double> m_Phi; + + private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) + G4ThreeVector t_Position; + unsigned int t_Index; + double t_Energy; + double t_Time; + unsigned int t_DetectorNbr; + unsigned int t_PixelFront; + unsigned int t_PixelBack; + double t_X; + double t_Y; + double t_Z; + double t_Theta; + double t_Phi; + + public: // information accessor + vector<unsigned int> GetIndexes(); + unsigned int GetPixelFront(unsigned int index); + unsigned int GetPixelBack(unsigned int index); + void GetARGBFront(unsigned int index,unsigned int& a,unsigned int& r,unsigned int& g,unsigned int& b); + void GetARGBBack(unsigned int index,unsigned int& a,unsigned int& r,unsigned int& g,unsigned int& b); + + double GetEnergy(unsigned int index); + double GetTime(unsigned int index); + unsigned int GetDetectorNbr(unsigned int index); + double GetX(unsigned int index); + double GetY(unsigned int index); + double GetZ(unsigned int index); + double GetTheta(unsigned int index); + double GetPhi(unsigned int index); + + }; + + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... struct RectangularOutput { G4double totalEnergy; @@ -48,46 +122,46 @@ namespace SILICONSCORERS { }; class PS_Silicon_Rectangle : public G4VPrimitiveScorer{ - - public: // with description - PS_Silicon_Rectangle(G4String name, G4int Level, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4int depth=0,G4String axis="xy"); - ~PS_Silicon_Rectangle(); - - protected: // with description - G4bool ProcessHits(G4Step*, G4TouchableHistory*); - - public: - void Initialize(G4HCofThisEvent*); - void EndOfEvent(G4HCofThisEvent*); - void clear(); - void DrawAll(); - void PrintAll(); - - private: // Geometry of the detector - G4double m_StripPlaneLength; - G4double m_StripPlaneWidth; - G4int m_NumberOfStripLength; - G4int m_NumberOfStripWidth; - G4double m_StripPitchLength; - G4double m_StripPitchWidth; - G4String m_Axis; - // Level at which to find the copy number linked to the detector number - G4int m_Level; - - private: // inherited from G4VPrimitiveScorer - G4int HCID; - NPS::HitsMap<G4double*>* EvtMap; - - private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) - G4ThreeVector m_Position ; - G4int m_DetectorNumber ; - G4int m_StripLengthNumber ; - G4int m_StripWidthNumber ; - G4int m_Index ; - + + public: // with description + PS_Silicon_Rectangle(G4String name, G4int Level, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripLength,G4int NumberOfStripWidth,G4int depth=0,G4String axis="xy"); + ~PS_Silicon_Rectangle(); + + protected: // with description + G4bool ProcessHits(G4Step*, G4TouchableHistory*); + + public: + void Initialize(G4HCofThisEvent*); + void EndOfEvent(G4HCofThisEvent*); + void clear(); + void DrawAll(); + void PrintAll(); + + private: // Geometry of the detector + G4double m_StripPlaneLength; + G4double m_StripPlaneWidth; + G4int m_NumberOfStripLength; + G4int m_NumberOfStripWidth; + G4double m_StripPitchLength; + G4double m_StripPitchWidth; + G4String m_Axis; + // Level at which to find the copy number linked to the detector number + G4int m_Level; + + private: // inherited from G4VPrimitiveScorer + G4int HCID; + NPS::HitsMap<G4double*>* EvtMap; + + private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) + G4ThreeVector m_Position ; + G4int m_DetectorNumber ; + G4int m_StripLengthNumber ; + G4int m_StripWidthNumber ; + G4int m_Index ; + }; - -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... struct AnnularOutput { G4double totalEnergy; @@ -104,51 +178,51 @@ namespace SILICONSCORERS { }; class PS_Silicon_Annular : public G4VPrimitiveScorer{ - - public: // with description - PS_Silicon_Annular(G4String name,G4int Level, G4double StripPlaneInnerRadius, G4double StripPlaneOuterRadius, G4double PhiStart,G4double PhiStop, G4int NumberOfStripRing,G4int NumberOfStripSector=1, G4int NumberOfQuadrant=1,G4int depth=0); - ~PS_Silicon_Annular(); - - protected: // with description - G4bool ProcessHits(G4Step*, G4TouchableHistory*); - - public: - void Initialize(G4HCofThisEvent*); - void EndOfEvent(G4HCofThisEvent*); - void clear(); - void DrawAll(); - void PrintAll(); - - private: // Geometry of the detector - G4double m_StripPlaneInnerRadius; - G4double m_StripPlaneOuterRadius; - G4double m_PhiStart; - G4double m_PhiStop; - G4int m_NumberOfStripRing; - G4int m_NumberOfStripSector; - G4int m_NumberOfStripQuadrant; - G4double m_StripPitchRing; - G4double m_StripPitchSector; - G4double m_StripPitchQuadrant; - G4String m_Axis; - // Level at which to find the copy number linked to the detector number - G4int m_Level; - - private: // inherited from G4VPrimitiveScorer - G4int HCID; - NPS::HitsMap<G4double*>* EvtMap; - - private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) - G4ThreeVector m_Position ; - G4ThreeVector m_uz ; - G4int m_DetectorNumber ; - G4int m_StripRingNumber ; - G4int m_StripSectorNumber ; - G4int m_StripQuadrantNumber ; - G4int m_Index ; - + + public: // with description + PS_Silicon_Annular(G4String name,G4int Level, G4double StripPlaneInnerRadius, G4double StripPlaneOuterRadius, G4double PhiStart,G4double PhiStop, G4int NumberOfStripRing,G4int NumberOfStripSector=1, G4int NumberOfQuadrant=1,G4int depth=0); + ~PS_Silicon_Annular(); + + protected: // with description + G4bool ProcessHits(G4Step*, G4TouchableHistory*); + + public: + void Initialize(G4HCofThisEvent*); + void EndOfEvent(G4HCofThisEvent*); + void clear(); + void DrawAll(); + void PrintAll(); + + private: // Geometry of the detector + G4double m_StripPlaneInnerRadius; + G4double m_StripPlaneOuterRadius; + G4double m_PhiStart; + G4double m_PhiStop; + G4int m_NumberOfStripRing; + G4int m_NumberOfStripSector; + G4int m_NumberOfStripQuadrant; + G4double m_StripPitchRing; + G4double m_StripPitchSector; + G4double m_StripPitchQuadrant; + G4String m_Axis; + // Level at which to find the copy number linked to the detector number + G4int m_Level; + + private: // inherited from G4VPrimitiveScorer + G4int HCID; + NPS::HitsMap<G4double*>* EvtMap; + + private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) + G4ThreeVector m_Position ; + G4ThreeVector m_uz ; + G4int m_DetectorNumber ; + G4int m_StripRingNumber ; + G4int m_StripSectorNumber ; + G4int m_StripQuadrantNumber ; + G4int m_Index ; + }; -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... struct ResistiveOutput { G4double upstreamEnergy; G4double downstreamEnergy; @@ -163,42 +237,42 @@ namespace SILICONSCORERS { }; class PS_Silicon_Resistive : public G4VPrimitiveScorer{ - - public: // with description - PS_Silicon_Resistive(G4String name, G4int Level, - G4double StripPlaneLength, G4double StripPlaneWidth, - G4int NumberOfStripWidth,G4int depth=0); - - ~PS_Silicon_Resistive(); - - protected: // with description - G4bool ProcessHits(G4Step*, G4TouchableHistory*); - - public: - void Initialize(G4HCofThisEvent*); - void EndOfEvent(G4HCofThisEvent*); - void clear(); - void DrawAll(); - void PrintAll(); - - private: // Geometry of the detector - G4double m_StripPlaneLength; - G4double m_StripPlaneWidth; - G4int m_NumberOfStripWidth; - G4double m_StripPitchWidth; - // Level at which to find the copy number linked to the detector number - G4int m_Level; - - private: // inherited from G4VPrimitiveScorer - G4int HCID; - NPS::HitsMap<G4double*>* EvtMap; - - private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) - G4ThreeVector m_Position ; - G4int m_DetectorNumber ; - G4int m_StripWidthNumber ; - G4int m_Index ; - + + public: // with description + PS_Silicon_Resistive(G4String name, G4int Level, + G4double StripPlaneLength, G4double StripPlaneWidth, + G4int NumberOfStripWidth,G4int depth=0); + + ~PS_Silicon_Resistive(); + + protected: // with description + G4bool ProcessHits(G4Step*, G4TouchableHistory*); + + public: + void Initialize(G4HCofThisEvent*); + void EndOfEvent(G4HCofThisEvent*); + void clear(); + void DrawAll(); + void PrintAll(); + + private: // Geometry of the detector + G4double m_StripPlaneLength; + G4double m_StripPlaneWidth; + G4int m_NumberOfStripWidth; + G4double m_StripPitchWidth; + // Level at which to find the copy number linked to the detector number + G4int m_Level; + + private: // inherited from G4VPrimitiveScorer + G4int HCID; + NPS::HitsMap<G4double*>* EvtMap; + + private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) + G4ThreeVector m_Position ; + G4int m_DetectorNumber ; + G4int m_StripWidthNumber ; + G4int m_Index ; + }; }