LauIsobarDynamics.hh

Go to the documentation of this file.

// Copyright University of Warwick 2005 - 2014.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

// Authors:
// Thomas Latham
// John Back
// Paul Harrison

#ifndef LAU_ISOBAR_DYNAMICS
#define LAU_ISOBAR_DYNAMICS

#include <set>
#include <vector>

#include "TString.h"

#include "LauAbsResonance.hh"
#include "LauComplex.hh"

class LauCacheData;
class LauDaughters;
class LauAbsEffModel;
class LauAbsIncohRes;
class LauFitDataTree;
class LauKMatrixPropagator;
class LauDPPartialIntegralInfo;
class LauKinematics;

class LauIsobarDynamics {

        public:
                typedef std::map<Int_t,LauAbsEffModel*> LauTagCatScfFractionModelMap;

                enum ToyMCStatus {
                        GenOK,        
                        MaxIterError, 
                        ASqMaxError   
                };


                LauIsobarDynamics(LauDaughters* daughters, LauAbsEffModel* effModel, LauAbsEffModel* scfFractionModel = 0);


                LauIsobarDynamics(LauDaughters* daughters, LauAbsEffModel* effModel, LauTagCatScfFractionModelMap scfFractionModel);

                virtual ~LauIsobarDynamics();


                void initialise(const std::vector<LauComplex>& coeffs);

                void recalculateNormalisation();


                inline void setIntFileName(const TString& fileName) {intFileName_ = fileName;}

                // Integration

                void setIntegralBinWidths(const Double_t m13BinWidth, const Double_t m23BinWidth,
                                          const Double_t mPrimeBinWidth = 0.001, const Double_t thPrimeBinWidth = 0.001);


                void setNarrowResonanceThreshold(const Double_t narrowWidth) { narrowWidth_ = narrowWidth; }


                void setIntegralBinningFactor(const Double_t binningFactor) { binningFactor_ = binningFactor; }


                void forceSymmetriseIntegration(const Bool_t force) { forceSymmetriseIntegration_ = force; }


                LauAbsResonance* addResonance(const TString& resName, const Int_t resPairAmpInt, const LauAbsResonance::LauResonanceModel resType, const LauBlattWeisskopfFactor::BlattWeisskopfCategory bwCategory = LauBlattWeisskopfFactor::Default);


                LauAbsResonance* addIncoherentResonance(const TString& resName, const Int_t resPairAmpInt, const LauAbsResonance::LauResonanceModel resType);


                void defineKMatrixPropagator(const TString& propName, const TString& paramFileName,
                                             Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles, Int_t rowIndex = 1);


                void addKMatrixProdPole(const TString& poleName, const TString& propName, Int_t poleIndex, Bool_t useProdAdler = kFALSE);


                void addKMatrixProdSVP(const TString& SVPName, const TString& propName, Int_t channelIndex, Bool_t useProdAdler = kFALSE);


                inline void setASqMaxValue(Double_t value) {aSqMaxSet_ = value;}


                inline Double_t getASqMaxSetValue() const { return aSqMaxSet_; }


                inline Double_t getASqMaxVarValue() const { return aSqMaxVar_; }


                Bool_t generate();


                ToyMCStatus checkToyMC(Bool_t printErrorMessages = kTRUE, Bool_t printInfoMessages = kFALSE);


                inline Int_t maxGenIterations() const {return iterationsMax_;}


                void calcLikelihoodInfo(const UInt_t iEvt);


                void calcLikelihoodInfo(const Double_t m13Sq, const Double_t m23Sq);


                void calcLikelihoodInfo(const Double_t m13Sq, const Double_t m23Sq, const Int_t tagCat);


                void calcExtraInfo(const Bool_t init = kFALSE);


                Bool_t gotReweightedEvent();


                Double_t getEventWeight();


                inline const LauComplex& getEvtDPAmp() const {return totAmp_;}


                inline Double_t getEvtm13Sq() const {return m13Sq_;}


                inline Double_t getEvtm23Sq() const {return m23Sq_;}


                inline Double_t getEvtmPrime() const {return mPrime_;}


                inline Double_t getEvtthPrime() const {return thPrime_;}


                inline Double_t getEvtEff() const {return eff_;}


                inline Double_t getEvtScfFraction() const {return scfFraction_;}


                inline Double_t getEvtJacobian() const {return jacobian_;}


                inline Double_t getEvtIntensity() const {return ASq_;}


                inline Double_t getEvtLikelihood() const {return evtLike_;}


                inline LauComplex getDynamicAmp(const Int_t resID) const {return ff_[resID].scale(fNorm_[resID]);}


                inline LauComplex getFullAmplitude(const Int_t resID) const {return Amp_[resID] * this->getDynamicAmp(resID);}


                inline const std::vector< std::vector<LauComplex> >& getFiFjSum() const {return fifjSum_;}


                inline const std::vector< std::vector<LauComplex> >& getFiFjEffSum() const {return fifjEffSum_;}


                inline const std::vector<Double_t>& getFNorm() const {return fNorm_;}


                void fillDataTree(const LauFitDataTree& fitDataTree);

                void modifyDataTree();


                Bool_t hasResonance(const TString& resName) const;


                Int_t resonanceIndex(const TString& resName) const;


                TString getConjResName(const TString& resName) const;


                const LauAbsResonance* findResonance(const TString& resName) const;


                const LauAbsResonance* getResonance(const UInt_t resIndex) const;


                void updateCoeffs(const std::vector<LauComplex>& coeffs);


                inline void flipHelicityForCPEigenstates(Bool_t boolean) {flipHelicity_ = boolean;}


                inline const LauParameter& getMeanEff() const {return meanDPEff_;}


                inline const LauParameter& getDPRate() const {return DPRate_;}


                inline const LauParArray& getFitFractions() const {return fitFrac_;}


                inline const LauParArray& getFitFractionsEfficiencyUncorrected() const {return fitFracEffUnCorr_;}


                inline UInt_t getnTotAmp() const {return nAmp_+nIncohAmp_;}


                inline UInt_t getnCohAmp() const {return nAmp_;}


                inline UInt_t getnIncohAmp() const {return nIncohAmp_;}


                inline Double_t getDPNorm() const {return DPNorm_;}


                inline const LauDaughters* getDaughters() const {return daughters_;}


                inline const LauKinematics* getKinematics() const {return kinematics_;}


                inline LauKinematics* getKinematics() {return kinematics_;}


                inline const LauAbsEffModel* getEffModel() const {return effModel_;}


                inline Bool_t usingScfModel() const { return ! scfFractionModel_.empty(); }


                inline const std::vector<LauParameter>& getExtraParameters() const {return extraParameters_;}


                inline std::vector<LauParameter*>& getFloatingParameters() {return resonancePars_;}

        protected:
                void initSummary();

                void initialiseVectors();

                void resetNormVectors();

                void calcDPNormalisation();


                std::vector< std::pair<Double_t,Double_t> > formGapsFromRegions(const std::vector< std::pair<Double_t,Double_t> >& regions, const Double_t min, const Double_t max) const;


                void cullNullRegions(std::vector<LauDPPartialIntegralInfo*>& regions) const;


                LauDPPartialIntegralInfo* newDPIntegrationRegion(const Double_t minm13, const Double_t maxm13,
                                                                 const Double_t minm23, const Double_t maxm23,
                                                                 const Double_t m13BinWidth, const Double_t m23BinWidth,
                                                                 const Double_t precision,
                                                                 const UInt_t nAmp,
                                                                 const UInt_t nIncohAmp) const;


                void correctDPOverlap(std::vector< std::pair<Double_t,Double_t> >& regions, const std::vector<Double_t>& binnings) const;


                std::vector<LauDPPartialIntegralInfo*> m23IntegrationRegions(const std::vector< std::pair<Double_t,Double_t> >& m13Regions,
                                                                             const std::vector< std::pair<Double_t,Double_t> >& m23Regions,
                                                                             const std::vector<Double_t>& m13Binnings,
                                                                             const std::vector<Double_t>& m23Binnings,
                                                                             const Double_t precision,
                                                                             const Double_t defaultBinning) const;


                std::vector<LauDPPartialIntegralInfo*> m13IntegrationRegions(const std::vector< std::pair<Double_t,Double_t> >& m13Regions,
                                                                             const std::vector< std::pair<Double_t,Double_t> >& m23Regions,
                                                                             const std::vector<Double_t>& m13Binnings,
                                                                             const Double_t precision,
                                                                             const Double_t defaultBinning) const;

                void calcDPNormalisationScheme();

                void findIntegralsToBeRecalculated();


                void calcDPPartialIntegral(LauDPPartialIntegralInfo* intInfo);

                void writeIntegralsFile();


                void setFFTerm(const UInt_t index, const Double_t realPart, const Double_t imagPart);


                void setIncohIntenTerm(const UInt_t index, const Double_t value);

                void calculateAmplitudes();


                void calculateAmplitudes( LauDPPartialIntegralInfo* intInfo, const UInt_t m13Point, const UInt_t m23Point );


                void addGridPointToIntegrals(const Double_t weight);


                void calcTotalAmp(const Bool_t useEff);


                Double_t retrieveEfficiency();


                Double_t retrieveScfFraction(Int_t tagCat);


                inline void setASqMaxVarValue(Double_t value) {aSqMaxVar_ = value;}


                Double_t calcSigDPNorm();


                LauComplex resAmp(const UInt_t index);


                Double_t incohResAmp(const UInt_t index);


                void setDataEventNo(UInt_t iEvt);


                LauAbsResonance* findResonance(const TString& resName);


                LauAbsResonance* getResonance(const UInt_t resIndex);


                void removeCharge(TString& string) const;


                Bool_t gotKMatrixMatch(UInt_t resAmpInt, const TString& propName) const;

        private:
                LauIsobarDynamics(const LauIsobarDynamics& rhs);

                LauIsobarDynamics& operator=(const LauIsobarDynamics& rhs);

                typedef std::map<TString, LauKMatrixPropagator*> KMPropMap;

                typedef std::map<TString, TString> KMStringMap;

                LauDaughters* daughters_;

                LauKinematics* kinematics_;

                LauAbsEffModel* effModel_;


                LauTagCatScfFractionModelMap scfFractionModel_;

                UInt_t nAmp_;

                UInt_t nIncohAmp_;

                std::vector<LauComplex> Amp_;

                Double_t DPNorm_;

                LauParArray fitFrac_;

                LauParArray fitFracEffUnCorr_;

                LauParameter DPRate_;

                LauParameter meanDPEff_;

                std::vector<LauCacheData*> data_;

                LauCacheData* currentEvent_;

                std::vector<LauParameter> extraParameters_;

                std::vector<LauAbsResonance*> sigResonances_;

                std::vector<LauAbsIncohRes*> sigIncohResonances_;

                KMPropMap kMatrixPropagators_;

                KMStringMap kMatrixPropSet_;

                std::vector<TString> resTypAmp_;

                std::vector<Int_t> resPairAmp_;

                std::vector<TString> incohResTypAmp_;

                std::vector<Int_t> incohResPairAmp_;

                std::vector<Int_t> typDaug_;

                Bool_t symmetricalDP_;

                Bool_t fullySymmetricDP_;

                Bool_t flavConjDP_;

                Bool_t integralsDone_;

                Bool_t normalizationSchemeDone_;

                Bool_t forceSymmetriseIntegration_;

                std::vector<LauDPPartialIntegralInfo*> dpPartialIntegralInfo_;

                TString intFileName_;

                Double_t m13BinWidth_;

                Double_t m23BinWidth_;

                Double_t mPrimeBinWidth_;

                Double_t thPrimeBinWidth_;

                Double_t narrowWidth_;

                Double_t binningFactor_;

                Double_t m13Sq_;

                Double_t m23Sq_;

                Double_t mPrime_;

                Double_t thPrime_;

                Int_t tagCat_;

                Double_t eff_;

                Double_t scfFraction_;

                Double_t jacobian_;

                Double_t ASq_;

                Double_t evtLike_;

                LauComplex totAmp_;


                std::vector< std::vector<LauComplex> > fifjEffSum_;


                std::vector< std::vector<LauComplex> > fifjSum_;

                std::vector<LauComplex> ff_;

                std::vector<Double_t> incohInten_;

                std::vector<Double_t> fSqSum_;

                std::vector<Double_t> fSqEffSum_;

                std::vector<Double_t> fNorm_;

                Int_t iterationsMax_;

                Int_t nSigGenLoop_;

                Double_t aSqMaxSet_;

                Double_t aSqMaxVar_;

                Bool_t flipHelicity_;

                Bool_t recalcNormalisation_;

                std::vector<LauParameter*> resonancePars_;

                std::vector<Double_t> resonanceParValues_;

                std::vector< std::vector<UInt_t> > resonanceParResIndex_;

                std::set<UInt_t> integralsToBeCalculated_;

                ClassDef(LauIsobarDynamics,0)
};

#endif