LauIsobarDynamics.hh

Go to the documentation of this file.

 
/*
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
*/
 
#ifndef LAU_ISOBAR_DYNAMICS
#define LAU_ISOBAR_DYNAMICS
 
#include "LauAbsResonance.hh"
#include "LauComplex.hh"
 
#include "TString.h"
 
#include <set>
#include <vector>
 
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 );
 
 
    LauKMatrixPropagator* 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;
        aSqMaxAuto_ = kFALSE;
    }
 
 
    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 );
 
 
    void collateResonanceParameters();
 
 
    inline void flipHelicityForCPEigenstates( const 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_; }
 
    inline void calculateRhoOmegaFitFractions( const Bool_t calcFF )
    {
        calculateRhoOmegaFitFractions_ = calcFF;
    }
 
  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 aSqMaxAuto_ { kTRUE };
 
    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_;
 
    Bool_t calculateRhoOmegaFitFractions_;
 
    ClassDef( LauIsobarDynamics, 0 )
};
 
#endif