diff --git a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc index f83c42af7..fbfc8720a 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,72 +82,41 @@ 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; - 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; - - 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 - 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); - } - } +void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const +{ + TLorentzVector beam( 0.0, 0.0, pKin[0][3], pKin[0][0]); + TVector3 eps(pKin[0][1], pKin[0][2], 0.); + // 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; @@ -158,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; + + 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 decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame + // 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]); @@ -175,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/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc index 956c633fc..d674cc690 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,15 +86,14 @@ 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 + // PlotGenerator works only with a 2-body vector decay - - // Fixed target + // Fixed target TLorentzVector target(0,0,0,0.938272); - + TVector3 eps; + //cout << "project event" << endl; TLorentzVector beam = kin->particle( 0 ); @@ -104,8 +102,8 @@ 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; @@ -120,25 +118,34 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ TLorentzVector proton_ps_b = proton + iso_daught1; TLorentzVector recoil_ps_b = recoil + iso_daught1; - + //Momentum transfer + double Mandt = fabs((target-recoil).M2()); - // 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 + + + // 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); - for(uint ioption=5; ioption 6) { + beam_polAngle = TMath::DegToRad()*parseValidatedNumber("polarization angle", args[6]); + 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); } - - - - //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); - + //forcibly put xy-components of the beam vector to zero + beam.SetX(0.); + beam.SetY(0.); + + + + + // 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..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 ); // 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 fbde50471..d4b59e8a7 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" @@ -310,6 +312,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() );