IsoPsRefl amplitude optionally reads the beam polarization (eps-vector) direclty from an input tree

src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc

@@ -1,4 +1,3 @@
-
 #include <cassert>
 #include <iostream>
 #include <string>
@@ -44,9 +43,8 @@ static double parseValidatedNumber(const string& label, const string& argInput){
 }
 
 
-Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) :
-UserAmplitude< Iso_ps_refl >( args ){
-
+Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) : UserAmplitude< Iso_ps_refl >( args )
+{
   // This function is only for the two-body Isobar decay: xi -> pipi
 
   // 6 possibilities to initialize this amplitude:
@@ -84,72 +82,41 @@ UserAmplitude< Iso_ps_refl >( args ){
   // Default polarization information stored in tree
   m_polInTree = true;
 
-  // Loop over any additional amplitude arguments to change defaults
-  for(uint ioption=6; ioption<args.size(); ioption++) {
-	  TString option = args[ioption].c_str();
-    // Polarization provided in configuration file
-    if(ioption==6){
-      m_polInTree = false;
-
-      polAngle = parseValidatedNumber("polarization angle", args[6]);    
-
-      TString polOption = args[7].c_str();
-      if(polOption.Contains(".root")){
-        polFraction = 0.0;
-        TFile* f = new TFile(polOption);
-        polFrac_vs_E = (TH1D*)f->Get(args[8].c_str());
-        if(polFrac_vs_E  != nullptr ){
-          throw std::runtime_error(
-            "Iso_ps_refl ERROR: Could not find histogram '" + args[8] +
-            "' in file " + std::string(polOption.Data()));
-        }
-      }
+  // Polarization information passed with additional arguments
+  if (args.size() > 6){
+
+    m_polInTree = false;
+    polAngle = parseValidatedNumber("polarization angle", args[6]);    
+      std::string pol_option = args[7];      
+
+      if (pol_option.find(".root") == std::string::npos)
+	polFraction = parseValidatedNumber("polarization fraction", args[7]);
       else{
-        polFraction = parseValidatedNumber("polarization fraction", args[7]);
-      }
-    }
+	polFraction = -1;
+        TFile* pol_file = new TFile(pol_option.c_str());
+        polFrac_vs_E = (TH1D*)pol_file->Get(args[8].c_str());
+        if(polFrac_vs_E  == nullptr )
+          throw std::runtime_error("Iso_ps_refl ERROR: Could not find histogram '" + args[8] + "' in file " + pol_option);
+      }	
   }  
-}
 
-void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{
+}
 
-  TLorentzVector beam;
-  TVector3 eps;
-  GDouble beam_polFraction;
-  GDouble beam_polAngle;
 
-  if(m_polInTree){
-    beam.SetPxPyPzE( 0.0, 0.0, pKin[0][0], pKin[0][0]);
-    eps.SetXYZ(pKin[0][1], pKin[0][2], 0.0); // beam polarization vector;
-
-    beam_polFraction = eps.Mag();
-    beam_polAngle = eps.Phi();
-  }
-  else{
-    beam.SetPxPyPzE( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] );
-    beam_polAngle = polAngle;
-
-    if(polFraction > 0.0){ // for fitting with fixed polarization
-	    beam_polFraction = polFraction;
-    }
-    else{ // for fitting with polarization vs E_gamma from input histogram 
-	    int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]);
-	    if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){
-		    beam_polFraction = 0.0;
-	    } else
-		    beam_polFraction = polFrac_vs_E->GetBinContent(bin);
-    }
-  }
+void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const
+{  
 
+  TLorentzVector beam( 0.0, 0.0, pKin[0][3], pKin[0][0]);
+  TVector3 eps(pKin[0][1], pKin[0][2], 0.);   
 
+
   // Fill in four-vectors for final state particles
   TLorentzVector target(0,0,0,0.938272);   // Fixed target
   TLorentzVector recoil( pKin[1][1]+pKin[5][1], pKin[1][2]+pKin[5][2], pKin[1][3]+pKin[5][3], pKin[1][0]+pKin[5][0] ); // Recoiling baryon
   TLorentzVector ps(pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0]);     // 1st after proton is always the pseudoscalar meson
 
 
-  // Compute isobar meson from its decay products (xi -> pipi)
-  // Make sure the order of daughters is correct in the config file!
+  // Fill in four-vectors for the isobar meson (xi -> pipi, make sure the order of daughters is correct in the config file!)
   TLorentzVector iso_daught1(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]);
   TLorentzVector iso_daught2(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]);
   TLorentzVector iso = iso_daught1 + iso_daught2;
@@ -158,13 +125,29 @@ void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{
   TLorentzVector X = iso + ps;
 
 
-  //Calculate production angle in the Gottfried-Jackson frame
-  GDouble prod_angle = getPhiProd(beam_polAngle, X, beam, target, 2, true);
+  //Define Pgamma and calculate BigPhi angle in the Gottfried-Jackson frame
+  GDouble beam_polFraction, beam_polAngle, prod_angle;
+
+  if(m_polInTree){    
+    beam_polFraction = eps.Mag();
+    prod_angle = getPhiProd(eps, X, beam, target, 2, true);
+  }
+  else{
+    beam_polAngle = TMath::DegToRad()*polAngle; //conversion factor, since angles in the config files are by default in degrees
+    prod_angle = getPhiProd(beam_polAngle, X, beam, target, 2, true);
+    beam_polFraction = polFraction;     
+
+    if(beam_polFraction == -1){ // extract the polarization fraction from its dependence on E_gamma 
+      int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]);
+      if (bin > 0 && bin <= polFrac_vs_E->GetXaxis()->GetNbins())
+	beam_polFraction = polFrac_vs_E->GetBinContent(bin);
+    }
+  }
+
 
-  // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame  
+  // Calculate two angles of X decay (we use the Gottfried-Jackson frame) and two angles of Isobar decay (we use the Helicity frame)
   vector <double> thetaPhiAnglesTwoStep;
   thetaPhiAnglesTwoStep = getTwoStepAngles(X, iso, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true);
-
 
 
   GDouble cosTheta = TMath::Cos(thetaPhiAnglesTwoStep[0]);
@@ -175,6 +158,7 @@ void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{
   GDouble m_iso = iso.M();
   GDouble m_ps = ps.M();
 
+  // Calculate amplitude 
   complex <GDouble> amplitude(0,0);
   complex <GDouble> i(0,1);
 

src/libraries/AMPTOOLS_AMPS/decayAngles.cc

@@ -44,6 +44,7 @@ vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inv
 	return xyz;
 }
 
+
 // Calculate production plane angle
 double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex)
 {
@@ -89,6 +90,50 @@ double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beam
 } 
 
 
+// Calculate production plane angle using eps-vector already defined in user's amplitude
+double getPhiProd(TVector3 eps, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex)
+{
+	// whichFrame = 1 for helicity, 2 for GJ
+	assert( whichFrame == 1 || whichFrame == 2 );
+
+	// Boost all P4 from lab to CM rest frame
+	TLorentzVector cmLab = beamLab + targetLab;
+	TVector3 boostCMToLab = cmLab.BoostVector();
+
+	TLorentzVector parentCM = parentLab;
+	TLorentzVector inverseLab;
+	if( upperVertex == false )
+		inverseLab = targetLab;
+	else
+		inverseLab = beamLab;
+	TLorentzVector inverseCM = inverseLab;			
+	parentCM.Boost( -1.0*boostCMToLab );
+	inverseCM.Boost( -1.0*boostCMToLab );
+
+	TVector3 boostParentToCM = parentCM.BoostVector();
+	TLorentzVector inverseParent = inverseCM;
+	inverseParent.Boost( -1.0*boostParentToCM );
+
+	// convert P4s to P3
+	TVector3 parentCM3 = parentCM.Vect();
+	TVector3 inverseCM3 = inverseCM.Vect();
+	TVector3 inverseParent3 = inverseParent.Vect();
+
+	vector< TVector3 > locxyz;
+	if( whichFrame == 1 )
+		locxyz = getHelicityAxes(parentCM3, inverseCM3);
+	else
+		locxyz = getGJAxes(parentCM3, inverseCM3, inverseParent3);
+
+	TVector3 y = locxyz[1];
+
+	double phiProd = atan2( y.Dot( eps ), inverseCM.Vect().Unit().Dot( eps.Cross( y ) ) );
+
+	return phiProd;	
+} 
+
+
+
 // Calculate angles for a one-step decay. If only one of the daughter particles has spin, that one should be used in this calculation 
 vector< double > getOneStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex)
 {

src/libraries/AMPTOOLS_AMPS/decayAngles.h

@@ -16,6 +16,7 @@ vector< TVector3 > getHelicityAxes(TVector3 parentCM, TVector3 inverseCM);
 vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inverseParent);
 
 double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex);
+double getPhiProd(TVector3 eps, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex);
 
 vector< double > getOneStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex = true);
 

src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc

@@ -62,11 +62,11 @@ PlotGenerator( )
 void IsoPsPlotGenerator::createHistograms( ) {
   cout << " calls to bookHistogram go here" << endl;
 
-   bookHistogram( kProd_Ang, new Histogram1D( 50, -180., 180., "ProdAng", "Production Angle [deg.]" ) );
+   bookHistogram( kProd_Ang, new Histogram1D( 50, -180., 180., "ProdAng", "#Phi [deg]" ) );
    bookHistogram( kCosTheta, new Histogram1D( 50, -1., 1., "CosTheta_GJ", "cos#theta^{[GJ]}" ) );
-   bookHistogram( kPhi, new Histogram1D( 50, -180., 180., "Phi_GJ", "#phi^{[GJ]} [deg.]" ) );
+   bookHistogram( kPhi, new Histogram1D( 50, -180., 180., "Phi_GJ", "#phi^{[GJ]} [deg]" ) );
    bookHistogram( kCosThetaH, new Histogram1D( 50, -1., 1., "CosTheta_HF", "cos#theta^{[HF]}" ) );
-   bookHistogram( kPhiH, new Histogram1D( 50, -180., 180., "Phi_HF", "#phi^{[HF]} [deg.]" ) );
+   bookHistogram( kPhiH, new Histogram1D( 50, -180., 180., "Phi_HF", "#phi^{[HF]} [deg]" ) );
 
    bookHistogram( kIsoMass, new Histogram1D( 200, 0., 3., "MIso", "m(2#pi)  [GeV]") );
    bookHistogram( kIsoPsMass, new Histogram1D( 200, 0.2, 3.2, "MIsoPs", "m(3#pi)  [GeV]") );
@@ -78,24 +78,22 @@ void IsoPsPlotGenerator::createHistograms( ) {
 
 }
 
-void
-IsoPsPlotGenerator::projectEvent( Kinematics* kin ){
+void IsoPsPlotGenerator::projectEvent( Kinematics* kin ){
 
   // this function will make this class backwards-compatible with older versions
   // (v0.10.x and prior) of AmpTools, but will not be able to properly obtain
   // the polariation plane in the lab when multiple orientations are used
   projectEvent( kin, "" );
 }
 
-void
-IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
+void IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
 
-   // We work only with a 2-body vector decay 
+   // PlotGenerator works only with a 2-body vector decay   
 
-
-  // Fixed target
+   // Fixed target
    TLorentzVector target(0,0,0,0.938272);
-
+   TVector3 eps;
+
 
    //cout << "project event" << endl;
    TLorentzVector beam   = kin->particle( 0 );
@@ -104,8 +102,8 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
    TLorentzVector iso_daught1 = kin->particle( 3 );
    TLorentzVector iso_daught2 = kin->particle( 4 );
    TLorentzVector piplusL = kin->particle( 5 );
+
 
-
    // Final state P4 momenta
    TLorentzVector X = iso_daught1 + iso_daught2 + bach;
    TLorentzVector recoil = proton + piplusL;
@@ -120,25 +118,34 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
    TLorentzVector proton_ps_b = proton + iso_daught1;
    TLorentzVector recoil_ps_b = recoil + iso_daught1;
 
-
+   //Momentum transfer
+   double Mandt = fabs((target-recoil).M2());
 
-   // Properly read polarization angle from config file if provided
-   double beam_polAngle=0;
-   // Check config file for optional parameters -- we assume here that the first amplitude in the list is a Iso_ps_refl amplitude
+
+
+   // Either read the diamond orientation angle directly from the config file (if provided), or read eps-vector from the (Beam_PX,Beam_PY)-components in the input tree   
+   // Then calculate polarization angle in the Gottfried-Jackson frame
+   // We assume here that the first amplitude in the list is always the Iso_ps_refl amplitude
    const vector< string > args = cfgInfo()->amplitudeList( reactionName, "", "" ).at(0)->factors().at(0);
-   for(uint ioption=5; ioption<args.size(); ioption++) {
-          TString option = args[ioption].c_str();
-	  if(ioption == 6) beam_polAngle = parseValidatedNumber("polarization angle", args[6]);
+   double beam_polAngle, prod_angle;
+
+   if (args.size() > 6) {    
+     beam_polAngle = TMath::DegToRad()*parseValidatedNumber("polarization angle", args[6]);    
+     prod_angle = TMath::RadToDeg()*getPhiProd(beam_polAngle, X, beam, target, 2, true);
+   }
+   else {
+     eps.SetXYZ(beam.X(),beam.Y(),0.);
+     prod_angle = TMath::RadToDeg()*getPhiProd(eps, X, beam, target, 2, true);
    }
-
-
-
-   //Momentum transfer
-   double Mandt = fabs((target-recoil).M2());
 
-   //Calculate production angle in the Gottfried-Jackson frame
-   double prod_angle = TMath::RadToDeg()*getPhiProd(beam_polAngle, X, beam, target, 2, true);
-
+   //forcibly put xy-components of the beam vector to zero     
+   beam.SetX(0.);
+   beam.SetY(0.);
+
+
+
+
+
    // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame  
    // Angles for the 1st permutation
    vector <double> thetaPhiAnglesTwoStep_a = getTwoStepAngles(X, iso_a, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true);

src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc

@@ -128,8 +128,8 @@ int main( int argc, char* argv[] ){
 	atiSetup();
         cout << "Plotgen results"<< endl;
 
-	//	IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // optional can be omitted
-	IsoPsPlotGenerator plotGen( results ); 
+	//	IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // 2nd argument is optional
+	IsoPsPlotGenerator plotGen( results );
 	cout << " Initialized ati and PlotGen" << endl;
 
 
@@ -409,7 +409,7 @@ int main( int argc, char* argv[] ){
     for(unsigned int j = i+1; j < fullamps.size(); j++){
 
       // leave only the Spring2017_PARA_0::ImagPosSign and Spring2017_PARA_0::ImagNegSign coherent sums
-      if (fullamps[i].find("Spring2017_PARA_0") == std::string::npos || fullamps[j].find("Spring2017_PARA_0") == std::string::npos) continue;      
+      if (fullamps[i].find("Spring2017") == std::string::npos || fullamps[j].find("Spring2017") == std::string::npos) continue;      
       if (fullamps[i].find("UniBG") != std::string::npos || fullamps[j].find("UniBG") != std::string::npos) continue;      
       if (fullamps[i].find("Real") != std::string::npos || fullamps[j].find("Real") != std::string::npos) continue;      
 

src/programs/Simulation/gen_amp_V2/gen_amp_v2.cc

@@ -28,6 +28,8 @@
 #include "AMPTOOLS_AMPS/Lambda1520Angles.h"
 #include "AMPTOOLS_AMPS/Lambda1520tdist.h"
 #include "AMPTOOLS_AMPS/Vec_ps_refl.h"
+#include "AMPTOOLS_AMPS/Iso_ps_refl.h"
+#include "AMPTOOLS_AMPS/PiPiSWaveAMPK.h"
 #include "AMPTOOLS_AMPS/Ylm.h"
 #include "AMPTOOLS_AMPS/Zlm.h"
 #include "AMPTOOLS_AMPS/DblRegge_FastEta.h"
@@ -310,6 +312,8 @@ int main( int argc, char* argv[] ){
   AmpToolsInterface::registerAmplitude( Lambda1520Angles() );
   AmpToolsInterface::registerAmplitude( Lambda1520tdist() );
   AmpToolsInterface::registerAmplitude( Vec_ps_refl() );
+  AmpToolsInterface::registerAmplitude( Iso_ps_refl() );
+  AmpToolsInterface::registerAmplitude( PiPiSWaveAMPK() );
   AmpToolsInterface::registerAmplitude( Ylm() );
   AmpToolsInterface::registerAmplitude( Zlm() );
   AmpToolsInterface::registerAmplitude( Hist2D() );