Zfluk.cxx

Go to the documentation of this file.

// $Author: zarko $ 
//
// $Revision: 1.10 $
// 
// $Name:  $
//
// $Id: Zfluk.cxx,v 1.10 2007/03/03 21:27:05 zarko Exp $
//
// $Log: Zfluk.cxx,v $
// Revision 1.10  2007/03/03 21:27:05  zarko
// Version used for PiMinus fits.
//
// Revision 1.9  2007/02/14 22:58:07  zarko
// In this version of zfluk GetWeight(ptype,pt,pz) returns a weight using
// 16 parameters (6 for pi+, 6 for K+, 2 for pi- and 2 for K-)
//
// Revision 1.8  2006/11/10 19:13:17  zarko
// Changed the way kaon weights are calculated so that K/pi ratio is
// preserved. Added K0L reweighting based on K+ and K- weights.
// To get these weights use SetParameters(vector<double>) and
// GetWeight(ptype, pT, pz) functions.
//
// Revision 1.7  2006/10/03 16:00:56  zarko
// Updated code to reweight numus and numubars
//
// Revision 1.6  2006/06/02 04:37:53  zarko
// Use arrays instead of histograms in GetPTshift to speed up the code.
//
// Revision 1.5  2006/04/14 02:37:31  zarko
// Zfluk now looks into SRT_PUBLIC_CONTEXT if the data cannot be found in SRT_PRIVATE_CONTEXT
//
// Revision 1.4  2006/03/13 23:58:32  zarko
// Updated SKZP parameterization. Now using 7 parameters.
// One for kaons and 6 to parameterize pt-xf weights.
//
// Revision 1.3  2006/02/21 00:04:24  zarko
// Modified Files:
//      Zfluk.cxx Zfluk.h
//
// Added function GetPTshift(par). This uses the fluka05ptxf.root file that should
// be in MCReweight/data. This function returns shift in mean pT for a given set
// of parameters, for a parameterization that was set with UseParameterization.
//
// Revision 1.2  2006/01/30 20:57:56  zarko
// Removed one parameter from SKZP (par[2])
// Redefined par[3] (low pT) and par[4] (kaons) in SKZP,
// and par[5] (kaons) in TV, so that they are just numbers.
// They are not multiplicative factors of sigma distortions any more.
//
//

#include "MCReweight/Zfluk.h"
#include "MessageService/MsgService.h"

#include <string>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include "TFile.h"
#include "TH2F.h"
#include "TSystem.h"
ClassImp(Zfluk)

CVSID("$Id: Zfluk.cxx,v 1.10 2007/03/03 21:27:05 zarko Exp $");

Zfluk::Zfluk(const char* dir) 
  :fPar(0)
{ 
  whichParameterization=0; //by default using SKZP parameterization

  // topDir is where the fluka05 pt-xf histogram file is
  std::string topDir=dir; // user may set location of input data
  if(topDir=="") { // by default, this code looks in a standard place
    topDir="MCReweight/data";
    std::string base="";
    base=getenv("SRT_PRIVATE_CONTEXT");
    if (base!="" && base!=".")
      {
        // check if directory exists in SRT_PRIVATE_CONTEXT
        std::string path = base + "/" + topDir;
        void *dir_ptr = gSystem -> OpenDirectory(path.c_str());
        if(!dir_ptr) base=getenv("SRT_PUBLIC_CONTEXT"); // if it doesn't exist then use SRT_PUBLIC_CONTEXT
      }
    else base=getenv("SRT_PUBLIC_CONTEXT");
    
    if(base=="") {
      MSG("MCReweight",Msg::kFatal)<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
      assert(false);
    }
    topDir = base+ "/" + topDir;
  }

  MSG("Zfluk",Msg::kDebug) <<"Zfluk reading data from: "<<topDir<<std::endl;
  std::string fileName =topDir+"/fluka05ptxf.root";

  TFile* hFile=new TFile(fileName.c_str());
  
  fPlist.push_back(kPiPlus);
  fPlist.push_back(kPiMinus);
  fPlist.push_back(kKPlus);
  fPlist.push_back(kKMinus);
  fPlist.push_back(kK0L);
  
  for (std::vector<ParticleType_t>::iterator itPlist=fPlist.begin();itPlist!=fPlist.end(); itPlist++)
    {
      TH2F* hist=dynamic_cast <TH2F*> (hFile->Get(Form("hF05ptxf%s",GetParticleName(*itPlist).c_str()))->Clone());
      hist->SetDirectory(0);
      TH2F* hist2=dynamic_cast <TH2F*> (hist->Clone(Form("hWeightedPTXF%s",GetParticleName(*itPlist).c_str())));
      hist2->SetDirectory(0);
      hist2->SetTitle(Form("%s weighted pt-pz",GetParticleName(*itPlist).c_str()));
      fWeightHist[*itPlist]=dynamic_cast<TH2F*> (hist->Clone(Form("hWeight%s",GetParticleName(*itPlist).c_str())));

      fWeightHist[*itPlist]->Divide(hist2);
      fWeightHist[*itPlist]->SetDirectory(0);

      std::pair<Zfluk::ParticleType_t, TH2F* > p(*itPlist, hist);
      std::pair<Zfluk::ParticleType_t, TH2F* > p2(*itPlist, hist2);
      fPTXF.insert(p);
      fWeightedPTXF.insert(p2);

      fPTXF[*itPlist]->SetDirectory(0);
      fWeightedPTXF[*itPlist]->SetDirectory(0);

      fMeanPT[*itPlist]=fPTXF[*itPlist]->ProjectionY()->GetMean()*1000.;
      fWeightedMeanPT[*itPlist]=fWeightedPTXF[*itPlist]->ProjectionY()->GetMean()*1000.;
            
      double N=fPTXF[*itPlist]->ProjectionY()->GetSumOfWeights();
      double wN=fWeightedPTXF[*itPlist]->ProjectionY()->GetSumOfWeights();

      fN[*itPlist]=N;
      fNWeighted[*itPlist]=wN;
    }
  hFile->Close();
  delete hFile;
  hFile=0;
}

Zfluk::~Zfluk() {
}

double Zfluk::GetWeight(int ptype, double pT, double pz, double* par)
{
  //this was used for PRL (par was an array of 7 parameters
  std::vector<double> parVec;
  int npar=0;

  if (whichParameterization==0) npar=7;
  else npar=6;
  
  for (int i=0;i<npar;i++) parVec.push_back(par[i]);
  //need 16 parameters to call GetWeight(ptype,pT,pz,parVec)
  // fill rest of the parameters so that weights for pi- and K- are 1
  for (int i=7;i<16;i++) parVec.push_back(double((i-1)%2));
  return GetWeight(ptype,pT,pz,parVec);
}

double Zfluk::GetWeight(int tpType, double pT, double tpz, std::vector<double> par)
{
  double xF;
  double weight=1.;
  if (tpz<0.) return weight;
  double A,B,C;
  double Ap,Bp,Cp;
  xF=tpz/120.;

  // This is SKZP parameterization
  if (whichParameterization==0)
    {
      if (tpType==211||tpType==8)
        {
          //Calculate weight for pions
          if (pT>=0.05)
            {
              // calculate the A, B and C in SKZP parameterization
              // A, B and C are best fit to Fluka 05
              
              A=-0.00761*pow((1.-xF),4.045)*(1.+9620.*xF)*pow(xF,-2.975);
              B=0.05465*pow((1.-xF),2.675)*(1.+69590.*xF)*pow(xF,-3.144);
              
              if (xF<0.22){
                C=-7.058/pow(xF,-0.1419)+9.188;
              }
              else{
                C = (3.008/exp((xF-0.1984)*3.577)) + 2.616*xF+0.1225;
              }
              // scale/skew A, B and C
              Ap=(par[0]+par[1]*xF)*A;
              Bp=(par[2]+par[3]*xF)*B;
              Cp=(par[4]+par[5]*xF)*C;
              
              weight=(Ap+Bp*pT)/(A+B*pT)*exp(-(Cp-C)*pow(pT,3./2.));
            }
          else weight=GetWeight(tpType,0.051,tpz,par);  // for low pT
        }
      else if (tpType==321||tpType==11)
        {
          //for kaons
          weight=par[6]+par[7]*xF;
        }
      else if (tpType==-211||tpType==9)
        {
          // calculate the A, B and C in SKZP parameterization
          // A, B and C are best fit to Fluka 05
          if (pT>=0.05)
            {
              A=-0.006327*pow((1.-xF),5.734)*(1.+13660.*xF)*pow(xF,-2.898);
              B=0.04605*pow((1.-xF),3.291)*(1.+58610.*xF)*pow(xF,-3.209);
              
              if (xF<0.22){
                C=-16.02/pow(xF,-0.06456)+17.61;
              }
              else{
                C = (2.963/exp((xF-0.1774)*2.265)) + 1.732*xF+0.03925;
              }
              
              // scale/skew A, B and C
              Ap=(par[8]+par[9]*xF)*A;
              Bp=(par[10]+par[11]*xF)*B;
              Cp=(par[12]+par[13]*xF)*C;
              
              weight=(Ap+Bp*pT)/(A+B*pT)*exp(-(Cp-C)*pow(pT,3./2.));
            }
          else weight=GetWeight(tpType,0.051,tpz,par);  // for low pT
          
        }
      else if (tpType==-321||tpType==12)
        {
          //for kaons
          weight=par[14]+par[15]*xF;
        }
    }
  else if (whichParameterization==1)
    {
      // This is TV parameterization
      if (tpType==211||tpType==8){
 
        weight=par[0]*(1.+par[1]*pT)/(1.+par[2]*pT)*
               exp(-par[3]*tpz*pow(pT,par[4]));
 
      }
      else if (tpType==321||tpType==11) {
        weight=par[5];
      }
      
    }
  else
    {
      MSG("MCReweight",Msg::kError) << "Please select one of the two hadron production"<<std::endl;
      MSG("MCReweight",Msg::kError) << " parameterizations (SKZP (default) or TV)"<<std::endl;
      return 0.;;
    }
  return weight;
}

double Zfluk::GetWeight(int ptype, double pT, double pz)
{
  double weight=1.;

  if (fPar.size()==0) 
    {
      MAXMSG("MCReweight",Msg::kWarning,10)
        <<"You need to set the parameters before calling "
        <<"Zfluk::GetWeight (use SetParameters(vector<double>))"<<std::endl
        <<"Returning weight = "<<weight<<std::endl;
      return weight;
    }

  double xF=pz/120.;
  double A,B,C;
  double Ap,Bp,Cp;
  Zfluk::ParticleType_t eptype=GetParticleEnum(ptype);

  // This is SKZP parameterization
  if (xF>1.||xF<0.) return weight;
  if (pT>1.||pT<0.) return weight;

  if (eptype==kPiPlus)
    {
      //Calculate weight for pions
      if (pT>=0.03)
        {
          // calculate the A, B and C in SKZP parameterization
          // A, B and C are best fit to Fluka 05
          
          A=-0.00761*pow((1.-xF),4.045)*(1.+9620.*xF)*pow(xF,-2.975);
          B=0.05465*pow((1.-xF),2.675)*(1.+69590.*xF)*pow(xF,-3.144);
          
          if (xF<0.22){
            C=-7.058/pow(xF,-0.1419)+9.188;
          }
          else{
            C = (3.008/exp((xF-0.1984)*3.577)) + 2.616*xF+0.1225;
          }
          
          // scale/skew A, B and C
          Ap=(fPar[0]+fPar[1]*xF)*A;
          Bp=(fPar[2]+fPar[3]*xF)*B;
          Cp=(fPar[4]+fPar[5]*xF)*C;
            
          weight=(Ap+Bp*pT)/(A+B*pT)*exp(-(Cp-C)*pow(pT,3./2.));
        }
      else weight=GetWeight(ptype,0.031,pz);  // for low pT
    }
  else if (eptype==kKPlus)
    {
      if (pT>=0.05)
        {
          // calculate the A, B and C in SKZP parameterization
          // A, B and C are best fit to Fluka 05
          
          A=-0.005187*pow((1.-xF),4.119)*(1.+2170.*xF)*pow(xF,-2.767);
          B=0.4918*pow((1.-xF),2.672)*(1.+1373.*xF)*pow(xF,-2.927);
          
          if (xF<0.22){
            C=-16.10/pow(xF,-0.04582)+17.92;
          }
          else{
            C = (6.905/exp((xF+0.163)*6.718)) - 0.4257*xF + 2.486;
          }
          
          // scale/skew A, B and C
          Ap=(fPar[6]+fPar[7]*xF)*A;
          Bp=(fPar[8]+fPar[9]*xF)*B;
          Cp=(fPar[10]+fPar[11]*xF)*C;

          weight=(Ap+Bp*pT)/(A+B*pT)*exp(-(Cp-C)*pow(pT,3./2.));
        }
      else 
        weight=GetWeight(ptype,0.051,pz);  // for low pT
    }
  else if (eptype==kPiMinus)
    {
      weight=GetWeight(kPiPlus,pT,pz);
      if ( pz > 4. ) weight *= ( fPar[12] + fPar[13] * xF );
    }
  else if (eptype==kKMinus)
    {
      weight=GetWeight(kKPlus,pT,pz)*(fPar[14]+fPar[15]*xF);
    }
  else if (eptype==kK0L)
    {
      //N(K0L) approximately given with (N(K+)+3*N(K-))/4
      weight= ((fNWeighted[kKPlus]+3.*fNWeighted[kKMinus])
               /(fN[kKPlus]+3.*fN[kKMinus]));
    }
  if (weight>10.) weight=10.;
  return weight;
}

void Zfluk::RecalculateWeights()
{
  
  for (std::vector<ParticleType_t>::iterator itPlist=fPlist.begin();
       itPlist!=fPlist.end(); itPlist++)
    {
      double xf,pt;
      double meanpt(0.);
      double sum(0.);
      for (int i=0;i<fPTXF[*itPlist]->GetNbinsX()+1;i++)
        {
          for (int j=0;j<fPTXF[*itPlist]->GetNbinsY()+1;j++)
            {
              xf=fPTXF[*itPlist]->GetXaxis()->GetBinCenter(i);
              pt=fPTXF[*itPlist]->GetYaxis()->GetBinCenter(j);
              fWeightedPTXF[*itPlist]
                ->SetBinContent(i,j,(fPTXF[*itPlist]->GetBinContent(i,j)
                                     *GetWeight(*itPlist,pt,xf)));
              fWeightHist[*itPlist]->SetBinContent(i,j,GetWeight(*itPlist,pt,xf));
              meanpt+=fWeightedPTXF[*itPlist]->GetBinContent(i,j)*pt;
              sum+=fWeightedPTXF[*itPlist]->GetBinContent(i,j);
            }
        }
      meanpt/=sum;
      fWeightedMeanPT[*itPlist]=meanpt*1000.; //GeV to MeV
      fNWeighted[*itPlist]=sum;
      meanpt=0.;
      sum=0.;
    }
}

double Zfluk::GetPTshift(int ptype, double pz_low, double pz_up)
{
  //returns the pT shift for ptype particle due to reweighting in xf slice
  //given by pz_low and pz_up
  //this function returns correct ptshift only if the parameters were set using SetParameters

  Zfluk::ParticleType_t eptype=GetParticleEnum(ptype);

  if (fPTXF.find(eptype)==fPTXF.end())
    {
      return 0.;
    }
  if (pz_low==0.&&pz_up==120.) return fWeightedMeanPT[eptype]-fMeanPT[eptype];
 
  int xlow=int(pz_low+1.);
  int xup=int(pz_up+1.);
 
  double pt;
  double meanpt(0.);
  double sum(0.);
  for (int i=xlow;i<xup;i++)
    {
      for (int j=1;j<51;j++)
        {
          pt=0.02*j-0.01;
          meanpt+=fWeightedPTXF[eptype]->GetBinContent(i,j)*pt;
          sum+=fWeightedPTXF[eptype]->GetBinContent(i,j);
        }
    }
  meanpt/=sum;
  meanpt*=1000.;  //GeV to MeV

  return meanpt-fMeanPT[eptype];

}
double Zfluk::GetNFrac(int ptype, double pz_low, double pz_up, double pt_low, double pt_up)
{
  //returns the fractional change in number of particles due to reweighting in xf-pt region 
  //given by pz_low, pz_up and pt_low, pt_up
  
  Zfluk::ParticleType_t eptype=GetParticleEnum(ptype);
  if (fPTXF.find(eptype)==fPTXF.end())
    {
      return 1.;
    }
  int xlow=int(pz_low+1.);
  int xup=int(pz_up+1.);
  int ylow=int(50.*pt_low+1.);
  int yup=int(50.*pt_up+1.);
  
  double wsum(0.);
  double sum(0.);
  for (int i=xlow;i<xup;i++)
    {
      for (int j=ylow;j<yup;j++)
        {
          sum+=fPTXF[eptype]->GetBinContent(i,j);
          wsum+=fWeightedPTXF[eptype]->GetBinContent(i,j);
        }
    }

  return wsum/sum;
}

std::string Zfluk::GetParticleName(Zfluk::ParticleType_t ptype)
{
switch (ptype)
   {
   case kPiPlus:          return "PiPlus";    break;
   case kPiMinus:         return "PiMinus";   break;
   case kKPlus:           return "KPlus";     break;
   case kKMinus:          return "KMinus";    break;
   case kK0L:             return "K0L";       break;
   case kUnknown:         return "Unknown";   break;
   default:               return "Unknown";   break;
   }
 
   return "Unknown";
}

Zfluk::ParticleType_t Zfluk::GetParticleEnum(int ptype)
{
  switch (ptype)
    {
    case 8   :       return kPiPlus;      break;
    case 211 :       return kPiPlus;      break;    
    case 9   :       return kPiMinus;     break;
    case -211:       return kPiMinus;     break;
    case 11  :       return kKPlus;       break;
    case 321 :       return kKPlus;       break;
    case 12  :       return kKMinus;      break;
    case -321:       return kKMinus;      break;
    case 10  :       return kK0L;         break;
    case 130 :       return kK0L;         break;
    default:         return kUnknown;     break;
    }
  return kUnknown;

}

---