Skip to content

Commit

Permalink
Merge pull request #173 from JeffersonLab/unstableLowerVertex
Browse files Browse the repository at this point in the history
Unstable lower vertex
  • Loading branch information
gleasonc authored Dec 3, 2020
2 parents 830a6f9 + 67460d4 commit 1f94234
Show file tree
Hide file tree
Showing 15 changed files with 1,139 additions and 111 deletions.
79 changes: 31 additions & 48 deletions src/libraries/AMPTOOLS_AMPS/omegapi_amplitude.cc
Original file line number Diff line number Diff line change
Expand Up @@ -30,16 +30,17 @@
omegapi_amplitude::omegapi_amplitude( const vector< string >& args ):
UserAmplitude< omegapi_amplitude >( args )
{
assert( args.size() == (7+4+2) || args.size() == (7+4+3) );
assert( args.size() == (6+4+2) || args.size() == (6+4+3) );

if(args.size() == (7+4+3)){
polAngle = atof(args[7+4+1].c_str() ); // azimuthal angle of the photon polarization vector in the lab measured in degrees.
polFraction = AmpParameter( args[7+4+2] ); // polarization fraction
if(args.size() == (6+4+3)){
polAngle = atof(args[6+4+1].c_str() ); // azimuthal angle of the photon polarization vector in the lab measured in degrees.
polFraction = AmpParameter( args[6+4+2] ); // polarization fraction
std::cout << "Fixed polarization fraction =" << polFraction << " and pol.angle= " << polAngle << " degrees." << std::endl;
}
else if (args.size() == (7+4+2)){//beam properties requires halld_sim
else if (args.size() == (6+4+2)){//beam properties requires halld_sim
// BeamProperties configuration file
TString beamConfigFile = args[7+4+1].c_str();
TString beamConfigFile = args[6+4+1].c_str();
cout<<beamConfigFile.Data()<<endl;
BeamProperties beamProp(beamConfigFile);
polFrac_vs_E = (TH1D*)beamProp.GetPolFrac();
polAngle = beamProp.GetPolAngle();
Expand All @@ -50,29 +51,21 @@ omegapi_amplitude::omegapi_amplitude( const vector< string >& args ):
//cout << polFrac_vs_E->GetBinContent(i) << endl;
}
}
else
assert(0);
else assert(0);

sign = atoi(args[0].c_str() );
lambda_gamma = atoi(args[1].c_str() );
spin = atoi(args[2].c_str() );
parity = atoi(args[3].c_str() );
spin_proj = atoi(args[4].c_str() );

c_0 = AmpParameter(args[5]);
registerParameter(c_0);

c_1 = AmpParameter(args[6]);
registerParameter(c_1);

c_2 = AmpParameter(args[7]);
registerParameter(c_2);

l = atoi(args[5].c_str() ); // partial wave l
nat_sign = atoi(args[6].c_str() ); // sign for amplitude given by naturality of exchange

//Dalitz Parameters
dalitz_alpha = AmpParameter(args[7+1]);
dalitz_beta = AmpParameter(args[7+2]);
dalitz_gamma = AmpParameter(args[7+3]);
dalitz_delta = AmpParameter(args[7+4]);
dalitz_alpha = AmpParameter(args[6+1]);
dalitz_beta = AmpParameter(args[6+2]);
dalitz_gamma = AmpParameter(args[6+3]);
dalitz_delta = AmpParameter(args[6+4]);

registerParameter(dalitz_alpha);
registerParameter(dalitz_beta);
Expand Down Expand Up @@ -109,7 +102,7 @@ omegapi_amplitude::calcUserVars( GDouble** pKin, GDouble* userVars ) const
GDouble Pgamma=polFraction;//fixed beam polarization fraction
if(polAngle == -1)
Pgamma = 0.;//if beam is amorphous set polarization fraction to 0
else if(polFrac_vs_E!=NULL){
else if(0) { //polFrac_vs_E!=NULL){
//This part causes seg fault with 34 amplitudes or more with gen_amp and gen_omegapi.
//Not needed for fixed beam pol angle and frac.
int bin = polFrac_vs_E->GetXaxis()->FindBin(beam.E());
Expand Down Expand Up @@ -140,9 +133,9 @@ omegapi_amplitude::calcUserVars( GDouble** pKin, GDouble* userVars ) const
double dalitz_s = rho.M2();//s=M2(pip pim)
double dalitz_t = (rhos_pip+omegas_pi).M2();//t=M2(pip pi0)
double dalitz_u = (rhos_pim+omegas_pi).M2();//u=M2(pim pi0)
double m3pi = (2*139.57018)+134.9766;
double m3pi = (2*0.13957018)+0.1349766;
double dalitz_d = 2*omega.M()*( omega.M() - m3pi);
double dalitz_sc = (1/3)*( omega.M2() - rhos_pip.M2() - rhos_pim.M2() - omegas_pi.M2());
double dalitz_sc = (1/3.)*( omega.M2() - rhos_pip.M2() - rhos_pim.M2() - omegas_pi.M2());
double dalitzx = sqrt(3)*(dalitz_t - dalitz_u)/dalitz_d;
double dalitzy = 3*(dalitz_sc - dalitz_s)/dalitz_d;
double dalitz_z = dalitzx*dalitzx + dalitzy*dalitzy;
Expand Down Expand Up @@ -171,36 +164,26 @@ omegapi_amplitude::calcAmplitude( GDouble** pKin, GDouble* userVars ) const
GDouble dalitz_z = userVars[uv_dalitz_z];
GDouble dalitz_sin3theta = userVars[uv_dalitz_sin3theta];

GDouble G = sqrt(1 + 2 * dalitz_alpha * dalitz_z + 2 * dalitz_beta * pow(dalitz_z,3/2) * dalitz_sin3theta
+ 2 * dalitz_gamma * pow(dalitz_z,2) + 2 * dalitz_delta * pow(dalitz_z,5/2) * dalitz_sin3theta );
GDouble G = sqrt(1 + 2 * dalitz_alpha * dalitz_z + 2 * dalitz_beta * pow(dalitz_z,3/2.) * dalitz_sin3theta
+ 2 * dalitz_gamma * pow(dalitz_z,2) + 2 * dalitz_delta * pow(dalitz_z,5/2.) * dalitz_sin3theta );

GDouble hel_c[3] = { c_0, c_1, c_2};

complex <GDouble> amplitude = CZero;
complex <GDouble> amplitude = CZero;

for (int lambda = -1; lambda <= 1; lambda++)//omega helicity
{
GDouble hel_amp = 0.0;

for(int l = 0; l <= 2; l++)//partial waves (l).
{//if ( (parity == -1 && l% 2 == 0) || (parity == 1 && l%2 != 0) ) continue;

hel_amp += hel_c[l] * clebschGordan(l, 1, 0, lambda, spin, lambda);
}//loop over l

amplitude += wignerD( spin, spin_proj, lambda, cosTheta, Phi ) * hel_amp * wignerD( 1, lambda, 0, cosThetaH, PhiH ) * G;
GDouble hel_amp = clebschGordan(l, 1, 0, lambda, spin, lambda);
amplitude += conj(wignerD( spin, spin_proj, lambda, cosTheta, Phi )) * hel_amp * conj(wignerD( 1, lambda, 0, cosThetaH, PhiH )) * G;
}//loop over lambda

// multiply by square root of photon spin density matrix in helicity basis
complex <GDouble> prefactor ( cos( lambda_gamma * prod_angle ), sin( lambda_gamma * prod_angle ));
if (sign == -1 && lambda_gamma == -1){amplitude *= -1*sqrt( ( 1 - (sign * polfrac) )/2 ) * prefactor;}
else{amplitude *= sqrt( ( 1 - (sign * polfrac) )/2 ) * prefactor;}

complex <GDouble> prefactor ( cos( lambda_gamma * prod_angle ), sin( lambda_gamma * prod_angle ));

if (sign == -1 && lambda_gamma == -1){amplitude *= -1*sqrt( ( 1 + (sign * polfrac) )/2 ) * prefactor;}
else{amplitude *= sqrt( ( 1 + (sign * polfrac) )/2 ) * prefactor;}

//Instead of the vertices "scale" in configuration file.
// if( (parity == +1 && ((spin+spin_proj)%2 == 0) ) || (parity == -1 && ((spin+spin_proj)%2 != 0) ) ){amplitude *= -1;}
// multiply amplitude by a sign given by the naturality of the exchange (set to +1 if unknown naturality)
amplitude *= nat_sign;

return amplitude;
return amplitude;
}

void omegapi_amplitude::updatePar( const AmpParameter& par ){
Expand All @@ -213,7 +196,7 @@ void
omegapi_amplitude::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const {

GPUomegapi_amplitude_exec( dimGrid, dimBlock, GPU_AMP_ARGS,
sign, lambda_gamma, spin, parity, spin_proj, c_0, c_1, c_2, dalitz_alpha, dalitz_beta, dalitz_gamma, dalitz_delta, polAngle, polFraction);
sign, lambda_gamma, spin, parity, spin_proj, l, dalitz_alpha, dalitz_beta, dalitz_gamma, dalitz_delta, polAngle, polFraction);
}
#endif //GPU_ACCELERATION

9 changes: 3 additions & 6 deletions src/libraries/AMPTOOLS_AMPS/omegapi_amplitude.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@

#ifdef GPU_ACCELERATION
void GPUomegapi_amplitude_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO,
int sign, int lambda_gamma, int spin, int parity, int spin_proj, GDouble c_0, GDouble c_1, GDouble c_2,
int sign, int lambda_gamma, int spin, int parity, int spin_proj, int l,
GDouble dalitz_alpha, GDouble dalitz_beta, GDouble dalitz_gamma, GDouble dalitz_delta, GDouble polAngle, GDouble polFraction);

#endif // GPU_ACCELERATION
Expand Down Expand Up @@ -79,11 +79,8 @@ complex< GDouble > calcAmplitude( GDouble** pKin, GDouble* userVars ) const;
int spin;
int parity;
int spin_proj;

AmpParameter c_0;
AmpParameter c_1;
AmpParameter c_2;

int l;
int nat_sign;
AmpParameter dalitz_alpha;
AmpParameter dalitz_beta;
AmpParameter dalitz_gamma;
Expand Down
111 changes: 111 additions & 0 deletions src/libraries/AMPTOOLS_DATAIO/OmegaPiPlotGenerator.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
//#include <ctime>
//#include <stdlib.h>
//#include <stdio.h>

//#include <cassert>
//#include <iostream>
//#include <string>
//#include <sstream>

#include "TLorentzVector.h"
#include "TLorentzRotation.h"

#include "AMPTOOLS_AMPS/omegapiAngles.h"

//#include <cmath>
//#include <complex>
//#include <vector>
//#include "TMath.h"

#include "AMPTOOLS_DATAIO/OmegaPiPlotGenerator.h"
#include "IUAmpTools/Histogram1D.h"
#include "IUAmpTools/Kinematics.h"

/* Constructor to display FitResults */
OmegaPiPlotGenerator::OmegaPiPlotGenerator( const FitResults& results ) :
PlotGenerator( results )
{
createHistograms();
}

/* Constructor for event generator (no FitResult) */
OmegaPiPlotGenerator::OmegaPiPlotGenerator( ) :
PlotGenerator( )
{
createHistograms();
}

void OmegaPiPlotGenerator::createHistograms( ) {
cout << " calls to bookHistogram go here" << endl;

bookHistogram( kOmegaPiMass, new Histogram1D( 200, 0.6, 2., "MOmegaPi", "Invariant Mass of #omega #pi") );
bookHistogram( kCosTheta, new Histogram1D( 50, -1., 1., "CosTheta", "cos#theta;cos#theta" ) );
bookHistogram( kPhi, new Histogram1D( 50, -1*PI, PI, "Phi", "#Phi; #Phi[rad.]" ) );
bookHistogram( kCosThetaH, new Histogram1D( 50, -1., 1., "CosTheta_H", "cos#theta_H;cos#theta_H" ) );
bookHistogram( kPhiH, new Histogram1D( 50, -1*PI, PI, "Phi_H", "#Phi_H; #Phi_H[rad.]" ) );
bookHistogram( kProd_Ang, new Histogram1D( 50, -1*PI, PI, "Prod_Ang", "Prod_Ang; Prod_Ang[rad.]" ) );
bookHistogram( kt, new Histogram1D( 100, 0, 2.0 , "t", "-t" ) );
bookHistogram( kRecoilMass, new Histogram1D( 100, 0.9, 1.9 , "MRecoil", "Invariant Mass of Recoil" ) );
bookHistogram( kTwoPiMass, new Histogram1D( 100, 0.25, 1.75, "MTwoPi", "Invariant Mass of 2 #pi" ) );

}

void
OmegaPiPlotGenerator::projectEvent( Kinematics* kin ){

//cout << "project event" << endl;
TLorentzVector beam = kin->particle( 0 );
TLorentzVector recoil = kin->particle( 1 );
TLorentzVector Xs_pi = kin->particle( 2 );//bachelor pi0
TLorentzVector omegas_pi = kin->particle( 3 );//omega's pi0
TLorentzVector rhos_pip = kin->particle( 4 );//pi-
TLorentzVector rhos_pim = kin->particle( 5 );//pi+

TLorentzVector two_pi = kin->particle( 2 );
for(uint i=6; i<kin->particleList().size(); i++) {
recoil += kin->particle(i);
two_pi += kin->particle(i);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
TLorentzVector rho = rhos_pip + rhos_pim;
TLorentzVector omega = rho + omegas_pi;
TLorentzVector X = omega + Xs_pi;

TLorentzVector target(0,0,0,0.938);
double polAngle = 0.0;
//cout << "Beam Polarization Angle Set to " << polAngle << endl;

double b1_mass = X.M();
double Mandt = fabs((target-recoil).M2());
double recoil_mass = recoil.M();

//////////////////////// Boost Particles and Get Angles//////////////////////////////////

//Helicity coordinate system
TLorentzVector Gammap = beam + target;

//Calculate decay angles in helicity frame
vector <double> locthetaphi = getomegapiAngles(polAngle, omega, X, beam, Gammap);

vector <double> locthetaphih = getomegapiAngles(rhos_pip, omega, X, Gammap, rhos_pim);

GDouble cosTheta = TMath::Cos(locthetaphi[0]);
GDouble Phi = locthetaphi[1];
GDouble cosThetaH = TMath::Cos(locthetaphih[0]);
GDouble PhiH = locthetaphih[1];
GDouble prod_angle = locthetaphi[2];

//cout << "calls to fillHistogram go here" << endl;
fillHistogram( kOmegaPiMass, b1_mass );
fillHistogram( kCosTheta, cosTheta );
fillHistogram( kPhi, Phi );
fillHistogram( kCosThetaH, cosThetaH );
fillHistogram( kPhiH, PhiH );
fillHistogram( kProd_Ang, prod_angle );
fillHistogram( kt, Mandt );
fillHistogram( kRecoilMass, recoil_mass );
fillHistogram( kTwoPiMass, two_pi.M() );

}

34 changes: 34 additions & 0 deletions src/libraries/AMPTOOLS_DATAIO/OmegaPiPlotGenerator.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
#if !(defined OMEGAPIPLOTGENERATOR)
#define OMEGAPIPLOTGENERATOR

#include <vector>
#include <string>

#include "IUAmpTools/PlotGenerator.h"

using namespace std;

class FitResults;
class Kinematics;

class OmegaPiPlotGenerator : public PlotGenerator
{

public:

// create an index for different histograms
enum { kOmegaPiMass = 0, kCosTheta = 1, kPhi = 2, kCosThetaH = 3, kPhiH = 4, kProd_Ang = 5, kt = 6, kRecoilMass = 7, kTwoPiMass = 8, kNumHists};

OmegaPiPlotGenerator( const FitResults& results );
OmegaPiPlotGenerator( );

void projectEvent( Kinematics* kin );

private:

void createHistograms( );

};

#endif

5 changes: 5 additions & 0 deletions src/libraries/AMPTOOLS_MCGEN/ProductionMechanism.cc
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,11 @@ ProductionMechanism::setGeneratorType( Type type ){
m_type = type;
}

void
ProductionMechanism::setRecoilMass( double recMass ){

m_recMass = recMass;
}

TLorentzVector
ProductionMechanism::produceResonance( const TLorentzVector& beam ){
Expand Down
1 change: 1 addition & 0 deletions src/libraries/AMPTOOLS_MCGEN/ProductionMechanism.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ class ProductionMechanism
void setMassRange( double low, double high );
void setTRange( double low, double high );
void setGeneratorType( Type type );
void setRecoilMass( double recMass );

TLorentzVector produceResonance( const TLorentzVector& beam );
TLorentzVector produceResonanceZ( const TLorentzVector& beam);
Expand Down
2 changes: 1 addition & 1 deletion src/programs/AmplitudeAnalysis/SConscript
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ import sbms

Import('*')

subdirs = ['fit', 'twopi_plotter', 'twopi_plotter_amp', 'twopi_plotter_mom', 'twopi_plotter_primakoff', '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']
subdirs = ['fit', 'twopi_plotter', 'twopi_plotter_amp', 'twopi_plotter_mom', 'twopi_plotter_primakoff', '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']

SConscript(dirs=subdirs, exports='env osname', duplicate=0)

2 changes: 2 additions & 0 deletions src/programs/AmplitudeAnalysis/fit/fit.cc
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
#include "AMPTOOLS_AMPS/Uniform.h"
#include "AMPTOOLS_AMPS/polCoef.h"
#include "AMPTOOLS_AMPS/dblRegge.h"
#include "AMPTOOLS_AMPS/omegapi_amplitude.h"

#include "MinuitInterface/MinuitMinimizationManager.h"
#include "IUAmpTools/AmpToolsInterface.h"
Expand Down Expand Up @@ -98,6 +99,7 @@ int main( int argc, char* argv[] ){
AmpToolsInterface::registerAmplitude( polCoef() );
AmpToolsInterface::registerAmplitude( Uniform() );
AmpToolsInterface::registerAmplitude( dblRegge() );
AmpToolsInterface::registerAmplitude( omegapi_amplitude() );

AmpToolsInterface::registerDataReader( ROOTDataReader() );
AmpToolsInterface::registerDataReader( ROOTDataReaderBootstrap() );
Expand Down
22 changes: 22 additions & 0 deletions src/programs/AmplitudeAnalysis/omegapi_plotter/SConscript
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@

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)
Loading

0 comments on commit 1f94234

Please sign in to comment.