BlastWave.cxx

Go to the documentation of this file.

#include "BlastWave.h"
#include <iostream>
 
namespace BlastWaveNamespace {
 
#include "../NumericalIntegration.h"
 
}
 
BlastWave::BlastWave(double mass,
                     int PDGID,
                     bool UseAcc,
                     double ymin,
                     double ymax,
                     double ycm,
                     double etam)
  : xlag32()
  , wlag32()
  , xleg32()
  , wleg32()
  , xlegeta()
  , wlegeta()
  , fMass(mass)
  , fPDGID(PDGID)
  , fUseAcceptance(UseAcc)
  , fYmin(ymin)
  , fYmax(ymax)
  , fYcm(ycm)
  , fEtaMax(etam)
  , fAcceptance()
  , fReconstructionEfficiency()
  , fTamt()
  , fNormT()
  , fNormT4pi()
//              fNorm(1.), fAcceptance(), fReconstructionEfficiency(), fTamt(), fNormT(), fNormT4pi()
{
  if (fabs(fEtaMax) < 1e-5) fEtaMax = 1.e-5;
  BlastWaveNamespace::GetCoefsIntegrateLegendre32(0., 3., xlag32, wlag32);
  BlastWaveNamespace::GetCoefsIntegrateLegendre32(fYmin, fYmax, xleg32, wleg32);
  BlastWaveNamespace::GetCoefsIntegrateLegendre32(
    -fEtaMax, fEtaMax, xlegeta, wlegeta);
  //if (fUseAcceptance)
  {
    TString work = getenv("VMCWORKDIR");
    TString dir =
      work + "/KF/KFModelParameters/common/";  //ThermalModel.mT-t.sts_v13.txt";
    char spdg[300];
    sprintf(spdg, "%d", fPDGID);
    if (!Acceptance::ReadAcceptanceFunction(fAcceptance,
                                            dir + "pty_acc_" + spdg + ".txt")) {
      fUseAcceptance = false;
      return;
    }
    fReconstructionEfficiency =
      Acceptance::ReconstructionEfficiencyFunction(0.99, 0.98, 0.135);
    std::vector<double> Tvec(0), mtacc(0), norm(0), norm4pi(0);
    double dT   = 0.002;
    double Tmax = 0.5;
    for (int i = 0;; ++i) {
      double tmpT = 0.01 + dT * (0.5 + i);
      if (tmpT > Tmax) break;
      Tvec.push_back(tmpT);
      mtacc.push_back(mtAv(tmpT));
      norm.push_back(Normalization(tmpT));
      norm4pi.push_back(Normalization4pi(tmpT));
    }
    fTamt     = TSpline3("fTamtacc", &mtacc[0], &Tvec[0], Tvec.size());
    fNormT    = TSpline3("fNormacc", &Tvec[0], &norm[0], Tvec.size());
    fNormT4pi = TSpline3("fNorm4pi", &Tvec[0], &norm4pi[0], Tvec.size());
  }
}
 
double BlastWave::mtAv(double T) {
  double ret1 = 0., ret2 = 0.;
  for (Int_t ie = 0; ie < 32; ie++) {
    for (Int_t i = 0; i < 32; i++) {
      for (Int_t j = 0; j < 32; j++) {
        double tmpf = 0.;
        double ty   = xleg32[j] - xlegeta[ie];
        double tp   = sqrt(xlag32[i] * xlag32[i] + fMass * fMass)
                    * TMath::CosH(xleg32[j] + fYcm);
        tp         = sqrt(tp * tp - fMass * fMass);
        double tmt = sqrt(xlag32[i] * xlag32[i] + fMass * fMass);
        tmpf       = wlegeta[ie] * wlag32[i] * wleg32[j] * xlag32[i] * tmt
               * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / T);
        if (fUseAcceptance)
          tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i])
                  * fReconstructionEfficiency.f(tp);
        ret1 += tmpf;
        ret2 += tmpf * tmt;
      }
    }
  }
  return ret2 / ret1;
}
 
double BlastWave::mtAv2(double T) {
  double ret = 0.;
  for (Int_t i = 0; i < 32; i++) {
    double tmt = sqrt(xlag32[i] * xlag32[i] + fMass * fMass);
    ret += wlag32[i] * xlag32[i] * fmt(tmt, T);
  }
  return ret;
}
 
double BlastWave::fmt(double mt, double T) {
  double ret = 0.;
  double tpt = sqrt(mt * mt - fMass * fMass);
  for (Int_t ie = 0; ie < 32; ie++) {
    for (Int_t j = 0; j < 32; j++) {
      double tmpf = 0.;
      double tp   = mt * TMath::CosH(xleg32[j] + fYcm);
      double ty   = xleg32[j] - xlegeta[ie];
      tp          = sqrt(tp * tp - fMass * fMass);
      tmpf        = wlegeta[ie] * wleg32[j] * tpt * mt * TMath::CosH(ty)
             * TMath::Exp(-mt * TMath::CosH(ty) / T);
      if (fUseAcceptance)
        tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, tpt)
                * fReconstructionEfficiency.f(tp);
      ret += tmpf;
    }
  }
  return ret * mt / tpt / fNormT.Eval(T);
}
 
double BlastWave::Normalization(double T) {
  double ret = 0.;
  {
    for (Int_t ie = 0; ie < 32; ie++) {
      for (Int_t i = 0; i < 32; i++) {
        for (Int_t j = 0; j < 32; j++) {
          double tmpf = 0.;
          double ty   = xleg32[j] - xlegeta[ie];
          double tp   = sqrt(xlag32[i] * xlag32[i] + fMass * fMass)
                      * TMath::CosH(xleg32[j] + fYcm);
          tp         = sqrt(tp * tp - fMass * fMass);
          double tmt = sqrt(xlag32[i] * xlag32[i] + fMass * fMass);
          tmpf       = wlegeta[ie] * wlag32[i] * wleg32[j] * xlag32[i] * tmt
                 * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / T);
          if (fUseAcceptance)
            tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i])
                    * fReconstructionEfficiency.f(tp);
          ret += tmpf;
        }
      }
    }
  }
  return ret;
}
 
double BlastWave::Normalization4pi(double T) {
  double ret = 0.;
  {
    for (Int_t ie = 0; ie < 32; ie++) {
      for (Int_t i = 0; i < 32; i++) {
        for (Int_t j = 0; j < 32; j++) {
          double tmpf = 0.;
          double ty   = xleg32[j] - xlegeta[ie];
          double tp   = sqrt(xlag32[i] * xlag32[i] + fMass * fMass)
                      * TMath::CosH(xleg32[j] + fYcm);
          tp         = sqrt(tp * tp - fMass * fMass);
          double tmt = sqrt(xlag32[i] * xlag32[i] + fMass * fMass);
          tmpf       = wlegeta[ie] * wlag32[i] * wleg32[j] * xlag32[i] * tmt
                 * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / T);
          ret += tmpf;
        }
      }
    }
  }
  return ret;
}