added abstract class for dynamics

PwaDynamics/Jamfile

@@ -5,7 +5,7 @@ lib PwaDynamics :
 	[ glob *.cc : *App.cc ] 
 	$(TOP)/Utils//Utils 
 	$(TOP)/qft++//qft++
-	$(TOP)/ErrLogger//ErrLogger 
+	$(TOP)/ErrLogger//ErrLogger
 	: 
 	: 
 	:  ;

PwaUtils/AbsDynamics.cc

+// AbsDynamics class definition file. -*- C++ -*-
+// Copyright 20123Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+
+#include "PwaUtils/AbsDynamics.hh"
+#include "ErrLogger/ErrLogger.hh"
+#include "Particle/Particle.hh"
+
+AbsDynamics::AbsDynamics(std::vector<Particle*>& fsParticles, Particle* mother) :
+  AbsParamHandler()
+  ,_fsParticles(fsParticles)
+  ,_mother(mother)
+{
+}
+
+AbsDynamics::~AbsDynamics()
+{
+}
+
+void AbsDynamics::cacheAmplitudes(){
+  _cacheAmps=true;
+}
+

PwaUtils/AbsDynamics.hh

+// AbsDynamics class definition file. -*- C++ -*-
+// Copyright 2013 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+#include <boost/shared_ptr.hpp>
+
+#include "PwaUtils/EvtDataBaseList.hh"
+#include "PwaUtils/FitParamsBase.hh"
+#include "PwaUtils/AbsParamHandler.hh"
+
+class Particle;
+
+class AbsDynamics : public AbsParamHandler{
+
+public:
+  AbsDynamics(std::vector<Particle*>& fsParticles, Particle* mother );
+  virtual ~AbsDynamics();
+
+  virtual complex<double> eval(EvtData* theData, Spin OrbMom=0)=0;
+  virtual void cacheAmplitudes();
+
+protected:
+  std::vector<Particle*> _fsParticles;
+  Particle* _mother;
+  std::map<int, complex<float> >  _cachedMap;
+
+private:
+
+};

PwaUtils/AbsXdecAmp.cc

@@ -11,6 +11,9 @@
 #include "qft++/relativistic-quantum-mechanics/Utils.hh"
 #include "ErrLogger/ErrLogger.hh"
 #include "Particle/Particle.hh"
+//#include "PwaUtils/AbsDynamics.hh"
+#include "PwaUtils/BreitWignerDynamics.hh"
+#include "PwaUtils/WoDynamics.hh"
 
 AbsXdecAmp::AbsXdecAmp(boost::shared_ptr<AbsDecay> theDec) :
   AbsParamHandler()
@@ -32,12 +35,17 @@ AbsXdecAmp::~AbsXdecAmp()
 }
 
 void AbsXdecAmp::initialize(){
+  std::vector<Particle*> fsParticles=_decay->finalStateParticles();
   if(_withDyn){
     if(!_decay->hasMother()){
       Alert << "no mother resonance; can not add dynamis" << endmsg;
       exit(1);
     }
     _massKey=_decay->massParKey();
+      _absDyn= boost::shared_ptr<AbsDynamics>(new BreitWignerDynamics(_massKey, fsParticles, _decay->motherPart()));
+  }
+  else{
+    _absDyn= boost::shared_ptr<AbsDynamics>(new WoDynamics(fsParticles, _decay->motherPart()));
   }
   
   if(!_daughter1IsStable){
@@ -63,6 +71,7 @@ complex<double> AbsXdecAmp::daughterAmp(Spin lam1, Spin lam2, EvtData* theData,
 
 void AbsXdecAmp::cacheAmplitudes(){
   _cacheAmps=true;
+  _absDyn->cacheAmplitudes();
   if(!_daughter1IsStable) _decAmpDaughter1->cacheAmplitudes();
   if(!_daughter2IsStable) _decAmpDaughter2->cacheAmplitudes();
 }

PwaUtils/AbsXdecAmp.hh

@@ -13,8 +13,10 @@
 #include "PwaUtils/EvtDataBaseList.hh"
 #include "PwaUtils/FitParamsBase.hh"
 #include "PwaUtils/AbsParamHandler.hh"
+#include "PwaUtils/AbsDynamics.hh"
 
 class AbsDecay;
+//class AbsDynamics;
 
 class AbsXdecAmp : public AbsParamHandler{
 
@@ -40,6 +42,7 @@ protected:
   boost::shared_ptr<AbsDecay> _decay; 
   const std::string _name;
   boost::shared_ptr<const jpcRes> _JPCPtr;
+  boost::shared_ptr<AbsDynamics> _absDyn;
   const std::vector<std::string> _hypVec;
   Spin _J_X;
   int _parity;
@@ -58,7 +61,7 @@ protected:
   double _currentXWidth;
   bool _daughter1IsStable;
   bool _daughter2IsStable;
-  std::map<int, std::map<Spin, std::map<Spin, complex<double> > > > _cachedAmpMap;
+  std::map<int, std::map<Spin, std::map<Spin, complex<float> > > > _cachedAmpMap;
 
   virtual void initialize();
 };

PwaUtils/BreitWignerDynamics.cc

+// BreitWignerDynamics class definition file. -*- C++ -*-
+// Copyright 20123Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+
+#include "PwaUtils/BreitWignerDynamics.hh"
+#include "ErrLogger/ErrLogger.hh"
+#include "Particle/Particle.hh"
+
+BreitWignerDynamics::BreitWignerDynamics(std::string& massKey, std::vector<Particle*>& fsParticles, Particle* mother) :
+  AbsDynamics(fsParticles, mother)
+  ,_massKey(massKey)
+{
+}
+
+BreitWignerDynamics::~BreitWignerDynamics()
+{
+}
+
+complex<double> BreitWignerDynamics::eval(EvtData* theData, Spin OrbMom){
+  int evtNo=theData->evtNo;
+  if ( _cacheAmps && !_recalculate){
+    return _cachedMap[evtNo];
+  }
+
+  Vector4<double> mass4Vec(0.,0.,0.,0.);
+  std::vector<Particle*>::iterator it;
+  for (it=_fsParticles.begin(); it != _fsParticles.end(); ++it){
+       mass4Vec+=theData->FourVecsString[(*it)->name()];
+    }
+
+  complex<double> result=BreitWigner(mass4Vec, _currentMass, _currentWidth);
+  
+  if ( _cacheAmps){
+#ifdef _OPENMP
+#pragma omp critical
+    {
+#endif
+      _cachedMap[evtNo]=result;
+#ifdef _OPENMP
+    }
+#endif
+  }  
+
+  return result;
+}
+
+void  BreitWignerDynamics::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
+    fitVal.Masses[_massKey]=_mother->mass();
+    fitErr.Masses[_massKey]=0.03;
+    fitVal.Widths[_massKey]=_mother->width();
+    fitErr.Widths[_massKey]=0.2*_mother->width();
+}
+
+bool BreitWignerDynamics::checkRecalculation(fitParams& theParamVal){
+  _recalculate=false;
+
+  double mass=theParamVal.Masses[_massKey];
+  if ( fabs(mass-_currentMass) > 1.e-10){
+    _recalculate=true;
+  }
+  double width=theParamVal.Widths[_massKey];
+  if ( fabs(width-_currentWidth) > 1.e-10){
+    _recalculate=true;
+  }
+
+  return _recalculate;
+}
+
+void BreitWignerDynamics::updateFitParams(fitParams& theParamVal){
+  _currentMass=theParamVal.Masses[_massKey];
+  _currentWidth=theParamVal.Widths[_massKey];
+}

PwaUtils/BreitWignerDynamics.hh

+// BreitWignerDynamics class definition file. -*- C++ -*-
+// Copyright 2013 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+#include <boost/shared_ptr.hpp>
+
+#include "PwaUtils/AbsDynamics.hh"
+
+class BreitWignerDynamics : public AbsDynamics{
+
+public:
+  BreitWignerDynamics(std::string& massKey, std::vector<Particle*>& fsParticles, Particle* mother);
+  virtual ~BreitWignerDynamics();
+
+  virtual complex<double> eval(EvtData* theData, Spin OrbMom=0);
+  
+  virtual void getDefaultParams(fitParams& fitVal, fitParams& fitErr);
+  virtual bool checkRecalculation(fitParams& theParamVal);
+  virtual void updateFitParams(fitParams& theParamVal);
+
+protected:
+  std::string _massKey;
+
+  double _currentMass;
+  double _currentWidth;
+  std::map<int, complex<double> >  _cachedMap; 
+
+private:
+
+};

PwaUtils/HeliDecAmps.cc

@@ -140,16 +140,18 @@ int evtNo=theData->evtNo;
     result+=amp*daughterAmp(lambda1, lambda2, theData, lamFs);
   }
 
-  if(_withDyn){
-    Vector4<double> mass4Vec(0.,0.,0.,0.);
-    std::vector<Particle*> fsParticleVec=_decay->finalStateParticles();
+  result*=_absDyn->eval(theData);
 
-    std::vector<Particle*>::iterator itPartVec;
-    for (itPartVec=fsParticleVec.begin(); itPartVec!=fsParticleVec.end(); ++itPartVec){
-       mass4Vec+=theData->FourVecsString[(*itPartVec)->name()];
-    }
-    result*=BreitWigner(mass4Vec, _currentXMass, _currentXWidth);
-  }
+  // if(_withDyn){
+  //   Vector4<double> mass4Vec(0.,0.,0.,0.);
+  //   std::vector<Particle*> fsParticleVec=_decay->finalStateParticles();
+
+  //   std::vector<Particle*>::iterator itPartVec;
+  //   for (itPartVec=fsParticleVec.begin(); itPartVec!=fsParticleVec.end(); ++itPartVec){
+  //      mass4Vec+=theData->FourVecsString[(*itPartVec)->name()];
+  //   }
+  //   result*=BreitWigner(mass4Vec, _currentXMass, _currentXWidth);
+  // }
 
   //  result*=sqrt((2.*_JPCPtr->J+1.)/12.56637);
   result*=sqrt(2.*_JPCPtr->J+1.);
@@ -187,12 +189,7 @@ void  HeliDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
   fitErr.MagLamLams[_key]=currentMagErrMap;
   fitErr.PhiLamLams[_key]=currentPhiErrMap;
 
-  if(_withDyn){
-    fitVal.Masses[_massKey]=_decay->motherPart()->mass();
-    fitErr.Masses[_massKey]=0.03;
-    fitVal.Widths[_massKey]=_decay->motherPart()->width();
-    fitErr.Widths[_massKey]=0.2*_decay->motherPart()->width();
-  }
+  _absDyn->getDefaultParams(fitVal, fitErr);
 
   if(!_daughter1IsStable) _decAmpDaughter1->getDefaultParams(fitVal, fitErr);
   if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);  
@@ -222,16 +219,7 @@ bool HeliDecAmps::checkRecalculation(fitParams& theParamVal){
      }
    }
 
-   if(_withDyn){
-     double mass=theParamVal.Masses[_massKey];
-     if ( fabs(mass-_currentXMass) > 1.e-10){
-       _recalculate=true;
-     }
-     double width=theParamVal.Widths[_massKey];
-     if ( fabs(width-_currentXWidth) > 1.e-10){
-       _recalculate=true;
-     }
-   }
+   if(_absDyn->checkRecalculation(theParamVal)) _recalculate=true; 
 
    if(!_daughter1IsStable) {
      if(_decAmpDaughter1->checkRecalculation(theParamVal)) _recalculate=true;
@@ -262,14 +250,10 @@ void  HeliDecAmps::updateFitParams(fitParams& theParamVal){
      }
    }
 
-   if (_withDyn){
-     double xMass=theParamVal.Masses[_massKey];
-     _currentXMass= xMass;
-     double xWidth=theParamVal.Widths[_massKey];
-     _currentXWidth=xWidth;
-   }
 
-  if(!_daughter1IsStable) _decAmpDaughter1->updateFitParams(theParamVal);
-  if(!_daughter2IsStable) _decAmpDaughter2->updateFitParams(theParamVal);
+   _absDyn->updateFitParams(theParamVal);
 
+   if(!_daughter1IsStable) _decAmpDaughter1->updateFitParams(theParamVal);
+   if(!_daughter2IsStable) _decAmpDaughter2->updateFitParams(theParamVal);
+   
 }

PwaUtils/Jamfile

@@ -5,7 +5,7 @@ lib PwaUtils :
 	[ glob *.cc : *App.cc ] 
 	$(TOP)/Utils//Utils 
 	$(TOP)/qft++//qft++
-	$(TOP)/ErrLogger//ErrLogger 
+	$(TOP)/ErrLogger//ErrLogger
 	: <use>$(TOP)//Minuit2
 	: 
 	: <library>$(TOP)//Minuit2 ;

PwaUtils/LSDecAmps.cc

@@ -99,10 +99,11 @@ complex<double> LSDecAmps::XdecPartAmp(Spin lamX, Spin lamDec, short fixDaughter
 
 
 complex<double> LSDecAmps::XdecAmp(Spin lamX, EvtData* theData, Spin lamFs){
-int evtNo=theData->evtNo;
-  
-  if ( _cacheAmps && !_recalculate){
-    complex<double> result(0.,0.);
+
+  int evtNo=theData->evtNo;
+  complex<double> result(0.,0.);  
+ 
+ if ( _cacheAmps && !_recalculate){
     result= _cachedAmpMap[evtNo][lamX][lamFs];
     return result;
   }
@@ -122,7 +123,6 @@ int evtNo=theData->evtNo;
     lam2Max=lamFs;
   }
 
-  complex<double> result(0.,0.);
   std::vector< boost::shared_ptr<const JPCLS> >::iterator it;
   for (it=_JPCLSs.begin(); it!=_JPCLSs.end(); ++it){
     if( fabs(lamX) > (*it)->J ) continue;
@@ -146,16 +146,7 @@ int evtNo=theData->evtNo;
     }
   }
 
-  if(_withDyn){
-    Vector4<double> mass4Vec(0.,0.,0.,0.);
-    std::vector<Particle*> fsParticleVec=_decay->finalStateParticles();
-
-    std::vector<Particle*>::iterator itPartVec;
-    for (itPartVec=fsParticleVec.begin(); itPartVec!=fsParticleVec.end(); ++itPartVec){
-       mass4Vec+=theData->FourVecsString[(*itPartVec)->name()];
-    }
-    result*=BreitWigner(mass4Vec, _currentXMass, _currentXWidth);
-  }
+  result*=_absDyn->eval(theData);
 
   if ( _cacheAmps){
 #ifdef _OPENMP
@@ -191,12 +182,8 @@ void  LSDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
   fitErr.Mags[_key]=currentMagErrMap;
   fitErr.Phis[_key]=currentPhiErrMap;
 
-  if(_withDyn){
-    fitVal.Masses[_massKey]=_decay->motherPart()->mass();
-    fitErr.Masses[_massKey]=0.03;
-    fitVal.Widths[_massKey]=_decay->motherPart()->width();
-    fitErr.Widths[_massKey]=0.2*_decay->motherPart()->width();
-  }
+  _absDyn->getDefaultParams(fitVal, fitErr);
+
 
   if(!_daughter1IsStable) _decAmpDaughter1->getDefaultParams(fitVal, fitErr);
   if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);  
@@ -226,16 +213,7 @@ bool LSDecAmps::checkRecalculation(fitParams& theParamVal){
      }
    }
 
-   if(_withDyn){
-     double mass=theParamVal.Masses[_massKey];
-     if ( fabs(mass-_currentXMass) > 1.e-10){
-       _recalculate=true;
-     }
-     double width=theParamVal.Widths[_massKey];
-     if ( fabs(width-_currentXWidth) > 1.e-10){
-       _recalculate=true;
-     }
-   }
+   if(_absDyn->checkRecalculation(theParamVal)) _recalculate=true; 
 
    if(!_daughter1IsStable) {
      if(_decAmpDaughter1->checkRecalculation(theParamVal)) _recalculate=true;
@@ -259,12 +237,7 @@ void  LSDecAmps::updateFitParams(fitParams& theParamVal){
      _currentParamPhis[*it]=thePhi;
    }
 
-   if (_withDyn){
-     double xMass=theParamVal.Masses[_massKey];
-     _currentXMass= xMass;
-     double xWidth=theParamVal.Widths[_massKey];
-     _currentXWidth=xWidth;
-   }
+   _absDyn->updateFitParams(theParamVal);
 
   if(!_daughter1IsStable) _decAmpDaughter1->updateFitParams(theParamVal);
   if(!_daughter2IsStable) _decAmpDaughter2->updateFitParams(theParamVal);

PwaUtils/WoDynamics.cc

+// WoDynamics class definition file. -*- C++ -*-
+// Copyright 20123Bertram Kopf
+
+#include <getopt.h>
+#include <fstream>
+#include <string>
+
+#include "PwaUtils/WoDynamics.hh"
+#include "ErrLogger/ErrLogger.hh"
+#include "Particle/Particle.hh"
+
+WoDynamics::WoDynamics(std::vector<Particle*>& fsParticles, Particle* mother) :
+  AbsDynamics(fsParticles, mother)
+{
+}
+
+WoDynamics::~WoDynamics()
+{
+}
+
+complex<double> WoDynamics::eval(EvtData* theData, Spin OrbMom){
+
+  complex<double> result(1.,0.);
+  return result;
+}
+
+void  WoDynamics::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
+}
+
+bool WoDynamics::checkRecalculation(fitParams& theParamVal){
+  _recalculate=false;
+  return _recalculate;
+}
+
+void WoDynamics::updateFitParams(fitParams& theParamVal){
+}

PwaUtils/WoDynamics.hh

+// WoDynamics class definition file. -*- C++ -*-
+// Copyright 2013 Bertram Kopf
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <complex>
+#include <map>
+#include <string>
+#include <boost/shared_ptr.hpp>
+
+#include "PwaUtils/AbsDynamics.hh"
+
+
+class WoDynamics : public AbsDynamics{
+
+public:
+  WoDynamics(std::vector<Particle*>& fsParticles, Particle* mother);
+  virtual ~WoDynamics();
+
+  virtual complex<double> eval(EvtData* theData, Spin OrbMom=0);
+  
+  virtual void getDefaultParams(fitParams& fitVal, fitParams& fitErr);
+  virtual bool checkRecalculation(fitParams& theParamVal);
+  virtual void updateFitParams(fitParams& theParamVal);
+
+protected:
+
+private:
+
+};