diff --git a/source/digits_hits/src/GateSpatialResolution.cc b/source/digits_hits/src/GateSpatialResolution.cc index 6c0b3573e..e1c1c762b 100644 --- a/source/digits_hits/src/GateSpatialResolution.cc +++ b/source/digits_hits/src/GateSpatialResolution.cc @@ -15,7 +15,7 @@ Includes functionalities from: GateSpblurring GateCC3DlocalSpblurring - GateDoIModels (TODO) + GateDoIModels Previous authors: Steven.Staelens@rug.ac.be(?), AE @@ -65,12 +65,12 @@ GateSpatialResolution::GateSpatialResolution(GateSinglesDigitizer *digitizer, G4 m_fwhmXDistrib(0), m_fwhmYDistrib(0), m_fwhmZDistrib(0), - m_nameAxis("YZ"), + m_nameAxis(""), m_fwhmXDistrib2D(0), m_fwhmYDistrib2D(0), m_fwhmZDistrib2D(0), m_IsConfined(true), - m_UseTruncatedGaussian(true), + m_UseTruncatedGaussian(false), m_Navigator(0), m_Touchable(0), m_systemDepth(-1), @@ -129,12 +129,30 @@ void GateSpatialResolution::SetSpatialResolutionParameters() { abort(); } + + if (m_nameAxis == "YZ" && (m_fwhmYDistrib2D ==0 || m_fwhmZDistrib2D ==0)) + { + G4cout << "***ERROR*** Spatial Resolution 2D distribution: you choose YZ axis but one or both parameters /fwhmYDistrib2D or /m_fwhmZDistrib2D are not defined" << G4endl; + abort(); + } + if (m_nameAxis == "XZ" && (m_fwhmXDistrib2D ==0 || m_fwhmZDistrib2D ==0)) + { + G4cout << "***ERROR*** Spatial Resolution 2D distribution: you choose XZ axis but one or both parameters /fwhmXDistrib2D or /m_fwhmZDistrib2D are not defined" << G4endl; + abort(); + } + + + + + + } void GateSpatialResolution::Digitize(){ + GateVSystem* m_system = ((GateSinglesDigitizer*)this->GetDigitizer())->GetSystem(); if (m_IsFirstEntrance) { @@ -190,6 +208,17 @@ void GateSpatialResolution::Digitize(){ GateDigi* inputDigi; + if(!m_IsConfined && !m_Navigator) + { + //Getting world Volume + // Do not use from TransportationManager as it is not recommended + G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4VPhysicalVolume *WorldVolume = navigator->GetWorldVolume(); + m_Navigator = new G4Navigator(); + m_Navigator->SetWorldVolume(WorldVolume); + } + + if (IDC) { G4int n_digi = IDC->entries(); @@ -208,39 +237,41 @@ void GateSpatialResolution::Digitize(){ G4double Pz = P.z(); G4double stddevX = 0., stddevY = 0., stddevZ = 0.; - // Use the configured axis pair (m_nameAxis) to evaluate Value2D. // Allowed pairs for PET: "XZ" or "YZ" (default "YZ"). - if (m_fwhmXDistrib2D || m_fwhmYDistrib2D || m_fwhmZDistrib2D) { - if (m_nameAxis == "XZ") { - if (m_fwhmXDistrib2D) stddevX = m_fwhmXDistrib2D->Value2D(P.x() * mm, P.z() * mm); - if (m_fwhmYDistrib2D) stddevY = m_fwhmYDistrib2D->Value2D(P.x() * mm, P.z() * mm); - if (m_fwhmZDistrib2D) stddevZ = m_fwhmZDistrib2D->Value2D(P.x() * mm, P.z() * mm); - } else { // YZ - if (m_fwhmXDistrib2D) stddevX = m_fwhmXDistrib2D->Value2D(P.y() * mm, P.z() * mm); - if (m_fwhmYDistrib2D) stddevY = m_fwhmYDistrib2D->Value2D(P.y() * mm, P.z() * mm); - if (m_fwhmZDistrib2D) stddevZ = m_fwhmZDistrib2D->Value2D(P.y() * mm, P.z() * mm); - } + if (m_nameAxis == "XZ") { //if 2D distributions are defined + stddevX = m_fwhmXDistrib2D->Value2D(P.x() * mm, P.z() * mm); + stddevZ = m_fwhmZDistrib2D->Value2D(P.x() * mm, P.z() * mm); } - else { - - if (m_fwhmXDistrib) stddevX = m_fwhmXDistrib->Value(P.x() * mm); - else if (fwhmX) stddevX = fwhmX / GateConstants::fwhm_to_sigma; - - if (m_fwhmYDistrib) stddevY = m_fwhmYDistrib->Value(P.y() * mm); - else if (fwhmY) stddevY = fwhmY / GateConstants::fwhm_to_sigma; - - if (m_fwhmZDistrib) stddevZ = m_fwhmZDistrib->Value(P.z() * mm); - else if (fwhmZ) stddevZ = fwhmZ / GateConstants::fwhm_to_sigma; - + if(m_nameAxis == "YZ") { // YZ + stddevY = m_fwhmYDistrib2D->Value2D(P.y() * mm, P.z() * mm); + stddevZ = m_fwhmZDistrib2D->Value2D(P.y() * mm, P.z() * mm); } - - - + + + if (m_fwhmXDistrib) stddevX = m_fwhmXDistrib->Value(P.x() * mm); + else if (fwhmX) + stddevX = fwhmX / GateConstants::fwhm_to_sigma; + + + if (m_fwhmYDistrib) stddevY = m_fwhmYDistrib->Value(P.y() * mm); + else if (fwhmY) stddevY = fwhmY / GateConstants::fwhm_to_sigma; + + if (m_fwhmZDistrib) stddevZ = m_fwhmZDistrib->Value(P.z() * mm); + else if (fwhmZ) stddevZ = fwhmZ / GateConstants::fwhm_to_sigma; + + // } + + + if(!stddevX || !stddevY || !stddevZ) + { + G4cout << " *** ERROR*** GateSpatialResolution::Digitize. Resolution for one of the axis is not set. If you really don't want to apply spatial resolution on one of the axis, please, contact kochebina [a] gmail (dot) com " << G4endl; + abort(); + } G4double PxNew ; G4double PyNew ; G4double PzNew ; -// store the computed stddevs into the digi for later ROOT output + // store the computed stddevs into the digi for later ROOT output m_outputDigi->SetSpatialRes2DStdDevX(stddevX); m_outputDigi->SetSpatialRes2DStdDevY(stddevY); m_outputDigi->SetSpatialRes2DStdDevZ(stddevZ); @@ -314,15 +345,11 @@ void GateSpatialResolution::Digitize(){ if(PxNew>Xmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; + m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC - m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC - //Getting world Volume - // Do not use from TransportationManager as it is not recommended - G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4VPhysicalVolume *WorldVolume = navigator->GetWorldVolume(); - m_Navigator = new G4Navigator(); - m_Navigator->SetWorldVolume(WorldVolume); + + m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC G4VPhysicalVolume* PV = m_Navigator->LocateGlobalPointAndSetup(m_outputDigi->GetGlobalPos()); m_Touchable = m_Navigator->CreateTouchableHistoryHandle(); @@ -332,6 +359,7 @@ void GateSpatialResolution::Digitize(){ { UpdateVolumeID(); } + m_OutputDigiCollection->insert(m_outputDigi); diff --git a/source/digits_hits/src/GateSpatialResolutionMessenger.cc b/source/digits_hits/src/GateSpatialResolutionMessenger.cc index 413b271e2..0197b5fab 100644 --- a/source/digits_hits/src/GateSpatialResolutionMessenger.cc +++ b/source/digits_hits/src/GateSpatialResolutionMessenger.cc @@ -57,6 +57,7 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol spresolutionZdistribCmd = new G4UIcmdWithAString(cmdName,this); spresolutionZdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmXDistrib2D"; spresolutionXDistrib2DCmd = new G4UIcmdWithAString(cmdName,this); spresolutionXDistrib2DCmd->SetGuidance("Set the 2D distribution for X axis (expects a 2D distribution object)");