InverseSlope.cxx

Go to the documentation of this file.

#include "InverseSlope.h"
#include <iostream>
 
namespace InverseSlopeNamespace {
 
#include "../NumericalIntegration.h"
 
}
 
InverseSlope::InverseSlope(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::GetCoefsIntegrateLegendre32(0., 3., xlag32, wlag32);
  InverseSlopeNamespace::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;
      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());
  }
}
 
double InverseSlope::mtAv(double T) {
  if (!fUseAcceptance) {
    return (fMass * fMass + 2. * fMass * T + 2. * T * T) / (fMass + T);
  } else {
    //double ret = 0.;
    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] * TMath::Exp(-tmt / T)
               * TMath::Exp(-xleg32[j] * xleg32[j] / 2. / fWidth / fWidth)
               * fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i])
               * fReconstructionEfficiency.f(tp);
        ret1 += tmpf;
        ret2 += tmpf * tmt;
      }
    }
    return ret2 / ret1;
  }
}
 
double InverseSlope::mtAv2(double T) {
  if (!fUseAcceptance) {
    return (fMass * fMass + 2. * fMass * T + 2. * T * T) / (fMass + T);
  } else {
    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 InverseSlope::fmt(double mt, double T) {
  if (!fUseAcceptance)
    return TMath::Exp((fMass - mt) / T) * mt / T / (fMass + T);
  else {
    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 * TMath::Exp(-mt / T)
             * TMath::Exp(-xleg32[j] * xleg32[j] / 2. / fWidth / fWidth)
             * fAcceptance.getAcceptance(xleg32[j] + fYcm, tpt)
             * fReconstructionEfficiency.f(tp);
      ret += tmpf;
    }
    return ret * mt / tpt / fNormT.Eval(T);
  }
}
 
double InverseSlope::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] * TMath::Exp(-tmt / T)
               * TMath::Exp(-xleg32[j] * xleg32[j] / 2. / fWidth / fWidth)
               * fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i])
               * fReconstructionEfficiency.f(tp);
        ret += tmpf;
      }
    }
  }
  return ret;
}
 
double InverseSlope::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] * TMath::Exp(-tmt / T)
               * TMath::Exp(-xleg32[j] * xleg32[j] / 2. / fWidth / fWidth);
        ret += tmpf;
      }
    }
  }
  return ret;
}