ThermalParticleSystem.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2014-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef THERMALPARTICLESYSTEM_H
#define THERMALPARTICLESYSTEM_H

#include <map>
#include <vector>
#include <set>
#include <fstream>

#include "HRGBase/ThermalParticle.h"

namespace thermalfist {

  class ThermalParticleSystem
  {
  public:
    
    ThermalParticleSystem(const std::string& InputFile = "", bool GenAntiP = true, double mcut = -1.) { Initialize(InputFile, std::string(""), GenAntiP, mcut); }

    ThermalParticleSystem(const std::string& InputFile, const std::string& DecayFile, bool GenAntiP = true, double mcut = -1.) { Initialize(InputFile, DecayFile, GenAntiP, mcut); }


    ThermalParticleSystem(
      const std::vector<std::string>& ListFiles,
      const std::vector<std::string>& DecayFiles = std::vector<std::string>(0),
      const std::set<std::string>& flags = std::set<std::string>(),
      double mcut = -1.);

    ~ThermalParticleSystem(void);

    ThermalParticle::ParticleDecaysVector GetDecaysFromAntiParticle(const ThermalParticle::ParticleDecaysVector& Decays);

    void ProcessDecays();

    void FillDecayProperties();
    
    void FillDecayThresholds();

    void FillResonanceDecays();

    void FillResonanceDecaysByFeeddown();


    typedef std::pair<double, int> SingleDecayContribution;

    typedef std::vector<SingleDecayContribution> DecayContributionsToParticle;

    typedef std::vector<DecayContributionsToParticle> DecayContributionsToAllParticles;

    typedef std::pair< std::vector<double>, int> SingleDecayCumulantsContribution;

    typedef std::vector<SingleDecayCumulantsContribution> DecayCumulantsContributionsToParticle;

    typedef std::vector<DecayCumulantsContributionsToParticle> DecayCumulantsContributionsToAllParticles;

    typedef std::pair< std::vector<double>, int> SingleDecayProbabilityDistribution;

    typedef std::vector<SingleDecayProbabilityDistribution> DecayProbabilityDistributionsToParticle;

    typedef std::vector<DecayProbabilityDistributionsToParticle> DecayProbabilityDistributionsToAllParticles;

    typedef std::vector< std::pair<double, std::vector<int> > > ResonanceFinalStatesDistribution;

    const std::vector<DecayContributionsToAllParticles>& DecayContributionsByFeeddown() const { return m_DecayContributionsByFeeddown; }

    const DecayCumulantsContributionsToAllParticles& DecayCumulants() const { return m_DecayCumulants; }

    const std::vector<ResonanceFinalStatesDistribution>& ResonanceFinalStatesDistributions() const { return m_ResonanceFinalStatesDistributions; }

    void LoadList(const std::vector<std::string>& ListFiles,
      const std::vector<std::string>& DecayFiles = std::vector<std::string>(0), 
      const std::set<std::string>& flags = std::set<std::string>(),
      double mcut = 1.e9);

    void LoadList(const std::string& InputFile, const std::string& DecayFile, bool GenerateAntiParticles = true, double mcut = 1.e9);

    void LoadList(const std::string& InputFile = "", bool GenerateAntiParticles = true, double mcut = -1.) { LoadList(InputFile, std::string(""), GenerateAntiParticles, mcut); }

    void AddParticlesToListFromFile(const std::string& InputFile = "", const std::set<std::string>& flags = std::set<std::string>(), double mcut = -1.);


    void LoadTable(const std::string& InputFile = "", bool GenerateAntiParticles = true, double mcut = -1.) { LoadList(InputFile, GenerateAntiParticles, mcut); }
    
    void SetTableFromVector(const std::vector<ThermalParticle> &part_in, bool GenerateAntiParticles = true);
    
    void WriteTableToFile(const std::string& OutputFile = "", bool WriteAntiParticles = false);

    void LoadDecays(const std::vector<std::string>& DecayFiles, const std::set<std::string>& flags = std::set<std::string>());

    void LoadDecays(const std::string& DecaysFile = "", bool GenerateAntiParticles = true);
    
    void WriteDecaysToFile(const std::string& OutputFile = "", bool WriteAntiParticles = false);

    void NormalizeBranchingRatios();

    void RestoreBranchingRatios();


    void SetCalculationType(IdealGasFunctions::QStatsCalculationType type);
    IdealGasFunctions::QStatsCalculationType QStatsCalculationType() const { return m_QStatsCalculationType; }

    void SetClusterExpansionOrder(int order);


    void SetResonanceWidthShape(ThermalParticle::ResonanceWidthShape shape);
    ThermalParticle::ResonanceWidthShape ResonanceWidthShape() const { return m_ResonanceWidthShape; }


    void SetResonanceWidthIntegrationType(ThermalParticle::ResonanceWidthIntegration type);
    ThermalParticle::ResonanceWidthIntegration ResonanceWidthIntegrationType() const { return m_ResonanceWidthIntegrationType; }

    std::string GetNameFromPDG(long long pdgid);

    bool hasBaryons() const { return (m_NumBaryons > 0); }

    bool hasCharged() const { return (m_NumCharged > 0); }

    bool hasStrange() const { return (m_NumStrange > 0); }

    bool hasCharmed() const { return (m_NumCharmed > 0); }

    int ComponentsNumber() const { return static_cast<int>(m_Particles.size()); }

    const std::vector<ThermalParticle>& Particles() const { return m_Particles; }


    const ThermalParticle& Particle(int id) const;
    ThermalParticle& Particle(int id);

    ThermalParticle& ParticleByPDG(long long pdgid);

    int  PdgToId(long long pdgid)    /*const*/ { return (m_PDGtoID.count(pdgid) > 0) ? m_PDGtoID[pdgid] : -1; }
    
    long long IdToPdg(int id)      const { return (id >= 0 && id < static_cast<int>(m_Particles.size())) ? m_Particles[id].PdgId() : 0; }

    void AddParticle(const ThermalParticle & part);

    void RemoveParticleAt(int ind);

    bool CheckDecayChargesConservation(int ind) const;

    bool CheckDecayChannelsAreSpecified() const;

    std::vector<int> CheckDecayChargesConservationVector(int ind) const;

    bool operator==(const ThermalParticleSystem &rhs) const;
    bool operator!=(const ThermalParticleSystem &rhs) const { return !(*this == rhs); }

    void FillPdgMap();

    enum SortModeType {
      SortByMass = 0,
      SortByMassAndPDG = 1,
      SortByBaryonAndMassAndPDG = 2
    };

    void SetSortMode(SortModeType type) { m_SortMode = type; FinalizeList(); ProcessDecays(); }

    SortModeType SortMode() const { return m_SortMode; }

    void FinalizeList();

    static ParticleDecayType::DecayType DecayTypeByParticleType(const ThermalParticle &part);

    static const std::string flag_no_antiparticles;
    static const std::string flag_nostrangeness;
    static const std::string flag_nocharm;
    static const std::string flag_nonuclei;
    static const std::string flag_noexcitednuclei;

  private:
    void GoResonance(int ind, int startind, double BR);

    void GoResonanceByFeeddown(int ind, int startind, double BR, Feeddown::Type feeddown);

    std::vector<double> GoResonanceDecayProbs(int ind, int goalind, bool firstdecay = false);

    std::vector<double> GoResonanceDecayProbsCharge(int ind, int nch, bool firstdecay = false);

    ResonanceFinalStatesDistribution GoResonanceDecayDistributions(int ind, bool firstdecay = false);

    bool AcceptParticle(const ThermalParticle& part, const std::set<std::string>& flags, double mcut = -1.) const;

    //void LoadTable_OldFormat(std::ifstream &fin, bool GenerateAntiParticles = true, double mcut = 1.e9);
    void LoadTable_OldFormat(std::ifstream& fin, const std::set<std::string>& flags = std::set<std::string>(), double mcut = 1.e9);

    //void LoadTable_NewFormat(std::ifstream &fin, bool GenerateAntiParticles = true, double mcut = 1.e9);
    void LoadTable_NewFormat(std::ifstream& fin, const std::set<std::string>& flags = std::set<std::string>(), double mcut = 1.e9);

    void ReadDecays_OldFormat(std::ifstream &fin);

    void ReadDecays_NewFormat(std::ifstream &fin);

    void Initialize(const std::vector<std::string>& ListFiles,
      const std::vector<std::string>& DecayFiles = std::vector<std::string>(0),
      const std::set<std::string>& flags = std::set<std::string>(),
      double mcut = -1.);

    void Initialize(const std::string& InputFile = "", const std::string& DecayFile = "", bool GenAntiP = true, double mcut = -1.);


  private:
    std::vector<ThermalParticle>    m_Particles;
    std::map<long long, int>              m_PDGtoID;
    int m_NumBaryons;
    int m_NumCharged;
    int m_NumStrange;
    int m_NumCharmed;

    int m_NumberOfParticles;

    ThermalParticle::ResonanceWidthIntegration m_ResonanceWidthIntegrationType;

    ThermalParticle::ResonanceWidthShape m_ResonanceWidthShape;

    IdealGasFunctions::QStatsCalculationType m_QStatsCalculationType;

    std::vector<DecayContributionsToAllParticles> m_DecayContributionsByFeeddown;

    DecayCumulantsContributionsToAllParticles m_DecayCumulants;

    DecayProbabilityDistributionsToAllParticles m_DecayProbabilities;

    std::vector<ResonanceFinalStatesDistribution> m_ResonanceFinalStatesDistributions;

    // Map for DP-based calculations of decay distributions
    std::vector<ResonanceFinalStatesDistribution> m_DecayDistributionsMap;

    SortModeType m_SortMode;
  };

  namespace CuteHRGHelper {
    std::vector<std::string> split(const std::string &s, char delim);
    void cutDecayDistributionsVector(std::vector<std::pair<double, std::vector<int> > > &vect, int maxsize = 1000);
  }

  namespace ExtraParticles {
    const ThermalParticle& Particle(int id);
    const ThermalParticle& ParticleByPdg(long long pdg);
    int PdgToId(long long pdg);
    bool Init();
    std::string NameByPdg(long long pdg);
  }

} // namespace thermalfist

#endif