diff --git a/src/libraries/AMPTOOLS_AMPS/GPUTwoPiAngles_Delta_DoubleSDMEs_kernel.cu b/src/libraries/AMPTOOLS_AMPS/GPUTwoPiAngles_Delta_DoubleSDMEs_kernel.cu new file mode 100755 index 000000000..eaff81ec0 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/GPUTwoPiAngles_Delta_DoubleSDMEs_kernel.cu @@ -0,0 +1,208 @@ + +#include + +#include "GPUManager/GPUCustomTypes.h" +#include "GPUManager/CUDA-Complex.cuh" + + +__global__ void +// +GPUTwoPiAngles_Delta_DoubleSDMEs_kernel( GPU_AMP_PROTO, + //order is like the appearing order in the paper + //Wmm + GDouble r00_33_0, + GDouble r11_33_0, GDouble r11_11_0, + GDouble r1m1_33_0, GDouble r1m1_11_0, + GDouble r10_33_0, GDouble r10_11_0, + + //Wbar + GDouble r00_31_0, GDouble r00_3m1_0, + GDouble r11_31_0, GDouble r11_3m1_0, + GDouble r1m1_31_0, GDouble r1m1_3m1_0, + GDouble r10_31_0, GDouble r10_3m1_0, + + //Wtil + GDouble rt1m1_31_0, GDouble rt1m1_3m1_0, + GDouble rt10_31_0, GDouble rt10_3m1_0, + // //alpha = 1 (same structure) + GDouble r00_33_1, GDouble r00_11_1, + GDouble r11_33_1, GDouble r11_11_1, + GDouble r1m1_33_1, GDouble r1m1_11_1, + GDouble r10_33_1, GDouble r10_11_1, + GDouble r00_31_1, GDouble r00_3m1_1, + GDouble r11_31_1, GDouble r11_3m1_1, + GDouble r1m1_31_1, GDouble r1m1_3m1_1, + GDouble r10_31_1, GDouble r10_3m1_1, + GDouble rt1m1_31_1, GDouble rt1m1_3m1_1, + GDouble rt10_31_1, GDouble rt10_3m1_1, + //alpha = 2 (different structure) + GDouble rt00_31_2, GDouble rt00_3m1_2, + GDouble rt11_31_2, GDouble rt11_3m1_2, + GDouble r1m1_33_2, GDouble r1m1_11_2, + GDouble r1m1_31_2, GDouble r1m1_3m1_2, + GDouble rt1m1_31_2, GDouble rt1m1_3m1_2, + GDouble r10_33_2, GDouble r10_11_2, + GDouble r10_31_2, GDouble r10_3m1_2, + GDouble rt10_31_2, GDouble rt10_3m1_2 +) + { + + int iEvent = GPU_THIS_EVENT; + + // here we need to be careful to index the user-defined + // data with the proper integer corresponding to the + // enumeration in the C++ header file + + + GDouble Pgamma = GPU_UVARS(0); + GDouble cosTheta_pim = GPU_UVARS(1); + GDouble sinSqTheta_pim = GPU_UVARS(2); + GDouble sin2Theta_pim = GPU_UVARS(3); + GDouble phi_pim = GPU_UVARS(4); + GDouble cosTheta_proton = GPU_UVARS(5); + GDouble sinSqTheta_proton = GPU_UVARS(6); + GDouble sin2Theta_proton = GPU_UVARS(7); + GDouble phi_proton = GPU_UVARS(8); + GDouble bigPhi = GPU_UVARS(9); + // flag for pol (0.0 = no pol, 1.0 = pol) + + GDouble sinSqTh = sinSqTheta_proton; + GDouble sin2Th = sin2Theta_proton; + GDouble cosSqTh = cosTheta_proton * cosTheta_proton; + + GDouble cphi = cos(phi_proton); + GDouble sphi = sin(phi_proton); + GDouble c2phi = cos(2.0 * phi_proton); + GDouble s2phi = sin(2.0 * phi_proton); + + GDouble sinSqThPi = sinSqTheta_pim; + GDouble sin2ThPi = sin2Theta_pim; + GDouble cosSqThPi = cosTheta_pim * cosTheta_pim; + + GDouble cphiPi = cos(phi_pim); + GDouble sphiPi = sin(phi_pim); + GDouble c2phiPi = cos(2.0 * phi_pim); + GDouble s2phiPi = sin(2.0 * phi_pim); + + GDouble r00_11_0 = 0.5 - r00_33_0 - 2.0 * r11_11_0 - 2.0 * r11_33_0; + + GDouble W00_0 = r00_33_0 * sinSqTh + r00_11_0 * (1.0/3.0 + cosSqTh); + GDouble Wb00_0 = (2.0 / sqrt(3.0)) * (r00_31_0 * cphi * sin2Th + r00_3m1_0 * c2phi * sinSqTh); + + GDouble W11_0 = r11_33_0 * sinSqTh + r11_11_0 * (1.0/3.0 + cosSqTh); + GDouble Wb11_0 = (2.0 / sqrt(3.0)) * (r11_31_0 * cphi * sin2Th + r11_3m1_0 * c2phi * sinSqTh); + + GDouble W1m1_0 = r1m1_33_0 * sinSqTh + r1m1_11_0 * (1.0/3.0 + cosSqTh); + GDouble Wb1m1_0 = (2.0 / sqrt(3.0)) * (r1m1_31_0 * cphi * sin2Th + r1m1_3m1_0 * c2phi * sinSqTh); + GDouble Wt1m1_0 = (2.0 / sqrt(3.0)) * (rt1m1_31_0 * sphi * sin2Th + rt1m1_3m1_0 * s2phi * sinSqTh); + + GDouble W10_0 = r10_33_0 * sinSqTh + r10_11_0 * (1.0/3.0 + cosSqTh); + GDouble Wb10_0 = (2.0 / sqrt(3.0)) * (r10_31_0 * cphi * sin2Th + r10_3m1_0 * c2phi * sinSqTh); + GDouble Wt10_0 = (2.0 / sqrt(3.0)) * (rt10_31_0 * sphi * sin2Th + rt10_3m1_0 * s2phi * sinSqTh); + GDouble W0 = cosSqThPi * (W00_0 - Wb00_0) + sinSqThPi * (W11_0 - Wb11_0) - sinSqThPi * (c2phiPi * (W1m1_0 - Wb1m1_0) + s2phiPi * Wt1m1_0) - sqrt(2.0) * sin2ThPi * (cphiPi * (W10_0 - Wb10_0) + sphiPi * Wt10_0); + + //alpha = 1 (same structure) + GDouble W00_1 = r00_33_1 * sinSqTh + r00_11_1 * (1.0/3.0 + cosSqTh); + GDouble Wb00_1 = (2.0 / sqrt(3.0)) * (r00_31_1 * cphi * sin2Th + r00_3m1_1 * c2phi * sinSqTh); + + GDouble W11_1 = r11_33_1 * sinSqTh + r11_11_1 * (1.0/3.0 + cosSqTh); + GDouble Wb11_1 = (2.0 / sqrt(3.0)) * (r11_31_1 * cphi * sin2Th + r11_3m1_1 * c2phi * sinSqTh); + + GDouble W1m1_1 = r1m1_33_1 * sinSqTh + r1m1_11_1 * (1.0/3.0 + cosSqTh); + GDouble Wb1m1_1 = (2.0 / sqrt(3.0)) * (r1m1_31_1 * cphi * sin2Th + r1m1_3m1_1 * c2phi * sinSqTh); + GDouble Wt1m1_1 = (2.0 / sqrt(3.0)) * (rt1m1_31_1 * sphi * sin2Th + rt1m1_3m1_1 * s2phi * sinSqTh); + + GDouble W10_1 = r10_33_1 * sinSqTh + r10_11_1 * (1.0/3.0 + cosSqTh); + GDouble Wb10_1 = (2.0 / sqrt(3.0)) * (r10_31_1 * cphi * sin2Th + r10_3m1_1 * c2phi * sinSqTh); + GDouble Wt10_1 = (2.0 / sqrt(3.0)) * (rt10_31_1 * sphi * sin2Th + rt10_3m1_1 * s2phi * sinSqTh); + + GDouble W1 =cosSqThPi * (W00_1 - Wb00_1) + sinSqThPi * (W11_1 - Wb11_1) - sinSqThPi * (c2phiPi * (W1m1_1 - Wb1m1_1) + s2phiPi * Wt1m1_1) - sqrt(2.0) * sin2ThPi * (cphiPi * (W10_1 - Wb10_1) + sphiPi * Wt10_1); + + // alpha = 2 (different structure) + GDouble Wt00_2 = (2.0 / sqrt(3.0)) * (rt00_31_2 * sphi * sin2Th + rt00_3m1_2 * s2phi * sinSqTh); + GDouble Wt11_2 = (2.0 / sqrt(3.0)) * (rt11_31_2 * sphi * sin2Th + rt11_3m1_2 * s2phi * sinSqTh); + + GDouble W1m1_2 = r1m1_33_2 * sinSqTh + r1m1_11_2 * (1.0/3.0 + cosSqTh); + GDouble Wb1m1_2 = (2.0 / sqrt(3.0)) * (r1m1_31_2 * cphi * sin2Th + r1m1_3m1_2 * c2phi * sinSqTh); + GDouble Wt1m1_2 = (2.0 / sqrt(3.0)) * (rt1m1_31_2 * sphi * sin2Th + rt1m1_3m1_2 * s2phi * sinSqTh); + + GDouble W10_2 = r10_33_2 * sinSqTh + r10_11_2 * (1.0/3.0 + cosSqTh); + GDouble Wb10_2 = (2.0 / sqrt(3.0)) * (r10_31_2 * cphi * sin2Th + r10_3m1_2 * c2phi * sinSqTh); + GDouble Wt10_2 = (2.0 / sqrt(3.0)) * (rt10_31_2 * sphi * sin2Th + rt10_3m1_2 * s2phi * sinSqTh); + GDouble W2 = cosSqThPi * Wt00_2+ sinSqThPi * Wt11_2+ sinSqThPi * (s2phiPi * (W1m1_2 - Wb1m1_2) - c2phiPi * Wt1m1_2) + sqrt(2.0) * sin2ThPi * (sphiPi * (W10_2 - Wb10_2) + cphiPi * Wt10_2); + + GDouble Wpol = W0 - Pgamma * cos(2.0 * bigPhi) * W1 - Pgamma * sin(2.0 * bigPhi) * W2; + + WCUComplex amp = { sqrt(fabs(Wpol)), 0.0 }; + pcDevAmp[iEvent] = amp; + +} + + +void +GPUTwoPiAngles_Delta_DoubleSDMEs_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble r00_33_0, + GDouble r11_33_0, GDouble r11_11_0, + GDouble r1m1_33_0, GDouble r1m1_11_0, + GDouble r10_33_0, GDouble r10_11_0, + GDouble r00_31_0, GDouble r00_3m1_0, + GDouble r11_31_0, GDouble r11_3m1_0, + GDouble r1m1_31_0, GDouble r1m1_3m1_0, + GDouble r10_31_0, GDouble r10_3m1_0, + GDouble rt1m1_31_0, GDouble rt1m1_3m1_0, + GDouble rt10_31_0, GDouble rt10_3m1_0, + GDouble r00_33_1, GDouble r00_11_1, + GDouble r11_33_1, GDouble r11_11_1, + GDouble r1m1_33_1, GDouble r1m1_11_1, + GDouble r10_33_1, GDouble r10_11_1, + GDouble r00_31_1, GDouble r00_3m1_1, + GDouble r11_31_1, GDouble r11_3m1_1, + GDouble r1m1_31_1, GDouble r1m1_3m1_1, + GDouble r10_31_1, GDouble r10_3m1_1, + GDouble rt1m1_31_1, GDouble rt1m1_3m1_1, + GDouble rt10_31_1, GDouble rt10_3m1_1, + GDouble rt00_31_2, GDouble rt00_3m1_2, + GDouble rt11_31_2, GDouble rt11_3m1_2, + GDouble r1m1_33_2, GDouble r1m1_11_2, + GDouble r1m1_31_2, GDouble r1m1_3m1_2, + GDouble rt1m1_31_2, GDouble rt1m1_3m1_2, + GDouble r10_33_2, GDouble r10_11_2, + GDouble r10_31_2, GDouble r10_3m1_2, + GDouble rt10_31_2, GDouble rt10_3m1_2 +) +{ + + GPUTwoPiAngles_Delta_DoubleSDMEs_kernel<<< dimGrid, dimBlock >>> + ( GPU_AMP_ARGS, + r00_33_0, + r11_33_0, r11_11_0, + r1m1_33_0, r1m1_11_0, + r10_33_0, r10_11_0, + r00_31_0, r00_3m1_0, + r11_31_0, r11_3m1_0, + r1m1_31_0, r1m1_3m1_0, + r10_31_0, r10_3m1_0, + rt1m1_31_0, rt1m1_3m1_0, + rt10_31_0, rt10_3m1_0, + r00_33_1, r00_11_1, + r11_33_1, r11_11_1, + r1m1_33_1, r1m1_11_1, + r10_33_1, r10_11_1, + r00_31_1, r00_3m1_1, + r11_31_1, r11_3m1_1, + r1m1_31_1, r1m1_3m1_1, + r10_31_1, r10_3m1_1, + rt1m1_31_1, rt1m1_3m1_1, + rt10_31_1, rt10_3m1_1, + rt00_31_2, rt00_3m1_2, + rt11_31_2, rt11_3m1_2, + r1m1_33_2, r1m1_11_2, + r1m1_31_2, r1m1_3m1_2, + rt1m1_31_2, rt1m1_3m1_2, + r10_33_2, r10_11_2, + r10_31_2, r10_3m1_2, + rt10_31_2, rt10_3m1_2 + + + ); +} diff --git a/src/libraries/AMPTOOLS_AMPS/GPUTwoPiAngles_Delta_factorized_kernel.cu b/src/libraries/AMPTOOLS_AMPS/GPUTwoPiAngles_Delta_factorized_kernel.cu new file mode 100755 index 000000000..e19482fc4 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/GPUTwoPiAngles_Delta_factorized_kernel.cu @@ -0,0 +1,96 @@ + +#include + +#include "GPUManager/GPUCustomTypes.h" +#include "GPUManager/CUDA-Complex.cuh" + + +__global__ void +// +GPUTwoPiAngles_Delta_factorized_kernel( GPU_AMP_PROTO, +GDouble rho000, +GDouble rho100, +GDouble rho1m10, +GDouble rho001, +GDouble rho111, +GDouble rho101, +GDouble rho1m11, +GDouble rho102, +GDouble rho1m12, +GDouble delta_rho011, +GDouble delta_rho031, +GDouble delta_rho03m1 +) + { + + int iEvent = GPU_THIS_EVENT; + + // here we need to be careful to index the user-defined + // data with the proper integer corresponding to the + // enumeration in the C++ header file + + + GDouble Pgamma = GPU_UVARS(0); + GDouble cosTheta_pim = GPU_UVARS(1); + GDouble sinSqTheta_pim = GPU_UVARS(2); + GDouble sin2Theta_pim = GPU_UVARS(3); + GDouble phi_pim = GPU_UVARS(4); + GDouble cosTheta_proton = GPU_UVARS(5); + GDouble sinSqTheta_proton = GPU_UVARS(6); + GDouble sin2Theta_proton = GPU_UVARS(7); + GDouble phi_proton = GPU_UVARS(8); + GDouble bigPhi = GPU_UVARS(9); + // flag for pol (0.0 = no pol, 1.0 = pol) + + GDouble W0_meson = 3.*(0.5*(1. - rho000) + 0.5*(3.*rho000 - 1.)*cosTheta_pim*cosTheta_pim - sqrt(2.)*rho100*sin2Theta_pim*cos(phi_pim) - rho1m10*sinSqTheta_pim*cos(2.*phi_pim)); + GDouble W1_meson= 3.*(rho111*sinSqTheta_pim + rho001*cosTheta_pim*cosTheta_pim - sqrt(2.)*rho101*sin2Theta_pim*cos(phi_pim) - rho1m11*sinSqTheta_pim*cos(2.*phi_pim)); + GDouble W2_meson = 3. * (sqrt(2.)*rho102*sin2Theta_pim*sin(phi_pim) + rho1m12*sinSqTheta_pim*sin(2.*phi_pim)); + + GDouble W0_baryon = 3.*(0.5 - delta_rho011)*sinSqTheta_proton + delta_rho011*(1.+3.*cosTheta_proton*cosTheta_proton) - 2.*sqrt(3.)*delta_rho031*cos(phi_proton)*sin2Theta_proton - 2.*sqrt(3.)*delta_rho03m1*cos(2.*phi_proton)*sinSqTheta_proton; + // GDouble W1_baryon = (3.*delta_rho133*sinSqTheta_piplus + delta_rho111*(1.+3.*cosTheta_piplus*cosTheta_piplus) - 2.*sqrt(3.)*delta_rho131*cos(phi_piplus)*sin2Theta_piplus - 2.*sqrt(3.)*delta_rho13m1*cos(2.*phi_piplus)*sinSqTheta_piplus); + // GDouble W2_baryon = (2.*sqrt(3.)*delta_rho231*sin(phi_piplus)*sin2Theta_piplus + 2.*sqrt(3.)*delta_rho23m1*sin(2.*phi_piplus)*sinSqTheta_piplus); + + GDouble Wpol = W0_meson*W0_baryon - Pgamma * cos(2.0 * bigPhi) *( W1_meson*W0_baryon) - Pgamma * sin(2.0 * bigPhi) *( W2_meson*W0_baryon); + + Wpol *= 1/(4.*PI); + + WCUComplex amp = { sqrt(fabs(Wpol)), 0.0 }; + pcDevAmp[iEvent] = amp; + +} + + +void +GPUTwoPiAngles_Delta_factorized_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble rho000, + GDouble rho100, + GDouble rho1m10, + GDouble rho001, + GDouble rho111, + GDouble rho101, + GDouble rho1m11, + GDouble rho102, + GDouble rho1m12, + GDouble delta_rho011, + GDouble delta_rho031, + GDouble delta_rho03m1 + +) +{ + + GPUTwoPiAngles_Delta_factorized_kernel<<< dimGrid, dimBlock >>> + ( GPU_AMP_ARGS, + rho000, + rho100, + rho1m10, + rho001, + rho111, + rho101, + rho1m11, + rho102, + rho1m12, + delta_rho011, + delta_rho031, + delta_rho03m1 + ); +} diff --git a/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.cc b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.cc new file mode 100755 index 000000000..6b022a82e --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.cc @@ -0,0 +1,458 @@ + +#include +#include +#include +#include +#include + +#include "TLorentzVector.h" +#include "TLorentzRotation.h" +#include "TFile.h" + +#include "IUAmpTools/Kinematics.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h" +#include "AMPTOOLS_AMPS/clebschGordan.h" +#include "AMPTOOLS_AMPS/wignerD.h" + +TwoPiAngles_Delta_DoubleSDMEs::TwoPiAngles_Delta_DoubleSDMEs( const vector< string >& args ) : + UserAmplitude< TwoPiAngles_Delta_DoubleSDMEs >( args ) +{ + +// assert(args.size() == 22 || args.size() == 23 || args.size() == 24 ); +frame = string( args[0] ); +polAngle = AmpParameter( args[1] ); + +polFraction = atof(args[2].c_str()); +cout << "Fitting with constant polarization " << polFraction << endl; +// r(t)m,m'_lambda,lambda'_alpha (rt corresponds to rtilde, m corresponds to rho helicity, lambda to Delta helicity, alpha to photon polarization) +r00_33_0 = AmpParameter( args[3] ); + +r11_33_0 = AmpParameter( args[4] ); +r11_11_0 = AmpParameter( args[5] ); + +r1m1_33_0 = AmpParameter( args[6] ); +r1m1_11_0 = AmpParameter( args[7] ); + +r10_33_0 = AmpParameter( args[8] ); +r10_11_0 = AmpParameter( args[9] ); + +r00_31_0 = AmpParameter( args[10] ); +r00_3m1_0 = AmpParameter( args[11] ); + +r11_31_0 = AmpParameter( args[12] ); +r11_3m1_0 = AmpParameter( args[13] ); + +r1m1_31_0 = AmpParameter( args[14] ); +r1m1_3m1_0 = AmpParameter( args[15] ); + +r10_31_0 = AmpParameter( args[16] ); +r10_3m1_0 = AmpParameter( args[17] ); + +rt1m1_31_0 = AmpParameter( args[18] ); +rt1m1_3m1_0 = AmpParameter( args[19] ); + +rt10_31_0 = AmpParameter( args[20] ); +rt10_3m1_0 = AmpParameter( args[21] ); + +//alpha = 1 (same structure) +r00_33_1 = AmpParameter( args[22] ); +r00_11_1 = AmpParameter( args[23] ); + +r11_33_1 = AmpParameter( args[24] ); +r11_11_1 = AmpParameter( args[25] ); + +r1m1_33_1 = AmpParameter( args[26] ); +r1m1_11_1 = AmpParameter( args[27] ); + +r10_33_1 = AmpParameter( args[28] ); +r10_11_1 = AmpParameter( args[29] ); + +r00_31_1 = AmpParameter( args[30] ); +r00_3m1_1 = AmpParameter( args[31] ); + +r11_31_1 = AmpParameter( args[32] ); +r11_3m1_1 = AmpParameter( args[33] ); + +r1m1_31_1 = AmpParameter( args[34] ); +r1m1_3m1_1 = AmpParameter( args[35] ); + +r10_31_1 = AmpParameter( args[36] ); +r10_3m1_1 = AmpParameter( args[37] ); + +rt1m1_31_1 = AmpParameter( args[38] ); +rt1m1_3m1_1 = AmpParameter( args[39] ); + +rt10_31_1 = AmpParameter( args[40] ); +rt10_3m1_1 = AmpParameter( args[41] ); + +//alpha = 2 (different structure) +rt00_31_2 = AmpParameter( args[42] ); +rt00_3m1_2 = AmpParameter( args[43] ); + +rt11_31_2 = AmpParameter( args[44] ); +rt11_3m1_2 = AmpParameter( args[45] ); + +r1m1_33_2 = AmpParameter( args[46] ); +r1m1_11_2 = AmpParameter( args[47] ); + +r1m1_31_2 = AmpParameter( args[48] ); +r1m1_3m1_2 = AmpParameter( args[49] ); + +rt1m1_31_2 = AmpParameter( args[50] ); +rt1m1_3m1_2 = AmpParameter( args[51] ); + +r10_33_2 = AmpParameter( args[52] ); +r10_11_2 = AmpParameter( args[53] ); + +r10_31_2 = AmpParameter( args[54] ); +r10_3m1_2 = AmpParameter( args[55] ); + +rt10_31_2 = AmpParameter( args[56] ); +rt10_3m1_2 = AmpParameter( args[57] ); + +// Parameter r00_11_0 was eliminated due to normalization constraint, so we don't register it); +registerParameter( r00_33_0 ); + +registerParameter( r11_33_0 ); +registerParameter( r11_11_0 ); + +registerParameter( r1m1_33_0 ); +registerParameter( r1m1_11_0 ); + +registerParameter( r10_33_0 ); +registerParameter( r10_11_0 ); + +registerParameter( r00_31_0 ); +registerParameter( r00_3m1_0 ); + +registerParameter( r11_31_0 ); +registerParameter( r11_3m1_0 ); + +registerParameter( r1m1_31_0 ); +registerParameter( r1m1_3m1_0 ); + +registerParameter( r10_31_0 ); +registerParameter( r10_3m1_0 ); + +registerParameter( rt1m1_31_0 ); +registerParameter( rt1m1_3m1_0 ); + +registerParameter( rt10_31_0 ); +registerParameter( rt10_3m1_0 ); + +registerParameter( r00_33_1 ); +registerParameter( r00_11_1 ); + +registerParameter( r11_33_1 ); +registerParameter( r11_11_1 ); + +registerParameter( r1m1_33_1 ); +registerParameter( r1m1_11_1 ); + +registerParameter( r10_33_1 ); +registerParameter( r10_11_1 ); + +registerParameter( r00_31_1 ); +registerParameter( r00_3m1_1 ); + +registerParameter( r11_31_1 ); +registerParameter( r11_3m1_1 ); + +registerParameter( r1m1_31_1 ); +registerParameter( r1m1_3m1_1 ); + +registerParameter( r10_31_1 ); +registerParameter( r10_3m1_1 ); + +registerParameter( rt1m1_31_1 ); +registerParameter( rt1m1_3m1_1 ); + +registerParameter( rt10_31_1 ); +registerParameter( rt10_3m1_1 ); + +//alpha = 2 (different structure) +registerParameter( rt00_31_2 ); +registerParameter( rt00_3m1_2 ); + +registerParameter( rt11_31_2 ); +registerParameter( rt11_3m1_2 ); + +registerParameter( r1m1_33_2 ); +registerParameter( r1m1_11_2 ); + +registerParameter( r1m1_31_2 ); +registerParameter( r1m1_3m1_2 ); + +registerParameter( rt1m1_31_2 ); +registerParameter( rt1m1_3m1_2 ); + +registerParameter( r10_33_2 ); +registerParameter( r10_11_2 ); + +registerParameter( r10_31_2 ); +registerParameter( r10_3m1_2 ); + +registerParameter( rt10_31_2 ); +registerParameter( rt10_3m1_2 ); + +registerParameter( polAngle ); + + // } + // Two possibilities to initialize this amplitude: + // 1: 11 arguments, fixed polarization + // Usage: amplitude :::: TwoPiAngles ... + + // 2: 12 arguments, read polarization from histogram in file + // Usage: amplitude :::: TwoPiAngles ... +// else if(args.size() == 13) { +// polFraction = 0.; +// TFile* f = new TFile( args[11].c_str() ); +// polFrac_vs_E = (TH1D*)f->Get( args[12].c_str() ); +// assert( polFrac_vs_E != NULL ); +// cout << "Fitting with polarization from " << polFrac_vs_E->GetName() << endl; +// } +// +} + + +complex< GDouble > +TwoPiAngles_Delta_DoubleSDMEs::calcAmplitude( GDouble** pKin, GDouble* userVars ) const { + + +//GDouble bigPhi = polAngle*0.017453293 + userVars[kBigPhi]; // rotate Phi (in rad) +GDouble bigPhi = userVars[kBigPhi]; +GDouble Pgamma = userVars[kPgamma]; +GDouble cosTheta_pim = userVars[kCosTheta_pim]; +GDouble sinSqTheta_pim = userVars[kSinSqTheta_pim]; +GDouble sin2Theta_pim = userVars[kSin2Theta_pim]; +GDouble phi_pim = userVars[kPhi_pim]; +GDouble cosTheta_proton = userVars[kCosTheta_proton]; +GDouble sinSqTheta_proton = userVars[kSinSqTheta_proton]; +GDouble sin2Theta_proton = userVars[kSin2Theta_proton]; +GDouble phi_proton = userVars[kPhi_proton]; + + +GDouble sinSqTh = sinSqTheta_proton; +GDouble sin2Th = sin2Theta_proton; +GDouble cosSqTh = cosTheta_proton * cosTheta_proton; + +GDouble cphi = cos(phi_proton); +GDouble sphi = sin(phi_proton); +GDouble c2phi = cos(2.0 * phi_proton); +GDouble s2phi = sin(2.0 * phi_proton); + +GDouble sinSqThPi = sinSqTheta_pim; +GDouble sin2ThPi = sin2Theta_pim; +GDouble cosSqThPi = cosTheta_pim * cosTheta_pim; + +GDouble cphiPi = cos(phi_pim); +GDouble sphiPi = sin(phi_pim); +GDouble c2phiPi = cos(2.0 * phi_pim); +GDouble s2phiPi = sin(2.0 * phi_pim); + +GDouble r00_11_0 = 0.5 - r00_33_0 - 2.0 * r11_11_0 - 2.0 * r11_33_0; + + +GDouble W00_0 = r00_33_0 * sinSqTh + r00_11_0 * (1.0/3.0 + cosSqTh); +GDouble Wb00_0 = (2.0 / sqrt(3.0)) * (r00_31_0 * cphi * sin2Th + r00_3m1_0 * c2phi * sinSqTh); + +GDouble W11_0 = r11_33_0 * sinSqTh + r11_11_0 * (1.0/3.0 + cosSqTh); +GDouble Wb11_0 = (2.0 / sqrt(3.0)) * (r11_31_0 * cphi * sin2Th + r11_3m1_0 * c2phi * sinSqTh); + +GDouble W1m1_0 = r1m1_33_0 * sinSqTh + r1m1_11_0 * (1.0/3.0 + cosSqTh); +GDouble Wb1m1_0 = (2.0 / sqrt(3.0)) * (r1m1_31_0 * cphi * sin2Th + r1m1_3m1_0 * c2phi * sinSqTh); +GDouble Wt1m1_0 = (2.0 / sqrt(3.0)) * (rt1m1_31_0 * sphi * sin2Th + rt1m1_3m1_0 * s2phi * sinSqTh); + +GDouble W10_0 = r10_33_0 * sinSqTh + r10_11_0 * (1.0/3.0 + cosSqTh); +GDouble Wb10_0 = (2.0 / sqrt(3.0)) * (r10_31_0 * cphi * sin2Th + r10_3m1_0 * c2phi * sinSqTh); +GDouble Wt10_0 = (2.0 / sqrt(3.0)) * (rt10_31_0 * sphi * sin2Th + rt10_3m1_0 * s2phi * sinSqTh); + +GDouble W0 = cosSqThPi * (W00_0 - Wb00_0) + sinSqThPi * (W11_0 - Wb11_0) - sinSqThPi * (c2phiPi * (W1m1_0 - Wb1m1_0) + s2phiPi * Wt1m1_0) - sqrt(2.0) * sin2ThPi * (cphiPi * (W10_0 - Wb10_0) + sphiPi * Wt10_0); +//alpha = 1 (same structure) +GDouble W00_1 = r00_33_1 * sinSqTh + r00_11_1 * (1.0/3.0 + cosSqTh); +GDouble Wb00_1 = (2.0 / sqrt(3.0)) * (r00_31_1 * cphi * sin2Th + r00_3m1_1 * c2phi * sinSqTh); + +GDouble W11_1 = r11_33_1 * sinSqTh + r11_11_1 * (1.0/3.0 + cosSqTh); +GDouble Wb11_1 = (2.0 / sqrt(3.0)) * (r11_31_1 * cphi * sin2Th + r11_3m1_1 * c2phi * sinSqTh); + +GDouble W1m1_1 = r1m1_33_1 * sinSqTh + r1m1_11_1 * (1.0/3.0 + cosSqTh); +GDouble Wb1m1_1 = (2.0 / sqrt(3.0)) * (r1m1_31_1 * cphi * sin2Th + r1m1_3m1_1 * c2phi * sinSqTh); +GDouble Wt1m1_1 = (2.0 / sqrt(3.0)) * (rt1m1_31_1 * sphi * sin2Th + rt1m1_3m1_1 * s2phi * sinSqTh); + +GDouble W10_1 = r10_33_1 * sinSqTh + r10_11_1 * (1.0/3.0 + cosSqTh); +GDouble Wb10_1 = (2.0 / sqrt(3.0)) * (r10_31_1 * cphi * sin2Th + r10_3m1_1 * c2phi * sinSqTh); +GDouble Wt10_1 = (2.0 / sqrt(3.0)) * (rt10_31_1 * sphi * sin2Th + rt10_3m1_1 * s2phi * sinSqTh); + +GDouble W1 =cosSqThPi * (W00_1 - Wb00_1) + sinSqThPi * (W11_1 - Wb11_1) - sinSqThPi * (c2phiPi * (W1m1_1 - Wb1m1_1) + s2phiPi * Wt1m1_1) - sqrt(2.0) * sin2ThPi * (cphiPi * (W10_1 - Wb10_1) + sphiPi * Wt10_1); + +// alpha = 2 (different structure) +GDouble Wt00_2 = (2.0 / sqrt(3.0)) * (rt00_31_2 * sphi * sin2Th + rt00_3m1_2 * s2phi * sinSqTh); +GDouble Wt11_2 = (2.0 / sqrt(3.0)) * (rt11_31_2 * sphi * sin2Th + rt11_3m1_2 * s2phi * sinSqTh); + +GDouble W1m1_2 = r1m1_33_2 * sinSqTh + r1m1_11_2 * (1.0/3.0 + cosSqTh); +GDouble Wb1m1_2 = (2.0 / sqrt(3.0)) * (r1m1_31_2 * cphi * sin2Th + r1m1_3m1_2 * c2phi * sinSqTh); +GDouble Wt1m1_2 = (2.0 / sqrt(3.0)) * (rt1m1_31_2 * sphi * sin2Th + rt1m1_3m1_2 * s2phi * sinSqTh); + +GDouble W10_2 = r10_33_2 * sinSqTh + r10_11_2 * (1.0/3.0 + cosSqTh); +GDouble Wb10_2 = (2.0 / sqrt(3.0)) * (r10_31_2 * cphi * sin2Th + r10_3m1_2 * c2phi * sinSqTh); +GDouble Wt10_2 = (2.0 / sqrt(3.0)) * (rt10_31_2 * sphi * sin2Th + rt10_3m1_2 * s2phi * sinSqTh); +GDouble W2 = cosSqThPi * Wt00_2+ sinSqThPi * Wt11_2+ sinSqThPi * (s2phiPi * (W1m1_2 - Wb1m1_2) - c2phiPi * Wt1m1_2) + sqrt(2.0) * sin2ThPi * (sphiPi * (W10_2 - Wb10_2) + cphiPi * Wt10_2); + +GDouble Wpol = W0 - Pgamma * cos(2.0 * bigPhi) * W1 - Pgamma * sin(2.0 * bigPhi) * W2; +return complex< GDouble > ( sqrt(fabs(Wpol)) ); +} + +void +TwoPiAngles_Delta_DoubleSDMEs::calcUserVars( GDouble** pKin, GDouble* userVars ) const { + + TLorentzVector beam ( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); + TLorentzVector proton ( pKin[1][1], pKin[1][2], pKin[1][3], pKin[1][0] ); + TLorentzVector p1 ( pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0] ); // pip + TLorentzVector p2 ( pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0] ); // pim + TLorentzVector p3 ( pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0] ); // pi0 + TLorentzVector target ( 0, 0, 0, 0.9382720813); + // determine boost vectors + TLorentzVector rho = p2 + p3; + TLorentzVector Delta = proton + p1; + TLorentzVector CM_motion_lab = beam + target; + TLorentzVector piminus = p2; + + TLorentzRotation rhoBoost( -rho.BoostVector() ); + TLorentzRotation CMBoost( -CM_motion_lab.BoostVector() ); + TLorentzRotation DeltaBoost( -Delta.BoostVector() ); + +//bost in rho restframe + TLorentzVector beam_resrho = rhoBoost * beam; + TLorentzVector Delta_resrho = rhoBoost * Delta; + TLorentzVector piminus_resrho = rhoBoost * piminus; + TLorentzVector proton_resrho = rhoBoost * proton; + +//boost in CM frame + TLorentzVector beam_cm = CMBoost * beam; + TLorentzVector Delta_cm = CMBoost * Delta; + TLorentzVector piminus_cm = CMBoost * piminus; + TLorentzVector rho_cm = CMBoost * rho; + +//boost in recoil frame + TLorentzVector target_resDelta = DeltaBoost * target; + TLorentzVector proton_resDelta = DeltaBoost * proton; + TLorentzVector rho_resDelta = DeltaBoost * rho; + +// Define coordinate systems: +// y equal in all frames + TVector3 y = (beam.Vect().Unit().Cross(-Delta.Vect().Unit())).Unit(); + + //rho restframe + + //Helicity frame + TVector3 z_hel_rho = -1. * Delta_resrho.Vect().Unit(); + TVector3 x_hel_rho = y.Cross(z_hel_rho).Unit(); + + // GJ frame + TVector3 z_GJ_rho = beam_resrho.Vect().Unit(); + TVector3 x_GJ_rho = y.Cross(z_GJ_rho).Unit(); + + //Delta system + // Helicity frame + TVector3 z_hel_Delta = -rho_resDelta.Vect().Unit(); + TVector3 x_hel_Delta = (y.Cross(z_hel_Delta)).Unit(); + + // GJ frame + TVector3 z_GJ_Delta = target_resDelta.Vect().Unit(); + TVector3 x_GJ_Delta = y.Cross(z_GJ_Delta).Unit(); + + + //Angle calculation + TVector3 angles_piminus_hel_rho( (piminus_resrho.Vect()).Dot(x_hel_rho), + (piminus_resrho.Vect()).Dot(y), + (piminus_resrho.Vect()).Dot(z_hel_rho) ); + + TVector3 angles_piminus_GJ_rho( (piminus_resrho.Vect()).Dot(x_GJ_rho), + (piminus_resrho.Vect()).Dot(y), + (piminus_resrho.Vect()).Dot(z_GJ_rho) ); + + TVector3 angles_proton_hel_Delta( (proton_resDelta.Vect()).Dot(x_hel_Delta), + (proton_resDelta.Vect()).Dot(y), + (proton_resDelta.Vect()).Dot(z_hel_Delta) ); + + TVector3 angles_proton_GJ_Delta( (proton_resDelta.Vect()).Dot(x_GJ_Delta), + (proton_resDelta.Vect()).Dot(y), + (proton_resDelta.Vect()).Dot(z_GJ_Delta) ); + + + if(frame=="Hel"){ + + userVars[kCosTheta_pim] = angles_piminus_hel_rho.CosTheta(); + userVars[kSinSqTheta_pim] = sin(angles_piminus_hel_rho.Theta())*sin(angles_piminus_hel_rho.Theta()); + userVars[kSin2Theta_pim] = sin(2.*angles_piminus_hel_rho.Theta()); + userVars[kPhi_pim] = angles_piminus_hel_rho.Phi(); + + userVars[kCosTheta_proton] = angles_proton_hel_Delta.CosTheta(); + userVars[kSinSqTheta_proton] = sin(angles_proton_hel_Delta.Theta())*sin(angles_proton_hel_Delta.Theta()); + userVars[kSin2Theta_proton] = sin(2.*angles_proton_hel_Delta.Theta()); + userVars[kPhi_proton] = angles_proton_hel_Delta.Phi(); + } + else{ + userVars[kCosTheta_pim] = angles_piminus_GJ_rho.CosTheta(); + userVars[kSinSqTheta_pim] = sin(angles_piminus_GJ_rho.Theta())*sin(angles_piminus_GJ_rho.Theta()); + userVars[kSin2Theta_pim] = sin(2.*angles_piminus_GJ_rho.Theta()); + userVars[kPhi_pim] = angles_piminus_GJ_rho.Phi(); + + userVars[kCosTheta_proton] = angles_proton_GJ_Delta.CosTheta(); + userVars[kSinSqTheta_proton] = sin(angles_proton_GJ_Delta.Theta())*sin(angles_proton_GJ_Delta.Theta()); + userVars[kSin2Theta_proton] = sin(2.*angles_proton_GJ_Delta.Theta()); + userVars[kPhi_proton] = angles_proton_GJ_Delta.Phi(); + } + TVector3 eps(cos(polAngle*TMath::DegToRad()), sin(polAngle*TMath::DegToRad()), 0.0); + userVars[kBigPhi] = atan2(y.Dot(eps), beam.Vect().Unit().Dot(eps.Cross(y))); + + // vector meson production from K. Schilling et. al. + GDouble Pgamma; + Pgamma = polFraction; + // else{ + // int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); + // if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){ + // Pgamma = 0.; + // } + // else Pgamma = polFrac_vs_E->GetBinContent(bin); + // } + userVars[kPgamma] = Pgamma; + +} + +#ifdef GPU_ACCELERATION +void +TwoPiAngles_Delta_DoubleSDMEs::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const { + + GPUTwoPiAngles_Delta_DoubleSDMEs_exec( dimGrid, dimBlock, GPU_AMP_ARGS, + r00_33_0, + r11_33_0, r11_11_0, + r1m1_33_0, r1m1_11_0, + r10_33_0, r10_11_0, + r00_31_0, r00_3m1_0, + r11_31_0, r11_3m1_0, + r1m1_31_0, r1m1_3m1_0, + r10_31_0, r10_3m1_0, + rt1m1_31_0, rt1m1_3m1_0, + rt10_31_0, rt10_3m1_0, + r00_33_1, r00_11_1, + r11_33_1, r11_11_1, + r1m1_33_1, r1m1_11_1, + r10_33_1, r10_11_1, + r00_31_1, r00_3m1_1, + r11_31_1, r11_3m1_1, + r1m1_31_1, r1m1_3m1_1, + r10_31_1, r10_3m1_1, + rt1m1_31_1, rt1m1_3m1_1, + rt10_31_1, rt10_3m1_1, + rt00_31_2, rt00_3m1_2, + rt11_31_2, rt11_3m1_2, + r1m1_33_2, r1m1_11_2, + r1m1_31_2, r1m1_3m1_2, + rt1m1_31_2, rt1m1_3m1_2, + r10_33_2, r10_11_2, + r10_31_2, r10_3m1_2, + rt10_31_2, rt10_3m1_2 +); +} +#endif //GPU_ACCELERATION \ No newline at end of file diff --git a/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h new file mode 100755 index 000000000..c58d57d58 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h @@ -0,0 +1,166 @@ +#if !defined(TWOPIANGLES_DELTA_DOUBLESDMES) +#define TWOPIANGLES_DELTA_DOUBLESDMES + +#include "IUAmpTools/Amplitude.h" +#include "IUAmpTools/UserAmplitude.h" +#include "IUAmpTools/AmpParameter.h" +#include "GPUManager/GPUCustomTypes.h" + +#include "TH1D.h" +#include +#include +#include + +#ifdef GPU_ACCELERATION +void +GPUTwoPiAngles_Delta_DoubleSDMEs_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, +GDouble r00_33_0, //GDouble r00_11_0 +GDouble r11_33_0, GDouble r11_11_0, +GDouble r1m1_33_0, GDouble r1m1_11_0, +GDouble r10_33_0, GDouble r10_11_0, + +GDouble r00_31_0, GDouble r00_3m1_0, +GDouble r11_31_0, GDouble r11_3m1_0, +GDouble r1m1_31_0, GDouble r1m1_3m1_0, +GDouble r10_31_0, GDouble r10_3m1_0, + +GDouble rt1m1_31_0, GDouble rt1m1_3m1_0, +GDouble rt10_31_0, GDouble rt10_3m1_0, +GDouble r00_33_1, GDouble r00_11_1, +GDouble r11_33_1, GDouble r11_11_1, +GDouble r1m1_33_1, GDouble r1m1_11_1, +GDouble r10_33_1, GDouble r10_11_1, +GDouble r00_31_1, GDouble r00_3m1_1, +GDouble r11_31_1, GDouble r11_3m1_1, +GDouble r1m1_31_1, GDouble r1m1_3m1_1, +GDouble r10_31_1, GDouble r10_3m1_1, +GDouble rt1m1_31_1, GDouble rt1m1_3m1_1, +GDouble rt10_31_1, GDouble rt10_3m1_1, +GDouble rt00_31_2, GDouble rt00_3m1_2, +GDouble rt11_31_2, GDouble rt11_3m1_2, +GDouble r1m1_33_2, GDouble r1m1_11_2, +GDouble r1m1_31_2, GDouble r1m1_3m1_2, +GDouble rt1m1_31_2, GDouble rt1m1_3m1_2, +GDouble r10_33_2, GDouble r10_11_2, +GDouble r10_31_2, GDouble r10_3m1_2, +GDouble rt10_31_2, GDouble rt10_3m1_2 +); + +#endif // GPU_ACCELERATION + +using std::complex; +using namespace std; + +class Kinematics; + +class TwoPiAngles_Delta_DoubleSDMEs : public UserAmplitude< TwoPiAngles_Delta_DoubleSDMEs > +{ + +public: + + TwoPiAngles_Delta_DoubleSDMEs() : UserAmplitude< TwoPiAngles_Delta_DoubleSDMEs >() { }; + TwoPiAngles_Delta_DoubleSDMEs( const vector< string >& args ); + + enum UserVars { kPgamma = 0, kCosTheta_pim, kSinSqTheta_pim, kSin2Theta_pim, kPhi_pim, kCosTheta_proton, kSinSqTheta_proton, kSin2Theta_proton, kPhi_proton, + kBigPhi, kNumUserVars }; + unsigned int numUserVars() const { return kNumUserVars; } + + string name() const { return "TwoPiAngles_Delta_DoubleSDMEs"; } + + complex< GDouble > calcAmplitude( GDouble** pKin, GDouble* userVars ) const; + void calcUserVars( GDouble** pKin, GDouble* userVars ) const; + + // we can calcualte everythign we need from userVars block so allow + // the framework to purge the four-vectors + bool needsUserVarsOnly() const { return true; } + +#ifdef GPU_ACCELERATION + + void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const; + bool isGPUEnabled() const { return true; } + +#endif // GPU_ACCELERATION + +private: + + AmpParameter r00_33_0; + // AmpParameter r00_11_0; + + AmpParameter r11_33_0; + AmpParameter r11_11_0; + + AmpParameter r1m1_33_0; + AmpParameter r1m1_11_0; + + AmpParameter r10_33_0; + AmpParameter r10_11_0; + + // \bar W^0 (31 / 3-1) + + AmpParameter r00_31_0; + AmpParameter r00_3m1_0; + + AmpParameter r11_31_0; + AmpParameter r11_3m1_0; + + AmpParameter r1m1_31_0; + AmpParameter r1m1_3m1_0; + + AmpParameter r10_31_0; + AmpParameter r10_3m1_0; + + // \tilde W^0 + + AmpParameter rt1m1_31_0; + AmpParameter rt1m1_3m1_0; + + AmpParameter rt10_31_0; + AmpParameter rt10_3m1_0; + + //alpha = 1 (same structure) + AmpParameter r00_33_1; + AmpParameter r00_11_1; + AmpParameter r11_33_1; + AmpParameter r11_11_1; + AmpParameter r1m1_33_1; + AmpParameter r1m1_11_1; + AmpParameter r10_33_1; + AmpParameter r10_11_1; + AmpParameter r00_31_1; + AmpParameter r00_3m1_1; + AmpParameter r11_31_1; + AmpParameter r11_3m1_1; + AmpParameter r1m1_31_1; + AmpParameter r1m1_3m1_1; + AmpParameter r10_31_1; + AmpParameter r10_3m1_1; + AmpParameter rt1m1_31_1; + AmpParameter rt1m1_3m1_1; + AmpParameter rt10_31_1; + AmpParameter rt10_3m1_1; + //alpha = 2 (different structure) + AmpParameter rt00_31_2; + AmpParameter rt00_3m1_2; + AmpParameter rt11_31_2; + AmpParameter rt11_3m1_2; + AmpParameter r1m1_33_2; + AmpParameter r1m1_11_2; + AmpParameter r1m1_31_2; + AmpParameter r1m1_3m1_2; + AmpParameter rt1m1_31_2; + AmpParameter rt1m1_3m1_2; + AmpParameter r10_33_2; + AmpParameter r10_11_2; + AmpParameter r10_31_2; + AmpParameter r10_3m1_2; + AmpParameter rt10_31_2; + AmpParameter rt10_3m1_2; + AmpParameter polAngle; + + double polFraction; + //TH1D *polFrac_vs_E; + string frame; + string AMO; +}; + +#endif diff --git a/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.cc b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.cc new file mode 100755 index 000000000..fbf8b3d84 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.cc @@ -0,0 +1,260 @@ +#include +#include +#include +#include +#include + +#include "TLorentzVector.h" +#include "TLorentzRotation.h" +#include "TFile.h" + +#include "IUAmpTools/Kinematics.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h" +#include "AMPTOOLS_AMPS/clebschGordan.h" +#include "AMPTOOLS_AMPS/wignerD.h" + +TwoPiAngles_Delta_factorized::TwoPiAngles_Delta_factorized( const vector< string >& args ) : + UserAmplitude< TwoPiAngles_Delta_factorized >( args ) +{ + +// assert(args.size() == 22 || args.size() == 23 || args.size() == 24 ); +frame = string( args[0] ); +polAngle = AmpParameter( args[1] ); + +polFraction = atof(args[2].c_str()); +cout << "Fitting with constant polarization " << polFraction << endl; + +rho000 = AmpParameter( args[3] ); +rho100 = AmpParameter( args[4] ); +rho1m10 = AmpParameter( args[5] ); +rho001 = AmpParameter( args[6] ); +rho111 = AmpParameter( args[7] ); +rho101 = AmpParameter( args[8] ); +rho1m11 = AmpParameter( args[9] ); +rho102 = AmpParameter( args[10] ); +rho1m12 = AmpParameter( args[11] ); + +delta_rho011 = AmpParameter( args[12] ); +delta_rho031 = AmpParameter( args[13] ); +delta_rho03m1 = AmpParameter( args[14] ); + +registerParameter( rho000 ); +registerParameter( rho100 ); +registerParameter( rho1m10 ); +registerParameter( rho001 ); +registerParameter( rho111 ); +registerParameter( rho101 ); +registerParameter( rho1m11 ); + +registerParameter( rho102 ); +registerParameter( rho1m12 ); +registerParameter( delta_rho011 ); +registerParameter( delta_rho031 ); +registerParameter( delta_rho03m1 ); + +registerParameter( polAngle ); + + // } + // Two possibilities to initialize this amplitude: + // 1: 11 arguments, fixed polarization + // Usage: amplitude :::: TwoPiAngles ... + + // 2: 12 arguments, read polarization from histogram in file + // Usage: amplitude :::: TwoPiAngles ... +// else if(args.size() == 13) { +// polFraction = 0.; +// TFile* f = new TFile( args[11].c_str() ); +// polFrac_vs_E = (TH1D*)f->Get( args[12].c_str() ); +// assert( polFrac_vs_E != NULL ); +// cout << "Fitting with polarization from " << polFrac_vs_E->GetName() << endl; +// } +// +} + + +complex< GDouble > +TwoPiAngles_Delta_factorized::calcAmplitude( GDouble** pKin, GDouble* userVars ) const { + + +//GDouble bigPhi = polAngle*0.017453293 + userVars[kBigPhi]; // rotate Phi (in rad) +GDouble bigPhi = userVars[kBigPhi]; +GDouble Pgamma = userVars[kPgamma]; +GDouble cosTheta_pim = userVars[kCosTheta_pim]; +GDouble sinSqTheta_pim = userVars[kSinSqTheta_pim]; +GDouble sin2Theta_pim = userVars[kSin2Theta_pim]; +GDouble phi_pim = userVars[kPhi_pim]; +GDouble cosTheta_proton = userVars[kCosTheta_proton]; +GDouble sinSqTheta_proton = userVars[kSinSqTheta_proton]; +GDouble sin2Theta_proton = userVars[kSin2Theta_proton]; +GDouble phi_proton = userVars[kPhi_proton]; + + + // vector meson production from K. Schilling et. al. last indize = upper + + + +// GDouble rho001_val; +// rho001_val = 2*(delta_rho111 + delta_rho133 - rho111); +// Use beam assy. to eliminate one parameter + +GDouble W0_meson = 3.*(0.5*(1. - rho000) + 0.5*(3.*rho000 - 1.)*cosTheta_pim*cosTheta_pim - sqrt(2.)*rho100*sin2Theta_pim*cos(phi_pim) - rho1m10*sinSqTheta_pim*cos(2.*phi_pim)); +GDouble W1_meson= 3.*(rho111*sinSqTheta_pim + rho001*cosTheta_pim*cosTheta_pim - sqrt(2.)*rho101*sin2Theta_pim*cos(phi_pim) - rho1m11*sinSqTheta_pim*cos(2.*phi_pim)); +GDouble W2_meson = 3. * (sqrt(2.)*rho102*sin2Theta_pim*sin(phi_pim) + rho1m12*sinSqTheta_pim*sin(2.*phi_pim)); + +GDouble W0_baryon = 3.*(0.5 - delta_rho011)*sinSqTheta_proton + delta_rho011*(1.+3.*cosTheta_proton*cosTheta_proton) - 2.*TMath::Sqrt(3.)*delta_rho031*cos(phi_proton)*sin2Theta_proton - 2.*TMath::Sqrt(3.)*delta_rho03m1*cos(2.*phi_proton)*sinSqTheta_proton; + +GDouble Wpol = W0_meson*W0_baryon - Pgamma * cos(2.0 * bigPhi) *( W1_meson*W0_baryon) - Pgamma * sin(2.0 * bigPhi) *( W2_meson*W0_baryon); +Wpol *= 1/(4.*PI); + + + +return complex< GDouble > ( sqrt(fabs(Wpol)) ); +} + +void +TwoPiAngles_Delta_factorized::calcUserVars( GDouble** pKin, GDouble* userVars ) const { + + TLorentzVector beam ( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); + TLorentzVector proton ( pKin[1][1], pKin[1][2], pKin[1][3], pKin[1][0] ); + TLorentzVector p1 ( pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0] ); // pip + TLorentzVector p2 ( pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0] ); // pim + TLorentzVector p3 ( pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0] ); // pi0 + TLorentzVector target ( 0, 0, 0, 0.9382720813); + // determine boost vectors + TLorentzVector recoil = proton + p1; + TLorentzVector rho = p2 + p3; + TLorentzVector Delta = proton + p1; + TLorentzVector CM_motion_lab = beam + target; + TLorentzVector piminus = p2; + + TLorentzRotation rhoBoost( -rho.BoostVector() ); + TLorentzRotation CMBoost( -CM_motion_lab.BoostVector() ); + TLorentzRotation DeltaBoost( -Delta.BoostVector() ); + +//bost in rho restframe + TLorentzVector beam_resrho = rhoBoost * beam; + TLorentzVector Delta_resrho = rhoBoost * Delta; + TLorentzVector piminus_resrho = rhoBoost * piminus; + TLorentzVector proton_resrho = rhoBoost * proton; +//boost in CM frame + TLorentzVector beam_cm = CMBoost * beam; + TLorentzVector Delta_cm = CMBoost * Delta; + TLorentzVector piminus_cm = CMBoost * piminus; + TLorentzVector rho_cm = CMBoost * rho; + +//boost in recoil frame + TLorentzVector target_resDelta = DeltaBoost * target; + TLorentzVector proton_resDelta = DeltaBoost * proton; + TLorentzVector piplus_resDelta = DeltaBoost * p1; + TLorentzVector rho_resDelta = DeltaBoost * rho; + + // Define coordinate systems: + // y equal in all frames + TVector3 y = (beam.Vect().Unit().Cross(-Delta.Vect().Unit())).Unit(); + + //rho system + + //Helicity frame + TVector3 z_hel_rho = -1. * Delta_resrho.Vect().Unit(); + TVector3 x_hel_rho = y.Cross(z_hel_rho).Unit(); + + // GJ frame + TVector3 z_GJ_rho = beam_resrho.Vect().Unit(); + TVector3 x_GJ_rho = y.Cross(z_GJ_rho).Unit(); + + //Delta system + + // Helicity frame + TVector3 z_hel_Delta = -rho_resDelta.Vect().Unit(); + TVector3 x_hel_Delta = (y.Cross(z_hel_Delta)).Unit(); + + // GJ frame + TVector3 z_GJ_Delta = target_resDelta.Vect().Unit(); + TVector3 x_GJ_Delta = y.Cross(z_GJ_Delta).Unit(); + + //Angle calculation + TVector3 angles_piminus_hel_rho( (piminus_resrho.Vect()).Dot(x_hel_rho), + (piminus_resrho.Vect()).Dot(y), + (piminus_resrho.Vect()).Dot(z_hel_rho) ); + + TVector3 angles_piminus_GJ_rho( (piminus_resrho.Vect()).Dot(x_GJ_rho), + (piminus_resrho.Vect()).Dot(y), + (piminus_resrho.Vect()).Dot(z_GJ_rho) ); + + + + TVector3 angles_proton_hel_Delta( (proton_resDelta.Vect()).Dot(x_hel_Delta), + (proton_resDelta.Vect()).Dot(y), + (proton_resDelta.Vect()).Dot(z_hel_Delta) ); + TVector3 angles_piplus_hel_Delta( (piplus_resDelta.Vect()).Dot(x_hel_Delta), + (piplus_resDelta.Vect()).Dot(y), + (piplus_resDelta.Vect()).Dot(z_hel_Delta) ); + + TVector3 angles_proton_GJ_Delta( (proton_resDelta.Vect()).Dot(x_GJ_Delta), + (proton_resDelta.Vect()).Dot(y), + (proton_resDelta.Vect()).Dot(z_GJ_Delta) ); + + TVector3 angles_piplus_GJ_Delta( (piplus_resDelta.Vect()).Dot(x_GJ_Delta), + (piplus_resDelta.Vect()).Dot(y), + (piplus_resDelta.Vect()).Dot(z_GJ_Delta) ); + + if(frame=="Hel"){ + + userVars[kCosTheta_pim] = angles_piminus_hel_rho.CosTheta(); + userVars[kSinSqTheta_pim] = sin(angles_piminus_hel_rho.Theta())*sin(angles_piminus_hel_rho.Theta()); + userVars[kSin2Theta_pim] = sin(2.*angles_piminus_hel_rho.Theta()); + userVars[kPhi_pim] = angles_piminus_hel_rho.Phi(); + + userVars[kCosTheta_proton] = angles_proton_hel_Delta.CosTheta(); + userVars[kSinSqTheta_proton] = sin(angles_proton_hel_Delta.Theta())*sin(angles_proton_hel_Delta.Theta()); + userVars[kSin2Theta_proton] = sin(2.*angles_proton_hel_Delta.Theta()); + userVars[kPhi_proton] = angles_proton_hel_Delta.Phi(); + } + else{ + userVars[kCosTheta_pim] = angles_piminus_GJ_rho.CosTheta(); + userVars[kSinSqTheta_pim] = sin(angles_piminus_GJ_rho.Theta())*sin(angles_piminus_GJ_rho.Theta()); + userVars[kSin2Theta_pim] = sin(2.*angles_piminus_GJ_rho.Theta()); + userVars[kPhi_pim] = angles_piminus_GJ_rho.Phi(); + + userVars[kCosTheta_proton] = angles_proton_GJ_Delta.CosTheta(); + userVars[kSinSqTheta_proton] = sin(angles_proton_GJ_Delta.Theta())*sin(angles_proton_GJ_Delta.Theta()); + userVars[kSin2Theta_proton] = sin(2.*angles_proton_GJ_Delta.Theta()); + userVars[kPhi_proton] = angles_proton_GJ_Delta.Phi(); + } + TVector3 eps(cos(polAngle*TMath::DegToRad()), sin(polAngle*TMath::DegToRad()), 0.0); + userVars[kBigPhi] = atan2(y.Dot(eps), beam.Vect().Unit().Dot(eps.Cross(y))); + + // vector meson production from K. Schilling et. al. + GDouble Pgamma; // for fitting with constant polarization + Pgamma = polFraction; + // else{ + // int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]); + // if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){ + // Pgamma = 0.; + // } + // else Pgamma = polFrac_vs_E->GetBinContent(bin); + // } + userVars[kPgamma] = Pgamma; + +} + +#ifdef GPU_ACCELERATION +void +TwoPiAngles_Delta_factorized::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const { + + GPUTwoPiAngles_Delta_factorized_exec( dimGrid, dimBlock, GPU_AMP_ARGS, + rho000, + rho100, + rho1m10, + rho001, + rho111, + rho101, + rho1m11, + rho102, + rho1m12, + delta_rho011, + delta_rho031, + delta_rho03m1 +); +} +#endif //GPU_ACCELERATION \ No newline at end of file diff --git a/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h new file mode 100755 index 000000000..3663e48b7 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h @@ -0,0 +1,90 @@ +#if !defined(TWOPIANGLES_DELTA_FACTORIZED) +#define TWOPIANGLES_DELTA_FACTORIZED + +#include "IUAmpTools/Amplitude.h" +#include "IUAmpTools/UserAmplitude.h" +#include "IUAmpTools/AmpParameter.h" +#include "GPUManager/GPUCustomTypes.h" + +#include "TH1D.h" +#include +#include +#include + +#ifdef GPU_ACCELERATION +void +GPUTwoPiAngles_Delta_factorized_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble rho000, + GDouble rho100, + GDouble rho1m10, + GDouble rho001, + GDouble rho111, + GDouble rho101, + GDouble rho1m11, + GDouble rho102, + GDouble rho1m12, + GDouble delta_rho011, + GDouble delta_rho031, + GDouble delta_rho03m1 +); + +#endif // GPU_ACCELERATION + +using std::complex; +using namespace std; + +class Kinematics; + +class TwoPiAngles_Delta_factorized : public UserAmplitude< TwoPiAngles_Delta_factorized > +{ + +public: + + TwoPiAngles_Delta_factorized() : UserAmplitude< TwoPiAngles_Delta_factorized >() { }; + TwoPiAngles_Delta_factorized( const vector< string >& args ); + + enum UserVars { kPgamma = 0, kCosTheta_pim, kSinSqTheta_pim, kSin2Theta_pim, kPhi_pim, kCosTheta_proton, kSinSqTheta_proton, kSin2Theta_proton, kPhi_proton, + kBigPhi, kNumUserVars }; + unsigned int numUserVars() const { return kNumUserVars; } + + string name() const { return "TwoPiAngles_Delta_factorized"; } + + complex< GDouble > calcAmplitude( GDouble** pKin, GDouble* userVars ) const; + void calcUserVars( GDouble** pKin, GDouble* userVars ) const; + + // we can calcualte everythign we need from userVars block so allow + // the framework to purge the four-vectors + bool needsUserVarsOnly() const { return true; } + +#ifdef GPU_ACCELERATION + + void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const; + bool isGPUEnabled() const { return true; } + +#endif // GPU_ACCELERATION + +private: + +AmpParameter rho000; +AmpParameter rho100; +AmpParameter rho1m10; + +AmpParameter rho001; +AmpParameter rho111; +AmpParameter rho101; +AmpParameter rho1m11; + +AmpParameter rho102; +AmpParameter rho1m12; +AmpParameter delta_rho011; +AmpParameter delta_rho031; +AmpParameter delta_rho03m1; + +AmpParameter polAngle; + + double polFraction; + //TH1D *polFrac_vs_E; + string frame; +}; + +#endif diff --git a/src/libraries/AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.cc new file mode 100755 index 000000000..22ddde407 --- /dev/null +++ b/src/libraries/AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.cc @@ -0,0 +1,303 @@ +#include "TLorentzVector.h" +#include "TLorentzRotation.h" +#include "TLorentzVector.h" +#include "TLorentzRotation.h" + +#include "AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.h" +#include "IUAmpTools/Histogram1D.h" +#include "IUAmpTools/Kinematics.h" +#include "IUAmpTools/FitResults.h" + +TwoPiDeltaPlotGenerator::TwoPiDeltaPlotGenerator( const FitResults& results ) : +PlotGenerator( results ) +{ + + vector< string > reactionVec = reactions(); + for( auto reac = reactionVec.begin(); reac != reactionVec.end(); ++reac ){ + + // obtain the polarization angle for this reaction by getting the list of amplitudes + // associated with this reaction -- we know all are SDME amplitudes + // // take the 1th argument of the first factor of the first amplitude in the first sum + // string ampArgument = { cfgInfo()->amplitudeList( *reac, "", "" ).at(0)->factors().at(0).at(1) }; + + // pick out the name of the parameter + // string parName = ampArgument.substr(1,ampArgument.length()-2); + + // // here results is of type FitResults which is passed into the constructor of the plot generator + // //cout << "Angle " << results.parValue( parName ) << endl; + // m_reactionAngleMap[*reac] = results.parValue( parName ); + + double polAngle = 0; + string angleString = cfgInfo()->amplitudeList( *reac, "", "" ).at(0)->factors().at(0).at(2); + + try{ polAngle = stod( angleString ); } + catch( ... ){ + + cout << "ERROR: Unable to get polarization angle from config file because converting\n" + << " the amplitude argument '" << angleString << "' to a double failed.\n" + << " Check the projectEvent method of the OmegaRadiativePlotGenerator to be sure\n" + << " the correct argument is specified. Also be sure the angle appears\n" + << " in the original fit config file as a number and not a parameter.\n" + << endl; + assert( false ); + } + m_reactionAngleMap[*reac] = polAngle; + + + } + + + +for(const auto& entry : m_reactionAngleMap){ + cout << "DEBUG: Reaction in map: " << entry.first + << " with angle: " << entry.second << endl; +} + + createHistograms(); +} + +TwoPiDeltaPlotGenerator::TwoPiDeltaPlotGenerator( ) : +PlotGenerator( ) +{ + createHistograms(); +} + +void TwoPiDeltaPlotGenerator::createHistograms() { + // calls to bookHistogram go here + //signal mass + bookHistogram( k2PiMass, new Histogram1D( 50, 0.6, 0.95, "M2pi", "M(#pi^{0} #pi^{-}) (GeV)") ); + bookHistogram( kPPipMass, new Histogram1D( 50, 1.1, 1.4, "Mppip", "M(p#pi^{+} )(GeV)") ); + //different combinations of invariant masses + bookHistogram( kPipPimMass, new Histogram1D( 100, 0.2, 2.5, "Mpippim", "M(#pi^{+}#pi^{-} )(GeV)") ); + bookHistogram( kPPi0Mass, new Histogram1D( 100, 1, 4.5, "Mppi0", "M(p#pi^{0}) (GeV)") ); + bookHistogram( kPipPi0Mass, new Histogram1D( 100, 0.1, 2.5, "Mpippi0", "M(#pi^{+}#pi^{0}) (GeV)") ); + + bookHistogram( kPPimMass, new Histogram1D( 100, 1, 4.5, "Mppim", "M( p#pi^{-}) (GeV)") ); + bookHistogram( kPi0PimMass, new Histogram1D( 100, 0.1, 1.5, "Mpi0pim", "M(#pi^{0}#pi^{-}) (GeV)") ); + bookHistogram( kPPipPi0Mass, new Histogram1D( 100, 1, 4.5, "Mppippi0", "M( p#pi^{+}#pi^{0}) (GeV)") ); + bookHistogram( kPPimPi0Mass, new Histogram1D( 100, 2.5, 4.5, "Mppimpi0", "M( p#pi^{-}#pi^{0}) (GeV)") ); + bookHistogram( kPipPi0PimMass, new Histogram1D( 100, 0.5, 4.5, "Mpippi0pim", "M( #pi^{-}#pi^{0}#pi^{+}) (GeV)") ); + bookHistogram( kPPipPimMass, new Histogram1D( 100, 1, 4.5, "Mppippim", "M( p#pi^{-}#pi^{+}) (GeV)") ); + + + + bookHistogram( kPhiPiPlus, new Histogram1D( 180, -180, 180, "PhiPiPlus", "#Phi_{#pi^{+}}" ) ); + bookHistogram( kPhiPiMinus, new Histogram1D( 180, -180, 180, "PhiPiMinus", "#Phi_{#pi^{-}}" ) ); + bookHistogram( kPhiProton, new Histogram1D( 180, -180, 180, "PhiProton", "#Phi_{p}" ) ); + + bookHistogram( kThetaPiPlus, new Histogram1D( 200, 0, 20, "ThetaPiPlus", "#Theta_{#pi^{+}}" ) ); + bookHistogram( kThetaPiMinus, new Histogram1D( 200, 0, 20, "ThetaPiMinus", "#Theta_{#pi^{-}}" ) ); + bookHistogram( kThetaDelta, new Histogram1D( 200, 50, 90, "ThetaDelta", "#Theta_{#Delta}" ) ); + + + + + // longitudinal momenta + bookHistogram( klongMomPiPlus, new Histogram1D( 60, 0, 0.5, "longMomPiPlus", "p_{z,#pi^{+}}" ) ); + bookHistogram( klongMomPiMinus, new Histogram1D( 60, 0, 9, "longMomPiMinus", "p_{z,#pi^{-}}" ) ); + bookHistogram( klongMomPi0, new Histogram1D( 60, 0, 9, "longMomPi0", "p_{z,#pi^{0}}" ) ); + bookHistogram( klongMomProton, new Histogram1D( 60, 0, 1.5, "longMomProton", "p_{z,p}" ) ); + + // momenta + bookHistogram( kMomPiPlus, new Histogram1D( 60, 0, 0.7, "MomPiPlus", "p_{#pi^{+}}" ) ); + bookHistogram( kMomPiMinus, new Histogram1D( 60, 0, 9, "MomPiMinus", "p_{#pi^{-}}" ) ); + bookHistogram( kMomPi0, new Histogram1D( 60, 0, 9, "MomPi0", "p_{#pi^{0}}" ) ); + bookHistogram( kMomProton, new Histogram1D( 180, 0, 1, "MomProton", "p_{p}" ) ); + + // Angles in different frames Phi + bookHistogram( kphi_PiMinus_hel_rho, new Histogram1D( 60, -180, 180, "phi_PiMinus_hel_rho", " #phi_{#pi^{-}} in hel (restframe #rho) " ) ); + bookHistogram( kphi_PiMinus_GJ_rho, new Histogram1D( 60, -180, 180, "phi_PiMinus_GJ_rho", "#phi_{#pi^{-}} in GJ (restframe #rho)" ) ); + + bookHistogram( kphi_PiPlus_hel_Delta, new Histogram1D( 60, -180, 180, "phi_PiPlus_hel_Delta", " #phi_{#pi^{+}} in hel (restframe #Delta)" ) ); + bookHistogram( kphi_PiPlus_GJ_Delta, new Histogram1D( 60, -180, 180, "phi_PiPlus_GJ_Delta", "#phi_{#pi^{+}} in GJ (restframe #Delta)" ) ); + + bookHistogram( kphi_Proton_hel_Delta, new Histogram1D( 60, -180, 180, "phi_Proton_hel_Delta", " #phi_{p} in hel (restframe #Delta)" ) ); + bookHistogram( kphi_Proton_GJ_Delta, new Histogram1D( 60, -180, 180, "phi_Proton_GJ_Delta", "#phi_{p} in GJ (restframe #Delta)" ) ); + + bookHistogram( kCosTheta_PiMinus_hel_rho, new Histogram1D( 60, -1., 1., "cosTheta_PiMinus_hel_rho", "cos( #theta_{#pi^{-}} ) in hel (restframe #rho)") ); + bookHistogram( kCosTheta_PiMinus_GJ_rho, new Histogram1D( 60, -1., 1., "cosTheta_PiMinus_GJ_rho", "cos( #theta_{#pi^{-}} ) in GJ (restframe #rho)") ); + + bookHistogram( kCosTheta_PiPlus_hel_Delta, new Histogram1D( 60, -1., 1., "cosTheta_PiPlus_hel_Delta", "cos( #theta_{#pi^{+}} ) in hel (restframe #Delta)") ); + bookHistogram( kCosTheta_PiPlus_GJ_Delta, new Histogram1D( 60, -1., 1., "cosTheta_PiPlus_GJ_Delta", "cos( #theta_{#pi^{+}} ) in GJ (restframe #Delta)") ); + + bookHistogram( kCosTheta_Proton_hel_Delta, new Histogram1D( 60, -1., 1., "cosTheta_Proton_hel_Delta", "cos( #theta_{p} ) in hel (restframe #Delta)") ); + bookHistogram( kCosTheta_Proton_GJ_Delta, new Histogram1D( 60, -1., 1., "cosTheta_Proton_GJ_Delta", "cos( #theta_{p} ) in GJ (restframe #Delta)") ); + + + bookHistogram( kPhi, new Histogram1D( 50, -180, 180, "Phi", "#Phi" ) ); + bookHistogram( kPsi, new Histogram1D( 20, -180, 180, "psi", "#psi" ) ); + bookHistogram( kt, new Histogram1D( 1000, 0, 1.4, "-t", "-t (GeV^2)" ) ); + +} + +void TwoPiDeltaPlotGenerator::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 + + projectEvent( kin, "" ); +} + +void TwoPiDeltaPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){ + + double polAngle = m_reactionAngleMap[ reactionName ]; + + TLorentzVector beam = kin->particle( 0 ); + TLorentzVector proton = kin->particle( 1 ); //proton + TLorentzVector p1 = kin->particle( 2 ); //pip + TLorentzVector p2 = kin->particle( 3 ); //pim + TLorentzVector p3 = kin->particle( 4 ); //pi0 + TLorentzVector target ( 0, 0, 0, 0.9382720813); + + TLorentzVector recoil = proton + p1; + TLorentzVector rho = p2 + p3; + TLorentzVector Delta = proton + p1; + TLorentzVector CM_motion_lab = beam + target; + TLorentzVector piminus = p2; + + TLorentzRotation rhoBoost( -rho.BoostVector() ); + TLorentzRotation CMBoost( -CM_motion_lab.BoostVector() ); + TLorentzRotation DeltaBoost( -Delta.BoostVector() ); + + +// choose helicity frame: z-axis opposite recoil Delta in rho rest frame +//bost in rho restframe + TLorentzVector beam_resrho = rhoBoost * beam; + TLorentzVector Delta_resrho = rhoBoost * Delta; + TLorentzVector piminus_resrho = rhoBoost * piminus; + TLorentzVector proton_resrho = rhoBoost * proton; +//boost in CM frame + TLorentzVector beam_cm = CMBoost * beam; + TLorentzVector Delta_cm = CMBoost * Delta; + TLorentzVector piminus_cm = CMBoost * piminus; + TLorentzVector rho_cm = CMBoost * rho; + +//boost in recoil frame + TLorentzVector target_resDelta = DeltaBoost * target; + TLorentzVector proton_resDelta = DeltaBoost * proton; + TLorentzVector piplus_resDelta = DeltaBoost * p1; + TLorentzVector rho_resDelta = DeltaBoost * rho; + + // Define coordinate systems: + // y equal in all frames + TVector3 y = (beam.Vect().Unit().Cross(-Delta.Vect().Unit())).Unit(); + + //rho system + + //Helicity frame + TVector3 z_hel_rho = -1. * Delta_resrho.Vect().Unit(); + TVector3 x_hel_rho = y.Cross(z_hel_rho).Unit(); + + // GJ frame + + TVector3 z_GJ_rho = beam_resrho.Vect().Unit(); + TVector3 x_GJ_rho = y.Cross(z_GJ_rho).Unit(); + + //Delta system + + // Helicity frame + TVector3 z_hel_Delta = -rho_resDelta.Vect().Unit(); + TVector3 x_hel_Delta = (y.Cross(z_hel_Delta)).Unit(); + + // GJ frame + TVector3 z_GJ_Delta = target_resDelta.Vect().Unit(); + TVector3 x_GJ_Delta = y.Cross(z_GJ_Delta).Unit(); + + //Angle calculation + TVector3 angles_piminus_hel_rho( (piminus_resrho.Vect()).Dot(x_hel_rho), + (piminus_resrho.Vect()).Dot(y), + (piminus_resrho.Vect()).Dot(z_hel_rho) ); + + TVector3 angles_piminus_GJ_rho( (piminus_resrho.Vect()).Dot(x_GJ_rho), + (piminus_resrho.Vect()).Dot(y), + (piminus_resrho.Vect()).Dot(z_GJ_rho) ); + + + TVector3 angles_proton_hel_Delta( (proton_resDelta.Vect()).Dot(x_hel_Delta), + (proton_resDelta.Vect()).Dot(y), + (proton_resDelta.Vect()).Dot(z_hel_Delta) ); + TVector3 angles_piplus_hel_Delta( (piplus_resDelta.Vect()).Dot(x_hel_Delta), + (piplus_resDelta.Vect()).Dot(y), + (piplus_resDelta.Vect()).Dot(z_hel_Delta) ); + + TVector3 angles_proton_GJ_Delta( (proton_resDelta.Vect()).Dot(x_GJ_Delta), + (proton_resDelta.Vect()).Dot(y), + (proton_resDelta.Vect()).Dot(z_GJ_Delta) ); + + TVector3 angles_piplus_GJ_Delta( (piplus_resDelta.Vect()).Dot(x_GJ_Delta), + (piplus_resDelta.Vect()).Dot(y), + (piplus_resDelta.Vect()).Dot(z_GJ_Delta) ); + + TVector3 eps(cos(polAngle*TMath::DegToRad()), sin(polAngle*TMath::DegToRad()), 0.0); // beam polarization vector + GDouble Phi = atan2(y.Dot(eps), beam.Vect().Unit().Dot(eps.Cross(y))); + + GDouble psi = angles_piminus_GJ_rho.Phi() - Phi; + if(psi < -1*PI) psi += 2*PI; + if(psi > PI) psi -= 2*PI; + + // compute invariant t + GDouble t = -(beam-p3-p2).M2(); + + // calls to fillHistogram go here + fillHistogram( k2PiMass, ( rho ).M() ); + fillHistogram( kPPipMass, ( proton+p1 ).M() ); + fillHistogram( kPPimMass, ( proton+p2 ).M() ); + fillHistogram( kPi0PimMass, ( p3+p2 ).M() ); + + + fillHistogram( kPipPimMass, (p1+p2 ).M()); + fillHistogram( kPPi0Mass, (proton+p3 ).M()); + fillHistogram( kPipPi0Mass, (p1+p3 ).M()); + + ; + fillHistogram( klongMomPiPlus, p1.Pz() ); + fillHistogram( klongMomPiMinus, p2.Pz() ); + fillHistogram( klongMomPi0, p3.Pz() ); + fillHistogram( klongMomProton, proton.Pz() ); + + fillHistogram( kMomPiPlus, p1.P() ); + fillHistogram( kMomPiMinus, p2.P() ); + fillHistogram( kMomPi0, p3.P() ); + fillHistogram( kMomProton, proton.P() ); + fillHistogram( kPPipPi0Mass, ( proton + p1 + p3 ).M() ); + fillHistogram( kPipPi0PimMass, ( p1 + p3 + p2 ).M() ); + fillHistogram( kPPipPimMass, ( proton + p1 + p2 ).M() ); + fillHistogram( kPPimPi0Mass, ( proton + p2 + p3 ).M() ); + + + //angles for fits + + // phi angles in different frames + fillHistogram( kphi_PiMinus_hel_rho, angles_piminus_hel_rho.Phi()*TMath::RadToDeg() ); + fillHistogram( kphi_PiMinus_GJ_rho, angles_piminus_GJ_rho.Phi()*TMath::RadToDeg()); + + fillHistogram( kphi_PiPlus_hel_Delta, angles_piplus_hel_Delta.Phi()*TMath::RadToDeg()); + fillHistogram( kphi_PiPlus_GJ_Delta, angles_piplus_GJ_Delta.Phi()*TMath::RadToDeg()); + + fillHistogram( kphi_Proton_hel_Delta, angles_proton_hel_Delta.Phi()*TMath::RadToDeg()); + fillHistogram( kphi_Proton_GJ_Delta, angles_proton_GJ_Delta.Phi()*TMath::RadToDeg()); + + // cosTheta angles in different frames + fillHistogram( kCosTheta_PiMinus_hel_rho, angles_piminus_hel_rho.CosTheta() ); + fillHistogram( kCosTheta_PiMinus_GJ_rho, angles_piminus_GJ_rho.CosTheta() ); + + fillHistogram( kCosTheta_PiPlus_hel_Delta, angles_piplus_hel_Delta.CosTheta() ); + fillHistogram( kCosTheta_PiPlus_GJ_Delta, angles_piplus_GJ_Delta.CosTheta() ); + + fillHistogram( kCosTheta_Proton_hel_Delta, angles_proton_hel_Delta.CosTheta() ); + fillHistogram( kCosTheta_Proton_GJ_Delta, angles_proton_GJ_Delta.CosTheta() ); + + //Different angles in lab frame + fillHistogram( kThetaPiPlus, p1.Theta()*TMath::RadToDeg() ); + fillHistogram( kThetaPiMinus, p2.Theta()*TMath::RadToDeg() ); + fillHistogram( kThetaDelta, recoil.Theta()*TMath::RadToDeg() ); + fillHistogram( kPhiPiPlus, p1.Phi()*TMath::RadToDeg() ); + fillHistogram( kPhiPiMinus, p2.Phi()*TMath::RadToDeg() ); + fillHistogram( kPhiProton, proton.Phi()*TMath::RadToDeg() ); + + //Polarisation angles + fillHistogram( kPhi, Phi*TMath::RadToDeg() ); + fillHistogram( kPsi, psi*TMath::RadToDeg() ); + fillHistogram( kt, t ); +} diff --git a/src/libraries/AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.h b/src/libraries/AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.h new file mode 100755 index 000000000..a97878a14 --- /dev/null +++ b/src/libraries/AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.h @@ -0,0 +1,77 @@ +#if !(defined TWOPIDELTAPLOTGENERATOR) +#define TWOPIDELTAPLOTGENERATOR + +#include +#include + +#include "IUAmpTools/PlotGenerator.h" + +using namespace std; + +class FitResults; +class Kinematics; + +class TwoPiDeltaPlotGenerator : public PlotGenerator +{ + +public: + + // create an index for different histograms + enum { + k2PiMass = 0, + kPPipMass, + kPPimMass, + kPi0PimMass, + kPPipPimMass, + kPPipPi0Mass, + kPipPi0PimMass, + kPPimPi0Mass, + kPipPimMass, + kPPi0Mass, + kPipPi0Mass, + kThetaDelta, + kThetaPiPlus, + kThetaPiMinus, + kMomPiPlus, + kMomPiMinus, + klongMomPiPlus, + klongMomPiMinus, + klongMomPi0, + klongMomProton, + kMomPi0, + kMomProton, + kPhiPiPlus, + kPhiPiMinus, + kPhiProton, + kphi_PiMinus_hel_rho, + kphi_PiMinus_GJ_rho, + kphi_PiPlus_hel_Delta, + kphi_PiPlus_GJ_Delta, + kphi_Proton_hel_Delta, + kphi_Proton_GJ_Delta, + kCosTheta_PiMinus_hel_rho, + kCosTheta_PiMinus_GJ_rho, + kCosTheta_PiPlus_hel_Delta, + kCosTheta_PiPlus_GJ_Delta, + kCosTheta_Proton_hel_Delta, + kCosTheta_Proton_GJ_Delta, + kPhi, + kPsi, + kt, + kNumHists + }; + + TwoPiDeltaPlotGenerator( const FitResults& results ); + TwoPiDeltaPlotGenerator( ); + + void projectEvent( Kinematics* kin ); + void projectEvent( Kinematics* kin, const string& reactionName ); + +private: + + void createHistograms(); + + map< string, double > m_reactionAngleMap; +}; + +#endif diff --git a/src/programs/AmplitudeAnalysis/SConscript b/src/programs/AmplitudeAnalysis/SConscript index 3ca97811b..043887f00 100644 --- a/src/programs/AmplitudeAnalysis/SConscript +++ b/src/programs/AmplitudeAnalysis/SConscript @@ -3,7 +3,8 @@ import sbms Import('*') -subdirs = ['fit', 'twopi_plotter', 'twopi_plotter_amp', 'twopi_plotter_mom', 'twopi_plotter_primakoff', 'twolepton_plotter', 'twoleptonGJ_plotter', 'split_mass', 'split_t', 'threepi_plotter_schilling', 'omega_radiative_plotter', 'project_moments', 'plot_etapi_delta', 'project_moments_polarized', 'Bootstrap_plot_etapi_delta_SPDG_allamps_mass_t_bins', 'Pol_moments_viafittedPW', 'project_moments_SPD_etapi0_posepsilon', 'omegapi_plotter', 'vecps_plotter', 'plot_etapi0', 'isops_plotter', 'convert_to_csv'] + +subdirs = ['fit', 'twopi_plotter', 'twopi_plotter_amp', 'twopi_plotter_mom', 'twopi_plotter_primakoff', 'twolepton_plotter', 'twoleptonGJ_plotter', 'split_mass', 'split_t', 'threepi_plotter_schilling', 'omega_radiative_plotter', 'project_moments', 'plot_etapi_delta', 'project_moments_polarized', 'Bootstrap_plot_etapi_delta_SPDG_allamps_mass_t_bins', 'Pol_moments_viafittedPW', 'project_moments_SPD_etapi0_posepsilon', 'omegapi_plotter', 'vecps_plotter', 'plot_etapi0', 'isops_plotter', 'convert_to_csv','twopidelta_plotter'] SConscript(dirs=subdirs, exports='env osname', duplicate=0) diff --git a/src/programs/AmplitudeAnalysis/fit/fit.cc b/src/programs/AmplitudeAnalysis/fit/fit.cc index b723de23d..4d7484674 100644 --- a/src/programs/AmplitudeAnalysis/fit/fit.cc +++ b/src/programs/AmplitudeAnalysis/fit/fit.cc @@ -22,6 +22,8 @@ #include "AMPTOOLS_AMPS/TwoPSHelicity.h" #include "AMPTOOLS_AMPS/TwoPiAngles.h" #include "AMPTOOLS_AMPS/TwoPiAngles_amp.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h" #include "AMPTOOLS_AMPS/TwoPiWt_primakoff.h" #include "AMPTOOLS_AMPS/TwoPiWt_sigma.h" #include "AMPTOOLS_AMPS/TwoPiW_brokenetas.h" @@ -439,6 +441,8 @@ int main( int argc, char* argv[] ){ AmpToolsInterface::registerAmplitude( TwoPSAngles() ); AmpToolsInterface::registerAmplitude( TwoPSHelicity() ); AmpToolsInterface::registerAmplitude( TwoPiAngles() ); + AmpToolsInterface::registerAmplitude( TwoPiAngles_Delta_DoubleSDMEs()); + AmpToolsInterface::registerAmplitude( TwoPiAngles_Delta_factorized()); AmpToolsInterface::registerAmplitude( TwoPiAngles_amp() ); AmpToolsInterface::registerAmplitude( TwoPiAngles_primakoff() ); AmpToolsInterface::registerAmplitude( TwoPiWt_primakoff() ); diff --git a/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc b/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc index 7656b0703..c60965657 100644 --- a/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc +++ b/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc @@ -20,6 +20,8 @@ #include "AMPTOOLS_AMPS/TwoPSHelicity.h" #include "AMPTOOLS_AMPS/TwoPiAngles.h" #include "AMPTOOLS_AMPS/TwoPiAngles_amp.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h" #include "AMPTOOLS_AMPS/TwoPiWt_primakoff.h" #include "AMPTOOLS_AMPS/TwoPiWt_sigma.h" #include "AMPTOOLS_AMPS/TwoPiW_brokenetas.h" @@ -414,6 +416,8 @@ int main( int argc, char* argv[] ){ AmpToolsInterface::registerAmplitude( TwoPSAngles() ); AmpToolsInterface::registerAmplitude( TwoPSHelicity() ); AmpToolsInterface::registerAmplitude( TwoPiAngles() ); + AmpToolsInterface::registerAmplitude( TwoPiAngles_Delta_DoubleSDMEs()); + AmpToolsInterface::registerAmplitude( TwoPiAngles_Delta_factorized()); AmpToolsInterface::registerAmplitude( TwoPiAngles_amp() ); AmpToolsInterface::registerAmplitude( TwoPiAngles_primakoff() ); AmpToolsInterface::registerAmplitude( TwoPiWt_primakoff() ); diff --git a/src/programs/AmplitudeAnalysis/twopidelta_plotter/SConscript b/src/programs/AmplitudeAnalysis/twopidelta_plotter/SConscript new file mode 100755 index 000000000..642f372be --- /dev/null +++ b/src/programs/AmplitudeAnalysis/twopidelta_plotter/SConscript @@ -0,0 +1,23 @@ + +import os +import sbms + +# get env object and clone it +Import('*') + +# Verify AMPTOOLS environment variable is set +if os.getenv('AMPTOOLS', 'nada')!='nada' and os.getenv('AMPPLOTTER', 'nada')!='nada': + + env = env.Clone() + + AMPTOOLS_LIBS = "AMPTOOLS_AMPS AMPTOOLS_DATAIO AMPTOOLS_MCGEN UTILITIES" + env.AppendUnique(LIBS = AMPTOOLS_LIBS.split()) + + sbms.AddHDDM(env) + sbms.AddAmpTools(env) + sbms.AddAmpPlotter(env) + sbms.AddUtilities(env) + sbms.AddROOT(env) + + sbms.executable(env) + diff --git a/src/programs/AmplitudeAnalysis/twopidelta_plotter/twopidelta_plotter.cc b/src/programs/AmplitudeAnalysis/twopidelta_plotter/twopidelta_plotter.cc new file mode 100755 index 000000000..76bb35a94 --- /dev/null +++ b/src/programs/AmplitudeAnalysis/twopidelta_plotter/twopidelta_plotter.cc @@ -0,0 +1,843 @@ +#include +#include +#include +#include + +#include "TClass.h" +#include "TApplication.h" +#include "TGClient.h" +#include "TROOT.h" +#include "TH1.h" +#include "TStyle.h" +#include "TClass.h" +#include "TFile.h" + +#include "IUAmpTools/AmpToolsInterface.h" +#include "IUAmpTools/FitResults.h" + +#include "AmpPlotter/PlotterMainWindow.h" +#include "AmpPlotter/PlotFactory.h" + +#include "AMPTOOLS_DATAIO/TwoPiDeltaPlotGenerator.h" +#include "AMPTOOLS_DATAIO/ROOTDataReader.h" +#include "AMPTOOLS_DATAIO/ROOTDataReaderHist.h" +#include "AMPTOOLS_AMPS/Zlm.h" +#include "AMPTOOLS_AMPS/TwoPiAngles.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_DoubleSDMEs.h" +#include "AMPTOOLS_AMPS/TwoPiAngles_Delta_factorized.h" +#include "AMPTOOLS_AMPS/BreitWigner.h" +#include "AMPTOOLS_AMPS/Uniform.h" + +typedef TwoPiDeltaPlotGenerator PlotGen; + +void atiSetup(){ + + AmpToolsInterface::registerAmplitude( Zlm() ); + + AmpToolsInterface::registerAmplitude( TwoPiAngles_Delta_DoubleSDMEs() ); + AmpToolsInterface::registerAmplitude( TwoPiAngles_Delta_factorized() ); + AmpToolsInterface::registerAmplitude( BreitWigner() ); + AmpToolsInterface::registerAmplitude( Uniform() ); + AmpToolsInterface::registerDataReader( ROOTDataReader() ); + AmpToolsInterface::registerDataReader( ROOTDataReaderHist() ); +} + +using namespace std; + +int main( int argc, char* argv[] ){ //char array for strings + + + // ************************ + // usage + // ************************ + + cout << endl << " *** Viewing Results Using AmpPlotter and writing root histograms *** " << endl << endl; + + if (argc < 2){ + cout << "Usage:" << endl << endl; + cout << "\ttwopidelta_plotter -o " << endl << endl; + return 0; + } + + bool showGui = false; + bool rootfile_creation = true; + string outName = "twopidelta_plot.root"; + string resultsName(argv[1]); + Double_t tmean; + Double_t terr; + double tmean_old = 0; + double terr_old = 0; + for (int i = 2; i < argc; i++){ + + string arg(argv[i]); + + if (arg == "-g"){ + showGui = true; + } + if (arg == "-o"){ + outName = argv[++i]; + } + if (arg == "-h"){ + cout << endl << " Usage for: " << argv[0] << endl << endl; + cout << "\t -o \t output file path" << endl; + cout << "\t -g \t show GUI" << endl; + exit(1); + } + if (arg == "-n"){ + rootfile_creation = false; + std::ifstream infile("twopidelta_fitPars.txt"); + std::string name; + double value, error; + while (infile >> name >> value >> error) { + if (name == "tmean") { + tmean_old = value; + terr_old = error; + break; + } + } +} + + } + + + // ************************ + // parse the command line parameters + // ************************ + + cout << "Fit results file name = " << resultsName << endl; + cout << "Output file name = " << outName << endl << endl; + + // ************************ + // load the results and display the configuration info + // ************************ + + FitResults results( resultsName ); //get fitresults + cout << "Fit results loaded" << endl; + if( !results.valid() ){ + + cout << "Invalid fit results in file: " << resultsName << endl; + exit( 1 ); + } + + // ************************ + // set up the plot generator + // ************************ + cout << "Setting up ATI" << endl; + atiSetup(); + cout << "ATI set up" << endl; + PlotGen plotGen( results ); + cout << "Plot generator set up" << endl; + + // ************************ + // set up an output ROOT file to store histograms + // ************************ + if (rootfile_creation){ + TFile* plotfile = new TFile( outName.c_str(), "recreate"); + TH1::AddDirectory(kFALSE); + + string reactionName = results.reactionList()[0]; + auto reactions = results.reactionList(); + cout << "Number of reactions: " << reactions.size() << endl; + for (size_t k = 0; k < reactions.size(); ++k) { + auto& reaction = reactions[k]; + string reactionName = reaction; + cout << "Reaction name: " << reactionName << endl; + plotGen.enableReaction( reactionName ); + vector sums = plotGen.uniqueSums(); + + + // loop over sum configurations (one for each of the individual contributions, and the combined sum of all) + + for (unsigned int isum = 0; isum <= sums.size(); isum++){ + + // turn on all sums by default + for (unsigned int i = 0; i < sums.size(); i++){ + plotGen.enableSum(i); + } + + // for individual contributions turn off all sums but the one of interest + if (isum < sums.size()){ + for (unsigned int i = 0; i < sums.size(); i++){ + if (i != isum) plotGen.disableSum(i); + } + } + + + // loop over data, accMC, and genMC + for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){ + if (isum < sums.size() && iplot == PlotGenerator::kData) continue; // only plot data once + + // loop over different variables + for (unsigned int ivar = 0; ivar < TwoPiDeltaPlotGenerator::kNumHists; ivar++){ + + // set unique histogram name for each plot (could put in directories...) + string histname = ""; + if (ivar == TwoPiDeltaPlotGenerator::k2PiMass) histname += "M2pi"; + else if (ivar == TwoPiDeltaPlotGenerator::kPPipMass) histname += "Mppip"; + else if (ivar == TwoPiDeltaPlotGenerator::kPPimMass) histname += "Mppim"; + else if (ivar == TwoPiDeltaPlotGenerator::kPi0PimMass) histname += "Mpi0pim"; + else if (ivar == TwoPiDeltaPlotGenerator::kPPipPimMass) histname += "Mppippim"; + else if (ivar == TwoPiDeltaPlotGenerator::kPPipPi0Mass) histname += "Mppippi0"; + else if (ivar == TwoPiDeltaPlotGenerator::kPPimPi0Mass) histname += "Mppimpi0"; + else if (ivar == TwoPiDeltaPlotGenerator::kPipPi0PimMass) histname += "Mpippi0pim"; + else if (ivar == TwoPiDeltaPlotGenerator::kPipPimMass) histname += "Mpippim"; + else if (ivar == TwoPiDeltaPlotGenerator::kPPi0Mass) histname += "Mppi0"; + else if (ivar == TwoPiDeltaPlotGenerator::kPipPi0Mass) histname += "Mpippi0"; + + else if (ivar == TwoPiDeltaPlotGenerator::kPhiPiPlus) histname += "PhiPiPlus"; + else if (ivar == TwoPiDeltaPlotGenerator::kPhiPiMinus) histname += "PhiPiMinus"; + else if (ivar == TwoPiDeltaPlotGenerator::kPhiProton) histname += "PhiProton"; + + else if (ivar == TwoPiDeltaPlotGenerator::kThetaPiPlus) histname += "ThetaPiPlus"; + else if (ivar == TwoPiDeltaPlotGenerator::kThetaPiMinus) histname += "ThetaPiMinus"; + else if (ivar == TwoPiDeltaPlotGenerator::kThetaDelta) histname += "ThetaDelta"; + + else if (ivar == TwoPiDeltaPlotGenerator::klongMomPiPlus) histname += "longMomPiPlus"; + else if (ivar == TwoPiDeltaPlotGenerator::klongMomPiMinus) histname += "longMomPiMinus"; + else if (ivar == TwoPiDeltaPlotGenerator::klongMomPi0) histname += "longMomPi0"; + else if (ivar == TwoPiDeltaPlotGenerator::klongMomProton) histname += "longMomProton"; + + else if (ivar == TwoPiDeltaPlotGenerator::kMomPiPlus) histname += "MomPiPlus"; + else if (ivar == TwoPiDeltaPlotGenerator::kMomPiMinus) histname += "MomPiMinus"; + else if (ivar == TwoPiDeltaPlotGenerator::kMomPi0) histname += "MomPi0"; + else if (ivar == TwoPiDeltaPlotGenerator::kMomProton) histname += "MomProton"; + + else if (ivar == TwoPiDeltaPlotGenerator::kphi_PiMinus_hel_rho) histname += "phi_PiMinus_hel_rho"; + else if (ivar == TwoPiDeltaPlotGenerator::kphi_PiMinus_GJ_rho) histname += "phi_PiMinus_GJ_rho"; + else if (ivar == TwoPiDeltaPlotGenerator::kphi_PiPlus_hel_Delta) histname += "phi_PiPlus_hel_Delta"; + else if (ivar == TwoPiDeltaPlotGenerator::kphi_PiPlus_GJ_Delta) histname += "phi_PiPlus_GJ_Delta"; + else if (ivar == TwoPiDeltaPlotGenerator::kphi_Proton_hel_Delta) histname += "phi_Proton_hel_Delta"; + else if (ivar == TwoPiDeltaPlotGenerator::kphi_Proton_GJ_Delta) histname += "phi_Proton_GJ_Delta"; + + else if (ivar == TwoPiDeltaPlotGenerator::kCosTheta_PiMinus_hel_rho) histname += "cosTheta_PiMinus_hel_rho"; + else if (ivar == TwoPiDeltaPlotGenerator::kCosTheta_PiMinus_GJ_rho) histname += "cosTheta_PiMinus_GJ_rho"; + else if (ivar == TwoPiDeltaPlotGenerator::kCosTheta_PiPlus_hel_Delta) histname += "cosTheta_PiPlus_hel_Delta"; + else if (ivar == TwoPiDeltaPlotGenerator::kCosTheta_PiPlus_GJ_Delta) histname += "cosTheta_PiPlus_GJ_Delta"; + else if (ivar == TwoPiDeltaPlotGenerator::kCosTheta_Proton_hel_Delta) histname += "cosTheta_Proton_hel_Delta"; + else if (ivar == TwoPiDeltaPlotGenerator::kCosTheta_Proton_GJ_Delta) histname += "cosTheta_Proton_GJ_Delta"; + + else if (ivar == TwoPiDeltaPlotGenerator::kPhi) histname += "Phi"; + else if (ivar == TwoPiDeltaPlotGenerator::kPsi) histname += "psi"; + else if (ivar == TwoPiDeltaPlotGenerator::kt) histname += "-t"; + + else continue; + + if (iplot == PlotGenerator::kData) histname += "dat"; + if (iplot == PlotGenerator::kAccMC) histname += "acc"; + if (iplot == PlotGenerator::kGenMC) histname += "gen"; + if (iplot == PlotGenerator::kBkgnd) histname += "_bkgnd"; + + if (isum < sums.size()){ + //ostringstream sdig; sdig << (isum + 1); + //histname += sdig.str(); + + // get name of sum for naming histogram + string sumName = sums[isum]; + histname += "_"; + histname += sumName; + } + ; + if (k==0){ + Histogram* hist = plotGen.projection(ivar, reactionName, iplot); + TH1* thist = hist->toRoot(); + thist->SetName(histname.c_str()); + + if (histname=="-tdat"&& k==reactions.size()-1){ + tmean = thist->GetMean(); + terr = thist->GetRMS(); + } + plotfile->cd(); + thist->Write(); + } + else + { + Histogram* hist = plotGen.projection(ivar, reactionName, iplot); + TH1* histSumReac = (TH1*) plotfile->Get(histname.c_str()); + TH1* thist = hist->toRoot(); + histSumReac->Add(thist); + if (histname=="-tdat" && k==reactions.size()-1){ + tmean = histSumReac->GetMean(); + terr = histSumReac->GetRMS(); + } + plotfile->cd(); + histSumReac->Write(histname.c_str(), TObject::kOverwrite); + + } + + // histname += "_"; + // histname += reactionName; + // Histogram* hist = plotGen.projection(ivar, reactionName, iplot); + // TH1* thist = hist->toRoot(); + // thist->SetName(histname.c_str()); + // plotfile->cd(); + // thist->Write(); + + } + } + } +} + plotfile->Close(); + } + // ************************ + // retrieve SDME parameters for plotting and asymmetry + // ************************ + + // parameters to check! be careful! for rho last number is upper case # Delta first number + vector< string > pars; + pars.push_back("rho000"); + pars.push_back("rho100"); + pars.push_back("rho1m10"); + + pars.push_back("rho111"); + pars.push_back("rho001"); + pars.push_back("rho101"); + pars.push_back("rho1m11"); + + pars.push_back("rho102"); + pars.push_back("rho1m12"); + pars.push_back("delta_rho011"); + pars.push_back("delta_rho031"); + pars.push_back("delta_rho03m1"); + pars.push_back("delta_rho111"); + pars.push_back("delta_rho133"); + pars.push_back("delta_rho131"); + pars.push_back("delta_rho13m1"); + pars.push_back("delta_rho231"); + pars.push_back("delta_rho23m1"); + pars.push_back("r00_33_0"); + pars.push_back("r00_11_0"); + pars.push_back("r11_33_0"); + pars.push_back("r11_11_0"); + pars.push_back("r1m1_33_0"); + pars.push_back("r1m1_11_0"); + pars.push_back("r10_33_0"); + pars.push_back("r10_11_0"); + pars.push_back("r00_31_0"); + pars.push_back("r00_3m1_0"); + pars.push_back("r11_31_0"); + pars.push_back("r11_3m1_0"); + pars.push_back("r1m1_31_0"); + pars.push_back("r1m1_3m1_0"); + pars.push_back("r10_31_0"); + pars.push_back("r10_3m1_0"); + pars.push_back("rt1m1_31_0"); + pars.push_back("rt1m1_3m1_0"); + pars.push_back("rt10_31_0"); + pars.push_back("rt10_3m1_0"); + // alpha = 1 + pars.push_back("r00_33_1"); + pars.push_back("r00_11_1"); + pars.push_back("r11_33_1"); + pars.push_back("r11_11_1"); + pars.push_back("r1m1_33_1"); + pars.push_back("r1m1_11_1"); + pars.push_back("r10_33_1"); + pars.push_back("r10_11_1"); + pars.push_back("r00_31_1"); + pars.push_back("r00_3m1_1"); + pars.push_back("r11_31_1"); + pars.push_back("r11_3m1_1"); + pars.push_back("r1m1_31_1"); + pars.push_back("r1m1_3m1_1"); + pars.push_back("r10_31_1"); + pars.push_back("r10_3m1_1"); + pars.push_back("rt1m1_31_1"); + pars.push_back("rt1m1_3m1_1"); + pars.push_back("rt10_31_1"); + pars.push_back("rt10_3m1_1"); + + // alpha = 2 + pars.push_back("rt00_31_2"); + pars.push_back("rt00_3m1_2"); + pars.push_back("rt11_31_2"); + pars.push_back("rt11_3m1_2"); + pars.push_back("r1m1_33_2"); + pars.push_back("r1m1_11_2"); + pars.push_back("r1m1_31_2"); + pars.push_back("r1m1_3m1_2"); + pars.push_back("rt1m1_31_2"); + pars.push_back("rt1m1_3m1_2"); + pars.push_back("r10_33_2"); + pars.push_back("r10_11_2"); + pars.push_back("r10_31_2"); + pars.push_back("r10_3m1_2"); + pars.push_back("rt10_31_2"); + pars.push_back("rt10_3m1_2"); + + + + + // file for writing parameters (later switch to putting in ROOT file) + ofstream outfile; + outfile.open( "twopidelta_fitPars.txt" ); + if (rootfile_creation){ + outfile << "tmean" << "\t" << tmean << "\t" << terr << "\t" << "\n"; + } + else{ + outfile << "tmean" << "\t" << tmean_old << "\t" << terr_old << "\t" << "\n"; + } + for(unsigned int i = 0; i > covMatrix; + covMatrix = results.errorMatrix(); + outfile << endl; + + // ************************ + // start the GUI + // ************************ + + if(showGui) { + + cout << ">> Plot generator ready, starting GUI..." << endl; + + int dummy_argc = 0; + char* dummy_argv[] = {}; + TApplication app( "app", &dummy_argc, dummy_argv ); + + gStyle->SetFillColor(10); + gStyle->SetCanvasColor(10); + gStyle->SetPadColor(10); + gStyle->SetFillStyle(1001); + gStyle->SetPalette(1); + gStyle->SetFrameFillColor(10); + gStyle->SetFrameFillStyle(1001); + + PlotFactory factory( plotGen ); + PlotterMainWindow mainFrame( gClient->GetRoot(), factory ); + + app.Run(); + } + + return 0; + +} +