From 41fc1a8d8e56f10a2e365497048c969358480d6f Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 11 Nov 2024 15:51:17 +0100 Subject: [PATCH 01/25] Update MCSource.cc --- src/cosima/src/MCSource.cc | 153 +++++++++++++++++++++++++++++++------ 1 file changed, 129 insertions(+), 24 deletions(-) diff --git a/src/cosima/src/MCSource.cc b/src/cosima/src/MCSource.cc index da549e8c..2f3495c1 100644 --- a/src/cosima/src/MCSource.cc +++ b/src/cosima/src/MCSource.cc @@ -14,8 +14,9 @@ * and all its terms. * */ - - + + + // Cosima: #include "MCCommon.hh" #include "MCSource.hh" @@ -116,7 +117,7 @@ const int MCSource::c_PolarizationAbsolute = 3; const int MCSource::c_PolarizationRelativeX = 4; const int MCSource::c_PolarizationRelativeY = 5; const int MCSource::c_PolarizationRelativeZ = 6; - +const int MCSource::c_PolarizationGalactic = 7; // Don't change this list because the ID is written to the Sim file const int MCSource::c_Gamma = 1; @@ -271,7 +272,7 @@ void MCSource::Initialize() m_PolarizationParam2 = c_Invalid; m_PolarizationParam3 = c_Invalid; m_PolarizationDegree = 0.0; - + m_UseFarFieldTransmissionProbability = false; m_UseEarthOccultation = false; m_ThetaMaxEarthOccultation = 0.0; @@ -430,7 +431,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) //Compute the theta max (angle between the particle direction and Earth-based zenith) //using the altitude from ori file - double R_Earth = 6378; //km + int R_Earth = 6378; //km ThetaMax = ( c_Pi - asin( R_Earth/(R_Earth+Alt/km) ) )/deg; } @@ -442,6 +443,16 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) // convert the Earth galactic l,b coordinate into cartesian representation G4ThreeVector EarthZenith(cos(Lat)*cos(Lon) , cos(Lat)*sin(Lon) , sin(Lat)); + + //compute zpointing of spacecraft + //double Xlat=0; + //double Xlon=0; + //double Zlat=0; + //double Zlon=0; + //Sky.GetOrientation(m_NextEmission,Xlat,Xlon,Zlat,Zlon); + //G4ThreeVector SpacecraftZaxis(cos(Zlat)*cos(Zlon) , cos(Zlat)*sin(Zlon) , sin(Zlat)); + //G4ThreeVector SpacecraftXaxis(cos(Xlat)*cos(Xlon) , cos(Xlat)*sin(Xlon) , sin(Xlat)); + //rotate the oriented coordinate system of the event into the local coordinate system //Sky.OrientDirectionInvers(m_NextEmission, Dir); Sky.OrientDirection(m_NextEmission, Dir); @@ -452,17 +463,20 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) //cout<< "particle Dir: "<GetCurrentSource()->GetEneDist()->SetMonoEnergy(0); } } else { - if (EarthZenith.angle(Dir)/deg >=ThetaMax ){ + if (abs(EarthZenith.angle(Dir)/deg) >=ThetaMax ){ ParticleGun->GetCurrentSource()->GetEneDist()->SetMonoEnergy(0); } } @@ -2256,6 +2270,7 @@ bool MCSource::SetPolarizationType(const int& PolarizationType) case c_PolarizationRelativeX: case c_PolarizationRelativeY: case c_PolarizationRelativeZ: + case c_PolarizationGalactic: m_PolarizationType = PolarizationType; return true; default: @@ -2292,6 +2307,9 @@ string MCSource::GetPolarizationTypeAsString() const case c_PolarizationRelativeZ: Name = "RelativeZ"; break; + case c_PolarizationGalactic: + Name = "Galactic"; + break; default: break; } @@ -2314,7 +2332,6 @@ bool MCSource::SetPolarization(double PolarizationParam1, return true; } - /****************************************************************************** * Set the degree of polarization 1.0 == 100% polarized */ @@ -3102,28 +3119,75 @@ bool MCSource::GeneratePosition(G4GeneralParticleSource* Gun) } else if (m_BeamType == c_FarFieldEarthOccultation) { // Determine a random position on the sphere between // theta min and theta max that are computed each time - // according to the pointing of the spacecraft in order - //to take into account the Earth Occultation + // according to the pointing of the spacecraft in order + //to take into account the Earth Occultation if (m_StartAreaType == c_StartAreaSphere) { - //determine the theta range + //determine the theta range - //get the Earth aspect information from the orientation file + //get the Earth aspect information from the orientation file const MCOrientation& Sky = MCRunManager::GetMCRunManager()->GetCurrentRun().GetSkyOrientationReference(); double Alt,Lat,Lon=0; - if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ - - + //double xLat,xLon,zLat,zLon=0; + if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ //&& Sky.GetOrientation(m_NextEmission,xLat,xLon,zLat,zLon)==true ){ + + double R_Earth = 6378; //km + + //compute theta max according to the current altitude + double ThetaMax = ( c_Pi - asin( R_Earth/(R_Earth+Alt/km) ) )/deg; + + // convert the Earth galactic l,b coordinate into cartesian representation + G4ThreeVector EarthZenith(cos(Lat)*cos(Lon) , cos(Lat)*sin(Lon) , sin(Lat)); + + // convert the spacecraft galacticz pointing l,b coordinate into cartesian representation + //G4ThreeVector SpacecraftZenith(cos(zLat)*cos(zLon) , cos(zLat)*sin(zLon) , sin(zLat)); + + //double ThetaRocking = EarthZenith.angle(SpacecraftZenith)/deg; + + //Theta upper lim should be Thetamax + ThetaRocking + //double ThetaUpper = ThetaMax + ThetaRocking; + + while(true){ + //sort theta and phi until it match requirement + + //Theta = acos(cos(0*deg) - CLHEP::RandFlat::shoot(1)*(cos(0*deg) - cos(ThetaUpper*deg))); + + // sort random theta and phi + Theta = acos(1 - CLHEP::RandFlat::shoot(1)*2); + Phi = CLHEP::RandFlat::shoot(1)*360*deg; + + // compute xyz of the particle + G4ThreeVector Dir(sin(Theta)*cos(Phi),sin(Theta)*sin(Phi),cos(Theta)); + + + //rotate the event into the galactic coordinate system + Sky.OrientDirection(m_NextEmission, Dir); + //For some reason Z of particle dir needs to be *-1 #reverseingeniering + Dir[2] = Dir[2]*-1; + + if (m_EarthOccultation_InverseCut){ + if (EarthZenith.angle(Dir)/deg >ThetaMax ){ + break; + } + + } else { + if (EarthZenith.angle(Dir)/deg <=ThetaMax ){ + break; + } + } + } + + + + + } else{ merr<<"the time is not in the ori time range : "<GetCurrentRun().GetSkyOrientationReference(); + + //celestial north pole in galactic coordinates l b is 122.93� and 27.13� + //see :https://lambda.gsfc.nasa.gov/product/about/pol_convention.html + //for some reason we need to multiply by -1 for Z #reverseingienering + G4ThreeVector CelestNorthPole(cos(27.13*deg)*cos(122.93*deg) , cos(27.13*deg)*sin(122.93*deg) , -1*sin(27.13*deg)); + + //convert the north pole in local coordinates + Sky.OrientDirectionInvers(m_NextEmission, CelestNorthPole); + cout<<"celest vector : "< Date: Mon, 11 Nov 2024 15:52:19 +0100 Subject: [PATCH 02/25] Update MCOrientation.cc --- src/cosima/src/MCOrientation.cc | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/cosima/src/MCOrientation.cc b/src/cosima/src/MCOrientation.cc index 61a917b9..d3a03f8d 100644 --- a/src/cosima/src/MCOrientation.cc +++ b/src/cosima/src/MCOrientation.cc @@ -236,10 +236,12 @@ bool MCOrientation::Parse(const MTokenizer& Tokenizer) } else if (Tokenizer.IsTokenAt(3, "File") == true) { - if (Tokenizer.GetNTokens() != 6) { + if (Tokenizer.GetNTokens() != 6 ) { mlog<<" *** Error *** You need exactly 6 tokens to for an orientation read from file"< 0) { m_IsOriented = true; @@ -392,6 +397,7 @@ bool MCOrientation::Read(MString FileName) } + return true; } From c171170f41d7c78f1bf459840aa9d8b90d4203a2 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 11 Nov 2024 15:52:52 +0100 Subject: [PATCH 03/25] Update MCParameterFile.cc --- src/cosima/src/MCParameterFile.cc | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/cosima/src/MCParameterFile.cc b/src/cosima/src/MCParameterFile.cc index 02e9f530..5a7748db 100644 --- a/src/cosima/src/MCParameterFile.cc +++ b/src/cosima/src/MCParameterFile.cc @@ -324,12 +324,8 @@ bool MCParameterFile::Parse() m_DecayMode = c_DecayModeIgnore; } else if (Mode == "activationbuildup") { m_DecayMode = c_DecayModeActivationBuildUp; - Typo(i, "Cannot use this activation option for this current MEGAlib branch. Please use the option Buildup"); - return false; } else if (Mode == "activationdelayeddecay") { m_DecayMode = c_DecayModeActivationDelayedDecay; - Typo(i, "Cannot use this activation option for this current MEGAlib branch. Please use the option Buildup"); - return false; } else { Typo(i, "Cannot parse token DecayMode correctly:" " Unknown Decay mode!"); @@ -865,14 +861,15 @@ bool MCParameterFile::Parse() if (O.Parse(*T) == false) { Typo(i, "Cannot parse token \"OrientationSky\" correctly"); return false; - } + } + // The sky can only be rotated if in Galactic coordinates - if (O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { + if ( O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { Typo(i, "\"OrientationSky\" can only have an orientation in Galactic coordinates!"); return false; } Run->SetSkyOrientation(O); - + } else { Typo(i, "Cannot parse token \"OrientationSky\" correctly: Number of tokens is not correct!"); return false; @@ -1449,14 +1446,16 @@ bool MCParameterFile::Parse() } } else if (Type == "farfieldearthoccultation" ) { - if (T->GetNTokens() >= 3) { + if (T->GetNTokens() >= 4) { Source->SetBeamType(MCSource::c_FarField, MCSource::c_FarFieldEarthOccultation); - + double Theta = 0.0; // here we don't use theta + bool InverseCut = T->GetTokenAtAsBoolean(2); + Source->SetEarthOccultation(Theta,InverseCut); mdebug<<"Using beam: "<GetNTokens() == 5) { + Source->SetPolarizationType(MCSource::c_PolarizationGalactic); + Source->SetPolarizationDegree(T->GetTokenAtAsDouble(3)); + Source->SetPolarization(T->GetTokenAtAsDouble(4)*deg); + } else { + Typo(i, "Cannot parse token \"Polarization galactic\" correctly: Number of tokens is not correct!"); + return false; + } + + } else if (Type == "relativex") { if (T->GetNTokens() == 5) { Source->SetPolarizationType(MCSource::c_PolarizationRelativeX); Source->SetPolarizationDegree(T->GetTokenAtAsDouble(3)); From 511fbf01504eb42326ed77ce483cf07c38ddd359 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 11 Nov 2024 15:53:29 +0100 Subject: [PATCH 04/25] Update MCSource.hh --- src/cosima/inc/MCSource.hh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/cosima/inc/MCSource.hh b/src/cosima/inc/MCSource.hh index 64655266..cf487f9b 100644 --- a/src/cosima/inc/MCSource.hh +++ b/src/cosima/inc/MCSource.hh @@ -184,6 +184,7 @@ public: bool SetPolarization(double PolarizationParam1 = c_Invalid, double PolarizationParam2 = c_Invalid, double PolarizationParam3 = c_Invalid); + /// Return the polarization vector G4ThreeVector GetPolarizationVector() const { return m_Polarization; } @@ -432,6 +433,8 @@ public: static const int c_PolarizationRandom; /// Id of the polarization being in absolute coordinates static const int c_PolarizationAbsolute; + /// Id of the polarization being in Galactic coordinates + static const int c_PolarizationGalactic; /// Id of the polarization being calculated relative to particle direction and x-axis static const int c_PolarizationRelativeX; /// Id of the polarization being calculated relative to particle direction and y-axis From 9d33573177553da462b002ee5e5ef094e3eea6c9 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 11 Nov 2024 16:02:00 +0100 Subject: [PATCH 05/25] Update MCSource.cc --- src/cosima/src/MCSource.cc | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/cosima/src/MCSource.cc b/src/cosima/src/MCSource.cc index 2f3495c1..aa32184e 100644 --- a/src/cosima/src/MCSource.cc +++ b/src/cosima/src/MCSource.cc @@ -14,9 +14,7 @@ * and all its terms. * */ - - - + // Cosima: #include "MCCommon.hh" #include "MCSource.hh" @@ -272,7 +270,6 @@ void MCSource::Initialize() m_PolarizationParam2 = c_Invalid; m_PolarizationParam3 = c_Invalid; m_PolarizationDegree = 0.0; - m_UseFarFieldTransmissionProbability = false; m_UseEarthOccultation = false; m_ThetaMaxEarthOccultation = 0.0; @@ -414,7 +411,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) - //get the Earth aspect information from the orientation file + //get the Earth aspect information from the orientation file const MCOrientation& Sky = MCRunManager::GetMCRunManager()->GetCurrentRun().GetSkyOrientationReference(); double Alt,Lat,Lon=0; if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ @@ -431,7 +428,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) //Compute the theta max (angle between the particle direction and Earth-based zenith) //using the altitude from ori file - int R_Earth = 6378; //km + double R_Earth = 6378; //km ThetaMax = ( c_Pi - asin( R_Earth/(R_Earth+Alt/km) ) )/deg; } @@ -3760,7 +3757,7 @@ bool MCSource::PerformOrientation(G4ThreeVector& Direction) // This reorientation can only happen is both are of the same coordinate system Sky.OrientDirectionInvers(m_NextEmission, Direction); } - } else if ( m_Orientation.GetCoordinateSystem() == MCOrientationCoordinateSystem::c_Galactic && Sky.GetCoordinateSystem() == MCOrientationCoordinateSystem::c_Galactic ) { + } else if ( m_Orientation.GetCoordinateSystem() == MCOrientationCoordinateSystem::c_Galactic && Sky.GetCoordinateSystem() == MCOrientationCoordinateSystem::c_Galactic) { if (m_CoordinateSystem != c_FarField) { mout< Date: Mon, 11 Nov 2024 16:04:35 +0100 Subject: [PATCH 06/25] Update MCParameterFile.cc --- src/cosima/src/MCParameterFile.cc | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/cosima/src/MCParameterFile.cc b/src/cosima/src/MCParameterFile.cc index 5a7748db..eb747092 100644 --- a/src/cosima/src/MCParameterFile.cc +++ b/src/cosima/src/MCParameterFile.cc @@ -324,8 +324,12 @@ bool MCParameterFile::Parse() m_DecayMode = c_DecayModeIgnore; } else if (Mode == "activationbuildup") { m_DecayMode = c_DecayModeActivationBuildUp; + Typo(i, "Cannot use this activation option for this current MEGAlib branch. Please use the option Buildup"); + return false; } else if (Mode == "activationdelayeddecay") { m_DecayMode = c_DecayModeActivationDelayedDecay; + Typo(i, "Cannot use this activation option for this current MEGAlib branch. Please use the option Buildup"); + return false; } else { Typo(i, "Cannot parse token DecayMode correctly:" " Unknown Decay mode!"); @@ -864,7 +868,7 @@ bool MCParameterFile::Parse() } // The sky can only be rotated if in Galactic coordinates - if ( O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { + if (O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { Typo(i, "\"OrientationSky\" can only have an orientation in Galactic coordinates!"); return false; } From 2b7c22a491de2e4d5919fb87920e2d018ff745cf Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 11 Nov 2024 16:05:38 +0100 Subject: [PATCH 07/25] Update MCOrientation.cc --- src/cosima/src/MCOrientation.cc | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/cosima/src/MCOrientation.cc b/src/cosima/src/MCOrientation.cc index d3a03f8d..61a917b9 100644 --- a/src/cosima/src/MCOrientation.cc +++ b/src/cosima/src/MCOrientation.cc @@ -236,12 +236,10 @@ bool MCOrientation::Parse(const MTokenizer& Tokenizer) } else if (Tokenizer.IsTokenAt(3, "File") == true) { - if (Tokenizer.GetNTokens() != 6 ) { + if (Tokenizer.GetNTokens() != 6) { mlog<<" *** Error *** You need exactly 6 tokens to for an orientation read from file"< 0) { m_IsOriented = true; @@ -397,7 +392,6 @@ bool MCOrientation::Read(MString FileName) } - return true; } From be7081026bd41073326256e8a546f4f96990f374 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 11 Nov 2024 16:07:02 +0100 Subject: [PATCH 08/25] Update MCParameterFile.cc --- src/cosima/src/MCParameterFile.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/cosima/src/MCParameterFile.cc b/src/cosima/src/MCParameterFile.cc index eb747092..89a0541f 100644 --- a/src/cosima/src/MCParameterFile.cc +++ b/src/cosima/src/MCParameterFile.cc @@ -865,10 +865,9 @@ bool MCParameterFile::Parse() if (O.Parse(*T) == false) { Typo(i, "Cannot parse token \"OrientationSky\" correctly"); return false; - } - + } // The sky can only be rotated if in Galactic coordinates - if (O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { + if (O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { Typo(i, "\"OrientationSky\" can only have an orientation in Galactic coordinates!"); return false; } From 2339d6c1cc8e505d95230d8a2b96348abd33b316 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 14:47:38 +0100 Subject: [PATCH 09/25] Update MCSource.cc --- src/cosima/src/MCSource.cc | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/cosima/src/MCSource.cc b/src/cosima/src/MCSource.cc index aa32184e..142be055 100644 --- a/src/cosima/src/MCSource.cc +++ b/src/cosima/src/MCSource.cc @@ -3660,12 +3660,12 @@ bool MCSource::GeneratePolarization(G4GeneralParticleSource* Gun) //celestial north pole in galactic coordinates l b is 122.93� and 27.13� //see :https://lambda.gsfc.nasa.gov/product/about/pol_convention.html - //for some reason we need to multiply by -1 for Z #reverseingienering + //because of left handed megalib convention we need to multiply by -1 for Z #reverseingienering G4ThreeVector CelestNorthPole(cos(27.13*deg)*cos(122.93*deg) , cos(27.13*deg)*sin(122.93*deg) , -1*sin(27.13*deg)); //convert the north pole in local coordinates Sky.OrientDirectionInvers(m_NextEmission, CelestNorthPole); - cout<<"celest vector : "< Will use zero polarization!"< Date: Mon, 6 Jan 2025 14:51:31 +0100 Subject: [PATCH 10/25] Update MCParameterFile.cc --- src/cosima/src/MCParameterFile.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosima/src/MCParameterFile.cc b/src/cosima/src/MCParameterFile.cc index 89a0541f..3e52b986 100644 --- a/src/cosima/src/MCParameterFile.cc +++ b/src/cosima/src/MCParameterFile.cc @@ -1453,7 +1453,7 @@ bool MCParameterFile::Parse() Source->SetBeamType(MCSource::c_FarField, MCSource::c_FarFieldEarthOccultation); double Theta = 0.0; // here we don't use theta - bool InverseCut = T->GetTokenAtAsBoolean(2); + bool InverseCut = T->GetTokenAtAsBoolean(3); Source->SetEarthOccultation(Theta,InverseCut); mdebug<<"Using beam: "< Date: Mon, 6 Jan 2025 14:56:27 +0100 Subject: [PATCH 11/25] Update MCSource.hh --- src/cosima/inc/MCSource.hh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/cosima/inc/MCSource.hh b/src/cosima/inc/MCSource.hh index cf487f9b..37eaf937 100644 --- a/src/cosima/inc/MCSource.hh +++ b/src/cosima/inc/MCSource.hh @@ -183,8 +183,7 @@ public: /// Return true, if the position vector could be set correctly bool SetPolarization(double PolarizationParam1 = c_Invalid, double PolarizationParam2 = c_Invalid, - double PolarizationParam3 = c_Invalid); - + double PolarizationParam3 = c_Invalid); /// Return the polarization vector G4ThreeVector GetPolarizationVector() const { return m_Polarization; } From 9654207f00b3e07641a252df013e568acd9f403c Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 14:57:14 +0100 Subject: [PATCH 12/25] Update MCSource.hh --- src/cosima/inc/MCSource.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosima/inc/MCSource.hh b/src/cosima/inc/MCSource.hh index 37eaf937..4aa15217 100644 --- a/src/cosima/inc/MCSource.hh +++ b/src/cosima/inc/MCSource.hh @@ -183,7 +183,7 @@ public: /// Return true, if the position vector could be set correctly bool SetPolarization(double PolarizationParam1 = c_Invalid, double PolarizationParam2 = c_Invalid, - double PolarizationParam3 = c_Invalid); + double PolarizationParam3 = c_Invalid); /// Return the polarization vector G4ThreeVector GetPolarizationVector() const { return m_Polarization; } From 54579c13f054dfcc99e59db517e43da4e4e99e93 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:00:05 +0100 Subject: [PATCH 13/25] Update MCParameterFile.cc --- src/cosima/src/MCParameterFile.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cosima/src/MCParameterFile.cc b/src/cosima/src/MCParameterFile.cc index 3e52b986..035f416b 100644 --- a/src/cosima/src/MCParameterFile.cc +++ b/src/cosima/src/MCParameterFile.cc @@ -329,7 +329,7 @@ bool MCParameterFile::Parse() } else if (Mode == "activationdelayeddecay") { m_DecayMode = c_DecayModeActivationDelayedDecay; Typo(i, "Cannot use this activation option for this current MEGAlib branch. Please use the option Buildup"); - return false; + return false; } else { Typo(i, "Cannot parse token DecayMode correctly:" " Unknown Decay mode!"); @@ -865,7 +865,7 @@ bool MCParameterFile::Parse() if (O.Parse(*T) == false) { Typo(i, "Cannot parse token \"OrientationSky\" correctly"); return false; - } + } // The sky can only be rotated if in Galactic coordinates if (O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { Typo(i, "\"OrientationSky\" can only have an orientation in Galactic coordinates!"); From 831d57173507cec9c5b613c4dd1a14552476a979 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:01:44 +0100 Subject: [PATCH 14/25] Update MCParameterFile.cc --- src/cosima/src/MCParameterFile.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cosima/src/MCParameterFile.cc b/src/cosima/src/MCParameterFile.cc index 035f416b..4cdb2ea6 100644 --- a/src/cosima/src/MCParameterFile.cc +++ b/src/cosima/src/MCParameterFile.cc @@ -865,14 +865,14 @@ bool MCParameterFile::Parse() if (O.Parse(*T) == false) { Typo(i, "Cannot parse token \"OrientationSky\" correctly"); return false; - } + } // The sky can only be rotated if in Galactic coordinates if (O.GetCoordinateSystem() != MCOrientationCoordinateSystem::c_Galactic && O.IsOriented() == true) { Typo(i, "\"OrientationSky\" can only have an orientation in Galactic coordinates!"); return false; } Run->SetSkyOrientation(O); - + } else { Typo(i, "Cannot parse token \"OrientationSky\" correctly: Number of tokens is not correct!"); return false; From c76abf7d31902b54b8ab9dd664459f128d516668 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:07:12 +0100 Subject: [PATCH 15/25] Update MCSource.cc --- src/cosima/src/MCSource.cc | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/src/cosima/src/MCSource.cc b/src/cosima/src/MCSource.cc index 142be055..fd554968 100644 --- a/src/cosima/src/MCSource.cc +++ b/src/cosima/src/MCSource.cc @@ -14,7 +14,8 @@ * and all its terms. * */ - + + // Cosima: #include "MCCommon.hh" #include "MCSource.hh" @@ -270,6 +271,7 @@ void MCSource::Initialize() m_PolarizationParam2 = c_Invalid; m_PolarizationParam3 = c_Invalid; m_PolarizationDegree = 0.0; + m_UseFarFieldTransmissionProbability = false; m_UseEarthOccultation = false; m_ThetaMaxEarthOccultation = 0.0; @@ -411,7 +413,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) - //get the Earth aspect information from the orientation file + //get the Earth aspect information from the orientation file const MCOrientation& Sky = MCRunManager::GetMCRunManager()->GetCurrentRun().GetSkyOrientationReference(); double Alt,Lat,Lon=0; if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ @@ -439,16 +441,6 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) // convert the Earth galactic l,b coordinate into cartesian representation G4ThreeVector EarthZenith(cos(Lat)*cos(Lon) , cos(Lat)*sin(Lon) , sin(Lat)); - - - //compute zpointing of spacecraft - //double Xlat=0; - //double Xlon=0; - //double Zlat=0; - //double Zlon=0; - //Sky.GetOrientation(m_NextEmission,Xlat,Xlon,Zlat,Zlon); - //G4ThreeVector SpacecraftZaxis(cos(Zlat)*cos(Zlon) , cos(Zlat)*sin(Zlon) , sin(Zlat)); - //G4ThreeVector SpacecraftXaxis(cos(Xlat)*cos(Xlon) , cos(Xlat)*sin(Xlon) , sin(Xlat)); //rotate the oriented coordinate system of the event into the local coordinate system //Sky.OrientDirectionInvers(m_NextEmission, Dir); @@ -460,9 +452,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) //cout<< "particle Dir: "< Date: Mon, 6 Jan 2025 15:17:00 +0100 Subject: [PATCH 16/25] Update MCSource.cc --- src/cosima/src/MCSource.cc | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/cosima/src/MCSource.cc b/src/cosima/src/MCSource.cc index fd554968..f17ad6b9 100644 --- a/src/cosima/src/MCSource.cc +++ b/src/cosima/src/MCSource.cc @@ -271,7 +271,6 @@ void MCSource::Initialize() m_PolarizationParam2 = c_Invalid; m_PolarizationParam3 = c_Invalid; m_PolarizationDegree = 0.0; - m_UseFarFieldTransmissionProbability = false; m_UseEarthOccultation = false; m_ThetaMaxEarthOccultation = 0.0; @@ -413,7 +412,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) - //get the Earth aspect information from the orientation file + //get the Earth aspect information from the orientation file const MCOrientation& Sky = MCRunManager::GetMCRunManager()->GetCurrentRun().GetSkyOrientationReference(); double Alt,Lat,Lon=0; if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ @@ -441,7 +440,7 @@ bool MCSource::GenerateParticles(G4GeneralParticleSource* ParticleGun) // convert the Earth galactic l,b coordinate into cartesian representation G4ThreeVector EarthZenith(cos(Lat)*cos(Lon) , cos(Lat)*sin(Lon) , sin(Lat)); - + //rotate the oriented coordinate system of the event into the local coordinate system //Sky.OrientDirectionInvers(m_NextEmission, Dir); Sky.OrientDirection(m_NextEmission, Dir); @@ -3117,8 +3116,8 @@ bool MCSource::GeneratePosition(G4GeneralParticleSource* Gun) //get the Earth aspect information from the orientation file const MCOrientation& Sky = MCRunManager::GetMCRunManager()->GetCurrentRun().GetSkyOrientationReference(); double Alt,Lat,Lon=0; - //double xLat,xLon,zLat,zLon=0; - if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ //&& Sky.GetOrientation(m_NextEmission,xLat,xLon,zLat,zLon)==true ){ + + if (Sky.GetEarthCoordinate(m_NextEmission, Alt, Lat, Lon)==true){ double R_Earth = 6378; //km @@ -3128,18 +3127,9 @@ bool MCSource::GeneratePosition(G4GeneralParticleSource* Gun) // convert the Earth galactic l,b coordinate into cartesian representation G4ThreeVector EarthZenith(cos(Lat)*cos(Lon) , cos(Lat)*sin(Lon) , sin(Lat)); - // convert the spacecraft galacticz pointing l,b coordinate into cartesian representation - //G4ThreeVector SpacecraftZenith(cos(zLat)*cos(zLon) , cos(zLat)*sin(zLon) , sin(zLat)); - - //double ThetaRocking = EarthZenith.angle(SpacecraftZenith)/deg; - - //Theta upper lim should be Thetamax + ThetaRocking - //double ThetaUpper = ThetaMax + ThetaRocking; while(true){ - //sort theta and phi until it match requirement - - //Theta = acos(cos(0*deg) - CLHEP::RandFlat::shoot(1)*(cos(0*deg) - cos(ThetaUpper*deg))); + //sort theta and phi until it match requirement // sort random theta and phi Theta = acos(1 - CLHEP::RandFlat::shoot(1)*2); @@ -3151,7 +3141,8 @@ bool MCSource::GeneratePosition(G4GeneralParticleSource* Gun) //rotate the event into the galactic coordinate system Sky.OrientDirection(m_NextEmission, Dir); - //For some reason Z of particle dir needs to be *-1 #reverseingeniering + //Due to the left-handed galactic system used in MEGAlib + //the Z of particle dir needs to be *-1 #reverseingeniering Dir[2] = Dir[2]*-1; if (m_EarthOccultation_InverseCut){ @@ -3649,7 +3640,7 @@ bool MCSource::GeneratePolarization(G4GeneralParticleSource* Gun) //get the Earth aspect information from the orientation file const MCOrientation& Sky = MCRunManager::GetMCRunManager()->GetCurrentRun().GetSkyOrientationReference(); - //celestial north pole in galactic coordinates l b is 122.93� and 27.13� + //celestial north pole in galactic coordinates l b is 122.93 deg and 27.13 deg //see :https://lambda.gsfc.nasa.gov/product/about/pol_convention.html //because of left handed megalib convention we need to multiply by -1 for Z #reverseingienering G4ThreeVector CelestNorthPole(cos(27.13*deg)*cos(122.93*deg) , cos(27.13*deg)*sin(122.93*deg) , -1*sin(27.13*deg)); @@ -3678,7 +3669,7 @@ bool MCSource::GeneratePolarization(G4GeneralParticleSource* Gun) cos(m_PolarizationParam1)*px[1]+sin(m_PolarizationParam1)*py[1], cos(m_PolarizationParam1)*px[2]+sin(m_PolarizationParam1)*py[2]); - //m_Polarization[2] = m_Polarization[2]*-1; + m_Polarization = m_Polarization.unit(); //cout<<"polarization vector : "< Will use zero polarization!"< Date: Mon, 6 Jan 2025 15:27:27 +0100 Subject: [PATCH 17/25] Create Polarization_IAUconvention --- resource/examples/advanced/Polarization_IAUconvention | 1 + 1 file changed, 1 insertion(+) create mode 100644 resource/examples/advanced/Polarization_IAUconvention diff --git a/resource/examples/advanced/Polarization_IAUconvention b/resource/examples/advanced/Polarization_IAUconvention new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/resource/examples/advanced/Polarization_IAUconvention @@ -0,0 +1 @@ + From 9f38f8d3ab3661e8298234e774b53b64a3e1d3e9 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:28:25 +0100 Subject: [PATCH 18/25] Delete resource/examples/advanced/Polarization_IAUconvention --- resource/examples/advanced/Polarization_IAUconvention | 1 - 1 file changed, 1 deletion(-) delete mode 100644 resource/examples/advanced/Polarization_IAUconvention diff --git a/resource/examples/advanced/Polarization_IAUconvention b/resource/examples/advanced/Polarization_IAUconvention deleted file mode 100644 index 8b137891..00000000 --- a/resource/examples/advanced/Polarization_IAUconvention +++ /dev/null @@ -1 +0,0 @@ - From 94ef288854594aec4fa6bd345318cc37d32b511e Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:29:23 +0100 Subject: [PATCH 19/25] Create Polarized.source --- .../advanced/Polarization_IAUconvention/Polarized.source | 1 + 1 file changed, 1 insertion(+) create mode 100644 resource/examples/advanced/Polarization_IAUconvention/Polarized.source diff --git a/resource/examples/advanced/Polarization_IAUconvention/Polarized.source b/resource/examples/advanced/Polarization_IAUconvention/Polarized.source new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/resource/examples/advanced/Polarization_IAUconvention/Polarized.source @@ -0,0 +1 @@ + From 0b121c34bfcad8f8a5ecf4f8f43f88a2f84dcf40 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:36:56 +0100 Subject: [PATCH 20/25] Update Polarized.source --- .../Polarized.source | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/resource/examples/advanced/Polarization_IAUconvention/Polarized.source b/resource/examples/advanced/Polarization_IAUconvention/Polarized.source index 8b137891..e2151459 100644 --- a/resource/examples/advanced/Polarization_IAUconvention/Polarized.source +++ b/resource/examples/advanced/Polarization_IAUconvention/Polarized.source @@ -1 +1,21 @@ +Geometry $(MEGALIB)/resource/examples/geomega/cosiballoon/COSIBalloon.9Detector.geo.setup +StoreSimulationInfo init-only +PhysicsListEM LivermorePol + +Run SpaceSim +SpaceSim.FileName PolarizedSource +SpaceSim.Time 10.0 +SpaceSim.OrientationSky Galactic File NoLoop test.ori + + +SpaceSim.Source AlbedoPhotons +AlbedoPhotons.ParticleType 1 +AlbedoPhotons.Beam FarFieldPointSource 0 0 +AlbedoPhotons.Orientation Galactic Fixed -48.90 176.96 +AlbedoPhotons.Spectrum Mono 511 +AlbedoPhotons.Flux 100 + +#fst number is the polarization fraction (0 to 1) and the second is the polarization angle (deg) +#following the IAU convention : https://lambda.gsfc.nasa.gov/product/about/pol_convention.html +AlbedoPhotons.Polarization galactic 1 45 From 8ad6d64ed6b9e0e86cea096832438e1722635df6 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:37:48 +0100 Subject: [PATCH 21/25] Create test.ori --- resource/examples/advanced/Polarization_IAUconvention/test.ori | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 resource/examples/advanced/Polarization_IAUconvention/test.ori diff --git a/resource/examples/advanced/Polarization_IAUconvention/test.ori b/resource/examples/advanced/Polarization_IAUconvention/test.ori new file mode 100644 index 00000000..2434eadc --- /dev/null +++ b/resource/examples/advanced/Polarization_IAUconvention/test.ori @@ -0,0 +1,3 @@ +Type OrientationsGalactic +OG 0 -60.18 96.33 27.13 122.93 +OG 10 -60.18 96.33 27.13 122.93 From 5d2e64486e54cac8158f90287062e2a482845cf9 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Thu, 9 Jan 2025 10:48:20 +0100 Subject: [PATCH 22/25] Rename EarthOccultation.source to EarthOccultation_pointsource.source --- ...arthOccultation.source => EarthOccultation_pointsource.source} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename resource/examples/advanced/EarthOccultation/{EarthOccultation.source => EarthOccultation_pointsource.source} (100%) diff --git a/resource/examples/advanced/EarthOccultation/EarthOccultation.source b/resource/examples/advanced/EarthOccultation/EarthOccultation_pointsource.source similarity index 100% rename from resource/examples/advanced/EarthOccultation/EarthOccultation.source rename to resource/examples/advanced/EarthOccultation/EarthOccultation_pointsource.source From c49c256a3683b53e7c434db0b01b4ef8949ec413 Mon Sep 17 00:00:00 2001 From: GallegoSav <123730578+GallegoSav@users.noreply.github.com> Date: Thu, 3 Jul 2025 13:40:39 +0200 Subject: [PATCH 23/25] Update SensitivityOptimizer.cxx --- src/addon/SensitivityOptimizer.cxx | 83 ++++++++++++++++++++++++++---- 1 file changed, 73 insertions(+), 10 deletions(-) diff --git a/src/addon/SensitivityOptimizer.cxx b/src/addon/SensitivityOptimizer.cxx index af8f4607..9b9af0c0 100644 --- a/src/addon/SensitivityOptimizer.cxx +++ b/src/addon/SensitivityOptimizer.cxx @@ -1530,7 +1530,7 @@ bool SensitivityOptimizer::ParseCommandLine(int argc, char** argv) m_MimrecSettings.Read(ConfigurationFile); m_EventSelector.SetSettings(&m_MimrecSettings); m_EventSelector.SetGeometry(&m_Geometry); - + cout<<"Coordinate system : "< Date: Wed, 29 Oct 2025 11:27:55 +0100 Subject: [PATCH 24/25] Update MResponseImagingBinnedMode.cxx --- src/response/src/MResponseImagingBinnedMode.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/response/src/MResponseImagingBinnedMode.cxx b/src/response/src/MResponseImagingBinnedMode.cxx index 43abee25..5d9bbc5f 100644 --- a/src/response/src/MResponseImagingBinnedMode.cxx +++ b/src/response/src/MResponseImagingBinnedMode.cxx @@ -119,8 +119,8 @@ MString MResponseImagingBinnedMode::Options() out<<" dmin: minimum distance (default: 0 cm)"< Date: Wed, 29 Oct 2025 11:28:41 +0100 Subject: [PATCH 25/25] Update MResponsePolarizationBinnedMode.cxx --- src/response/src/MResponsePolarizationBinnedMode.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/response/src/MResponsePolarizationBinnedMode.cxx b/src/response/src/MResponsePolarizationBinnedMode.cxx index 9a525920..7eebec11 100644 --- a/src/response/src/MResponsePolarizationBinnedMode.cxx +++ b/src/response/src/MResponsePolarizationBinnedMode.cxx @@ -124,8 +124,8 @@ MString MResponsePolarizationBinnedMode::Options() out<<" dbins: number of distance bins between min and max distance (default: 1)"<