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,10 +43,9 @@ static double parseValidatedNumber(const string& label, const string& argInput){
 }
 
 
-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
+Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) : UserAmplitude< Iso_ps_refl >( args )
+{
+  // This function applies only to the two-body Isobar decay, such as xi -> pipi
 
   // 6 possibilities to initialize this amplitude:
   // <S>: isobar spin
@@ -57,6 +55,14 @@ UserAmplitude< Iso_ps_refl >( args ){
   // <r>: +1/-1 for real/imaginary part; 
   // <s>: +1/-1 sign in P_gamma term
 
+  // By default, beam polarization info is read from the input tree  
+  // 3 additional arguments to pass the beam polarization pars in this amplitude
+  // <Phi0>: diamond orientation angle in degrees
+  // <Pgamma> or <*.root>: either Pgamma value, or root-file with a TH1D Pgamma_vs_E histogram 
+  // <histname>: name of the histogram 
+
+
+
   // Isobar spin S
   m_s = static_cast<int>(parseValidatedNumber("S",args[0]));    
   // Resonance spin J
@@ -84,72 +90,43 @@ 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()));
-        }
-      }
-      else{
-        polFraction = parseValidatedNumber("polarization fraction", args[7]);
-      }
-    }
-  }  
-}
+  // Polarization information passed with additional arguments
+  if (args.size() > 6){
 
-void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{
+    m_polInTree = false;
+    polAngle = parseValidatedNumber("polarization angle", args[6]);    
+    std::string pol_option = args[7];      
 
-  TLorentzVector beam;
-  TVector3 eps;
-  GDouble beam_polFraction;
-  GDouble beam_polAngle;
+      if (pol_option.find(".root") == std::string::npos)
+	polFraction = parseValidatedNumber("polarization fraction", args[7]);
+      else{
+	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);
+	polFrac_vs_E->SetDirectory(0); 	// Detach histogram from the file, so it persists after the file is closed
+	pol_file->Close();
+      }	
+  }  
 
-  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 +135,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 +168,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

@@ -23,7 +23,7 @@
 #include "TMath.h"
 
 // Calculate the x, y, and z axes in the helicity reference frame
-vector< TVector3 > getHelicityAxes(TVector3 parentCM, TVector3 inverseCM)
+vector< TVector3 > getHelicityAxes(const TVector3 &parentCM, const TVector3 &inverseCM)
 {
 	TVector3 z = parentCM.Unit();
 	TVector3 y = ( inverseCM.Unit() ).Cross( parentCM.Unit() ).Unit();
@@ -34,7 +34,7 @@ vector< TVector3 > getHelicityAxes(TVector3 parentCM, TVector3 inverseCM)
 } 
 
 // Calculate the x, y, and z axes in the Gottfried-Jackson reference frame
-vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inverseParent)
+vector< TVector3 > getGJAxes(const TVector3 &parentCM, const TVector3 &inverseCM, const TVector3 &inverseParent)
 {
 	TVector3 z = inverseParent.Unit();
 	TVector3 y = ( inverseCM.Unit() ).Cross( parentCM.Unit() ).Unit();
@@ -44,8 +44,9 @@ 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)
+double getPhiProd(double polAngle, const TLorentzVector &parentLab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex)
 {
 	// whichFrame = 1 for helicity, 2 for GJ
 	assert( whichFrame == 1 || whichFrame == 2 );
@@ -89,8 +90,52 @@ double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beam
 } 
 
 
+// Calculate production plane angle using already defined eps-vector
+double getPhiProd(const TVector3 &eps, const TLorentzVector &parentLab, const TLorentzVector &beamLab, const 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)
+vector< double > getOneStepAngles(const TLorentzVector &parentLab, const TLorentzVector &daughterLab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex)
 {
 	// whichFrame = 1 for helicity, 2 for GJ
 	assert( whichFrame == 1 || whichFrame == 2 );
@@ -148,7 +193,7 @@ vector< double > getOneStepAngles(TLorentzVector parentLab, TLorentzVector daugh
 
 
 // Calculate angles for a two-step decay. This has been tested for photoproduction of a resonance decaying to omega and pi, with the omega subequently decaying to three pions. If the second step of the decay is a two-body decay, the input for granddaughter2Lab should be TLorentzVector(0,0,0,0)  
-vector< double > getTwoStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector granddaughter1Lab, TLorentzVector granddaughter2Lab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex)
+vector< double > getTwoStepAngles(const TLorentzVector &parentLab, const TLorentzVector &daughterLab, const TLorentzVector &granddaughter1Lab, const TLorentzVector &granddaughter2Lab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex)
 {
 	// whichFrame = 1 for helicity, 2 for GJ
 	assert( whichFrame == 1 || whichFrame == 2 );

src/libraries/AMPTOOLS_AMPS/decayAngles.h

@@ -11,14 +11,15 @@
 using std::complex;
 using namespace std;
 
-vector< TVector3 > getHelicityAxes(TVector3 parentCM, TVector3 inverseCM);
+vector< TVector3 > getHelicityAxes(const TVector3 &parentCM, const TVector3 &inverseCM);
 
-vector< TVector3 > getGJAxes(TVector3 parentCM, TVector3 inverseCM, TVector3 inverseParent);
+vector< TVector3 > getGJAxes(const TVector3 &parentCM, const TVector3 &inverseCM, const TVector3 &inverseParent);
 
-double getPhiProd(double polAngle, TLorentzVector parentLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex);
+double getPhiProd(double polAngle, const TLorentzVector &parentLab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex);
+double getPhiProd(const TVector3 &eps, const TLorentzVector &parentLab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex);
 
-vector< double > getOneStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex = true);
+vector< double > getOneStepAngles(const TLorentzVector &parentLab, const TLorentzVector &daughterLab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex = true);
 
-vector< double > getTwoStepAngles(TLorentzVector parentLab, TLorentzVector daughterLab, TLorentzVector granddaughter1Lab, TLorentzVector granddaughter2Lab, TLorentzVector beamLab, TLorentzVector targetLab, int whichFrame, bool upperVertex);
+vector< double > getTwoStepAngles(const TLorentzVector &parentLab, const TLorentzVector &daughterLab, const TLorentzVector &granddaughter1Lab, const TLorentzVector &granddaughter2Lab, const TLorentzVector &beamLab, const TLorentzVector &targetLab, int whichFrame, bool upperVertex);
 
 #endif

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,36 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
    TLorentzVector proton_ps_b = proton + iso_daught1;
    TLorentzVector recoil_ps_b = recoil + iso_daught1;
 
-
-
-   // Properly read polarization angle from config file if provided
-   double beam_polAngle=0;
-   // Check config file for optional parameters -- we assume here that the first amplitude in the list is a Iso_ps_refl amplitude
-   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]);
-   }
-
-
-
    //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);
+
+   // 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);
+   double beam_polAngle, prod_angle;
 
+
+   //Initialize eps-vector and forcibly set xy-components of the beam vector to zero     
+   eps.SetXYZ(beam.X(),beam.Y(),0.);
+   beam.SetX(0.);
+   beam.SetY(0.);
+
+
+   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 {
+     prod_angle = TMath::RadToDeg()*getPhiProd(eps, X, beam, target, 2, true);
+   }
+
+
+
+
+
    // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame  
    // Angles for the 1st permutation
    vector <double> thetaPhiAnglesTwoStep_a = getTwoStepAngles(X, iso_a, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true);