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,244 @@ + +#include +#include +#include + +#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< [nExpt = 1] [firstExpt = 0]"< [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 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 as the are dealt with internally + std::vector 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 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); + + //Set magnitude and phase for the knots (including the end points here) + ((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) ); // Still need a total mag/phase for the MIPW shape + 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::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 );