BoltzmannDistribution.cxx

Go to the documentation of this file.

#include "BoltzmannDistribution.h"
#include <iostream>
 
namespace BoltzmannDistributionNamespace {
 
#include "../NumericalIntegration.h"
 
}
 
BoltzmannDistribution::BoltzmannDistribution(double mass,
                                             int PDGID,
                                             bool UseAcc,
                                             double ymin,
                                             double ymax,
                                             double ycm,
                                             double width)
  : xlag32()
  , wlag32()
  , xleg32()
  , wleg32()
  , fMass(mass)
  , fPDGID(PDGID)
  , fUseAcceptance(UseAcc)
  , fYmin(ymin)
  , fYmax(ymax)
  , fYcm(ycm)
  , fWidth(width)
  , fAcceptance()
  , fReconstructionEfficiency()
  , fTamt()
  , fNormT()
  , fNormT4pi()
//              fNorm(1.), fAcceptance(), fReconstructionEfficiency(), fTamt(), fNormT(), fNormT4pi()
{
  //InverseSlopeNamespace::GetCoefsIntegrateLaguerre32(xlag32, wlag32);
  fWidth = 1.;
  BoltzmannDistributionNamespace::GetCoefsIntegrateLegendre32(
    0., 3., xlag32, wlag32);
  BoltzmannDistributionNamespace::GetCoefsIntegrateLegendre32(
    fYmin, fYmax, xleg32, wleg32);
  //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;
      //std::cout << dir + "pty_acc_" + spdg + ".txt" << "\nAcceptance not OK!\n" << "\n";
      return;
    }
    //std::cout << "Acceptance OK!\n" << "\n";
    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;
    //std::cout << ymin << "<" << "y" << "<" << ymax << "\n";
    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));
      //std::cout << "Calc: " << Tvec[i] << " " << mtacc[i] << " " << norm[i] << "\n";
    }
    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());
    //for(int i=0;i<Tvec.size();++i) {
    //  std::cout << "ACalc: " << fTamt.Eval(mtacc[i]) << " " << mtacc[i] << " " << fNormT.Eval(Tvec[i]) << "\n";
    //  std::cout << "Avs: " << mtAv(Tvec[i]) << " " << mtAv2(Tvec[i]) << "\n";
    //}
  }
}
 
double BoltzmannDistribution::mtAv(double T) {
  double ret1 = 0., ret2 = 0.;
  for (Int_t i = 0; i < 32; i++) {
    for (Int_t j = 0; j < 32; j++) {
      double tmpf = 0.;
      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 = wlag32[i] * wleg32[j] * xlag32[i] * tmt * TMath::CosH(xleg32[j])
             * TMath::Exp(-tmt * TMath::CosH(xleg32[j]) / T);
      if (fUseAcceptance)
        tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i])
                * fReconstructionEfficiency.f(tp);
      ret1 += tmpf;
      ret2 += tmpf * tmt;
    }
  }
  return ret2 / ret1;
}
 
double BoltzmannDistribution::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 BoltzmannDistribution::fmt(double mt, double T) {
  double ret = 0.;
  double tpt = sqrt(mt * mt - fMass * fMass);
  for (Int_t j = 0; j < 32; j++) {
    double tmpf = 0.;
    double tp   = mt * TMath::CosH(xleg32[j] + fYcm);
    tp          = sqrt(tp * tp - fMass * fMass);
    tmpf        = wleg32[j] * tpt * mt * TMath::CosH(xleg32[j])
           * TMath::Exp(-mt * TMath::CosH(xleg32[j]) / T);
    if (fUseAcceptance)
      tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, tpt)
              * fReconstructionEfficiency.f(tp);
    ret += tmpf;
  }
  return ret * mt / tpt / fNormT.Eval(T);
}
 
double BoltzmannDistribution::Normalization(double T) {
  double ret = 0.;
  //if (fUseAcceptance)
  {
    for (Int_t i = 0; i < 32; i++) {
      for (Int_t j = 0; j < 32; j++) {
        double tmpf = 0.;
        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 = wlag32[i] * wleg32[j] * xlag32[i] * tmt * TMath::CosH(xleg32[j])
               * TMath::Exp(-tmt * TMath::CosH(xleg32[j]) / T);
        if (fUseAcceptance)
          tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i])
                  * fReconstructionEfficiency.f(tp);
        ret += tmpf;
      }
    }
  }
  return ret;
}
 
double BoltzmannDistribution::Normalization4pi(double T) {
  double ret = 0.;
  //if (fUseAcceptance)
  {
    for (Int_t i = 0; i < 32; i++) {
      for (Int_t j = 0; j < 32; j++) {
        double tmpf = 0.;
        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 = wlag32[i] * wleg32[j] * xlag32[i] * tmt * TMath::CosH(xleg32[j])
               * TMath::Exp(-tmt * TMath::CosH(xleg32[j]) / T);
        ret += tmpf;
      }
    }
  }
  return ret;
}
 
double BoltzmannDistribution::dndy(double y, double A, double T) {
  double ret = 0.;
  for (Int_t i = 0; i < 32; i++) {
    //for(Int_t j = 0 ; j < 32 ; j++){
    double tmpf = 0.;
    double tp =
      sqrt(xlag32[i] * xlag32[i] + fMass * fMass) * TMath::CosH(y + fYcm);
    tp         = sqrt(tp * tp - fMass * fMass);
    double tmt = sqrt(xlag32[i] * xlag32[i] + fMass * fMass);
    tmpf       = wlag32[i] * xlag32[i] * tmt * TMath::CosH(y)
           * TMath::Exp(-tmt * TMath::CosH(y) / T);
    ret += tmpf;
    //}
  }
  return ret / fNormT4pi.Eval(T) * A;
}