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); }

    bool hasMultiBaryons() const { return (m_MaxAbsBaryonNumber > 1); }

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

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


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

    const ThermalParticle& ParticleByPDG(long long pdgid) const;
    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;

    bool CheckAbsoluteQuarkNumbers() 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;

     std::vector<double> GetConservedChargesVector(ConservedCharge::Name charge);
 
  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, const std::set<std::string>& flags = std::set<std::string>(), 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.);
 
    void FinalizeListLoad();
 
    void FinalizeDecaysLoad();

    bool CheckListIsiSS(const std::string &filename);

    void LoadListiSS(const std::string& filename, const std::set<std::string>& flags = std::set<std::string>(), double mcut = 1.e9);
 
  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_MaxAbsBaryonNumber;
 
    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 DecayLifetimes {
    // In units of ctau (fm)
    double GetLifetime(long long pdg);
  }
 
} // namespace thermalfist
 
#endif