Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: single-hadron SIDIS kinematics #295

Merged
merged 8 commits into from
Oct 23, 2024
Merged

feat: single-hadron SIDIS kinematics #295

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
0e8659b
refactor: move common physics tools to their own source files
c-dilks Oct 23, 2024
b65cc01
feat: algorithm and validator
c-dilks Oct 23, 2024
3286db2
fix: algorith name in `Config.yaml`
c-dilks Oct 23, 2024
86ebb81
fix: set result bank rows and filter
c-dilks Oct 23, 2024
451d1f8
fix: only build `physics::tools` if ROOT is found
c-dilks Oct 23, 2024
bececa1
fix: force `pdg` to be an `int`
c-dilks Oct 23, 2024
9343b03
fix: yet another `hipo::structure` ubsan suppression
c-dilks Oct 23, 2024
e7bf918
fix(?): set `xi_h` range to [-1,1]
c-dilks Oct 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions CODEOWNERS
View file
Open in desktop
Validating CODEOWNERS rules …
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,14 @@
#################################################################################

* @c-dilks
src/iguana/algorithms/clas12/EventBuilderFilter/* @c-dilks
src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar
src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3
src/iguana/algorithms/clas12/MomentumCorrection/* @RichCap @c-dilks
src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12
src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3
src/iguana/algorithms/clas12/SectorFinder/* @rtysonCLAS12
src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar
src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12
src/iguana/algorithms/physics/DihadronKinematics/* @c-dilks
src/iguana/algorithms/physics/InclusiveKinematics/* @c-dilks
src/iguana/algorithms/physics/SingleHadronKinematics/* @c-dilks
src/iguana/services/YAMLReader.* @rtysonCLAS12 @c-dilks
src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3
src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3
1 change: 1 addition & 0 deletions meson/ubsan.supp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,4 @@ alignment:hipo::structure::getFloatAt
alignment:hipo::structure::getShortAt
alignment:hipo::structure::getDoubleAt
alignment:hipo::structure::putDoubleAt
alignment:hipo::structure::putIntAt
26 changes: 24 additions & 2 deletions src/iguana/algorithms/meson.build
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,15 @@ algo_dict = [
'algorithm_needs_ROOT': true,
'test_args': {'banks': [ 'REC::Particle', 'RUN::config' ]},
},
{
'name': 'physics::SingleHadronKinematics',
'algorithm_needs_ROOT': true,
'has_action_yaml': false, # FIXME: needs a vector action function
'test_args': {
'banks': [ 'REC::Particle', 'RUN::config' ],
'prerequisites': [ 'physics::InclusiveKinematics' ],
},
},
{
'name': 'physics::DihadronKinematics',
'algorithm_needs_ROOT': true,
Expand All @@ -92,8 +101,21 @@ algo_dict = [
]

# make lists of objects to build; inclusion depends on whether ROOT is needed or not, and if we have ROOT
algo_sources = [ 'Algorithm.cc', 'AlgorithmFactory.cc', 'AlgorithmSequence.cc' ]
algo_headers = [ 'Algorithm.h', 'AlgorithmBoilerplate.h', 'TypeDefs.h', 'AlgorithmSequence.h' ]
algo_sources = [
'Algorithm.cc',
'AlgorithmFactory.cc',
'AlgorithmSequence.cc',
]
algo_headers = [
'Algorithm.h',
'AlgorithmBoilerplate.h',
'TypeDefs.h',
'AlgorithmSequence.h',
]
if ROOT_dep.found()
algo_sources += [ 'physics/Tools.cc' ]
algo_headers += [ 'physics/Tools.h' ]
endif
vdor_sources = [ 'Validator.cc' ]
vdor_headers = [ 'Validator.h' ]
algo_configs = []
Expand Down
71 changes: 9 additions & 62 deletions src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,41 +145,41 @@ namespace iguana::physics {
double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R());

// calculate phiH
double phiH = PlaneAngle(
double phiH = tools::PlaneAngle(
p_q.Vect(),
p_beam.Vect(),
p_q.Vect(),
p_Ph.Vect()).value_or(UNDEF);
p_Ph.Vect()).value_or(tools::UNDEF);

// calculate PhiR
double phiR = UNDEF;
double phiR = tools::UNDEF;
switch(m_phi_r_method) {
case e_RT_via_covariant_kT:
{
for(auto& had : {&had_a, &had_b}) {
had->z = p_target.Dot(had->p) / p_target.Dot(p_q);
had->p_perp = RejectVector(had->p.Vect(), p_q.Vect());
had->p_perp = tools::RejectVector(had->p.Vect(), p_q.Vect());
}
if(had_a.p_perp.has_value() && had_b.p_perp.has_value()) {
auto RT = (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()) / (had_a.z + had_b.z);
phiR = PlaneAngle(
phiR = tools::PlaneAngle(
p_q.Vect(),
p_beam.Vect(),
p_q.Vect(),
RT).value_or(UNDEF);
RT).value_or(tools::UNDEF);
}
break;
}
}

// calculate theta
double theta = UNDEF;
double theta = tools::UNDEF;
switch(m_theta_method) {
case e_hadron_a:
{
theta = VectorAngle(
theta = tools::VectorAngle(
boost__dih(had_a.p).Vect(),
p_Ph.Vect()).value_or(UNDEF);
p_Ph.Vect()).value_or(tools::UNDEF);
break;
}
}
Expand Down Expand Up @@ -238,59 +238,6 @@ namespace iguana::physics {

///////////////////////////////////////////////////////////////////////////////

std::optional<double> DihadronKinematics::PlaneAngle(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b,
ROOT::Math::XYZVector const v_c,
ROOT::Math::XYZVector const v_d)
{
auto cross_ab = v_a.Cross(v_b); // A x B
auto cross_cd = v_c.Cross(v_d); // C x D

auto sgn = cross_ab.Dot(v_d); // (A x B) . D
if(!(std::abs(sgn) > 0))
return std::nullopt;
sgn /= std::abs(sgn); // sign of (A x B) . D

auto numer = cross_ab.Dot(cross_cd); // (A x B) . (C x D)
auto denom = cross_ab.R() * cross_cd.R(); // |A x B| * |C x D|
if(!(std::abs(denom) > 0))
return std::nullopt;
return sgn * std::acos(numer / denom);
}

std::optional<ROOT::Math::XYZVector> DihadronKinematics::ProjectVector(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b)
{
auto denom = v_b.Dot(v_b);
if(!(std::abs(denom) > 0))
return std::nullopt;
return v_b * ( v_a.Dot(v_b) / denom );
}

std::optional<ROOT::Math::XYZVector> DihadronKinematics::RejectVector(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b)
{
auto v_c = ProjectVector(v_a, v_b);
if(v_c.has_value())
return v_a - v_c.value();
return std::nullopt;
}

std::optional<double> DihadronKinematics::VectorAngle(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b)
{
double m = v_a.R() * v_b.R();
if(m > 0)
return std::acos(v_a.Dot(v_b) / m);
return std::nullopt;
}

///////////////////////////////////////////////////////////////////////////////

void DihadronKinematics::Stop()
{
}
Expand Down
52 changes: 6 additions & 46 deletions src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#pragma once

#include "iguana/algorithms/Algorithm.h"
#include "iguana/algorithms/physics/Tools.h"
#include <Math/Vector3D.h>
#include <Math/Vector4D.h>

Expand All @@ -18,27 +19,27 @@ namespace iguana::physics {
int pdg_b;
/// @brief @latex{M_h}: Invariant mass of the dihadron
double Mh;
/// @brief @latex{z}: Fraction of energy of fragmenting parton carried by the dihadron
/// @brief @latex{z}: Momentum fraction of the fragmenting parton carried by the dihadron
double z;
/// @brief @latex{M_X(ehhX)}: Missing mass of the dihadron
double MX;
/// @brief @latex{x_F}: Feynman-x of the dihadron
double xF;
/// @brief @latex{\phi_h}: @latex{q}-azimuthal angle between the lepton-scattering plane and the @latex{\vec{q}\times\vec{P}_h} plane;
/// if the value is `DihadronKinematics::UNDEF`, the calculation failed
/// if the value is `tools::UNDEF`, the calculation failed
double phiH;
/// @brief @latex{\phi_R}: @latex{q}-azimuthal angle between the lepton-scattering plane and dihadron plane;
/// if the value is `DihadronKinematics::UNDEF`, the calculation failed
/// if the value is `tools::UNDEF`, the calculation failed
double phiR;
/// @brief @latex{\theta}: the "decay" angle of hadron A in the dihadron rest frame, with respect;
/// @brief @latex{\theta}: The "decay" angle of hadron A in the dihadron rest frame, with respect;
/// to the dihadron momentum direction
double theta;
};

/// @brief_algo Calculate semi-inclusive dihadron kinematic quantities defined in `iguana::physics::DihadronKinematicsVars`
///
/// @begin_doc_algo{physics::DihadronKinematics | Creator}
/// @input_banks{REC::Particle, physics::InclusiveKinematics}
/// @input_banks{REC::Particle, %physics::InclusiveKinematics}
/// @output_banks{%physics::DihadronKinematics}
/// @end_doc
///
Expand Down Expand Up @@ -75,52 +76,11 @@ namespace iguana::physics {
void Run(hipo::banklist& banks) const override;
void Stop() override;

/// a value used when some calculation fails
double const UNDEF{-10000};

/// @brief form dihadrons by pairing hadrons
/// @param particle_bank the particle bank (`REC::Particle`)
/// @returns a list of pairs of hadron rows
std::vector<std::pair<int,int>> PairHadrons(hipo::bank const& particle_bank) const;

/// @brief calculate the angle between two planes
///
/// The two planes are transverse to @latex{\vec{v}_a\times\vec{v}_b} and @latex{\vec{v}_c\times\vec{v}_d}
/// @param v_a vector @latex{\vec{v}_a}
/// @param v_b vector @latex{\vec{v}_b}
/// @param v_c vector @latex{\vec{v}_c}
/// @param v_d vector @latex{\vec{v}_d}
/// @returns the angle between the planes, in radians, if the calculation is successful
static std::optional<double> PlaneAngle(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b,
ROOT::Math::XYZVector const v_c,
ROOT::Math::XYZVector const v_d);

/// @brief projection of one vector onto another
/// @param v_a vector @latex{\vec{v}_a}
/// @param v_b vector @latex{\vec{v}_b}
/// @returns the vector @latex{\vec{v}_a} projected onto vector @latex{\vec{v}_b}, if the calculation is successful
static std::optional<ROOT::Math::XYZVector> ProjectVector(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b);

/// @brief projection of one vector onto the plane transverse to another vector
/// @param v_a vector @latex{\vec{v}_a}
/// @param v_b vector @latex{\vec{v}_b}
/// @returns the vector @latex{\vec{v}_a} projected onto the plane transverse to @latex{\vec{v}_b}, if the calculation is successful
static std::optional<ROOT::Math::XYZVector> RejectVector(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b);

/// @brief calculate the angle between two vectors
/// @param v_a vector @latex{\vec{v}_a}
/// @param v_b vector @latex{\vec{v}_b}
/// @returns the angle between @latex{\vec{v}_a} and @latex{\vec{v}_b}, if the calculation is successful
static std::optional<double> VectorAngle(
ROOT::Math::XYZVector const v_a,
ROOT::Math::XYZVector const v_b);

private:

// banklist indices
Expand Down
2 changes: 1 addition & 1 deletion src/iguana/algorithms/physics/DihadronKinematics/Validator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ namespace iguana::physics {
[](auto const& b, auto const r) { return b.getDouble("z", r); }
},
{
new TH1D("MX_dist", "missing mass M_X [GeV];", n_bins, 0, 4),
new TH1D("MX_dist", "missing mass M_{X} [GeV];", n_bins, 0, 4),
[](auto const& b, auto const r) { return b.getDouble("MX", r); }
},
{
Expand Down
Loading
Loading