From 6af99941f7d8ec29dae6726abcd7fa49ed0ab5e0 Mon Sep 17 00:00:00 2001
From: Bertram Kopf <bertram@pc14.ep1.rub.de>
Date: Thu, 30 Oct 2014 17:16:41 +0100
Subject: [PATCH] modified caching for LS and heli amplitudes for improving
 performance

---
 PwaUtils/AbsDecay.cc                          |  86 +++++---
 PwaUtils/AbsDecay.hh                          |   8 +-
 .../{HeliDecAmps.cc => AbsHeliDecAmps.cc}     | 130 +-----------
 .../{HeliDecAmps.hh => AbsHeliDecAmps.hh}     |  14 +-
 PwaUtils/{LSDecAmps.cc => AbsLSDecAmps.cc}    | 169 ++-------------
 PwaUtils/{LSDecAmps.hh => AbsLSDecAmps.hh}    |  22 +-
 PwaUtils/HeliDecNonRefAmps.cc                 | 156 ++++++++++++++
 PwaUtils/HeliDecNonRefAmps.hh                 |  72 +++++++
 PwaUtils/HeliDecRefAmps.cc                    | 163 +++++++++++++++
 PwaUtils/HeliDecRefAmps.hh                    |  72 +++++++
 PwaUtils/LSDecNonRefAmps.cc                   | 189 +++++++++++++++++
 PwaUtils/LSDecNonRefAmps.hh                   |  74 +++++++
 PwaUtils/LSDecRefAmps.cc                      | 193 ++++++++++++++++++
 PwaUtils/LSDecRefAmps.hh                      |  74 +++++++
 PwaUtils/LSOmegaTo3PiDecAmps.cc               |   2 +-
 PwaUtils/LSOmegaTo3PiDecAmps.hh               |   4 +-
 PwaUtils/XdecAmpRegistry.cc                   |  12 +-
 pbarpUtils/PbarpChannelEnv.cc                 |   8 +
 pbarpUtils/pbarpBaseLh.cc                     |   1 -
 pbarpUtils/pbarpCanoLh.cc                     |   1 -
 pbarpUtils/pbarpHeliLh.cc                     |   1 -
 pbarpUtils/pbarpTensorLh.cc                   |   1 -
 22 files changed, 1121 insertions(+), 331 deletions(-)
 rename PwaUtils/{HeliDecAmps.cc => AbsHeliDecAmps.cc} (58%)
 rename PwaUtils/{HeliDecAmps.hh => AbsHeliDecAmps.hh} (82%)
 rename PwaUtils/{LSDecAmps.cc => AbsLSDecAmps.cc} (51%)
 rename PwaUtils/{LSDecAmps.hh => AbsLSDecAmps.hh} (78%)
 create mode 100644 PwaUtils/HeliDecNonRefAmps.cc
 create mode 100644 PwaUtils/HeliDecNonRefAmps.hh
 create mode 100644 PwaUtils/HeliDecRefAmps.cc
 create mode 100644 PwaUtils/HeliDecRefAmps.hh
 create mode 100644 PwaUtils/LSDecNonRefAmps.cc
 create mode 100644 PwaUtils/LSDecNonRefAmps.hh
 create mode 100644 PwaUtils/LSDecRefAmps.cc
 create mode 100644 PwaUtils/LSDecRefAmps.hh

diff --git a/PwaUtils/AbsDecay.cc b/PwaUtils/AbsDecay.cc
index 7c5fa873..fd4f7279 100644
--- a/PwaUtils/AbsDecay.cc
+++ b/PwaUtils/AbsDecay.cc
@@ -73,6 +73,7 @@ AbsDecay::AbsDecay(Particle* mother, Particle* daughter1, Particle* daughter2, C
    ,_useProdBarrier(false)
    ,_massSumFsParticles(0.)
    ,_Lmin(0)
+   ,_decLevel(decayLevel::noLevel)
 {
   _absDecDaughter1=GlobalEnv::instance()->Channel(_channelId)->absDecayList()->decay(_daughter1);
   if(0 != _absDecDaughter1){
@@ -157,6 +158,7 @@ AbsDecay::AbsDecay(std::shared_ptr<const IGJPC> motherIGJPCPtr, Particle* daught
    ,_useProdBarrier(false)
    ,_massSumFsParticles(0.)
    ,_Lmin(0)
+   ,_decLevel(decayLevel::noLevel)
 {
   _absDecDaughter1=GlobalEnv::instance()->Channel(_channelId)->absDecayList()->decay(_daughter1);
 
@@ -327,46 +329,63 @@ void AbsDecay::fillWignerDs(std::map<std::string, Vector4<double> >& fsMap, Vect
     }
   }
 
-  for (Spin lamMother=-lamMotherMax; lamMother<=lamMotherMax; ++lamMother){
-    for (Spin lam12=-lam12Max; lam12<=lam12Max; ++lam12){
-      double thePhi=0.;
-      if(_hasMotherPart) thePhi=daughter2HelMother.Phi();
-      Id3StringType IdSpinMotherLamMotherLam12=FunctionUtils::spin3Index(spinMother, lamMother, lam12);
-
-      std::map<Id3StringType, complex<double> >::iterator found = evtData->WignerDStringStringId[_wignerDKey][refKey].find(IdSpinMotherLamMotherLam12);
-      if(found != evtData->WignerDStringStringId[_wignerDKey][refKey].end()){
-	continue;
-      } 
-      evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12]=Wigner_D(thePhi,daughter2HelMother.Theta(),0,spinMother,lamMother,lam12);
-
-      if(evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12] != evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12]){
-	DebugMsg << "WignerDStringString nan:\t" << evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12] << endmsg;
-    	DebugMsg << "daughter2HelMother.Theta():\t" << daughter2HelMother.Theta() << endmsg;
-	DebugMsg << "thePhi:\t" << thePhi << endmsg;
-	DebugMsg << "daughter2_4Vec:\t" << daughter2_4Vec << endmsg;
-	DebugMsg << "daughter2HelMother:\t" << daughter2HelMother << endmsg;
-	DebugMsg << "mother4Vec:\t" << mother4Vec << "\tP:" << mother4Vec.P() << endmsg;
+  if (whichDecayLevel()==decayLevel::isProdAmp || whichDecayLevel()==decayLevel::firstLevel){
+    for (Spin lamMother=-lamMotherMax; lamMother<=lamMotherMax; ++lamMother){
+      for (Spin lam12=-lam12Max; lam12<=lam12Max; ++lam12){
+	double thePhi=0.;
+	if(_hasMotherPart) thePhi=daughter2HelMother.Phi();
+	Id3StringType IdSpinMotherLamMotherLam12=FunctionUtils::spin3Index(spinMother, lamMother, lam12); 
+	evtData->WignerDStringId[_wignerDKey][IdSpinMotherLamMotherLam12]=Wigner_D(thePhi,daughter2HelMother.Theta(),0,spinMother,lamMother,lam12);
+	if(evtData->WignerDStringId[_wignerDKey][IdSpinMotherLamMotherLam12] != evtData->WignerDStringId[_wignerDKey][IdSpinMotherLamMotherLam12]){
+	  DebugMsg << "WignerDsString nan:\t" << evtData->WignerDsString[_wignerDKey][spinMother][lamMother][lam12] << endmsg;
+	  DebugMsg << "daughter2HelMother.Theta():\t" << daughter2HelMother.Theta() << endmsg;
+	  DebugMsg << "thePhi:\t" << thePhi << endmsg;
+	  DebugMsg << "daughter2_4Vec:\t" << daughter2_4Vec << endmsg;
+	  DebugMsg << "daughter2HelMother:\t" << daughter2HelMother << endmsg;
+	  DebugMsg << "mother4Vec:\t" << mother4Vec << "\tP:" << mother4Vec.P() << endmsg;
+	}
+      }
+    }
+  }
+  else{
+    for (Spin lamMother=-lamMotherMax; lamMother<=lamMotherMax; ++lamMother){
+      for (Spin lam12=-lam12Max; lam12<=lam12Max; ++lam12){
+	double thePhi=0.;
+	if(_hasMotherPart) thePhi=daughter2HelMother.Phi();
+	Id3StringType IdSpinMotherLamMotherLam12=FunctionUtils::spin3Index(spinMother, lamMother, lam12);
+	
+	std::map<Id3StringType, complex<double> >::iterator found = evtData->WignerDStringStringId[_wignerDKey][refKey].find(IdSpinMotherLamMotherLam12);
+	if(found != evtData->WignerDStringStringId[_wignerDKey][refKey].end()){
+	  continue;
+	} 
+	evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12]=Wigner_D(thePhi,daughter2HelMother.Theta(),0,spinMother,lamMother,lam12);
+	
+	if(evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12] != evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12]){
+	  DebugMsg << "WignerDStringString nan:\t" << evtData->WignerDStringStringId[_wignerDKey][refKey][IdSpinMotherLamMotherLam12] << endmsg;
+	  DebugMsg << "daughter2HelMother.Theta():\t" << daughter2HelMother.Theta() << endmsg;
+	  DebugMsg << "thePhi:\t" << thePhi << endmsg;
+	  DebugMsg << "daughter2_4Vec:\t" << daughter2_4Vec << endmsg;
+	  DebugMsg << "daughter2HelMother:\t" << daughter2HelMother << endmsg;
+	  DebugMsg << "mother4Vec:\t" << mother4Vec << "\tP:" << mother4Vec.P() << endmsg;
+	}
       }
     }
   }
 
+
   if(_isProdAmp){
      if (!_daughter1IsStable){
-       //       _absDecDaughter1->fillWignerDs(fsMap, daughter1_4Vec, evtData);
        _absDecDaughter1->fillWignerDs(fsMap, mother4Vec, evtData, _refKey); 
      }
      if (!_daughter2IsStable){
-       //       _absDecDaughter2->fillWignerDs(fsMap, daughter2_4Vec, evtData);
        _absDecDaughter2->fillWignerDs(fsMap, mother4Vec, evtData, _refKey);
      }
   }
    else{
      if (!_daughter1IsStable){
-       // _absDecDaughter1->fillWignerDs(fsMap, prodParticle4Vec, evtData);
        _absDecDaughter1->fillWignerDs(fsMap, mother4Vec, evtData, _refKey);
      }
      if (!_daughter2IsStable){
-       //       _absDecDaughter2->fillWignerDs(fsMap, prodParticle4Vec, evtData);
        _absDecDaughter2->fillWignerDs(fsMap, mother4Vec, evtData, _refKey);
      }
    }
@@ -475,3 +494,24 @@ void AbsDecay::extractLmin(){
   }
    
 }
+
+void AbsDecay::setDecayLevel(decLevel theLevel){
+  if (_decLevel==decayLevel::noLevel){
+    _decLevel=theLevel;
+  }
+  else if(_decLevel != theLevel) _decLevel=decayLevel::severalLevels;
+  else return;
+  Info << name() << " set decay level to " << _decLevel << endmsg; 
+}
+
+void AbsDecay::setDecayLevelTree(decLevel theLevel){
+
+  setDecayLevel(theLevel);
+
+  if (!_daughter1IsStable){
+    _absDecDaughter1->setDecayLevelTree(decayLevel(theLevel+1));
+  }
+  if (!_daughter2IsStable){
+    _absDecDaughter2->setDecayLevelTree(decayLevel(theLevel+1));
+  }
+}
diff --git a/PwaUtils/AbsDecay.hh b/PwaUtils/AbsDecay.hh
index 76927c85..aa736067 100644
--- a/PwaUtils/AbsDecay.hh
+++ b/PwaUtils/AbsDecay.hh
@@ -46,6 +46,8 @@ class AbsDynamics;
 class AbsDecay : public std::enable_shared_from_this<AbsDecay>{
 
 public:
+  typedef enum decayLevel {noLevel=-1, isProdAmp, firstLevel, secondLevel, thirdLevel, fourthLevel, fifthLevel, severalLevels} decLevel;
+  
   AbsDecay(Particle* mother, Particle* daughter1, Particle* daughter2, ChannelID channelId);
   AbsDecay(std::shared_ptr<const IGJPC> motherIGJPCPtr, Particle* daughter1, Particle* daughter2, std::string motherName, ChannelID channelId);
   virtual ~AbsDecay();
@@ -67,7 +69,10 @@ public:
   bool hasMother() {return _hasMotherPart;}
   bool isDaughter1Stable() {return _daughter1IsStable;}
   bool isDaughter2Stable() {return _daughter2IsStable;}
-  virtual bool isTensorAmp() {return false;} 
+  virtual bool isTensorAmp() {return false;}
+  virtual void setDecayLevel(decLevel theLevel);
+  virtual void setDecayLevelTree(decLevel theLevel); 
+  virtual decLevel whichDecayLevel() {return _decLevel;}
 
   std::vector<Particle*> finalStateParticles() {return _finalStateParticles;}
   std::vector<Particle*> finalStateParticlesDaughter2() {return _finalStateParticlesDaughter2;}
@@ -169,4 +174,5 @@ protected:
   double _massSumFsParticles;
 
   int _Lmin;
+  decLevel _decLevel;
 };
diff --git a/PwaUtils/HeliDecAmps.cc b/PwaUtils/AbsHeliDecAmps.cc
similarity index 58%
rename from PwaUtils/HeliDecAmps.cc
rename to PwaUtils/AbsHeliDecAmps.cc
index d7ac16fd..13a6df16 100644
--- a/PwaUtils/HeliDecAmps.cc
+++ b/PwaUtils/AbsHeliDecAmps.cc
@@ -1,6 +1,6 @@
 //************************************************************************//
 //									  //
-//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)			  //
+//  Copyright 2014 Bertram Kopf (bertram@ep1.rub.de)			  //
 //  	      	   Julian Pychy (julian@ep1.rub.de)			  //
 //          	   - Ruhr-Universität Bochum 				  //
 //									  //
@@ -21,15 +21,15 @@
 //									  //
 //************************************************************************//
 
-// HeliDecAmps class definition file. -*- C++ -*-
-// Copyright 2012 Bertram Kopf
+// AbsHeliDecAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
 
 #include <getopt.h>
 #include <fstream>
 #include <string>
 #include <mutex>
 
-#include "PwaUtils/HeliDecAmps.hh"
+#include "PwaUtils/AbsHeliDecAmps.hh"
 #include "qft++/relativistic-quantum-mechanics/Utils.hh"
 #include "PwaUtils/DataUtils.hh"
 #include "PwaUtils/AbsChannelEnv.hh"
@@ -39,7 +39,7 @@
 #include "Particle/Particle.hh"
 #include "ErrLogger/ErrLogger.hh"
 
-HeliDecAmps::HeliDecAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID) :
+AbsHeliDecAmps::AbsHeliDecAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID) :
   AbsXdecAmp(theDec, channelID)
   ,_JPClamlams(theDec->JPClamlamAmps())
   ,_factorMag(1.)
@@ -68,7 +68,7 @@ HeliDecAmps::HeliDecAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID chan
   }
 }
 
-HeliDecAmps::HeliDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+AbsHeliDecAmps::AbsHeliDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
   AbsXdecAmp(theDec, channelID)
 {
   Particle* daughter1=_decay->daughter1Part();
@@ -77,119 +77,11 @@ HeliDecAmps::HeliDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID)
   Info << "_parityFactor=\t" << _parityFactor << endmsg;
 }
 
-HeliDecAmps::~HeliDecAmps()
+AbsHeliDecAmps::~AbsHeliDecAmps()
 {
 }
 
-
-complex<double> HeliDecAmps::XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
-  complex<double> result(0.,0.);
-
-  std::string refKey=_refKey;
-  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
-
-
-  bool lamFs_daughter1=false;
-  if( _daughter1IsStable && _Jdaughter1>0) lamFs_daughter1=true;
-
-  bool lamFs_daughter2=false;
-  if( _daughter2IsStable && _Jdaughter2>0) lamFs_daughter2=true;
-
-  std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >::iterator it;
-
-  for(it=_currentParamMagLamLams.begin(); it!=_currentParamMagLamLams.end(); ++it){
-    std::shared_ptr<const JPClamlam> currentJPClamlam=it->first;
-    if( fabs(lamX) > currentJPClamlam->J) continue;
-
-    Spin lambda1= currentJPClamlam->lam1;
-    Spin lambda2= currentJPClamlam->lam2;
-    Spin lambda = lambda1-lambda2;
-    if( fabs(lambda) > currentJPClamlam->J) continue;
-    if(lamFs_daughter1 && lamFs!=lambda1) continue;
-    if(lamFs_daughter2 && lamFs!=lambda2) continue;
-    if(fixDaughterNr==1 && lamDec!=lambda1) continue;
-    if(fixDaughterNr==2 && lamDec!=lambda2) continue;
-
-    double theMag=it->second;
-    double thePhi=_currentParamPhiLamLams[currentJPClamlam];
-    complex<double> expi(cos(thePhi), sin(thePhi));
-    Id3StringType IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
-    complex<double> amp = currentJPClamlam->parityFactor*theMag*expi*conj(theData->WignerDStringStringId.at(_wignerDKey).at(refKey).at(IdJLamXLam12));
-    result+=amp;
-  }
-
-  result*=_preFactor*_isospinCG*sqrt(2.*_JPCPtr->J+1.);
-  return result;
-}
-
-
-
-
-complex<double> HeliDecAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
-
-  complex<double> result(0.,0.);
-
-  std::string refKey=_refKey;
-  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
-
-  if( fabs(lamX) > _JPCPtr->J) return result;
-
-  int evtNo=theData->evtNo;
-  Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs);
-
-  if ( _cacheAmps && !_recalculate){
-    result=_cachedAmpMapNew.at(evtNo).at(refKey).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
-    //    result*=_absDyn->eval(theData, grandmaAmp);
-    if(result.real()!=result.real()) DebugMsg << "result:\t" << result << endmsg;
-    return result;
-  }
-
-  std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >::iterator it;
-
-  for(it=_currentParamMagLamLams.begin(); it!=_currentParamMagLamLams.end(); ++it){
-
-    Spin lambda1= it->first->lam1;
-    Spin lambda2= it->first->lam2;
-    Spin lambda = lambda1-lambda2;
-    if( fabs(lambda) > it->first->J) continue;
-
-    if(_enabledlamFsDaughter1 && lamFs!=lambda1) continue;
-    if(_enabledlamFsDaughter2 && lamFs!=lambda2) continue;
-
-    double theMag=it->second;
-    double thePhi=_currentParamPhiLamLams.at(it->first);
-    complex<double> expi(cos(thePhi), sin(thePhi));
-    unsigned int IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
-
-    complex<double> amp = it->first->parityFactor*theMag*expi*conj( theData->WignerDStringStringId.at(_wignerDKey).at(refKey).at(IdJLamXLam12));
-    result+=amp*daughterAmp(lambda1, lambda2, theData, lamFs);
-  }
-
-  result*=_preFactor*_isospinCG*sqrt(2.*_JPCPtr->J+1.);
-
-  // if(absDec()->useProdBarrier()){
-  //   result *= BarrierFactor::BlattWeisskopf(absDec()->orbMomMin(), theData->DoubleString.at(_wignerDKey), BarrierFactor::qRDefault) /
-  //     BarrierFactor::BlattWeisskopf(absDec()->orbMomMin(), theData->DoubleString.at(_wignerDKey + "qNorm"), BarrierFactor::qRDefault);
-  // }
-  // else result*=_absDyn->eval(theData, grandmaAmp, absDec()->orbMomMin());
-
-  result*=_absDyn->eval(theData, grandmaAmp, absDec()->orbMomMin());
-  
-  if(result.real()!=result.real()){
-    Alert << "result:\t" << result << endmsg;
-    exit(0);
-  }
-
-  if ( _cacheAmps){
-    theMutex.lock();
-    _cachedAmpMapNew[evtNo][refKey][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
-    theMutex.unlock();
-  }
-
-  return result;
-}
-
-void  HeliDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
+void  AbsHeliDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
 
   std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess > currentMagValMap;
   std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess > currentPhiValMap;
@@ -215,12 +107,12 @@ void  HeliDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
   if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);
 }
 
-void HeliDecAmps::print(std::ostream& os) const{
+void AbsHeliDecAmps::print(std::ostream& os) const{
   return; //dummy
 }
 
 
-bool HeliDecAmps::checkRecalculation(fitParams& theParamVal){
+bool AbsHeliDecAmps::checkRecalculation(fitParams& theParamVal){
   _recalculate=false;
 
    if(_absDyn->checkRecalculation(theParamVal)) _recalculate=true;
@@ -254,7 +146,7 @@ bool HeliDecAmps::checkRecalculation(fitParams& theParamVal){
 }
 
 
-void  HeliDecAmps::updateFitParams(fitParams& theParamVal){
+void  AbsHeliDecAmps::updateFitParams(fitParams& theParamVal){
    std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >& magMap=theParamVal.MagLamLams[_key];
    std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >& phiMap=theParamVal.PhiLamLams[_key];
 
diff --git a/PwaUtils/HeliDecAmps.hh b/PwaUtils/AbsHeliDecAmps.hh
similarity index 82%
rename from PwaUtils/HeliDecAmps.hh
rename to PwaUtils/AbsHeliDecAmps.hh
index a4e88e6e..9140f88f 100644
--- a/PwaUtils/HeliDecAmps.hh
+++ b/PwaUtils/AbsHeliDecAmps.hh
@@ -21,7 +21,7 @@
 //									  //
 //************************************************************************//
 
-// HeliDecAmps class definition file. -*- C++ -*-
+// AbsHeliDecAmps class definition file. -*- C++ -*-
 // Copyright 2012 Bertram Kopf
 
 #pragma once
@@ -40,25 +40,21 @@
 class IsobarHeliDecay;
 class AbsDecay;
 
-class HeliDecAmps : public AbsXdecAmp{
+class AbsHeliDecAmps : public AbsXdecAmp{
 
 public:
 
   // create/copy/destroy:
 
   ///Constructor
-  HeliDecAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID);
-  HeliDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
+  AbsHeliDecAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID);
+  AbsHeliDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
   /** Destructor */
-  virtual ~HeliDecAmps();
+  virtual ~AbsHeliDecAmps();
 
 
   // Getters:
 
-  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
-  virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
-				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
-
   virtual void getDefaultParams(fitParams& fitVal, fitParams& fitErr);
   virtual void print(std::ostream& os) const;
   virtual bool checkRecalculation(fitParams& theParamVal);
diff --git a/PwaUtils/LSDecAmps.cc b/PwaUtils/AbsLSDecAmps.cc
similarity index 51%
rename from PwaUtils/LSDecAmps.cc
rename to PwaUtils/AbsLSDecAmps.cc
index 29eda954..01700afb 100644
--- a/PwaUtils/LSDecAmps.cc
+++ b/PwaUtils/AbsLSDecAmps.cc
@@ -1,6 +1,6 @@
 //************************************************************************//
 //									  //
-//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)			  //
+//  Copyright 2014 Bertram Kopf (bertram@ep1.rub.de)			  //
 //  	      	   Julian Pychy (julian@ep1.rub.de)			  //
 //          	   - Ruhr-Universität Bochum 				  //
 //									  //
@@ -21,15 +21,15 @@
 //									  //
 //************************************************************************//
 
-// LSDecAmps class definition file. -*- C++ -*-
-// Copyright 2012 Bertram Kopf
+// AbsLSDecAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
 
 #include <getopt.h>
 #include <fstream>
 #include <string>
 #include <mutex>
 
-#include "PwaUtils/LSDecAmps.hh"
+#include "PwaUtils/AbsLSDecAmps.hh"
 #include "qft++/relativistic-quantum-mechanics/Utils.hh"
 #include "PwaUtils/DataUtils.hh"
 #include "PwaUtils/GlobalEnv.hh"
@@ -39,7 +39,7 @@
 #include "Particle/Particle.hh"
 #include "ErrLogger/ErrLogger.hh"
 
-LSDecAmps::LSDecAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID) :
+AbsLSDecAmps::AbsLSDecAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID) :
   AbsXdecAmp(theDec, channelID)
   ,_LSs(theDec->LSAmps())
   ,_factorMag(1.)
@@ -52,7 +52,7 @@ LSDecAmps::LSDecAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID)
   fillCgPreFactor();
 }
 
-LSDecAmps::LSDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+AbsLSDecAmps::AbsLSDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
   AbsXdecAmp(theDec, channelID)
 {
   Particle* daughter1=_decay->daughter1Part();
@@ -62,156 +62,11 @@ LSDecAmps::LSDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
   fillCgPreFactor();
 }
 
-LSDecAmps::~LSDecAmps()
+AbsLSDecAmps::~AbsLSDecAmps()
 {
 }
 
-
-complex<double> LSDecAmps::XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
-
-  complex<double> result(0.,0.);
-
-  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;
-  }
-
-  result=lsLoop( grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, false);
-
-  return result;
-}
-
-
-
-
-complex<double> LSDecAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
-
-  std::string refKey=_refKey;
-  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
-  
-  // Info <<"\nlamX: " << lamX << "\tlamFs: " << lamFs << endmsg;
-  complex<double> result(0.,0.);
-  if( fabs(lamX) > _JPCPtr->J) return result;
-
-  int evtNo=theData->evtNo;
-
-  Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs); 
-  //  unsigned short currentSpinIndex=lamX.ToIndex()*100+lamFs.ToIndex();
-
-  if ( _cacheAmps && !_recalculate){
-    result=_cachedAmpMapNew.at(evtNo).at(refKey).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
-    return result;
-  }
-
-  //  Spin lam1Min=-_Jdaughter1;
-  Spin lam1Min=-_Jdaughter1;
-  Spin lam1Max= _Jdaughter1;
-  Spin lam2Min=-_Jdaughter2;
-  Spin lam2Max=_Jdaughter2;
-
-  if(_enabledlamFsDaughter1){
-    lam1Min=lamFs;
-    lam1Max=lamFs;
-  }
-  else if(_enabledlamFsDaughter2){
-    lam2Min=lamFs;
-    lam2Max=lamFs;
-  }
-
-
-  result=lsLoop(grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, true, lamFs);
-
-  if ( _cacheAmps){
-     theMutex.lock();
-     _cachedAmpMapNew[evtNo][refKey][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
-     theMutex.unlock();
-  }
-
-  if(result.real()!=result.real()){
-    Info << "dyn name: " << _absDyn->name() 
-	 << "\nname(): " << name()
-	 << endmsg;
-    Alert << "result:\t" << result << endmsg;
-    exit(0);
-  }
-  return result;
-}
-
-
-complex<double> LSDecAmps::lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs ){
-  std::string refKey=_refKey;
-  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
-  // Info << "\n_JPCPtr->J: " << _JPCPtr->J << "\tlamX: " << lamX << "\tlam1Min: " << lam1Min << "\tlam2Min: " << lam2Min << "\tlam1Max: " << lam1Max << "\tlam2Max: " << lam2Max << "\tlamFs: " <<lamFs << endmsg;
- 
-  complex<double> result(0.,0.);
-
-  //  map<Spin,complex<double> >& currentWignerDsMap=theData->WignerDsString.at(_wignerDKey).at(_JPCPtr->J).at(lamX);
-  //  Spin currentJ=_JPCPtr->J;
-  std::map<Id3StringType, complex<double> >& currentWignerDMap=theData->WignerDStringStringId.at(_wignerDKey).at(refKey);
-
-  std::vector< std::shared_ptr<const LScomb> >::iterator it;
-  for (it=_LSs.begin(); it!=_LSs.end(); ++it){
-
-    map<Spin,map<Spin, double > >& currentCgFactor=_cgPreFactor.at(*it);
-
-    double theMag=_currentParamMags.at(*it);
-    double thePhi=_currentParamPhis.at(*it);
-    complex<double> expi(cos(thePhi), sin(thePhi));
-
-    complex<double> tmpResult(0.,0.);
-    for(Spin lambda1=lam1Min; lambda1<=lam1Max; ++lambda1){
-      for(Spin lambda2=lam2Min; lambda2<=lam2Max; ++lambda2){
-	Spin lambda = lambda1-lambda2;
-	if( fabs(lambda)>_JPCPtr->J || fabs(lambda)>(*it)->S) continue;
-	Id3StringType IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
-	complex<double> amp = theMag*expi*currentCgFactor.at(lambda1).at(lambda2)*conj(currentWignerDMap.at(IdJLamXLam12));
-	//	complex<double> amp = theMag*expi*currentCgFactor.at(lambda1).at(lambda2)*conj(currentWignerDsMap.at(lambda));
-      	if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData, lamFs);
-	tmpResult+=amp;
-      }
-    }
-
-    // if(absDec()->useProdBarrier()){
-    //   //      tmpResult *= BarrierFactor::BlattWeisskopf((*it)->L, theData->DoubleString[_wignerDKey], 0.197);
-    //   tmpResult *= BarrierFactor::BlattWeisskopf((*it)->L, theData->DoubleString.at(_wignerDKey), BarrierFactor::qRDefault) /
-    //   	            BarrierFactor::BlattWeisskopf((*it)->L, theData->DoubleString.at(_wignerDKey + "qNorm"), BarrierFactor::qRDefault);
-    // }
-    
-    // else tmpResult*=_absDyn->eval(theData, grandmaAmp, (*it)->L);
-    tmpResult*=_absDyn->eval(theData, grandmaAmp, (*it)->L);
-
-    result+=tmpResult; 
-  }
-
-  result*=_preFactor*_isospinCG;
-  if(result.real()!=result.real()){
-    Alert << "result:\t" << result << endmsg;
-    exit(0);
-  }
-  return result;
-}
-
-
-void  LSDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
+void  AbsLSDecAmps::getDefaultParams(fitParams& fitVal, fitParams& 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;
@@ -238,12 +93,12 @@ void  LSDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
   if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);
 }
 
-void LSDecAmps::print(std::ostream& os) const{
+void AbsLSDecAmps::print(std::ostream& os) const{
   return; //dummy
 }
 
 
-bool LSDecAmps::checkRecalculation(fitParams& theParamVal){
+bool AbsLSDecAmps::checkRecalculation(fitParams& theParamVal){
   _recalculate=false;
 
    if(_absDyn->checkRecalculation(theParamVal)) _recalculate=true;
@@ -279,7 +134,7 @@ bool LSDecAmps::checkRecalculation(fitParams& theParamVal){
 }
 
 
-void  LSDecAmps::updateFitParams(fitParams& theParamVal){
+void  AbsLSDecAmps::updateFitParams(fitParams& theParamVal){
    std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess >& magMap=theParamVal.MagsLS[_key];
    std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess >& phiMap=theParamVal.PhisLS[_key];
 
@@ -298,7 +153,7 @@ void  LSDecAmps::updateFitParams(fitParams& theParamVal){
 
 }
 
-void  LSDecAmps::fillCgPreFactor(){
+void  AbsLSDecAmps::fillCgPreFactor(){
 
   std::vector< std::shared_ptr<const LScomb> >::iterator it;
   for (it=_LSs.begin(); it!=_LSs.end(); ++it){
diff --git a/PwaUtils/LSDecAmps.hh b/PwaUtils/AbsLSDecAmps.hh
similarity index 78%
rename from PwaUtils/LSDecAmps.hh
rename to PwaUtils/AbsLSDecAmps.hh
index 5f64fe65..56df8bd0 100644
--- a/PwaUtils/LSDecAmps.hh
+++ b/PwaUtils/AbsLSDecAmps.hh
@@ -1,6 +1,6 @@
 //************************************************************************//
 //									  //
-//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)			  //
+//  Copyright 2014 Bertram Kopf (bertram@ep1.rub.de)			  //
 //  	      	   Julian Pychy (julian@ep1.rub.de)			  //
 //          	   - Ruhr-Universität Bochum 				  //
 //									  //
@@ -21,8 +21,8 @@
 //									  //
 //************************************************************************//
 
-// LSDecAmps class definition file. -*- C++ -*-
-// Copyright 2012 Bertram Kopf
+// AbsLSDecAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
 
 #pragma once
 
@@ -40,24 +40,24 @@
 class IsobarLSDecay;
 class AbsDecay;
 
-class LSDecAmps : public AbsXdecAmp{
+class AbsLSDecAmps : public AbsXdecAmp{
 
 public:
 
   // create/copy/destroy:
 
   ///Constructor
-  LSDecAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID);
-  LSDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
+  AbsLSDecAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID);
+  AbsLSDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
   /** Destructor */
-  virtual ~LSDecAmps();
+  virtual ~AbsLSDecAmps();
 
 
   // Getters:
 
-  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
-  virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
-				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+  //  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+  // virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
+  //				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
 
   virtual void print(std::ostream& os) const;
   std::vector< std::shared_ptr<const LScomb> >& lsVec() {return _LSs;}
@@ -76,7 +76,7 @@ protected:
   std::map< std::shared_ptr<const LScomb>, double, pawian::Collection::SharedPtrLess > _currentParamPhis;
 
   void  fillCgPreFactor();
-  virtual complex<double> lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs=0 );
+  virtual complex<double> lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs=0 )=0;
 private:
 
 
diff --git a/PwaUtils/HeliDecNonRefAmps.cc b/PwaUtils/HeliDecNonRefAmps.cc
new file mode 100644
index 00000000..5e657c46
--- /dev/null
+++ b/PwaUtils/HeliDecNonRefAmps.cc
@@ -0,0 +1,156 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// HeliDecNonRefAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+#include <mutex>
+
+#include "PwaUtils/HeliDecNonRefAmps.hh"
+#include "qft++/relativistic-quantum-mechanics/Utils.hh"
+#include "PwaUtils/DataUtils.hh"
+#include "PwaUtils/AbsChannelEnv.hh"
+#include "PwaUtils/IsobarHeliDecay.hh"
+#include "PwaDynamics/BarrierFactor.hh"
+#include "Utils/FunctionUtils.hh"
+#include "Particle/Particle.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+HeliDecNonRefAmps::HeliDecNonRefAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID) :
+  AbsHeliDecAmps(theDec, channelID)
+{
+}
+
+HeliDecNonRefAmps::HeliDecNonRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+  AbsHeliDecAmps(theDec, channelID)
+{
+}
+
+HeliDecNonRefAmps::~HeliDecNonRefAmps()
+{
+}
+
+
+complex<double> HeliDecNonRefAmps::XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+  complex<double> result(0.,0.);
+
+  bool lamFs_daughter1=false;
+  if( _daughter1IsStable && _Jdaughter1>0) lamFs_daughter1=true;
+
+  bool lamFs_daughter2=false;
+  if( _daughter2IsStable && _Jdaughter2>0) lamFs_daughter2=true;
+
+  std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >::iterator it;
+
+  for(it=_currentParamMagLamLams.begin(); it!=_currentParamMagLamLams.end(); ++it){
+    std::shared_ptr<const JPClamlam> currentJPClamlam=it->first;
+    if( fabs(lamX) > currentJPClamlam->J) continue;
+
+    Spin lambda1= currentJPClamlam->lam1;
+    Spin lambda2= currentJPClamlam->lam2;
+    Spin lambda = lambda1-lambda2;
+    if( fabs(lambda) > currentJPClamlam->J) continue;
+    if(lamFs_daughter1 && lamFs!=lambda1) continue;
+    if(lamFs_daughter2 && lamFs!=lambda2) continue;
+    if(fixDaughterNr==1 && lamDec!=lambda1) continue;
+    if(fixDaughterNr==2 && lamDec!=lambda2) continue;
+
+    double theMag=it->second;
+    double thePhi=_currentParamPhiLamLams[currentJPClamlam];
+    complex<double> expi(cos(thePhi), sin(thePhi));
+    Id3StringType IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
+    complex<double> amp = currentJPClamlam->parityFactor*theMag*expi*conj(theData->WignerDStringId.at(_wignerDKey).at(IdJLamXLam12));
+    result+=amp;
+  }
+
+  result*=_preFactor*_isospinCG*sqrt(2.*_JPCPtr->J+1.);
+  return result;
+}
+
+
+
+
+complex<double> HeliDecNonRefAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+
+  complex<double> result(0.,0.);
+
+  if( fabs(lamX) > _JPCPtr->J) return result;
+
+  int evtNo=theData->evtNo;
+  Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs);
+
+  if ( _cacheAmps && !_recalculate){
+    result=_cachedAmpMap.at(evtNo).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
+    //    result*=_absDyn->eval(theData, grandmaAmp);
+    if(result.real()!=result.real()) DebugMsg << "result:\t" << result << endmsg;
+    return result;
+  }
+
+  std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >::iterator it;
+
+  for(it=_currentParamMagLamLams.begin(); it!=_currentParamMagLamLams.end(); ++it){
+
+    Spin lambda1= it->first->lam1;
+    Spin lambda2= it->first->lam2;
+    Spin lambda = lambda1-lambda2;
+    if( fabs(lambda) > it->first->J) continue;
+
+    if(_enabledlamFsDaughter1 && lamFs!=lambda1) continue;
+    if(_enabledlamFsDaughter2 && lamFs!=lambda2) continue;
+
+    double theMag=it->second;
+    double thePhi=_currentParamPhiLamLams.at(it->first);
+    complex<double> expi(cos(thePhi), sin(thePhi));
+    unsigned int IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
+
+    complex<double> amp = it->first->parityFactor*theMag*expi*conj( theData->WignerDStringId.at(_wignerDKey).at(IdJLamXLam12));
+    result+=amp*daughterAmp(lambda1, lambda2, theData, lamFs);
+  }
+
+  result*=_preFactor*_isospinCG*sqrt(2.*_JPCPtr->J+1.);
+
+  // if(absDec()->useProdBarrier()){
+  //   result *= BarrierFactor::BlattWeisskopf(absDec()->orbMomMin(), theData->DoubleString.at(_wignerDKey), BarrierFactor::qRDefault) /
+  //     BarrierFactor::BlattWeisskopf(absDec()->orbMomMin(), theData->DoubleString.at(_wignerDKey + "qNorm"), BarrierFactor::qRDefault);
+  // }
+  // else result*=_absDyn->eval(theData, grandmaAmp, absDec()->orbMomMin());
+
+  result*=_absDyn->eval(theData, grandmaAmp, absDec()->orbMomMin());
+  
+  if(result.real()!=result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
+
+  if ( _cacheAmps){
+    theMutex.lock();
+    _cachedAmpMap[evtNo][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
+    theMutex.unlock();
+  }
+
+  return result;
+}
+
diff --git a/PwaUtils/HeliDecNonRefAmps.hh b/PwaUtils/HeliDecNonRefAmps.hh
new file mode 100644
index 00000000..5c82b4b7
--- /dev/null
+++ b/PwaUtils/HeliDecNonRefAmps.hh
@@ -0,0 +1,72 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// HeliDecNonRefAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+
+#include <cassert>
+#include <memory>
+
+#include "PwaUtils/AbsHeliDecAmps.hh"
+
+class IsobarHeliDecay;
+class AbsDecay;
+
+class HeliDecNonRefAmps : public AbsHeliDecAmps{
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor
+  HeliDecNonRefAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID);
+  HeliDecNonRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
+  /** Destructor */
+  virtual ~HeliDecNonRefAmps();
+
+
+  // Getters:
+
+  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+  virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
+				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+
+protected:
+
+private:
+
+
+
+
+};
+
+
+
diff --git a/PwaUtils/HeliDecRefAmps.cc b/PwaUtils/HeliDecRefAmps.cc
new file mode 100644
index 00000000..5ceb2c59
--- /dev/null
+++ b/PwaUtils/HeliDecRefAmps.cc
@@ -0,0 +1,163 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// HeliDecRefAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+#include <mutex>
+
+#include "PwaUtils/HeliDecRefAmps.hh"
+#include "qft++/relativistic-quantum-mechanics/Utils.hh"
+#include "PwaUtils/DataUtils.hh"
+#include "PwaUtils/AbsChannelEnv.hh"
+#include "PwaUtils/IsobarHeliDecay.hh"
+#include "PwaDynamics/BarrierFactor.hh"
+#include "Utils/FunctionUtils.hh"
+#include "Particle/Particle.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+HeliDecRefAmps::HeliDecRefAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID) :
+  AbsHeliDecAmps(theDec, channelID)
+{
+}
+
+HeliDecRefAmps::HeliDecRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+  AbsHeliDecAmps(theDec, channelID)
+{
+}
+
+HeliDecRefAmps::~HeliDecRefAmps()
+{
+}
+
+
+complex<double> HeliDecRefAmps::XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+  complex<double> result(0.,0.);
+
+  std::string refKey=_refKey;
+  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
+
+
+  bool lamFs_daughter1=false;
+  if( _daughter1IsStable && _Jdaughter1>0) lamFs_daughter1=true;
+
+  bool lamFs_daughter2=false;
+  if( _daughter2IsStable && _Jdaughter2>0) lamFs_daughter2=true;
+
+  std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >::iterator it;
+
+  for(it=_currentParamMagLamLams.begin(); it!=_currentParamMagLamLams.end(); ++it){
+    std::shared_ptr<const JPClamlam> currentJPClamlam=it->first;
+    if( fabs(lamX) > currentJPClamlam->J) continue;
+
+    Spin lambda1= currentJPClamlam->lam1;
+    Spin lambda2= currentJPClamlam->lam2;
+    Spin lambda = lambda1-lambda2;
+    if( fabs(lambda) > currentJPClamlam->J) continue;
+    if(lamFs_daughter1 && lamFs!=lambda1) continue;
+    if(lamFs_daughter2 && lamFs!=lambda2) continue;
+    if(fixDaughterNr==1 && lamDec!=lambda1) continue;
+    if(fixDaughterNr==2 && lamDec!=lambda2) continue;
+
+    double theMag=it->second;
+    double thePhi=_currentParamPhiLamLams[currentJPClamlam];
+    complex<double> expi(cos(thePhi), sin(thePhi));
+    Id3StringType IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
+    complex<double> amp = currentJPClamlam->parityFactor*theMag*expi*conj(theData->WignerDStringStringId.at(_wignerDKey).at(refKey).at(IdJLamXLam12));
+    result+=amp;
+  }
+
+  result*=_preFactor*_isospinCG*sqrt(2.*_JPCPtr->J+1.);
+  return result;
+}
+
+
+
+
+complex<double> HeliDecRefAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+
+  complex<double> result(0.,0.);
+
+  std::string refKey=_refKey;
+  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
+
+  if( fabs(lamX) > _JPCPtr->J) return result;
+
+  int evtNo=theData->evtNo;
+  Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs);
+
+  if ( _cacheAmps && !_recalculate){
+    result=_cachedAmpMapNew.at(evtNo).at(refKey).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
+    //    result*=_absDyn->eval(theData, grandmaAmp);
+    if(result.real()!=result.real()) DebugMsg << "result:\t" << result << endmsg;
+    return result;
+  }
+
+  std::map< std::shared_ptr<const JPClamlam>, double, pawian::Collection::SharedPtrLess >::iterator it;
+
+  for(it=_currentParamMagLamLams.begin(); it!=_currentParamMagLamLams.end(); ++it){
+
+    Spin lambda1= it->first->lam1;
+    Spin lambda2= it->first->lam2;
+    Spin lambda = lambda1-lambda2;
+    if( fabs(lambda) > it->first->J) continue;
+
+    if(_enabledlamFsDaughter1 && lamFs!=lambda1) continue;
+    if(_enabledlamFsDaughter2 && lamFs!=lambda2) continue;
+
+    double theMag=it->second;
+    double thePhi=_currentParamPhiLamLams.at(it->first);
+    complex<double> expi(cos(thePhi), sin(thePhi));
+    unsigned int IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
+
+    complex<double> amp = it->first->parityFactor*theMag*expi*conj( theData->WignerDStringStringId.at(_wignerDKey).at(refKey).at(IdJLamXLam12));
+    result+=amp*daughterAmp(lambda1, lambda2, theData, lamFs);
+  }
+
+  result*=_preFactor*_isospinCG*sqrt(2.*_JPCPtr->J+1.);
+
+  // if(absDec()->useProdBarrier()){
+  //   result *= BarrierFactor::BlattWeisskopf(absDec()->orbMomMin(), theData->DoubleString.at(_wignerDKey), BarrierFactor::qRDefault) /
+  //     BarrierFactor::BlattWeisskopf(absDec()->orbMomMin(), theData->DoubleString.at(_wignerDKey + "qNorm"), BarrierFactor::qRDefault);
+  // }
+  // else result*=_absDyn->eval(theData, grandmaAmp, absDec()->orbMomMin());
+
+  result*=_absDyn->eval(theData, grandmaAmp, absDec()->orbMomMin());
+  
+  if(result.real()!=result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
+
+  if ( _cacheAmps){
+    theMutex.lock();
+    _cachedAmpMapNew[evtNo][refKey][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
+    theMutex.unlock();
+  }
+
+  return result;
+}
+
diff --git a/PwaUtils/HeliDecRefAmps.hh b/PwaUtils/HeliDecRefAmps.hh
new file mode 100644
index 00000000..71406afd
--- /dev/null
+++ b/PwaUtils/HeliDecRefAmps.hh
@@ -0,0 +1,72 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// HeliDecRefAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+
+#include <cassert>
+#include <memory>
+
+#include "PwaUtils/AbsHeliDecAmps.hh"
+
+class IsobarHeliDecay;
+class AbsDecay;
+
+class HeliDecRefAmps : public AbsHeliDecAmps{
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor
+  HeliDecRefAmps(std::shared_ptr<IsobarHeliDecay> theDec, ChannelID channelID);
+  HeliDecRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
+  /** Destructor */
+  virtual ~HeliDecRefAmps();
+
+
+  // Getters:
+
+  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+  virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
+				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+
+protected:
+
+private:
+
+
+
+
+};
+
+
+
diff --git a/PwaUtils/LSDecNonRefAmps.cc b/PwaUtils/LSDecNonRefAmps.cc
new file mode 100644
index 00000000..801f629a
--- /dev/null
+++ b/PwaUtils/LSDecNonRefAmps.cc
@@ -0,0 +1,189 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// LSDecNonRefAmps class definition file. -*- C++ -*-
+// Copyright 2012 Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+#include <mutex>
+
+#include "PwaUtils/LSDecNonRefAmps.hh"
+#include "qft++/relativistic-quantum-mechanics/Utils.hh"
+#include "PwaUtils/DataUtils.hh"
+#include "PwaUtils/GlobalEnv.hh"
+#include "PwaUtils/IsobarLSDecay.hh"
+#include "PwaDynamics/BarrierFactor.hh"
+#include "Utils/FunctionUtils.hh"
+#include "Particle/Particle.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+LSDecNonRefAmps::LSDecNonRefAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID) :
+  AbsLSDecAmps(theDec, channelID)
+{
+}
+
+LSDecNonRefAmps::LSDecNonRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+  AbsLSDecAmps(theDec, channelID)
+{
+}
+
+LSDecNonRefAmps::~LSDecNonRefAmps()
+{
+}
+
+
+complex<double> LSDecNonRefAmps::XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+
+  complex<double> result(0.,0.);
+
+  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;
+  }
+
+  result=lsLoop( grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, false);
+
+  return result;
+}
+
+
+
+
+complex<double> LSDecNonRefAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+
+  // Info <<"\nlamX: " << lamX << "\tlamFs: " << lamFs << endmsg;
+  complex<double> result(0.,0.);
+  if( fabs(lamX) > _JPCPtr->J) return result;
+
+  int evtNo=theData->evtNo;
+
+  Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs); 
+  //  unsigned short currentSpinIndex=lamX.ToIndex()*100+lamFs.ToIndex();
+
+  if ( _cacheAmps && !_recalculate){
+    result=_cachedAmpMap.at(evtNo).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
+    return result;
+  }
+
+  //  Spin lam1Min=-_Jdaughter1;
+  Spin lam1Min=-_Jdaughter1;
+  Spin lam1Max= _Jdaughter1;
+  Spin lam2Min=-_Jdaughter2;
+  Spin lam2Max=_Jdaughter2;
+
+  if(_enabledlamFsDaughter1){
+    lam1Min=lamFs;
+    lam1Max=lamFs;
+  }
+  else if(_enabledlamFsDaughter2){
+    lam2Min=lamFs;
+    lam2Max=lamFs;
+  }
+
+
+  result=lsLoop(grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, true, lamFs);
+
+  if ( _cacheAmps){
+     theMutex.lock();
+     _cachedAmpMap[evtNo][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
+     theMutex.unlock();
+  }
+
+  if(result.real()!=result.real()){
+    Info << "dyn name: " << _absDyn->name() 
+	 << "\nname(): " << name()
+	 << endmsg;
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
+  return result;
+}
+
+
+complex<double> LSDecNonRefAmps::lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs ){
+ 
+  complex<double> result(0.,0.);
+
+  //  map<Spin,complex<double> >& currentWignerDsMap=theData->WignerDsString.at(_wignerDKey).at(_JPCPtr->J).at(lamX);
+  //  Spin currentJ=_JPCPtr->J;
+  //  std::map<Id3StringType, complex<double> >& currentWignerDMap=theData->WignerDStringStringId.at(_wignerDKey).at(refKey);
+
+  std::map<Id3StringType, complex<double> >& currentWignerDMap=theData->WignerDStringId.at(_wignerDKey);
+
+  std::vector< std::shared_ptr<const LScomb> >::iterator it;
+  for (it=_LSs.begin(); it!=_LSs.end(); ++it){
+
+    map<Spin,map<Spin, double > >& currentCgFactor=_cgPreFactor.at(*it);
+
+    double theMag=_currentParamMags.at(*it);
+    double thePhi=_currentParamPhis.at(*it);
+    complex<double> expi(cos(thePhi), sin(thePhi));
+
+    complex<double> tmpResult(0.,0.);
+    for(Spin lambda1=lam1Min; lambda1<=lam1Max; ++lambda1){
+      for(Spin lambda2=lam2Min; lambda2<=lam2Max; ++lambda2){
+	Spin lambda = lambda1-lambda2;
+	if( fabs(lambda)>_JPCPtr->J || fabs(lambda)>(*it)->S) continue;
+	Id3StringType IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
+	complex<double> amp = theMag*expi*currentCgFactor.at(lambda1).at(lambda2)*conj(currentWignerDMap.at(IdJLamXLam12));
+	//	complex<double> amp = theMag*expi*currentCgFactor.at(lambda1).at(lambda2)*conj(currentWignerDsMap.at(lambda));
+      	if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData, lamFs);
+	tmpResult+=amp;
+      }
+    }
+
+    tmpResult*=_absDyn->eval(theData, grandmaAmp, (*it)->L);
+
+    result+=tmpResult; 
+  }
+
+  result*=_preFactor*_isospinCG;
+  if(result.real()!=result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
+  return result;
+}
+
+
diff --git a/PwaUtils/LSDecNonRefAmps.hh b/PwaUtils/LSDecNonRefAmps.hh
new file mode 100644
index 00000000..a00d314f
--- /dev/null
+++ b/PwaUtils/LSDecNonRefAmps.hh
@@ -0,0 +1,74 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// LSDecNonRefAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+
+#include <cassert>
+#include <memory>
+
+#include "PwaUtils/AbsLSDecAmps.hh"
+
+class IsobarLSDecay;
+class AbsDecay;
+
+class LSDecNonRefAmps : public AbsLSDecAmps{
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor
+  LSDecNonRefAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID);
+  LSDecNonRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
+  /** Destructor */
+  virtual ~LSDecNonRefAmps();
+
+
+  // Getters:
+
+  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+  virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
+				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+
+
+protected:
+  virtual complex<double> lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs=0 );
+
+private:
+
+
+
+
+};
+
+
+
diff --git a/PwaUtils/LSDecRefAmps.cc b/PwaUtils/LSDecRefAmps.cc
new file mode 100644
index 00000000..41e1f818
--- /dev/null
+++ b/PwaUtils/LSDecRefAmps.cc
@@ -0,0 +1,193 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// LSDecRefAmps class definition file. -*- C++ -*-
+// Copyright 2012 Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+#include <mutex>
+
+#include "PwaUtils/LSDecRefAmps.hh"
+#include "qft++/relativistic-quantum-mechanics/Utils.hh"
+#include "PwaUtils/DataUtils.hh"
+#include "PwaUtils/GlobalEnv.hh"
+#include "PwaUtils/IsobarLSDecay.hh"
+#include "PwaDynamics/BarrierFactor.hh"
+#include "Utils/FunctionUtils.hh"
+#include "Particle/Particle.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+LSDecRefAmps::LSDecRefAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID) :
+  AbsLSDecAmps(theDec, channelID)
+{
+}
+
+LSDecRefAmps::LSDecRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
+  AbsLSDecAmps(theDec, channelID)
+{
+}
+
+LSDecRefAmps::~LSDecRefAmps()
+{
+}
+
+
+complex<double> LSDecRefAmps::XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+
+  complex<double> result(0.,0.);
+
+  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;
+  }
+
+  result=lsLoop( grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, false);
+
+  return result;
+}
+
+
+
+
+complex<double> LSDecRefAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
+
+  std::string refKey=_refKey;
+  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
+  
+  // Info <<"\nlamX: " << lamX << "\tlamFs: " << lamFs << endmsg;
+  complex<double> result(0.,0.);
+  if( fabs(lamX) > _JPCPtr->J) return result;
+
+  int evtNo=theData->evtNo;
+
+  Id2StringType currentSpinIndex=FunctionUtils::spin2Index(lamX,lamFs); 
+  //  unsigned short currentSpinIndex=lamX.ToIndex()*100+lamFs.ToIndex();
+
+  if ( _cacheAmps && !_recalculate){
+    result=_cachedAmpMapNew.at(evtNo).at(refKey).at(_absDyn->grandMaKey(grandmaAmp)).at(currentSpinIndex);
+    return result;
+  }
+
+  //  Spin lam1Min=-_Jdaughter1;
+  Spin lam1Min=-_Jdaughter1;
+  Spin lam1Max= _Jdaughter1;
+  Spin lam2Min=-_Jdaughter2;
+  Spin lam2Max=_Jdaughter2;
+
+  if(_enabledlamFsDaughter1){
+    lam1Min=lamFs;
+    lam1Max=lamFs;
+  }
+  else if(_enabledlamFsDaughter2){
+    lam2Min=lamFs;
+    lam2Max=lamFs;
+  }
+
+
+  result=lsLoop(grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, true, lamFs);
+
+  if ( _cacheAmps){
+     theMutex.lock();
+     _cachedAmpMapNew[evtNo][refKey][_absDyn->grandMaKey(grandmaAmp)][currentSpinIndex]=result;
+     theMutex.unlock();
+  }
+
+  if(result.real()!=result.real()){
+    Info << "dyn name: " << _absDyn->name() 
+	 << "\nname(): " << name()
+	 << endmsg;
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
+  return result;
+}
+
+
+complex<double> LSDecRefAmps::lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs ){
+  std::string refKey=_refKey;
+  if (0!=grandmaAmp) refKey=grandmaAmp->refKey();
+  // Info << "\n_JPCPtr->J: " << _JPCPtr->J << "\tlamX: " << lamX << "\tlam1Min: " << lam1Min << "\tlam2Min: " << lam2Min << "\tlam1Max: " << lam1Max << "\tlam2Max: " << lam2Max << "\tlamFs: " <<lamFs << endmsg;
+ 
+  complex<double> result(0.,0.);
+
+  //  map<Spin,complex<double> >& currentWignerDsMap=theData->WignerDsString.at(_wignerDKey).at(_JPCPtr->J).at(lamX);
+  //  Spin currentJ=_JPCPtr->J;
+  std::map<Id3StringType, complex<double> >& currentWignerDMap=theData->WignerDStringStringId.at(_wignerDKey).at(refKey);
+
+  std::vector< std::shared_ptr<const LScomb> >::iterator it;
+  for (it=_LSs.begin(); it!=_LSs.end(); ++it){
+
+    map<Spin,map<Spin, double > >& currentCgFactor=_cgPreFactor.at(*it);
+
+    double theMag=_currentParamMags.at(*it);
+    double thePhi=_currentParamPhis.at(*it);
+    complex<double> expi(cos(thePhi), sin(thePhi));
+
+    complex<double> tmpResult(0.,0.);
+    for(Spin lambda1=lam1Min; lambda1<=lam1Max; ++lambda1){
+      for(Spin lambda2=lam2Min; lambda2<=lam2Max; ++lambda2){
+	Spin lambda = lambda1-lambda2;
+	if( fabs(lambda)>_JPCPtr->J || fabs(lambda)>(*it)->S) continue;
+	Id3StringType IdJLamXLam12=FunctionUtils::spin3Index(_J, lamX, lambda);
+	complex<double> amp = theMag*expi*currentCgFactor.at(lambda1).at(lambda2)*conj(currentWignerDMap.at(IdJLamXLam12));
+	//	complex<double> amp = theMag*expi*currentCgFactor.at(lambda1).at(lambda2)*conj(currentWignerDsMap.at(lambda));
+      	if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData, lamFs);
+	tmpResult+=amp;
+      }
+    }
+
+    tmpResult*=_absDyn->eval(theData, grandmaAmp, (*it)->L);
+
+    result+=tmpResult; 
+  }
+
+  result*=_preFactor*_isospinCG;
+  if(result.real()!=result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
+  return result;
+}
+
+
diff --git a/PwaUtils/LSDecRefAmps.hh b/PwaUtils/LSDecRefAmps.hh
new file mode 100644
index 00000000..a846835c
--- /dev/null
+++ b/PwaUtils/LSDecRefAmps.hh
@@ -0,0 +1,74 @@
+//************************************************************************//
+//									  //
+//  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/>.	  //
+//									  //
+//************************************************************************//
+
+// LSDecRefAmps class definition file. -*- C++ -*-
+// Copyright 2014 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+
+#include <cassert>
+#include <memory>
+
+#include "PwaUtils/AbsLSDecAmps.hh"
+
+class IsobarLSDecay;
+class AbsDecay;
+
+class LSDecRefAmps : public AbsLSDecAmps{
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor
+  LSDecRefAmps(std::shared_ptr<IsobarLSDecay> theDec, ChannelID channelID);
+  LSDecRefAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID);
+  /** Destructor */
+  virtual ~LSDecRefAmps();
+
+
+  // Getters:
+
+  virtual complex<double> XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+  virtual complex<double> XdecPartAmp(Spin& lamX, Spin& lamDec, short fixDaughterNr,
+				      EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp);
+
+
+protected:
+  virtual complex<double> lsLoop(AbsXdecAmp* grandmaAmp, Spin& lamX, EvtData* theData, Spin& lam1Min, Spin& lam1Max, Spin& lam2Min, Spin& lam2Max, bool withDecs, Spin lamFs=0 );
+
+private:
+
+
+
+
+};
+
+
+
diff --git a/PwaUtils/LSOmegaTo3PiDecAmps.cc b/PwaUtils/LSOmegaTo3PiDecAmps.cc
index 405accbb..ab423271 100644
--- a/PwaUtils/LSOmegaTo3PiDecAmps.cc
+++ b/PwaUtils/LSOmegaTo3PiDecAmps.cc
@@ -39,7 +39,7 @@
 
 
 LSOmegaTo3PiDecAmps::LSOmegaTo3PiDecAmps(std::shared_ptr<OmegaTo3PiLSDecay> theDec, ChannelID channelID) :
-  LSDecAmps( (std::shared_ptr<AbsDecay>) theDec, channelID)
+  LSDecRefAmps( (std::shared_ptr<AbsDecay>) theDec, channelID)
   ,_lambdaDecKey(theDec->lambdaDecKey())
 {
   _LSs=theDec->LSAmps();
diff --git a/PwaUtils/LSOmegaTo3PiDecAmps.hh b/PwaUtils/LSOmegaTo3PiDecAmps.hh
index 8f7b379e..68ad3039 100644
--- a/PwaUtils/LSOmegaTo3PiDecAmps.hh
+++ b/PwaUtils/LSOmegaTo3PiDecAmps.hh
@@ -35,11 +35,11 @@
 #include <cassert>
 #include <memory>
 
-#include "PwaUtils/LSDecAmps.hh"
+#include "PwaUtils/LSDecRefAmps.hh"
 #include "PwaUtils/OmegaTo3PiLSDecay.hh"
 //class OmegaTo3PiLSDecay;
 
-class LSOmegaTo3PiDecAmps : public LSDecAmps{
+class LSOmegaTo3PiDecAmps : public LSDecRefAmps{
 
 public:
 
diff --git a/PwaUtils/XdecAmpRegistry.cc b/PwaUtils/XdecAmpRegistry.cc
index 5e41ab14..31b8c4b7 100644
--- a/PwaUtils/XdecAmpRegistry.cc
+++ b/PwaUtils/XdecAmpRegistry.cc
@@ -37,8 +37,10 @@
 #include "PwaUtils/OmegaTo3PiLSDecay.hh"
 #include "PwaUtils/OmegaTo3PiTensorDecay.hh"
 #include "PwaUtils/AbsXdecAmp.hh"
-#include "PwaUtils/LSDecAmps.hh"
-#include "PwaUtils/HeliDecAmps.hh"
+#include "PwaUtils/LSDecNonRefAmps.hh"
+#include "PwaUtils/LSDecRefAmps.hh"
+#include "PwaUtils/HeliDecNonRefAmps.hh"
+#include "PwaUtils/HeliDecRefAmps.hh"
 #include "PwaUtils/TensorDecAmps.hh"
 #include "PwaUtils/TensorPsiToGamXDecAmps.hh"
 #include "PwaUtils/LSOmegaTo3PiDecAmps.hh"
@@ -72,11 +74,13 @@ std::shared_ptr<AbsXdecAmp> XdecAmpRegistry::getXdecAmp(short channelID, std::sh
   else{
     if(theAbsXDec->type()=="IsobarLSDecay"){
       std::shared_ptr<IsobarLSDecay> decLS =  std::dynamic_pointer_cast<IsobarLSDecay>(theAbsXDec);
-      result=std::shared_ptr<AbsXdecAmp>(new LSDecAmps(decLS, channelID));
+      if(theAbsXDec->whichDecayLevel()==AbsDecay::decayLevel::isProdAmp || theAbsXDec->whichDecayLevel()==AbsDecay::decayLevel::firstLevel) result=std::shared_ptr<AbsXdecAmp>(new LSDecNonRefAmps(decLS, channelID)); 
+      else result=std::shared_ptr<AbsXdecAmp>(new LSDecRefAmps(decLS, channelID));
     }
     else if(theAbsXDec->type()=="IsobarHeliDecay"){
       std::shared_ptr<IsobarHeliDecay> decLamLam =  std::dynamic_pointer_cast<IsobarHeliDecay>(theAbsXDec);
-      result=std::shared_ptr<AbsXdecAmp>(new HeliDecAmps(decLamLam, channelID));
+      if(theAbsXDec->whichDecayLevel()==AbsDecay::decayLevel::isProdAmp || theAbsXDec->whichDecayLevel()==AbsDecay::decayLevel::firstLevel) result=std::shared_ptr<AbsXdecAmp>(new HeliDecNonRefAmps(decLamLam, channelID));
+      else result=std::shared_ptr<AbsXdecAmp>(new HeliDecRefAmps(decLamLam, channelID));
     }
     else if(theAbsXDec->type()=="IsobarTensorDecay"){
       std::shared_ptr<IsobarTensorDecay> decTensor =  std::dynamic_pointer_cast<IsobarTensorDecay>(theAbsXDec);
diff --git a/pbarpUtils/PbarpChannelEnv.cc b/pbarpUtils/PbarpChannelEnv.cc
index b687bf78..7f386e9b 100644
--- a/pbarpUtils/PbarpChannelEnv.cc
+++ b/pbarpUtils/PbarpChannelEnv.cc
@@ -223,6 +223,14 @@ void PbarpChannelEnv::setup(ChannelID id){
       }
    }
 
+
+   //set decay levels
+   std::vector<std::shared_ptr<AbsDecay> > prodDecList= _prodDecList->getList();
+   std::vector<std::shared_ptr<AbsDecay> >::iterator itProdDecList;
+   for (itProdDecList=prodDecList.begin(); itProdDecList!=prodDecList.end(); ++itProdDecList){
+     (*itProdDecList)->setDecayLevelTree(AbsDecay::decayLevel::isProdAmp);    
+   } 
+
    // spin density particles
    _spinDensity = _theParser->spinDensityNames();
 
diff --git a/pbarpUtils/pbarpBaseLh.cc b/pbarpUtils/pbarpBaseLh.cc
index fbe2d683..077ec52d 100644
--- a/pbarpUtils/pbarpBaseLh.cc
+++ b/pbarpUtils/pbarpBaseLh.cc
@@ -32,7 +32,6 @@
 #include "pbarpUtils/pbarpReaction.hh"
 #include "pbarpUtils/PbarpChannelEnv.hh"
 #include "PwaUtils/GlobalEnv.hh"
-#include "PwaUtils/LSDecAmps.hh"
 #include "PwaUtils/EvtDataBaseList.hh"
 #include "PwaUtils/AbsXdecAmp.hh"
 #include "PwaUtils/AbsDecay.hh"
diff --git a/pbarpUtils/pbarpCanoLh.cc b/pbarpUtils/pbarpCanoLh.cc
index 96828e68..fc3c4e28 100644
--- a/pbarpUtils/pbarpCanoLh.cc
+++ b/pbarpUtils/pbarpCanoLh.cc
@@ -30,7 +30,6 @@
 
 #include "pbarpUtils/pbarpCanoLh.hh"
 #include "pbarpUtils/pbarpReaction.hh"
-#include "PwaUtils/LSDecAmps.hh"
 #include "PwaUtils/EvtDataBaseList.hh"
 #include "PwaUtils/AbsXdecAmp.hh"
 #include "PwaUtils/AbsDecay.hh"
diff --git a/pbarpUtils/pbarpHeliLh.cc b/pbarpUtils/pbarpHeliLh.cc
index bafd7b16..a833a5a0 100644
--- a/pbarpUtils/pbarpHeliLh.cc
+++ b/pbarpUtils/pbarpHeliLh.cc
@@ -31,7 +31,6 @@
 #include "pbarpUtils/pbarpHeliLh.hh"
 #include "pbarpUtils/pbarpReaction.hh"
 #include "PwaUtils/HeliDecAmps.hh"
-#include "PwaUtils/LSDecAmps.hh"
 #include "PwaUtils/EvtDataBaseList.hh"
 #include "PwaUtils/AbsXdecAmp.hh"
 #include "PwaUtils/AbsDecay.hh"
diff --git a/pbarpUtils/pbarpTensorLh.cc b/pbarpUtils/pbarpTensorLh.cc
index bf32a47c..a7acb68e 100644
--- a/pbarpUtils/pbarpTensorLh.cc
+++ b/pbarpUtils/pbarpTensorLh.cc
@@ -32,7 +32,6 @@
 #include "pbarpUtils/pbarpReaction.hh"
 #include "pbarpUtils/PbarpChannelEnv.hh"
 #include "PwaUtils/GlobalEnv.hh"
-#include "PwaUtils/LSDecAmps.hh"
 #include "PwaUtils/EvtDataBaseList.hh"
 #include "PwaUtils/AbsXdecAmp.hh"
 #include "PwaUtils/AbsDecay.hh"
-- 
GitLab