Version 1.2, beta now means beta_s in cylindrical blast wave

CMakeLists.txt

@@ -7,8 +7,8 @@ project (ThermalFIST)
 
 # The version number.
 set (ThermalFIST_VERSION_MAJOR 1)
-set (ThermalFIST_VERSION_MINOR 1)
-set (ThermalFIST_VERSION_DEVEL 6)
+set (ThermalFIST_VERSION_MINOR 2)
+set (ThermalFIST_VERSION_DEVEL 0)
 
 # configure a header file to pass some of the CMake settings
 # to the source code

include/HRGEventGenerator/CylindricalBlastWaveEventGenerator.h

@@ -29,31 +29,40 @@ namespace thermalfist {
      * \param TPS    A pointer to the particle list
      * \param config Event generator configuration
      * \param Tkin   The kinetic freeze-out temperature (in GeV) 
-     * \param beta   The transverse flow velocity
+     * \param betas  The transverse flow velocity at the surface
      * \param etamax The longitudinal space-time rapidity cut-off
      * \param npow   The power in the transverse flow profile function
      */
     CylindricalBlastWaveEventGenerator(ThermalParticleSystem *TPS, 
                                        const EventGeneratorConfiguration& config, 
                                        double T = 0.120, 
-                                       double beta = 0.5, 
+                                       double betas = 0.5, 
                                        double etamax = 0.5, 
                                        double npow = 1.);
 
     /// \deprecated
     /// \brief Old constructor. Included for backward compatibility.
-    CylindricalBlastWaveEventGenerator(ThermalModelBase *THM, double T = 0.120, double beta = 0.5, double etamax = 0.5, double npow = 1., bool onlyStable = false, EventGeneratorConfiguration::ModelType EV = EventGeneratorConfiguration::PointParticle, ThermalModelBase *THMEVVDW = NULL);
+    CylindricalBlastWaveEventGenerator(ThermalModelBase *THM, double T = 0.120, double betas = 0.5, double etamax = 0.5, double npow = 1., bool onlyStable = false, EventGeneratorConfiguration::ModelType EV = EventGeneratorConfiguration::PointParticle, ThermalModelBase *THMEVVDW = NULL);
 
     ~CylindricalBlastWaveEventGenerator() { }
 
     /// Sets the momentum distribution parameters
-    void SetParameters(double T, double beta, double etamax, double npow = 1.);
+    void SetParameters(double T, double betas, double etamax, double npow = 1.);
+
+    /**
+     * \brief Set the mean transverse flow velocity.
+     *
+     * Surface velocity is calculated as \\beta_s = \\langle \\beta_T \\rangle (2+n)/2 
+     * 
+     * \param betaT  The mean transverse flow velocity
+     */
+    void SetMeanBetaT(double betaT);
 
     /// Sets up the random generators of particle momenta
     /// and resonances masses
     void SetMomentumGenerators();
   private:
-    double m_T, m_Beta, m_EtaMax, m_n;
+    double m_T, m_BetaS, m_EtaMax, m_n;
   };
 
   /// For backward compatibility

include/HRGEventGenerator/MomentumDistribution.h

@@ -151,14 +151,14 @@ namespace thermalfist {
     /**
      * \copydoc MomentumDistributionBase()
      * \param T      The kinetic temperature (in GeV)
-     * \param beta   The transverse flow velocity
+     * \param beta   The transverse flow velocity at the surface
      * \param etamax The longitudinal space-time rapidity cut-off
      * \param npow   The power in the transverse flow profile function
      * \param norm   Whether the momentum distribution should be normalized to unity 
      */
-    SSHDistribution(int pdgid = 0, double mass = 0., double T = 0.100, double beta = 0.5, double etamax = 0.5, double npow = 1., bool norm = false) :
+    SSHDistribution(int pdgid = 0, double mass = 0., double T = 0.100, double betas = 0.5, double etamax = 0.5, double npow = 1., bool norm = false) :
       MomentumDistributionBase(pdgid, mass),
-      m_T(T), m_Beta(beta), m_EtaMax(etamax), m_n(npow)
+      m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow)
     {
       m_NormY = m_NormPt = m_Norm = 1.;
       if (norm) Normalize();
@@ -172,16 +172,16 @@ namespace thermalfist {
      * \brief Set the parameters of the longitudinal blast-wave distribution
      * 
      * \param T      The kinetic temperature (in GeV)
-     * \param beta   The transverse flow velocity
+     * \param betas  The transverse flow velocity at the surface
      * \param etamax The longitudinal space-time rapidity cut-off
      * \param npow   The power in the transverse flow profile function
      * \param mass   Particle mass (in GeV)
      * \param pdgid  Particle PDG code
      * \param norm   Whether the momentum distribution should be normalized to unity 
      */
-    void SetParameters(double T, double beta, double etamax, double npow, double mass, int pdgid = 0, bool norm = true) {
+    void SetParameters(double T, double betas, double etamax, double npow, double mass, int pdgid = 0, bool norm = true) {
       m_T = T;
-      m_Beta = beta;
+      m_BetaS = betas;
       m_EtaMax = etamax;
       m_n = npow;
       m_Mass = mass;
@@ -192,6 +192,17 @@ namespace thermalfist {
       else Initialize();
     }
 
+    /**
+     * \brief Set the mean transverse flow velocity.
+     *
+     * Surface velocity is calculated as \\beta_s = \\langle \\beta_T \\rangle (2+n)/2 
+     * 
+     * \param betaT  The mean transverse flow velocity
+     */
+    void SetMeanBetaT(double betaT) {
+      m_BetaS = (2. + m_n) / 2. * betaT;
+    }
+
     void Normalize();
 
     /// Rapidity distribution at fixed pT
@@ -239,11 +250,11 @@ namespace thermalfist {
 
     double betar(double r) const {
       if (m_n == 1.)
-        return m_Beta * r;
+        return m_BetaS * r;
       else if (m_n == 2.)
-        return m_Beta * r * m_Beta * r;
+        return m_BetaS * r * r;
       else
-        return pow(m_Beta * r, m_n);
+        return m_BetaS * pow(r, m_n);
     }
 
     double rho(double r) const { return atanh(betar(r)); }
@@ -253,7 +264,7 @@ namespace thermalfist {
     double y2Av() const;
 
     double m_T;
-    double m_Beta, m_EtaMax;
+    double m_BetaS, m_EtaMax;
     double m_NormY, m_NormPt, m_Norm;
     double m_n;
     std::vector<double> m_xlag, m_wlag;

include/HRGEventGenerator/RandomGenerators.h

@@ -204,13 +204,13 @@ namespace thermalfist {
       /**
        * \brief Construct a new SSHGenerator object
        * \param T      The kinetic temperature (in GeV)
-       * \param beta   The transverse flow velocity
+       * \param betas  The transverse flow velocity at the surface
        * \param etamax The longitudinal space-time rapidity cut-off
        * \param npow   The power in the transverse flow profile function
        * \param mass   Particle mass (in GeV) 
        */
-      SSHMomentumGenerator(double T, double beta, double etamax, double npow, double mass) :m_T(T), m_Beta(beta), m_EtaMax(etamax), m_n(npow), m_Mass(mass) {
-        m_distr = SSHDistribution(0, m_Mass, m_T, m_Beta, m_EtaMax, m_n, false);
+      SSHMomentumGenerator(double T, double betas, double etamax, double npow, double mass) :m_T(T), m_BetaS(betas), m_EtaMax(etamax), m_n(npow), m_Mass(mass) {
+        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n, false);
         m_dPt = 0.02;
         m_dy = 0.05;
         FixParameters2();
@@ -223,18 +223,33 @@ namespace thermalfist {
        * \brief Sets the parameters of the distribution
        * 
        * \param T      The kinetic temperature (in GeV)
-       * \param beta   The transverse flow velocity
+       * \param betas  The transverse flow velocity at the surface
        * \param etamax The longitudinal space-time rapidity cut-off
        * \param npow   The power in the transverse flow profile function
        * \param mass   Particle mass (in GeV) 
        */
-      void SetParameters(double T, double beta, double etamax, double npow, double mass) {
+      void SetParameters(double T, double betas, double etamax, double npow, double mass) {
         m_T = T;
-        m_Beta = beta;
+        m_BetaS = betas;
         m_EtaMax = etamax;
         m_Mass = mass;
         m_n = npow;
-        m_distr = SSHDistribution(0, m_Mass, m_T, m_Beta, m_EtaMax, m_n);
+        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n);
+        m_dPt = 0.02;
+        m_dy = 0.05;
+        FixParameters2();
+      }
+
+      /**
+        * \brief Set the mean transverse flow velocity.
+        *
+        * Surface velocity is calculated as \\beta_s = \\langle \\beta_T \\rangle (2+n)/2 
+        * 
+        * \param betaT  The mean transverse flow velocity
+        */
+      void SetMeanBetaT(double betaT) {
+        m_BetaS = (2. + m_n) / 2. * betaT;
+        m_distr = SSHDistribution(0, m_Mass, m_T, m_BetaS, m_EtaMax, m_n);
         m_dPt = 0.02;
         m_dy = 0.05;
         FixParameters2();
@@ -275,7 +290,7 @@ namespace thermalfist {
 
       std::pair<double, double> GetRandom2() const;
 
-      double m_T, m_Beta, m_EtaMax, m_n, m_Mass;
+      double m_T, m_BetaS, m_EtaMax, m_n, m_Mass;
       double m_MaxY, m_MaxPt;
       SSHDistribution m_distr;
       SplineFunction m_dndpt;

src/gui/QtThermalFIST/SpectralFunctionDialog.cpp

@@ -111,6 +111,8 @@ SpectralFunctionDialog::SpectralFunctionDialog(QWidget *parent, ThermalParticle
   setWindowTitle(QString(particle->Name().c_str()) + tr(" spectral function"));
 
   calculate();
+
+  cpath = QApplication::applicationDirPath();
 }
 
 void SpectralFunctionDialog::replot()
@@ -185,7 +187,9 @@ void SpectralFunctionDialog::replotSpectralFunctions()
     plot->graph(1)->addData(BWm[i], mnozh * BWTHval[i]);
   }
 
-
+  plot->setContextMenuPolicy(Qt::CustomContextMenu);
+  connect(plot, SIGNAL(customContextMenuRequested(QPoint)), this, SLOT(contextMenuRequest(QPoint)));
+
 
   plot->replot();
 }
@@ -476,4 +480,113 @@ void SpectralFunctionDialog::calculate()
   replot();
 }
 
+void SpectralFunctionDialog::contextMenuRequest(QPoint pos)
+{
+  QMenu *menu = new QMenu(this);
+  menu->setAttribute(Qt::WA_DeleteOnClose);
+
+  menu->addAction("Save as pdf", this, SLOT(saveAsPdf()));
+  menu->addAction("Save as png", this, SLOT(saveAsPng()));
+  menu->addAction("Write computed values to file", this, SLOT(saveAsAscii()));
+
+  menu->popup(plot->mapToGlobal(pos));
+}
+
+void SpectralFunctionDialog::saveAsPdf()
+{
+  saveAs(0);
+}
+
+void SpectralFunctionDialog::saveAsPng()
+{
+  saveAs(1);
+}
+
+void SpectralFunctionDialog::saveAsAscii()
+{
+  saveAs(2);
+}
+
+void SpectralFunctionDialog::saveAs(int type)
+{
+  QString tname = QString::fromStdString(particle->Name());
+  if (comboView->currentIndex() == 0)
+    tname += "_spectral";
+  else if (comboView->currentIndex() == 1)
+    tname += "_BRs";
+  else
+    tname += "_widths";
+  for (int i = 0; i < tname.size(); ++i) {
+    if (tname[i] == '/')
+      tname[i] = '.';
+  }
+
+  QVector<QString> exts;
+  exts.push_back("pdf");
+  exts.push_back("png");
+  exts.push_back("dat");
+
+  QString listpathprefix = cpath + "/" + tname + "." + exts[type];
+  QString path = QFileDialog::getSaveFileName(this, "Save plot as " + exts[type], listpathprefix, "*." + exts[type]);
+  if (path.length()>0)
+  {
+    if (type == 0)
+      plot->savePdf(path, plot->width(), plot->height());
+    else if (type == 1)
+      plot->savePng(path, plot->width(), plot->height());
+    else {
+      QFile fout(path);
+
+      if (fout.open(QIODevice::WriteOnly | QIODevice::Text)) {
+        QTextStream out(&fout);
+        if (comboView->currentIndex() == 0) {
+          out.setFieldWidth(20);
+          out.setFieldAlignment(QTextStream::AlignLeft);
+          out << "M[GeV]";
+          out << "rho[GeV^-1]_vac";
+          out << "rho[GeV^-1]_th";
+          out << qSetFieldWidth(0) << endl << qSetFieldWidth(20);
+          for (int i = 0; i < plot->graph(0)->dataCount(); ++i) {
+            out.setFieldWidth(20);
+            out.setFieldAlignment(QTextStream::AlignLeft);
+            out << plot->graph(0)->dataMainKey(i);
+            out << plot->graph(0)->dataMainValue(i);
+            out << plot->graph(1)->dataMainValue(i);
+            out << qSetFieldWidth(0) << endl << qSetFieldWidth(20);
+          }
+        }
+        else  {
+          out.setFieldWidth(30);
+          out.setFieldAlignment(QTextStream::AlignLeft);
+          out << "M[GeV]";
+          for (int ch = 0; ch < plot->graphCount(); ++ch) {
+            QString tnamech = "";
+            const QString& name = plot->graph(ch)->name();
+            for (int ic = 0; ic < name.size(); ++ic)
+              if (name[ic] != QChar(' '))
+                tnamech += name[ic];
+            if (comboView->currentIndex() == 1)
+              out << "BR_" + tnamech;
+            else
+              out << "Gamm_" + tnamech + "[GeV]";
+          }
+          out << qSetFieldWidth(0) << endl << qSetFieldWidth(30);
+          for (int i = 0; i < plot->graph(0)->dataCount(); ++i) {
+            out.setFieldWidth(30);
+            out.setFieldAlignment(QTextStream::AlignLeft);
+            out << plot->graph(0)->dataMainKey(i);
+            for (int ch = 0; ch < plot->graphCount(); ++ch) {
+              out << plot->graph(ch)->dataMainValue(i);
+            }
+            out << qSetFieldWidth(0) << endl << qSetFieldWidth(30);
+          }
+        }
+      }
+    }
+
+    QFileInfo saved(path);
+    cpath = saved.absolutePath();
+  }
+}
+
 

src/gui/QtThermalFIST/SpectralFunctionDialog.h

@@ -63,6 +63,8 @@ class SpectralFunctionDialog : public QDialog
     //QCPColorMap *colormap;
     //QCPColorScale *colorScale;
 
+    QString cpath;
+
     QComboBox *comboScheme;
 
     QDoubleSpinBox *spinT;
@@ -86,6 +88,13 @@ public slots:
     void replotSpectralFunctions();
     void replotBranchingRatios();
     void replotPartialWidths();
+
+    void contextMenuRequest(QPoint pos);
+    void saveAsPdf();
+    void saveAsPng();
+    void saveAsAscii();
+    // saveAs type: 0 - pdf, 1 - png, 2 - ascii data
+    void saveAs(int type);
 };
 
 #endif // SPECTRALFUNCTIONDIALOG_H

src/gui/QtThermalFIST/eventgeneratortab.cpp

@@ -374,7 +374,7 @@ EventGeneratorTab::EventGeneratorTab(QWidget *parent, ThermalModelBase *modelop)
     spinBeta->setDecimals(3);
     spinBeta->setValue(0.5);
     spinBeta->setToolTip(tr("Radial flow velocity"));
-    QLabel *labelBetat = new QLabel(tr("β<sub>T</sub>:"));
+    QLabel *labelBetat = new QLabel(tr("⟨β<sub>T</sub>⟩:"));
     spinBetat = new QDoubleSpinBox();
     spinBetat->setMinimum(0.);
     spinBetat->setMaximum(1.);
@@ -1067,10 +1067,11 @@ void EventGeneratorTab::generateEvents(const ThermalModelConfig & config)
 
     model->CalculateDensitiesGCE();
 
-
+    // Convert the mean transverse velocity into the one at the surface
+    double betaS = (2. + spinn->value()) / 2. * spinBetat->value();
 
     if (radioSR->isChecked()) spectra->Reset(model, spinTkin->value() * 1.e-3, spinBeta->value());
-    else if (radioSSH->isChecked()) spectra->Reset(model, spinTkin->value() * 1.e-3, spinBetat->value(), 1, spinEtaMax->value(), spinn->value());
+    else if (radioSSH->isChecked()) spectra->Reset(model, spinTkin->value() * 1.e-3, betaS, 1, spinEtaMax->value(), spinn->value());
 
     int id = getCurrentRow();
     if (id<0) id = 0;
@@ -1127,7 +1128,7 @@ void EventGeneratorTab::generateEvents(const ThermalModelConfig & config)
     if (radioSR->isChecked())
       generator = new SphericalBlastWaveEventGenerator(model->TPS(), configMC, spinTkin->value() * 1.e-3, spinBeta->value());
     else 
-      generator = new CylindricalBlastWaveEventGenerator(model->TPS(), configMC, spinTkin->value() * 1.e-3, spinBetat->value(), spinEtaMax->value(), spinn->value());
+      generator = new CylindricalBlastWaveEventGenerator(model->TPS(), configMC, spinTkin->value() * 1.e-3, betaS, spinEtaMax->value(), spinn->value());
 
 
     //if (radioSR->isChecked()) generator = new SphericalBlastWaveEventGenerator(model, spinTkin->value() * 1.e-3, spinBeta->value(), false, EV, modelEVVDW);

src/library/HRGBase/Broyden.cpp

@@ -45,6 +45,7 @@ namespace thermalfist {
     for (size_t i = 0; i < x.size(); ++i) {
       h[i] = m_dx*abs(h[i]);
       if (h[i] == 0.0) h[i] = m_dx;
+      if (h[i] < 1.e-10) h[i] = 1.e-10;
       //h[i] = max(m_dx, h[i]);
     }