added K-matrix with background terms

ConfigParser/KMatrixParser.cc

@@ -38,6 +38,10 @@ KMatrixParser::KMatrixParser(std::string& path)
       , _projection("")
       , _useBarrierFactors(false)
       , _orbitalMom(0)
+      ,_orderKMatrixBackground(-1)
+      ,_useAdler0(false)
+      ,_s0Adler(1.)
+      ,_snormAdler(1.)
       ,_config(new po::options_description("Configuration file options"))
 
      {
@@ -51,6 +55,10 @@ KMatrixParser::KMatrixParser(std::string& path)
       ("projection",po::value< string>(&_projection),  "projection of the P-vector via pair of decay particles")
       ("useBarrierFactors",po::value<bool>(&_useBarrierFactors),  "calculation with or without barrier factors")
       ("orbitalMomentum",po::value<unsigned int>(&_orbitalMom),  "orbital momentum of the decay")
+      ("orderKMatrixBackground",po::value<int>(&_orderKMatrixBackground),  "order of the K-matrix background")
+      ("useAdler0",po::value<bool>(&_useAdler0),  "use adler0 term")
+      ("s0Adler",po::value<double>(&_s0Adler),  "s0Adler parameter")
+      ("snormAdler",po::value<double>(&_snormAdler),  "snormAdler parameter")
        ;
 
 
@@ -86,7 +94,11 @@ bool KMatrixParser::parseCommandLine()
 	      << "number of poles: " << _noOfPoles  << "\n\n"
 	      << "\nprojection: " << _projection << "\n\n"
 	      << "with barrier factors: " << _useBarrierFactors  << "\n\n"
-	      << "orbital momentum: " << _orbitalMom  << "\n\n" 
+	      << "orbital momentum: " << _orbitalMom  << "\n\n"
+	      << "order of the K-matrix background: " << _orderKMatrixBackground << "\n\n"
+	      << "use Adler0 term: " << _useAdler0 << "\n\n"
+	      << "s0Adler: " << _s0Adler << "\n\n"
+	      << "snormAdler: " << _snormAdler << "\n\n" 
 	      << endl;
     
     std::cout << "g-factors are defined as follows:" << std::endl;

ConfigParser/KMatrixParser.hh

@@ -58,6 +58,10 @@ public:
   const std::string projection() const { return _projection; }
   const int useBarrierFactors() const {return _useBarrierFactors;}
   const unsigned int orbitalMom() const {return _orbitalMom;}
+  const int orderBg() const {return _orderKMatrixBackground;}
+  const bool useAdler() const {return _useAdler0;}
+  const double s0Adler() const {return _s0Adler;}
+  const double snormAdler() const {return _snormAdler;}
  
 protected:
   //  virtual bool parseCommandLine(int argc,char **argv); 
@@ -70,6 +74,10 @@ protected:
   std::string _projection;
   bool _useBarrierFactors;
   unsigned int _orbitalMom;
+  int _orderKMatrixBackground;
+  bool _useAdler0;
+  double _s0Adler;
+  double _snormAdler;
  
   po::options_description* _config;
 };

Examples/Tutorial/LineShapes/KPiSWaveTMatrix.cc

@@ -80,8 +80,9 @@ KPiSWaveTMatrix::KPiSWaveTMatrix() :
     theTMatrix32->evalMatrix(mass);        
 
     complex<double> currentVal=(*theTMatrix12)(0,0)+0.5*(*theTMatrix32)(0,0);
-    _KPiAmpRealH1->Fill(mass, sqrt(thePhpVecs[0]->factor(mass4Vec.M()).real())*currentVal.real());    
-    _KPiAmpImagH1->Fill(mass, sqrt(thePhpVecs[0]->factor(mass4Vec.M()).real())*currentVal.imag());   
+    //    complex<double> currentVal=(*theTMatrix12)(0,0);
+       _KPiAmpRealH1->Fill(mass, sqrt(thePhpVecs[0]->factor(mass4Vec.M()).real())*currentVal.real());    
+       _KPiAmpImagH1->Fill(mass, sqrt(thePhpVecs[0]->factor(mass4Vec.M()).real())*currentVal.imag());
   }
 }
 

Examples/Tutorial/LineShapes/TMatrixGeneral.cc

@@ -36,6 +36,7 @@
 #include "PwaDynamics/KPole.hh"
 #include "PwaDynamics/KPoleBarrier.hh"
 #include "PwaDynamics/KMatrixRel.hh"
+#include "PwaDynamics/KMatrixRelBg.hh"
 #include "PwaDynamics/AbsPhaseSpace.hh"
 #include "PwaDynamics/PhaseSpaceIsobar.hh"
 #include "ConfigParser/KMatrixParser.hh"
@@ -116,6 +117,12 @@ TMatrixGeneral::TMatrixGeneral(std::string pathToConfigParser) :
     currentSqrT11H1->SetYTitle("|T_{11}|^{2}");
     currentSqrT11H1->SetXTitle("mass/GeV");
     _SqrT11H1Vec.push_back(currentSqrT11H1);
+
+    std::string currentphpKey="phase space factor"+key;
+    TH1F* currentphpH1=new TH1F(currentphpKey.c_str(), currentphpKey.c_str(), _noOfSteps-1, _massMin, _massMax);
+    currentphpH1->SetYTitle("#rho");
+    currentphpH1->SetXTitle("mass/GeV");
+    _phpH1Vec.push_back(currentphpH1);
   }
 
   DebugMsg << "_massMin: "<< _massMin
@@ -130,8 +137,8 @@ TMatrixGeneral::TMatrixGeneral(std::string pathToConfigParser) :
     for(unsigned int i=0; i<_gFactorNames.size(); ++i){
       complex<double> currentRho=_phpVecs.at(i)->factor(mass);
       //      std::cout << mass << "(" << i << ")\t"<< (*_tMatr)(i,i) << std::endl;
-      _AmpRealH1Vec.at(i)->Fill(mass, (*_tMatr)(i,i).real());
-      _AmpImagH1Vec.at(i)->Fill(mass, (*_tMatr)(i,i).imag());
+      _AmpRealH1Vec.at(i)->Fill(mass, sqrt(currentRho.real())*(*_tMatr)(i,i).real());
+      _AmpImagH1Vec.at(i)->Fill(mass, sqrt(currentRho.real())*(*_tMatr)(i,i).imag());
       
       _ArgandH2Vec.at(i)->Fill( currentRho.real()*(*_tMatr)(i,i).real(), currentRho.real()*(*_tMatr)(i,i).imag());
       double currentphase=360.*atan2((*_tMatr)(i,i).imag(),(*_tMatr)(i,i).real()) / 3.1415;
@@ -141,7 +148,9 @@ TMatrixGeneral::TMatrixGeneral(std::string pathToConfigParser) :
       complex<double> S00_rel=complex<double>(1.,0.)+2.*complex<double>(0.,1.)*currentRho.real()*(*_tMatr)(i,i);
       _ElasticityH1Vec.at(i)->Fill(mass, sqrt(norm(S00_rel)));
 
-      _SqrT11H1Vec.at(i)->Fill(mass,currentRho.real()*norm((*_tMatr)(i,i))); 
+      _SqrT11H1Vec.at(i)->Fill(mass,currentRho.real()*norm((*_tMatr)(i,i)));
+
+      _phpH1Vec.at(i)->Fill(mass, sqrt(norm(currentRho))); 
     }    
   }
 }
@@ -243,7 +252,30 @@ void TMatrixGeneral::init(){
     _kPoles.push_back(currentPole);
   }
 
-  _kMatr=std::shared_ptr<KMatrixRel>(new KMatrixRel(_kPoles,_phpVecs ));
+  int orderBg=_kMatrixParser->orderBg();
+  if(orderBg<0) _kMatr=std::shared_ptr<KMatrixRel>(new KMatrixRel(_kPoles,_phpVecs ));
+  else{
+    bool withAdler=_kMatrixParser->useAdler();
+    _kMatr=std::shared_ptr<KMatrixRel>(new KMatrixRelBg(_kPoles,_phpVecs, orderBg, withAdler));
+    _kMatr->updateBgTerms(0, 0, 0, 0.79299);
+    _kMatr->updateBgTerms(0, 0, 1, 0.15040);
+    _kMatr->updateBgTerms(0, 1, 1, 0.17054);
+    if(orderBg>0){
+      _kMatr->updateBgTerms(1, 0, 0, -0.15099);
+      _kMatr->updateBgTerms(1, 0, 1, -0.038266);
+      _kMatr->updateBgTerms(1, 1, 1, -0.0219);
+    }
+    if(orderBg>1){
+      _kMatr->updateBgTerms(2, 0, 0, 0.00811);
+      _kMatr->updateBgTerms(2, 0, 1, 0.0022596);
+      _kMatr->updateBgTerms(2, 1, 1, 0.00085655);
+    }
+
+    if(withAdler){
+      _kMatr->updates0Adler(_kMatrixParser->s0Adler());
+      _kMatr->updatesnormAdler(_kMatrixParser->snormAdler());
+    }
+  }
   _tMatr=std::shared_ptr<TMatrixRel>(new TMatrixRel(_kMatr));
 
 

Examples/Tutorial/LineShapes/TMatrixGeneral.hh

@@ -90,6 +90,7 @@ private:
   std::vector<TH2F*> _PhaseH2Vec;
   std::vector<TH1F*> _ElasticityH1Vec;
   std::vector<TH1F*> _SqrT11H1Vec;
+  std::vector<TH1F*> _phpH1Vec;
 
   void init();
 };

PwaDynamics/KMatrixBase.cc

@@ -25,52 +25,52 @@
 #include "PwaDynamics/KPole.hh"
 #include "qft++/relativistic-quantum-mechanics/Utils.hh"
 #include "qft++/matrix/IdentityMatrix.hh"
+#include "ErrLogger/ErrLogger.hh"
 
 KMatrixBase::KMatrixBase(vector<std::shared_ptr<KPole> > Kpoles, vector<std::shared_ptr<AbsPhaseSpace> > phpVecs) :
   Matrix< complex<double> >::Matrix(int(phpVecs.size()), int(phpVecs.size()))
   ,_KPoles(Kpoles)
   ,_phpVecs(phpVecs)
+  ,_orderBg(0)
+  ,_s0Adler(0.)
+  ,_snormAdler(1.)
  {
-   _bgTerms.resize(phpVecs.size());
-   for(unsigned int i=0; i<phpVecs.size(); ++i){
-     std::vector<double> currentbgVec;
-     currentbgVec.resize(phpVecs.size());
-     for(unsigned int j=0; j<phpVecs.size(); ++j){
-       currentbgVec[j]=0.;
-     }
-     _bgTerms[i]=currentbgVec;
-   }
  }
 
 KMatrixBase::KMatrixBase(vector<std::shared_ptr<AbsPhaseSpace> > phpVecs, int numCols, int numRows) :
   Matrix< complex<double> >::Matrix(numCols, numRows)
   ,_phpVecs(phpVecs)
+  ,_orderBg(0)
+  ,_s0Adler(0.)
+  ,_snormAdler(1.)
  {
-   _bgTerms.resize(numRows);
-   for(int i=0; i<numRows; ++i){
-     std::vector<double> currentbgVec;
-     currentbgVec.resize(numCols);
-     for(int j=0; j<numCols; ++j){
-       currentbgVec[j]=0.;
-     }
-     _bgTerms[i]=currentbgVec;
-   }
  }
 
 KMatrixBase::KMatrixBase(int numCols, int numRows) :
   Matrix< complex<double> >::Matrix(numCols, numRows)
+  ,_orderBg(0)
  {
-   _bgTerms.resize(numRows);
-   for(int i=0; i<numRows; ++i){
-     std::vector<double> currentbgVec;
-     currentbgVec.resize(numCols);
-     for(int j=0; j<numCols; ++j){
-       currentbgVec[j]=0.;
-     }
-     _bgTerms[i]=currentbgVec;
-   }
  }
 
 KMatrixBase::~KMatrixBase(){
 }
 
+void KMatrixBase::updateBgTerms(unsigned int order, unsigned int row,  unsigned int column, double theVal){
+  if(row>column){
+    Alert << "K-matrix is symmetric; for update row<=column" << endmsg;
+    exit(0);
+  }
+  if(order>_orderBg){
+    Alert << "background parameter for order " << order << " not available!!!" << endmsg;
+    exit(0);
+  }
+  if(row>=NumRows() || column>=NumCols()){
+    Alert << "row " << row << " or column " << column 
+	  << " >= NumRows " << NumRows() << " or  >= NumCols " << NumCols() << endmsg;
+    exit(0);
+  }
+
+  _bgTerms.at(order).at(row).at(column)=theVal;
+  _bgTerms.at(order).at(column).at(row)=theVal;
+}
+

PwaDynamics/KMatrixBase.hh

@@ -57,11 +57,17 @@ public:
   virtual void evalMatrix(const double mass) {return;}
   virtual vector<std::shared_ptr<AbsPhaseSpace> > phaseSpaceVec() {return _phpVecs;}
   virtual vector<std::shared_ptr<KPole> > kpoles() {return _KPoles;}
+  virtual void updateBgTerms(unsigned int order, unsigned int row,  unsigned int column, double theVal);
+  virtual void updates0Adler(double s0Adler) {_s0Adler=s0Adler;}
+  virtual void updatesnormAdler(double snormAdler) {_snormAdler=snormAdler;}
 
 protected:
   vector<std::shared_ptr<KPole> > _KPoles;
   vector<std::shared_ptr<AbsPhaseSpace> > _phpVecs;
-  vector< vector<double> > _bgTerms; 
+  vector< vector< vector<double> > > _bgTerms;
+  unsigned int _orderBg;
+  double _s0Adler;
+  double _snormAdler;
 };
 //_____________________________________________________________________________
 

PwaDynamics/KMatrixKPiSFocus.cc

@@ -67,7 +67,7 @@ KMatrixKPiSFocus::~KMatrixKPiSFocus(){
 }
 
 void KMatrixKPiSFocus::evalMatrix(const double mass){
-
+ 
   Matrix< complex<double> > theKMatrix(NumRows(), NumRows());
   vector<std::shared_ptr<KPole> >::iterator it;
   for (it =_KPoles.begin(); it != _KPoles.end(); ++it){

PwaDynamics/KMatrixRel.cc

@@ -47,7 +47,7 @@ void KMatrixRel::evalMatrix(const double mass){
 
   for (int i=0; i<theKMatrix.NumRows(); ++i){
     for (int j=0; j<theKMatrix.NumCols(); ++j){
-      this->operator()(i,j)=theKMatrix(i,j)+_bgTerms.at(i).at(j);
+      this->operator()(i,j)=theKMatrix(i,j);
     }
   }
 }

PwaDynamics/KMatrixRelBg.cc

+//************************************************************************//
+//									  //
+//  Copyright 2013 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/>.	  //
+//									  //
+//************************************************************************//
+
+#include "PwaDynamics/KMatrixRelBg.hh"
+#include "PwaDynamics/KPole.hh"
+#include "PwaDynamics/AbsPhaseSpace.hh"
+#include "qft++/relativistic-quantum-mechanics/Utils.hh"
+#include "qft++/matrix/IdentityMatrix.hh"
+
+KMatrixRelBg::KMatrixRelBg(vector<std::shared_ptr<KPole> > Kpoles, vector<std::shared_ptr<AbsPhaseSpace> > phpVecs, unsigned int orderBg, bool withAdler) :
+  KMatrixRel(Kpoles, phpVecs)
+  ,_withAdler(withAdler)
+ {
+   _orderBg=orderBg;
+   _bgTerms.resize(_orderBg+1);
+   for(unsigned int i=0; i<=_orderBg; ++i){
+     _bgTerms.at(i).resize(phpVecs.size());
+     for(unsigned int j=0; j<phpVecs.size(); ++j){
+       _bgTerms.at(i).at(j).resize(phpVecs.size());
+       for(unsigned int k=0; k<phpVecs.size(); ++k){
+	 _bgTerms.at(i).at(j).at(k)=0.;
+       }
+     }
+   }
+ }
+
+
+KMatrixRelBg::~KMatrixRelBg(){
+}
+
+void KMatrixRelBg::evalMatrix(const double mass){
+
+  double adlerTerm=1.;
+  double s_hat=mass*mass;
+  if(_withAdler){
+    adlerTerm=(mass*mass-_s0Adler)/_snormAdler;
+    s_hat=mass*mass/_snormAdler-1.;
+  }
+
+  Matrix< complex<double> > theKMatrix(NumRows(), NumRows());
+  vector<std::shared_ptr<KPole> >::iterator it;
+  for (it =_KPoles.begin(); it != _KPoles.end(); ++it){
+    (*it)->evalMatrix(mass);
+    theKMatrix += *(*it);
+  }
+
+  for (int i=0; i<NumRows(); ++i){
+    for (int j=i; j<NumRows(); ++j){
+      complex<double> currentBg(0.,0.);
+      for (unsigned int k=0; k<=_orderBg; ++k){
+	currentBg+=complex<double>(_bgTerms.at(k).at(i).at(j)*pow(s_hat,k), 0.);
+      }
+      this->operator()(i,j)=theKMatrix(i,j)+currentBg;
+      this->operator()(i,j)*=adlerTerm;
+      this->operator()(j,i)=this->operator()(i,j);
+    }
+  }
+}

PwaDynamics/KMatrixRelBg.hh

+//************************************************************************//
+//									  //
+//  Copyright 2013 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/>.	  //
+//									  //
+//************************************************************************//
+
+// KMatrixRelBg class definition file. -*- C++ -*-
+// Copyright 2013 Bertram Kopf
+
+#pragma once 
+
+//_____________________________________________________________________________
+// @file KMatrixRelBg.h
+//_____________________________________________________________________________
+
+#include "PwaDynamics/KMatrixRel.hh"
+#include <iostream>
+#include <vector>
+#include <memory>
+
+class KPole;
+class AbsPhaseSpace;
+
+using namespace std;
+//_____________________________________________________________________________
+//_____________________________________________________________________________
+
+class KMatrixRelBg : public KMatrixRel {
+
+public:
+
+  /// Constructor 
+  KMatrixRelBg(vector<std::shared_ptr<KPole> > Kpoles, vector<std::shared_ptr<AbsPhaseSpace> > phpVecs, unsigned int orderBg, bool withAdler); 
+
+  /// Destructor
+  virtual ~KMatrixRelBg();
+
+  virtual void evalMatrix(const double mass);
+
+protected:
+  bool _withAdler;
+};
+//_____________________________________________________________________________
+
+