diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -7,6 +7,7 @@
GenFit3KS
GenFit3pi
GenFitBelleCPKpipi
+ GenFitDpipi
GenFitDs2KKpi
GenFitEFKLLM
GenFitKpipi
diff --git a/examples/GenFitDpipi.cc b/examples/GenFitDpipi.cc
new file mode 100644
--- /dev/null
+++ b/examples/GenFitDpipi.cc
@@ -0,0 +1,243 @@
+
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+#include "TFile.h"
+#include "TH2.h"
+#include "TString.h"
+#include "TTree.h"
+
+#include "LauSimpleFitModel.hh"
+#include "LauResonanceMaker.hh"
+#include "LauBkgndDPModel.hh"
+#include "LauDaughters.hh"
+#include "LauEffModel.hh"
+#include "LauIsobarDynamics.hh"
+#include "LauMagPhaseCoeffSet.hh"
+#include "LauRealImagCoeffSet.hh"
+#include "LauRandom.hh"
+#include "LauVetoes.hh"
+#include "LauAbsModIndPartWave.hh"
+#include "LauModIndPartWaveRealImag.hh"
+#include "LauModIndPartWaveMagPhase.hh"
+
+void usage( std::ostream& out, const TString& progName )
+{
+ out<<"Usage:\n";
+ out<<progName<<" gen [nExpt = 1] [firstExpt = 0]\n";
+ out<<"or\n";
+ out<<progName<<" fit <iFit> [nExpt = 1] [firstExpt = 0]"<<std::endl;
+ out<<"or\n";
+ out<<progName<<" plot"<<std::endl;
+}
+
+int main( int argc, char** argv )
+{
+ // Process command-line arguments
+ // Usage:
+ // ./GenFit3pi gen [nExpt = 1] [firstExpt = 0]
+ // or
+ // ./GenFit3pi fit <iFit> [nExpt = 1] [firstExpt = 0]
+ if ( argc < 2 ) {
+ usage( std::cerr, argv[0] );
+ return EXIT_FAILURE;
+ }
+
+ TString command = argv[1];
+ command.ToLower();
+ Int_t iFit(0);
+ Int_t nExpt(1);
+ Int_t firstExpt(0);
+ if ( command == "gen" ) {
+ if ( argc > 2 ) {
+ nExpt = atoi( argv[2] );
+ if ( argc > 3 ) {
+ firstExpt = atoi( argv[3] );
+ }
+ }
+ } else if ( command == "fit" ) {
+ if ( argc < 3 ) {
+ usage( std::cerr, argv[0] );
+ return EXIT_FAILURE;
+ }
+ iFit = atoi( argv[2] );
+ if ( argc > 3 ) {
+ nExpt = atoi( argv[3] );
+ if ( argc > 4 ) {
+ firstExpt = atoi( argv[4] );
+ }
+ }
+ } else {
+ usage( std::cerr, argv[0] );
+ return EXIT_FAILURE;
+ }
+
+ LauRandom::setSeed(0);
+
+ // If you want to use square DP histograms for efficiency,
+ // backgrounds or you just want the square DP co-ordinates
+ // stored in the toy MC ntuple then set this to kTRUE
+ Bool_t squareDP = kTRUE;
+
+ // This defines the DP => decay is B+ -> pi+ pi+ D-
+ // Particle 1 = pi+
+ // Particle 2 = pi+
+ // Particle 3 = D-
+ // The DP is defined in terms of m13Sq and m23Sq
+ LauDaughters* daughters = new LauDaughters("B+", "pi+", "pi+", "D-", squareDP);
+
+ // Optionally apply some vetoes to the DP
+ LauVetoes* vetoes = new LauVetoes();
+
+ LauEffModel* effModel = new LauEffModel(daughters, vetoes);
+
+ // Set the values of the Blatt-Weisskopf barrier radii and whether they are fixed or floating
+ LauResonanceMaker& resMaker = LauResonanceMaker::get();
+ resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 4.0 );
+ resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Charm, 4.0 );
+ resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 );
+ resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Beauty, 4.0 );
+ resMaker.fixBWRadius( LauBlattWeisskopfFactor::Parent, kTRUE);
+ resMaker.fixBWRadius( LauBlattWeisskopfFactor::Charm, kTRUE);
+ resMaker.fixBWRadius( LauBlattWeisskopfFactor::Light, kTRUE);
+ resMaker.fixBWRadius( LauBlattWeisskopfFactor::Beauty, kTRUE);
+ // Create the isobar model
+
+ std::vector<LauAbsCoeffSet*> coeffset;
+ LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel);
+ LauAbsResonance* reson(0);
+
+ reson = sigModel->addResonance("D*0", 2, LauAbsResonance::RelBW);
+ reson = sigModel->addResonance("D*0_0", 2, LauAbsResonance::MIPW_MagPhase);
+
+ //vector of knot masses - ignore ends and let Lauura deal with them
+ std::vector<double> knot_mass;
+
+ knot_mass.push_back(2.100);
+ knot_mass.push_back(2.200);
+ knot_mass.push_back(2.300);
+ knot_mass.push_back(2.400);
+ knot_mass.push_back(2.500);
+ knot_mass.push_back(2.600);
+ knot_mass.push_back(2.700);
+ knot_mass.push_back(2.800);
+ knot_mass.push_back(2.900);
+ knot_mass.push_back(3.100);
+ knot_mass.push_back(4.100);
+
+ //find knot at 2.4 within vector - this knot will be fixed
+ //fill set with values from vector of knots
+ int fixed_knot=-1;
+ std::set<double> knot_mass_;
+ for (Size_t knot(0); knot < knot_mass.size(); knot++){
+ knot_mass_.insert(knot_mass.at(knot));
+ if (knot_mass.at(knot) == 2.400) {
+ std::cout << knot << " " << knot_mass.at(knot) << std::endl;
+ fixed_knot = knot+1;
+ }
+ }
+
+ ((LauModIndPartWaveMagPhase*)reson)->defineKnots(knot_mass_);
+ int nKnots = (int)(((LauModIndPartWaveMagPhase*)reson)->nKnots());
+ std::cout << "there are " << nKnots << " knots." << std::endl;
+
+ ((LauModIndPartWaveMagPhase*)reson)->floatKnotsSecondStage(kFALSE);
+
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(0, 0.12, -2.82,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(1, 0.58, -1.56,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(2, 0.73, -1.00,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(3, 0.68, -0.42,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(4, 0.5, 0.0,kTRUE,kTRUE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(5, 0.23, -0.00,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(6, 0.23, -0.42,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(7, 0.15, -0.31,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(8, 0.17, -0.63,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(9, 0.20, -0.87,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(10, 0.14, -1.16,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(11, 0.08, 1.02,kFALSE,kFALSE);
+ ((LauModIndPartWaveMagPhase*)reson)->setKnotAmp(12, 0.0, 0.0,kTRUE, kTRUE);
+
+ reson = sigModel->addResonance("D*0_2", 2, LauAbsResonance::RelBW);
+ reson = sigModel->addResonance("D*0_1(2680)", 2, LauAbsResonance::RelBW);
+ reson = sigModel->addResonance("B*0", 2, LauAbsResonance::RelBW);
+ reson = sigModel->addResonance("D*0_3(2760)", 2, LauAbsResonance::RelBW);
+ reson = sigModel->addResonance("D0(3000)", 2, LauAbsResonance::RelBW);
+
+ sigModel->setASqMaxValue(1.0);
+
+ // Create the fit model
+ LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel);
+
+ coeffset.push_back( new LauMagPhaseCoeffSet("D*0", 0.55, -0.38, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("D*0_0", 1.26, -0.28, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("D*0_2", 1.00, 0.00, kTRUE, kTRUE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("D*0_1(2680)", 0.48, 2.47, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("B*0", 0.27, 0.14, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("D*0_3(2760)", 0.17, 0.01, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("D0(3000)", 0.08, -0.84, kFALSE, kFALSE) );
+
+ for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
+ fitModel->setAmpCoeffSet(*iter);
+ }
+
+ double nSig(50000);
+
+ TString sigEventsName = "signalEvents";
+ LauParameter* nSigEvents = new LauParameter(sigEventsName,nSig,-2.0*nSig,2.0*nSig,kTRUE);
+ fitModel->setNSigEvents(nSigEvents);
+
+ // Set the number of experiments to generate or fit and which
+ // experiment to start with
+ fitModel->setNExpts( nExpt, firstExpt );
+
+ // Switch on/off calculation of asymmetric errors.
+ fitModel->useAsymmFitErrors(kFALSE);
+
+ // Randomise initial fit values for the signal mode
+ fitModel->useRandomInitFitPars(kTRUE);
+
+ // Switch on/off Poissonian smearing of total number of events
+ fitModel->doPoissonSmearing(kTRUE);
+
+ // Switch on/off Extended ML Fit option
+ fitModel->doEMLFit(kFALSE);
+
+ // Switch on the two-stage fit (for the resonance parameters)
+ fitModel->twoStageFit(kFALSE);
+
+ // Generate toys from the fit model
+ TString toyName = "fitToyDpipi"; toyName += ".root";
+ //fitModel->compareFitData(100,toyName,"",kTRUE);
+
+ // Set the names of the files to read/write
+ TString dataFile("Dpipi_data");
+
+ TString treeName("DecayTree");
+ TString rootFileName("");
+ TString tableFileName("");
+ TString fitToyFileName("fitToyMC_");
+
+ if (command.Contains("fit")) {
+ rootFileName += command; rootFileName += iFit;
+ rootFileName += "_expt_"; rootFileName += firstExpt;
+ rootFileName += "-"; rootFileName += (firstExpt+nExpt-1);
+ rootFileName += ".root";
+ tableFileName = "fitResults_"; tableFileName += "_"; tableFileName += iFit;
+ fitToyFileName += iFit;
+ fitToyFileName += ".root";
+ } else {
+ rootFileName = "dummy.root";
+ tableFileName = "genResults";
+ }
+
+ std::cout << rootFileName << std::endl;
+
+ // Execute the generation/fit
+ if (command=="gen") {
+ fitModel->run( command, toyName, treeName, rootFileName, tableFileName );
+ }
+ else fitModel->run( command, toyName, treeName, rootFileName, tableFileName );
+
+ return EXIT_SUCCESS;
+}
diff --git a/src/LauResonanceMaker.cc b/src/LauResonanceMaker.cc
--- a/src/LauResonanceMaker.cc
+++ b/src/LauResonanceMaker.cc
@@ -382,6 +382,9 @@
resInfo_.push_back( neutral );
resInfo_.push_back( positve );
resInfo_.push_back( negatve );
+ // D(2680)
+ neutral = new LauResonanceInfo("D*0_1(2680)", 2.6811, 0.1867, 1, 0, LauBlattWeisskopfFactor::Charm );
+ resInfo_.push_back( neutral );
// D(2760)
//OLD-- neutral = new LauResonanceInfo("D0(2760)", 2.7633, 0.061, 1, 0 );
//OLD-- positve = new LauResonanceInfo("D+(2760)", 2.7697, 0.061, 1, 1 );
@@ -391,8 +394,10 @@
resInfo_.push_back( neutral );
resInfo_.push_back( positve );
resInfo_.push_back( negatve );
+ neutral = new LauResonanceInfo("D*0_3(2760)", 2.7755, 0.0953, 3, 0, LauBlattWeisskopfFactor::Charm );
+ resInfo_.push_back( neutral );
// D(2900)
- neutral = new LauResonanceInfo("D0(3000)", 3.214, 0.186, 0, 0, LauBlattWeisskopfFactor::Charm );
+ neutral = new LauResonanceInfo("D0(3000)", 3.214, 0.186, 0, 0, LauBlattWeisskopfFactor::Charm );
resInfo_.push_back( neutral );
// D(3400)
neutral = new LauResonanceInfo("D0(3400)", 3.4, 0.15, 0, 0, LauBlattWeisskopfFactor::Charm );