added missing class TMatrixDynamics

PwaUtils/TMatrixDynamics.cc

+//************************************************************************//
+//									  //
+//  Copyright 2017 Bertram Kopf (bertram@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/>.	  //
+//									  //
+//************************************************************************//
+
+// TMatrixDynamics class definition file. -*- C++ -*-
+// Copyright 2017 Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+
+#include "PwaUtils/TMatrixDynamics.hh"
+#include "PwaUtils/XdecAmpRegistry.hh"
+#include "PwaUtils/AbsDecay.hh"
+#include "PwaUtils/AbsXdecAmp.hh"
+#include "PwaUtils/GlobalEnv.hh"
+#include "ErrLogger/ErrLogger.hh"
+#include "Particle/Particle.hh"
+#include "Particle/ParticleTable.hh"
+#include "ConfigParser/KMatrixParser.hh"
+#include "PwaDynamics/KMatrixRel.hh"
+#include "PwaDynamics/KMatrixRelBg.hh"
+#include "PwaDynamics/TMatrixRel.hh"
+#include "PwaDynamics/KPole.hh"
+#include "PwaDynamics/KPoleBarrier.hh"
+#include "PwaDynamics/PPole.hh"
+#include "PwaDynamics/PPoleBarrier.hh"
+#include "PwaDynamics/AbsPhaseSpace.hh"
+#include "PwaDynamics/PhaseSpaceFactory.hh"
+#include "FitParams/AbsPawianParameters.hh"
+
+TMatrixDynamics::TMatrixDynamics(std::string& name, std::vector<Particle*>& fsParticles, Particle* mother, std::string& pathToConfigParser) :
+  AbsDynamics(name, fsParticles, mother)
+  ,_kMatName("")
+  ,_prodProjectionIndex(0)
+  ,_decProjectionIndex(0)
+  , _orderKMatBg(-1)
+  ,_withKMatAdler(false)
+  ,_currentAdler0(0.)
+  ,_kMatrixParser(new KMatrixParser(pathToConfigParser))
+{
+  init();
+  _isLdependent=true;
+}
+
+TMatrixDynamics::~TMatrixDynamics()
+{
+}
+
+complex<double> TMatrixDynamics::eval(EvtData* theData, AbsXdecAmp* grandmaAmp, Spin OrbMom){
+  _tMatr->evalMatrix(theData->DoubleMassId.at(_dynId));
+  return (*_tMatr)(_prodProjectionIndex,_decProjectionIndex);
+}
+
+void TMatrixDynamics::fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar){
+  //pole positions
+  for(unsigned int i=0; i<_currentPoleMasses.size(); ++i){
+    double valMass=_currentPoleMasses.at(i);
+    double errMass=0.02;
+    double minMass=valMass-5.*errMass;
+    if(minMass<0.) minMass=0.;
+    double maxMass=valMass+5.*errMass;
+    fitPar->Add(_poleNames.at(i)+"Mass", valMass, errMass);
+    fitPar->SetLimits(_poleNames.at(i)+"Mass", minMass, maxMass);
+  }
+
+  //g-factors
+  for(unsigned int i=0; i<_poleNames.size(); ++i){
+    std::vector<double> currentgFactorVec=_currentgFactorMap.at(i);
+    for(unsigned int j=0; j<currentgFactorVec.size(); ++j){
+      std::string currentName=_poleNames.at(i)+_gFactorNames.at(j)+"gFactor";
+      double currentError=currentgFactorVec.at(j)/3.;
+      if (currentError<1.e-5) currentError=0.01;
+      fitPar->Add(currentName, currentgFactorVec.at(j), currentError);
+    }
+  }
+
+  //k-matrix bg-terms
+  if(_orderKMatBg>=0){
+    for(unsigned int i=0; i<=fabs(_orderKMatBg); ++i){
+      for(unsigned int j=0; j<_phpVecs.size(); ++j){
+  	for(unsigned int k=j; k<_phpVecs.size(); ++k){
+  	  std::string currentName=_bgTermNames.at(i).at(j).at(k);
+  	  fitPar->Add(currentName, _currentBgTerms.at(i).at(j).at(k), fabs(_currentBgTerms.at(i).at(j).at(k))+0.3);
+  	}
+      }
+    }
+
+    //Adler-term
+    if(_withKMatAdler){
+      fitPar->Add("s0"+_kMatName, _currentAdler0, fabs(_currentAdler0)+0.2);
+    }
+  }
+}
+
+void TMatrixDynamics::fillParamNameList(){
+  _paramNameList.clear();
+
+  std::map<std::string, std::vector<std::string> >::iterator itNameMap;  
+  //pole positions
+  for(unsigned int i=0; i<_currentPoleMasses.size(); ++i){
+    _paramNameList.push_back(_poleNames.at(i)+"Mass");
+    for(itNameMap=_paramNameListMap.begin(); itNameMap!=_paramNameListMap.end(); ++itNameMap){
+      itNameMap->second.push_back(_poleNames.at(i)+"Mass");
+    }
+  }
+  
+  //g-factors
+  for(unsigned int i=0; i<_poleNames.size(); ++i){
+    std::vector<double> currentgFactorVec=_currentgFactorMap.at(i);
+    for(unsigned int j=0; j<currentgFactorVec.size(); ++j){
+      std::string currentName=_poleNames.at(i)+_gFactorNames.at(j)+"gFactor";
+      _paramNameList.push_back(currentName);
+      for(itNameMap=_paramNameListMap.begin(); itNameMap!=_paramNameListMap.end(); ++itNameMap){
+	itNameMap->second.push_back(currentName);
+      }
+    }
+  }
+  
+  //k-matrix bg-terms
+  if(_orderKMatBg>=0){
+    for(unsigned int i=0; i<=fabs(_orderKMatBg); ++i){
+      for(unsigned int j=0; j<_phpVecs.size(); ++j){
+  	for(unsigned int k=j; k<_phpVecs.size(); ++k){
+  	  std::string currentName=_bgTermNames.at(i).at(j).at(k);
+	  _paramNameList.push_back(currentName);
+          for(itNameMap=_paramNameListMap.begin(); itNameMap!=_paramNameListMap.end(); ++itNameMap){
+	    itNameMap->second.push_back(currentName);
+	  }
+  	}
+      }
+    }
+    //Adler-term
+    if(_withKMatAdler){
+      _paramNameList.push_back("s0"+_kMatName);
+      for(itNameMap=_paramNameListMap.begin(); itNameMap!=_paramNameListMap.end(); ++itNameMap){
+	itNameMap->second.push_back("s0"+_kMatName);
+      }
+    }
+  }
+}
+
+bool TMatrixDynamics::checkRecalculation(std::shared_ptr<AbsPawianParameters> fitParNew, std::shared_ptr<AbsPawianParameters> fitParOld){
+  return true;
+}
+
+void TMatrixDynamics::updateFitParams(std::shared_ptr<AbsPawianParameters> fitPar){
+
+  //pole positions
+  for(unsigned int i=0; i<_currentPoleMasses.size(); ++i){
+    std::string currentPoleName=_poleNames.at(i)+"Mass";
+    double currentPoleMass=fitPar->Value(currentPoleName);
+    _currentPoleMasses.at(i)=currentPoleMass;
+    _kPoles.at(i)->updatePoleMass(currentPoleMass);
+  }
+ 
+ //g-factors
+  for(unsigned int i=0; i<_poleNames.size(); ++i){
+    std::vector<double>& currentgFactorVec=_currentgFactorMap.at(i);
+    for(unsigned int j=0; j<currentgFactorVec.size(); ++j){
+      std::string currentName=_poleNames.at(i)+_gFactorNames.at(j)+"gFactor";
+      currentgFactorVec.at(j)=fitPar->Value(currentName);
+    }
+    _kPoles.at(i)->updategFactors(currentgFactorVec);
+  }
+    
+  //k-matrix bg-terms
+  if(_orderKMatBg>=0){
+    for(unsigned int i=0; i<=fabs(_orderKMatBg); ++i){
+      for(unsigned int j=0; j<_phpVecs.size(); ++j){
+	for(unsigned int k=j; k<_phpVecs.size(); ++k){
+	  double newVal=fitPar->Value(_bgTermNames.at(i).at(j).at(k));
+	  _currentBgTerms.at(i).at(j).at(k)=newVal;
+	  _kMatr->updateBgTerms(i,j,k,newVal);
+	}
+      }
+    }
+
+    //Adler-term
+    if(_withKMatAdler){
+      _currentAdler0=fitPar->Value("s0"+_kMatName);
+      _kMatr->updates0Adler(_currentAdler0);
+    }
+  }
+}
+
+
+
+void TMatrixDynamics::init(){
+  _kMatName=_kMatrixParser->keyName();
+  _orderKMatBg=_kMatrixParser->orderBg();
+  _withKMatAdler=_kMatrixParser->useAdler();  
+
+  std::vector<std::string> poleNameAndMassVecs=_kMatrixParser->poles();
+  std::vector<std::string>::iterator itString;
+  for (itString=poleNameAndMassVecs.begin(); itString!=poleNameAndMassVecs.end(); ++itString){
+    std::istringstream poleIString(*itString);
+    std::string currentPoleName;
+    std::string currentPoleMassStr;
+    poleIString >> currentPoleName >> currentPoleMassStr;
+
+    std::istringstream currentPoleMassiStr(currentPoleMassStr);
+    double currentValue;
+    if(!(currentPoleMassiStr >> currentValue)){
+      Alert << "cannot convert " << currentPoleMassStr << " to a double value" << endmsg;
+      exit(0);
+    }
+
+    _currentPoleMasses.push_back(currentValue);
+    _poleNames.push_back(currentPoleName);
+  }
+
+  if(_currentPoleMasses.size()!= _kMatrixParser->noOfPoles()){
+    Alert << "number of poles != number of pole masses" << endmsg;
+    exit(0);
+  }
+
+  const std::vector<std::string> gFacStringVec=_kMatrixParser->gFactors();
+  std::map<std::pair<std::string, std::string>, std::string> phpDescriptionVec=_kMatrixParser->phpDescriptionMap();
+
+  std::vector<std::string>::const_iterator itStrConst;
+  for(itStrConst=gFacStringVec.begin(); itStrConst!=gFacStringVec.end(); ++itStrConst){
+    std::istringstream particles(*itStrConst);
+    std::string firstParticleName;
+    std::string secondParticleName;
+    particles >> firstParticleName >> secondParticleName;
+
+    Particle* firstParticle = GlobalEnv::instance()->particleTable()->particle(firstParticleName);
+    Particle* secondParticle = GlobalEnv::instance()->particleTable()->particle(secondParticleName);
+    if(0==firstParticle || 0==secondParticle){
+      Alert << "particle with name: " << firstParticleName <<" or " << secondParticleName << " doesn't exist in pdg-table" << endmsg;
+      exit(0);
+    }
+
+  std::pair<std::string, std::string> currentParticlePair=make_pair(firstParticleName, secondParticleName);
+    std::string currentPhpDescription = phpDescriptionVec.at(currentParticlePair);
+
+    std::vector<double> fsParticleMasses;
+    fsParticleMasses.push_back(firstParticle->mass());
+    fsParticleMasses.push_back(secondParticle->mass()); 
+
+    std::shared_ptr<AbsPhaseSpace> currentPhp=PhaseSpaceFactory::instance()->getPhpPointer(currentPhpDescription, fsParticleMasses);
+    _phpVecs.push_back(currentPhp);
+
+    std::string gFactorKey=firstParticleName+secondParticleName;    
+    _gFactorNames.push_back(gFactorKey);   
+
+    for (int i=0; i<int(_kMatrixParser->noOfPoles()); ++i){
+      std::string currentPoleName=_poleNames[i];
+      std::string currentgValueStr;
+      if(!(particles >> currentgValueStr)){
+	Alert << "g-factors for pole " << currentPoleName << " does not exist!" << endmsg;
+	exit(0);
+      }
+      std::istringstream currentgValueiStr(currentgValueStr);
+      double currentGValue;
+      if (!(currentgValueiStr >> currentGValue)){
+	Alert << "cannot convert " << currentgValueStr << " to a double value" << endmsg;
+	exit(0);
+      }
+      _currentgFactorMap[i].push_back(currentGValue);
+    }
+  }
+
+
+  std::map<int, std::vector<double> >::iterator itgFac;
+  for (itgFac=_currentgFactorMap.begin(); itgFac!=_currentgFactorMap.end(); ++itgFac){
+    std::vector<double> currentgVector=itgFac->second;
+    std::shared_ptr<KPole> currentPole;
+    if (_kMatrixParser->useBarrierFactors()) currentPole=std::shared_ptr<KPole>(new KPoleBarrier(currentgVector, _currentPoleMasses.at(itgFac->first), _phpVecs, _kMatrixParser->orbitalMom(), _kMatrixParser->useTruncatedBarrierFactors()));
+    else currentPole=std::shared_ptr<KPole>(new KPole(currentgVector, _currentPoleMasses.at(itgFac->first)));
+    _kPoles.push_back(currentPole);
+  }
+
+  if(_orderKMatBg<0) _kMatr=std::shared_ptr<KMatrixRel>(new KMatrixRel(_kPoles,_phpVecs ));
+  else {
+    _currentBgTerms.resize(_orderKMatBg+1);
+    _bgTermNames.resize(_orderKMatBg+1);
+    for(unsigned int i=0; i<= fabs(_orderKMatBg); ++i){
+      _currentBgTerms.at(i).resize(_phpVecs.size());
+      _bgTermNames.at(i).resize(_phpVecs.size());
+      for(unsigned int j=0; j<_phpVecs.size(); ++j){
+	_currentBgTerms.at(i).at(j).resize(_phpVecs.size());
+	_bgTermNames.at(i).at(j).resize(_phpVecs.size());
+	for(unsigned int k=j; k<_phpVecs.size(); ++k){
+	  _currentBgTerms.at(i).at(j).at(k)=0.;
+
+	  std::stringstream keyOrderStrStr;
+	  keyOrderStrStr << i << j << k;
+	  std::string keyOrder=keyOrderStrStr.str();
+	  std::string currentName="bg"+keyOrder+_kMatName;
+	  _bgTermNames.at(i).at(j).at(k)=currentName;
+	}
+      }
+    }
+
+    _kMatr=std::shared_ptr<KMatrixRel>(new KMatrixRelBg(_kPoles,_phpVecs, _orderKMatBg, _withKMatAdler ));
+    if(_withKMatAdler){
+      _currentAdler0=_kMatrixParser->s0Adler();
+      _kMatr->updates0Adler(_currentAdler0);
+      _kMatr->updatesnormAdler(_kMatrixParser->snormAdler());
+    }
+  }
+ 
+  const std::string porjectionParticleNames=_kMatrixParser->projection();
+  std::istringstream projParticles(porjectionParticleNames);
+  std::string firstProjParticleName;
+  std::string secondProjParticleName;
+  projParticles >> firstProjParticleName >> secondProjParticleName;
+  std::string projKey=firstProjParticleName+secondProjParticleName;    
+
+  bool found=false;
+  for(unsigned int i=0; i<_gFactorNames.size();++i){
+    if(projKey==_gFactorNames[i]){
+      _decProjectionIndex=i;
+      found=true;
+    }
+  }
+  if (!found){
+    Alert << "decay projection index for key " << projKey << " not found" << endmsg;
+    exit(0);
+  }
+
+  _tMatr=std::shared_ptr<TMatrixRel>(new TMatrixRel(_kMatr));  
+}
+
+

PwaUtils/TMatrixDynamics.hh

+//************************************************************************//
+//									  //
+//  Copyright 2017 Bertram Kopf (bertram@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/>.	  //
+//									  //
+//************************************************************************//
+
+// TMatrixDynamics class definition file. -*- C++ -*-
+// Copyright 2017 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+#include <memory>
+
+#include "PwaUtils/AbsDynamics.hh"
+
+class AbsXdecAmp;
+class KMatrixRel;
+class TMatrixRel;
+class KPole;
+class KMatrixParser;
+class AbsPhaseSpace;
+class AbsPawianParameters;
+
+class TMatrixDynamics : public AbsDynamics{
+
+public:
+  TMatrixDynamics(std::string& name, std::vector<Particle*>& fsParticles, Particle* mother, std::string& pathToConfigParser);
+  virtual ~TMatrixDynamics();
+
+  virtual std::string type() {return "TMatrixDynamics";}
+  virtual complex<double> eval(EvtData* theData, AbsXdecAmp* grandmaAmp, Spin OrbMom=0);
+  
+  virtual void fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar);
+  virtual void fillParamNameList();
+  virtual bool checkRecalculation(std::shared_ptr<AbsPawianParameters> fitParNew, std::shared_ptr<AbsPawianParameters> fitParOld); 
+  virtual void updateFitParams(std::shared_ptr<AbsPawianParameters> fitPar);
+ 
+protected:
+  std::string _kMatName;
+  int _prodProjectionIndex;
+  int _decProjectionIndex;
+  int _orderKMatBg;
+  bool _withKMatAdler;
+  
+  std::shared_ptr<KMatrixRel> _kMatr;
+  std::shared_ptr<TMatrixRel> _tMatr;
+  std::vector< std::shared_ptr<KPole> > _kPoles;
+  std::vector<std::shared_ptr<AbsPhaseSpace> > _phpVecs;
+  std::vector< std::string> _poleNames;
+  std::vector< std::string> _gFactorNames;
+
+  std::vector<double> _currentPoleMasses;
+  std::map<int, std::vector<double> > _currentgFactorMap;
+  std::vector< std::vector< std::vector<double> > > _currentBgTerms;
+  std::vector< std::vector< std::vector<std::string> > > _bgTermNames;
+  double _currentAdler0;
+
+  std::map<std::string, bool > _recalcMap;  
+  std::shared_ptr<KMatrixParser> _kMatrixParser;
+  std::map<std::string, std::vector<std::string> > _paramNameListMap;
+
+  virtual void init();
+
+private:
+
+};