From 388866883748f43203e4a019225bdb788310869b Mon Sep 17 00:00:00 2001 From: edbarriga Date: Wed, 8 Jul 2026 08:02:23 -0400 Subject: [PATCH 1/4] Updates to Vec_ps_refl and VecPSPlotGenerator to better handle pol info in the beam P4 --- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc | 518 ++++++++++++++---- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h | 12 +- .../AMPTOOLS_DATAIO/VecPsPlotGenerator.cc | 267 ++++++--- .../AMPTOOLS_DATAIO/VecPsPlotGenerator.h | 7 +- src/programs/AmplitudeAnalysis/fit/fit.cc | 156 +++++- 5 files changed, 731 insertions(+), 229 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc index 0394d7074..549c7ab31 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc @@ -1,9 +1,120 @@ - -#include -#include +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// Here are some useful notes to understand how this amplitude works and how to +// use it in your AmpTools configuration file: +// +// PURPOSE: +// Calculates the lineary polarized photoproduction amplitude for a system +// composed of a vector (V) and pseudoscalar (ps) +// +// PHYSICS SUMMARY: +// The amplitude uses the helicity formalism (see [GlueX-doc-4858]): +// a resonance X with total spin J, spin projection M, and partial wave L +// decays into V and ps. The angular distribution is described by +// Wigner D-functions and Clebsch-Gordan coefficients. A barrier factor +// accounts for the centrifugal suppression near threshold. +// +// Beam polarization information is required in two places: +// - The fraction linear polarization. +// - The angle between the polarization angle and the production plane +// +// The default behavior is that the vector decays into two pseudoscalars +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// CONSTRUCTOR ARGUMENTS (for the AmpTools config file): +// We support two formats for optional arguments: a positional input format +// and a key=value format +// +// REQUIRED--the first 5 arguments are required and positional: +// [0] m_j - Total spin (J) of the resonance X (integer >= 0) +// [1] m_m - Spin projection (m) of X; must satisfy |m| <= J +// [2] m_l - Orbital angular momentum (l) between V and ps (integer >= 0) +// [3] m_r - +1 (r)eal part or -1 imaginary part of the amplitude +// [4] m_s - (S)ign of the polarization fraction +1 for (1 + P_gamma) or +// -1 for (1 - P_gamma) +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// OPTIONAL--positional arguments: +// [5] polAngle Beam polarization angle. Degrees will be converted +// to radians +// [6] polFraction Fixed polarization fraction [0, 1], +// OR a path to a .root file containing a TH1D histogram +// of polarization vs beam energy +// [7] histName (required if [6] is a .root file) - Name of the TH1D +// inside the root file +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// OPTIONAL--key=value arguments: +// polAngle= Beam polarization angle +// polFraction= Fixed polarization fraction [0, 1] +// polFile= path to a .root file containing a TH1D histogram +// of polarization vs beam energy +// polHist= Must be used if polFile is used +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// Note: If NO Optional arguments are given, the code assumes that you have +// stored the polarization information in the Px and Py components of the +// beam. This is useful since the acceptance within a run period is assumed +// to be constant and reduces the data set by a factor of 4. +// When constructing your tree save the beam 4-vector as follows: +// Px = polFraction*cos(polAngle) *Make sure the angle is in radians +// Py = polFraction*sin(polAngle) *Make sure the angle is in radians +// Pz = Regular beam Pz +// E = Regular beam energy *Note that Pz = E +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// The following are keywords to turn-on features on the amplitude. They must +// be placed at the end of optional arguments if using the positional format, +// or can be placed anywhere if using the key=value format: +// [X] "omega3pi" V decays will be described as omega -> pi+ pi- pi0 +// [Y] "omegagpi0" V decays will be described as omega -> gamma pi0 +// [Y+1] gHelicity When omegagpi0 is used, one must specify the +// helicity of the bachelor photon (+1 or -1) +// gHelicity=<+1|-1> Alternative key=value format for the photon helicity +/// "nobarrier" Turn off the barrier factors +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// Examples of usage in config file: +// For a 1S m=-1 wave from an omega->3pi that belongs to the imaginary part of +// the amplitude sum and has a negative sign in front of the polarization term, +// you would write: +// If polarization info stored in the beam P4: +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 -1 0 -1 -1 omega3pi +// If polarization angle and fraction are fixed: +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 -1 0 -1 -1 90 .35 omega3pi +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 -1 0 -1 -1 polFraction=0.35 polAngle=90 omega3pi +// If polarization vs E_gamma from histogram: +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 -1 0 -1 -1 90 pathToFile.root histName omega3pi +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 -1 0 -1 -1 polAngle=90 polFile=pathToFile.root polHist=histName omega3pi +// For a 2P m=0 wave from an omega->gpi that belongs to the real part of +// the amplitude sum,has a negative sign in front of the polarization term, and +// a photon helicity -1, you would write: +// you would write: +// If polarization info stored in the beam P4: +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 omegagpi -1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 omegagpi gHelicity=-1 +// If polarization angle and fraction are fixed: +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 90 .35 omegagpi -1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 polFraction=0.35 polAngle=90 omegagpi gHelicity=-1 +// If polarization vs E_gamma from histogram: +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 90 pathToFile.root histName omegagpi -1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 polAngle=90 polFile=pathToFile.root polHist=histName omegagpi gHelicity=-1 +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// When the reaction is specified in the config file, the following format is +// assumed: +// reaction reactionName Beam Proton Ps VecDaught1 VecDaught2 [VecDaught3] +// +// This order should be reflected in the way the particles are saved in the +// input file. Mistakes deviating from this assumption will not trigger +// an error and will lead to incorrect results. +// Assumed ordering of particles in the array from the input file +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #include #include #include +#include #include "TLorentzVector.h" #include "TLorentzRotation.h" @@ -19,8 +130,38 @@ #include "UTILITIES/BeamProperties.h" -// Function to check if a string is a valid number -static bool isValidNumber(const string& argInput, double &value){ +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........ Structs ........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +struct Vec_ps_refl_Args { + // Required positional args + int j; // Total spin + int m; // Spin projection + int l; // Relative orbital angular momentum + int r; // +1/-1 for (r)eal/imaginary part of the amplitude + int s; // +1/-1 (s)ign of the P_gamma term + + // Polarization information + bool polInfoInPhotonP4 = true; + double polAngle = -2.0; // -2 = not set + double polFraction = -2.0; // -2 = not set + // Amorphous sometimes is referred as having a polAngle -1, so we use -2 to + // avoid possible confusion for both polAngle and polFraction + std::string polFile = ""; // path to .root file with polarization vs E_beam + std::string polHist = ""; // name of TH1D inside polFile + + // Default is 2-body vector decay, but if omega is used as vector: + // Omega decay mode + bool omega3pi = false; + bool omegagpi0 = false; + int gHelicity = 0; // +1/-1 helicity of the bachelor photon in omega->g pi0 + + bool noBarrier = false; // Turn off barrier factors +}; +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........ Helper Functions ........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +static bool isValidNumber(const std::string& argInput, double &value){ char* end = nullptr; errno = 0; // reset global error value = std::strtod(argInput.c_str(), &end); @@ -32,124 +173,239 @@ static bool isValidNumber(const string& argInput, double &value){ if(end == argInput.c_str() || *end != '\0' || errno != 0) { return false; // not a valid number } - // If end points to the end of string, it's fully numeric + // If "end" points to the end of string, it's fully numeric return true; } -static double parseValidatedNumber(const string& label, const string& argInput){ +static double parseValidatedNumber(const std::string& label, const std::string& argInput){ double tmpValue = 0.0; if(!isValidNumber(argInput, tmpValue)){ - throw std::invalid_argument("Vec_ps_refl: invalid " + label + ": " + argInput); + throw std::invalid_argument("[ Vec_ps_refl ]: invalid " + label + ": " + argInput); } return tmpValue; } +static Vec_ps_refl_Args parsedArgs( const std::vector& args ){ + Vec_ps_refl_Args inputArgs; + + // Required positional args + inputArgs.j = static_cast( parseValidatedNumber( "J", args[0] ) ); + inputArgs.m = static_cast( parseValidatedNumber( "M", args[1] ) ); + inputArgs.l = static_cast( parseValidatedNumber( "L", args[2] ) ); + inputArgs.r = static_cast( parseValidatedNumber( "Re/Im", args[3] ) ); + inputArgs.s = static_cast( parseValidatedNumber( "P_gamma sign", args[4] ) ); + + // Physics validation of required args + if( abs(inputArgs.m) > inputArgs.j ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: |M| must be <= J, got J=" + args[0] + " M=" + args[1] ); + if( abs(inputArgs.r) != 1 ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: Re/Im flag must be +1 or -1, got " + args[3] ); + if( abs(inputArgs.s) != 1 ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: P_gamma sign must be +1 or -1, got " + args[4] ); + + // If no optional args, return now + if( args.size() <= 5 ) return inputArgs; + + // Detect format by checking if args[5] contains '=' + bool isKeyValue = std::any_of( args.begin() + 5, args.end(), + []( const std::string& s ){ return s.find('=') != std::string::npos; } ); + + // Lambda function to detect bare flags (valid in both formats) + auto isBareFlag = []( const std::string& s ){ + return s == "omega3pi" || s == "omegagpi0" || s == "nobarrier"; + }; + + if( !isKeyValue ){ + size_t i = 5; + + // arg[5]: polarization angle + if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", args[i++] ); + inputArgs.polInfoInPhotonP4 = false; + } -Vec_ps_refl::Vec_ps_refl( const vector< string >& args ) : -UserAmplitude< Vec_ps_refl >( args ){ + // arg[6]: fixed fraction or .root file path + if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ + if( args[i].find(".root") != std::string::npos ){ + inputArgs.polFile = args[i++]; + if( i >= args.size() || isBareFlag(args[i]) ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: .root polarization file requires a histogram name immediately after it" ); + inputArgs.polHist = args[i++]; + } else { + inputArgs.polFraction = parseValidatedNumber( "polFraction", args[i++] ); + } + } - // 5 possibilities to initialize this amplitude: - // : total spin, - // : spin projection, - // : partial wave, - // : +1/-1 for real/imaginary part; - // : +1/-1 sign in P_gamma term - - // Resonance spin J - m_j = static_cast(parseValidatedNumber("J", args[0])); - // Spin projection (Lambda) - m_m = static_cast(parseValidatedNumber("M", args[1])); - // Partial wave L - m_l = static_cast(parseValidatedNumber("L", args[2])); - // Real (+1) or imaginary (-1) - m_r = static_cast(parseValidatedNumber("Re/Im", args[3])); - // Sign for polarization in amplitude - m_s = static_cast(parseValidatedNumber("P_gamma sign", args[4])); - - // make sure values are reasonable - assert( abs( m_m ) <= m_j ); - // m_r = +1 for real - // m_r = -1 for imag - assert( abs( m_r ) == 1 ); - // m_s = +1 for 1 + Pgamma - // m_s = -1 for 1 - Pgamma - assert( abs( m_s ) == 1 ); - - // Default polarization information stored in tree - m_polInTree = true; - - // Default is 2-body vector decay (set flag in config file for omega->3pi) - m_3pi = false; - - // Loop over any additional amplitude arguments to change defaults - for(uint ioption=5; ioptionGet(args[7].c_str()); - if(polFrac_vs_E == nullptr ){ - throw std::runtime_error( - "Vec_ps_refl ERROR: Could not find histogram '" + args[7] + - "' in file " + std::string(polOption.Data())); - } + // Remaining: bare flags, and gHelicity bare integer after omegagpi0 + while( i < args.size() ){ + if( args[i] == "omega3pi" ){ + inputArgs.omega3pi = true; + i++; + } + else if( args[i] == "omegagpi0" ){ + inputArgs.omegagpi0 = true; + if( i + 1 >= args.size() ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: omegagpi0 requires a photon helicity value (+1 or -1) immediately after it" ); + inputArgs.gHelicity = static_cast( parseValidatedNumber( "gHelicity", args[i+1] ) ); + i += 2; + } + else if( args[i] == "nobarrier" ){ + inputArgs.noBarrier = true; + i++; } else{ - polFraction = parseValidatedNumber("polarization fraction", args[6]); - + throw std::invalid_argument( + "[ Vec_ps_refl ]: unrecognized positional argument '" + args[i] + "'" ); } } - // Check for omega->3pi option - if(option.EqualTo("omega3pi")) m_3pi = true; - } + } + else { // key=value format + for( size_t i = 5; i < args.size(); i++ ){ + const std::string& arg = args[i]; + + if( arg == "omega3pi" ){ inputArgs.omega3pi = true; continue; } + if( arg == "omegagpi0" ){ inputArgs.omegagpi0 = true; continue; } + if( arg == "nobarrier" ){ inputArgs.noBarrier = true; continue; } + auto sep = arg.find('='); // '=' is used as separator + if( sep == std::string::npos ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: unrecognized argument '" + arg + + "' (expected key=value pair or known flag)" ); + + const std::string key = arg.substr(0, sep); + const std::string value = arg.substr(sep + 1); + + if( key == "polAngle" ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", value ); + inputArgs.polInfoInPhotonP4 = false; + } + else if( key == "polFraction" ){ + inputArgs.polFraction = parseValidatedNumber( "polFraction", value ); + } + else if( key == "polFile" ){ inputArgs.polFile = value; } + else if( key == "polHist" ){ inputArgs.polHist = value; } + else if( key == "gHelicity" ){ + inputArgs.gHelicity = static_cast( parseValidatedNumber( "gHelicity", value ) ); + } + else{ + throw std::invalid_argument( + "[ Vec_ps_refl ]: unrecognized key '" + key + "'" ); + } + } + } + + // Post-parse validation + if( inputArgs.polFile != "" && inputArgs.polHist == "" ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: polFile='" + inputArgs.polFile + + "' requires polHist= to also be specified" ); + + if( inputArgs.polHist != "" && inputArgs.polFile == "" ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: polHist='" + inputArgs.polHist + + "' was specified without polFile=" ); + + if( inputArgs.polFraction >= 0.0 && inputArgs.polFile != "" ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: polFraction and polFile are mutually exclusive; specify only one" ); + + if( inputArgs.omegagpi0 && abs(inputArgs.gHelicity) != 1 ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: omegagpi0 requires gHelicity=+1 or gHelicity=-1" ); + + if( inputArgs.gHelicity != 0 && !inputArgs.omegagpi0 ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: gHelicity was specified but omegagpi0 flag is missing" ); + + if( inputArgs.omega3pi && inputArgs.omegagpi0 ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: omega3pi and omegagpi0 are mutually exclusive; specify only one" ); + + if( !inputArgs.polInfoInPhotonP4 && inputArgs.polFraction < 0.0 && inputArgs.polFile == "" ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: polAngle was given but neither polFraction nor polFile+polHist were provided" ); + + return inputArgs; } +Vec_ps_refl::Vec_ps_refl( const vector< string >& args ) : +UserAmplitude< Vec_ps_refl >( args ){ + Vec_ps_refl_Args inputArgs = parsedArgs( args ); + m_j = inputArgs.j; + m_m = inputArgs.m; + m_l = inputArgs.l; + m_r = inputArgs.r; + m_s = inputArgs.s; + + m_polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; + polAngle = inputArgs.polAngle; + polFraction = inputArgs.polFraction; + polFrac_vs_E = nullptr; + + m_3pi = inputArgs.omega3pi; + m_gpi0 = inputArgs.omegagpi0; + m_ghel = inputArgs.gHelicity; + m_noBarrier = inputArgs.noBarrier; + + + // Open polarization file and retrieve histogram + if( inputArgs.polFile != "" ){ + TFile* f = new TFile( inputArgs.polFile.c_str() ); + if( !f || f->IsZombie() ) + throw std::runtime_error( + "[ Vec_ps_refl ]: could not open polarization file '" + inputArgs.polFile + "'" ); + polFrac_vs_E = (TH1D*)f->Get( inputArgs.polHist.c_str() ); + if( !polFrac_vs_E ) + throw std::runtime_error( + "[ Vec_ps_refl ]: histogram '" + inputArgs.polHist + + "' not found in '" + inputArgs.polFile + "'" ); + // Detach histogram from file so it persists after file is closed + polFrac_vs_E->SetDirectory(0); + f->Close(); + delete f; + } +} + +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........ Main Functions ........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + void Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ TLorentzVector beam; - TVector3 eps; + TVector3 eps; // beam polarization vector (eps)ilon GDouble beam_polFraction; GDouble beam_polAngle; - if(m_polInTree){ - beam.SetPxPyPzE( 0.0, 0.0, pKin[0][0], pKin[0][0]); - eps.SetXYZ(pKin[0][1], pKin[0][2], 0.0); // beam polarization vector; - + if(m_polInfoInPhotonP4){ + // When pol info is stored in the photon 4-vector in the tree + // the energy (pKin[0][0]) and the pz (pKin[0][3]) are used as normal + beam.SetPxPyPzE( 0.0, 0.0, pKin[0][3], pKin[0][0]); + // while the px and py components store the pol info + // The values should be stored as px = polFraction*cos(polAngle) + // and py = polFraction*sin(polAngle) + eps.SetXYZ(pKin[0][1], pKin[0][2], 0.0); beam_polFraction = eps.Mag(); beam_polAngle = eps.Phi(); } else{ beam.SetPxPyPzE( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); beam_polAngle = polAngle; - - if(polFraction > 0.0){ // for fitting with fixed polarization + // for fixed polarization fraction + if(polFrac_vs_E == nullptr) beam_polFraction = polFraction; - } else{ // for fitting with polarization vs E_gamma from input histogram int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){ - beam_polFraction = 0.0; - // Warn if using events outside histogram range - static bool warned_hist = false; // only warn once - if(!warned_hist){ - report( NOTICE ) - << " Energy for polFraction: " << pKin[0][0] - << " appears to be outside histogram range." - << " the value will be set to zero." - << " This message will only be shown once." - << endl; - warned_hist = true; - } + throw std::runtime_error( + "[ Vec_ps_refl ]: energy " + std::to_string(pKin[0][0]) + + " outside histogram range" ); } else beam_polFraction = polFrac_vs_E->GetBinContent(bin); } @@ -164,10 +420,11 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ // Make sure the order of daughters is correct in the config file! TLorentzVector vec, vec_daught1, vec_daught2; - if(m_3pi){ - // Omega ps proton, omega -> 3pi (6 particles) - // Omega pi- Delta++, omega -> 3pi (7 particles) + // Omega ps proton, omega -> 3pi (6 particles): + // beam proton ps pi0 pip pim + // Omega pi- Delta++, omega -> 3pi (6 particles): + // beam delta ps pi0 pip pim TLorentzVector pi0(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); TLorentzVector pip(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); TLorentzVector pim(pKin[5][1], pKin[5][2], pKin[5][3], pKin[5][0]); @@ -176,13 +433,17 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ vec_daught2 = pim; } else{ - // Omega ps proton, omega -> pi0 g (4 particles) - // Omega pi- Delta++, omega -> pi0 g (5 particles) - - // (vec 2-body) ps proton, vec 2-body -> pipi, KK (5 particles) - // (vec 2-body) pi- Delta++, vec 2-body -> pipi, KK (6 particles) - // (vec 2-body) K+ Lambda, vec 2-body -> Kpi (6 particles) - vec_daught1 = TLorentzVector(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); + // Omega ps proton, omega -> gpi0 (5 particles): + // beam proton ps pi0 gamma + // Omega pi- Delta++, omega -> gpi (5 particles): + // beam delta ps pi0 gamma + // Vec(KK) ps proton, (5 particles) + // beam proton ps K K + // Vec(pipi) pi- Delta++, (5 particles) + // beam proton ps pi pi + // Vec(Kpi) K+ Lambda, (5 particles) + // beam proton ps K pi + vec_daught1 = TLorentzVector(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); vec_daught2 = TLorentzVector(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); vec = vec_daught1 + vec_daught2; } @@ -192,23 +453,22 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ ///////////////// Boost Particles and Get Angles///////////////////// - TLorentzVector target(0,0,0,0.938); - // Helicity coordinate system - TLorentzVector beamP = beam + target; + TLorentzVector target(0,0,0,0.938); // proton at rest + TLorentzVector beamTarget = beam + target; // Calculate decay angles for X in helicity frame (same for all vectors) // Change getXDecayAngles to get Gottfried-Jackson angles if needed // Note: it also calculates the production angle - vector xDecayAngles = getXDecayAngles( beam_polAngle, beam, beamP, X, vec); + vector xDecayAngles = getXDecayAngles( beam_polAngle, beam, beamTarget, X, vec); // Calculate vector decay angles (unique for each vector) vector vectorDecayAngles; if(m_3pi){ - vectorDecayAngles = getVectorDecayAngles( beamP, X, vec, + vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, vec_daught1, vec_daught2); } else{ - vectorDecayAngles = getVectorDecayAngles( beamP, X, vec, + vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, vec_daught1, TLorentzVector(0,0,0,0)); } @@ -222,17 +482,39 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ GDouble m_ps = ps.M(); complex amplitude(0,0); - complex i(0,1); - - for (int lambda = -1; lambda <= 1; lambda++) { // sum over vector helicity - GDouble hel_amp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); - amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * - hel_amp * conj(wignerD(1, lambda, 0, cosThetaH, phiH)); + static const complex i(0,1); + + if(m_gpi0){ // radiative omega decay requires a handling of the photon helicity + for (int lambda = -1; lambda <= 1; lambda++) { // sum over vector helicity + GDouble hel_amp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); + if(lambda==0){ + amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * + hel_amp * conj(wignerD(1, lambda, m_ghel, cosThetaH, phiH)) * + (m_ghel*1.0); // m_ghel is int + } + else{ + amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * + hel_amp * conj(wignerD(1, -1*lambda, m_ghel, cosThetaH, phiH)) * + -1.0; + // the power of lambda is irrelevant for lambda =/= 0: (-1)^lambda = -1 + } + } + } + else{ // for any other vector decay + for (int lambda = -1; lambda <= 1; lambda++) { // sum over vector helicity + GDouble hel_amp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); + amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * + hel_amp * conj(wignerD(1, lambda, 0, cosThetaH, phiH)); + } } + // The amplitude is multiplied by a factor, either sqrt(1 + -P_gamma) or + // sqrt(1 + P_gamma) depending on the what sum is being calculated GDouble factor = sqrt(1 + m_s * beam_polFraction); + // The result of the function that depends on the angles is stored in zjm complex zjm = 0; - // - -> + in prod_angle + // A - sign translates to + in the prod_angle + // This is because reflectivity conventions introduce exp(-iPhi) complex rotateY = polar((GDouble)1., (GDouble)(-1. * prod_angle )); if (m_r == 1) @@ -240,18 +522,20 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ if (m_r == -1) zjm = i*imag(amplitude * rotateY); - // E852 Nozar thesis has sqrt(2*s+1)*sqrt(2*l+1)*F_l(p_omega)*sqrt(omega) - double kinFactor = barrierFactor(m_X, m_l, m_vec, m_ps); - //kinFactor *= sqrt(3.) * sqrt(2.*m_l + 1.); + if( !m_noBarrier ){ + GDouble kinFactor = barrierFactor(m_X, m_l, m_vec, m_ps); factor *= kinFactor; - + // Alternatives: + // E852 Nozar thesis has sqrt(2*s+1)*sqrt(2*l+1)*F_l(p_omega)*sqrt(omega) + // kinFactor *= sqrt(3.) * sqrt(2.*m_l + 1.); + } + userVars[uv_ampRe] = ( factor * zjm ).real(); userVars[uv_ampIm] = ( factor * zjm ).imag(); return; } - /////////////////////// Amplitude Calculation ////////////////////////// complex< GDouble > diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h index b583113a9..052f184c6 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h @@ -12,7 +12,6 @@ #include using std::complex; -using namespace std; #ifdef GPU_ACCELERATION void @@ -81,13 +80,18 @@ class Vec_ps_refl : public UserAmplitude< Vec_ps_refl > int m_l; int m_r; int m_s; - int m_3pi; - //AmpParameter polAngle; + bool m_polInfoInPhotonP4; GDouble polFraction; GDouble polAngle; - bool m_polInTree; TH1D *polFrac_vs_E; + + bool m_3pi; + bool m_gpi0; + int m_ghel; + bool m_noBarrier; + + }; #endif diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc index 15ab13dd5..35fd3d499 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc @@ -1,13 +1,33 @@ +#include + #include "TLorentzVector.h" #include "TLorentzRotation.h" #include "AMPTOOLS_AMPS/vecPsAngles.h" - #include "AMPTOOLS_DATAIO/VecPsPlotGenerator.h" + #include "IUAmpTools/Histogram1D.h" #include "IUAmpTools/Kinematics.h" -// Function to check if a string is a valid number +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........ Structs ........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +struct Vec_ps_plot_Args { + bool isMoment = false; // default Vec_ps_refl, true if Vec_ps_moment detected + + // Polarization information + bool polInfoInPhotonP4 = true; + double polAngle = -2.0; // -2 = not set + // Amorphous sometimes is referred as having a polAngle -1, so we use -2 to + // avoid possible confusion for both polAngle + + // Default is 2-body vector decay, but if omega is used as vector: + // Omega decay mode + bool omega3pi = false; +}; +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +//....oooOO0OOooo........ Helper Functions ........oooOO0OOooo..... +//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ static bool isValidNumber(const string& argInput, double &value){ char* end = nullptr; errno = 0; // reset global error @@ -27,11 +47,73 @@ static bool isValidNumber(const string& argInput, double &value){ static double parseValidatedNumber(const string& label, const string& argInput){ double tmpValue = 0.0; if(!isValidNumber(argInput, tmpValue)){ - throw std::invalid_argument("Vec_ps_refl: invalid " + label + ": " + argInput); + throw std::invalid_argument("[ VecPsPlotGenerator ]: invalid " + label + ": " + argInput); } return tmpValue; } +static Vec_ps_plot_Args parsedArgs( const std::vector& args ){ + Vec_ps_plot_Args inputArgs; + + inputArgs.isMoment = ( args[0] == "Vec_ps_moment" ); + if( inputArgs.isMoment ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", args[1] ); + inputArgs.omega3pi = true; + inputArgs.polInfoInPhotonP4 = false; + return inputArgs; + } + + // If no optional args, return now + if( args.size() <= 6 ) return inputArgs; + + // Detect format by checking if args[6] contains '=' + bool isKeyValue = std::any_of( args.begin() + 6, args.end(), + []( const std::string& s ){ return s.find('=') != std::string::npos; } ); + + // Lambda function to detect bare flags (valid in both formats) + auto isBareFlag = []( const std::string& s ){ + return s == "omega3pi" || s == "omegagpi0" || s == "nobarrier"; + }; + + if( !isKeyValue ){ + size_t i = 6; + + // arg[6]: polarization angle + if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", args[i++] ); + inputArgs.polInfoInPhotonP4 = false; + } + + // Remaining: bare flags, and gHelicity bare integer after omegagpi0 + for( size_t i = 6; i < args.size(); i++ ){ + if( args[i] == "omega3pi" ){ + inputArgs.omega3pi = true; + } + } + } + else { // key=value format + for( size_t i = 6; i < args.size(); i++ ){ + const std::string& arg = args[i]; + + if( arg == "omega3pi" ){ inputArgs.omega3pi = true; continue; } + auto sep = arg.find('='); // '=' is used as separator + if( sep == std::string::npos ) + throw std::invalid_argument( + "[ Vec_ps_refl ]: unrecognized argument '" + arg + + "' (expected key=value pair or known flag)" ); + + const std::string key = arg.substr(0, sep); + const std::string value = arg.substr(sep + 1); + + if( key == "polAngle" ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", value ); + inputArgs.polInfoInPhotonP4 = false; + } + } + } + return inputArgs; +} + /* Constructor to display FitResults */ VecPsPlotGenerator::VecPsPlotGenerator( const FitResults& results, Option opt ) : PlotGenerator( results, opt ) @@ -84,102 +166,108 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin ){ // This function will make this class backwards-compatible with older versions // (v0.10.x and prior) of AmpTools, but will not be able to properly obtain - // the polariation plane in the lab when multiple orientations are used + // the polarization plane in the lab when multiple orientations are used projectEvent( kin, "" ); } void VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ - //cout << "project event" << endl; - TLorentzVector beam = kin->particle( 0 ); - TLorentzVector recoil = kin->particle( 1 ); - // 1st after proton is always the pseudoscalar (bachelor) meson - TLorentzVector bach = kin->particle( 2 ); - - // Default is 2-body vector decay (set flag in config file for omega->3pi) - bool m_3pi = false; - - // Min particle index for recoil sum - int min_recoil = 5; - - - // Properly read polarization angle from config file if provided - double beam_polAngle=0; - // Check config file for optional parameters -- we assume here that the - // first amplitude in the list is a Vec_ps_refl amplitude or a Vec_ps_moment - const vector args = - cfgInfo()->amplitudeList(reactionName, "", "").at(0)->factors().at(0); - for(uint ioption = 0; ioption < args.size(); ioption++){ - TString option = args[ioption].c_str(); - if(option.EqualTo("Vec_ps_moment")){ - // Force the 3pi decay as it's currently the only one handled by it - m_3pi = true; - break; - } - else{ - if(ioption == 6){ - beam_polAngle = parseValidatedNumber("polarization angle", args[6]); - } - if(option.EqualTo("omega3pi")){ - m_3pi = true; - } - } - } + // The first argument in the list is Vec_ps_refl or Vec_ps_moment + const vector args = + cfgInfo()->amplitudeList(reactionName, "", "").at(0)->factors().at(0); + Vec_ps_plot_Args inputArgs = parsedArgs( args ); + + bool m_polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; + double beam_polAngle = inputArgs.polAngle; + + bool m_3pi = inputArgs.omega3pi; + + //cout << "project event" << endl; + TLorentzVector readBeam = kin->particle( 0 ); + TLorentzVector recoil = kin->particle( 1 ); + // 1st after proton is always the pseudoscalar (bachelor) meson + TLorentzVector ps = kin->particle( 2 ); + + TVector3 eps; // beam polarization vector (eps)ilon + TLorentzVector beam; + if(m_polInfoInPhotonP4){ + // When pol info is stored in the photon 4-vector in the tree + // the energy (pKin[0][0]) and the pz (pKin[0][3]) are used as normal + beam.SetPxPyPzE( 0.0, 0.0, readBeam.Pz(), readBeam.E() ); + // while the px and py components store the pol info + // The values should be stored as px = polFraction*cos(polAngle) + // and py = polFraction*sin(polAngle) + eps.SetXYZ(readBeam.Px(), readBeam.Py(), 0.0); + beam_polAngle = eps.Phi(); + } + else{ + beam = readBeam; + } - // Compute vector meson from its decay products - // Make sure the order of daughters is correct in the config file! - TLorentzVector vec, vec_daught1, vec_daught2; - double dalitz_s, dalitz_t, dalitz_u, dalitz_d, dalitz_sc; - double dalitz_x = 0.0; - double dalitz_y = 0.0; - - if(m_3pi) { - // Omega ps proton, omega -> 3pi (6 particles) - // Omega pi- Delta++, omega -> 3pi (7 particles) - TLorentzVector pi0 = kin->particle( 3 ); - TLorentzVector pip = kin->particle( 4 ); - TLorentzVector pim = kin->particle( 5 ); - vec = pi0 + pip + pim; - vec_daught1 = pip; - vec_daught2 = pim; - min_recoil = 6; - - // Dalitz variables + // Min particle index for recoil sum + int min_recoil = 5; + + // Compute vector meson from its decay products + // Make sure the order of daughters is correct in the config file! + TLorentzVector vec, vec_daught1, vec_daught2; + double dalitz_s, dalitz_t, dalitz_u, dalitz_d, dalitz_sc; + double dalitz_x = 0.0; + double dalitz_y = 0.0; + + if(m_3pi) { + // Omega ps proton, omega -> 3pi (6 particles): + // beam proton ps pi0 pip pim + // Omega pi- Delta++, omega -> 3pi (6 particles): + // beam delta ps pi0 pip pim + TLorentzVector pi0 = kin->particle( 3 ); + TLorentzVector pip = kin->particle( 4 ); + TLorentzVector pim = kin->particle( 5 ); + vec = pi0 + pip + pim; + vec_daught1 = pip; + vec_daught2 = pim; + min_recoil = 6; + + // Dalitz variables are included here. + // They amplitude tends to be multiply in the config file but that's not + // somthing we know here (can we?) const auto& p2 = pi0; const auto& p3 = pip; const auto& p4 = pim; - double pdg_m_pi0 = 0.1349766; - double pdg_m_chargedPi = 0.13957018; - dalitz_s = (p3 + p4).M2(); // s=M(pip pim) - dalitz_t = (p2 + p3).M2(); // s=M(pip pi0) - dalitz_u = (p2 + p4).M2(); // s=M(pim pi0) - dalitz_d = 2 * (p2 + p3 + p4).M() * - ((p2 + p3 + p4).M() - ((2 * pdg_m_chargedPi) + pdg_m_pi0)); - dalitz_sc = (1.0 / 3.0) * ((p2 + p3 + p4).M2() + - ((2 * (pdg_m_chargedPi * pdg_m_chargedPi)) + - (pdg_m_pi0 * pdg_m_pi0))); - dalitz_x = sqrt(3.0) * (dalitz_t - dalitz_u) / dalitz_d; - dalitz_y = 3.0 * (dalitz_sc - dalitz_s) / dalitz_d; - - } + double pdg_m_pi0 = 0.1349766; + double pdg_m_chargedPi = 0.13957018; + dalitz_s = (p3 + p4).M2(); // s=M(pip pim) + dalitz_t = (p2 + p3).M2(); // s=M(pip pi0) + dalitz_u = (p2 + p4).M2(); // s=M(pim pi0) + dalitz_d = 2 * (p2 + p3 + p4).M() * + ((p2 + p3 + p4).M() - ((2 * pdg_m_chargedPi) + pdg_m_pi0)); + dalitz_sc = (1.0 / 3.0) * ((p2 + p3 + p4).M2() + + ((2 * (pdg_m_chargedPi * pdg_m_chargedPi)) + + (pdg_m_pi0 * pdg_m_pi0))); + dalitz_x = sqrt(3.0) * (dalitz_t - dalitz_u) / dalitz_d; + dalitz_y = 3.0 * (dalitz_sc - dalitz_s) / dalitz_d; + } else { - // omega ps proton, omega -> pi0 g (4 particles) - // omega pi- Delta++, omega -> pi0 g (5 particles) - - // (vec 2-body) ps proton, vec 2-body -> pipi, KK (5 particles) - // (vec 2-body) pi- Delta++, vec 2-body -> pipi, KK (6 particles) - // (vec 2-body) K+ Lambda, vec 2-body -> Kpi (6 particles) + // Omega ps proton, omega -> gpi0 (5 particles): + // beam proton ps pi0 gamma + // Omega pi- Delta++, omega -> gpi (5 particles): + // beam delta ps pi0 gamma + // Vec(KK) ps proton, (5 particles) + // beam proton ps K K + // Vec(pipi) pi- Delta++, (5 particles) + // beam proton ps pi pi + // Vec(Kpi) K+ Lambda, (5 particles) + // beam proton ps K pi vec_daught1 = kin->particle( 3 ); - vec_daught2 = kin->particle( 4 ); - vec = vec_daught1 + vec_daught2; + vec_daught2 = kin->particle( 4 ); + vec = vec_daught1 + vec_daught2; min_recoil = 5; } // Final meson system P4 - TLorentzVector X = vec + bach; + TLorentzVector X = vec + ps; - TLorentzVector proton_ps = recoil + bach; + TLorentzVector proton_ps = recoil + ps; TLorentzVector recoil_ps = proton_ps; for(uint i=min_recoil; iparticleList().size(); i++){ // Add mesons to recoil system (e.g. Delta or Lambda) @@ -187,30 +275,29 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ recoil_ps += kin->particle(i); } - TLorentzVector target(0,0,0,0.938); - // Helicity coordinate system - TLorentzVector beamP = beam + target; + TLorentzVector target(0,0,0,0.938); // proton at rest + TLorentzVector beamTarget = beam + target; // Calculate decay angles for X in helicity frame (same for all vectors) // Change getXDecayAngles to get Gottfried-Jackson angles if needed // Note: it also calculates the production angle - vector xDecayAngles = getXDecayAngles( beam_polAngle, beam, beamP, X, vec); + vector xDecayAngles = getXDecayAngles( beam_polAngle, beam, beamTarget, X, vec); // set polarization angle to zero to see shift in Phi_Prod distributions - vector xDecayAngles_offset = getXDecayAngles( 0, beam, beamP, X, vec); + vector xDecayAngles_offset = getXDecayAngles( 0, beam, beamTarget, X, vec); // Calculate vector decay angles (unique for each vector) vector vectorDecayAngles; if(m_3pi){ - vectorDecayAngles = getVectorDecayAngles( beamP, X, vec, + vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, vec_daught1, vec_daught2); } else{ - vectorDecayAngles = getVectorDecayAngles( beamP, X, vec, + vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, vec_daught1, TLorentzVector(0,0,0,0)); } double Mandt = fabs((target-recoil).M2()); - double recoil_mass = recoil.M(); + double recoil_mass = recoil.M2(); GDouble cosTheta = TMath::Cos(xDecayAngles[0]); GDouble phi = xDecayAngles[1]; diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h index 4f54c2a5b..cd45365c6 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h @@ -1,4 +1,4 @@ -#if !(defined VCPSPLOTGENERATOR) +#if !(defined VECPSPLOTGENERATOR) #define VECPSPLOTGENERATOR #include @@ -34,10 +34,13 @@ class VecPsPlotGenerator : public PlotGenerator case VecPsPlotGenerator::kCosThetaH: return "CosTheta_H"; case VecPsPlotGenerator::kPhiH: return "Phi_H"; case VecPsPlotGenerator::kProd_Ang: return "Prod_Ang"; + case VecPsPlotGenerator::kProdOffset: return "ProdOffset"; case VecPsPlotGenerator::kt: return "t"; + case VecPsPlotGenerator::kRecoilMass: return "MRecoil"; + case VecPsPlotGenerator::kProtonPsMass: return "MProtonPs"; case VecPsPlotGenerator::kRecoilPsMass: return "MRecoilPs"; case VecPsPlotGenerator::kLambda: return "Lambda"; - case VecPsPlotGenerator::kProdOffset: return "ProdOffset"; + case VecPsPlotGenerator::kDalitz: return "Dalitz"; case VecPsPlotGenerator::kPhi_ProdVsPhi: return "Phi_ProdVsPhi"; case VecPsPlotGenerator::kPhiOffsetVsPhi: return "PhiOffsetVsPhi"; // Add more variables here if needed diff --git a/src/programs/AmplitudeAnalysis/fit/fit.cc b/src/programs/AmplitudeAnalysis/fit/fit.cc index b723de23d..6adfc70a0 100644 --- a/src/programs/AmplitudeAnalysis/fit/fit.cc +++ b/src/programs/AmplitudeAnalysis/fit/fit.cc @@ -71,23 +71,146 @@ using std::complex; using namespace std; -void summarizeFits(std::vector>& fitLLs) { - if(fitLLs.size() == 0) return; - std::sort(fitLLs.begin(), fitLLs.end(), [](const std::tuple& a, const std::tuple& b){ return std::get<4>(a) < std::get<4>(b); }); - std::ofstream fout("fit_ranking.txt"); - if(!fout) std::cerr << "Error: cannot open fit_ranking.txt\n"; - auto print_header = [](std::ostream& os){ os << "\nSUMMARY OF ALL FITS:\n" << std::left << std::setw(6) << "#" << std::setw(9) << "Success" << std::setw(11) << "FitStatus" << std::setw(9) << "eMatrix" << std::setw(14) << "LogL" << '\n'; }; - auto print_row = [&](std::ostream& os, size_t i){ os << std::left << std::setw(6) << std::get<0>(fitLLs[i]) << std::setw(9) << (std::get<1>(fitLLs[i]) ? "Y" : "N") << std::setw(11) << std::get<2>(fitLLs[i]) << std::setw(9) << std::get<3>(fitLLs[i]) << std::setw(14) << (std::ostringstream() << std::setprecision(8) << std::defaultfloat << std::get<4>(fitLLs[i])).str() << '\n'; }; - print_header(std::cout); - if(fout) print_header(fout); - for(size_t i=0; i>& fitLLs, bool csvOutput = false, + double tieEps = 1e-3) { + // tieEps: LogL difference below which two converged fits are considered "tied". + // MINUIT's default EDM convergence criterion is on the order of 1e-3 (for UP=1, + // i.e. -logL directly), so two independent runs that land on the same true + // minimum will typically agree to within about this tolerance. Tighten this + // (e.g. 1e-4) if you called MIGRAD with a smaller "tolerance" arg / SetErrorDef, + // or loosen it if you want a "statistically indistinguishable" cut instead + // (e.g. ~2.0, akin to a likelihood-ratio / AIC-difference threshold). + + if(fitLLs.empty()) return; + + // sort: converged (success) fits first, ranked by LogL ascending; + // failed fits pushed to the bottom (also LogL-sorted among themselves, but unranked) + std::sort(fitLLs.begin(), fitLLs.end(), + [](const auto& a, const auto& b){ + bool sa = std::get<1>(a), sb = std::get<1>(b); + if (sa != sb) return sa; // successes before failures + return std::get<4>(a) < std::get<4>(b); + }); + + // collect all *successful* fits sharing the minimum LogL (ties, within tieEps) + auto tie_min = [&](int idx) { + using FitType = std::tuple; + double minVal = std::numeric_limits::max(); + for (const auto& v : fitLLs) { + if (!std::get<1>(v)) continue; // only consider converged fits + double val = (idx == 4) ? std::get<4>(v) : (double)std::get<3>(v); + if (val < minVal) minVal = val; + } + std::vector out; + for (const auto& v : fitLLs) { + if (!std::get<1>(v)) continue; + double val = (idx == 4) ? std::get<4>(v) : (double)std::get<3>(v); + if (std::fabs(val - minVal) < tieEps) out.push_back(v); + } + return out; + }; + + // success rate + size_t succ = 0; + for (const auto& f : fitLLs) if (std::get<1>(f)) ++succ; + double pct = 100.0 * static_cast(succ) / static_cast(fitLLs.size()); + + // best successful fit (lowest LogL among successes only) + int bestSuccIdx = -1; + double bestSuccLogL = std::numeric_limits::max(); + for (const auto& f : fitLLs) { + if (std::get<1>(f) && std::get<4>(f) < bestSuccLogL) { + bestSuccLogL = std::get<4>(f); + bestSuccIdx = std::get<0>(f); + } } - size_t succ = 0; // print overall success rate - for (size_t i = 0; i < fitLLs.size(); ++i) if (std::get<1>(fitLLs[i])) ++succ; - double pct = fitLLs.empty() ? 0.0 : (100.0 * static_cast(succ) / static_cast(fitLLs.size())); - std::cout << "\nSuccess: " << succ << " / " << fitLLs.size() << " (" << std::setprecision(3) << std::fixed << pct << "%)\n"; + + auto bestByLogL = tie_min(4); + + // reference LogL for DeltaLogL column: best converged fit if one exists, + // otherwise just fall back to the first row so the column still prints + double refLogL = (bestSuccIdx >= 0) ? bestSuccLogL : std::get<4>(fitLLs[0]); + + // ---- write to both stdout and file ---- + std::string outfile = csvOutput ? "fit_ranking.csv" : "fit_ranking.txt"; + std::ofstream fout(outfile); + if (!fout) std::cerr << "Error: cannot open " << outfile << "\n"; + + auto write = [&](std::ostream& os) { + + if (csvOutput) { + os << "Rank,FitIndex,Success,FitStatus,eMatrix,LogL,DeltaLogL\n"; + int rank = 0; + for (size_t i = 0; i < fitLLs.size(); ++i) { + const auto& f = fitLLs[i]; + bool isSucc = std::get<1>(f); + if (isSucc) ++rank; + os << (isSucc ? std::to_string(rank) : std::string("-")) << "," + << std::get<0>(f) << "," + << (isSucc ? 1:0) << "," + << std::get<2>(f) << "," + << std::get<3>(f) << "," + << std::setprecision(10) << std::fixed << std::get<4>(f) << "," + << std::setprecision(6) << std::fixed << (std::get<4>(f) - refLogL) + << "\n"; + } + } else { + os << "\n===== SUMMARY OF ALL FITS (converged fits ranked by LogL; failed fits listed unranked) =====\n"; + os << std::left + << std::setw(6) << "Rank" + << std::setw(6) << "#" + << std::setw(9) << "Success" + << std::setw(11) << "FitStatus" + << std::setw(9) << "eMatrix" + << std::setw(18) << "LogL" + << std::setw(14) << "DeltaLogL" + << "\n"; + int rank = 0; + for (size_t i = 0; i < fitLLs.size(); ++i) { + const auto& f = fitLLs[i]; + bool isSucc = std::get<1>(f); + if (isSucc) ++rank; + std::string rankStr = isSucc ? std::to_string(rank) : std::string("-"); + os << std::left + << std::setw(6) << rankStr + << std::setw(6) << std::get<0>(f) + << std::setw(9) << (isSucc ? "Y" : "N") + << std::setw(11) << std::get<2>(f) + << std::setw(9) << std::get<3>(f) + << std::setw(18) << (std::ostringstream() << std::setprecision(10) << std::fixed << std::get<4>(f)).str() + << std::setw(14) << (std::ostringstream() << std::setprecision(6) << std::fixed << (std::get<4>(f) - refLogL)).str() + << "\n"; + } + } + + // best fit callout (converged fits only) + if (!bestByLogL.empty()) { + os << "\nBEST FIT(S) by LogL among converged fits (" << bestByLogL.size() << " tied):\n"; + for (const auto& f : bestByLogL) + os << " Fit #" << std::get<0>(f) + << " success=" << (std::get<1>(f) ? "Y":"N") + << " status=" << std::get<2>(f) + << " eMatrix=" << std::get<3>(f) + << " LogL=" << std::setprecision(10) << std::fixed << std::get<4>(f) + << "\n"; + } else { + os << "\nBEST FIT(S) by LogL among converged fits: none (no fits converged)\n"; + } + + if (bestSuccIdx >= 0) + os << "BEST SUCCESSFUL FIT: #" << bestSuccIdx + << " LogL=" << std::setprecision(10) << std::fixed << bestSuccLogL << "\n"; + else + os << "BEST SUCCESSFUL FIT: none (all fits failed)\n"; + + // success rate + os << "\nSuccess: " << succ << " / " << fitLLs.size() + << " (" << std::setprecision(1) << std::fixed << pct << "%)\n"; + os << "=================================================\n"; + }; + + write(std::cout); + if (fout) write(fout); } double runSingleFit(ConfigurationInfo* cfgInfo, bool useMinos, bool hesse, int maxIter, string seedfile, int eMatrixRequirement) { @@ -236,6 +359,7 @@ void runRndFits(ConfigurationInfo* cfgInfo, bool useMinos, bool hesse, int maxIt gSystem->Exec(Form("cp %s_%d%s %s", seedfileBaseName.data(), minFitTag, seedfileExtension.data(), seedfile.data())); } // print summary of all fits + // summarizeFits(fitLLs,true); // write summary to CSV file summarizeFits(fitLLs); } } From 786be2d39b9c50a89c937fa9e22cf9501ae1c050 Mon Sep 17 00:00:00 2001 From: edbarriga Date: Wed, 8 Jul 2026 12:41:39 -0400 Subject: [PATCH 2/4] Updated name of variables. Class member data start with "m_" and local variables are camelCase --- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc | 80 +++++++------- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h | 6 +- .../AMPTOOLS_DATAIO/VecPsPlotGenerator.cc | 104 +++++++++--------- 3 files changed, 95 insertions(+), 95 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc index 549c7ab31..1cf5ab623 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc @@ -343,9 +343,9 @@ UserAmplitude< Vec_ps_refl >( args ){ m_s = inputArgs.s; m_polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; - polAngle = inputArgs.polAngle; - polFraction = inputArgs.polFraction; - polFrac_vs_E = nullptr; + m_polAngle = inputArgs.polAngle; + m_polFraction = inputArgs.polFraction; + m_polFracVsE = nullptr; m_3pi = inputArgs.omega3pi; m_gpi0 = inputArgs.omegagpi0; @@ -359,13 +359,13 @@ UserAmplitude< Vec_ps_refl >( args ){ if( !f || f->IsZombie() ) throw std::runtime_error( "[ Vec_ps_refl ]: could not open polarization file '" + inputArgs.polFile + "'" ); - polFrac_vs_E = (TH1D*)f->Get( inputArgs.polHist.c_str() ); - if( !polFrac_vs_E ) + m_polFracVsE = (TH1D*)f->Get( inputArgs.polHist.c_str() ); + if( !m_polFracVsE ) throw std::runtime_error( "[ Vec_ps_refl ]: histogram '" + inputArgs.polHist + "' not found in '" + inputArgs.polFile + "'" ); // Detach histogram from file so it persists after file is closed - polFrac_vs_E->SetDirectory(0); + m_polFracVsE->SetDirectory(0); f->Close(); delete f; } @@ -380,8 +380,8 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ TLorentzVector beam; TVector3 eps; // beam polarization vector (eps)ilon - GDouble beam_polFraction; - GDouble beam_polAngle; + GDouble beamPolFraction; + GDouble beamPolAngle; if(m_polInfoInPhotonP4){ // When pol info is stored in the photon 4-vector in the tree @@ -391,23 +391,23 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ // The values should be stored as px = polFraction*cos(polAngle) // and py = polFraction*sin(polAngle) eps.SetXYZ(pKin[0][1], pKin[0][2], 0.0); - beam_polFraction = eps.Mag(); - beam_polAngle = eps.Phi(); + beamPolFraction = eps.Mag(); + beamPolAngle = eps.Phi(); } else{ beam.SetPxPyPzE( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); - beam_polAngle = polAngle; + beamPolAngle = m_polAngle; // for fixed polarization fraction - if(polFrac_vs_E == nullptr) - beam_polFraction = polFraction; + if(m_polFracVsE == nullptr) + beamPolFraction = m_polFraction; else{ // for fitting with polarization vs E_gamma from input histogram - int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); - if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){ + int bin = m_polFracVsE->GetXaxis()->FindBin(pKin[0][0]); + if (bin == 0 || bin > m_polFracVsE->GetXaxis()->GetNbins()){ throw std::runtime_error( "[ Vec_ps_refl ]: energy " + std::to_string(pKin[0][0]) + " outside histogram range" ); } else - beam_polFraction = polFrac_vs_E->GetBinContent(bin); + beamPolFraction = m_polFracVsE->GetBinContent(bin); } } @@ -418,7 +418,7 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ TLorentzVector ps(pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0]); // Compute vector meson from its decay products // Make sure the order of daughters is correct in the config file! - TLorentzVector vec, vec_daught1, vec_daught2; + TLorentzVector vec, vecDaught1, vecDaught2; if(m_3pi){ // Omega ps proton, omega -> 3pi (6 particles): @@ -429,8 +429,8 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ TLorentzVector pip(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); TLorentzVector pim(pKin[5][1], pKin[5][2], pKin[5][3], pKin[5][0]); vec = pi0 + pip + pim; - vec_daught1 = pip; - vec_daught2 = pim; + vecDaught1 = pip; + vecDaught2 = pim; } else{ // Omega ps proton, omega -> gpi0 (5 particles): @@ -443,13 +443,13 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ // beam proton ps pi pi // Vec(Kpi) K+ Lambda, (5 particles) // beam proton ps K pi - vec_daught1 = TLorentzVector(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); - vec_daught2 = TLorentzVector(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); - vec = vec_daught1 + vec_daught2; + vecDaught1 = TLorentzVector(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); + vecDaught2 = TLorentzVector(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); + vec = vecDaught1 + vecDaught2; } // Final meson system P4 - TLorentzVector X = vec + ps; + TLorentzVector xMeson = vec + ps; ///////////////// Boost Particles and Get Angles///////////////////// @@ -459,42 +459,42 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ // Calculate decay angles for X in helicity frame (same for all vectors) // Change getXDecayAngles to get Gottfried-Jackson angles if needed // Note: it also calculates the production angle - vector xDecayAngles = getXDecayAngles( beam_polAngle, beam, beamTarget, X, vec); + vector xDecayAngles = getXDecayAngles( beamPolAngle, beam, beamTarget, xMeson, vec); // Calculate vector decay angles (unique for each vector) vector vectorDecayAngles; if(m_3pi){ - vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, - vec_daught1, vec_daught2); + vectorDecayAngles = getVectorDecayAngles( beamTarget, xMeson, vec, + vecDaught1, vecDaught2); } else{ - vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, - vec_daught1, TLorentzVector(0,0,0,0)); + vectorDecayAngles = getVectorDecayAngles( beamTarget, xMeson, vec, + vecDaught1, TLorentzVector(0,0,0,0)); } GDouble cosTheta = TMath::Cos(xDecayAngles[0]); GDouble phi = xDecayAngles[1]; - GDouble prod_angle = xDecayAngles[2]; // bigPhi + GDouble prodAngle = xDecayAngles[2]; // bigPhi GDouble cosThetaH = TMath::Cos(vectorDecayAngles[0]); GDouble phiH = vectorDecayAngles[1]; - GDouble m_X = X.M(); - GDouble m_vec = vec.M(); - GDouble m_ps = ps.M(); + GDouble xMesonMass = xMeson.M(); + GDouble vecMass = vec.M(); + GDouble psMass = ps.M(); complex amplitude(0,0); static const complex i(0,1); if(m_gpi0){ // radiative omega decay requires a handling of the photon helicity for (int lambda = -1; lambda <= 1; lambda++) { // sum over vector helicity - GDouble hel_amp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); + GDouble helAmp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); if(lambda==0){ amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * - hel_amp * conj(wignerD(1, lambda, m_ghel, cosThetaH, phiH)) * + helAmp * conj(wignerD(1, lambda, m_ghel, cosThetaH, phiH)) * (m_ghel*1.0); // m_ghel is int } else{ amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * - hel_amp * conj(wignerD(1, -1*lambda, m_ghel, cosThetaH, phiH)) * + helAmp * conj(wignerD(1, -1*lambda, m_ghel, cosThetaH, phiH)) * -1.0; // the power of lambda is irrelevant for lambda =/= 0: (-1)^lambda = -1 } @@ -502,20 +502,20 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ } else{ // for any other vector decay for (int lambda = -1; lambda <= 1; lambda++) { // sum over vector helicity - GDouble hel_amp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); + GDouble helAmp = clebschGordan(m_l, 1, 0, lambda, m_j, lambda); amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi)) * - hel_amp * conj(wignerD(1, lambda, 0, cosThetaH, phiH)); + helAmp * conj(wignerD(1, lambda, 0, cosThetaH, phiH)); } } // The amplitude is multiplied by a factor, either sqrt(1 + -P_gamma) or // sqrt(1 + P_gamma) depending on the what sum is being calculated - GDouble factor = sqrt(1 + m_s * beam_polFraction); + GDouble factor = sqrt(1 + m_s * beamPolFraction); // The result of the function that depends on the angles is stored in zjm complex zjm = 0; // A - sign translates to + in the prod_angle // This is because reflectivity conventions introduce exp(-iPhi) - complex rotateY = polar((GDouble)1., (GDouble)(-1. * prod_angle )); + complex rotateY = polar((GDouble)1., (GDouble)(-1. * prodAngle )); if (m_r == 1) zjm = real(amplitude * rotateY); @@ -523,7 +523,7 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ zjm = i*imag(amplitude * rotateY); if( !m_noBarrier ){ - GDouble kinFactor = barrierFactor(m_X, m_l, m_vec, m_ps); + GDouble kinFactor = barrierFactor(xMesonMass, m_l, vecMass, psMass); factor *= kinFactor; // Alternatives: // E852 Nozar thesis has sqrt(2*s+1)*sqrt(2*l+1)*F_l(p_omega)*sqrt(omega) diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h index 052f184c6..07e822e0a 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h @@ -82,9 +82,9 @@ class Vec_ps_refl : public UserAmplitude< Vec_ps_refl > int m_s; bool m_polInfoInPhotonP4; - GDouble polFraction; - GDouble polAngle; - TH1D *polFrac_vs_E; + GDouble m_polFraction; + GDouble m_polAngle; + TH1D *m_polFracVsE; bool m_3pi; bool m_gpi0; diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc index 35fd3d499..2926c6e04 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc @@ -141,7 +141,7 @@ void VecPsPlotGenerator::createHistograms( ) { bookHistogram(kPhiH, new Histogram1D(50, -PI, PI, "Phi_H", ";#phi_H [rad.]")); bookHistogram(kProd_Ang, new Histogram1D(50, -PI, PI, "Prod_Ang", ";Prod_Ang (#Phi) [rad.]")); - bookHistogram(kProdOffset, new Histogram1D(50, -PI, PI, "ProdOffset", + bookHistogram(kProdOffset, new Histogram1D(50, -PI, PI, "ProdOffset", ";Prod_Ang (#Phi) Uncorrected [rad.]" ) ); bookHistogram(kt, new Histogram1D(100, 0, 2.0, "t", ";-t")); bookHistogram(kRecoilMass, new Histogram1D(100, 0.9, 1.9, "MRecoil", @@ -178,10 +178,10 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ cfgInfo()->amplitudeList(reactionName, "", "").at(0)->factors().at(0); Vec_ps_plot_Args inputArgs = parsedArgs( args ); - bool m_polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; - double beam_polAngle = inputArgs.polAngle; + bool polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; + double beamPolAngle = inputArgs.polAngle; - bool m_3pi = inputArgs.omega3pi; + bool omega3pi = inputArgs.omega3pi; //cout << "project event" << endl; TLorentzVector readBeam = kin->particle( 0 ); @@ -191,7 +191,7 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ TVector3 eps; // beam polarization vector (eps)ilon TLorentzVector beam; - if(m_polInfoInPhotonP4){ + if(polInfoInPhotonP4){ // When pol info is stored in the photon 4-vector in the tree // the energy (pKin[0][0]) and the pz (pKin[0][3]) are used as normal beam.SetPxPyPzE( 0.0, 0.0, readBeam.Pz(), readBeam.E() ); @@ -199,23 +199,23 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ // The values should be stored as px = polFraction*cos(polAngle) // and py = polFraction*sin(polAngle) eps.SetXYZ(readBeam.Px(), readBeam.Py(), 0.0); - beam_polAngle = eps.Phi(); + beamPolAngle = eps.Phi(); } else{ beam = readBeam; } // Min particle index for recoil sum - int min_recoil = 5; + int minRecoil = 5; // Compute vector meson from its decay products // Make sure the order of daughters is correct in the config file! - TLorentzVector vec, vec_daught1, vec_daught2; - double dalitz_s, dalitz_t, dalitz_u, dalitz_d, dalitz_sc; - double dalitz_x = 0.0; - double dalitz_y = 0.0; + TLorentzVector vec, vecDaught1, vecDaught2; + double dalitzS, dalitzT, dalitzU, dalitzD, dalitzSC; + double dalitzX = 0.0; + double dalitzY = 0.0; - if(m_3pi) { + if(omega3pi) { // Omega ps proton, omega -> 3pi (6 particles): // beam proton ps pi0 pip pim // Omega pi- Delta++, omega -> 3pi (6 particles): @@ -224,9 +224,9 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ TLorentzVector pip = kin->particle( 4 ); TLorentzVector pim = kin->particle( 5 ); vec = pi0 + pip + pim; - vec_daught1 = pip; - vec_daught2 = pim; - min_recoil = 6; + vecDaught1 = pip; + vecDaught2 = pim; + minRecoil = 6; // Dalitz variables are included here. // They amplitude tends to be multiply in the config file but that's not @@ -234,18 +234,18 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ const auto& p2 = pi0; const auto& p3 = pip; const auto& p4 = pim; - double pdg_m_pi0 = 0.1349766; - double pdg_m_chargedPi = 0.13957018; - dalitz_s = (p3 + p4).M2(); // s=M(pip pim) - dalitz_t = (p2 + p3).M2(); // s=M(pip pi0) - dalitz_u = (p2 + p4).M2(); // s=M(pim pi0) - dalitz_d = 2 * (p2 + p3 + p4).M() * - ((p2 + p3 + p4).M() - ((2 * pdg_m_chargedPi) + pdg_m_pi0)); - dalitz_sc = (1.0 / 3.0) * ((p2 + p3 + p4).M2() + - ((2 * (pdg_m_chargedPi * pdg_m_chargedPi)) + - (pdg_m_pi0 * pdg_m_pi0))); - dalitz_x = sqrt(3.0) * (dalitz_t - dalitz_u) / dalitz_d; - dalitz_y = 3.0 * (dalitz_sc - dalitz_s) / dalitz_d; + double pdgMassPi0 = 0.1349766; + double pdgMassPi = 0.13957018; + dalitzS = (p3 + p4).M2(); // s=M(pip pim) + dalitzT = (p2 + p3).M2(); // s=M(pip pi0) + dalitzU = (p2 + p4).M2(); // s=M(pim pi0) + dalitzD = 2 * (p2 + p3 + p4).M() * + ((p2 + p3 + p4).M() - ((2 * pdgMassPi) + pdgMassPi0)); + dalitzSC = (1.0 / 3.0) * ((p2 + p3 + p4).M2() + + ((2 * (pdgMassPi * pdgMassPi)) + + (pdgMassPi0 * pdgMassPi0))); + dalitzX = sqrt(3.0) * (dalitzT - dalitzU) / dalitzD; + dalitzY = 3.0 * (dalitzSC - dalitzS) / dalitzD; } else { // Omega ps proton, omega -> gpi0 (5 particles): @@ -258,21 +258,21 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ // beam proton ps pi pi // Vec(Kpi) K+ Lambda, (5 particles) // beam proton ps K pi - vec_daught1 = kin->particle( 3 ); - vec_daught2 = kin->particle( 4 ); - vec = vec_daught1 + vec_daught2; - min_recoil = 5; + vecDaught1 = kin->particle( 3 ); + vecDaught2 = kin->particle( 4 ); + vec = vecDaught1 + vecDaught2; + minRecoil = 5; } // Final meson system P4 TLorentzVector X = vec + ps; - TLorentzVector proton_ps = recoil + ps; - TLorentzVector recoil_ps = proton_ps; - for(uint i=min_recoil; iparticleList().size(); i++){ + TLorentzVector protonPs = recoil + ps; + TLorentzVector recoilPs = protonPs; + for(uint i=minRecoil; iparticleList().size(); i++){ // Add mesons to recoil system (e.g. Delta or Lambda) recoil += kin->particle(i); - recoil_ps += kin->particle(i); + recoilPs += kin->particle(i); } TLorentzVector target(0,0,0,0.938); // proton at rest @@ -281,28 +281,28 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ // Calculate decay angles for X in helicity frame (same for all vectors) // Change getXDecayAngles to get Gottfried-Jackson angles if needed // Note: it also calculates the production angle - vector xDecayAngles = getXDecayAngles( beam_polAngle, beam, beamTarget, X, vec); + vector xDecayAngles = getXDecayAngles( beamPolAngle, beam, beamTarget, X, vec); // set polarization angle to zero to see shift in Phi_Prod distributions - vector xDecayAngles_offset = getXDecayAngles( 0, beam, beamTarget, X, vec); + vector xDecayAnglesOffset = getXDecayAngles( 0, beam, beamTarget, X, vec); // Calculate vector decay angles (unique for each vector) vector vectorDecayAngles; - if(m_3pi){ + if(omega3pi){ vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, - vec_daught1, vec_daught2); + vecDaught1, vecDaught2); } else{ vectorDecayAngles = getVectorDecayAngles( beamTarget, X, vec, - vec_daught1, TLorentzVector(0,0,0,0)); + vecDaught1, TLorentzVector(0,0,0,0)); } - double Mandt = fabs((target-recoil).M2()); - double recoil_mass = recoil.M2(); + double momentumTransfer = fabs((target-recoil).M2()); + double recoilMass = recoil.M2(); GDouble cosTheta = TMath::Cos(xDecayAngles[0]); GDouble phi = xDecayAngles[1]; - GDouble prod_angle = xDecayAngles[2]; // bigPhi - GDouble prod_angle_offset = xDecayAngles_offset[2]; // bigPhi not corrected + GDouble prodAngle = xDecayAngles[2]; // bigPhi + GDouble prodAngleOffset = xDecayAnglesOffset[2]; // bigPhi not corrected GDouble cosThetaH = TMath::Cos(vectorDecayAngles[0]); GDouble phiH = vectorDecayAngles[1]; GDouble lambda = vectorDecayAngles[2]; @@ -313,13 +313,13 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ fillHistogram( kPhi, phi ); fillHistogram( kCosThetaH, cosThetaH ); fillHistogram( kPhiH, phiH ); - fillHistogram( kProd_Ang, prod_angle ); - fillHistogram( kt, Mandt ); - fillHistogram( kRecoilMass, recoil_mass ); - fillHistogram( kProtonPsMass, proton_ps.M() ); - fillHistogram( kRecoilPsMass, recoil_ps.M() ); + fillHistogram( kProd_Ang, prodAngle ); + fillHistogram( kt, momentumTransfer ); + fillHistogram( kRecoilMass, recoilMass ); + fillHistogram( kProtonPsMass, protonPs.M() ); + fillHistogram( kRecoilPsMass, recoilPs.M() ); fillHistogram( kLambda, lambda ); - fillHistogram( kDalitz, dalitz_x, dalitz_y ); - fillHistogram( kPhi_ProdVsPhi, phi, prod_angle ); - fillHistogram( kPhiOffsetVsPhi, phi, prod_angle_offset ); + fillHistogram( kDalitz, dalitzX, dalitzY ); + fillHistogram( kPhi_ProdVsPhi, phi, prodAngle ); + fillHistogram( kPhiOffsetVsPhi, phi, prodAngleOffset ); } From 563d643229e8dcc42dd5ba56009a3d8f02893d85 Mon Sep 17 00:00:00 2001 From: edbarriga Date: Fri, 10 Jul 2026 09:28:46 -0400 Subject: [PATCH 3/4] The plotter class now does parsing in the constructor and re-uses parsing functions from the amplitude class --- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc | 360 ++++++++---------- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h | 63 +++ .../AMPTOOLS_DATAIO/VecPsPlotGenerator.cc | 134 ++----- .../AMPTOOLS_DATAIO/VecPsPlotGenerator.h | 6 + 4 files changed, 255 insertions(+), 308 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc index 1cf5ab623..c3befdcf0 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc @@ -130,212 +130,162 @@ #include "UTILITIES/BeamProperties.h" -//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -//....oooOO0OOooo........ Structs ........oooOO0OOooo..... -//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -struct Vec_ps_refl_Args { - // Required positional args - int j; // Total spin - int m; // Spin projection - int l; // Relative orbital angular momentum - int r; // +1/-1 for (r)eal/imaginary part of the amplitude - int s; // +1/-1 (s)ign of the P_gamma term - - // Polarization information - bool polInfoInPhotonP4 = true; - double polAngle = -2.0; // -2 = not set - double polFraction = -2.0; // -2 = not set - // Amorphous sometimes is referred as having a polAngle -1, so we use -2 to - // avoid possible confusion for both polAngle and polFraction - std::string polFile = ""; // path to .root file with polarization vs E_beam - std::string polHist = ""; // name of TH1D inside polFile - - // Default is 2-body vector decay, but if omega is used as vector: - // Omega decay mode - bool omega3pi = false; - bool omegagpi0 = false; - int gHelicity = 0; // +1/-1 helicity of the bachelor photon in omega->g pi0 - - bool noBarrier = false; // Turn off barrier factors -}; -//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -//....oooOO0OOooo........ Helper Functions ........oooOO0OOooo..... -//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -static bool isValidNumber(const std::string& argInput, double &value){ - char* end = nullptr; - errno = 0; // reset global error - value = std::strtod(argInput.c_str(), &end); - - // Check if - // (1) no conversion was performed - // (2) there are leftover characters - // (3) an overflow/underflow occurred - if(end == argInput.c_str() || *end != '\0' || errno != 0) { - return false; // not a valid number - } - // If "end" points to the end of string, it's fully numeric - return true; -} - -static double parseValidatedNumber(const std::string& label, const std::string& argInput){ - double tmpValue = 0.0; - if(!isValidNumber(argInput, tmpValue)){ - throw std::invalid_argument("[ Vec_ps_refl ]: invalid " + label + ": " + argInput); - } - return tmpValue; -} - -static Vec_ps_refl_Args parsedArgs( const std::vector& args ){ - Vec_ps_refl_Args inputArgs; - - // Required positional args - inputArgs.j = static_cast( parseValidatedNumber( "J", args[0] ) ); - inputArgs.m = static_cast( parseValidatedNumber( "M", args[1] ) ); - inputArgs.l = static_cast( parseValidatedNumber( "L", args[2] ) ); - inputArgs.r = static_cast( parseValidatedNumber( "Re/Im", args[3] ) ); - inputArgs.s = static_cast( parseValidatedNumber( "P_gamma sign", args[4] ) ); - - // Physics validation of required args - if( abs(inputArgs.m) > inputArgs.j ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: |M| must be <= J, got J=" + args[0] + " M=" + args[1] ); - if( abs(inputArgs.r) != 1 ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: Re/Im flag must be +1 or -1, got " + args[3] ); - if( abs(inputArgs.s) != 1 ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: P_gamma sign must be +1 or -1, got " + args[4] ); - - // If no optional args, return now - if( args.size() <= 5 ) return inputArgs; - - // Detect format by checking if args[5] contains '=' - bool isKeyValue = std::any_of( args.begin() + 5, args.end(), - []( const std::string& s ){ return s.find('=') != std::string::npos; } ); - - // Lambda function to detect bare flags (valid in both formats) - auto isBareFlag = []( const std::string& s ){ - return s == "omega3pi" || s == "omegagpi0" || s == "nobarrier"; - }; - - if( !isKeyValue ){ - size_t i = 5; - - // arg[5]: polarization angle - if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ - inputArgs.polAngle = parseValidatedNumber( "polAngle", args[i++] ); - inputArgs.polInfoInPhotonP4 = false; - } - - // arg[6]: fixed fraction or .root file path - if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ - if( args[i].find(".root") != std::string::npos ){ - inputArgs.polFile = args[i++]; - if( i >= args.size() || isBareFlag(args[i]) ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: .root polarization file requires a histogram name immediately after it" ); - inputArgs.polHist = args[i++]; - } else { - inputArgs.polFraction = parseValidatedNumber( "polFraction", args[i++] ); - } - } - - // Remaining: bare flags, and gHelicity bare integer after omegagpi0 - while( i < args.size() ){ - if( args[i] == "omega3pi" ){ - inputArgs.omega3pi = true; - i++; - } - else if( args[i] == "omegagpi0" ){ - inputArgs.omegagpi0 = true; - if( i + 1 >= args.size() ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: omegagpi0 requires a photon helicity value (+1 or -1) immediately after it" ); - inputArgs.gHelicity = static_cast( parseValidatedNumber( "gHelicity", args[i+1] ) ); - i += 2; - } - else if( args[i] == "nobarrier" ){ - inputArgs.noBarrier = true; - i++; - } - else{ - throw std::invalid_argument( - "[ Vec_ps_refl ]: unrecognized positional argument '" + args[i] + "'" ); - } - } - } - else { // key=value format - for( size_t i = 5; i < args.size(); i++ ){ - const std::string& arg = args[i]; - - if( arg == "omega3pi" ){ inputArgs.omega3pi = true; continue; } - if( arg == "omegagpi0" ){ inputArgs.omegagpi0 = true; continue; } - if( arg == "nobarrier" ){ inputArgs.noBarrier = true; continue; } - auto sep = arg.find('='); // '=' is used as separator - if( sep == std::string::npos ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: unrecognized argument '" + arg + - "' (expected key=value pair or known flag)" ); - - const std::string key = arg.substr(0, sep); - const std::string value = arg.substr(sep + 1); - - if( key == "polAngle" ){ - inputArgs.polAngle = parseValidatedNumber( "polAngle", value ); - inputArgs.polInfoInPhotonP4 = false; - } - else if( key == "polFraction" ){ - inputArgs.polFraction = parseValidatedNumber( "polFraction", value ); - } - else if( key == "polFile" ){ inputArgs.polFile = value; } - else if( key == "polHist" ){ inputArgs.polHist = value; } - else if( key == "gHelicity" ){ - inputArgs.gHelicity = static_cast( parseValidatedNumber( "gHelicity", value ) ); - } - else{ - throw std::invalid_argument( - "[ Vec_ps_refl ]: unrecognized key '" + key + "'" ); - } - } - } - - // Post-parse validation - if( inputArgs.polFile != "" && inputArgs.polHist == "" ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: polFile='" + inputArgs.polFile + - "' requires polHist= to also be specified" ); - - if( inputArgs.polHist != "" && inputArgs.polFile == "" ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: polHist='" + inputArgs.polHist + - "' was specified without polFile=" ); - - if( inputArgs.polFraction >= 0.0 && inputArgs.polFile != "" ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: polFraction and polFile are mutually exclusive; specify only one" ); - - if( inputArgs.omegagpi0 && abs(inputArgs.gHelicity) != 1 ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: omegagpi0 requires gHelicity=+1 or gHelicity=-1" ); - - if( inputArgs.gHelicity != 0 && !inputArgs.omegagpi0 ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: gHelicity was specified but omegagpi0 flag is missing" ); - - if( inputArgs.omega3pi && inputArgs.omegagpi0 ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: omega3pi and omegagpi0 are mutually exclusive; specify only one" ); - - if( !inputArgs.polInfoInPhotonP4 && inputArgs.polFraction < 0.0 && inputArgs.polFile == "" ) - throw std::invalid_argument( - "[ Vec_ps_refl ]: polAngle was given but neither polFraction nor polFile+polHist were provided" ); - - return inputArgs; -} + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + //....oooOO0OOooo........ Helper Functions ........oooOO0OOooo..... + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Vec_ps_refl::VecPsReflArgs Vec_ps_refl::parsedArgs( + const std::vector& args, + const std::string& context){ + VecPsReflArgs inputArgs; + + // Required positional args + inputArgs.j = static_cast( parseValidatedNumber( "J", args[0], context ) ); + inputArgs.m = static_cast( parseValidatedNumber( "M", args[1], context ) ); + inputArgs.l = static_cast( parseValidatedNumber( "L", args[2], context ) ); + inputArgs.r = static_cast( parseValidatedNumber( "Re/Im", args[3], context ) ); + inputArgs.s = static_cast( parseValidatedNumber( "P_gamma sign", args[4], context ) ); + + // Physics validation of required args + if( abs(inputArgs.m) > inputArgs.j ) + throw std::invalid_argument( + "[ " + context + " ] |M| must be <= J, got J=" + args[0] + " M=" + args[1] ); + if( abs(inputArgs.r) != 1 ) + throw std::invalid_argument( + "[ " + context + " ] Re/Im flag must be +1 or -1, got " + args[3] ); + if( abs(inputArgs.s) != 1 ) + throw std::invalid_argument( + "[ " + context + " ] P_gamma sign must be +1 or -1, got " + args[4] ); + + // If no optional args, return now + if( args.size() <= 5 ) return inputArgs; + + // Detect format by checking if args[5] contains '=' + bool isKeyValue = std::any_of( args.begin() + 5, args.end(), + []( const std::string& s ){ return s.find('=') != std::string::npos; } ); + + // Lambda function to detect bare flags (valid in both formats) + auto isBareFlag = []( const std::string& s ){ + return s == "omega3pi" || s == "omegagpi0" || s == "nobarrier"; + }; + + if( !isKeyValue ){ + size_t i = 5; + + // arg[5]: polarization angle + if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", args[i++], context ); + inputArgs.polInfoInPhotonP4 = false; + } + + // arg[6]: fixed fraction or .root file path + if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ + if( args[i].find(".root") != std::string::npos ){ + inputArgs.polFile = args[i++]; + if( i >= args.size() || isBareFlag(args[i]) ) + throw std::invalid_argument( + "[ " + context + " ] .root polarization file requires a histogram name immediately after it" ); + inputArgs.polHist = args[i++]; + } else { + inputArgs.polFraction = parseValidatedNumber( "polFraction", args[i++], context ); + } + } + + // Remaining: bare flags, and gHelicity bare integer after omegagpi0 + while( i < args.size() ){ + if( args[i] == "omega3pi" ){ + inputArgs.omega3pi = true; + i++; + } + else if( args[i] == "omegagpi0" ){ + inputArgs.omegagpi0 = true; + if( i + 1 >= args.size() ) + throw std::invalid_argument( + "[ " + context + " ] omegagpi0 requires a photon helicity value (+1 or -1) immediately after it" ); + inputArgs.gHelicity = static_cast( parseValidatedNumber( "gHelicity", args[i+1], context ) ); + i += 2; + } + else if( args[i] == "nobarrier" ){ + inputArgs.noBarrier = true; + i++; + } + else{ + throw std::invalid_argument( + "[ " + context + " ] unrecognized positional argument '" + args[i] + "'" ); + } + } + } + else { // key=value format + for( size_t i = 5; i < args.size(); i++ ){ + const std::string& arg = args[i]; + + if( arg == "omega3pi" ){ inputArgs.omega3pi = true; continue; } + if( arg == "omegagpi0" ){ inputArgs.omegagpi0 = true; continue; } + if( arg == "nobarrier" ){ inputArgs.noBarrier = true; continue; } + auto sep = arg.find('='); // '=' is used as separator + if( sep == std::string::npos ) + throw std::invalid_argument( + "[ " + context + " ] unrecognized argument '" + arg + + "' (expected key=value pair or known flag)" ); + + const std::string key = arg.substr(0, sep); + const std::string value = arg.substr(sep + 1); + + if( key == "polAngle" ){ + inputArgs.polAngle = parseValidatedNumber( "polAngle", value, context ); + inputArgs.polInfoInPhotonP4 = false; + } + else if( key == "polFraction" ){ + inputArgs.polFraction = parseValidatedNumber( "polFraction", value, context ); + } + else if( key == "polFile" ){ inputArgs.polFile = value; } + else if( key == "polHist" ){ inputArgs.polHist = value; } + else if( key == "gHelicity" ){ + inputArgs.gHelicity = static_cast( parseValidatedNumber( "gHelicity", value, context ) ); + } + else{ + throw std::invalid_argument( + "[ " + context + " ] unrecognized key '" + key + "'" ); + } + } + } + + // Post-parse validation + if( inputArgs.polFile != "" && inputArgs.polHist == "" ) + throw std::invalid_argument( + "[ " + context + " ] polFile='" + inputArgs.polFile + + "' requires polHist= to also be specified" ); + + if( inputArgs.polHist != "" && inputArgs.polFile == "" ) + throw std::invalid_argument( + "[ " + context + " ] polHist='" + inputArgs.polHist + + "' was specified without polFile=" ); + + if( inputArgs.polFraction >= 0.0 && inputArgs.polFile != "" ) + throw std::invalid_argument( + "[ " + context + " ] polFraction and polFile are mutually exclusive; specify only one" ); + + if( inputArgs.omegagpi0 && abs(inputArgs.gHelicity) != 1 ) + throw std::invalid_argument( + "[ " + context + " ] omegagpi0 requires gHelicity=+1 or gHelicity=-1" ); + + if( inputArgs.gHelicity != 0 && !inputArgs.omegagpi0 ) + throw std::invalid_argument( + "[ " + context + " ] gHelicity was specified but omegagpi0 flag is missing" ); + + if( inputArgs.omega3pi && inputArgs.omegagpi0 ) + throw std::invalid_argument( + "[ " + context + " ] omega3pi and omegagpi0 are mutually exclusive; specify only one" ); + + if( !inputArgs.polInfoInPhotonP4 && inputArgs.polFraction < 0.0 && inputArgs.polFile == "" ) + throw std::invalid_argument( + "[ " + context + " ] polAngle was given but neither polFraction nor polFile+polHist were provided" ); + + return inputArgs; + } Vec_ps_refl::Vec_ps_refl( const vector< string >& args ) : UserAmplitude< Vec_ps_refl >( args ){ - Vec_ps_refl_Args inputArgs = parsedArgs( args ); + VecPsReflArgs inputArgs = parsedArgs( args ); m_j = inputArgs.j; m_m = inputArgs.m; m_l = inputArgs.l; @@ -343,8 +293,8 @@ UserAmplitude< Vec_ps_refl >( args ){ m_s = inputArgs.s; m_polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; - m_polAngle = inputArgs.polAngle; - m_polFraction = inputArgs.polFraction; + m_polAngle = inputArgs.polAngle; + m_polFraction = inputArgs.polFraction; m_polFracVsE = nullptr; m_3pi = inputArgs.omega3pi; diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h index 07e822e0a..756444f3d 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.h @@ -65,6 +65,69 @@ class Vec_ps_refl : public UserAmplitude< Vec_ps_refl > void updatePar( const AmpParameter& par ); + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + //....oooOO0OOooo........ Structs ........oooOO0OOooo..... + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + struct VecPsReflArgs { + // Required positional args + int j; // Total spin + int m; // Spin projection + int l; // Relative orbital angular momentum + int r; // +1/-1 for (r)eal/imaginary part of the amplitude + int s; // +1/-1 (s)ign of the P_gamma term + + // Polarization information + bool polInfoInPhotonP4 = true; + double polAngle = -2.0; // -2 = not set + double polFraction = -2.0; // -2 = not set + // Amorphous sometimes is referred as having a polAngle -1, so we use -2 to + // avoid possible confusion for both polAngle and polFraction + std::string polFile = ""; // path to .root file with polarization vs E_beam + std::string polHist = ""; // name of TH1D inside polFile + + // Default is 2-body vector decay, but if omega is used as vector: + // Omega decay mode + bool omega3pi = false; + bool omegagpi0 = false; + int gHelicity = 0; // +1/-1 helicity of the bachelor photon in omega->g pi0 + + bool noBarrier = false; // Turn off barrier factors + }; + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + //....oooOO0OOooo........ Helper Functions ........oooOO0OOooo..... + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + static bool isValidNumber(const std::string& argInput, double &value){ + char* end = nullptr; + errno = 0; // reset global error + value = std::strtod(argInput.c_str(), &end); + + // Check if + // (1) no conversion was performed + // (2) there are leftover characters + // (3) an overflow/underflow occurred + if(end == argInput.c_str() || *end != '\0' || errno != 0) { + return false; // not a valid number + } + // If "end" points to the end of string, it's fully numeric + return true; + } + + static double parseValidatedNumber(const std::string& label, + const std::string& argInput, + const std::string& context = "Vec_ps_refl"){ + double tmpValue = 0.0; + if(!isValidNumber(argInput, tmpValue)){ + throw std::invalid_argument("[ " + context + " ] invalid " + + label + ": " + argInput); + } + return tmpValue; + } + + static VecPsReflArgs parsedArgs(const std::vector& args, + const std::string& context = "Vec_ps_refl"); + #ifdef GPU_ACCELERATION void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const; diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc index 2926c6e04..ae8e8b43b 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc @@ -9,115 +9,47 @@ #include "IUAmpTools/Histogram1D.h" #include "IUAmpTools/Kinematics.h" -//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -//....oooOO0OOooo........ Structs ........oooOO0OOooo..... -//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -struct Vec_ps_plot_Args { - bool isMoment = false; // default Vec_ps_refl, true if Vec_ps_moment detected - - // Polarization information - bool polInfoInPhotonP4 = true; - double polAngle = -2.0; // -2 = not set - // Amorphous sometimes is referred as having a polAngle -1, so we use -2 to - // avoid possible confusion for both polAngle - - // Default is 2-body vector decay, but if omega is used as vector: - // Omega decay mode - bool omega3pi = false; -}; //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ //....oooOO0OOooo........ Helper Functions ........oooOO0OOooo..... //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -static bool isValidNumber(const string& argInput, double &value){ - char* end = nullptr; - errno = 0; // reset global error - value = std::strtod(argInput.c_str(), &end); - - // Check if - // (1) no conversion was performed - // (2) there are leftover characters - // (3) an overflow/underflow occurred - if(end == argInput.c_str() || *end != '\0' || errno != 0) { - return false; // not a valid number - } - // If end points to the end of string, it's fully numeric - return true; -} - -static double parseValidatedNumber(const string& label, const string& argInput){ - double tmpValue = 0.0; - if(!isValidNumber(argInput, tmpValue)){ - throw std::invalid_argument("[ VecPsPlotGenerator ]: invalid " + label + ": " + argInput); - } - return tmpValue; -} - -static Vec_ps_plot_Args parsedArgs( const std::vector& args ){ - Vec_ps_plot_Args inputArgs; - - inputArgs.isMoment = ( args[0] == "Vec_ps_moment" ); - if( inputArgs.isMoment ){ - inputArgs.polAngle = parseValidatedNumber( "polAngle", args[1] ); - inputArgs.omega3pi = true; - inputArgs.polInfoInPhotonP4 = false; - return inputArgs; - } - - // If no optional args, return now - if( args.size() <= 6 ) return inputArgs; - - // Detect format by checking if args[6] contains '=' - bool isKeyValue = std::any_of( args.begin() + 6, args.end(), - []( const std::string& s ){ return s.find('=') != std::string::npos; } ); - - // Lambda function to detect bare flags (valid in both formats) - auto isBareFlag = []( const std::string& s ){ - return s == "omega3pi" || s == "omegagpi0" || s == "nobarrier"; - }; - - if( !isKeyValue ){ - size_t i = 6; - - // arg[6]: polarization angle - if( i < args.size() && !isBareFlag(args[i]) && args[i].find('=') == std::string::npos ){ - inputArgs.polAngle = parseValidatedNumber( "polAngle", args[i++] ); - inputArgs.polInfoInPhotonP4 = false; - } - - // Remaining: bare flags, and gHelicity bare integer after omegagpi0 - for( size_t i = 6; i < args.size(); i++ ){ - if( args[i] == "omega3pi" ){ +void VecPsPlotGenerator::cacheArgs(){ + + const auto& reactionList = cfgInfo()->reactionList(); + for( const auto& reaction : reactionList ){ + const std::string reactionName = reaction->reactionName(); + const std::vector& fullAmplitudeArgs = + cfgInfo()->amplitudeList(reactionName, "", "").at(0)->factors().at(0); + + // Remove the amplitude name since the Vec_ps_refl parser expects only + // the arguments of the amplitude + std::vector args(fullAmplitudeArgs.begin() + 1, + fullAmplitudeArgs.end()); + + Vec_ps_refl::VecPsReflArgs inputArgs; + const std::string context = "VecPsPlotGenerator"; + // The first argument in the list is Vec_ps_refl or Vec_ps_moment + if( fullAmplitudeArgs[0] == "Vec_ps_moment" ){ + inputArgs.polAngle = Vec_ps_refl::parseValidatedNumber( "polAngle", args[1], context ); inputArgs.omega3pi = true; + inputArgs.polInfoInPhotonP4 = false; } - } - } - else { // key=value format - for( size_t i = 6; i < args.size(); i++ ){ - const std::string& arg = args[i]; - - if( arg == "omega3pi" ){ inputArgs.omega3pi = true; continue; } - auto sep = arg.find('='); // '=' is used as separator - if( sep == std::string::npos ) + else if( fullAmplitudeArgs[0] == "Vec_ps_refl" ){ + inputArgs = Vec_ps_refl::parsedArgs( args , context ); + } + else{ throw std::invalid_argument( - "[ Vec_ps_refl ]: unrecognized argument '" + arg + - "' (expected key=value pair or known flag)" ); - - const std::string key = arg.substr(0, sep); - const std::string value = arg.substr(sep + 1); - - if( key == "polAngle" ){ - inputArgs.polAngle = parseValidatedNumber( "polAngle", value ); - inputArgs.polInfoInPhotonP4 = false; + "[ " + context + " ] unrecognized amplitude '" + fullAmplitudeArgs[0] + + "' (expected 'Vec_ps_refl' or 'Vec_ps_moment')" ); } + m_args.emplace(reactionName, inputArgs); } } - return inputArgs; -} /* Constructor to display FitResults */ VecPsPlotGenerator::VecPsPlotGenerator( const FitResults& results, Option opt ) : PlotGenerator( results, opt ) { + cacheArgs(); createHistograms(); } @@ -173,15 +105,11 @@ VecPsPlotGenerator::projectEvent( Kinematics* kin ){ void VecPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ - // The first argument in the list is Vec_ps_refl or Vec_ps_moment - const vector args = - cfgInfo()->amplitudeList(reactionName, "", "").at(0)->factors().at(0); - Vec_ps_plot_Args inputArgs = parsedArgs( args ); + const auto& inputArgs = m_args.at(reactionName); - bool polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; - double beamPolAngle = inputArgs.polAngle; - - bool omega3pi = inputArgs.omega3pi; + bool polInfoInPhotonP4 = inputArgs.polInfoInPhotonP4; + double beamPolAngle = inputArgs.polAngle; + bool omega3pi = inputArgs.omega3pi; //cout << "project event" << endl; TLorentzVector readBeam = kin->particle( 0 ); diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h index cd45365c6..280286f30 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h @@ -4,6 +4,8 @@ #include #include +#include "AMPTOOLS_AMPS/Vec_ps_refl.h" + #include "IUAmpTools/PlotGenerator.h" using namespace std; @@ -54,6 +56,10 @@ class VecPsPlotGenerator : public PlotGenerator void projectEvent( Kinematics* kin, const string& reactionName ); void createHistograms( ); + + std::map m_args; + + void cacheArgs(); }; From 036881f961ca42ab1263bacf007f895cd85a5d90 Mon Sep 17 00:00:00 2001 From: edbarriga Date: Thu, 16 Jul 2026 13:43:09 -0400 Subject: [PATCH 4/4] Corrected the usage example --- src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc index c3befdcf0..d6aa7225a 100644 --- a/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc +++ b/src/libraries/AMPTOOLS_AMPS/Vec_ps_refl.cc @@ -87,16 +87,15 @@ // For a 2P m=0 wave from an omega->gpi that belongs to the real part of // the amplitude sum,has a negative sign in front of the polarization term, and // a photon helicity -1, you would write: -// you would write: // If polarization info stored in the beam P4: -// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 omegagpi -1 -// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 omegagpi gHelicity=-1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 2 0 1 1 -1 omegagpi -1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 2 0 1 1 -1 omegagpi gHelicity=-1 // If polarization angle and fraction are fixed: -// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 90 .35 omegagpi -1 -// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 polFraction=0.35 polAngle=90 omegagpi gHelicity=-1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 2 0 1 1 -1 90 .35 omegagpi -1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 2 0 1 1 -1 polFraction=0.35 polAngle=90 omegagpi gHelicity=-1 // If polarization vs E_gamma from histogram: -// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 90 pathToFile.root histName omegagpi -1 -// amplitude reactionName::SumName::WaveName Vec_ps_refl 1 1 2 -1 -1 polAngle=90 polFile=pathToFile.root polHist=histName omegagpi gHelicity=-1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 2 0 1 1 -1 90 pathToFile.root histName omegagpi -1 +// amplitude reactionName::SumName::WaveName Vec_ps_refl 2 0 1 1 -1 polAngle=90 polFile=pathToFile.root polHist=histName omegagpi gHelicity=-1 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -389,11 +388,11 @@ Vec_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ // beam delta ps pi0 gamma // Vec(KK) ps proton, (5 particles) // beam proton ps K K - // Vec(pipi) pi- Delta++, (5 particles) + // Vec(pipi) pi- Delta++, (5 particles) // beam proton ps pi pi - // Vec(Kpi) K+ Lambda, (5 particles) + // Vec(Kpi) K+ Lambda, (5 particles) // beam proton ps K pi - vecDaught1 = TLorentzVector(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); + vecDaught1 = TLorentzVector(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); vecDaught2 = TLorentzVector(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); vec = vecDaught1 + vecDaught2; }