new 4pi phase space factor; update of PiPiSWaveTMatrix.cc

Examples/Tutorial/LineShapes/PiPiSWaveTMatrix.cc

@@ -37,6 +37,7 @@
 //#include "PwaDynamics/KPoleNonRel.hh"
 #include "PwaDynamics/KMatrixSlowAdlerCorRel.hh"
 #include "PwaDynamics/KMatrixNonRel.hh"
+#include "Utils/PawianConstants.hh"
 
 #include "TFile.h"
 #include "TH1F.h"
@@ -76,12 +77,12 @@ PiPiSWaveTMatrix::PiPiSWaveTMatrix() :
   _argandH2=new TH2F("_argandH2","Argand plot K matrix",301, -1., 1., 301, 0., 1.3);
   _argandH2->SetYTitle("Im(T_{00})");
   _argandH2->SetXTitle("Re(T_{00})");
-  _argandRelH2=new TH2F("_argandRelH2","Argand plot K matrix (rel)",301, -1., 1., 301, -1.3, 1.3);
+  _argandRelH2=new TH2F("_argandRelH2","Argand plot K matrix (rel)",301, -1., 1., 301, 0., 1.3);
   _argandRelH2->SetYTitle("Im(#rho_{00}#hat{T}_{00})");
   _argandRelH2->SetXTitle("Re(#rho_{00}#hat{T}_{00})");
-  _phaseShiftH2=new TH2F("_phaseShiftH2", "phase shift",301, massMin, massMax, 301, 0., 3.1415);
-  _phaseShiftRelH2=new TH2F("_phaseShiftRelH2", "phase shift (rel)",301, massMin, massMax, 301, 0., 3.1415);
-
+  _phaseShiftDegH2=new TH2F("_phaseShiftDegH2", "phase shift deg",301, massMin, massMax, 301, 0., 360.);
+  _phaseShiftDegRelH2=new TH2F("_phaseShiftDegRelH2", "phase shift deg rel",301, massMin, massMax, 301, 0., 360.);
+  _elasticityH1=new TH1F("_elasticityH1", "elasticity", size, massMin, massMax);
 
 
   std::shared_ptr<KMatrixSlowAdlerCorRel> theKMatrix(new KMatrixPiPiS());
@@ -104,34 +105,67 @@ PiPiSWaveTMatrix::PiPiSWaveTMatrix() :
   std::shared_ptr<KMatrixSlowAdlerCorRel> theKMatrixSigmaPole(new KMatrixSlowAdlerCorRel(sigmaPoleVec, thePhpVecs, fScatProd, soScat));
   std::shared_ptr<TMatrixBase> theTMatrixSigmaPole(new TMatrixRel(theKMatrixSigmaPole));
 
+  double oldT00Real=1.;
+  int n180Shift(0);
+  double oldT00RelReal=1.;
+  int n180ShiftRel(0);
+
   for (double mass=massMin; mass<massMax; mass+=stepSize){
     Vector4<double> mass4Vec(mass, 0.,0.,0.);
     
     theTMatrixNonRel->evalMatrix(mass);
         
-    complex<double> currentValLow=(*theTMatrixNonRel)(0,0);
-    _invPiPiMassH1->Fill(mass4Vec.M(), norm(currentValLow) );    
-    _argandH2->Fill(currentValLow.real(),currentValLow.imag());
-    _phaseShiftH2->Fill(mass, atan2(currentValLow.imag(), currentValLow.real())); 
-
+    complex<double> currentT00=(*theTMatrixNonRel)(0,0);
+    _invPiPiMassH1->Fill(mass4Vec.M(), norm(currentT00) );    
+    _argandH2->Fill(currentT00.real(),currentT00.imag());
+
+    // T - E = 0.5*i, where E = inelasticity vector, pointing to T from (0,i/2)
+    double currentReE = currentT00.real();
+    double currentImE = currentT00.imag()- 0.5;;
+
+    // Find the phase shift angle, delta
+    double delta = 0.5*atan2(currentImE, fabs(currentReE))*PawianConstants::radToDeg + 45.0;
+    if (currentT00.real() < 0.0) {delta = 180.0 - delta;}
+    // Have we gone through 180 deg (or 2*delta through 360 deg)?
+    if (oldT00Real < 0.0 && currentT00.real() > 0.0) {n180Shift += 1;}
+    // Add 180 if required
+    delta += 180.0*n180Shift;
+    _phaseShiftDegH2->Fill(mass, delta);
 
     theTMatrix->evalMatrix(mass);    
-    complex<double> currentValLowRel=(*theTMatrix)(0,0);
+    complex<double> currentT00Rel=(*theTMatrix)(0,0);
     complex<double> S00_rel=complex<double>(1.,0.)+2.*complex<double>(0.,1.)*thePhpVecs[0]->factor(mass4Vec.M()).real()*(*theTMatrix)(0,0);
 
-    _invPiPiMassRelH1->Fill(mass4Vec.M(), norm(sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentValLowRel*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() ) ) );
-    _absT00RelH1->Fill(mass4Vec.M(), sqrt(norm(sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentValLowRel*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() ) )) );
+    _invPiPiMassRelH1->Fill(mass4Vec.M(), norm(sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentT00Rel*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() ) ) );
+    _absT00RelH1->Fill(mass4Vec.M(), sqrt(norm(sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentT00Rel*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() ) )) );
     _absS00RelH1->Fill(mass4Vec.M(), sqrt(norm(S00_rel))); 
        
-    _argandRelH2->Fill( sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentValLowRel.real()*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() ),sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentValLowRel.imag()*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )); 
-    _phaseShiftRelH2->Fill(mass, atan2(currentValLowRel.imag(), currentValLowRel.real()));
+    _argandRelH2->Fill( sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentT00Rel.real()*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() ),sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )*currentT00Rel.imag()*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real() )); 
 
     _pipiPhaseSpaceFactorH2->Fill(mass, thePhpVecs[0]->factor(mass4Vec.M()).real() );
     _pipipipiPhaseSpaceFactorH2->Fill(mass, thePhpVecs[2]->factor(mass4Vec.M()).real() );
 
     theTMatrixSigmaPole->evalMatrix(mass);
     complex<double> currentValRelSigmaPole=(*theTMatrixSigmaPole)(0,0);
-    _sqrT00RelSigmaPoleH1->Fill( mass4Vec.M(), norm( sqrt((thePhpVecs[0]->factor(mass4Vec.M())).real())*currentValRelSigmaPole*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real())) );     
+    _sqrT00RelSigmaPoleH1->Fill( mass4Vec.M(), norm( sqrt((thePhpVecs[0]->factor(mass4Vec.M())).real())*currentValRelSigmaPole*sqrt( (thePhpVecs[0]->factor(mass4Vec.M())).real())) );
+
+    _elasticityH1->Fill(mass4Vec.M(), sqrt(norm(S00_rel)));
+
+    // phase shift relat. matrix
+    complex<double> currentT00Rel_rho=(*theTMatrix)(0,0)*thePhpVecs[0]->factor(mass4Vec.M());
+    double currentReERel = currentT00Rel_rho.real();
+    double currentImERel = currentT00Rel_rho.imag() - 0.5;;
+
+    // Find the phase shift angle, delta
+    double deltaRel = 0.5*atan2(currentImERel, fabs(currentReERel))*PawianConstants::radToDeg + 45.0;
+    if (currentT00Rel_rho.real()  < 0.0) {deltaRel = 180.0 - deltaRel;}
+    // Have we gone through 180 deg (or 2*delta through 360 deg)?
+    if (oldT00RelReal < 0.0 && currentT00Rel_rho.real() > 0.0) {n180ShiftRel += 1;}
+    // Add 180 if required
+    deltaRel += 180.0*n180ShiftRel;
+    _phaseShiftDegRelH2->Fill(mass, deltaRel);
+    oldT00Real=currentT00.real();
+    oldT00RelReal=currentT00Rel_rho.real();     
   }
 }
 

Examples/Tutorial/LineShapes/PiPiSWaveTMatrix.hh

@@ -73,8 +73,9 @@ private:
   TH2F* _pipipipiPhaseSpaceFactorH2;
   TH2F*  _argandH2;
   TH2F*  _argandRelH2;
-  TH2F* _phaseShiftH2;
-  TH2F* _phaseShiftRelH2;
+  TH2F* _phaseShiftDegH2;
+  TH2F* _phaseShiftDegRelH2;
+  TH1F* _elasticityH1;
 };
 
 

PwaDynamics/PhaseSpace4Pi.cc

@@ -21,13 +21,45 @@
 //									  //
 //************************************************************************//
 
+// this paramererizaton of this 4pi phasespce factor has been taken over from 
+// the Laura++ software package:
+// 
+// Copyright University of Warwick 2008 - 2013.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+//
+// Calculate the 4pi phase space factor. This uses a 6th-order polynomial 
+// parameterisation that approximates the multi-body phase space double integral
+// defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the 
+// BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the
+// exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2.
+// Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2 
+// for s1 and s2, the invariant energy squared of each of the di-pion states,
+// with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and
+// s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be
+// 0.775 GeV and the energy-dependent width of the 4pi system 
+// Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is 
+// the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV 
+// (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state).
+// The normalisation term rho_0 is found by ensuring that the phase space integral
+// at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this 
+// factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5.
+// Plotting the value of this double integral as a function of s can then be fitted
+// to a 6th-order polynomial (for s < 1), which is the result used below
+
+
 #include "PwaDynamics/PhaseSpace4Pi.hh"
+#include "Utils/PawianConstants.hh"
 #include "qft++/relativistic-quantum-mechanics/Utils.hh"
 #include "ErrLogger/ErrLogger.hh"
 
 PhaseSpace4Pi::PhaseSpace4Pi():
   AbsPhaseSpace()
-  , _piMass(0.13957)
+  ,_fourPiFactor1(16.0*PawianConstants::mPiSq)
 {
 }
 
@@ -38,14 +70,20 @@ PhaseSpace4Pi::~PhaseSpace4Pi(){
 complex<double> PhaseSpace4Pi::factor(const double mass){
   double mass_sqr=mass*mass;
   complex<double> result(0.,0.);
-  if( mass_sqr <= 1 ){
-      double real = 1.2274 + .00370909 / ( mass_sqr * mass_sqr ) - .111203 / mass_sqr - 6.39017 * mass_sqr + 
-	 16.8358*mass_sqr*mass_sqr - 21.8845*mass_sqr*mass_sqr*mass_sqr + 11.3153*mass_sqr*mass_sqr*mass_sqr*mass_sqr;
 
-      double cont32 = sqrt(1.0-(16.0*_piMass*_piMass));
-      result = complex<double>( cont32 * real, 0 );
-    }
-  else result = complex<double>( sqrt( 1. - 16. * _piMass * _piMass / mass_sqr ), 0. );
+  if (fabs(mass_sqr) < 1e-10) {return result;}
+
+  if (mass_sqr <= 1.0) {
+    double rhoTerm = ((1.07885*mass_sqr + 0.13655)*mass_sqr - 0.29744)*mass_sqr - 0.20840;
+    rhoTerm = ((rhoTerm*mass_sqr + 0.13851)*mass_sqr - 0.01933)*mass_sqr + 0.00051;
+    // For some values of s (below 2*mpi), this term is a very small 
+    // negative number. Check for this and set the rho term to zero.
+    if (rhoTerm < 0.0) {rhoTerm = 0.0;}
+    result= complex<double>( rhoTerm, 0. );
+  } else {
+    result= complex<double>( sqrt(1.0 - (_fourPiFactor1/mass_sqr)), 0. );
+  }  
+
   return result;
 }
 
@@ -53,20 +91,13 @@ complex<double> PhaseSpace4Pi::factor(const complex<double> mass){
   complex<double> mass_sqr=mass*mass;
   complex<double> result(0.,0.);
   if( norm(mass_sqr) <= 1. ){
-    complex<double> real = 1.2274 + .00370909 / ( mass_sqr * mass_sqr ) - .111203 / mass_sqr - 6.39017 * mass_sqr + 
-	 16.8358*mass_sqr*mass_sqr - 21.8845*mass_sqr*mass_sqr*mass_sqr + 11.3153*mass_sqr*mass_sqr*mass_sqr*mass_sqr;
-
-      double cont32 = sqrt(1.0-(16.0*_piMass*_piMass));
-      result = cont32*real;
-      //      result = complex<double>( cont32 * real, 0 );
+    complex<double> rhoTerm = ((1.07885*mass_sqr + 0.13655)*mass_sqr - 0.29744)*mass_sqr - 0.20840;
+    rhoTerm = ((rhoTerm*mass_sqr + 0.13851)*mass_sqr - 0.01933)*mass_sqr + 0.00051;
+    result = rhoTerm;
     }
-  //  else result = complex<double>( sqrt( 1. - 16. * _piMass * _piMass / mass_sqr ), 0. );
- else result = sqrt( 1. - 16. * _piMass * _piMass / mass_sqr );
+  else result = sqrt(1.0 - (_fourPiFactor1/mass_sqr));
   CorrectForChosenSign(result);
   return result;
-   // Alert << "PhaseSpace4Pi does currently not support complex masses" << endmsg;
-   // exit(EXIT_FAILURE);
-   // return 0;
 }
 
 complex<double> PhaseSpace4Pi::breakUpMom(const double mass){

PwaDynamics/PhaseSpace4Pi.hh

@@ -61,7 +61,7 @@ public:
 protected:
 
 private:
-  const double _piMass;
+  const double _fourPiFactor1;
 };
 //_____________________________________________________________________________
 

PwaUtils/GlobalEnv.hh

@@ -41,7 +41,6 @@ class GlobalEnv
 
 public:
    static GlobalEnv* instance();
-   GlobalEnv();
    ~GlobalEnv();
 
    void setup(ParserBase* theParser);
@@ -66,6 +65,7 @@ public:
    std::vector<std::string> fixedParams();
 
 private:
+  GlobalEnv();
    static GlobalEnv* _instance;
    bool _alreadySetUp;
    bool _channelEnvsAlredySetup;

Scripts/pawianfixparametersEpEm.pl

+#!/usr/bin/perl
+
+#************************************************************************#
+#                                                                        #
+#  Copyright 2014 Bertram Kopf (bertram@ep1.rub.de)                      #
+#                 Julian Pychy (julian@ep1.rub.de)                       #
+#                 - Ruhr-Universität Bochum                              #
+#                                                                        #
+#  This file is part of Pawian.                                          #
+#                                                                        #
+#  Pawian is free software: you can redistribute it and/or modify        #
+#  it under the terms of the GNU General Public License as published by  #
+#  the Free Software Foundation, either version 3 of the License, or     #
+#  (at your option) any later version.                                   #
+#                                                                        #
+#  Pawian is distributed in the hope that it will be useful,             #
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of        #
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         #
+#  GNU General Public License for more details.                          #
+#                                                                        #
+#  You should have received a copy of the GNU General Public License     #
+#  along with Pawian.  If not, see <http:#www.gnu.org/licenses/>.        #
+#                                                                        #
+#************************************************************************#
+
+#######################################################
+##                                                   ##
+##  pawianfixparametersEpEm                          ##
+##                                                   ##
+##  This script generated the mnParFix list to be    ##
+##  used in the Pawian configuration file for        ##
+##  e+ e- reactions.                                 ##
+##                                                   ##
+##  - Bertram Kopf     bertram@ep1.rub.de            ##
+##                                                   ##
+#######################################################
+
+use strict;
+use warnings;
+use Getopt::Long;
+
+my $showHelp=0;
+my $filename="";
+my $prio="";
+
+GetOptions ("help|h" => \$showHelp,
+	    "file=s" => \$filename,
+	    "prio=s" => \$prio)  or die "Bad arguments!";
+
+if($showHelp){
+    Help();
+    exit;
+}
+
+if($filename eq "" || !(-e $filename)){
+    die "File not found.";
+}
+
+
+# Read the parameter file
+open FILE, "<", $filename or die $!;
+
+my @linesInFile;
+while (<FILE>) { 
+    push(@linesInFile, $_);
+}
+
+close FILE;
+
+
+my @fixedproddec = FixProductionsAndDecays(\@linesInFile);
+my @fixedphase = Fixepem(\@linesInFile);
+my @fixedpipiswave = FixPiPiSWave(\@linesInFile);
+my @fixedkpiswave = FixKPiSWaveI12(\@linesInFile);
+push(@fixedkpiswave, FixKPiSWaveI32(\@linesInFile));
+my @fixedgenerickmatrix = FixGenericKMatrix(\@linesInFile);
+
+my $numfixedproddec = @fixedproddec;
+my $numfixedphase = @fixedphase;
+my $numfixedpipiswave = @fixedpipiswave;
+my $numfixedkpiswave = @fixedkpiswave;
+my $numfixedgenerickmatrix = @fixedgenerickmatrix;
+
+my $numfixed = $numfixedproddec + $numfixedphase + $numfixedgenerickmatrix + 
+               $numfixedpipiswave + $numfixedkpiswave;
+
+print "\n\n";
+print "######################################################\n";
+print "##                                                  ##\n";
+print "##  Output from the pawianfixparametersEpEm script  ##\n";
+print "##                                                  ##\n";
+print "######################################################\n";
+print "#\n"; 
+print "# Fixing $numfixed parameters\n";
+print "#\n";
+print "# !! The following list is given without any        !!\n";
+print "# !! guarantee and has to be checked by the analyst !!\n";
+print "#\n\n";
+
+my $fixedline;
+
+
+print "############################################\n";
+print "# Fixing one phase of the e+ e- amplidudes #\n";
+print "############################################\n";
+foreach $fixedline (@fixedphase){
+    print "mnParFix = ".$fixedline."\n";
+}
+
+
+print "\n";
+print "###########################################################\n";
+print "# Fixing $numfixedproddec production and decay parameters #\n";
+print "###########################################################\n";
+foreach $fixedline (@fixedproddec){
+    print "mnParFix = ".$fixedline."\n";
+}
+
+print "\n";
+print "##################################################\n";
+print "# Fixing $numfixedpipiswave PiPiS-wave parameters\n";
+print "##################################################\n";
+foreach $fixedline (@fixedpipiswave){
+    print "mnParFix = ".$fixedline."\n";
+}
+
+print "\n";
+print "##################################################\n";
+print "# Fixing $numfixedkpiswave KPiS-wave parameters\n";
+print "##################################################\n";
+foreach $fixedline (@fixedkpiswave){
+    print "mnParFix = ".$fixedline."\n";
+}
+
+print "\n";
+print "##################################################\n";
+print "# Fixing $numfixedgenerickmatrix generic k-matrix parameters\n";
+print "##################################################\n";
+foreach $fixedline (@fixedgenerickmatrix){
+    print "mnParFix = ".$fixedline."\n";
+}
+
+sub FixProductionsAndDecays
+{
+    my @linesInFile = @{$_[0]};
+    my @fixedParams;
+    my @fixedSources;
+
+    foreach my $currentLine (@linesInFile){
+
+	if($currentLine =~ /Matrix/){
+	    next;
+	}
+	if($currentLine =~ /pipiS/){
+	    next;
+	}
+	if($currentLine =~ /KPiS/){
+	    next;
+	}
+	if($currentLine =~ /(.*_J.*P.*C.*To.*).*/){
+	    next;
+	} 
+
+	if($currentLine =~ /(.*_(.*)To(.*))Mag/ ){
+	    my $currentSource = $2;
+	    my $currentDest = $3;
+	    my $found = 0;
+
+	    if(grep $_ eq $currentSource, @fixedSources){
+		$found = 1;
+	    }
+
+	    # Fix the priority system instead of this one
+	    elsif(!($currentDest eq $prio)){	    
+		foreach my $alternative(@linesInFile){
+		    if( $alternative =~ /(.*_(.*)To(.*))Mag.*/ ){
+			if(($currentSource eq $2) && ($prio eq $3)){
+			    $found = 1;
+			    last;
+			}
+		    }
+		}
+	    }
+	    if($found == 0){
+		push(@fixedSources, $currentSource);
+		push(@fixedParams, $1."Mag");
+		push(@fixedParams, $1."Phi");
+	    }
+	}
+    }
+
+    return @fixedParams;
+}
+
+
+
+sub Fixepem
+{
+    my @linesInFile = @{$_[0]};
+    my @fixedParams;
+    my @fixedSources;
+
+    my $epemIsFixed=0;
+
+    foreach my $currentLine (@linesInFile){
+
+        # Fix epem phases
+        if($currentLine =~ /(.*_J.*P.*C.*To.*Phi).*/ && $epemIsFixed ==0){
+	    push(@fixedParams, $1);
+	    $epemIsFixed = 1;
+        }
+    }
+    
+
+    if(!$epemIsFixed){
+        die "Failed to fix epem phases.";
+    }
+
+    return @fixedParams;
+}
+
+sub FixPiPiSWave
+{
+    my @linesInFile = @{$_[0]};
+    my @fixedParams;
+    my @fixedSources;
+
+    my $poleIsFixed=0;
+    my $decayIsFixed=0;
+
+    foreach my $currentLine (@linesInFile){
+	
+	if( $currentLine =~ /(pipiS.*b_pole1.*)Mag.*/ ){
+	    my $currentSource = $1;
+
+	    if(!(grep $_ eq $currentSource, @fixedSources) && $poleIsFixed == 0 ){
+		push(@fixedSources, $currentSource);
+		push(@fixedParams, $1."Mag");
+		push(@fixedParams, $1."Phi");
+		$poleIsFixed = 1;
+	    }
+
+	}
+	elsif( $currentLine =~ /(pipiS.*S0).*/){
+	    push(@fixedParams, $1);
+	}
+	elsif( $currentLine =~ /(.*pipiS.*To.*)Mag/ && $decayIsFixed==0){
+	    push(@fixedParams, $1."Mag");
+	    push(@fixedParams, $1."Phi");
+	    $decayIsFixed=1;
+	}
+    }
+
+    return @fixedParams;
+}
+
+
+
+
+sub FixKPiSWaveI12
+{
+    my @linesInFile = @{$_[0]};
+    my @fixedParams;
+
+    my $poleIsFixed=0;
+    foreach my $currentLine (@linesInFile){
+
+	if( $currentLine =~ /(KPiS.*b_pole1.*)Mag.*/ && $poleIsFixed == 0 ){
+	    push(@fixedParams, $1."Mag");
+	    push(@fixedParams, $1."Phi");
+	    $poleIsFixed = 1;
+	}
+    }
+    
+    return @fixedParams;
+}
+
+
+
+sub FixKPiSWaveI32
+{
+    my @linesInFile = @{$_[0]};
+    my @fixedParams;
+
+    foreach my $currentLine (@linesInFile){
+
+	if($currentLine =~ /(.*_(.*)ToKpiS32_.*)Mag.*/ ){
+	    my $currentSource = $2;
+	    push(@fixedParams, $1."Mag");
+	    push(@fixedParams, $1."Phi");
+	}
+    }
+    
+    return @fixedParams;
+}
+
+
+
+sub FixGenericKMatrix
+{
+    my @linesInFile = @{$_[0]};
+    my @fixedParams;
+    my @fixedJPCKMat;
+    my $poleIsFixed=0;
+    my @fixedKMatrixPoles; 
+
+    foreach my $currentLine (@linesInFile){
+
+	if( $currentLine =~ /((.*Matrix.*b_).*)Mag/){
+	    if(grep $_ eq $2, @fixedKMatrixPoles){
+		next;
+	       }
+    	    push(@fixedParams, $1."Mag");
+	    push(@fixedParams, $1."Phi");
+	    push(@fixedKMatrixPoles, $2);	    
+	}
+
+#	elsif($currentLine =~ /(.*_(.*To.*Matrix.*))Mag/){
+	elsif($currentLine =~ /(.*(_.*Matrix.*To.*))Mag/){
+	    if(grep $_ eq $2, @fixedJPCKMat){
+		next;
+	    }	    
+
+	    push(@fixedParams, $1."Mag");
+	    push(@fixedParams, $1."Phi");
+	    push(@fixedJPCKMat, $2);
+	}
+	
+    }
+    return @fixedParams;
+}
+
+
+
+sub Help
+{
+
+    print STDERR <<INLINE_LITERAL_TEXT;
+
+NAME
+    pawianfixparameters - automatically fix parameters for PAWIAN Partial Wave Analyses
+
+SYNOPSIS
+    pawianfixparameters -f <file> [-p <priority>]
+
+OPTIONS
+
+    -h | --help     Show this help
+    -f | --file     Select the file containing the fit parameters
+    -p | --prio     Set a priority production system that becomes fixed
+                    rather than another system
+
+EXAMPLES
+
+    pawianfixparameters -f defaultparams.dat
+
+    pawianfixparameters -f myparams.dat -prio "K*K" 
+    #Fixes K*K amplitudes instead of K*(1430)K amplitudes
+
+
+INLINE_LITERAL_TEXT
+}

Utils/PawianConstants.hh

+#pragma once
+
+#include <math.h>
+
+namespace PawianConstants {
+  //mass of pi+- (GeV/c^2)                      
+  const double mPi              = 0.13957018; 
+  //mass of pi0 (GeV/c^2)
+  const double mPi0             = 0.1349766;
+
+  //square of pi+- mass 
+  const double mPiSq            = mPi*mPi;
+  //square of pi0 mass 
+  const double mPi0Sq           = mPi0*mPi0;
+
+  //square of pi 
+  const double pi           = M_PI;
+
+  //rad to deg 
+  const double radToDeg = 180.0/M_PI;
+}