initial import of subdir Examples/Tutorial

Examples/Jamfile

@@ -6,3 +6,4 @@ build-project Psi2STo2K2PiGam ;
 build-project pbarp ;
 build-project SeqDecHeliAmps ; 
 build-project pbarpToOmegaPiLS ;
+build-project Tutorial ;

Examples/Tutorial/DfuncClebschG/ClebschGApp.cc

+#include <iostream>
+#include <cstring>
+#include <string>
+#include <sstream>
+
+
+#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+int main(int __argc,char *__argv[]){
+
+  ErrLogger::instance()->setLevel(log4cpp::Priority::DEBUG);
+  
+  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
+                    || strcmp( __argv[1], "--help" ) == 0 ) ){
+
+    Info << "\nThis is a test application calculates the Clebsch-Gordan coefficients C(j1, lambda_j1, j2, lambda_j2, J, lambda_J) with lambda_ji= -ji, -ji+1, ... , ji\n"
+	 << "The switches are:\n\n"
+         << "-2j1  (default 1)\n\n" 
+         << "-2j2  (default 3)\n\n"
+         << "-2J   (default 4)"
+         << endmsg;
+    return 0;
+  }
+
+  std::string twoj1Str="1";
+  std::string twoj2Str="3";
+  std::string twoJStr="4";
+  
+  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
+    bool found=false;
+    std::string sw = __argv[optind];
+    if (sw=="-2j1"){
+      optind++;
+      twoj1Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-2j2"){
+      optind++;
+      twoj2Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-2J"){
+      optind++;
+      twoJStr = __argv[optind];
+      found=true;
+    }
+    if (!found){
+      Warning << "Unknown switch: " 
+            << __argv[optind] << endmsg;
+      optind++;
+    }
+    
+    while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
+    }
+
+  std::stringstream twoj1StrStr(twoj1Str);
+  unsigned int twoj1=0;
+  twoj1StrStr >> twoj1;
+
+  std::stringstream twoj2StrStr(twoj2Str);
+  unsigned int twoj2=0;
+  twoj2StrStr >> twoj2;
+
+  std::stringstream twoJStrStr(twoJStr);
+  unsigned int twoJ=0;
+  twoJStrStr >> twoJ;
+
+  Spin j1( double(twoj1)/2. );
+  Spin j2( double(twoj2)/2. );
+  Spin J( double(twoJ)/2. );
+
+  if (J> j1+j2){
+    Alert << "J< j1+j2 !!!\t" << J << " < " << j1+j2 << endmsg;
+    exit(0);
+  }
+  
+  if (J< fabs(j1-j2)){
+    Alert << "J< |j1-j2| \t" << J << " < " << fabs(j1-j2) << endmsg;
+    exit(0);
+  }
+  if ( (twoj1+twoj2)%2 != twoJ%2){
+    Alert << "j1= "<< j1 << " j2= " << j2 << "  cannot couple to J= " << J << endmsg;
+    exit(0);
+  }
+
+
+  for (Spin lamj1=-j1;lamj1<=j1; ++lamj1){
+    for (Spin lamj2=-j2;lamj2<=j2; ++lamj2){
+      Spin lamJ=lamj1+lamj2;
+      if (fabs(lamJ) > J) continue;
+      Info << "C(j1=" << j1 << ", lamj1=" << lamj1 
+	   << ", j2=" << j2 << ", lamj2=" << lamj2
+	   << ", J=" << J << ", lamJ=" << lamJ << ")\t"
+	   << Clebsch(j1, lamj1, j2, lamj2, J, lamJ)
+	   << "\t  -> sqr= " 
+	   << Clebsch(j1, lamj1, j2, lamj2, J, lamJ)*Clebsch(j1, lamj1, j2, lamj2, J, lamJ)
+	   << endmsg;
+
+    }
+  }
+
+
+  return 0;
+}
+

Examples/Tutorial/DfuncClebschG/JPCdecays.cc

+#include <getopt.h>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include "Examples/Tutorial/DfuncClebschG/JPCdecays.hh"
+#include "PwaUtils/DataUtils.hh"
+
+
+
+#include "ErrLogger/ErrLogger.hh"
+
+JPCdecays::JPCdecays(boost::shared_ptr<const jpcRes> motherJPC, boost::shared_ptr<const jpcRes> daughter1JPC, boost::shared_ptr<const jpcRes> daughter2JPC, bool sameDaughters) :
+  _motherJPC(motherJPC)
+  ,_daughter1JPC(daughter1JPC)
+  ,_daughter2JPC(daughter2JPC)
+  ,_sameDaughters(sameDaughters)
+{
+  _allLSs=GetValidLS(_motherJPC->J, _motherJPC->P, _daughter1JPC->J,_daughter1JPC->P, _daughter2JPC->J, _daughter2JPC->P);
+
+
+  //check C parity
+  bool allowedC=true;
+  if (!sameDaughters){
+    if (_motherJPC->C != (daughter1JPC->C * daughter2JPC->C)) allowedC=false;
+  }
+
+  std::vector<LS>::const_iterator itLS;
+
+  for (itLS=_allLSs.begin(); itLS!=_allLSs.end(); ++itLS){
+    if (sameDaughters){
+      bool bosons=true;
+      if (int(2*daughter1JPC->J)%2 !=0) bosons=false;
+      
+      Spin theL=itLS->L;
+      Spin theS=itLS->S;
+      int cparity=0;
+
+      if (bosons) cparity=pow(-1,int(theL));
+      else cparity=pow(-1,int(theL+theS)); //daughters are fermions	
+
+      if(_motherJPC->C != cparity ) allowedC=false;
+    }
+
+    if (allowedC) _allowedCParityLSs.push_back(*itLS);
+    else _forbiddenCParityLSs.push_back(*itLS);
+  }
+  
+}
+
+JPCdecays::~JPCdecays()
+{
+
+}
+
+void JPCdecays::print(std::ostream& os) const{
+
+  std::cout << "the following LS combinations are allowed for the decay: ";
+  _motherJPC->print(os);
+  std::cout << " -> ";
+  _daughter1JPC->print(os);
+  std::cout << " + ";
+  _daughter2JPC->print(os);
+
+  std::cout << "\n __________LS combinations with allowed C parity______________" << std::endl;
+
+  std::vector<LS>::const_iterator itLS;
+
+  for (itLS=_allowedCParityLSs.begin(); itLS!=_allowedCParityLSs.end(); ++itLS){
+    std::cout << "L: " << itLS->L << "\tS: " << itLS->S << std::endl;
+  }
+
+ std::cout << "\n __________LS combinations with forbidden C parity______________" << std::endl;
+
+  for (itLS=_forbiddenCParityLSs.begin(); itLS!=_forbiddenCParityLSs.end(); ++itLS){
+    std::cout << "L: " << itLS->L << "\tS: " << itLS->S << std::endl;
+  }
+}

Examples/Tutorial/DfuncClebschG/JPCdecays.hh

+#ifndef _JPCdecays_H
+#define _JPCdecays_H
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include <cassert>
+
+#include <boost/shared_ptr.hpp>
+
+#include "TROOT.h"
+// #include <TSystem.h>
+#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
+
+class jpcRes;
+
+class JPCdecays {
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor 
+  JPCdecays(boost::shared_ptr<const jpcRes> motherJPC, boost::shared_ptr<const jpcRes> daughter1JPC, boost::shared_ptr<const jpcRes> daughter2JPC, bool sameDaughters=false);
+
+
+  /** Destructor */
+  virtual ~JPCdecays();
+
+  // Getters:
+  void print(std::ostream& os) const; 
+
+protected:
+
+
+private:
+  boost::shared_ptr<const jpcRes> _motherJPC;
+  boost::shared_ptr<const jpcRes> _daughter1JPC;
+  boost::shared_ptr<const jpcRes> _daughter2JPC;
+  std::vector<LS> _allLSs;
+  std::vector<LS> _allowedCParityLSs;
+  std::vector<LS> _forbiddenCParityLSs;
+  bool _sameDaughters;
+};
+
+#endif

Examples/Tutorial/DfuncClebschG/JPCdecaysApp.cc

+#include <iostream>
+#include <cstring>
+#include <string>
+#include <sstream>
+
+#include "Examples/Tutorial/DfuncClebschG/JPCdecays.hh"
+#include "PwaUtils/DataUtils.hh"
+#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+bool findPCAndRemove(std::string& theString, int& theP, int& theC){
+
+  bool result=true;
+
+  std::string plusplusStr("++");
+  std::string plusminusStr("+-");
+  std::string minusplusStr("-+");
+  std::string minusminusStr("--");
+
+  size_t strfound=0;
+  if ( (strfound=theString.find(plusplusStr)) != std::string::npos){
+   theString.erase (strfound, strfound+2);
+   theP=1;
+   theC=1;
+  }
+  else  if ( (strfound=theString.find(plusminusStr)) != std::string::npos){
+   theString.erase (strfound, strfound+2);
+   theP=1;
+   theC=-1;
+  }
+  else if ( (strfound=theString.find(minusplusStr)) != std::string::npos){
+   theString.erase (strfound, strfound+2);
+   theP=-1;
+   theC=1;
+  }
+  else if ( (strfound=theString.find(minusminusStr)) != std::string::npos){
+   theString.erase (strfound, strfound+2);
+   theP=-1;
+   theC=-1;
+  }
+  else result=false;
+
+  return result;
+}
+
+
+
+int main(int __argc,char *__argv[]){
+
+  ErrLogger::instance()->setLevel(log4cpp::Priority::DEBUG);
+  
+  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
+                    || strcmp( __argv[1], "--help" ) == 0 ) ){
+
+    Info << "\nThis is a test application which calculates all allowed LS combinations for the decay a(JPC) -> b(JPC) + c(JPC)\n"
+	 << "The switches are:\n\n"
+         << "-mother  (default 1+-)\n\n" 
+         << "-daughter1  (default 2--)\n\n"
+         << "-daughter2   (default 1-+)"
+         << endmsg;
+    return 0;
+  }
+
+  std::string motherStr="1+-";
+  std::string daughter1Str="2--";
+  std::string daughter2Str="1-+";
+  
+  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
+    bool found=false;
+    std::string sw = __argv[optind];
+    if (sw=="-mother"){
+      optind++;
+      motherStr = __argv[optind];
+      found=true;
+    }
+    if (sw=="-daughter1"){
+      optind++;
+      daughter1Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-daughter2"){
+      optind++;
+      daughter2Str = __argv[optind];
+      found=true;
+    }
+    if (!found){
+      Warning << "Unknown switch: " 
+            << __argv[optind] << endmsg;
+      optind++;
+    }
+    
+    while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
+    }
+
+  Info << "the JPC combinations of the particles are:\n\n"
+       << "mother:\t" << motherStr
+       << "\ndaughter1:\t" << daughter1Str
+       << "\ndaughter2:\t" << daughter2Str << endmsg;
+
+
+  int motherP=0;
+  int motherC=0;
+  int daughter1P=0;
+  int daughter1C=0;
+  int daughter2P=0;
+  int daughter2C=0;
+
+  bool found=findPCAndRemove(motherStr, motherP, motherC);
+
+  if (!found) {
+    Alert << "JPCmother: " << motherStr << " not allowed!!!" << endmsg;
+    exit(0);
+  }
+
+  std::stringstream motherJStrStr(motherStr);
+  double motherJ=0.;
+  motherJStrStr >> motherJ;
+  Info << "mother:\tJ: " << motherJ << "\tP: " << motherP << "\tC: " << motherC << endmsg;
+
+  found=findPCAndRemove(daughter1Str, daughter1P, daughter1C);
+  if (!found) {
+    Alert << "JPCdaughter1: " << daughter1Str << " not allowed!!!" << endmsg;
+    exit(0);
+  }
+  std::stringstream daughter1JStrStr(daughter1Str);
+  double daughter1J=0.;
+  daughter1JStrStr >> daughter1J;
+  Info << "daughter1:\tJ: " << daughter1J << "\tP: " << daughter1P << "\tC: " << daughter1C << endmsg;
+
+  found=findPCAndRemove(daughter2Str, daughter2P, daughter2C);
+  if (!found) {
+    Alert << "JPCdaughter2: " << daughter2Str << " not allowed!!!" << endmsg;
+    exit(0);
+  }
+  std::stringstream daughter2JStrStr(daughter2Str);
+  double daughter2J=0.;
+  daughter2JStrStr >> daughter2J;
+  Info << "daughter2:\tJ: " << daughter2J << "\tP: " << daughter2P << "\tC: " << daughter2C << endmsg;
+
+
+  
+  boost::shared_ptr<const jpcRes> motherJPC(new jpcRes(Spin(motherJ), motherP, motherC));
+  boost::shared_ptr<const jpcRes> daughter1JPC(new jpcRes(Spin(daughter1J), daughter1P, daughter1C)); 
+  boost::shared_ptr<const jpcRes> daughter2JPC(new jpcRes(Spin(daughter2J), daughter2P, daughter2C));
+
+  JPCdecays theDecays(motherJPC, daughter1JPC, daughter2JPC);
+  theDecays.print(std::cout);
+  return 0;
+}
+

Examples/Tutorial/DfuncClebschG/Jamfile

+project : 
+        ;
+
+lib DfuncClebschG : 
+        [ glob *.cc : *App.cc ]  
+        $(TOP)/qft++//qft++ 
+        $(TOP)/ErrLogger//ErrLogger
+	: <use>$(TOP)//Minuit2 <use>$(TOP)//Geneva 
+	: 
+	: <library>$(TOP)//Minuit2 <library>$(TOP)//Geneva ;
+
+exe ClebschGApp : ClebschGApp.cc DfuncClebschG : ;
+exe JPCdecaysApp : JPCdecaysApp.cc DfuncClebschG : ;

Examples/Tutorial/Jamfile

+
+build-project LineShapes ;
+build-project DfuncClebschG ;

Examples/Tutorial/LineShapes/BwShape.cc

+#include <getopt.h>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include "Examples/Tutorial/LineShapes/BwShape.hh"
+
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TMath.h"
+
+
+#include "ErrLogger/ErrLogger.hh"
+
+BwShape::BwShape(double MassRes, double MassWidth, double MassDec1, double MassDec2, unsigned int Lmax, double deltaMass) :
+  _theTFile(0)
+{
+
+  if (MassRes <= MassDec1+MassDec2){
+    Alert << "Mass of the resonance smaller than the masses of the decay particles\n" 
+	  << MassRes << " < " << MassDec1+MassDec2 << endmsg;
+    exit(0); 
+  }
+
+   std::stringstream Lmaxstrstr;
+   Lmaxstrstr << Lmax;
+   std::string rootFileName="./BwShapeLmax"+Lmaxstrstr.str()+".root";
+   _theTFile=new TFile(rootFileName.c_str(),"recreate");
+ 
+  for (unsigned int i=0; i<=Lmax; ++i){
+
+    std::stringstream Lstrstr;
+    Lstrstr << i;
+    std::string histName="BreitWigner_L"+Lstrstr.str();
+
+
+   _histMap[i]= new TH1F(histName.c_str(),histName.c_str(),301, MassRes-deltaMass, MassRes+deltaMass);
+
+    histName="Argand_L"+Lstrstr.str();
+    TH2F* currentArgandHist=new TH2F(histName.c_str(),histName.c_str(),301, -1., 1., 301, 0., 1.3);
+    currentArgandHist->SetXTitle("Re(Bw)");
+    currentArgandHist->SetYTitle("Im(Bw)");
+    _argandHistMap[i]=currentArgandHist;
+    
+    
+    //initialize here the phase shift histogramms _phaseHistMap
+  }
+  
+  double stepSize=2.*deltaMass/300;
+  
+  for (unsigned int lIt=0; lIt<=Lmax; ++lIt){
+    
+    TH1F* currentHist=_histMap[lIt];
+    TH2F* currentArgandHist=_argandHistMap[lIt];
+
+    for (double massIt=MassRes-deltaMass; massIt<MassRes+deltaMass; massIt+=stepSize){
+      
+      Vector4<double>  res4V(massIt, 0., 0., 0.);
+      
+      complex<double> currentBW=BreitWignerBlattW(res4V, MassDec1, MassDec2, MassRes, MassWidth, lIt);
+      double weight=norm(currentBW);
+      currentHist->Fill(massIt,weight);
+      
+      currentArgandHist->Fill(currentBW.real(),currentBW.imag());
+      
+      //fill here the phase shift histogramms: hint "atan2(imag,real)"
+
+    }
+
+  }
+}
+
+BwShape::~BwShape()
+{
+   _theTFile->Write();
+   _theTFile->Close();
+}
+

Examples/Tutorial/LineShapes/BwShape.hh

+#ifndef _BwShape_H
+#define _BwShape_H
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <map>
+#include <cassert>
+
+#include <boost/shared_ptr.hpp>
+
+#include "TROOT.h"
+#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
+
+class TFile;
+class TH1F;
+class TH2F;
+
+class BwShape {
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor 
+  BwShape(double MassRes, double MassWidth, double MassDec1, double MassDec2, unsigned int Lmax, double deltaMass);
+
+
+  /** Destructor */
+  virtual ~BwShape();
+
+  // Getters:
+ 
+protected:
+
+
+private:
+  TFile* _theTFile;
+  std::map <unsigned int, TH1F* > _histMap;
+  std::map <unsigned int, TH2F* > _argandHistMap;
+  std::map <unsigned int, TH2F* > _phaseHistMap;
+};
+
+#endif

Examples/Tutorial/LineShapes/BwShapeApp.cc

+#include <iostream>
+#include <cstring>
+#include <string>
+#include <sstream>
+#include <vector>
+#include <map>
+
+#include <boost/shared_ptr.hpp>
+
+#include "Examples/Tutorial/LineShapes/BwShape.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+int main(int __argc,char *__argv[]){
+  ErrLogger::instance()->setLevel(log4cpp::Priority::DEBUG);
+  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
+                    || strcmp( __argv[1], "--help" ) == 0 ) ){
+
+    Info << "\nThis is a test application for histogramming the Argand plot, the phase shift and the mass shape of the relativistic Breit-Wigner function with Blatt-Weisskopf barrier factors\n"
+	 << "The switches are:\n\n"
+         << "-mass (mass of the resonance;default 1.318)\n\n" 
+         << "-width (width of the resonance;default 0.1)\n\n"
+         << "-mDec1 (mass of the first decay particle; default 0.13498\n\n"
+         << "-mDec2 (mass of the second decay particle; default 0.547853\n\n"
+         << "-Lmax (histogams of mass shapes; default 3\n\n"
+         << "-DeltaM (range of histograms +-DeltaM; default 0.4\n\n"  
+         << endmsg;
+    return 0;
+  }
+  std::string massStr="1.318";
+  std::string widthStr="0.1";
+  std::string mDec1Str="0.13498";
+  std::string mDec2Str="0.547853";
+  std::string LmaxStr="3";
+  std::string DeltaMStr="0.4";
+
+  
+  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
+    bool found=false;
+    std::string sw = __argv[optind];
+    if (sw=="-mass"){
+      optind++;
+      massStr = __argv[optind];
+      found=true;
+    }
+    if (sw=="-width"){
+      optind++;
+      widthStr = __argv[optind];
+      found=true;
+    }
+    if (sw=="-mDec1"){
+      optind++;
+      mDec1Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-mDec2"){
+      optind++;
+      mDec2Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-Lmax"){
+      optind++;
+      LmaxStr = __argv[optind];
+      found=true;
+    }
+    if (sw=="-DeltaM"){
+      optind++;
+      DeltaMStr = __argv[optind];
+      found=true;
+    }
+
+    if (!found){
+      Warning << "Unknown switch: " 
+            << __argv[optind] << endmsg;
+      optind++;
+    }
+    
+    while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
+    }
+
+
+  std::stringstream massStrStr(massStr);
+  double mass=0.;
+  massStrStr >> mass;
+
+  std::stringstream widthStrStr(widthStr);
+  double width=0.;
+  widthStrStr >> width;
+
+  std::stringstream mDec1StrStr(mDec1Str);
+  double mDec1=0.;
+  mDec1StrStr >> mDec1;
+
+  std::stringstream mDec2StrStr(mDec2Str);
+  double mDec2=0.;
+  mDec2StrStr >> mDec2;
+
+  std::stringstream LmaxStrStr(LmaxStr);
+  unsigned int Lmax=0;
+  LmaxStrStr >> Lmax;
+
+  std::stringstream DeltaMStrStr(DeltaMStr);
+  double DeltaM=0;
+  DeltaMStrStr >> DeltaM;
+
+  BwShape bwShape(mass, width, mDec1, mDec2, Lmax, DeltaM);
+
+  return 0;
+}
+

Examples/Tutorial/LineShapes/FlatteShape.cc

+#include <getopt.h>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include "Examples/Tutorial/LineShapes/FlatteShape.hh"
+
+#include "TFile.h"
+#include "TH1F.h"
+#include "TMath.h"
+
+
+#include "ErrLogger/ErrLogger.hh"
+
+FlatteShape::FlatteShape(std::string ptype, double g1, double g2) :
+  _theTFile(0)
+  ,_histShapeLow(0)
+  ,_histShapeHigh(0)
+  ,_pType(ptype)
+{
+  if(_pType!="a0" && _pType!="f0"){
+    Alert <<"only a0 and f0 are supported" << endmsg;
+    exit(1);
+  }
+  if(g1<0.  || g2<0.){
+    Alert <<"requirement: g1>0. and g2>0." << endmsg;
+    exit(1);
+  }
+
+  std::string rootFileName="./FlatteShape"+_pType+".root";
+  _theTFile=new TFile(rootFileName.c_str(),"recreate");
+
+
+  double mass1=0.;
+  double mass2=0.;
+
+  if(_pType=="a0"){
+    mass1 = 0.1349766;
+    mass2 = 0.547853;
+  }
+  else if(_pType=="f0"){
+    mass1 = 0.1349766;
+    mass2= 0.1349766;
+  }
+
+  std::pair <const double, const double> decPairLow=make_pair(mass1, mass2); 
+  const double KplusMass = 0.493677;
+  const double K0Mass = 0.497614;
+  std::pair <const double, const double> decPairHigh=make_pair(KplusMass,K0Mass);
+
+
+  int size=700;
+  double massMin=0.35;
+  double massMax=1.85;
+
+//   double massMin=1.;
+//   double massMax=1.004;
+  double stepSize=(massMax-massMin)/size;
+
+  _histShapeLow= new TH1F("_histShapeLow","hist low",size+1, massMin, massMax);
+  _histShapeHigh= new TH1F("_histShapeHigh","hist high",size+1, massMin, massMax);
+
+  for (double mass=massMin; mass<massMax; mass+=stepSize){
+    Vector4<double> mass4Vec(mass, 0.,0.,0.);
+
+    _histShapeLow->Fill(mass4Vec.M(), norm(Flatte(mass4Vec, decPairLow, decPairHigh, 0.981, g1, g2)));
+    _histShapeHigh->Fill(mass4Vec.M(), norm(Flatte(mass4Vec, decPairHigh, decPairLow, 0.981, g2, g1))); 
+  }
+
+}
+
+FlatteShape::~FlatteShape()
+{
+  _theTFile->Write();
+  _theTFile->Close();
+}
+

Examples/Tutorial/LineShapes/FlatteShape.hh

+#ifndef _FlatteShape_H
+#define _FlatteShape_H
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include <cassert>
+
+#include <boost/shared_ptr.hpp>
+
+#include "TROOT.h"
+// #include <TSystem.h>
+#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
+
+class TFile;
+class TH1F;
+
+class FlatteShape {
+
+public:
+
+  // create/copy/destroy:
+
+  ///Constructor 
+  FlatteShape(std::string ptype, double g1, double g2);
+
+
+  /** Destructor */
+  virtual ~FlatteShape();
+
+  // Getters:
+ 
+protected:
+
+
+private:
+  TFile* _theTFile;
+  TH1F* _histShapeLow;
+  TH1F* _histShapeHigh;
+
+  std::string _pType;
+
+};
+
+#endif

Examples/Tutorial/LineShapes/FlatteShapeApp.cc

+#include <iostream>
+#include <cstring>
+#include <string>
+#include <sstream>
+#include <vector>
+#include <map>
+
+#include <boost/shared_ptr.hpp>
+
+#include "Examples/Tutorial/LineShapes/FlatteShape.hh"
+#include "ErrLogger/ErrLogger.hh"
+
+int main(int __argc,char *__argv[]){
+  ErrLogger::instance()->setLevel(log4cpp::Priority::DEBUG);
+  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
+                    || strcmp( __argv[1], "--help" ) == 0 ) ){
+
+    Info << "\nThis is a test application for histogramming the Flatte shape\n"
+	 << "The switches are:\n\n"
+         << "-ptype particle type (so far: a0, f0) (default a0)\n\n" 
+         << "-g11 g(pi,eta) for a0; or g(pi,pi) for f0 (default 0.5857)\n\n"
+         << "-g22 g(K,K) for a0 and f0 (default 1.86243)\n\n"
+         << "-mass for a0 or f0  (default 0.982)\n"  
+         << endmsg;
+    return 0;
+  }
+
+  std::string ptypeStr="a0";
+  std::string g11Str="0.5857";
+  std::string g22Str="1.86243";
+  std::string massStr="0.982";
+  
+  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
+    bool found=false;
+    std::string sw = __argv[optind];
+    if (sw=="-ptype"){
+      optind++;
+      ptypeStr = __argv[optind];
+      found=true;
+    }
+    if (sw=="-g11"){
+      optind++;
+      g11Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-g22"){
+      optind++;
+      g22Str = __argv[optind];
+      found=true;
+    }
+    if (sw=="-mass"){
+      optind++;
+      massStr = __argv[optind];
+      found=true;
+    }
+    if (!found){
+      Warning << "Unknown switch: " 
+            << __argv[optind] << endmsg;
+      optind++;
+    }
+    
+    while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
+    }
+
+  std::stringstream g11StrStr(g11Str);
+  double g11=0.;
+  g11StrStr >> g11;
+
+  std::stringstream g22StrStr(g22Str);
+  double g22=0.;
+  g22StrStr >> g22;
+
+  std::stringstream massStrStr(massStr);
+  double mass=0.;
+  massStrStr >> mass;
+  
+
+//   std::string resStr="a0";
+//  std::string resStr="f0";
+
+
+  FlatteShape flatteShape(ptypeStr, g11, g22);
+
+  return 0;
+}
+

Examples/Tutorial/LineShapes/Jamfile

+project : 
+        ;
+
+lib LineShapes : 
+        [ glob *.cc : *App.cc ]  
+        $(TOP)/qft++//qft++ 
+        $(TOP)/ErrLogger//ErrLogger
+	: <use>$(TOP)//Minuit2 <use>$(TOP)//Geneva 
+	: 
+	: <library>$(TOP)//Minuit2 <library>$(TOP)//Geneva ;
+
+exe FlatteShapeApp : FlatteShapeApp.cc LineShapes : ;
+exe BwShapeApp : BwShapeApp.cc LineShapes : ;

qft++/relativistic-quantum-mechanics/BlattWeisskopf.cc

@@ -79,10 +79,12 @@ double BlattWeisskopf::compute(double p) const
       if(1 == _LL) return sqrt(x/(1.0+x));
       else
 	if(2 == _LL) return sqrt((13.*x*x)/((x-3.)*(x-3.)+9.*x));
-	else {
-	  std::cout << "Angular momentum " << _LL << " not implemented" << std::endl;
-	  assert(0);
-	}
+	else
+	  if(3 == _LL) return sqrt((277.*x*x*x)/(x*(x-15.)*(x-15.)+9.*(2.*x-5)*(2.*x-5)));
+	  else {
+	    std::cout << "Angular momentum " << _LL << " not implemented" << std::endl;
+	    assert(0);
+	  }
   }
 
 }