From d63bdc41a948cb8ee6532973616681b45d2c138c Mon Sep 17 00:00:00 2001
From: Orestis Afedulidis <oafedulidis@pc68.ep1.rub.de>
Date: Tue, 16 Jan 2024 11:19:25 +0100
Subject: [PATCH] reversed lsLoop for loop swap in TensorDecAmps and added
preCaching for daughterAmps(lam1,lam2)
---
PwaUtils/.nfs0000000003d2692b00000057 | 263 ++++++++++++++++++++++++++
PwaUtils/IsobarTensorDecay.cc | 4 +-
PwaUtils/TensorDecAmps.cc | 89 ++++++---
PwaUtils/TensorDecAmps.hh | 8 +-
4 files changed, 337 insertions(+), 27 deletions(-)
create mode 100644 PwaUtils/.nfs0000000003d2692b00000057
diff --git a/PwaUtils/.nfs0000000003d2692b00000057 b/PwaUtils/.nfs0000000003d2692b00000057
new file mode 100644
index 00000000..346fa658
--- /dev/null
+++ b/PwaUtils/.nfs0000000003d2692b00000057
@@ -0,0 +1,263 @@
+//************************************************************************//
+// //
+// 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/>. //
+// //
+//************************************************************************//
+
+// TensorDecAmps class definition file. -*- C++ -*-
+// Copyright 2012 Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+#include <mutex>
+
+#include "PwaUtils/TensorDecAmps.hh"
+#include "qft++/relativistic-quantum-mechanics/Utils.hh"
+#include "PwaUtils/DataUtils.hh"
+#include "PwaUtils/IsobarTensorDecay.hh"
+#include "Particle/Particle.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+TensorDecAmps::TensorDecAmps(std::shared_ptr<IsobarTensorDecay> theDec, ChannelID channelID) :
+ AbsXdecAmp(theDec, channelID)
+ ,_LSs(theDec->LSAmps())
+ ,_factorMag(1.)
+{
+ initialize();
+ if(_LSs.size()>0) _factorMag/=sqrt(_LSs.size());
+}
+
+TensorDecAmps::TensorDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+ AbsXdecAmp(theDec, channelID)
+ ,_factorMag(1.)
+{
+ initialize();
+ if(_LSs.size()>0) _factorMag/=sqrt(_LSs.size());
+}
+
+TensorDecAmps::~TensorDecAmps()
+{
+}
+
+
+complex<double> TensorDecAmps::XdecPartAmp(const Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs,AbsXdecAmp* grandmaAmp){
+
+ Spin lam1Min=-_Jdaughter1;
+ Spin lam1Max= _Jdaughter1;
+ Spin lam2Min=-_Jdaughter2;
+ Spin lam2Max=_Jdaughter2;
+
+ if(fixDaughterNr == 1){
+ lam1Min = lam1Max = lamDec;
+ }
+ else if(fixDaughterNr == 2){
+ lam2Min = lam2Max = lamDec;
+ }
+ else{
+ Alert << "Invalid fixDaughterNr in XdecPartAmp." << endmsg;
+ }
+
+ if(_enabledlamFsDaughter1){
+ lam1Min=lamFs;
+ lam1Max=lamFs;
+ }
+ else if(_enabledlamFsDaughter2){
+ lam2Min=lamFs;
+ lam2Max=lamFs;
+ }
+
+ complex<double> result=lsLoop(grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, false);
+
+ return result;
+}
+
+
+
+
+complex<double> TensorDecAmps::XdecAmp(const Spin& lamX, EvtData* theData, AbsXdecAmp* grandmaAmp){
+
+ complex<double> result(0.,0.);
+ if( fabs(lamX) > _JPCPtr->J) return result;
+
+ short currentSpinIndex=FunctionUtils::spin1IdIndex(_projId,lamX);
+
+ if (!_recalculate){
+ return _cachedAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex);
+ }
+
+ result=lsLoop(grandmaAmp, lamX, theData, _lam1MinProj,
+ _lam1MaxProj,
+ _lam2MinProj,
+ _lam2MaxProj, true);
+
+ if ( _cacheAmps){
+ // _cachedAmpMap[evtNo][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
+ _cachedAmpIdMap[theData->evtNo][_absDyn->grandMaId(grandmaAmp)][currentSpinIndex]=result;
+ }
+ return result;
+}
+
+
+complex<double> TensorDecAmps::lsLoop(AbsXdecAmp* grandmaAmp, Spin lamX, EvtData* theData, Spin lam1Min,
+ Spin lam1Max, Spin lam2Min, Spin lam2Max, bool withDecs) {
+ complex<double> result(0.,0.);
+
+ map<unsigned short, map<Id3StringType, complex<double> > >& currentLS3SpinMap=theData->ComplexLS3Spin.at(_decay->nameId());
+ std::vector< std::shared_ptr<const LScomb> >::iterator it;
+ for (it=_LSs.begin(); it!=_LSs.end(); ++it){
+ map<Id3StringType, complex<double> >& current3SpinMap = currentLS3SpinMap.at((*it)->idnumberLS);
+ complex<double> theMagExpi=_currentParamMagExpi.at(*it);
+
+ complex<double> tmpResult(0.,0.);
+ for(Spin lambda1=lam1Min; lambda1<=lam1Max; ++lambda1){
+ for(Spin lambda2=lam2Min; lambda2<=lam2Max; ++lambda2){
+ Id3StringType IdLamXLam1Lam2=FunctionUtils::spin3Index(lamX, lambda1, lambda2);
+ complex<double> amp = theMagExpi*current3SpinMap.at(IdLamXLam1Lam2);
+ if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData);
+ tmpResult+=amp;
+ }
+ }
+ if (_absDyn->isLdependent()){
+ tmpResult*=_cachedDynIdLSMap.at((*it)->L).at(_absDyn->grandMaId(grandmaAmp));
+ }
+ result+=tmpResult;
+ }
+
+ if (!_absDyn->isLdependent()) result *=_cachedDynIdMap.at(_absDyn->grandMaId(grandmaAmp));
+
+ result*=_isospinCG;
+ return result;
+}
+
+
+// void TensorDecAmps::getDefaultParams(fitParCol& fitVal, fitParCol& fitErr){
+
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentMagValMap;
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentPhiValMap;
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentMagErrMap;
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentPhiErrMap;
+
+// std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+// for(itLS=_LSs.begin(); itLS!=_LSs.end(); ++itLS){
+// currentMagValMap[*itLS]=_factorMag;
+// currentPhiValMap[*itLS]=0.;
+// currentMagErrMap[*itLS]=_factorMag/3.;
+// currentPhiErrMap[*itLS]=0.3;
+// }
+
+// fitVal.MagsLS[_key]=currentMagValMap;
+// fitVal.PhisLS[_key]=currentPhiValMap;
+// fitErr.MagsLS[_key]=currentMagErrMap;
+// fitErr.PhisLS[_key]=currentPhiErrMap;
+
+// _absDyn->getDefaultParams(fitVal, fitErr);
+
+
+// if(!_daughter1IsStable) _decAmpDaughter1->getDefaultParams(fitVal, fitErr);
+// if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);
+// }
+
+
+void TensorDecAmps::fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar){
+
+ std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+ for(itLS=_LSs.begin(); itLS!=_LSs.end(); ++itLS){
+ //fill magnitude
+ std::string magName=(*itLS)->name()+_key+"Mag";
+ double valMag=_factorMag;
+ double errMag=_factorMag/2.;
+ // double minMag=0.;
+ // double maxMag=_factorMag+30.*errMag;
+
+ fitPar->Add(magName, valMag, errMag);
+ // fitPar->SetLimits(magName, minMag, maxMag);
+
+ std::string phiName=(*itLS)->name()+_key+"Phi";
+ double valPhi=0.;
+ double errPhi=0.2;
+ //no limits for phi parameter
+ fitPar->Add(phiName, valPhi, errPhi);
+ }
+
+ _absDyn->fillDefaultParams(fitPar);
+
+ if(!_daughter1IsStable) _decAmpDaughter1->fillDefaultParams(fitPar);
+ if(!_daughter2IsStable) _decAmpDaughter2->fillDefaultParams(fitPar);
+}
+
+void TensorDecAmps::fillParamNameList(){
+ _paramNameList.clear();
+ std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+ for(itLS=_LSs.begin(); itLS!=_LSs.end(); ++itLS){
+ std::string magName=(*itLS)->name()+_key+"Mag";
+ _paramNameList.push_back(magName);
+
+ std::string phiName=(*itLS)->name()+_key+"Phi";
+ _paramNameList.push_back(phiName);
+ }
+}
+
+void TensorDecAmps::print(std::ostream& os) const{
+ return; //dummy
+}
+
+void TensorDecAmps::updateFitParams(std::shared_ptr<AbsPawianParameters> fitPar){
+
+ std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+ for(itLS=_LSs.begin(); itLS!=_LSs.end(); ++itLS){
+ //magnitude
+ std::string magName=(*itLS)->name()+_key+"Mag";
+ std::string phiName=(*itLS)->name()+_key+"Phi";
+ double theMag=fabs(fitPar->Value(magName));
+ double thePhi=fitPar->Value(phiName);
+
+ _currentParamMags[*itLS]=theMag;
+ _currentParamPhis[*itLS]=thePhi;
+
+ complex<double> expi(cos(thePhi), sin(thePhi));
+ _currentParamMagExpi[*itLS]=theMag*expi;
+ }
+
+ _absDyn->updateFitParams(fitPar);
+ if(!_daughter1IsStable) _decAmpDaughter1->updateFitParams(fitPar);
+ if(!_daughter2IsStable) _decAmpDaughter2->updateFitParams(fitPar);
+}
+
+
+void TensorDecAmps::calcDynamics(EvtData* theData, AbsXdecAmp* grandmaAmp){
+ if(!_recalculate) return;
+
+ if(!_absDyn->isLdependent()){
+ AbsXdecAmp::calcDynamics(theData, grandmaAmp);
+ return;
+ }
+
+ std::vector< std::shared_ptr<const LScomb> >::iterator it;
+ for (it=_LSs.begin(); it!=_LSs.end(); ++it){
+ _cachedDynIdLSMap[(*it)->L][_absDyn->grandMaId(grandmaAmp)]=_absDyn->eval(theData, grandmaAmp, (*it)->L);
+ }
+
+ if(!_daughter1IsStable) _decAmpDaughter1->calcDynamics(theData, this);
+ if(!_daughter2IsStable) _decAmpDaughter2->calcDynamics(theData, this);
+ return;
+}
+
+
diff --git a/PwaUtils/IsobarTensorDecay.cc b/PwaUtils/IsobarTensorDecay.cc
index 1a1692cf..9bdb772f 100644
--- a/PwaUtils/IsobarTensorDecay.cc
+++ b/PwaUtils/IsobarTensorDecay.cc
@@ -295,8 +295,8 @@ void IsobarTensorDecay::fillWignerDs(std::map<std::string, Vector4<double> >& fs
}
Id3StringType IdLamXLam1Lam2=FunctionUtils::spin3Index(lamMother, lamDaughter1, lamDaughter2);
unsigned short currentLSId=(*itJPCLS)->idnumberLS;
- // evtData->ComplexLS3Spin[_nameId][currentLSId][IdLamXLam1Lam2]=result(0);
- evtData->Complex3SpinLS[_nameId][IdLamXLam1Lam2][currentLSId]=result(0);
+ evtData->ComplexLS3Spin[_nameId][currentLSId][IdLamXLam1Lam2]=result(0);
+ //evtData->Complex3SpinLS[_nameId][IdLamXLam1Lam2][currentLSId]=result(0);
}
}
}
diff --git a/PwaUtils/TensorDecAmps.cc b/PwaUtils/TensorDecAmps.cc
index 2931e096..61836dfb 100644
--- a/PwaUtils/TensorDecAmps.cc
+++ b/PwaUtils/TensorDecAmps.cc
@@ -95,9 +95,9 @@ complex<double> TensorDecAmps::XdecPartAmp(const Spin& lamX, Spin& lamDec, short
complex<double> TensorDecAmps::XdecAmp(const Spin& lamX, EvtData* theData, AbsXdecAmp* grandmaAmp){
complex<double> result(0.,0.);
- if( std::abs(lamX) > _JPCPtr->J) return result;
+ if( fabs(lamX) > _JPCPtr->J) return result;
- short currentSpinIndex=FunctionUtils::spin1IdIndex(_projId,lamX);
+ short currentSpinIndex=FunctionUtils::spin1IdIndex(_projId,lamX);
if (!_recalculate){
return _cachedAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex);
@@ -109,6 +109,7 @@ complex<double> TensorDecAmps::XdecAmp(const Spin& lamX, EvtData* theData, AbsXd
_lam2MaxProj, true);
if ( _cacheAmps){
+ // _cachedAmpMap[evtNo][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
_cachedAmpIdMap[theData->evtNo][_absDyn->grandMaId(grandmaAmp)][currentSpinIndex]=result;
}
return result;
@@ -117,38 +118,75 @@ complex<double> TensorDecAmps::XdecAmp(const Spin& lamX, EvtData* theData, AbsXd
complex<double> TensorDecAmps::lsLoop(AbsXdecAmp* grandmaAmp, Spin lamX, EvtData* theData, Spin lam1Min,
Spin lam1Max, Spin lam2Min, Spin lam2Max, bool withDecs) {
+
complex<double> result(0.,0.);
- // map<unsigned short, map<Id3StringType, complex<double> > >& currentLS3SpinMap=theData->ComplexLS3Spin.at(_decay->nameId());
- map<Id3StringType, map<unsigned short, complex<double> > >& current3SpinLSMap=theData->Complex3SpinLS.at(_decay->nameId());
- std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+ map<unsigned short, map<Id3StringType, complex<double> > >& currentLS3SpinMap=theData->ComplexLS3Spin.at(_decay->nameId());
+ std::vector< std::shared_ptr<const LScomb> >::iterator it;
- for(Spin lambda1=-_Jdaughter1; lambda1<=_Jdaughter1; ++lambda1){
- if(lambda1<lam1Min || lambda1>lam1Max) continue;
- for(Spin lambda2=-_Jdaughter2; lambda2<=_Jdaughter2; ++lambda2){
- if(lambda2<lam2Min || lambda2>lam2Max) continue;
- Id3StringType IdLamXLam1Lam2=FunctionUtils::spin3Index(lamX, lambda1, lambda2);
- map<unsigned short, complex<double> >& currentLSMap = current3SpinLSMap.at(IdLamXLam1Lam2);
- complex<double> amp(0.,0.);
- for(itLS=_LSs.begin(); itLS!=_LSs.end(); ++itLS){
- //complex<double> tmpamp=_currentParamMagExpi.at(*itLS)*currentLS3SpinMap.at((*itLS)->idnumberLS).at(IdLamXLam1Lam2);
- complex<double> tmpamp=_currentParamMagExpi.at(*itLS)*currentLSMap.at((*itLS)->idnumberLS);
- if (_absDyn->isLdependent()){
- tmpamp *=_cachedDynIdLSMap.at((*itLS)->L).at(_absDyn->grandMaId(grandmaAmp));
- }
- amp+=tmpamp;
+ std::vector<std::vector<complex<double>>> dAmps(lam1Max - lam1Min + 1, std::vector<complex<double>>(lam2Max - lam2Min + 1));
+ if (withDecs){
+ for (Spin l1 = lam1Min; l1<=lam1Max; l1++){
+ for (Spin l2 = lam2Min; l2 <=lam2Max; l2++){
+ dAmps.at(l1 - lam1Min).at(l2 - lam2Min) = daughterAmp(l1, l2, theData);
}
- if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData);
- result+=amp;
}
}
+ for (it=_LSs.begin(); it!=_LSs.end(); ++it){
+ map<Id3StringType, complex<double> >& current3SpinMap = currentLS3SpinMap.at((*it)->idnumberLS);
+ complex<double> theMagExpi=_currentParamMagExpi.at(*it);
+
+ complex<double> tmpResult(0.,0.);
+ for(Spin lambda1=lam1Min; lambda1<=lam1Max; ++lambda1){
+ for(Spin lambda2=lam2Min; lambda2<=lam2Max; ++lambda2){
+ Id3StringType IdLamXLam1Lam2=FunctionUtils::spin3Index(lamX, lambda1, lambda2);
+ complex<double> amp = theMagExpi*current3SpinMap.at(IdLamXLam1Lam2);
+ if(withDecs) amp *=dAmps.at(lambda1 - lam1Min).at(lambda2 - lam2Min);
+ tmpResult+=amp;
+ }
+ }
+ if (_absDyn->isLdependent()){
+ tmpResult*=_cachedDynIdLSMap.at((*it)->L).at(_absDyn->grandMaId(grandmaAmp));
+ }
+ result+=tmpResult;
+ }
+
if (!_absDyn->isLdependent()) result *=_cachedDynIdMap.at(_absDyn->grandMaId(grandmaAmp));
+
result*=_isospinCG;
return result;
}
+// void TensorDecAmps::getDefaultParams(fitParCol& fitVal, fitParCol& fitErr){
+
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentMagValMap;
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentPhiValMap;
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentMagErrMap;
+// std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > currentPhiErrMap;
+
+// std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+// for(itLS=_LSs.begin(); itLS!=_LSs.end(); ++itLS){
+// currentMagValMap[*itLS]=_factorMag;
+// currentPhiValMap[*itLS]=0.;
+// currentMagErrMap[*itLS]=_factorMag/3.;
+// currentPhiErrMap[*itLS]=0.3;
+// }
+
+// fitVal.MagsLS[_key]=currentMagValMap;
+// fitVal.PhisLS[_key]=currentPhiValMap;
+// fitErr.MagsLS[_key]=currentMagErrMap;
+// fitErr.PhisLS[_key]=currentPhiErrMap;
+
+// _absDyn->getDefaultParams(fitVal, fitErr);
+
+
+// if(!_daughter1IsStable) _decAmpDaughter1->getDefaultParams(fitVal, fitErr);
+// if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);
+// }
+
+
void TensorDecAmps::fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar){
std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
@@ -199,10 +237,14 @@ void TensorDecAmps::updateFitParams(std::shared_ptr<AbsPawianParameters> fitPar)
//magnitude
std::string magName=(*itLS)->name()+_key+"Mag";
std::string phiName=(*itLS)->name()+_key+"Phi";
- double theMag=std::abs(fitPar->Value(magName));
+ double theMag=fabs(fitPar->Value(magName));
double thePhi=fitPar->Value(phiName);
- _currentParamMagExpi[*itLS]=std::polar(theMag, thePhi);
+ _currentParamMags[*itLS]=theMag;
+ _currentParamPhis[*itLS]=thePhi;
+
+ complex<double> expi(cos(thePhi), sin(thePhi));
+ _currentParamMagExpi[*itLS]=theMag*expi;
}
_absDyn->updateFitParams(fitPar);
@@ -230,4 +272,3 @@ void TensorDecAmps::calcDynamics(EvtData* theData, AbsXdecAmp* grandmaAmp){
}
-
diff --git a/PwaUtils/TensorDecAmps.hh b/PwaUtils/TensorDecAmps.hh
index e8408d58..f41d9abd 100644
--- a/PwaUtils/TensorDecAmps.hh
+++ b/PwaUtils/TensorDecAmps.hh
@@ -74,12 +74,18 @@ public:
protected:
std::vector< std::shared_ptr<const LScomb> > _LSs;
double _factorMag;
+ std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > _currentParamMags;
+ std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > _currentParamPhis;
std::map< std::shared_ptr<const LScomb>, complex<double>, pawian::Collection::SharedPtrLess > _currentParamMagExpi;
+ // std::map<std::thread::id, std::map<Spin, complex<double> > > _cachedDynLSMap;
std::map<unsigned short, std::map<Spin, complex<double> > > _cachedDynIdLSMap;
virtual complex<double> lsLoop(AbsXdecAmp* grandmaAmp, Spin lamX, EvtData* theData, Spin lam1Min, Spin lam1Max, Spin lam2Min, Spin lam2Max, bool withDecs);
-
private:
+
+
+
+
};
--
GitLab