Ilbelov 3pi fits

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,60 +82,46 @@ 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;
+void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const
+{  
 
-  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;
+  GDouble beam_polFraction, beam_polAngle;
+  TLorentzVector beam( 0.0, 0.0, pKin[0][3], pKin[0][0]);
 
-    beam_polFraction = eps.Mag();
-    beam_polAngle = eps.Phi();
+
+  if(m_polInTree){    
+    beam_polAngle = TMath::DegToRad()*pKin[0][1]; //Px
+    beam_polFraction = pKin[0][2];  //Py  
   }
   else{
-    beam.SetPxPyPzE( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] );
-    beam_polAngle = polAngle;
+    beam_polAngle = TMath::DegToRad()*polAngle;
+    beam_polFraction = polFraction;     
 
-    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);
+    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);
     }
   }
 

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,22 +78,20 @@ 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 
+   // We work only with a 2-body vector decay   
 
-  
-  // Fixed target
+
+   // Fixed target
    TLorentzVector target(0,0,0,0.938272);
 
 
@@ -104,8 +102,9 @@ 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;
@@ -121,24 +120,29 @@ IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
    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
+
+   // Read polarization angle either from config file (if provided), or from the X-component of p_beam
+   // 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;
+
+   if (args.size() > 6)     
+     beam_polAngle = TMath::DegToRad()*parseValidatedNumber("polarization angle", args[6]);    
+   else 
+     beam_polAngle = TMath::DegToRad()*beam.X(); //Px should be the 0th position in TLorentzVector
+
+   //forcibly put xy-components of the beam vector to zero     
+   beam.SetX(0.);
+   beam.SetY(0.);
+
 
    //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);
-
+
+
    // 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"
@@ -308,6 +310,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() );