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_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; 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 163d3eff6..e51b4f410 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" @@ -308,6 +310,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() );