From b9e19512ce7a3a8b7dd57d7e41146c07664757c4 Mon Sep 17 00:00:00 2001
From: Bertram Kopf <bertram@ep1.rub.de>
Date: Fri, 12 Jan 2024 14:44:36 +0100
Subject: [PATCH] improved performance for tensor amp calculations

---
 PwaUtils/TensorDecAmps.cc           | 74 ++++++++---------------------
 PwaUtils/TensorDecAmps.hh           |  8 +---
 PwaUtils/TensorOmegaTo3PiDecAmps.cc | 19 ++++----
 PwaUtils/TensorOmegaTo3PiDecAmps.hh |  3 +-
 PwaUtils/TensorPsiToGamXDecAmps.cc  | 19 ++------
 PwaUtils/TensorPsiToGamXDecAmps.hh  |  4 +-
 6 files changed, 39 insertions(+), 88 deletions(-)

diff --git a/PwaUtils/TensorDecAmps.cc b/PwaUtils/TensorDecAmps.cc
index 346fa658..0a31f355 100644
--- a/PwaUtils/TensorDecAmps.cc
+++ b/PwaUtils/TensorDecAmps.cc
@@ -95,7 +95,7 @@ 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( fabs(lamX) > _JPCPtr->J) return result;
+  if( std::abs(lamX) > _JPCPtr->J) return result;
 
   short currentSpinIndex=FunctionUtils::spin1IdIndex(_projId,lamX);  
   
@@ -109,7 +109,6 @@ 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;
@@ -121,61 +120,31 @@ complex<double> TensorDecAmps::lsLoop(AbsXdecAmp* grandmaAmp, Spin lamX, EvtData
   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;
+  std::vector< std::shared_ptr<const LScomb> >::const_iterator itLS;
+
+  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);
+      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);
+	if (_absDyn->isLdependent()){
+          tmpamp *=_cachedDynIdLSMap.at((*itLS)->L).at(_absDyn->grandMaId(grandmaAmp));
+        }
+        amp+=tmpamp;
       }
+      if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData);
+      result+=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;
@@ -226,14 +195,10 @@ 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=fabs(fitPar->Value(magName));
+    double theMag=std::abs(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;
+    _currentParamMagExpi[*itLS]=std::polar(theMag, thePhi);    
   }
 
   _absDyn->updateFitParams(fitPar);
@@ -261,3 +226,4 @@ void TensorDecAmps::calcDynamics(EvtData* theData, AbsXdecAmp* grandmaAmp){
 }
 
 
+
diff --git a/PwaUtils/TensorDecAmps.hh b/PwaUtils/TensorDecAmps.hh
index f41d9abd..e8408d58 100644
--- a/PwaUtils/TensorDecAmps.hh
+++ b/PwaUtils/TensorDecAmps.hh
@@ -74,18 +74,12 @@ 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:
-
-
-
 
+private:
 };
 
 
diff --git a/PwaUtils/TensorOmegaTo3PiDecAmps.cc b/PwaUtils/TensorOmegaTo3PiDecAmps.cc
index 69c086d4..728677e5 100644
--- a/PwaUtils/TensorOmegaTo3PiDecAmps.cc
+++ b/PwaUtils/TensorOmegaTo3PiDecAmps.cc
@@ -31,7 +31,7 @@
 
 #include "PwaUtils/TensorOmegaTo3PiDecAmps.hh"
 #include "qft++/relativistic-quantum-mechanics/Utils.hh"
-#include "PwaUtils/DataUtils.hh"
+//#include "PwaUtils/DataUtils.hh"
 #include "PwaUtils/OmegaTo3PiTensorDecay.hh"
 #include "Utils/FunctionUtils.hh"
 #include "Particle/Particle.hh"
@@ -67,31 +67,34 @@ complex<double> TensorOmegaTo3PiDecAmps::XdecPartAmp(const Spin& lamX, Spin& lam
 
 complex<double> TensorOmegaTo3PiDecAmps::XdecAmp(const Spin& lamX, EvtData* theData, AbsXdecAmp* grandmaAmp){
 
-  complex<double> result(0.,0.);
+  //  complex<double> result(0.,0.);
 
   //  int evtNo=theData->evtNo;
   // Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs);
   Id1StringType currentSpinIndex=FunctionUtils::spin1Index(lamX);
   
   if (!_recalculate){
-    result=_cachedLocalAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex);
-    return result;
+    //  result=_cachedLocalAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex);
+    // return result;
+    return _cachedLocalAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex); 
   }
 
+  complex<double> result(0.,0.);
   Spin motherJ(_JPCPtr->J);
   std::vector< std::shared_ptr<const LScomb> >::iterator it;
   for (it=_LSs.begin(); it!=_LSs.end(); ++it){
     Id2StringType IdJLamX=FunctionUtils::spin2Index(motherJ,lamX);
-    complex<double> amp = _currentMagExpi[*it]*theData->Complex2Spin.at(_decay->wignerDId()).at(IdJLamX);
-    result+=amp;
+    //complex<double> amp = _currentMagExpi[*it]*theData->Complex2Spin.at(_decay->wignerDId()).at(IdJLamX);
+    //result+=amp;
+    result+=_currentMagExpi[*it]*theData->Complex2Spin.at(_decay->wignerDId()).at(IdJLamX); 
   }
 
-  result*=100.;
+  result*=100.*_absDyn->eval(theData, grandmaAmp);
   if ( _cacheAmps){
      _cachedLocalAmpIdMap[theData->evtNo][_absDyn->grandMaId(grandmaAmp)][currentSpinIndex]=result;
   }
 
-  result*=_absDyn->eval(theData, grandmaAmp);
+  // result*=_absDyn->eval(theData, grandmaAmp);
   //  InfoMsg <<"TensorOmegaTo3PiDecAmps result: " << result << endmsg; 
   return result;
 }
diff --git a/PwaUtils/TensorOmegaTo3PiDecAmps.hh b/PwaUtils/TensorOmegaTo3PiDecAmps.hh
index 3874d2aa..4cd0e2c3 100644
--- a/PwaUtils/TensorOmegaTo3PiDecAmps.hh
+++ b/PwaUtils/TensorOmegaTo3PiDecAmps.hh
@@ -36,7 +36,7 @@
 #include <memory>
 
 #include "PwaUtils/AbsXdecAmp.hh"
-
+#include "PwaUtils/DataUtils.hh"
 class OmegaTo3PiTensorDecay;
 class Particle;
 class AbsPawianParameters;
@@ -80,6 +80,7 @@ protected:
   Particle* _daughter2;
   Particle* _daughter3;
   tensorOmegaTo3PiCachedIdMap _cachedLocalAmpIdMap;
+  
 private:
 
 
diff --git a/PwaUtils/TensorPsiToGamXDecAmps.cc b/PwaUtils/TensorPsiToGamXDecAmps.cc
index 212a25d0..be9b24b7 100644
--- a/PwaUtils/TensorPsiToGamXDecAmps.cc
+++ b/PwaUtils/TensorPsiToGamXDecAmps.cc
@@ -119,21 +119,13 @@ complex<double> TensorPsiToGamXDecAmps::XdecAmp(const Spin& lamX, EvtData* theDa
   short currentSpinIndex=FunctionUtils::spin1IdIndex(_projId,lamX);
   
   if (!_recalculate){
-    //    result=_cachedAmpMap.at(evtNo).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
     return _cachedAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex);
-    //    result*=_absDyn->eval(theData, grandmaAmp);
-    //    return result;
   }
 
   for(Spin lambda2=_lam2MinProj; lambda2<=_lam2MaxProj; ++lambda2){
     Id3StringType IdLamMotherLamGamLamX=FunctionUtils::spin3Index(lamX, _lam1MinProj, lambda2);
     complex<double> tmpResult(0.,0.);
     for(int i=0; i<_noOfAmps; ++i){
-      // double theMag=_currentParamLocalMags.at(i);
-      // double thePhi=_currentParamLocalPhis.at(i);
-      // complex<double> expi(cos(thePhi), sin(thePhi));
-      // tmpResult+=theMag*expi*theData->ComplexDoubleInt3SpinString.at(_name).at(i).at(lamX).at(lamFs).at(lambda2)*_absDyn->eval(theData, grandmaAmp,_ampLMap.at(i));
-      //      tmpResult+=_currentParamLocalMagExpi.at(i)*theData->ComplexDoubleInt3SpinString.at(_name).at(i).at(lamX).at(lamFs).at(lambda2)*_absDyn->eval(theData, grandmaAmp,_ampLMap.at(i));
       tmpResult+=_currentParamLocalMagExpi.at(i)*theData->ComplexN3Spin.at(_decay->nameId()).at(i).at(IdLamMotherLamGamLamX)*_absDyn->eval(theData, grandmaAmp,_ampLMap.at(i));
     }
 
@@ -141,7 +133,6 @@ complex<double> TensorPsiToGamXDecAmps::XdecAmp(const Spin& lamX, EvtData* theDa
   }
 
   if ( _cacheAmps){
-     //_cachedAmpMap[evtNo][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
      _cachedAmpIdMap[theData->evtNo][_absDyn->grandMaId(grandmaAmp)][currentSpinIndex]=result;
   }
 
@@ -192,15 +183,11 @@ void TensorPsiToGamXDecAmps::updateFitParams(std::shared_ptr<AbsPawianParameters
 
   for (int i=0; i<_noOfAmps; ++i){
     std::string magName=_MagParamNames.at(i);
-    double theMag=fabs(fitPar->Value(magName));
-    _currentParamLocalMags[i]=theMag;
-
+    double theMag=std::abs(fitPar->Value(magName));
+    
     std::string phiName=_PhiParamNames.at(i);
     double thePhi=fitPar->Value(phiName);
-    _currentParamLocalPhis[i]=thePhi;
-
-    complex<double> expi(cos(thePhi), sin(thePhi));
-    _currentParamLocalMagExpi[i]=theMag*expi;
+    _currentParamLocalMagExpi[i]=std::polar(theMag, thePhi);
   }
 
   _absDyn->updateFitParams(fitPar);
diff --git a/PwaUtils/TensorPsiToGamXDecAmps.hh b/PwaUtils/TensorPsiToGamXDecAmps.hh
index cbae7935..19b642e3 100644
--- a/PwaUtils/TensorPsiToGamXDecAmps.hh
+++ b/PwaUtils/TensorPsiToGamXDecAmps.hh
@@ -67,8 +67,8 @@ public:
   virtual void fillParamNameList();
 
 protected:
-  std::map< int, double> _currentParamMags;
-  std::map< int, double> _currentParamPhis;
+  // std::map< int, double> _currentParamMags;
+  // std::map< int, double> _currentParamPhis;
   short _noOfAmps;
   std::vector<std::string> _MagParamNames;
   std::vector<std::string> _PhiParamNames;
-- 
GitLab