diff --git a/README.md b/README.md index 904ae6d..2bf2383 100644 --- a/README.md +++ b/README.md @@ -1,474 +1,475 @@ # Laura++ README ## Introduction This package is a C++ development of LAURA, the FORTRAN "Likelihood Analysis Unofficial RooFitDalitz Alternative" code written by Paul Harrison for performing Dalitz plot analyses of 3-body decays of B mesons at BaBar. This C++ version was originally developed by John Back and is now mainly developed by Tom Latham with continuting contributions from John Back and others. ## Structure of the Package The package consists of the following directories and files: ``` laura.git ├── CMakeLists.txt - the top-level CMake configuration file ├── LICENSE-2.0 - the Apache License, Version 2.0 ├── README.md - this README file ├── cmake - directory containing custom CMake modules and templates │   ├── Modules │   │   ├── LauraCompilerFlags.cmake │   │   └── LauraExternalDependencies.cmake │   └── Templates │   └── LauraConfig.cmake.in ├── doc - directory related to building Doxygen documentation │   ├── CMakeLists.txt - CMake configuration file │   ├── Doxyfile.in - Doxygen configuration file │   ├── mainpage.dox.in - front page of the Doxygen documents │   └── ReleaseNotes.md - history of commits to the package ├── examples - directory containing code examples and their associated data files and job scripts │   ├── B3piKMNoAdler.dat │   ├── B3piKMPoles.dat │   ├── B3piKMSVPs.dat │   ├── B3piKMatrixCoeff.dat │   ├── B3piKMatrixMassProj.cc │   ├── B3piKMatrixPlots.cc │   ├── CMakeLists.txt │   ├── CalcChiSq.cc │   ├── DToKspipiKMatrixCoeff.dat │   ├── FOCUSD3pi.dat │   ├── GenFit3K.cc │   ├── GenFit3KS.cc │   ├── GenFit3pi.cc │   ├── GenFit3pi.py │   ├── GenFitBelleCPKpipi.cc │   ├── GenFitDs2KKpi.cc │   ├── GenFitEFKLLM.cc │   ├── GenFitKpipi.cc │   ├── GenFitNoDP.cc │   ├── GenFitNoDPMultiDim.cc │   ├── KMatrixDto3pi.cc │   ├── KMatrixExample.cc │   ├── LauKMatrixCoeff.dat │   ├── LauKMatrixCoeff2.dat │   ├── MergeDataFiles.cc │   ├── PlotKMatrixTAmp.cc │   ├── PlotResults.cc │   ├── ResultsExtractor.cc │   ├── SimFitCoordinator.cc │   ├── SimFitTask.cc │   ├── SimFitTaskRooFit.cc │   ├── chiSqInput.txt │   ├── mixedSampleTest.cc │   ├── point2PointTestSample.cc │   ├── runCoordinatorRooFitTask.sh │   ├── runCoordinatorTask.sh │   ├── runPoint2PointTest.sh │   └── usfactor.dat ├── inc - directory containing the header files for the Laura++ library │   ├── CMakeLists.txt │   ├── Lau1DCubicSpline.hh │   ├── Lau1DHistPdf.hh │   ├── Lau2DAbsDP.hh │   ├── Lau2DAbsDPPdf.hh │   ├── Lau2DAbsHistDP.hh │   ├── Lau2DAbsHistDPPdf.hh │   ├── Lau2DCubicSpline.hh │   ├── Lau2DHistDP.hh │   ├── Lau2DHistDPPdf.hh │   ├── Lau2DHistPdf.hh │   ├── Lau2DSplineDP.hh │   ├── Lau2DSplineDPPdf.hh │   ├── LauAbsBkgndDPModel.hh │   ├── LauAbsCoeffSet.hh │   ├── LauAbsEffModel.hh │   ├── LauAbsFitModel.hh │   ├── LauAbsFitter.hh │   ├── LauAbsIncohRes.hh │   ├── LauAbsModIndPartWave.hh │   ├── LauAbsPdf.hh │   ├── LauAbsRValue.hh │   ├── LauAbsResonance.hh │   ├── LauArgusPdf.hh │   ├── LauAsymmCalc.hh │   ├── LauBelleCPCoeffSet.hh │   ├── LauBelleNR.hh │   ├── LauBelleSymNR.hh │   ├── LauBifurcatedGaussPdf.hh │   ├── LauBkgndDPModel.hh │   ├── LauBlattWeisskopfFactor.hh │   ├── LauBlind.hh │   ├── LauBreitWignerRes.hh │   ├── LauCPFitModel.hh │   ├── LauCacheData.hh │   ├── LauCalcChiSq.hh │   ├── LauCartesianCPCoeffSet.hh │   ├── LauCartesianGammaCPCoeffSet.hh │   ├── LauChebychevPdf.hh │   ├── LauCleoCPCoeffSet.hh │   ├── LauComplex.hh │   ├── LauConstants.hh │   ├── LauCruijffPdf.hh │   ├── LauCrystalBallPdf.hh │   ├── LauDPDepBifurGaussPdf.hh │   ├── LauDPDepCruijffPdf.hh │   ├── LauDPDepGaussPdf.hh │   ├── LauDPDepMapPdf.hh │   ├── LauDPDepSumPdf.hh │   ├── LauDPPartialIntegralInfo.hh │   ├── LauDabbaRes.hh │   ├── LauDatabasePDG.hh │   ├── LauDaughters.hh │   ├── LauEFKLLMRes.hh │   ├── LauEffModel.hh │   ├── LauEmbeddedData.hh │   ├── LauExponentialPdf.hh │   ├── LauFitDataTree.hh │   ├── LauFitNtuple.hh │   ├── LauFitObject.hh │   ├── LauFitter.hh │   ├── LauFlatNR.hh │   ├── LauFlatteRes.hh │   ├── LauFormulaPar.hh │   ├── LauGaussIncohRes.hh │   ├── LauGaussPdf.hh │   ├── LauGenNtuple.hh │   ├── LauGounarisSakuraiRes.hh │   ├── LauIntegrals.hh │   ├── LauIsobarDynamics.hh │   ├── LauKMatrixProdPole.hh │   ├── LauKMatrixProdSVP.hh │   ├── LauKMatrixPropFactory.hh │   ├── LauKMatrixPropagator.hh │   ├── LauKappaRes.hh │   ├── LauKinematics.hh │   ├── LauLASSBWRes.hh │   ├── LauLASSNRRes.hh │   ├── LauLASSRes.hh │   ├── LauLinearPdf.hh │   ├── LauMagPhaseCPCoeffSet.hh │   ├── LauMagPhaseCoeffSet.hh │   ├── LauMergeDataFiles.hh │   ├── LauMinuit.hh │   ├── LauModIndPartWaveMagPhase.hh │   ├── LauModIndPartWaveRealImag.hh │   ├── LauNRAmplitude.hh │   ├── LauNovosibirskPdf.hh │   ├── LauParamFixed.hh │   ├── LauParameter.hh │   ├── LauParametricStepFuncPdf.hh │   ├── LauParticlePDG.hh │   ├── LauPolNR.hh │   ├── LauPolarFormFactorNR.hh │   ├── LauPolarFormFactorSymNR.hh │   ├── LauPolarGammaCPCoeffSet.hh │   ├── LauPoleRes.hh │   ├── LauPrint.hh │   ├── LauRandom.hh │   ├── LauRealImagCPCoeffSet.hh │   ├── LauRealImagCoeffSet.hh │   ├── LauRealImagGammaCPCoeffSet.hh │   ├── LauRelBreitWignerRes.hh │   ├── LauRescattering2Res.hh │   ├── LauRescatteringRes.hh │   ├── LauResonanceInfo.hh │   ├── LauResonanceMaker.hh │   ├── LauResultsExtractor.hh │   ├── LauRhoOmegaMix.hh │   ├── LauRooFitTask.hh │   ├── LauSPlot.hh │   ├── LauScfMap.hh │   ├── LauSigmaRes.hh │   ├── LauSigmoidPdf.hh │   ├── LauSimFitCoordinator.hh │   ├── LauSimFitTask.hh │   ├── LauSimpleFitModel.hh │   ├── LauString.hh │   ├── LauSumPdf.hh │   ├── LauTextFileParser.hh │   ├── LauVetoes.hh │   ├── LauWeightedSumEffModel.hh │   └── Laura++_LinkDef.h ├── macros - directory containing some ROOT macro files for plotting results, etc. │   ├── makeFitFile.C │   ├── plotComplexVars.C │   ├── plotCorrs.C │   ├── plotCorrsFromToy.C │   └── plotDataIsobars.C ├── src - directory containing the source files for the Laura++ library │   ├── CMakeLists.txt │   ├── Lau1DCubicSpline.cc │   ├── Lau1DHistPdf.cc │   ├── Lau2DAbsHistDP.cc │   ├── Lau2DAbsHistDPPdf.cc │   ├── Lau2DCubicSpline.cc │   ├── Lau2DHistDP.cc │   ├── Lau2DHistDPPdf.cc │   ├── Lau2DHistPdf.cc │   ├── Lau2DSplineDP.cc │   ├── Lau2DSplineDPPdf.cc │   ├── LauAbsBkgndDPModel.cc │   ├── LauAbsCoeffSet.cc │   ├── LauAbsFitModel.cc │   ├── LauAbsIncohRes.cc │   ├── LauAbsModIndPartWave.cc │   ├── LauAbsPdf.cc │   ├── LauAbsResonance.cc │   ├── LauArgusPdf.cc │   ├── LauAsymmCalc.cc │   ├── LauBelleCPCoeffSet.cc │   ├── LauBelleNR.cc │   ├── LauBelleSymNR.cc │   ├── LauBifurcatedGaussPdf.cc │   ├── LauBkgndDPModel.cc │   ├── LauBlattWeisskopfFactor.cc │   ├── LauBlind.cc │   ├── LauBreitWignerRes.cc │   ├── LauCPFitModel.cc │   ├── LauCacheData.cc │   ├── LauCalcChiSq.cc │   ├── LauCartesianCPCoeffSet.cc │   ├── LauCartesianGammaCPCoeffSet.cc │   ├── LauChebychevPdf.cc │   ├── LauCleoCPCoeffSet.cc │   ├── LauComplex.cc │   ├── LauCruijffPdf.cc │   ├── LauCrystalBallPdf.cc │   ├── LauDPDepBifurGaussPdf.cc │   ├── LauDPDepCruijffPdf.cc │   ├── LauDPDepGaussPdf.cc │   ├── LauDPDepMapPdf.cc │   ├── LauDPDepSumPdf.cc │   ├── LauDPPartialIntegralInfo.cc │   ├── LauDabbaRes.cc │   ├── LauDatabasePDG.cc │   ├── LauDaughters.cc │   ├── LauEFKLLMRes.cc │   ├── LauEffModel.cc │   ├── LauEmbeddedData.cc │   ├── LauExponentialPdf.cc │   ├── LauFitDataTree.cc │   ├── LauFitNtuple.cc │   ├── LauFitObject.cc │   ├── LauFitter.cc │   ├── LauFlatNR.cc │   ├── LauFlatteRes.cc │   ├── LauFormulaPar.cc │   ├── LauGaussIncohRes.cc │   ├── LauGaussPdf.cc │   ├── LauGenNtuple.cc │   ├── LauGounarisSakuraiRes.cc │   ├── LauIntegrals.cc │   ├── LauIsobarDynamics.cc │   ├── LauKMatrixProdPole.cc │   ├── LauKMatrixProdSVP.cc │   ├── LauKMatrixPropFactory.cc │   ├── LauKMatrixPropagator.cc │   ├── LauKappaRes.cc │   ├── LauKinematics.cc │   ├── LauLASSBWRes.cc │   ├── LauLASSNRRes.cc │   ├── LauLASSRes.cc │   ├── LauLinearPdf.cc │   ├── LauMagPhaseCPCoeffSet.cc │   ├── LauMagPhaseCoeffSet.cc │   ├── LauMergeDataFiles.cc │   ├── LauMinuit.cc │   ├── LauModIndPartWaveMagPhase.cc │   ├── LauModIndPartWaveRealImag.cc │   ├── LauNRAmplitude.cc │   ├── LauNovosibirskPdf.cc │   ├── LauParameter.cc │   ├── LauParametricStepFuncPdf.cc │   ├── LauParticlePDG.cc │   ├── LauPolNR.cc │   ├── LauPolarFormFactorNR.cc │   ├── LauPolarFormFactorSymNR.cc │   ├── LauPolarGammaCPCoeffSet.cc │   ├── LauPoleRes.cc │   ├── LauPrint.cc │   ├── LauRandom.cc │   ├── LauRealImagCPCoeffSet.cc │   ├── LauRealImagCoeffSet.cc │   ├── LauRealImagGammaCPCoeffSet.cc │   ├── LauRelBreitWignerRes.cc │   ├── LauRescattering2Res.cc │   ├── LauRescatteringRes.cc │   ├── LauResonanceInfo.cc │   ├── LauResonanceMaker.cc │   ├── LauResultsExtractor.cc │   ├── LauRhoOmegaMix.cc │   ├── LauRooFitTask.cc │   ├── LauSPlot.cc │   ├── LauScfMap.cc │   ├── LauSigmaRes.cc │   ├── LauSigmoidPdf.cc │   ├── LauSimFitCoordinator.cc │   ├── LauSimFitTask.cc │   ├── LauSimpleFitModel.cc │   ├── LauString.cc │   ├── LauSumPdf.cc │   ├── LauTextFileParser.cc │   ├── LauVetoes.cc │   └── LauWeightedSumEffModel.cc └── test - directory containing code for some test executables ├── CMakeLists.txt ├── TestCovariant.cc ├── TestCovariant2.cc └── TestNewKinematicsMethods.cc ``` ## Building the library ### Prerequisites Compilation of Laura++ on a UNIX operating system requires: * the [CMake](https://www.cmake.org) build tool * a C++14 compatible compiler ([GCC](https://gcc.gnu.org) 8 or better, [Clang](https://clang.llvm.org) 8 or better are recommended) * [GNU Make](https://www.gnu.org/software/make/) or [Ninja](https://ninja-build.org) (not currently tested) The package depends only on [ROOT](https://root.cern.ch). Before building the code, it is necessary that the ROOT package be findable by CMake, which can be achieved by doing one of the following: - the `ROOTSYS` environment variable is set to the directory of the ROOT package - the directory containing the `root-config` program is in the `PATH` environment variable - the ROOT cmake directory is in the `CMAKE_PREFIX_PATH` environment variable In order to setup an environment that satifies all of the above requirements we highly recommend use of the [LCG views](https://lcginfo.cern.ch). For example, on a machine with the [CVMFS client](https://cvmfs.readthedocs.io/en/stable/) installed, one can do: ``` source /cvmfs/sft.cern.ch/lcg/views/setupViews.(c)sh ``` for example: ``` source /cvmfs/sft.cern.ch/lcg/views/setupViews.sh LCG_96b x86_64-centos7-gcc9-opt ``` ### Build procedure To build from a clone of this repository, open a terminal window and change directory into that holding this README file. Create a build directory in which to run `cmake` and the build, and change into it by doing the following: ``` $ mkdir ../laura.build $ cd ../laura.build ``` It is best to install the package after building. You can either install into an existing location or you can create an install directory specifically for Laura++ (we assume the latter in the following instructions): ``` $ mkdir ../laura.install ``` Run `cmake` in this directory, pointing it to the directory holding this README, and consequently the top level CMake script for the project: ``` $ cmake ../ -DCMAKE_INSTALL_PREFIX=../laura.install ... system specific output ... -- Configuring done -- Generating done -- Build files have been written to: ... your build dir path ... $ ``` The exact output will depend on your system, compiler and build directory location, but you should not see any errors. CMake will generate Makefiles by default on UNIX platforms, so to build, simply run `make` in the build directory: ``` $ make ... verbose output ... [100%] Built target Laura++ ... $ make install Install the project... ... verbose output ... $ ``` Again, the exact output will be system specific, but you should see the `Laura++` target built without error and that the installation was successful. If compiling with gcc 5, see [the corresponding note](#notes). ## Examples and documentation Example code is included in the examples directory. To enable the building and installation of the example executables you need to supply the following option when running `cmake`: ``` -DLAURA_BUILD_EXAMPLES=ON ``` After building and installing, the example executables will be in the `bin` directory in the install location. To build the Doxygen documentation you need to supply the following option when running `cmake`: ``` -DLAURA_BUILD_DOCS=ON ``` After building and installing, you can load the `share/doc/Laura++/html/index.html` file (found in the install location) into your web browser. ## Building your own code If you wish to build your own code that uses the Laura++ library you should add the Laura++ installation area to the `CMAKE_PREFIX_PATH` environment variable. You can then add lines like the following to your own `CMakeLists.txt` file: ```cmake # Find ROOT and Laura++ find_package(Laura REQUIRED) find_package(ROOT REQUIRED) # Now build the executable add_executable(MyExe MyExe.cc) target_link_libraries(MyExe PRIVATE Laura::Laura++) ``` ## Authors and contributors As mentioned in the Introduction above, the authors of the package are: Thomas Latham John Back Paul Harrison The authors would also like to thank the following people for their invaluable contributions: Sian Morgan Tim Gershon Pablo del Amo Sanchez Jelena Ilic Eugenia Puccio Mark Whitehead Daniel Craik Rafael Coutinho Charlotte Wallace Juan Otalora Wenbin Qian Daniel O'Hanlon +Daniel Johnson Many thanks also go to numerous members of the BaBar, Belle, and LHCb collaborations for their helpful input. Contact: Tom Latham (laura@projects.hepforge.org) ## Notes {#notes} * [It has been reported](https://laura.hepforge.org/trac/ticket/67) that when compiling on Ubuntu 16.04 LTS with gcc 5.4.0 the build will fail due to, apparently spurious, "variable defined but not used" errors related to the constants defined in the `LauConstants.hh` header file. The workaround is to add the option `-Wno-error=unused-variable` to the `CMAKE_CXX_FLAGS` variable in the file: cmake/Modules/LauraCompilerFlags.cmake diff --git a/inc/LauBlattWeisskopfFactor.hh b/inc/LauBlattWeisskopfFactor.hh index f5a71ce..1f9d325 100644 --- a/inc/LauBlattWeisskopfFactor.hh +++ b/inc/LauBlattWeisskopfFactor.hh @@ -1,145 +1,146 @@ /* Copyright 2014 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 */ /*! \file LauBlattWeisskopfFactor.hh \brief File containing declaration of LauBlattWeisskopfFactor class. */ /*! \class LauBlattWeisskopfFactor \brief Class that implements the Blatt-Weisskopf barrier factor This factor is used in resonance models such as the Relativistic Breit-Wigner. */ #ifndef LAU_BLATTWEISSKOPFFACTOR #define LAU_BLATTWEISSKOPFFACTOR #include "Rtypes.h" #include "TString.h" class LauParameter; class LauResonanceInfo; class LauBlattWeisskopfFactor final { public: //! Define the allowed types of barrier factors enum BarrierType { BWBarrier, /*!< Blatt-Weisskopf barrier factor (for use when momentum terms not used in angular term) */ BWPrimeBarrier, /*!< Blatt-Weisskopf barrier factor (for use when momentum terms are used in angular term) - the default */ ExpBarrier /*!< expoential barrier factor (mostly used for virtual contributions) */ }; //! Define the rest frame in which the momentum should be calculated (only relevant for bachelor) enum RestFrame { Covariant, /*!< use expression from Covariant spin factor */ ParentFrame, /*!< momentum calculated in parent rest frame */ ResonanceFrame /*!< momentum calculated in resonance rest frame */ }; //! Define resonance categories that will share common barrier factor radii enum BlattWeisskopfCategory { Default, //*!< indicates that LauResonanceMaker should use the appropriate category for the given resonance */ Parent, //*!< indicates that this is the factor for the decay of the parent particle */ Indep, //*!< indicates that this resonance should be independent of all others */ Light, //*!< default category for light unflavoured states, e.g. rho(77), f0(980), etc. */ Kstar, //*!< default category for K* states */ Charm, //*!< default category for D* states */ StrangeCharm, //*!< default category for Ds* states */ Charmonium, //*!< default category for ccbar states */ Beauty, //*!< default category for B* states */ StrangeBeauty, //*!< default category for Bs* states */ CharmBeauty, //*!< default category for Bc* states */ Custom1, //*!< user-defined category */ Custom2, //*!< user-defined category */ Custom3, //*!< user-defined category */ Custom4 //*!< user-defined category */ }; //! Constructor LauBlattWeisskopfFactor( const LauResonanceInfo& resInfo, const BarrierType barrierType, const RestFrame restFrame, const BlattWeisskopfCategory category ); //! Constructor LauBlattWeisskopfFactor( const LauResonanceInfo& resInfo, const Double_t resRadius, const BarrierType barrierType, const RestFrame restFrame, const BlattWeisskopfCategory category ); //! Constructor LauBlattWeisskopfFactor( const Int_t spin, const Double_t resRadius, const BarrierType barrierType, const RestFrame restFrame, const BlattWeisskopfCategory category ); //! Method to create a new factor with cloned radius parameter /*! \param newSpin the value of the spin to use for the created instance + \param newBarrierType the BarrierType to use for the created instance */ LauBlattWeisskopfFactor* createClone( const UInt_t newSpin , const BarrierType newBarrierType ); //! Retrieve the radius parameter const LauParameter* getRadiusParameter() const { return radius_; } //! Retrieve the radius parameter LauParameter* getRadiusParameter() { return radius_; } //! Retrieve the barrier type BarrierType getBarrierType() const { return barrierType_; } //! Retrieve the rest frame information RestFrame getRestFrame() const { return restFrame_; } //! Calculate form factor value /*! \param p the value of the momentum */ Double_t calcFormFactor( const Double_t p ) const; private: //! Copy constructor LauBlattWeisskopfFactor( const LauBlattWeisskopfFactor& other, const UInt_t newSpin, const BarrierType newBarrierType ); //! Copy assignment operator (not implemented) LauBlattWeisskopfFactor& operator=( const LauBlattWeisskopfFactor& other ); //! Set the name of the radius parameter TString setRadiusName( const LauResonanceInfo& resInfo, const BlattWeisskopfCategory category ); //! Set the name of the radius parameter TString setRadiusName( const BlattWeisskopfCategory category ); //! Set the name of the radius parameter TString setRadiusName( const TString& categoryName ); //! Resonance spin const UInt_t spin_; //! Radius parameter LauParameter* radius_; //! Barrier type const BarrierType barrierType_; //! Rest frame const RestFrame restFrame_; ClassDef(LauBlattWeisskopfFactor, 0) }; #endif diff --git a/inc/LauKMatrixPropagator.hh b/inc/LauKMatrixPropagator.hh index 46619a5..3c11acd 100644 --- a/inc/LauKMatrixPropagator.hh +++ b/inc/LauKMatrixPropagator.hh @@ -1,637 +1,637 @@ /* 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 */ /*! \file LauKMatrixPropagator.hh \brief File containing declaration of LauKMatrixPropagator class. */ /*! \class LauKMatrixPropagator \brief Class for defining a K-matrix propagator. Class used to define a K-matrix propagator. See the following papers for info: hep-ph/0204328, hep-ex/0312040, [hep-ex]0804.2089 and hep-ph/9705401. */ #ifndef LAU_KMATRIX_PROPAGATOR #define LAU_KMATRIX_PROPAGATOR #include "LauConstants.hh" #include "LauResonanceMaker.hh" #include "LauResonanceInfo.hh" #include "TMatrixD.h" #include "TString.h" #include #include class LauParameter; class LauKinematics; class LauComplex; class LauKMatrixPropagator { public: //! Constructor /*! \param [in] name name of the propagator \param [in] paramFileName the parameter file name \param [in] resPairAmpInt the number of the daughter not produced by the resonance \param [in] nChannels the number of channels \param [in] nPoles the number of poles \param [in] rowIndex this specifies which row of the propagator should be used when summing over the amplitude channels */ LauKMatrixPropagator( const TString& name, const TString& paramFileName, const Int_t resPairAmpInt, const Int_t nChannels, const Int_t nPoles, const Int_t rowIndex = 1 ); //! Destructor virtual ~LauKMatrixPropagator(); //! Calculate the K-matrix propagator for the given s value /*! \param [in] s the invariant mass squared */ void updatePropagator(const Double_t s); //! Read an input file to set parameters /*! \param [in] inputFile name of the input file */ void setParameters(const TString& inputFile); //! Set flag to ignore Blatt-Weisskopf-like barrier factor void ignoreBWBarrierFactor() {includeBWBarrierFactor_=kFALSE;} //! Get the scattering K matrix /*! \return the real, symmetric scattering K matrix */ TMatrixD getKMatrix() const {return ScattKMatrix_;} //! Get the real part of the propagator full matrix /*! \return the real part of the propagator full matrix */ TMatrixD getRealPropMatrix() const {return realProp_;} //! Get the negative imaginary part of the full propagator matrix /*! \return the negative imaginary part of the full propagator matrix */ TMatrixD getNegImagPropMatrix() const {return negImagProp_;} //! Get the real part of the term of the propagator /*! \param [in] channelIndex the channel number \return the real part of the propagator term */ Double_t getRealPropTerm(const Int_t channelIndex) const; //! Get the imaginary part of the term of the propagator /*! \param [in] channelIndex the channel number \return the imaginiary part of the propagator term */ Double_t getImagPropTerm(const Int_t channelIndex) const; //! Get the 1/(m_pole^2 -s) terms for the scattering and production K-matrix formulae /*! \param [in] poleIndex the number of the pole required \return the value of 1/(m_pole^2 -s) */ Double_t getPoleDenomTerm(const Int_t poleIndex) const; //! Get spin of K-matrix /*! \param [in] iChannel the index of the channel \return the value of the orbital angular momentum, L_, for this channel */ Int_t getL(const Int_t iChannel) const {return L_[iChannel];} //! Get index of final channel /*! \return the index of the channel into which the scattering happens */ Int_t getIndex() const {return index_;}; //! Get pole mass parameters, set according to the input file /*! \param [in] poleIndex number of the required pole - \param [in] channelIndex number of the required channel \return the parameter of the pole mass */ LauParameter& getPoleMassSqParameter(const Int_t poleIndex); //! Get coupling constants that were loaded from the input file /*! \param [in] poleIndex number of the required pole \param [in] channelIndex number of the required channel \return the value of the coupling constant */ Double_t getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const; //! Get coupling parameters, set according to the input file /*! \param [in] poleIndex number of the required pole \param [in] channelIndex number of the required channel \return the parameter of the coupling constant */ LauParameter& getCouplingParameter(const Int_t poleIndex, const Int_t channelIndex); //! Get scattering constants that were loaded from the input file /*! \param [in] channel1Index number of the first channel index \param [in] channel2Index number of the second channel index \return the value of the scattering constant */ Double_t getScatteringConstant(const Int_t channel1Index, const Int_t channel2Index) const; //! Get scattering parameters, set according to the input file /*! \param [in] channel1Index number of the first channel index \param [in] channel2Index number of the second channel index \return the parameter of the scattering constant */ LauParameter& getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index); //! Get mSq0 production parameter /*! \return the mSq0 parameter */ LauParameter& getmSq0() {return mSq0_;} //! Get s0Scatt production parameter /*! \return the s0Scatt parameter */ LauParameter& gets0Scatt() {return s0Scatt_;} //! Get s0 production parameter /*! \return the s0Prod parameter */ LauParameter& gets0Prod() {return s0Prod_;} //! Get sA production parameter /*! \return the sA parameter */ LauParameter& getsA() {return sA_;} //! Get sA0 production parameter /*! \return the sA0 parameter */ LauParameter& getsA0() {return sA0_;} //! Get the "slowly-varying part" term of the amplitude /*! \return the svp term */ Double_t getProdSVPTerm() const {return prodSVP_;} //! Get the full complex propagator term for a given channel /*! \param [in] channelIndex the number of the required channel \return the complex propagator term */ LauComplex getPropTerm(const Int_t channelIndex) const; //! Get the DP axis identifier /*! \return the value to identify the DP axis in question */ Int_t getResPairAmpInt() const {return resPairAmpInt_;} //! Get the number of channels /*! \return the number of channels */ Int_t getNChannels() const {return nChannels_;} //! Get the number of poles /*! \return the number of poles */ Int_t getNPoles() const {return nPoles_;} //! Get the propagator name /*! \return the name of the propagator */ TString getName() const {return name_;} //! Get the unitary transition amplitude for the given channel /*! \param [in] s The invariant mass squared \param [in] channel The index number of the channel process \return the complex amplitude T */ LauComplex getTransitionAmp(const Double_t s, const Int_t channel); //! Get the complex phase space term for the given channel and invariant mass squared /*! \param [in] s The invariant mass squared \param [in] channel The index number of the channel process \return the complex phase space term rho(channel, channel) */ LauComplex getPhaseSpaceTerm(const Double_t s, const Int_t channel); //! Get the Adler zero factor, which is set when updatePropagator is called /*! \return the Adler zero factor */ Double_t getAdlerZero() const {return adlerZeroFactor_;} //! Get the THat amplitude for the given s and channel number /*! \param [in] s The invariant mass squared \param [in] channel The index number of the channel process \return the complex THat amplitude */ LauComplex getTHat(const Double_t s, const Int_t channel); protected: // Integers to specify the allowed channels for the phase space calculations. // Please keep Zero at the start and leave TotChannels at the end // whenever more channels are added to this. //! Integers to specify the allowed channels for the phase space calculations enum class KMatrixChannels {Zero, PiPi, KK, FourPi, EtaEta, EtaEtaP, KPi, KEtaP, KThreePi, D0K, Dstar0K, TotChannels}; //! Calculate the scattering K-matrix for the given value of s /*! \param [in] s the invariant mass squared */ void calcScattKMatrix(const Double_t s); //! Calculate the real and imaginary part of the phase space density diagonal matrix /*! \param [in] s the invariant mass squared */ void calcRhoMatrix(const Double_t s); //! Retrieve the complex phasespace density for a given channel /*! \param [in] s the invariant mass squared \param [in] phaseSpaceIndex the phasespace index of the channel \return the complex phasespace density */ LauComplex getRho(const Double_t s, const LauKMatrixPropagator::KMatrixChannels) const; //! Calculate the (real) gamma angular-momentum barrier matrix /*! \param [in] s the invariant mass squared */ void calcGammaMatrix(const Double_t s); //! Calculate the gamma angular-momentum barrier /*! + \param [in] iCh the channel index \param [in] s the invariant mass squared \return the centrifugal barrier factor for L=0,1, or 2 */ - Double_t calcGamma(const Int_t iCh, const Double_t s, KMatrixChannels phaseSpaceIndex) const; + Double_t calcGamma(const Int_t iCh, const Double_t s) const; //! Calulate the term 1/(m_pole^2 - s) for the scattering and production K-matrix formulae /*! \param [in] s the invariant mass squared */ void calcPoleDenomVect(const Double_t s); //! Calculate the D0K+ phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcD0KRho(const Double_t s) const; //! Calculate the D*0K+ phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcDstar0KRho(const Double_t s) const; //! Calculate the pipi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcPiPiRho(const Double_t s) const; //! Calculate the KK phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKKRho(const Double_t s) const; //! Calculate the 4 pi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcFourPiRho(const Double_t s) const; //! Calculate the eta-eta phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcEtaEtaRho(const Double_t s) const; //! Calculate the eta-eta' phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcEtaEtaPRho(const Double_t s) const; //! Calculate the Kpi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKPiRho(const Double_t s) const; //! Calculate the K-eta' phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKEtaPRho(const Double_t s) const; //! Calculate the Kpipipi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKThreePiRho(const Double_t s) const; //! Calculate the "slow-varying part" /*! \param [in] s the invariant mass squared \param [in] s0 the invariant mass squared at the Adler zero \return the SVP term */ Double_t calcSVPTerm(const Double_t s, const Double_t s0) const; //! Update the scattering "slowly-varying part" /*! \param [in] s the invariant mass squared */ void updateScattSVPTerm(const Double_t s); //! Update the production "slowly-varying part" /*! \param [in] s the invariant mass squared */ void updateProdSVPTerm(const Double_t s); //! Calculate the multiplicative factor containing severa Adler zero constants /*! \param [in] s the invariant mass squared */ void updateAdlerZeroFactor(const Double_t s); //! Check the phase space factors that need to be used /*! \param [in] phaseSpaceInt phase space types \return true of false */ Bool_t checkPhaseSpaceType(const Int_t phaseSpaceInt) const; //! Get the unitary transition amplitude matrix for the given kinematics /*! \param [in] kinematics The pointer to the constant kinematics */ void getTMatrix(const LauKinematics* kinematics); //! Get the unitary transition amplitude matrix for the given kinematics /*! \param [in] s The invariant mass squared of the system */ void getTMatrix(const Double_t s); //! Get the square root of the phase space matrix void getSqrtRhoMatrix(); private: //! Copy constructor (not implemented) LauKMatrixPropagator(const LauKMatrixPropagator& rhs)=delete; //! Copy assignment operator (not implemented) LauKMatrixPropagator& operator=(const LauKMatrixPropagator& rhs)=delete; //! Initialise and set the dimensions for the internal matrices and parameter arrays void initialiseMatrices(); //! Store the (phase space) channel indices from a line in the parameter file /*! \param [in] theLine Vector of strings corresponding to the line from the parameter file */ void storeChannels(const std::vector& theLine); //! Store the pole mass and couplings from a line in the parameter file /*! \param [in] theLine Vector of strings corresponding to the line from the parameter file */ void storePole(const std::vector& theLine); //! Store the scattering coefficients from a line in the parameter file /*! \param [in] theLine Vector of strings corresponding to the line from the parameter file */ void storeScattering(const std::vector& theLine); //! Store the channels' characteristic radii from a line in the parameter file /*! \param [in] theLine Vector of strings corresponding to the line from the parameter file */ void storeRadii(const std::vector& theLine); //! Store the channels' angular momenta from a line in the parameter file /*! \param [in] theLine Vector of strings corresponding to the line from the parameter file \param [in] a Vector of integers corresponding to parameter in the propagator denominator */ void storeOrbitalAngularMomenta(const std::vector& theLine, std::vector& a); //! Store the barrier-factor parameter from a line in the parameter file /*! \param [in] theLine Vector of strings corresponding to the line from the parameter file \param [in] a Vector of integers corresponding to parameter in the propagator denominator */ void storeBarrierFactorParameter(const std::vector& theLine, std::vector& a); //! Store miscelleanous parameters from a line in the parameter file /*! \param [in] keyword the name of the parameter to be set \param [in] parString the string containing the value of the parameter */ void storeParameter(const TString& keyword, const TString& parString); //! String to store the propagator name TString name_; //! Name of the input parameter file TString paramFileName_; //! Number to identify the DP axis in question Int_t resPairAmpInt_; //! Row index - 1 Int_t index_; //! s value of the previous pole Double_t previousS_{0.0}; //! "slowly-varying part" for the scattering K-matrix Double_t scattSVP_{0.0}; //! "slowly-varying part" for the production K-matrix Double_t prodSVP_{0.0}; //! Real part of the propagator matrix TMatrixD realProp_; //! Imaginary part of the propagator matrix TMatrixD negImagProp_; //! Scattering K-matrix TMatrixD ScattKMatrix_; //! Real part of the phase space density diagonal matrix TMatrixD ReRhoMatrix_; //! Imaginary part of the phase space density diagonal matrix TMatrixD ImRhoMatrix_; //! Gamma angular-momentum barrier matrix TMatrixD GammaMatrix_; //! Identity matrix TMatrixD IMatrix_; //! Null matrix TMatrixD zeroMatrix_; //! Real part of the square root of the phase space density diagonal matrix TMatrixD ReSqrtRhoMatrix_; //! Imaginary part of the square root of the phase space density diagonal matrix TMatrixD ImSqrtRhoMatrix_; //! Real part of the unitary T matrix TMatrixD ReTMatrix_; //! Imaginary part of the unitary T matrix TMatrixD ImTMatrix_; //! Number of channels Int_t nChannels_; //! Number of poles Int_t nPoles_; //! Vector of orbital angular momenta std::vector L_; //! Boolean to indicate whether storeBarrierFactorParameter has been called Bool_t haveCalled_storeBarrierFactorParameter{kFALSE}; //! Boolean to dictate whether to include Blatt-Weisskopf-like denominator in K-matrix centrifugal barrier factor Bool_t includeBWBarrierFactor_{kTRUE}; //! Vector of squared pole masses std::vector mSqPoles_; //! Array of coupling constants LauParArray gCouplings_; //! Array of scattering SVP values LauParArray fScattering_; //! Vector of characteristic radii std::vector radii_; //! Vector of ratio a/R^2 std::vector gamAInvRadSq_; //! Vector of phase space types std::vector phaseSpaceTypes_; //! Vector of squared masses std::vector mSumSq_; //! Vector of mass differences std::vector mDiffSq_; //! Vector of 1/(m_pole^2 - s) terms for scattering and production K-matrix formulae std::vector poleDenomVect_; //! Constant from input file LauParameter mSq0_; //! Constant from input file LauParameter s0Scatt_; //! Constant from input file LauParameter s0Prod_; //! Constant from input file LauParameter sA_; //! Constant from input file LauParameter sA0_; //! Defined as 0.5*sA*mPi*mPi Double_t sAConst_{0.0}; //! Defined as 4*mPi*mPi const Double_t m2piSq_{4.0*LauConstants::mPiSq}; //! Defined as 4*mK*mK const Double_t m2KSq_{4.0*LauConstants::mKSq}; //! Defined as 4*mEta*mEta const Double_t m2EtaSq_{4.0*LauConstants::mEtaSq}; //! Defined as (mEta+mEta')^2 const Double_t mEtaEtaPSumSq_{(LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)}; //! Defined as (mEta-mEta')^2 const Double_t mEtaEtaPDiffSq_{(LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)}; //! Defined as (mK+mPi)^2 const Double_t mKpiSumSq_{(LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)}; //! Defined as (mK-mPi)^2 const Double_t mKpiDiffSq_{(LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)}; //! Defined as (mK+mEta')^2 const Double_t mKEtaPSumSq_{(LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)}; //! Defined as (mK-mEta')^2 const Double_t mKEtaPDiffSq_{(LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)}; //! Defined as (mK-3*mPi)^2 const Double_t mK3piDiffSq_{(LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)}; //! Factor used to calculate the Kpipipi phase space term const Double_t k3piFactor_{TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)}; //! Factor used to calculate the pipipipi phase space term const Double_t fourPiFactor1_{16.0*LauConstants::mPiSq}; //! Factor used to calculate the pipipipi phase space term const Double_t fourPiFactor2_{TMath::Sqrt(1.0 - fourPiFactor1_)}; //! Defined as (mD0+mK)^2 const Double_t mD0KSumSq_{(LauConstants::mD0 + LauConstants::mK)*(LauConstants::mD0 + LauConstants::mK)}; //! Defined as (mD0-mK)^2 const Double_t mD0KDiffSq_{(LauConstants::mD0 - LauConstants::mK)*(LauConstants::mD0 - LauConstants::mK)}; //! Defined as (mD*0+mK)^2 const Double_t mDstar0KSumSq_{(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)}; //! Defined as (mD*0-mK)^2 const Double_t mDstar0KDiffSq_{(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)}; //! Multiplicative factor containing various Adler zero constants Double_t adlerZeroFactor_{0.0}; //! Tracks if all params have been set Bool_t parametersSet_{kFALSE}; //! Control the output of the functions static constexpr Bool_t verbose_{kFALSE}; //! Control if scattering constants are channel symmetric: f_ji = f_ij Bool_t scattSymmetry_{kFALSE}; ClassDef(LauKMatrixPropagator,0) // K-matrix amplitude model }; #endif diff --git a/src/LauKMatrixPropagator.cc b/src/LauKMatrixPropagator.cc index b45d833..0955961 100644 --- a/src/LauKMatrixPropagator.cc +++ b/src/LauKMatrixPropagator.cc @@ -1,1425 +1,1425 @@ /* 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 */ /*! \file LauKMatrixPropagator.cc \brief File containing implementation of LauKMatrixPropagator class. */ #include "LauKMatrixPropagator.hh" #include "LauTextFileParser.hh" #include "LauKinematics.hh" #include "LauComplex.hh" #include "TMath.h" #include "TSystem.h" #include #include #include #include using std::cout; using std::endl; using std::cerr; ClassImp(LauKMatrixPropagator) LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile, Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles, Int_t rowIndex) : name_(name), paramFileName_(paramFile), resPairAmpInt_(resPairAmpInt), index_(rowIndex - 1), nChannels_(nChannels), nPoles_(nPoles) { // Constructor // Check that the index is OK if (index_ < 0 || index_ >= nChannels_) { std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to " << rowIndex << ", must be between 1 and the number of channels " << nChannels_ << std::endl; gSystem->Exit(EXIT_FAILURE); } this->setParameters(paramFile); } LauKMatrixPropagator::~LauKMatrixPropagator() { // Destructor realProp_.Clear(); negImagProp_.Clear(); ScattKMatrix_.Clear(); ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear(); GammaMatrix_.Clear(); ReTMatrix_.Clear(); ImTMatrix_.Clear(); IMatrix_.Clear(); zeroMatrix_.Clear(); } LauComplex LauKMatrixPropagator::getPropTerm(const Int_t channel) const { // Get the (i,j) = (index_, channel) term of the propagator // matrix. This allows us not to return the full propagator matrix. Double_t realTerm = this->getRealPropTerm(channel); Double_t imagTerm = this->getImagPropTerm(channel); LauComplex propTerm(realTerm, imagTerm); return propTerm; } Double_t LauKMatrixPropagator::getRealPropTerm(const Int_t channel) const { // Get the real part of the (i,j) = (index_, channel) term of the propagator // matrix. This allows us not to return the full propagator matrix. if (parametersSet_ == kFALSE) {return 0.0;} Double_t propTerm = realProp_[index_][channel]; return propTerm; } Double_t LauKMatrixPropagator::getImagPropTerm(const Int_t channel) const { // Get the imaginary part of the (i,j) = (index_, channel) term of the propagator // matrix. This allows us not to return the full propagator matrix. if (parametersSet_ == kFALSE) {return 0.0;} Double_t imagTerm = -1.0*negImagProp_[index_][channel]; return imagTerm; } void LauKMatrixPropagator::updatePropagator(const Double_t s) { // Calculate the K-matrix propagator for the given s value. // The K-matrix amplitude is given by // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix. // i = index for the state (e.g. S-wave index = 0). // Here, we only find the (I - iK*rho)^-1 matrix part. // Check if we have almost the same s value as before. If so, don't re-calculate // the propagator nor any of the pole mass summation terms. if (TMath::Abs(s - previousS_) < 1e-6*s) { //cout<<"Already got propagator for s = "<calcPoleDenomVect(s); this->updateAdlerZeroFactor(s); // Calculate the scattering K-matrix (real and symmetric) this->calcScattKMatrix(s); // Calculate the phase space density matrix, which is diagonal, but can be complex // if the quantity s is below various threshold values (analytic continuation). this->calcRhoMatrix(s); // Calculate the angular momentum barrier matrix, which is real and diagonal this->calcGammaMatrix(s); // Calculate K*rho*(gamma^2) (real and imaginary parts, since rho can be complex) TMatrixD GammaMatrixSq = (GammaMatrix_*GammaMatrix_); TMatrixD K_realRhoGammaSq(ScattKMatrix_); K_realRhoGammaSq *= ReRhoMatrix_; K_realRhoGammaSq *= GammaMatrixSq; TMatrixD K_imagRhoGammaSq(ScattKMatrix_); K_imagRhoGammaSq *= ImRhoMatrix_; K_imagRhoGammaSq *= GammaMatrixSq; // A = I + K*Imag(rho)Gamma^2, B = -K*Real(Rho)Gamma^2 // Calculate C and D matrices such that (A + iB)*(C + iD) = I, // ie. C + iD = (I - i K*rhoGamma^2)^-1, separated into real and imaginary parts. // realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C TMatrixD A(IMatrix_); A += K_imagRhoGammaSq; TMatrixD B(zeroMatrix_); B -= K_realRhoGammaSq; TMatrixD invA(TMatrixD::kInverted, A); TMatrixD invA_B(invA); invA_B *= B; TMatrixD B_invA_B(B); B_invA_B *= invA_B; TMatrixD invC(A); invC += B_invA_B; // Set the real part of the propagator matrix ("C") realProp_ = TMatrixD(TMatrixD::kInverted, invC); // Set the (negative) imaginary part of the propagator matrix ("-D") TMatrixD BC(B); BC *= realProp_; negImagProp_ = TMatrixD(invA); negImagProp_ *= BC; // Pre-multiply by the Gamma matrix: realProp_ = GammaMatrix_ * realProp_; negImagProp_ = GammaMatrix_ * negImagProp_; if(verbose_) { std::cout << "In LauKMatrixPropagator::updatePropagator(s). D[1-iKrhoD^2]^-1: " << std::endl; TString realOutput("Real part:"), imagOutput("Imag part:"); for (int iChannel = 0; iChannel < nChannels_; iChannel++) { for (int jChannel = 0; jChannel < nChannels_; jChannel++) { realOutput += Form("\t%.6f",realProp_[iChannel][jChannel]); imagOutput += Form("\t%.6f",-1*negImagProp_[iChannel][jChannel]); } realOutput += "\n "; imagOutput += "\n "; } std::cout << realOutput << std::endl; std::cout << imagOutput << std::endl; } // Also calculate the production SVP term, since this uses Adler-zero parameters // defined in the parameter file. this->updateProdSVPTerm(s); // Finally, keep track of the value of s we just used. previousS_ = s; } void LauKMatrixPropagator::setParameters(const TString& inputFile) { // Read an input file that specifies the values of the coupling constants g_i for // the given number of poles and their (bare) masses. Also provided are the f_{ab} // slow-varying constants. The input file should also provide the Adler zero // constants s_0, s_A and s_A0. parametersSet_ = kFALSE; cout<<"Initialising K-matrix propagator "<Exit(EXIT_FAILURE); } UInt_t iLine(0); for (iLine = 1; iLine <= nTotLines; iLine++) { // Get the line of strings std::vector theLine = readFile.getLine(iLine); // There should always be at least two strings: a keyword and at least 1 value if (theLine.size() < 2) {continue;} TString keyword(theLine[0].c_str()); keyword.ToLower(); // Use lowercase if (!keyword.CompareTo("channels")) { // Channel indices for phase-space factors this->storeChannels(theLine); } else if (!keyword.CompareTo("pole")) { // Pole terms this->storePole(theLine); } else if (!keyword.CompareTo("scatt")) { // Scattering terms this->storeScattering(theLine); } else if (!keyword.CompareTo("angularmomentum")) { // Orbital angular momentum state for each channel & set default a values if called before storeBarrierFactorParameter this->storeOrbitalAngularMomenta(theLine, a); } else if (!keyword.CompareTo("barrierfactorparameter")) { // Value of parameter "a" in denominator of centrifugal barrier factor, gamma this->storeBarrierFactorParameter(theLine, a); } else if (!keyword.CompareTo("radii")) { // Values of characteristic radius this->storeRadii(theLine); } else { // Usually Adler-zero constants TString parString(theLine[1].c_str()); this->storeParameter(keyword, parString); } } sAConst_ = 0.5*sA_.unblindValue()*LauConstants::mPiSq; // Symmetrise scattering parameters if enabled if (scattSymmetry_ == kTRUE) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { for (Int_t jChannel = iChannel; jChannel < nChannels_; jChannel++) { LauParameter fPar = fScattering_[iChannel][jChannel]; fScattering_[jChannel][iChannel] = LauParameter(fPar); } } } // Now that radii and barrier-factor-denominator parameters have been set, cache the value of "a/(R*R)" for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { gamAInvRadSq_[iChannel] = a[iChannel]/(radii_[iChannel]*radii_[iChannel]); } // All required parameters have been set parametersSet_ = kTRUE; cout<<"Finished initialising K-matrix propagator "< >. // Set their sizes using the number of poles and channels defined in the constructor gCouplings_.clear(); gCouplings_.resize(nPoles_); for (Int_t iPole = 0; iPole < nPoles_; iPole++) { gCouplings_[iPole].resize(nChannels_); } fScattering_.clear(); fScattering_.resize(nChannels_); for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { fScattering_[iChannel].resize(nChannels_); } // Clear other vectors phaseSpaceTypes_.clear(); phaseSpaceTypes_.resize(nChannels_); mSqPoles_.clear(); mSqPoles_.resize(nPoles_); haveCalled_storeBarrierFactorParameter = kFALSE; } void LauKMatrixPropagator::storeChannels(const std::vector& theLine) { // Get the list of channel indices to specify what phase space factors should be used // e.g. pipi, Kpi, eta eta', 4pi etc.. // Check that the line has nChannels_+1 strings Int_t nTypes = static_cast(theLine.size()) - 1; if (nTypes != nChannels_) { cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines " <checkPhaseSpaceType(phaseSpaceInt); if (checkChannel == kTRUE) { cout<<"Adding phase space channel index "<(phaseSpaceInt); } else { cerr<<"Phase space channel index "<(LauKMatrixPropagator::KMatrixChannels::TotChannels)-1<& theLine) { // Store the pole mass and its coupling constants for each channel. // Each line will contain: Pole poleNumber poleMass poleCouplingsPerChannel // Check that the line has nChannels_ + 3 strings Int_t nWords = static_cast(theLine.size()); Int_t nExpect = nChannels_ + 3; if (nWords == nExpect) { Int_t poleIndex = std::atoi(theLine[1].c_str()) - 1; if (poleIndex >= 0 && poleIndex < nPoles_) { Double_t poleMass = std::atof(theLine[2].c_str()); Double_t poleMassSq = poleMass*poleMass; LauParameter mPoleParam(Form("KM_%s_poleMassSq_%i",name_.Data(),poleIndex),poleMassSq); mSqPoles_[poleIndex] = mPoleParam; cout<<"Added bare pole mass "<& theLine) { // Store the scattering constants (along one of the channel rows). // Each line will contain: Scatt ScattIndex ScattConstantsPerChannel // Check that the line has nChannels_ + 2 strings Int_t nWords = static_cast(theLine.size()); Int_t nExpect = nChannels_ + 2; if (nWords == nExpect) { Int_t scattIndex = std::atoi(theLine[1].c_str()) - 1; if (scattIndex >= 0 && scattIndex < nChannels_) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { Double_t scattConst = std::atof(theLine[iChannel+2].c_str()); LauParameter scattParam(Form("KM_%s_fScatteringConst_%i_%i",name_.Data(),scattIndex,iChannel),scattConst); fScattering_[scattIndex][iChannel] = scattParam; cout<<"Added scattering parameter f("<& theLine, std::vector& a) { // Store the orbital angular momentum for each channel // Each line will contain: angularmomentum OrbitalAngularMomentumPerChannel // Check that the line has nChannels_ + 1 strings Int_t nWords = static_cast(theLine.size()); Int_t nExpect = nChannels_ + 1; if (nWords == nExpect) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { Int_t angularMomentum = std::atoi(theLine[iChannel+1].c_str()); L_[iChannel] = angularMomentum; cout<<"Defined K-matrix orbital angular momentum "<Exit(EXIT_FAILURE); } } } } void LauKMatrixPropagator::storeRadii(const std::vector& theLine) { // Store the characteristic radii (measured in GeV^{-1}) // Each line will contain: Radii RadiusConstantsPerChannel // Check that the line has nChannels_ + 1 strings Int_t nWords = static_cast(theLine.size()); Int_t nExpect = nChannels_ + 1; if (nWords == nExpect) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { Double_t radiusConst = std::atof(theLine[iChannel+1].c_str()); radii_[iChannel] = radiusConst; cout<<"Added K-matrix radius "<& theLine, std::vector& a) { // Store the parameter of the barrier factor // Each line will contain: barrierfactorparameter ParameterValuePerchannel // Check that the line has nChannels_ + 1 strings Int_t nWords = static_cast(theLine.size()); Int_t nExpect = nChannels_ + 1; if (nWords == nExpect) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { Double_t parameterValue = std::atof(theLine[iChannel+1].c_str()); a[iChannel] = parameterValue; cout<<"Added K-matrix barrier factor parameter value "<updateScattSVPTerm(s); // Now loop over iChannel, jChannel to calculate Kij = Kji. for (iChannel = 0; iChannel < nChannels_; iChannel++) { // Scattering matrix is real and symmetric. Start j loop from i. for (jChannel = iChannel; jChannel < nChannels_; jChannel++) { Double_t Kij(0.0); // Calculate pole mass summation term for (iPole = 0; iPole < nPoles_; iPole++) { Double_t g_i = this->getCouplingConstant(iPole, iChannel); Double_t g_j = this->getCouplingConstant(iPole, jChannel); Kij += poleDenomVect_[iPole]*g_i*g_j; if (verbose_) {cout<<"1: Kij for i = "<getScatteringConstant(iChannel, jChannel); Kij += fij*scattSVP_; Kij *= adlerZeroFactor_; if (verbose_) {cout<<"2: Kij for i = "< 1.0e-6) {invPoleTerm = 1.0/poleTerm;} poleDenomVect_.push_back(invPoleTerm); } } Double_t LauKMatrixPropagator::getPoleDenomTerm(const Int_t poleIndex) const { if (parametersSet_ == kFALSE) {return 0.0;} Double_t poleDenom(0.0); poleDenom = poleDenomVect_[poleIndex]; return poleDenom; } LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(const Int_t poleIndex) { if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) ) { std::cerr << "ERROR from LauKMatrixPropagator::getPoleMassSqParameter(). Invalid pole." << std::endl; gSystem->Exit(EXIT_FAILURE); } return mSqPoles_[poleIndex]; } Double_t LauKMatrixPropagator::getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const { if (parametersSet_ == kFALSE) {return 0.0;} if (poleIndex < 0 || poleIndex >= nPoles_) {return 0.0;} if (channelIndex < 0 || channelIndex >= nChannels_) {return 0.0;} Double_t couplingConst = gCouplings_[poleIndex][channelIndex].unblindValue(); return couplingConst; } LauParameter& LauKMatrixPropagator::getCouplingParameter(const Int_t poleIndex, const Int_t channelIndex) { if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) || (channelIndex < 0 || channelIndex >= nChannels_) ) { std::cerr << "ERROR from LauKMatrixPropagator::getCouplingParameter(). Invalid coupling." << std::endl; gSystem->Exit(EXIT_FAILURE); } //std::cout << "Minvalue + range for " << poleIndex << ", " << channelIndex << ": " << gCouplings_[poleIndex][channelIndex].minValue() << " => + " << gCouplings_[poleIndex][channelIndex].range() << // " and init value: " << gCouplings_[poleIndex][channelIndex].initValue() << std::endl; return gCouplings_[poleIndex][channelIndex]; } Double_t LauKMatrixPropagator::getScatteringConstant(const Int_t channel1Index, const Int_t channel2Index) const { if (parametersSet_ == kFALSE) {return 0.0;} if (channel1Index < 0 || channel1Index >= nChannels_) {return 0.0;} if (channel2Index < 0 || channel2Index >= nChannels_) {return 0.0;} Double_t scatteringConst = fScattering_[channel1Index][channel2Index].unblindValue(); return scatteringConst; } LauParameter& LauKMatrixPropagator::getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index) { if ( (parametersSet_ == kFALSE) || (channel1Index < 0 || channel1Index >= nChannels_) || (channel2Index < 0 || channel2Index >= nChannels_) ) { std::cerr << "ERROR from LauKMatrixPropagator::getScatteringParameter(). Invalid chanel index." << std::endl; gSystem->Exit(EXIT_FAILURE); } return fScattering_[channel1Index][channel2Index]; } Double_t LauKMatrixPropagator::calcSVPTerm(const Double_t s, const Double_t s0) const { if (parametersSet_ == kFALSE) {return 0.0;} // Calculate the "slowly-varying part" (SVP) Double_t result(0.0); Double_t deltaS = s - s0; if (TMath::Abs(deltaS) > 1.0e-6) { result = (mSq0_.unblindValue() - s0)/deltaS; } return result; } void LauKMatrixPropagator::updateScattSVPTerm(const Double_t s) { // Update the scattering "slowly-varying part" (SVP) Double_t s0Scatt = s0Scatt_.unblindValue(); scattSVP_ = this->calcSVPTerm(s, s0Scatt); } void LauKMatrixPropagator::updateProdSVPTerm(const Double_t s) { // Update the production "slowly-varying part" (SVP) Double_t s0Prod = s0Prod_.unblindValue(); prodSVP_ = this->calcSVPTerm(s, s0Prod); } void LauKMatrixPropagator::updateAdlerZeroFactor(const Double_t s) { // Calculate the multiplicative factor containing various Adler zero // constants. adlerZeroFactor_ = 0.0; Double_t sA0Val = sA0_.unblindValue(); Double_t deltaS = s - sA0Val; if (TMath::Abs(deltaS) > 1e-6) { adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS; } } void LauKMatrixPropagator::calcGammaMatrix(const Double_t s) { // Calculate the gamma angular momentum barrier matrix // for the given invariant mass squared quantity, s. // Initialise all entries to zero GammaMatrix_.Zero(); Double_t gamma(0.0); for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { - LauKMatrixPropagator::KMatrixChannels phaseSpaceIndex = phaseSpaceTypes_[iChannel]; - if ( L_[iChannel] != 0 ) { - gamma = this->calcGamma(iChannel,s,phaseSpaceIndex); + gamma = this->calcGamma(iChannel,s); } else { gamma = 1.0; // S-wave } if (verbose_) { cout<<"GammaMatrix("<Exit(EXIT_FAILURE); } return rho; } LauComplex LauKMatrixPropagator::calcD0KRho(const Double_t s) const { // Calculate the D0K+ phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mD0KSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mD0KDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcDstar0KRho(const Double_t s) const { // Calculate the Dstar0K+ phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mDstar0KSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mDstar0KDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcPiPiRho(const Double_t s) const { // Calculate the pipi phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-m2piSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKKRho(const Double_t s) const { // Calculate the KK phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-m2KSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcFourPiRho(const Double_t s) const { // Calculate the 4pi phase space factor. This uses a 6th-order polynomial // parameterisation that approximates the multi-body phase space double integral // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2. // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2 // for s1 and s2, the invariant energy squared of each of the di-pion states, // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be // 0.775 GeV and the energy-dependent width of the 4pi system // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state). // The normalisation term rho_0 is found by ensuring that the phase space integral // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5. // Plotting the value of this double integral as a function of s can then be fitted // to a 6th-order polynomial (for s < 1), which is the result used below LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} if (s <= 1.0) { Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840; rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051; // For some values of s (below 2*mpi), this term is a very small // negative number. Check for this and set the rho term to zero. if (rhoTerm < 0.0) {rhoTerm = 0.0;} rho.setRealPart( rhoTerm ); } else { rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) ); } return rho; } LauComplex LauKMatrixPropagator::calcEtaEtaRho(const Double_t s) const { // Calculate the eta-eta phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-m2EtaSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcEtaEtaPRho(const Double_t s) const { // Calculate the eta-eta' phase space factor. Note that the // mass difference term m_eta - m_eta' is not included, // since this corresponds to a "t or u-channel crossing", // which means that we cannot simply analytically continue // this part of the phase space factor below threshold, which // we can do for s-channel contributions. This is actually an // unsolved problem, e.g. see Guo et al 1409.8652, and // Danilkin et al 1409.7708. Anisovich and Sarantsev in // hep-ph/0204328 "solve" this issue by setting the mass // difference term to unity, which is what we do here... LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-mEtaEtaPSumSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKPiRho(const Double_t s) const { // Calculate the K-pi phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mKpiSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKEtaPRho(const Double_t s) const { // Calculate the K-eta' phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKThreePiRho(const Double_t s) const { // Calculate the Kpipipi + multimeson phase space factor. // Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form // used for the rest of that paper (thankfully, the amplitude does not depend // significantly on the form used for the K3pi phase space factor). LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} if (s < 1.44) { Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0; if (powerTerm < 0.0) { rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) ); } else { rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) ); } } else { rho.setRealPart( 1.0 ); } return rho; } Bool_t LauKMatrixPropagator::checkPhaseSpaceType(const Int_t phaseSpaceInt) const { Bool_t passed(kFALSE); if (phaseSpaceInt >= 1 && phaseSpaceInt < static_cast(LauKMatrixPropagator::KMatrixChannels::TotChannels)) { passed = kTRUE; } return passed; } LauComplex LauKMatrixPropagator::getTransitionAmp(const Double_t s, const Int_t channel) { // Get the complex (unitary) transition amplitude T for the given channel LauComplex TAmp(0.0, 0.0); if (channel <= 0 || channel > nChannels_) {return TAmp;} this->getTMatrix(s); TAmp.setRealPart(ReTMatrix_[index_][channel-1]); TAmp.setImagPart(ImTMatrix_[index_][channel-1]); return TAmp; } LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(const Double_t s, const Int_t channel) { // Get the complex (unitary) transition amplitude T for the given channel LauComplex rho(0.0, 0.0); if (channel <= 0 || channel > nChannels_) {return rho;} // If s has changed from the previous value, recalculate rho if (TMath::Abs(s - previousS_) > 1e-6*s) { this->calcRhoMatrix(s); } rho.setRealPart(ReRhoMatrix_[channel][channel-1]); rho.setImagPart(ImRhoMatrix_[channel][channel-1]); return rho; } void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) { // Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho), // and T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix, // which has phase-space factors included (rho). This function is not // needed to calculate the K-matrix amplitudes, but allows us // to check the variation of T as a function of s (kinematics) if (!kinematics) {return;} // Get the invariant mass squared (s) from the kinematics object. // Use the resPairAmpInt to find which mass-squared combination to use. Double_t s(0.0); if (resPairAmpInt_ == 1) { s = kinematics->getm23Sq(); } else if (resPairAmpInt_ == 2) { s = kinematics->getm13Sq(); } else if (resPairAmpInt_ == 3) { s = kinematics->getm12Sq(); } this->getTMatrix(s); } void LauKMatrixPropagator::getTMatrix(const Double_t s) { // Find the unitary transition T matrix, where // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and // T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix, // which has phase-space factors included (rho). Note that the first // sqrt of the rho matrix is complex conjugated. // This function is not needed to calculate the K-matrix amplitudes, but // allows us to check the variation of T as a function of s (kinematics) // Initialse the real and imaginary parts of the T matrix to zero ReTMatrix_.Zero(); ImTMatrix_.Zero(); if (parametersSet_ == kFALSE) {return;} // Update K, rho and the propagator (I - i K rho)^-1 this->updatePropagator(s); // Find the real and imaginary T_hat matrices TMatrixD THatReal = realProp_*ScattKMatrix_; TMatrixD THatImag(zeroMatrix_); THatImag -= negImagProp_*ScattKMatrix_; // Find the square-root of the phase space matrix this->getSqrtRhoMatrix(); // Let sqrt(rho) = A + iB and T_hat = C + iD // => T = A(CA-DB) + B(DA+CB) + i[A(DA+CB) + B(DB-CA)] TMatrixD CA(THatReal); CA *= ReSqrtRhoMatrix_; TMatrixD DA(THatImag); DA *= ReSqrtRhoMatrix_; TMatrixD CB(THatReal); CB *= ImSqrtRhoMatrix_; TMatrixD DB(THatImag); DB *= ImSqrtRhoMatrix_; TMatrixD CAmDB(CA); CAmDB -= DB; TMatrixD DApCB(DA); DApCB += CB; TMatrixD DBmCA(DB); DBmCA -= CA; // Find the real and imaginary parts of the transition matrix T ReTMatrix_ = ReSqrtRhoMatrix_*CAmDB + ImSqrtRhoMatrix_*DApCB; ImTMatrix_ = ReSqrtRhoMatrix_*DApCB + ImSqrtRhoMatrix_*DBmCA; } void LauKMatrixPropagator::getSqrtRhoMatrix() { // Find the square root of the (current) phase space matrix so that // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first // sqrt of rho matrix is complex conjugated // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2). // Since rho is diagonal, then the square root of rho will also be diagonal, // with its real and imaginary matrix elements equal to c and d, respectively // Initialise the real and imaginary parts of the square root of // the rho matrix to zero ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero(); for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { Double_t realRho = ReRhoMatrix_[iChannel][iChannel]; Double_t imagRho = ImRhoMatrix_[iChannel][iChannel]; Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho); Double_t rhoSum = rhoMag + realRho; Double_t rhoDiff = rhoMag - realRho; Double_t reSqrtRho(0.0), imSqrtRho(0.0); if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);} if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);} ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho; ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho; } } LauComplex LauKMatrixPropagator::getTHat(const Double_t s, const Int_t channel) { LauComplex THat(0.0, 0.0); if (channel <= 0 || channel > nChannels_) {return THat;} this->updatePropagator(s); // Find the real and imaginary T_hat matrices TMatrixD THatReal = realProp_*ScattKMatrix_; TMatrixD THatImag(zeroMatrix_); THatImag -= negImagProp_*ScattKMatrix_; // Return the specific THat component THat.setRealPart(THatReal[index_][channel-1]); THat.setImagPart(THatImag[index_][channel-1]); return THat; }