diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 66531ac..78a398a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,50 +1,79 @@ add_subdirectory(ProgOpts) +# List of files (without extensions) that need pre-processing before compiling +list(APPEND EXAMPLE_SOURCES_CONFIGURE + KMatrixDto3pi + ) + +# List of files (without extensions) that do not need pre-processing before compiling list(APPEND EXAMPLE_SOURCES B3piKMatrixMassProj B3piKMatrixPlots CalcChiSq GenFit3K GenFit3KS GenFit3pi GenFitBelleCPKpipi GenFitDpipi GenFitDs2KKpi GenFitEFKLLM GenFitIncoherent_Bs2KSKpi GenFitKpipi GenFitNoDP GenFitNoDPMultiDim GenFitTimeDep GenFitTimeDep_Bs2KSKpi GenFitTimeDep_Bd2Dpipi - KMatrixDto3pi KMatrixExample MergeDataFiles mixedSampleTest PlotKMatrixTAmp PlotResults point2PointTestSample QuasiFlatSqDalitz QuasiFlatSqDalitz-CustomMasses ResultsExtractor SimFitCoordinator SimFitTask SimFitTaskRooFit ) if(NOT LAURA_BUILD_ROOFIT_TASK) list(REMOVE_ITEM EXAMPLE_SOURCES SimFitTaskRooFit) endif() +# Configure (if appropriate), build and install examples +foreach( _example ${EXAMPLE_SOURCES_CONFIGURE}) + configure_file(${_example}.cc.in ${_example}.cc @ONLY) + add_executable(${_example} ${CMAKE_CURRENT_BINARY_DIR}/${_example}.cc) + target_link_libraries(${_example} PRIVATE Laura++ ProgramOptions) + install(TARGETS ${_example} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) +endforeach() + foreach( _example ${EXAMPLE_SOURCES}) add_executable(${_example} ${_example}.cc) target_link_libraries(${_example} PRIVATE Laura++ ProgramOptions) install(TARGETS ${_example} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) endforeach() # Also install the python script version of GenFit3pi configure_file(GenFit3pi.py.in GenFit3pi.py @ONLY) install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/GenFit3pi.py DESTINATION ${CMAKE_INSTALL_BINDIR}) +# Also install the files to configure the models, coeffs, etc. +# Some need pre-processing +list(APPEND EXAMPLE_CONF_CONFIGURE + KMatrixDto3pi.json + ) +foreach( _json ${EXAMPLE_CONF_CONFIGURE}) + configure_file(${_json}.in ${_json} @ONLY) +endforeach() +list(APPEND EXAMPLE_CONF + GenFit3pi.json + FOCUSD3pi.dat + ) +list(TRANSFORM EXAMPLE_CONF PREPEND ${CMAKE_CURRENT_SOURCE_DIR}/) +list(TRANSFORM EXAMPLE_CONF_CONFIGURE PREPEND ${CMAKE_CURRENT_BINARY_DIR}/) +list(APPEND EXAMPLE_CONF ${EXAMPLE_CONF_CONFIGURE}) +install(FILES ${EXAMPLE_CONF} DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}) diff --git a/examples/GenFit3pi.json b/examples/GenFit3pi.json new file mode 100644 index 0000000..284082c --- /dev/null +++ b/examples/GenFit3pi.json @@ -0,0 +1,203 @@ +{ + "coeffs": { + "coeffs": [ + { + "Mag": 1.0, + "MagBlind": false, + "MagFixed": true, + "MagSecondStage": false, + "Phase": 0.0, + "PhaseBlind": false, + "PhaseFixed": true, + "PhaseSecondStage": false, + "clone": false, + "name": "rho0(770)", + "type": "MagPhase" + }, + { + "Mag": 0.37, + "MagBlind": false, + "MagFixed": false, + "MagSecondStage": false, + "Phase": 1.99, + "PhaseBlind": false, + "PhaseFixed": false, + "PhaseSecondStage": false, + "clone": false, + "name": "rho0(1450)", + "type": "MagPhase" + }, + { + "Mag": 0.27, + "MagBlind": false, + "MagFixed": false, + "MagSecondStage": false, + "Phase": -1.59, + "PhaseBlind": false, + "PhaseFixed": false, + "PhaseSecondStage": false, + "clone": false, + "name": "f_0(980)", + "type": "MagPhase" + }, + { + "Mag": 0.53, + "MagBlind": false, + "MagFixed": false, + "MagSecondStage": false, + "Phase": 1.39, + "PhaseBlind": false, + "PhaseFixed": false, + "PhaseSecondStage": false, + "clone": false, + "name": "f_2(1270)", + "type": "MagPhase" + }, + { + "Mag": 0.54, + "MagBlind": false, + "MagFixed": false, + "MagSecondStage": false, + "Phase": -0.84, + "PhaseBlind": false, + "PhaseFixed": false, + "PhaseSecondStage": false, + "clone": false, + "name": "BelleNR_Swave", + "type": "MagPhase" + } + ], + "nCoeffs": 5 + }, + "model": { + "commonSettings": { + "disableGSScalingOfWidthByBWFactors": false, + "setBWBachelorRestFrame": "ResonanceFrame", + "setBWType": "BWPrimeBarrier", + "setDefaultBWRadius": [ + { + "category": "Parent", + "fix": true, + "value": 5.0 + }, + { + "category": "Light", + "fix": true, + "value": 4.0 + } + ], + "setSpinFormalism": "Zemach_P" + }, + "resonances": [ + { + "changeResonance": { + "mass": { + "fix": true, + "value": 0.77526 + }, + "spin": { + "value": 1 + }, + "width": { + "fix": true, + "value": 0.1478 + } + }, + "resName": "rho0(770)", + "resPairAmpInt": 1, + "resType": "GS" + }, + { + "changeResonance": { + "mass": { + "fix": true, + "value": 1.465 + }, + "spin": { + "value": 1 + }, + "width": { + "fix": true, + "value": 0.4 + } + }, + "resName": "rho0(1450)", + "resPairAmpInt": 1, + "resType": "RelBW" + }, + { + "changeResonance": { + "mass": { + "fix": true, + "value": 0.965 + }, + "spin": { + "value": 0 + }, + "width": { + "fix": true, + "value": 0.07 + } + }, + "parameters": [ + { + "float": false, + "name": "g1", + "value": 0.2 + }, + { + "float": false, + "name": "g2", + "value": 1.0 + } + ], + "resName": "f_0(980)", + "resPairAmpInt": 1, + "resType": "Flatte" + }, + { + "changeResonance": { + "mass": { + "fix": true, + "value": 1.2755 + }, + "spin": { + "value": 2 + }, + "width": { + "fix": true, + "value": 0.1867 + } + }, + "resName": "f_2(1270)", + "resPairAmpInt": 1, + "resType": "RelBW" + }, + { + "changeResonance": { + "mass": { + "fix": true, + "value": 0.0 + }, + "spin": { + "value": 0 + }, + "width": { + "fix": true, + "value": 0.0 + } + }, + "parameters": [ + { + "float": false, + "name": "alpha", + "value": 0.2 + } + ], + "resName": "BelleNR_Swave", + "resPairAmpInt": 0, + "resType": "BelleSymNRNoInter" + } + ] + } +} diff --git a/examples/GenFit3pi.py.in b/examples/GenFit3pi.py.in index 9b7785f..eab1735 100644 --- a/examples/GenFit3pi.py.in +++ b/examples/GenFit3pi.py.in @@ -1,228 +1,211 @@ #!/usr/bin/env python """ Copyright 2017 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. Laura++ package authors: John Back Paul Harrison Thomas Latham """ import sys # Process the command-line arguments def usage( progName ) : print('Usage:') - print('%s gen [nExpt = 1] [firstExpt = 0]' % progName) - print('%s fit [nExpt = 1] [firstExpt = 0]' % progName) + print(f'{progName} gen [nExpt = 1] [firstExpt = 0] [jsonFileName]') + print(f'{progName} fit [nExpt = 1] [firstExpt = 0] [jsonFileName]') if len(sys.argv) < 2 : usage( sys.argv[0] ) sys.exit(1) command = sys.argv[1] command = command.lower() iFit = 0 nExpt = 1 firstExpt = 0 + +jsonFileName = '@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/GenFit3pi.json' + if command == 'gen' : if len(sys.argv) > 2 : nExpt = int( sys.argv[2] ) if len(sys.argv) > 3 : firstExpt = int( sys.argv[3] ) + if len(sys.argv) > 4 : + jsonFileName = sys.argv[4] elif command == 'fit' : if len(sys.argv) < 3 : usage( sys.argv[0] ) sys.exit(1) iFit = int( sys.argv[2] ) if len(sys.argv) > 3 : nExpt = int( sys.argv[3] ) if len(sys.argv) > 4 : firstExpt = int( sys.argv[4] ) + if len(sys.argv) > 5 : + jsonFileName = sys.argv[5] else : usage( sys.argv[0] ) sys.exit(1) # Import ROOT and load the appropriate libraries import ROOT ROOT.gSystem.AddDynamicPath('@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_LIBDIR@') ROOT.gSystem.Load('libEG') ROOT.gSystem.Load('libHist') ROOT.gSystem.Load('libMatrix') ROOT.gSystem.Load('libTree') ROOT.gSystem.Load('libTreePlayer') ROOT.gSystem.Load('libMinuit') ROOT.gSystem.Load('libLaura++') # 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 True squareDP = False # This defines the DP => decay is B+ -> pi+ pi+ pi- # Particle 1 = pi+ # Particle 2 = pi+ # Particle 3 = pi- # The DP is defined in terms of m13Sq and m23Sq daughters = ROOT.LauDaughters('B+', 'pi+', 'pi+', 'pi-', squareDP) # Optionally apply some vetoes to the DP # (example syntax given but commented-out) vetoes = ROOT.LauVetoes() #DMin = 1.70 #DMax = 1.925 #JpsiMin = 3.051 #JpsiMax = 3.222 #psi2SMin = 3.676 #psi2SMax = 3.866 # D0 veto, m23 (and automatically m13) #vetoes.addMassVeto(1, DMin, DMax) # J/psi veto, m23 (and automatically m13) #vetoes.addMassVeto(1, JpsiMin, JpsiMax) # psi(2S) veto, m23 (and automatically m13) #vetoes.addMassVeto(1, psi2SMin, psi2SMax) # Define the efficiency model (defaults to unity everywhere) effModel = ROOT.LauEffModel(daughters, vetoes) # Can optionally provide a histogram to model variation over DP # (example syntax given but commented-out) #effHistFile = ROOT.TFile.Open('histoFiles/B3piNRDPEff.root', 'read') #effHist = effHistFile.Get('effHist') #useInterpolation = True #fluctuateBins = False #useUpperHalf = True #effModel.setEffHisto(effHist, useInterpolation, fluctuateBins, 0.0, 0.0, useUpperHalf, squareDP) -# Create the isobar model - -# Set the values of the Blatt-Weisskopf barrier radii and whether they are fixed or floating -resMaker = ROOT.LauResonanceMaker.get() -resMaker.setDefaultBWRadius( ROOT.LauBlattWeisskopfFactor.Category.Parent, 5.0 ) -resMaker.setDefaultBWRadius( ROOT.LauBlattWeisskopfFactor.Category.Light, 4.0 ) -resMaker.fixBWRadius( ROOT.LauBlattWeisskopfFactor.Category.Parent, True ) -resMaker.fixBWRadius( ROOT.LauBlattWeisskopfFactor.Category.Light, True ) - +# Create the isobar model, loading everything from a JSON file sigModel = ROOT.LauIsobarDynamics(daughters, effModel) - -# Add various components to the isobar model, -# modifying some resonance parameters -# addResonance arguments: resName, resPairAmpInt, resType -reson = sigModel.addResonance('rho0(770)', 1, ROOT.LauAbsResonance.GS) -reson = sigModel.addResonance('rho0(1450)', 1, ROOT.LauAbsResonance.RelBW) -reson = sigModel.addResonance('f_0(980)', 1, ROOT.LauAbsResonance.Flatte) -reson.setResonanceParameter('g1',0.2) -reson.setResonanceParameter('g2',1.0) -reson = sigModel.addResonance('f_2(1270)', 1, ROOT.LauAbsResonance.RelBW) -reson = sigModel.addResonance('BelleNR_Swave', 0, ROOT.LauAbsResonance.BelleSymNRNoInter) -reson.setResonanceParameter("alpha", 0.2) +sigModel.constructModelFromJson(jsonFileName, 'model') # Reset the maximum signal DP ASq value # This will be automatically adjusted to avoid bias or extreme # inefficiency if you get the value wrong but best to set this by # hand once you've found the right value through some trial and error. sigModel.setASqMaxValue(0.35) # Create the fit model fitModel = ROOT.LauSimpleFitModel(sigModel) # Create the complex coefficients for the isobar model -# Here we're using the magnitude and phase form: -# c_j = a_j exp(i*delta_j) -coeffset = [] -coeffset.append( ROOT.LauMagPhaseCoeffSet('rho0(770)', 1.00, 0.00, True, True) ) -coeffset.append( ROOT.LauMagPhaseCoeffSet("rho0(1450)", 0.37, 1.99, False, False) ) -coeffset.append( ROOT.LauMagPhaseCoeffSet("f_0(980)", 0.27, -1.59, False, False) ) -coeffset.append( ROOT.LauMagPhaseCoeffSet("f_2(1270)", 0.53, 1.39, False, False) ) -coeffset.append( ROOT.LauMagPhaseCoeffSet('BelleNR_Swave', 0.54, -0.84, False, False) ) -for i in coeffset : - fitModel.setAmpCoeffSet(i) +fitModel.setAmpCoeffSets( ROOT.LauAbsCoeffSet.readFromJson(jsonFileName, 'coeffs') ) # Set the signal yield and define whether it is fixed or floated nSigEvents = 1500.0 signalEvents = ROOT.LauParameter('signalEvents',nSigEvents,-1.0*nSigEvents,2.0*nSigEvents,False) fitModel.setNSigEvents(signalEvents) # Set the number of experiments to generate or fit and which experiment to start with fitModel.setNExpts( nExpt, firstExpt ) # Set up a background model # First declare the names of the background class(es) bkgndNames = ROOT.std.vector(ROOT.TString)() bkgndNames.push_back(ROOT.TString('qqbar')) fitModel.setBkgndClassNames( bkgndNames ) # Define and set the yield parameter for the background nBkg = 1250.0 nBkgndEvents = ROOT.LauParameter('qqbar',nBkg,-2.0*nBkg,2.0*nBkg,False) fitModel.setNBkgndEvents( nBkgndEvents ) # Create the background DP model qqbarModel = ROOT.LauBkgndDPModel(daughters, vetoes) # Load in background DP model histogram # (example syntax given but commented-out - the background will be treated as being uniform in the DP in the absence of a histogram) #qqFileName = 'histoFiles/offResDP.root' #qqFile = TFile.Open(qqFileName, 'read') #qqDP = qqFile.Get('AllmTheta') #qqbarModel.setBkgndHisto(qqDP, useInterpolation, fluctuateBins, useUpperHalf, squareDP) # Add the background DP model into the fit model fitModel.setBkgndDPModel( 'qqbar', qqbarModel ) # Configure various fit options # Switch on/off calculation of asymmetric errors. fitModel.useAsymmFitErrors(False) # Randomise initial fit values for the signal mode fitModel.useRandomInitFitPars(True) haveBkgnds = ( fitModel.nBkgndClasses() > 0 ) # Switch on/off Poissonian smearing of total number of events fitModel.doPoissonSmearing(haveBkgnds) # Switch on/off Extended ML Fit option fitModel.doEMLFit(haveBkgnds) # Generate toy from the fitted parameters -#fitToyFileName = 'fitToyMC_3pi_%d.root' % iFit +#fitToyFileName = f'fitToyMC_3pi_{iFit}.root' #fitModel.compareFitData(100, fitToyFileName) # Write out per-event likelihoods and sWeights -#splotFileName = 'splot_3pi_%d.root' % iFit +#splotFileName = f'splot_3pi_{iFit}.root' #fitModel.writeSPlotData(splotFileName, 'splot', False) # Set the names of the files to read/write dataFile = 'gen-3pi.root' treeName = 'genResults' rootFileName = '' tableFileName = '' if command == 'fit' : - rootFileName = 'fit3pi_%d_expt_%d-%d.root' % ( iFit, firstExpt, firstExpt+nExpt-1 ) - tableFileName = 'fit3piResults_%d' % iFit + rootFileName = f'fit3pi_{iFit}_expt_{firstExpt}-{firstExpt+nExpt-1}.root' + tableFileName = f'fit3piResults_{iFit}' else : rootFileName = 'dummy.root' tableFileName = 'gen3piResults' # Execute the generation/fit fitModel.run( command, dataFile, treeName, rootFileName, tableFileName ) +# Write out the model and coeffs to JSON for comparison with input +outputJsonFile = f'GenFit3pi-{command}.json' +sigModel.writeModelToJson( outputJsonFile, 'model' ) +ROOT.LauAbsCoeffSet.writeToJson( fitModel.getAmpCoeffs(), outputJsonFile, 'coeffs', True ) diff --git a/examples/GenFitBelleCPKpipi.cc b/examples/GenFitBelleCPKpipi.cc index 305b8ec..3341201 100644 --- a/examples/GenFitBelleCPKpipi.cc +++ b/examples/GenFitBelleCPKpipi.cc @@ -1,251 +1,257 @@ /* Copyright 2005 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ #include #include #include #include #include "TFile.h" #include "TH2.h" #include "TString.h" #include "TTree.h" #include "LauBkgndDPModel.hh" #include "LauCPFitModel.hh" #include "LauDaughters.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauBelleCPCoeffSet.hh" #include "LauVetoes.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; } // 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 = kFALSE; // These define the DP => decay is B+ -> K+ pi+ pi- and it's charge conjugate // The DP is defined in terms of m13Sq and m23Sq LauDaughters* negDaughters = new LauDaughters("B-", "K-", "pi-", "pi+", squareDP); LauDaughters* posDaughters = new LauDaughters("B+", "K+", "pi+", "pi-", squareDP); // Create the isobar models LauAbsResonance* res(0); using ResonanceModel = LauAbsResonance::ResonanceModel; LauIsobarDynamics* negSigModel = new LauIsobarDynamics(negDaughters, 0); res = negSigModel->addResonance("K*0(892)", 2, ResonanceModel::RelBW); // resPairAmpInt = 2 => resonance mass is m13. res = negSigModel->addResonance("K*0_0(1430)", 2, ResonanceModel::LASS_BW); res->floatResonanceParameter("a"); res->floatResonanceParameter("r"); res = negSigModel->addResonance("LASSNR0", 2, ResonanceModel::LASS_NR); res = negSigModel->addResonance("rho0(770)", 1, ResonanceModel::RelBW); // resPairAmpInt = 1 => resonance mass is m23. res = negSigModel->addResonance("f_0(980)", 1, ResonanceModel::Flatte); res = negSigModel->addResonance("chi_c0", 1, ResonanceModel::RelBW); res = negSigModel->addResonance("NonReson", 1, ResonanceModel::BelleNR); res->setResonanceParameter("alpha", 0.50); res->floatResonanceParameter("alpha"); LauIsobarDynamics* posSigModel = new LauIsobarDynamics(posDaughters, 0); res = posSigModel->addResonance("K*0(892)", 2, ResonanceModel::RelBW); res = posSigModel->addResonance("K*0_0(1430)", 2, ResonanceModel::LASS_BW); res->floatResonanceParameter("a"); res->floatResonanceParameter("r"); res = posSigModel->addResonance("LASSNR0", 2, ResonanceModel::LASS_NR); res = posSigModel->addResonance("rho0(770)", 1, ResonanceModel::RelBW); res = posSigModel->addResonance("f_0(980)", 1, ResonanceModel::Flatte); res = posSigModel->addResonance("chi_c0", 1, ResonanceModel::RelBW); res = posSigModel->addResonance("NonReson", 1, ResonanceModel::BelleNR); res->setResonanceParameter("alpha", 0.50); res->floatResonanceParameter("alpha"); // Set the file names for the integrals information (can be useful for debugging) negSigModel->setIntFileName("integ_neg.dat"); posSigModel->setIntFileName("integ_pos.dat"); // Reset the maximum signal DP ASq value // This will be automatically adjusted to avoid bias or extreme // inefficiency if you get the value wrong but best to set this by // hand once you've found the right value through some trial and // error. Double_t aSqMaxValue = 1.62; negSigModel->setASqMaxValue(aSqMaxValue); posSigModel->setASqMaxValue(aSqMaxValue); // Create the fit model, giving it both isobar models LauCPFitModel* fitModel = new LauCPFitModel( negSigModel, posSigModel ); // Create the complex coefficients for the isobar model // Here we're using the "Belle" form: // The amplitude has the form a * exp(i*delta) * ( 1 +/- b * exp(i*phi) ) where // a is a CP conserving magnitude, // b is a CP violating magnitude, // delta is the strong phase // and phi is the weak phase. // The last two arguments indicate that the CP-violating parameters // should only be floated in the second-stage of the two-stage fit auto a0CoeffSet = std::make_unique("K*0(892)", 1.00, 0.00, 0.00, 0.00, kTRUE, kTRUE, kFALSE, kTRUE, kTRUE, kTRUE); // Here we show how to use the coefficient cloning mechanism to clone // only the CP-violating parameters. This means that the K*0_0(1430) // component will have identical CPV parameters to the K*0(892). // But the CP-conserving parameters will be independent and as such we // need to set their values and, in this case because we're cloning // from the reference amplitude, also set them to float. //auto a1CoeffSet = std::make_unique("K*0_0(1430)", 2.00, 3.00, 0.00, 0.00, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE); auto a1CoeffSet = a0CoeffSet->createClone("K*0_0(1430)", LauAbsCoeffSet::CloneOption::TieCPPars); a1CoeffSet->setParameterValue("A",2.0,kTRUE); a1CoeffSet->setParameterValue("Delta",3.0,kTRUE); a1CoeffSet->floatParameter("A"); a1CoeffSet->floatParameter("Delta"); auto a2CoeffSet = a1CoeffSet->createClone("LASSNR0", LauAbsCoeffSet::CloneOption::TieCPPars); auto a3CoeffSet = std::make_unique("rho0(770)", 0.66, 1.00, 0.00, 0.00, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE); auto a4CoeffSet = std::make_unique("f_0(980)", 1.00, -1.00, 0.00, 0.00, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE); auto a5CoeffSet = std::make_unique("chi_c0", 0.33, 0.50, 0.00, 0.00, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE); auto a6CoeffSet = std::make_unique("NonReson", 0.50, 1.50, 0.00, 0.00, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE); std::vector> coeffset; coeffset.push_back( std::move(a0CoeffSet) ); coeffset.push_back( std::move(a1CoeffSet) ); coeffset.push_back( std::move(a2CoeffSet) ); coeffset.push_back( std::move(a3CoeffSet) ); coeffset.push_back( std::move(a4CoeffSet) ); coeffset.push_back( std::move(a5CoeffSet) ); coeffset.push_back( std::move(a6CoeffSet) ); for ( auto& coeff : coeffset ) { fitModel->setAmpCoeffSet( std::move(coeff) ); } // Set the signal yield and define whether it is fixed or floated const Double_t nSigEvents = 5000.0; Bool_t fixNSigEvents = kFALSE; LauParameter * signalEvents = new LauParameter("signalEvents",nSigEvents,-1.0*nSigEvents,2.0*nSigEvents,fixNSigEvents); fitModel->setNSigEvents(signalEvents); // Set the number of experiments to generate or fit and which // experiment to start with fitModel->setNExpts( nExpt, firstExpt ); // Configure various fit options // Switch on/off calculation of asymmetric errors. fitModel->useAsymmFitErrors(kFALSE); // Randomise initial fit values for the signal mode fitModel->useRandomInitFitPars(kTRUE); const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 ); // Switch on/off Poissonian smearing of total number of events fitModel->doPoissonSmearing(haveBkgnds); // Switch on/off Extended ML Fit option fitModel->doEMLFit(haveBkgnds); // Switch on/off two-stage fit fitModel->twoStageFit(kTRUE); // Generate toy from the fitted parameters //TString fitToyFileName("fitToyMC_BelleCPKpipi_"); //fitToyFileName += iFit; //fitToyFileName += ".root"; //fitModel->compareFitData(100, fitToyFileName); // Write out per-event likelihoods and sWeights //TString splotFileName("splot_BelleCPKpipi_"); //splotFileName += iFit; //splotFileName += ".root"; //fitModel->writeSPlotData(splotFileName, "splot", kFALSE); // Set the names of the files to read/write TString dataFile("gen-BelleCPKpipi.root"); TString treeName("genResults"); TString rootFileName(""); TString tableFileName(""); if (command == "fit") { rootFileName = "fitBelleCPKpipi_"; rootFileName += iFit; rootFileName += "_expt_"; rootFileName += firstExpt; rootFileName += "-"; rootFileName += (firstExpt+nExpt-1); rootFileName += ".root"; tableFileName = "fitBelleCPKpipiResults_"; tableFileName += iFit; } else { rootFileName = "dummy.root"; tableFileName = "genBelleCPKpipiResults"; } // Execute the generation/fit fitModel->run( command, dataFile, treeName, rootFileName, tableFileName ); + // Write out the DP models and coeffs to JSON + const TString jsonFileName {"GenFitBelleCPKpipi-"+command+"-model.json"}; + negSigModel->writeModelToJson( jsonFileName, "model_Bm", false ); + posSigModel->writeModelToJson( jsonFileName, "model_Bp", true ); + LauAbsCoeffSet::writeToJson( fitModel->getAmpCoeffs(), jsonFileName, "coeffs", true ); + return EXIT_SUCCESS; } diff --git a/examples/KMatrixDto3pi.cc b/examples/KMatrixDto3pi.cc.in similarity index 53% rename from examples/KMatrixDto3pi.cc rename to examples/KMatrixDto3pi.cc.in index 8aa03c6..8a43149 100644 --- a/examples/KMatrixDto3pi.cc +++ b/examples/KMatrixDto3pi.cc.in @@ -1,183 +1,158 @@ /* Copyright 2008 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ #include using std::cout; using std::cerr; using std::endl; #include #include #include "LauDaughters.hh" #include "LauVetoes.hh" #include "LauIsobarDynamics.hh" #include "LauSimpleFitModel.hh" #include "LauAbsCoeffSet.hh" #include "LauCartesianCPCoeffSet.hh" #include "LauEffModel.hh" #include "TString.h" void usage( std::ostream& out, const TString& progName ) { out<<"Usage:\n"; - out< [nExpt = 1] [firstExpt = 0]"< [nExpt = 1] [firstExpt = 0] [jsonFileName]"< [nExpt = 1] [firstExpt = 0] + // ./KMatrixDto3pi fit [nExpt = 1] [firstExpt = 0] [jsonFileName] 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); + + TString jsonFileName { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/KMatrixDto3pi.json" }; + if ( command == "gen" ) { if ( argc > 2 ) { nExpt = atoi( argv[2] ); if ( argc > 3 ) { firstExpt = atoi( argv[3] ); + if ( argc > 4 ) { + jsonFileName = argv[4]; + } } } } 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] ); + if ( argc > 5 ) { + jsonFileName = argv[5]; + } } } } else { usage( std::cerr, argv[0] ); return EXIT_FAILURE; } // Generate K-matrix amplitude distribution across the DP LauDaughters* daughters = new LauDaughters("D+", "pi+", "pi+", "pi-", kFALSE); LauVetoes* vetoes = new LauVetoes(); LauEffModel* effModel = new LauEffModel(daughters, vetoes); LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel); - - // Define the "S-wave" K-matrix propagator - Int_t nChannels = 5; - Int_t nPoles = 5; - Int_t resPairInt = 1; - Int_t KMatrixIndex = 1; // for S-wave - sigModel->defineKMatrixPropagator("KMSWave", "FOCUSD3pi.dat", resPairInt, nChannels, nPoles, KMatrixIndex); - // Now define the pole and "slowly-varying part" (SVP) amplitudes for the K-matrix. - // Each part will have their own complex coefficients which can (should) be fitted. - sigModel->addKMatrixProdPole("KMPole1", "KMSWave", 1); - sigModel->addKMatrixProdPole("KMPole2", "KMSWave", 2); - sigModel->addKMatrixProdPole("KMPole3", "KMSWave", 3); - //sigModel->addKMatrixProdPole("KMPole4", "KMSWave", 4); - //sigModel->addKMatrixProdPole("KMPole5", "KMSWave", 5); - - sigModel->addKMatrixProdSVP("KMSVP1", "KMSWave", 1); - sigModel->addKMatrixProdSVP("KMSVP2", "KMSWave", 2); - //sigModel->addKMatrixProdSVP("KMSVP3", "KMSWave", 3); - //sigModel->addKMatrixProdSVP("KMSVP4", "KMSWave", 4); - //sigModel->addKMatrixProdSVP("KMSVP5", "KMSWave", 5); - - sigModel->addResonance("rho0(770)", 1, LauAbsResonance::ResonanceModel::RelBW); - sigModel->addResonance("f_2(1270)", 1, LauAbsResonance::ResonanceModel::RelBW); - + sigModel->constructModelFromJson( jsonFileName, "model" ); sigModel->setASqMaxValue(2.1); LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel); - - std::vector> coeffset; - - coeffset.push_back(std::make_unique("KMPole1", 0.85, -0.55, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - coeffset.push_back(std::make_unique("KMPole2", -1.10, -0.30, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - coeffset.push_back(std::make_unique("KMPole3", -0.84, 0.35, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - //coeffset.push_back(std::make_unique("KMPole4", -0.004, .46, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - //coeffset.push_back(std::make_unique("KMPole5", 0.00, 0.00, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - - coeffset.push_back(std::make_unique("KMSVP1", 0.230, 0.352, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - coeffset.push_back(std::make_unique("KMSVP2", 0.330, 0.42, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - //coeffset.push_back(std::make_unique("KMSVP3", -0.48, 1.25, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - //coeffset.push_back(std::make_unique("KMSVP4", 0.30, 1.32, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - //coeffset.push_back(std::make_unique("KMSVP5", 0.00, 0.00, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - - coeffset.push_back(std::make_unique("rho0(770)", 1.0, 0.00, 0.0, 0.0, kTRUE, kTRUE, kTRUE, kTRUE)); - coeffset.push_back(std::make_unique("f_2(1270)", 0.45, 0.35, 0.0, 0.0, kFALSE, kFALSE, kTRUE, kTRUE)); - - for ( auto& coeff : coeffset ) { - fitModel->setAmpCoeffSet( std::move(coeff) ); - } + auto coeffs { LauAbsCoeffSet::readFromJson( jsonFileName, "coeffs" ) }; + fitModel->setAmpCoeffSets( coeffs ); Double_t nSigEvents = 20000; Bool_t fixSigEvents = kTRUE; LauParameter * nSig = new LauParameter("signalEvents",nSigEvents,-2.0*nSigEvents,2.0*nSigEvents,fixSigEvents); fitModel->setNSigEvents(nSig); fitModel->setNExpts(nExpt, firstExpt); // Do not calculate asymmetric errors fitModel->useAsymmFitErrors(kFALSE); // Randomise initial fit values for the signal mode fitModel->useRandomInitFitPars(kTRUE); // Switch off Poissonian smearing of total number of events fitModel->doPoissonSmearing(kFALSE); // Switch on Extended ML Fit option - only really use if Poisson smearing is also enabled fitModel->doEMLFit(kFALSE); // Set the names of the files to read/write TString dataFile("gen-KMatrixDto3pi.root"); TString treeName("genResults"); TString rootFileName(""); TString tableFileName(""); if (command == "fit") { rootFileName = "fitKMatrixDto3pi_"; rootFileName += iFit; rootFileName += "_expt_"; rootFileName += firstExpt; rootFileName += "-"; rootFileName += (firstExpt+nExpt-1); rootFileName += ".root"; tableFileName = "fitKMatrixDto3piResults_"; tableFileName += iFit; } else { rootFileName = "dummy.root"; tableFileName = "genKMatrixDto3piResults"; } // Execute the generation/fit fitModel->run( command, dataFile, treeName, rootFileName, tableFileName ); + + // Write out the DP models and coeffs to JSON (for comparison with input) + const TString outputJsonFileName {"KMatrixDto3pi-"+command+".json"}; + sigModel->writeModelToJson( outputJsonFileName, "model", false ); + LauAbsCoeffSet::writeToJson( fitModel->getAmpCoeffs(), outputJsonFileName, "coeffs", true ); + + return EXIT_SUCCESS; } diff --git a/examples/KMatrixDto3pi.json.in b/examples/KMatrixDto3pi.json.in new file mode 100644 index 0000000..bda8b80 --- /dev/null +++ b/examples/KMatrixDto3pi.json.in @@ -0,0 +1,246 @@ +{ + "coeffs": { + "coeffs": [ + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": 1.0, + "XBlind": false, + "XFixed": true, + "XSecondStage": false, + "Y": 0.0, + "YBlind": false, + "YFixed": true, + "YSecondStage": false, + "clone": false, + "name": "rho0(770)", + "type": "CartesianCP" + }, + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": 0.45, + "XBlind": false, + "XFixed": false, + "XSecondStage": false, + "Y": 0.35, + "YBlind": false, + "YFixed": false, + "YSecondStage": false, + "clone": false, + "name": "f_2(1270)", + "type": "CartesianCP" + }, + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": 0.85, + "XBlind": false, + "XFixed": false, + "XSecondStage": false, + "Y": -0.55, + "YBlind": false, + "YFixed": false, + "YSecondStage": false, + "clone": false, + "name": "KMPole1", + "type": "CartesianCP" + }, + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": -1.1, + "XBlind": false, + "XFixed": false, + "XSecondStage": false, + "Y": -0.3, + "YBlind": false, + "YFixed": false, + "YSecondStage": false, + "clone": false, + "name": "KMPole2", + "type": "CartesianCP" + }, + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": -0.84, + "XBlind": false, + "XFixed": false, + "XSecondStage": false, + "Y": 0.35, + "YBlind": false, + "YFixed": false, + "YSecondStage": false, + "clone": false, + "name": "KMPole3", + "type": "CartesianCP" + }, + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": 0.23, + "XBlind": false, + "XFixed": false, + "XSecondStage": false, + "Y": 0.352, + "YBlind": false, + "YFixed": false, + "YSecondStage": false, + "clone": false, + "name": "KMSVP1", + "type": "CartesianCP" + }, + { + "DeltaX": 0.0, + "DeltaXBlind": false, + "DeltaXFixed": true, + "DeltaXSecondStage": false, + "DeltaY": 0.0, + "DeltaYBlind": false, + "DeltaYFixed": true, + "DeltaYSecondStage": false, + "X": 0.33, + "XBlind": false, + "XFixed": false, + "XSecondStage": false, + "Y": 0.42, + "YBlind": false, + "YFixed": false, + "YSecondStage": false, + "clone": false, + "name": "KMSVP2", + "type": "CartesianCP" + } + ], + "nCoeffs": 7 + }, + "model": { + "commonSettings": { + "setBWBachelorRestFrame": "ResonanceFrame", + "setBWType": "BWPrimeBarrier", + "setDefaultBWRadius": [ + { + "category": "Parent", + "fix": true, + "value": 4.0 + }, + { + "category": "Light", + "fix": true, + "value": 5.3 + } + ], + "setSpinFormalism": "Zemach_P" + }, + "kmatrix": [ + { + "poles": [ + { + "poleIndex": 1, + "poleName": "KMPole1" + }, + { + "poleIndex": 2, + "poleName": "KMPole2" + }, + { + "poleIndex": 3, + "poleName": "KMPole3" + } + ], + "propagator": { + "nChannels": 5, + "nPoles": 5, + "paramFileName": "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/FOCUSD3pi.dat", + "propName": "KMSWave", + "resPairAmpInt": 1 + }, + "svps": [ + { + "channelIndex": 1, + "svpName": "KMSVP1" + }, + { + "channelIndex": 2, + "svpName": "KMSVP2" + } + ] + } + ], + "resonances": [ + { + "changeResonance": { + "mass": { + "fix": true, + "value": 0.77526 + }, + "spin": { + "value": 1 + }, + "width": { + "fix": true, + "value": 0.1478 + } + }, + "resName": "rho0(770)", + "resPairAmpInt": 1, + "resType": "RelBW" + }, + { + "changeResonance": { + "mass": { + "fix": true, + "value": 1.2755 + }, + "spin": { + "value": 2 + }, + "width": { + "fix": true, + "value": 0.1867 + } + }, + "resName": "f_2(1270)", + "resPairAmpInt": 1, + "resType": "RelBW" + } + ] + } +}