From 42edbfbc149ce5e15632e11117e288d490537316 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Thu, 24 Sep 2026 16:55:26 +0200 Subject: [PATCH 1/7] Amplitudes will work on trees with merged diamond orientations --- src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc | 80 ++++++++----------- src/libraries/AMPTOOLS_AMPS/decayAngles.cc | 6 +- .../AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc | 48 ++++++----- .../isops_plotter/isops_plotter.cc | 6 +- 4 files changed, 64 insertions(+), 76 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc index f83c42af7..6a931bc26 100644 --- a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc +++ b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc @@ -1,4 +1,3 @@ - #include #include #include @@ -44,9 +43,8 @@ static double parseValidatedNumber(const string& label, const string& argInput){ } -Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) : -UserAmplitude< Iso_ps_refl >( args ){ - +Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) : UserAmplitude< Iso_ps_refl >( args ) +{ // This function is only for the two-body Isobar decay: xi -> pipi // 6 possibilities to initialize this amplitude: @@ -84,60 +82,46 @@ UserAmplitude< Iso_ps_refl >( args ){ // Default polarization information stored in tree m_polInTree = true; - // Loop over any additional amplitude arguments to change defaults - for(uint ioption=6; ioptionGet(args[8].c_str()); - if(polFrac_vs_E != nullptr ){ - throw std::runtime_error( - "Iso_ps_refl ERROR: Could not find histogram '" + args[8] + - "' in file " + std::string(polOption.Data())); - } - } + // Polarization information passed with additional arguments + if (args.size() > 6){ + + m_polInTree = false; + polAngle = parseValidatedNumber("polarization angle", args[6]); + std::string pol_option = args[7]; + + if (pol_option.find(".root") == std::string::npos) + polFraction = parseValidatedNumber("polarization fraction", args[7]); else{ - polFraction = parseValidatedNumber("polarization fraction", args[7]); - } - } + polFraction = -1; + TFile* pol_file = new TFile(pol_option.c_str()); + polFrac_vs_E = (TH1D*)pol_file->Get(args[8].c_str()); + if(polFrac_vs_E == nullptr ) + throw std::runtime_error("Iso_ps_refl ERROR: Could not find histogram '" + args[8] + "' in file " + pol_option); + } } + } -void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{ - TLorentzVector beam; - TVector3 eps; - GDouble beam_polFraction; - GDouble beam_polAngle; +void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const +{ - 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; + GDouble beam_polFraction, beam_polAngle; + TLorentzVector beam( 0.0, 0.0, pKin[0][3], pKin[0][0]); - beam_polFraction = eps.Mag(); - beam_polAngle = eps.Phi(); + + if(m_polInTree){ + beam_polAngle = TMath::DegToRad()*pKin[0][1]; //Px + beam_polFraction = pKin[0][2]; //Py } else{ - beam.SetPxPyPzE( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); - beam_polAngle = polAngle; + beam_polAngle = TMath::DegToRad()*polAngle; + beam_polFraction = polFraction; - if(polFraction > 0.0){ // for fitting with fixed polarization - 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; - } else - beam_polFraction = polFrac_vs_E->GetBinContent(bin); + if(beam_polFraction == -1){ // extract the polarization fraction from its dependence on E_gamma + int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); + if (bin > 0 && bin <= polFrac_vs_E->GetXaxis()->GetNbins()) + beam_polFraction = polFrac_vs_E->GetBinContent(bin); } } diff --git a/src/libraries/AMPTOOLS_AMPS/decayAngles.cc b/src/libraries/AMPTOOLS_AMPS/decayAngles.cc index 7e0087ef7..a21f5f6fa 100644 --- a/src/libraries/AMPTOOLS_AMPS/decayAngles.cc +++ b/src/libraries/AMPTOOLS_AMPS/decayAngles.cc @@ -46,8 +46,8 @@ vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inv // Calculate production plane angle double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex) -{ - // whichFrame = 1 for helicity, 2 for GJ +{ + // whichFrame = 1 for helicity, 2 for GJ assert( whichFrame == 1 || whichFrame == 2 ); // Boost all P4 from lab to CM rest frame @@ -82,7 +82,7 @@ double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beam TVector3 y = locxyz[1]; TVector3 eps( cos( polAngle ), sin( polAngle ), 0.0 ); - + double phiProd = atan2( y.Dot( eps ), inverseCM.Vect().Unit().Dot( eps.Cross( y ) ) ); return phiProd; diff --git a/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc index 956c633fc..0266269b9 100644 --- a/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc +++ b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc @@ -62,11 +62,11 @@ PlotGenerator( ) void IsoPsPlotGenerator::createHistograms( ) { cout << " calls to bookHistogram go here" << endl; - bookHistogram( kProd_Ang, new Histogram1D( 50, -180., 180., "ProdAng", "Production Angle [deg.]" ) ); + bookHistogram( kProd_Ang, new Histogram1D( 50, -180., 180., "ProdAng", "#Phi [deg]" ) ); bookHistogram( kCosTheta, new Histogram1D( 50, -1., 1., "CosTheta_GJ", "cos#theta^{[GJ]}" ) ); - bookHistogram( kPhi, new Histogram1D( 50, -180., 180., "Phi_GJ", "#phi^{[GJ]} [deg.]" ) ); + bookHistogram( kPhi, new Histogram1D( 50, -180., 180., "Phi_GJ", "#phi^{[GJ]} [deg]" ) ); bookHistogram( kCosThetaH, new Histogram1D( 50, -1., 1., "CosTheta_HF", "cos#theta^{[HF]}" ) ); - bookHistogram( kPhiH, new Histogram1D( 50, -180., 180., "Phi_HF", "#phi^{[HF]} [deg.]" ) ); + bookHistogram( kPhiH, new Histogram1D( 50, -180., 180., "Phi_HF", "#phi^{[HF]} [deg]" ) ); bookHistogram( kIsoMass, new Histogram1D( 200, 0., 3., "MIso", "m(2#pi) [GeV]") ); bookHistogram( kIsoPsMass, new Histogram1D( 200, 0.2, 3.2, "MIsoPs", "m(3#pi) [GeV]") ); @@ -78,8 +78,7 @@ void IsoPsPlotGenerator::createHistograms( ) { } -void -IsoPsPlotGenerator::projectEvent( Kinematics* kin ){ +void IsoPsPlotGenerator::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 @@ -87,13 +86,12 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin ){ projectEvent( kin, "" ); } -void -IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ +void IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ - // We work only with a 2-body vector decay + // We work only with a 2-body vector decay - - // Fixed target + + // Fixed target TLorentzVector target(0,0,0,0.938272); @@ -104,8 +102,9 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ TLorentzVector iso_daught1 = kin->particle( 3 ); TLorentzVector iso_daught2 = kin->particle( 4 ); TLorentzVector piplusL = kin->particle( 5 ); - + + // Final state P4 momenta TLorentzVector X = iso_daught1 + iso_daught2 + bach; TLorentzVector recoil = proton + piplusL; @@ -121,24 +120,29 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ TLorentzVector recoil_ps_b = recoil + iso_daught1; - - // 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 Iso_ps_refl amplitude + + // Read polarization angle either from config file (if provided), or from the X-component of p_beam + // We assume here that the first amplitude in the list is always the Iso_ps_refl amplitude const vector< string > args = cfgInfo()->amplitudeList( reactionName, "", "" ).at(0)->factors().at(0); - for(uint ioption=5; ioption 6) + beam_polAngle = TMath::DegToRad()*parseValidatedNumber("polarization angle", args[6]); + else + beam_polAngle = TMath::DegToRad()*beam.X(); //Px should be the 0th position in TLorentzVector + + //forcibly put xy-components of the beam vector to zero + beam.SetX(0.); + beam.SetY(0.); + //Momentum transfer double Mandt = fabs((target-recoil).M2()); //Calculate production angle in the Gottfried-Jackson frame double prod_angle = TMath::RadToDeg()*getPhiProd(beam_polAngle, X, beam, target, 2, true); - + + // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame // Angles for the 1st permutation vector thetaPhiAnglesTwoStep_a = getTwoStepAngles(X, iso_a, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true); diff --git a/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc b/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc index 29ecb4a65..c836e20e5 100644 --- a/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc +++ b/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc @@ -128,8 +128,8 @@ int main( int argc, char* argv[] ){ atiSetup(); cout << "Plotgen results"<< endl; - // IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // optional can be omitted - IsoPsPlotGenerator plotGen( results ); + IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // 2nd argument is optional + // IsoPsPlotGenerator plotGen( results ); cout << " Initialized ati and PlotGen" << endl; @@ -409,7 +409,7 @@ int main( int argc, char* argv[] ){ for(unsigned int j = i+1; j < fullamps.size(); j++){ // leave only the Spring2017_PARA_0::ImagPosSign and Spring2017_PARA_0::ImagNegSign coherent sums - if (fullamps[i].find("Spring2017_PARA_0") == std::string::npos || fullamps[j].find("Spring2017_PARA_0") == std::string::npos) continue; + if (fullamps[i].find("Spring2017") == std::string::npos || fullamps[j].find("Spring2017") == std::string::npos) continue; if (fullamps[i].find("UniBG") != std::string::npos || fullamps[j].find("UniBG") != std::string::npos) continue; if (fullamps[i].find("Real") != std::string::npos || fullamps[j].find("Real") != std::string::npos) continue; From 0ce124c36eb23283586138026eaedd25e6846af6 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Fri, 25 Sep 2026 14:54:23 +0200 Subject: [PATCH 2/7] Register my amplituides in gen_amp_V2 --- src/programs/Simulation/gen_amp_V2/gen_amp_v2.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/programs/Simulation/gen_amp_V2/gen_amp_v2.cc b/src/programs/Simulation/gen_amp_V2/gen_amp_v2.cc index f5ba52fe7..82b9265d0 100644 --- a/src/programs/Simulation/gen_amp_V2/gen_amp_v2.cc +++ b/src/programs/Simulation/gen_amp_V2/gen_amp_v2.cc @@ -28,6 +28,8 @@ #include "AMPTOOLS_AMPS/Lambda1520Angles.h" #include "AMPTOOLS_AMPS/Lambda1520tdist.h" #include "AMPTOOLS_AMPS/Vec_ps_refl.h" +#include "AMPTOOLS_AMPS/Iso_ps_refl.h" +#include "AMPTOOLS_AMPS/PiPiSWaveAMPK.h" #include "AMPTOOLS_AMPS/Ylm.h" #include "AMPTOOLS_AMPS/Zlm.h" #include "AMPTOOLS_AMPS/DblRegge_FastEta.h" @@ -296,6 +298,8 @@ int main( int argc, char* argv[] ){ AmpToolsInterface::registerAmplitude( Lambda1520Angles() ); AmpToolsInterface::registerAmplitude( Lambda1520tdist() ); AmpToolsInterface::registerAmplitude( Vec_ps_refl() ); + AmpToolsInterface::registerAmplitude( Iso_ps_refl() ); + AmpToolsInterface::registerAmplitude( PiPiSWaveAMPK() ); AmpToolsInterface::registerAmplitude( Ylm() ); AmpToolsInterface::registerAmplitude( Zlm() ); AmpToolsInterface::registerAmplitude( Hist2D() ); From 9ff9ba2febe7783717b60aeadc3609cf2d95c028 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Mon, 28 Sep 2026 16:51:27 +0200 Subject: [PATCH 3/7] Restore decayAngles.cc --- src/libraries/AMPTOOLS_AMPS/decayAngles.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/decayAngles.cc b/src/libraries/AMPTOOLS_AMPS/decayAngles.cc index a21f5f6fa..7e0087ef7 100644 --- a/src/libraries/AMPTOOLS_AMPS/decayAngles.cc +++ b/src/libraries/AMPTOOLS_AMPS/decayAngles.cc @@ -46,8 +46,8 @@ vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inv // Calculate production plane angle double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex) -{ - // whichFrame = 1 for helicity, 2 for GJ +{ + // whichFrame = 1 for helicity, 2 for GJ assert( whichFrame == 1 || whichFrame == 2 ); // Boost all P4 from lab to CM rest frame @@ -82,7 +82,7 @@ double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beam TVector3 y = locxyz[1]; TVector3 eps( cos( polAngle ), sin( polAngle ), 0.0 ); - + double phiProd = atan2( y.Dot( eps ), inverseCM.Vect().Unit().Dot( eps.Cross( y ) ) ); return phiProd; From c3a87d09a0a76dd1e287ff2d5193d1a6755ab3e2 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Tue, 13 Oct 2026 17:13:32 +0200 Subject: [PATCH 4/7] Amplitude works on the trees with polarization info stored in the (Beam_Px, Beam_Py). Now, the format of epsilon-vector is implemented. --- src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc | 46 +++++++++---------- src/libraries/AMPTOOLS_AMPS/decayAngles.cc | 45 ++++++++++++++++++ src/libraries/AMPTOOLS_AMPS/decayAngles.h | 1 + .../isops_plotter/isops_plotter.cc | 4 +- .../Simulation/mcsmear/mcsmear_config.cc | 10 +++- 5 files changed, 80 insertions(+), 26 deletions(-) diff --git a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc index 6a931bc26..fbfc8720a 100644 --- a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc +++ b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc @@ -106,34 +106,17 @@ Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) : UserAmplitude< Iso_ps void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const { - GDouble beam_polFraction, beam_polAngle; TLorentzVector beam( 0.0, 0.0, pKin[0][3], pKin[0][0]); + TVector3 eps(pKin[0][1], pKin[0][2], 0.); - if(m_polInTree){ - beam_polAngle = TMath::DegToRad()*pKin[0][1]; //Px - beam_polFraction = pKin[0][2]; //Py - } - else{ - beam_polAngle = TMath::DegToRad()*polAngle; - beam_polFraction = polFraction; - - if(beam_polFraction == -1){ // extract the polarization fraction from its dependence on E_gamma - int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); - if (bin > 0 && bin <= polFrac_vs_E->GetXaxis()->GetNbins()) - beam_polFraction = polFrac_vs_E->GetBinContent(bin); - } - } - - // Fill in four-vectors for final state particles TLorentzVector target(0,0,0,0.938272); // Fixed target TLorentzVector recoil( pKin[1][1]+pKin[5][1], pKin[1][2]+pKin[5][2], pKin[1][3]+pKin[5][3], pKin[1][0]+pKin[5][0] ); // Recoiling baryon TLorentzVector ps(pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0]); // 1st after proton is always the pseudoscalar meson - // Compute isobar meson from its decay products (xi -> pipi) - // Make sure the order of daughters is correct in the config file! + // Fill in four-vectors for the isobar meson (xi -> pipi, make sure the order of daughters is correct in the config file!) TLorentzVector iso_daught1(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); TLorentzVector iso_daught2(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]); TLorentzVector iso = iso_daught1 + iso_daught2; @@ -142,13 +125,29 @@ void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const TLorentzVector X = iso + ps; - //Calculate production angle in the Gottfried-Jackson frame - GDouble prod_angle = getPhiProd(beam_polAngle, X, beam, target, 2, true); + //Define Pgamma and calculate BigPhi angle in the Gottfried-Jackson frame + GDouble beam_polFraction, beam_polAngle, prod_angle; + + if(m_polInTree){ + beam_polFraction = eps.Mag(); + prod_angle = getPhiProd(eps, X, beam, target, 2, true); + } + else{ + beam_polAngle = TMath::DegToRad()*polAngle; //conversion factor, since angles in the config files are by default in degrees + prod_angle = getPhiProd(beam_polAngle, X, beam, target, 2, true); + beam_polFraction = polFraction; - // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame + if(beam_polFraction == -1){ // extract the polarization fraction from its dependence on E_gamma + int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); + if (bin > 0 && bin <= polFrac_vs_E->GetXaxis()->GetNbins()) + beam_polFraction = polFrac_vs_E->GetBinContent(bin); + } + } + + + // Calculate two angles of X decay (we use the Gottfried-Jackson frame) and two angles of Isobar decay (we use the Helicity frame) vector thetaPhiAnglesTwoStep; thetaPhiAnglesTwoStep = getTwoStepAngles(X, iso, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true); - GDouble cosTheta = TMath::Cos(thetaPhiAnglesTwoStep[0]); @@ -159,6 +158,7 @@ void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const GDouble m_iso = iso.M(); GDouble m_ps = ps.M(); + // Calculate amplitude complex amplitude(0,0); complex i(0,1); diff --git a/src/libraries/AMPTOOLS_AMPS/decayAngles.cc b/src/libraries/AMPTOOLS_AMPS/decayAngles.cc index 7e0087ef7..733287d6e 100644 --- a/src/libraries/AMPTOOLS_AMPS/decayAngles.cc +++ b/src/libraries/AMPTOOLS_AMPS/decayAngles.cc @@ -44,6 +44,7 @@ vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inv return xyz; } + // Calculate production plane angle double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex) { @@ -89,6 +90,50 @@ double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beam } +// Calculate production plane angle using eps-vector already defined in user's amplitude +double getPhiProd(TVector3 eps, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex) +{ + // whichFrame = 1 for helicity, 2 for GJ + assert( whichFrame == 1 || whichFrame == 2 ); + + // Boost all P4 from lab to CM rest frame + TLorentzVector cmLab = beamLab + targetLab; + TVector3 boostCMToLab = cmLab.BoostVector(); + + TLorentzVector parentCM = parentLab; + TLorentzVector inverseLab; + if( upperVertex == false ) + inverseLab = targetLab; + else + inverseLab = beamLab; + TLorentzVector inverseCM = inverseLab; + parentCM.Boost( -1.0*boostCMToLab ); + inverseCM.Boost( -1.0*boostCMToLab ); + + TVector3 boostParentToCM = parentCM.BoostVector(); + TLorentzVector inverseParent = inverseCM; + inverseParent.Boost( -1.0*boostParentToCM ); + + // convert P4s to P3 + TVector3 parentCM3 = parentCM.Vect(); + TVector3 inverseCM3 = inverseCM.Vect(); + TVector3 inverseParent3 = inverseParent.Vect(); + + vector< TVector3 > locxyz; + if( whichFrame == 1 ) + locxyz = getHelicityAxes(parentCM3, inverseCM3); + else + locxyz = getGJAxes(parentCM3, inverseCM3, inverseParent3); + + TVector3 y = locxyz[1]; + + double phiProd = atan2( y.Dot( eps ), inverseCM.Vect().Unit().Dot( eps.Cross( y ) ) ); + + return phiProd; +} + + + // Calculate angles for a one-step decay. If only one of the daughter particles has spin, that one should be used in this calculation vector< double > getOneStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex) { diff --git a/src/libraries/AMPTOOLS_AMPS/decayAngles.h b/src/libraries/AMPTOOLS_AMPS/decayAngles.h index ced397cc1..268b4ca8b 100644 --- a/src/libraries/AMPTOOLS_AMPS/decayAngles.h +++ b/src/libraries/AMPTOOLS_AMPS/decayAngles.h @@ -16,6 +16,7 @@ vector< TVector3 > getHelicityAxes(TVector3 parentCM, TVector3 inverseCM); vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inverseParent); double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex); +double getPhiProd(TVector3 eps, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex); vector< double > getOneStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex = true); diff --git a/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc b/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc index c836e20e5..03f52b5ab 100644 --- a/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc +++ b/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc @@ -128,8 +128,8 @@ int main( int argc, char* argv[] ){ atiSetup(); cout << "Plotgen results"<< endl; - IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // 2nd argument is optional - // IsoPsPlotGenerator plotGen( results ); + // IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // 2nd argument is optional + IsoPsPlotGenerator plotGen( results ); cout << " Initialized ati and PlotGen" << endl; diff --git a/src/programs/Simulation/mcsmear/mcsmear_config.cc b/src/programs/Simulation/mcsmear/mcsmear_config.cc index ba03889ec..a92d6c052 100644 --- a/src/programs/Simulation/mcsmear/mcsmear_config.cc +++ b/src/programs/Simulation/mcsmear/mcsmear_config.cc @@ -134,8 +134,16 @@ bool mcsmear_config_t::ParseRCDBConfigFile(int runNumber) } auto json = rtvsCondition->ToJsonDocument(); // The CODA rtvs is serialized as JSon dictionary. - string fileName(json.at("%(config)").get_string()); // The file name is stored in '%(config)' + if(!json.IsObject() || !json.HasMember("%(config)") || !json["%(config)"].IsString()) { + jerr << "RCDB: 'rtvs' condition for run " << runNumber + << " has no string member '%(config)'" << endl; + return false; + } + + string fileName(json["%(config)"].GetString()); // The file name is stored in '%(config)' + + // Get file out of RCDB (indexed by run number and name) auto file = rcdb_connection->GetFile(runNumber, fileName); From 7c48d476995449162ebc470de8109ec890c003f6 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Wed, 14 Oct 2026 11:30:14 +0200 Subject: [PATCH 5/7] Plot generator adjusted accordingly --- .../AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc | 31 ++++++++++--------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc index 0266269b9..d674cc690 100644 --- a/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc +++ b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc @@ -88,12 +88,12 @@ void IsoPsPlotGenerator::projectEvent( Kinematics* kin ){ void IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ - // We work only with a 2-body vector decay - + // PlotGenerator works only with a 2-body vector decay // Fixed target TLorentzVector target(0,0,0,0.938272); - + TVector3 eps; + //cout << "project event" << endl; TLorentzVector beam = kin->particle( 0 ); @@ -103,7 +103,6 @@ void IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionNa TLorentzVector iso_daught2 = kin->particle( 4 ); TLorentzVector piplusL = kin->particle( 5 ); - // Final state P4 momenta TLorentzVector X = iso_daught1 + iso_daught2 + bach; @@ -119,28 +118,32 @@ void IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionNa TLorentzVector proton_ps_b = proton + iso_daught1; TLorentzVector recoil_ps_b = recoil + iso_daught1; + //Momentum transfer + double Mandt = fabs((target-recoil).M2()); + - // Read polarization angle either from config file (if provided), or from the X-component of p_beam + // Either read the diamond orientation angle directly from the config file (if provided), or read eps-vector from the (Beam_PX,Beam_PY)-components in the input tree + // Then calculate polarization angle in the Gottfried-Jackson frame // We assume here that the first amplitude in the list is always the Iso_ps_refl amplitude const vector< string > args = cfgInfo()->amplitudeList( reactionName, "", "" ).at(0)->factors().at(0); - double beam_polAngle; + double beam_polAngle, prod_angle; - if (args.size() > 6) + if (args.size() > 6) { beam_polAngle = TMath::DegToRad()*parseValidatedNumber("polarization angle", args[6]); - else - beam_polAngle = TMath::DegToRad()*beam.X(); //Px should be the 0th position in TLorentzVector + prod_angle = TMath::RadToDeg()*getPhiProd(beam_polAngle, X, beam, target, 2, true); + } + else { + eps.SetXYZ(beam.X(),beam.Y(),0.); + prod_angle = TMath::RadToDeg()*getPhiProd(eps, X, beam, target, 2, true); + } //forcibly put xy-components of the beam vector to zero beam.SetX(0.); beam.SetY(0.); - - //Momentum transfer - double Mandt = fabs((target-recoil).M2()); - //Calculate production angle in the Gottfried-Jackson frame - double prod_angle = TMath::RadToDeg()*getPhiProd(beam_polAngle, X, beam, target, 2, true); + // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame From ea49fd39b4afb29f70cb439445c074340745ff08 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Wed, 14 Oct 2026 11:55:13 +0200 Subject: [PATCH 6/7] Exclude changes in the mcsmear_config.cc file --- .../Simulation/mcsmear/mcsmear_config.cc | 300 ------------------ 1 file changed, 300 deletions(-) delete mode 100644 src/programs/Simulation/mcsmear/mcsmear_config.cc diff --git a/src/programs/Simulation/mcsmear/mcsmear_config.cc b/src/programs/Simulation/mcsmear/mcsmear_config.cc deleted file mode 100644 index a92d6c052..000000000 --- a/src/programs/Simulation/mcsmear/mcsmear_config.cc +++ /dev/null @@ -1,300 +0,0 @@ - -#include "mcsmear_config.h" - - -#include -#include - -#ifdef HAVE_RCDB -#include -#include "RCDB/ConfigParser.h" -string RCDB_CONNECTION; -rcdb::Connection *rcdb_connection; -#endif // HAVE_RCDB - - -//----------- -// mcsmear_config_t (Constructor) -//----------- -mcsmear_config_t::mcsmear_config_t() -{ - // default values - DROP_TRUTH_HITS = false; - ADD_NOISE = false; - SMEAR_HITS = true; - //SMEAR_BCAL = true; - IGNORE_SEEDS = false; - DUMP_RCDB_CONFIG = false; - APPLY_EFFICIENCY_CORRECTIONS = true; - APPLY_HITS_TRUNCATION = true; - - FCAL_ADD_LIGHTGUIDE_HITS = true; - FCAL_LIGHTGUIDE_SCALE_FACTOR = -1.; - FCAL_NEW_TIME_SMEAR = false; - - SKIP_READING_RCDB = false; - MERGE_TAGGER_HITS = true; - - BCAL_NO_T_SMEAR = false; - BCAL_NO_DARK_PULSES = false; - BCAL_NO_SAMPLING_FLUCTUATIONS = false; - BCAL_NO_SAMPLING_FLOOR_TERM = false; - BCAL_NO_POISSON_STATISTICS = false; - BCAL_NO_FADC_SATURATION = false; - BCAL_NO_SIPM_SATURATION = false; - - TRIGGER_LOOKBACK_TIME = -100; // ns - -#ifdef HAVE_RCDB - // RCDB configuration - // first determine which database to connect to - if( getenv("RCDB_CONNECTION")!= NULL ) - RCDB_CONNECTION = getenv("RCDB_CONNECTION"); - else - RCDB_CONNECTION = "mysql://rcdb@hallddb.jlab.org/rcdb"; // default to outward-facing MySQL DB - - // load RCDB later, so that the DApplication interface is initialized - rcdb_connection = NULL; -#endif //HAVE_RCDB -} - -//----------- -// mcsmear_config_t (Destructor) -//----------- -mcsmear_config_t::~mcsmear_config_t() { -#ifdef HAVE_RCDB - if (rcdb_connection != NULL) - delete rcdb_connection; -#endif //HAVE_RCDB -} - -//----------- -// SetSeeds -//----------- -void mcsmear_config_t::SetSeeds(const char *vals) -{ - /// This is called from the command line parser to - /// set the initial seeds based on user input from - /// the command line. - // - std::stringstream ss(vals); - Int_t seed1, seed2, seed3; - ss >> seed1 >> seed2 >> seed3; - UInt_t *useed1 = reinterpret_cast(&seed1); - UInt_t *useed2 = reinterpret_cast(&seed2); - UInt_t *useed3 = reinterpret_cast(&seed3); - gDRandom.SetSeeds(*useed1, *useed2, *useed3); - - cout << "Seeds set from command line. Any random number" << endl; - cout << "seeds found in the input file will be ignored!" << endl; - IGNORE_SEEDS = true; -} - - -#ifdef HAVE_RCDB - -//----------- -// LoadRCDBConnection -//----------- -void mcsmear_config_t::LoadRCDBConnection() -{ - // We want to connect to RCDB as late as possible for two reasons: - // 1) No need to connect to the database unless we are actually using this information, - // which should speed things up - // 2) To use the JANA command line parameter interface, we need to make sure that the - // DApplication is initialized, and it might not be when the mcsmear_config_t - // constructor is called - - // if we're already connected, then stop now - if(rcdb_connection != NULL) - return; - japp->SetDefaultParameter("RCDB_CONNECTION", RCDB_CONNECTION, "URL used to access RCDB."); - - // load connection to RCDB - rcdb_connection = new rcdb::Connection(RCDB_CONNECTION); -} - -//----------- -// ParseRCDBConfigFile -//----------- -bool mcsmear_config_t::ParseRCDBConfigFile(int runNumber) -{ - // To get a lot of the configuration parameters we need, we need to parse the CODA configuration files - // which are saved in RCDB in the JSON format - - // Lazily connect to RCDB - LoadRCDBConnection(); - - // The "rtvs" condition contains the file name of the CODA configuration file - // What else does this contain?? - auto rtvsCondition = rcdb_connection->GetCondition(runNumber, "rtvs"); // Get condition by run and name - if(!rtvsCondition) { - jerr << "RCDB: 'rtvs' condition is not set for run " << runNumber << endl; - return false; - } - - auto json = rtvsCondition->ToJsonDocument(); // The CODA rtvs is serialized as JSon dictionary. - - if(!json.IsObject() || !json.HasMember("%(config)") || !json["%(config)"].IsString()) { - jerr << "RCDB: 'rtvs' condition for run " << runNumber - << " has no string member '%(config)'" << endl; - return false; - } - - string fileName(json["%(config)"].GetString()); // The file name is stored in '%(config)' - - - - // Get file out of RCDB (indexed by run number and name) - auto file = rcdb_connection->GetFile(runNumber, fileName); - if(!file) { // If there is no such file, null is returned - jerr << "File with name \"" << fileName - << "\" doesn't exist for run "<< runNumber << endl; - return false; - } - - // Parse CODA config file - vector SectionNames = {"TRIGGER", "GLOBAL", "FCAL", "ECAL", "CCAL", "BCAL", "TOF", "ST", "TAGH", - "TAGM", "PS", "PSC", "TPOL", "CDC", "FDC"}; - string fileContent = file->GetContent(); // Get file content - auto result = rcdb::ConfigParser::Parse(fileContent, SectionNames); // Parse it! - - // EXAMPLE -// double CDC_FADC125_DAC = stod(result.Sections["CDC"].NameValues["FADC125_DAC"]); -// double CDC_FADC125_THR = stod(result.Sections["CDC"].NameValues["FADC125_THR"]); - - if(DUMP_RCDB_CONFIG) { - //// DEBUG //// - //cout << "CODA config file contents:" << endl; - std::ofstream coda_ofile("rcdb_coda.config"); - coda_ofile << "Full file: " << endl; - coda_ofile << fileContent << endl; - - - coda_ofile << endl << "Parsed:" << endl; - for(auto sectionName: SectionNames) { - coda_ofile << "Section: " << sectionName << endl; - auto sectionData = result.Sections[sectionName]; - - for(auto value : sectionData.NameValues) - coda_ofile << value.first << " = " << value.second << endl; - for(auto data_vec : sectionData.NameVectors) { - coda_ofile << data_vec.first << " = "; - for(auto value: data_vec.second) - coda_ofile << value << " "; - coda_ofile << endl; - } - - coda_ofile << endl; - } - coda_ofile.close(); - } - - // then we get stuff out of it - - double dvalue; - std::stringstream deco; - vector::iterator iter; - vector fadc250_sys = {"FCAL", "ECAL", "CCAL", "BCAL", "TOF", "ST", "TAGH", - "TAGM", "PS", "PSC", "TPOL"}; - for (iter = fadc250_sys.begin(); iter != fadc250_sys.end(); ++iter) { - - // Make sure that the sub-detector exists in the configuration file - auto section_length = result.Sections[*iter].Rows; - if(section_length.size() == 0) continue; - - - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC250_NPEAK"]); - deco >> dvalue; - readout[*iter]["NPEAK"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC250_NSA"]); - deco >> dvalue; - readout[*iter]["NSA"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC250_NSB"]); - deco >> dvalue; - readout[*iter]["NSB"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC250_W_WIDTH"]); - deco >> dvalue; - readout[*iter]["WINDOW"] = dvalue; - } - - vector fadc125_sys = {"FDC", "CDC"}; - for (iter = fadc125_sys.begin(); iter != fadc125_sys.end(); ++iter) { - - auto section_length = result.Sections[*iter].Rows; - if(section_length.size() == 0) continue; - - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC125_NPEAK"]); - deco >> dvalue; - readout[*iter]["NPEAK"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC125_IE"]); - deco >> dvalue; - readout[*iter]["IE"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC125_PG"]); - deco >> dvalue; - readout[*iter]["PG"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["FADC125_W_WIDTH"]); - deco >> dvalue; - readout[*iter]["NW"] = dvalue; - } - - vector f1tdc_sys = {"BCAL", "ST", "TAGH", "TAGM", "PSC", "FDC"}; - for (iter = f1tdc_sys.begin(); iter != f1tdc_sys.end(); ++iter) { - - auto section_length = result.Sections[*iter].Rows; - if(section_length.size() == 0) continue; - - readout[*iter]["NHITS"] = 8.; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["F1TDC_WINDOW"]); - deco >> dvalue; - readout[*iter]["WINDOW"] = dvalue; - deco.clear(); - if (result.Sections.at(*iter).NameVectors["DSC2_WIDTH"].size() > 0) { - deco.str(result.Sections.at(*iter).NameVectors["DSC2_WIDTH"][0]); - deco >> dvalue; - } - else { - dvalue = 50; // default value; - } - readout[*iter]["WIDTH"] = dvalue; - } - - vector tdc1290_sys = {"TOF"}; - for (iter = tdc1290_sys.begin(); iter != tdc1290_sys.end(); ++iter) { - - auto section_length = result.Sections[*iter].Rows; - if(section_length.size() == 0) continue; - - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["TDC1290_N_HITS"]); - deco >> dvalue; - readout[*iter]["NHITS"] = dvalue; - deco.clear(); - deco.str(result.Sections.at(*iter).NameValues["TDC1290_W_WIDTH"]); - deco >> dvalue; - readout[*iter]["WINDOW"] = dvalue; - deco.clear(); - if (result.Sections.at(*iter).NameVectors["DSC2_WIDTH"].size() > 0) { - deco.str(result.Sections.at(*iter).NameVectors["DSC2_WIDTH"][0]); - deco >> dvalue; - } - else { - dvalue = 50; // default value - } - readout[*iter]["WIDTH"] = dvalue; - } - - return true; -} - -#endif // HAVE_RCDB - From 1651329da7756cfad6bdec61750fb5e3dc52e5b9 Mon Sep 17 00:00:00 2001 From: Ilia Belov Date: Wed, 14 Oct 2026 12:02:42 +0200 Subject: [PATCH 7/7] Retain mcsmear_config.cc from master --- .../Simulation/mcsmear/mcsmear_config.cc | 292 ++++++++++++++++++ 1 file changed, 292 insertions(+) create mode 100644 src/programs/Simulation/mcsmear/mcsmear_config.cc diff --git a/src/programs/Simulation/mcsmear/mcsmear_config.cc b/src/programs/Simulation/mcsmear/mcsmear_config.cc new file mode 100644 index 000000000..ba03889ec --- /dev/null +++ b/src/programs/Simulation/mcsmear/mcsmear_config.cc @@ -0,0 +1,292 @@ + +#include "mcsmear_config.h" + + +#include +#include + +#ifdef HAVE_RCDB +#include +#include "RCDB/ConfigParser.h" +string RCDB_CONNECTION; +rcdb::Connection *rcdb_connection; +#endif // HAVE_RCDB + + +//----------- +// mcsmear_config_t (Constructor) +//----------- +mcsmear_config_t::mcsmear_config_t() +{ + // default values + DROP_TRUTH_HITS = false; + ADD_NOISE = false; + SMEAR_HITS = true; + //SMEAR_BCAL = true; + IGNORE_SEEDS = false; + DUMP_RCDB_CONFIG = false; + APPLY_EFFICIENCY_CORRECTIONS = true; + APPLY_HITS_TRUNCATION = true; + + FCAL_ADD_LIGHTGUIDE_HITS = true; + FCAL_LIGHTGUIDE_SCALE_FACTOR = -1.; + FCAL_NEW_TIME_SMEAR = false; + + SKIP_READING_RCDB = false; + MERGE_TAGGER_HITS = true; + + BCAL_NO_T_SMEAR = false; + BCAL_NO_DARK_PULSES = false; + BCAL_NO_SAMPLING_FLUCTUATIONS = false; + BCAL_NO_SAMPLING_FLOOR_TERM = false; + BCAL_NO_POISSON_STATISTICS = false; + BCAL_NO_FADC_SATURATION = false; + BCAL_NO_SIPM_SATURATION = false; + + TRIGGER_LOOKBACK_TIME = -100; // ns + +#ifdef HAVE_RCDB + // RCDB configuration + // first determine which database to connect to + if( getenv("RCDB_CONNECTION")!= NULL ) + RCDB_CONNECTION = getenv("RCDB_CONNECTION"); + else + RCDB_CONNECTION = "mysql://rcdb@hallddb.jlab.org/rcdb"; // default to outward-facing MySQL DB + + // load RCDB later, so that the DApplication interface is initialized + rcdb_connection = NULL; +#endif //HAVE_RCDB +} + +//----------- +// mcsmear_config_t (Destructor) +//----------- +mcsmear_config_t::~mcsmear_config_t() { +#ifdef HAVE_RCDB + if (rcdb_connection != NULL) + delete rcdb_connection; +#endif //HAVE_RCDB +} + +//----------- +// SetSeeds +//----------- +void mcsmear_config_t::SetSeeds(const char *vals) +{ + /// This is called from the command line parser to + /// set the initial seeds based on user input from + /// the command line. + // + std::stringstream ss(vals); + Int_t seed1, seed2, seed3; + ss >> seed1 >> seed2 >> seed3; + UInt_t *useed1 = reinterpret_cast(&seed1); + UInt_t *useed2 = reinterpret_cast(&seed2); + UInt_t *useed3 = reinterpret_cast(&seed3); + gDRandom.SetSeeds(*useed1, *useed2, *useed3); + + cout << "Seeds set from command line. Any random number" << endl; + cout << "seeds found in the input file will be ignored!" << endl; + IGNORE_SEEDS = true; +} + + +#ifdef HAVE_RCDB + +//----------- +// LoadRCDBConnection +//----------- +void mcsmear_config_t::LoadRCDBConnection() +{ + // We want to connect to RCDB as late as possible for two reasons: + // 1) No need to connect to the database unless we are actually using this information, + // which should speed things up + // 2) To use the JANA command line parameter interface, we need to make sure that the + // DApplication is initialized, and it might not be when the mcsmear_config_t + // constructor is called + + // if we're already connected, then stop now + if(rcdb_connection != NULL) + return; + japp->SetDefaultParameter("RCDB_CONNECTION", RCDB_CONNECTION, "URL used to access RCDB."); + + // load connection to RCDB + rcdb_connection = new rcdb::Connection(RCDB_CONNECTION); +} + +//----------- +// ParseRCDBConfigFile +//----------- +bool mcsmear_config_t::ParseRCDBConfigFile(int runNumber) +{ + // To get a lot of the configuration parameters we need, we need to parse the CODA configuration files + // which are saved in RCDB in the JSON format + + // Lazily connect to RCDB + LoadRCDBConnection(); + + // The "rtvs" condition contains the file name of the CODA configuration file + // What else does this contain?? + auto rtvsCondition = rcdb_connection->GetCondition(runNumber, "rtvs"); // Get condition by run and name + if(!rtvsCondition) { + jerr << "RCDB: 'rtvs' condition is not set for run " << runNumber << endl; + return false; + } + + auto json = rtvsCondition->ToJsonDocument(); // The CODA rtvs is serialized as JSon dictionary. + string fileName(json.at("%(config)").get_string()); // The file name is stored in '%(config)' + + + // Get file out of RCDB (indexed by run number and name) + auto file = rcdb_connection->GetFile(runNumber, fileName); + if(!file) { // If there is no such file, null is returned + jerr << "File with name \"" << fileName + << "\" doesn't exist for run "<< runNumber << endl; + return false; + } + + // Parse CODA config file + vector SectionNames = {"TRIGGER", "GLOBAL", "FCAL", "ECAL", "CCAL", "BCAL", "TOF", "ST", "TAGH", + "TAGM", "PS", "PSC", "TPOL", "CDC", "FDC"}; + string fileContent = file->GetContent(); // Get file content + auto result = rcdb::ConfigParser::Parse(fileContent, SectionNames); // Parse it! + + // EXAMPLE +// double CDC_FADC125_DAC = stod(result.Sections["CDC"].NameValues["FADC125_DAC"]); +// double CDC_FADC125_THR = stod(result.Sections["CDC"].NameValues["FADC125_THR"]); + + if(DUMP_RCDB_CONFIG) { + //// DEBUG //// + //cout << "CODA config file contents:" << endl; + std::ofstream coda_ofile("rcdb_coda.config"); + coda_ofile << "Full file: " << endl; + coda_ofile << fileContent << endl; + + + coda_ofile << endl << "Parsed:" << endl; + for(auto sectionName: SectionNames) { + coda_ofile << "Section: " << sectionName << endl; + auto sectionData = result.Sections[sectionName]; + + for(auto value : sectionData.NameValues) + coda_ofile << value.first << " = " << value.second << endl; + for(auto data_vec : sectionData.NameVectors) { + coda_ofile << data_vec.first << " = "; + for(auto value: data_vec.second) + coda_ofile << value << " "; + coda_ofile << endl; + } + + coda_ofile << endl; + } + coda_ofile.close(); + } + + // then we get stuff out of it + + double dvalue; + std::stringstream deco; + vector::iterator iter; + vector fadc250_sys = {"FCAL", "ECAL", "CCAL", "BCAL", "TOF", "ST", "TAGH", + "TAGM", "PS", "PSC", "TPOL"}; + for (iter = fadc250_sys.begin(); iter != fadc250_sys.end(); ++iter) { + + // Make sure that the sub-detector exists in the configuration file + auto section_length = result.Sections[*iter].Rows; + if(section_length.size() == 0) continue; + + + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC250_NPEAK"]); + deco >> dvalue; + readout[*iter]["NPEAK"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC250_NSA"]); + deco >> dvalue; + readout[*iter]["NSA"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC250_NSB"]); + deco >> dvalue; + readout[*iter]["NSB"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC250_W_WIDTH"]); + deco >> dvalue; + readout[*iter]["WINDOW"] = dvalue; + } + + vector fadc125_sys = {"FDC", "CDC"}; + for (iter = fadc125_sys.begin(); iter != fadc125_sys.end(); ++iter) { + + auto section_length = result.Sections[*iter].Rows; + if(section_length.size() == 0) continue; + + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC125_NPEAK"]); + deco >> dvalue; + readout[*iter]["NPEAK"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC125_IE"]); + deco >> dvalue; + readout[*iter]["IE"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC125_PG"]); + deco >> dvalue; + readout[*iter]["PG"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["FADC125_W_WIDTH"]); + deco >> dvalue; + readout[*iter]["NW"] = dvalue; + } + + vector f1tdc_sys = {"BCAL", "ST", "TAGH", "TAGM", "PSC", "FDC"}; + for (iter = f1tdc_sys.begin(); iter != f1tdc_sys.end(); ++iter) { + + auto section_length = result.Sections[*iter].Rows; + if(section_length.size() == 0) continue; + + readout[*iter]["NHITS"] = 8.; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["F1TDC_WINDOW"]); + deco >> dvalue; + readout[*iter]["WINDOW"] = dvalue; + deco.clear(); + if (result.Sections.at(*iter).NameVectors["DSC2_WIDTH"].size() > 0) { + deco.str(result.Sections.at(*iter).NameVectors["DSC2_WIDTH"][0]); + deco >> dvalue; + } + else { + dvalue = 50; // default value; + } + readout[*iter]["WIDTH"] = dvalue; + } + + vector tdc1290_sys = {"TOF"}; + for (iter = tdc1290_sys.begin(); iter != tdc1290_sys.end(); ++iter) { + + auto section_length = result.Sections[*iter].Rows; + if(section_length.size() == 0) continue; + + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["TDC1290_N_HITS"]); + deco >> dvalue; + readout[*iter]["NHITS"] = dvalue; + deco.clear(); + deco.str(result.Sections.at(*iter).NameValues["TDC1290_W_WIDTH"]); + deco >> dvalue; + readout[*iter]["WINDOW"] = dvalue; + deco.clear(); + if (result.Sections.at(*iter).NameVectors["DSC2_WIDTH"].size() > 0) { + deco.str(result.Sections.at(*iter).NameVectors["DSC2_WIDTH"][0]); + deco >> dvalue; + } + else { + dvalue = 50; // default value + } + readout[*iter]["WIDTH"] = dvalue; + } + + return true; +} + +#endif // HAVE_RCDB +