diff --git a/src/cosima/inc/MCSource.hh b/src/cosima/inc/MCSource.hh index 35b18cb5..32a8c2a2 100644 --- a/src/cosima/inc/MCSource.hh +++ b/src/cosima/inc/MCSource.hh @@ -567,6 +567,20 @@ public: /// Id of an Omega particle static const int c_Omega; + /// Normalize angles between 0 and π (0 to 180 degrees) + double NormalizeAngle(double angle) const; + + /// Orthographic projection: returns px and py vectors + std::pair OrthographicProjection(const G4ThreeVector& vec, const G4ThreeVector& ref_vector) const; + + /// Stereographic projection: returns px and py vectors + std::pair StereographicProjection(const G4ThreeVector& vec, const G4ThreeVector& ref_vector) const; + + /// Transform the polarization angle between conventions + double TransformPolarizationAngle(double pa_1, const std::pair& projection1, const std::pair& projection2) const; + + /// Example function to demonstrate the projections and angle transformations + bool ProcessPolarizationTransformation(); // protected methods: protected: @@ -592,7 +606,6 @@ protected: /// Perform an orientation of the vector from local into oriented coordinate system bool PerformOrientation(G4ThreeVector& Position, G4ThreeVector& Direction); - // protected members: protected: diff --git a/src/cosima/src/MCSource.cc b/src/cosima/src/MCSource.cc index 73a5f62a..b2fee445 100644 --- a/src/cosima/src/MCSource.cc +++ b/src/cosima/src/MCSource.cc @@ -3626,6 +3626,104 @@ double MCSource::Comptonized(const double Energy, double Alpha, double Epeak) co return pow(Energy, Alpha)*exp(-(Alpha+2)*Energy/Epeak); } +/****************************************************************************** + +* Polarization conventions +*/ + + +// Constants for angle normalization +const double PI = 3.14159265358979323846; + + +// Normalize angles between 0 and π (0 to 180 degrees) +double MCSource::NormalizeAngle(double angle) const { + if (angle < 0) { + angle += PI; + } + return angle; +} + +// Orthographic projection: returns px and py vectors +std::pair MCSource::OrthographicProjection(const G4ThreeVector& vec, const G4ThreeVector& ref_vector) const { + double dot_product = vec.dot(ref_vector); + G4ThreeVector px = ref_vector - dot_product * vec; + px = px.unit(); // Normalize px + G4ThreeVector py = vec.cross(px); + return std::make_pair(px, py); +} + +// Stereographic projection: returns px and py vectors +std::pair MCSource::StereographicProjection(const G4ThreeVector& vec, const G4ThreeVector& ref_vector) const { + G4ThreeVector norm_ref = ref_vector.unit(); + double denom = (vec.z() + 1.0); + G4ThreeVector px(1.0 - (vec.x() * vec.x() - vec.y() * vec.y()) / (denom * denom), + -2.0 * vec.x() * vec.y() / (denom * denom), + -2.0 * vec.x() / denom); + px = px.unit(); // Normalize px + G4ThreeVector py = vec.cross(px); + return std::make_pair(px, py); +} + +// Transform the polarization angle between conventions +ddouble MCSource::TransformPolarizationAngle(double pa_1, const std::pair& projection1, const std::pair& projection2) const { + double cos_pa_1 = std::cos(pa_1); + double sin_pa_1 = std::sin(pa_1); + + // Defining the first projection's basis vectors + G4ThreeVector px1 = projection1.first; + G4ThreeVector py1 = projection1.second; + + // Defining the second projection's basis vectors + G4ThreeVector px2 = projection2.first; + G4ThreeVector py2 = projection2.second; + + // Compute the polarization vector in the original coordinate system + G4ThreeVector pol_vec = px1 * cos_pa_1 + py1 * sin_pa_1; + + // Compute dot products with the new coordinate system + double a = pol_vec.dot(px2); + double b = pol_vec.dot(py2); + + // Calculate the new polarization angle + double pa_2 = std::atan2(b, a); + return NormalizeAngle(pa_2); +} + + +// Example function to demonstrate the projections and angle transformations +bool MCSource::ProcessPolarizationTransformation() { + // Example input: polarization angle in radians (45 degrees) + double pa_1 = PI / 4.0; + + // Reference vector: Z-axis (default) + G4ThreeVector ref_vector(0, 0, 1); + + // Example source direction vector + G4ThreeVector vec(std::sin(PI / 6) * std::cos(PI / 4), + std::sin(PI / 6) * std::sin(PI / 4), + std::cos(PI / 6)); + + // Perform orthographic projection + auto ortho_projection = OrthographicProjection(vec, ref_vector); + + // Perform stereographic projection + auto stereo_projection = StereographicProjection(vec, ref_vector); + + // Transform the polarization angle from orthographic to stereographic + double pa_2_ortho_to_stereo = TransformPolarizationAngle(pa_1, ortho_projection, stereo_projection); + + // Transform the polarization angle from stereographic back to orthographic + double pa_2_stereo_to_ortho = TransformPolarizationAngle(pa_2_ortho_to_stereo, stereo_projection, ortho_projection); + + // Output results (replace this with the appropriate output mechanism for your framework) + mout << "PA_1: " << pa_1 << " radians" << endl; + mout << "PA_2 (Ortho to Stereo): " << pa_2_ortho_to_stereo << " radians" << endl; + mout << "PA_2 (Stereo to Ortho): " << pa_2_stereo_to_ortho << " radians" << endl; + + return true; +} + /* * MCSource.cc: the end... ******************************************************************************/