Lfluk.cxx

Go to the documentation of this file.

#include "MCReweight/Lfluk.h"
#include "MessageService/MsgService.h"
#include "Registry/Registry.h"

#include <string>
#include <iostream>
#include <iomanip>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <sstream>
#include "TFile.h"
#include "TH2F.h"
#include "TF1.h"
#include "TSystem.h"
ClassImp(Lfluk)

CVSID("$Id");

Lfluk::Lfluk(/*const char* dir*/) 
   :fParam("skzp"),
    fSkew("lin"),
    fNA49Fit("Poly3"),
    fPath(""),
    fMinW(-999.0),
    fMaxW(999.0),
    fKPiRatScale(1.0),
    fWghtFromHist(true),
    fUseFlukPi(false),
    fPar(0),
    fFlukKPiRat(0)
{ 
  
}

//-------------------------------------------------------------------------------------
Lfluk::~Lfluk()
{
}

//------------------------------------------------------------------------------------------
void Lfluk::Init()
{
   cout << "Lfluk::Init" << endl;
   
   std::string dir = "";
    TH1::AddDirectory(kFALSE);

   // 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("Lfluk",Msg::kDebug) <<"Lfluk reading data from: "<<topDir<<std::endl;
   std::string fileName =topDir+"/fluka05ptxf.root";
   
   TFile* hFile=new TFile(fileName.c_str(), "READ");
   if(!(hFile->IsOpen()))
   {
       hFile = new TFile("/minos/scratch/loiacono/Data/ForMonFit/fluka05ptxf.root", "READ");
       if(!hFile || !(hFile->IsOpen()))
       {
          cout << "Can't open /minos/scratch/loiacono/Data/ForMonFit/fluka05ptxf.root" << endl;
          cout << "   Failed to open file: " << fileName << endl;
          cout << "   **************I am going to crash***************"<< endl;
       }
    
   }
   else
   {
      cout << "   Reading Fluka PTXF from: " << fileName << endl;
   }

   
   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()));

      const int nbinsx = 120;
      const int nbinsy = hist -> GetNbinsY();
      //cout << "nbinsx = " << nbinsx << "nbinsy = " << nbinsy << endl;
      //fWeightHist[*itPlist]=dynamic_cast<TH2D*> (hist->Clone(Form("hWeight%s",GetParticleName(*itPlist).c_str())));
      //
      //The bins need to be the same as the bins that the pt-pz plots for the horn current scans!!!!!!!!!!!!!!!!!!!!!!!
      //
      fWeightHist[*itPlist] = new TH2D(Form("hWeightedPTXF%s",GetParticleName(*itPlist).c_str()),GetParticleName(*itPlist).c_str(), nbinsx, 0, 120, nbinsy, 0, 1);
      for(int ibinx = 0; ibinx <= nbinsx+1; ++ibinx)
      {
         for(int ibiny = 0; ibiny <= nbinsy+1; ++ibiny)
         {
            //cout << "ibinx = " << ibinx << "ibiny = " << ibiny << endl;
            fWeightHist[*itPlist] -> SetBinContent(ibinx, ibiny, 1.0);
         }
      }
      //TH2D* temp = dynamic_cast<TH2D*>(hist->Clone());
      //fWeightHist[*itPlist]->Divide(temp);
      fWeightHist[*itPlist]->SetDirectory(0);

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

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

      //
      // Projection DOES include overflow bins
      // Get Mean DOES NOT
      //
      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;

   //
   //Get NA49 Ratio from Fit.
   //
   std::string path = string(fPath + "na49Ratio.root");
   TFile *fitfile = new TFile(path.c_str(), "READ");
   if(!(fitfile->IsOpen()))
   {
      cout << "   Failed to open file: " << path << endl;
      cout << "   **************I am going to crash***************"<< endl;
   }
   else
   {
      cout << "   Reading NA49 Pi Ratio from: " << path << endl;
   }

   stringstream name;
   name << "NA49PipPimRatio_" << fNA49Fit << "_fit";
   if(!(fitfile -> Get(name.str().c_str())))
   {
      cout << "   Failed to get " << name.str() << " from " << path << endl;
   }
   fNA49pirat = dynamic_cast<TF1 *>( (fitfile -> Get(name.str().c_str())) -> Clone());
   if(!fNA49pirat)
   {
      cout << "   Failed to get " << name.str() << " from " << path << endl;
      fNA49pirat = 0;
   }
   fitfile -> Close();
   delete fitfile;


   /*
   //
   //added stuff above xf=0.5.
   //see MuonPhysics/Fit/Misc/NA49PiRatio.cxx
   //

   //na49 data thru xf = 0.5
   double xfna49[] ={0.01,   0.02,   0.03,   0.04,   0.05,
                     0.075,  0.1,    0.125,  0.15,   0.2,
                     0.25,   0.3,    0.4,    0.5,    70.5,
                     80.5,   90.5,  100.5,  119.5};
   double pipna49[]={25.756, 20.903, 16.683, 13.389, 10.981,
                     7.067,  4.904,  3.531,  2.648,  1.5499,
                     0.9734, 0.5996, 0.2455, 0.1081, 3.69741,
                     4.36113,5.3943, 5.89819,7.60617};
   double pimna49[]={24.424, 19.247, 14.938, 11.757, 9.516,
                     5.790,  3.803,  2.561,  1.826,  0.9555,
                     0.5475, 0.3190, 0.1041, 0.0336, 1.0,
                     1.0,    1.0,    1.0,    1.0};
   //double piratna49[14];
   double piratna49[19];
   
   for (int i=0;i<19;i++) piratna49[i]=pipna49[i]/pimna49[i];
   for (int i=0;i<14;i++) xfna49[i]*=120.;
   TGraph* pina49=new TGraph(19,xfna49,piratna49);                                                                                    
   fNA49pirat = new TMVA::TSpline1("sp0", pina49);
*/
   
   fPiRatio.resize(122);
   fFluka05KPiRatio.resize(122);
   fFluka05PiRatio.resize(122);
   fKPiRatio.resize(122);
   TH1D* piplus=dynamic_cast<TH1D*> (GetPTXF(8)
                                     ->ProjectionX()->Clone("piplus"));
   TH1D* piminus=dynamic_cast<TH1D*> (GetPTXF(9)
                                      ->ProjectionX()->Clone("piplus"));
   TH1D* kplus=dynamic_cast<TH1D*> (GetPTXF(11)
                                    ->ProjectionX()->Clone("kplus"));
   
   kplus->Divide(piplus);
   piplus->Divide(piminus);
   
   fFlukKPiRat = dynamic_cast<TH1D*>( kplus->Clone("FlukKPiRatio") );
   fFlukPiRat = dynamic_cast<TH1D*>( piplus->Clone("FlukPiRatio") );
  

   
   for (int i=1;i<121;i++)
   {
      if (kplus->GetBinContent(i)!=0.) 
         fFluka05KPiRatio[i] = kplus->GetBinContent(i);
      else if (i<119) 
         fFluka05KPiRatio[i] = (fFluka05KPiRatio[i-2]+fFluka05KPiRatio[i-1]
                                +kplus->GetBinContent(i+1)+kplus->GetBinContent(i+2))/4.;
      else
         fFluka05KPiRatio[i] = fFluka05KPiRatio[i-1];
      
      fFlukKPiRat -> SetBinContent(i, fFluka05KPiRatio[i]);
   }

  

   for (int i=1;i<121;i++)
   {
      if (piplus->GetBinContent(i)!=0.) 
         fFluka05PiRatio[i] = piplus->GetBinContent(i);
      else if (i<119) 
         fFluka05PiRatio[i] = (fFluka05PiRatio[i-2]+fFluka05PiRatio[i-1]
                                +piplus->GetBinContent(i+1)+piplus->GetBinContent(i+2))/4.;
      else
         fFluka05PiRatio[i] = fFluka05PiRatio[i-1];
      
      fFlukPiRat -> SetBinContent(i, fFluka05PiRatio[i]);
   }


   isFirst=true;


}
//------------------------------------------------------------------------------------------
void Lfluk::SetConfig(const std::string param, const std::string skew, const std::string na49fit, const double wmin, const double wmax)
{
   fParam = param;
   fSkew = skew;
   fNA49Fit = na49fit;
   fMinW = wmin;
   fMaxW = wmax;

}

//------------------------------------------------------------------------------------------
void Lfluk::Config(const Registry &reg)
{
   reg.Get("WEIGHT_MIN", fMinW);
   reg.Get("WEIGHT_MAX", fMaxW);
   reg.Get("KPiRatioScale", fKPiRatScale);

   if(fMaxW < 990.0 && fMinW > -990.0  && fMaxW <= fMinW)
   {
      fMaxW = fMinW + 10.0;
      cout << "Lfluk::Config - Weight Max <= Weight Min, setting fMaxW = " << fMaxW << " fMinW = " << fMinW << endl; 
   }

   const char *value_char = 0;

   if(reg.Get("ZRatTuneParam", value_char) && value_char)
   {
      stringstream strstrm;
      strstrm << value_char;
      fParam = strstrm.str();
   }

   value_char = 0;
   if(reg.Get("ZRatTuneSkew", value_char) && value_char)
   {
      stringstream strstrm;
      strstrm << value_char;
      fSkew = strstrm.str();
   }

   value_char = 0;
   if(reg.Get("NA49FIT", value_char) && value_char)
   {
      stringstream strstrm;
      strstrm << value_char;
      fNA49Fit = strstrm.str();
   }

   value_char = 0;
   if(reg.Get("ROOTFILES", value_char) && value_char)
   {
      stringstream strstrm;
      strstrm << value_char;
      fPath = strstrm.str();
   }

   value_char = 0;
   if(reg.Get("WGHTFROMHIST", value_char) && value_char)
   {
      if(std::strcmp(value_char, "no") == 0) 
      {
         fWghtFromHist = false;
      }
   }

   value_char = 0;
   if(reg.Get("UseFlukaPiRatio", value_char) && value_char)
   {
      if(std::strcmp(value_char, "yes") == 0) 
      {
         fUseFlukPi = true;
      }
   }

     
   cout << "Lfluk::Config" << endl;
   cout << "   Root File Path = " << fPath << endl
        << "   Get Weight From Hist = " <<  fWghtFromHist << endl
        << "   Parameterization = " << fParam << endl
        << "   Skew Config = " << fSkew << endl
        << "   Use Fluka Pi Ratio = " << fUseFlukPi << endl
        << "   NA49 Ratio Fit = " << fNA49Fit << endl
        << "   KPi Ratio Scale = " << fKPiRatScale << endl
        << "   Weight Min = " << fMinW << endl
        << "   Weight Max = " << fMaxW << endl;
   
   
}

//------------------------------------------------------------------------------------------
void Lfluk::SetParameters(std::vector<double> par) 
{
   fPar=par;
   isFirst=true;

   if(fParam == "skzp" && fSkew == "fluka")
   {
      cout << " Pars = ";
      for(unsigned int ip = 0; ip < fPar.size(); ++ip)
      {
         cout << fPar[ip] << " ";
      }
      cout << "end pars " << endl;

      /*
      Lfluk::ParticleType_t eptype=GetParticleEnum(8);
      
      double xF; bool print = false;
      for (int i=0;i<fPTXF[eptype]->GetNbinsX()+1;i++)
      {
         xF=fPTXF[eptype]->GetXaxis()->GetBinCenter(i);
         xF = xF/120.0;
         if (xF>1.||xF<0.) continue;
         
         if( fabs((int)(xF*10) - xF*10) < 0.1) print = true;
         
         if(print) cout << "xF = " << xF;
         
         double A = Lfluk::GetSKZPA(eptype, xF);
         if(print) cout << " A = " << A;
         
         double B = Lfluk::GetSKZPB(eptype, xF);
         if(print) cout << " B = " << B;
         
         double C = Lfluk::GetSKZPC(eptype, xF);
         if(print) cout << " C = " << C;
         
         long double Ap= fPar[0]*pow((1.-xF),fPar[1])*(1.+fPar[2]*xF)*pow(xF,fPar[3]);
         if(print) cout << " Ap = " << Ap;
         
         long double Bp= fPar[4]*pow((1.-xF),fPar[5])*(1.+fPar[6]*xF)*pow(xF,fPar[7]);
         if(print) cout << " Bp = " << Bp;
         
         long double Cp;   
         if (xF<0.22)
         {
            Cp= fPar[8]/pow(xF,fPar[9])+fPar[10];
            if(Cp < -100.0) Cp = -100.0;
            if(Cp > 100.0) Cp = 100.0;
            if(print) cout << " Cp = " << Cp;
         }
         else
         {
             double arg22 = (0.22-fPar[11])*fPar[12];
            if(arg22 < -100.0) arg22 = -100.0;
            if(arg22 > 100.0) arg22 = 100.0;
            double match = (fPar[8]/pow(0.22,fPar[9])+fPar[10]-
                            fPar[13]*0.22-fPar[14])*exp(arg22);
            double exparg = (xF-fPar[11])*fPar[12];
            if(exparg < -100.0) exparg = -100.0;
            if(exparg > 100.0) exparg = 100.0;
            Cp = ( match / exp(exparg)) + fPar[13]*xF+fPar[14];
            
            if(Cp < -100.0) Cp = -100.0;
            if(Cp > 100.0) Cp = 100.0;
            if(print) cout << " Cp = " << Cp;
         }
         
         
         if(print)
         {
            cout << " exp(-(Cp-C)) = " << exp(-(Cp-C)) << endl;
            
            cout << "   pT = 0.01, (Ap+Bp*pT)/(A+B*pT) = " << (Ap+Bp*0.01)/(A+B*0.01);
            double weight=(Ap+Bp*0.01)/(A+B*0.01)*exp(-(Cp-C)*pow(0.01,3./2.));
            cout << " weight = " << weight << endl;

            cout << "   pT = 0.95, (Ap+Bp*pT)/(A+B*pT) = " << (Ap+Bp*0.95)/(A+B*0.95);
            weight=(Ap+Bp*0.95)/(A+B*0.95)*exp(-(Cp-C)*pow(0.95,3./2.));
            cout << " weight = " << weight << endl;

         }
         print = false;
      
      }
      */

   cout << "Done Setting pars...recalc weights" << endl;
   cout << endl;

   }




   
   RecalculateWeights();
}

//-------------------------------------------------------------------------------------
TH2* Lfluk::GetWeightHist(const std::string &ptype)
{
   TH2* hist = 0;
   hist = dynamic_cast <TH2*> (fWeightHist[GetParticleEnum(ptype)]->Clone());
   hist->SetDirectory(0);
   
   return hist;

}

//-------------------------------------------------------------------------------------
double Lfluk::GetWeight(int ptype, double pT, double pz)
{
   
   
   if(fWghtFromHist)
   {
            
      double weight = 1.0;

      Lfluk::ParticleType_t eptype=GetParticleEnum(ptype);
      string name = GetParticleName(eptype);  
      
      if(eptype != kPiPlus && eptype != kKPlus && eptype != kPiMinus && eptype != kPiMinus)
      {
         return 1.0;
      }
      
      if (pz > 120.0 || pz < 0.0 || pT < 0.0)
      {
         return weight;
      }
      
      const int nbinsx = fWeightHist[eptype] -> GetNbinsX();
      const int nbinsy = fWeightHist[eptype] -> GetNbinsY();
      
      int binx = (fWeightHist[eptype] -> GetXaxis()) -> FindFixBin(pz);
      int biny = (fWeightHist[eptype] -> GetYaxis()) -> FindFixBin(pT);
      
      if (pT > 1.0) biny = nbinsy+1;
      
      if(!(binx < 0) && !(biny < 0) && !(binx > nbinsx+1) && !(biny > nbinsy+1))
      {
         weight = fWeightHist[eptype] -> GetBinContent(binx, biny);
      }

      /*
      double wt = CalcWeight(eptype,pT,pz);
      string name = GetParticleName(eptype);            
      //if(name == "PiMinus" && pz > 117 && pz < 119 && weight-wt != 0)
      if(weight - wt != 0)
      {
         cout << " Lfluk::GetWeight(): ptype = " << name << " pz = " << pz << " pt = " << pT << " weight from calc = "
              << setprecision(10) << wt << " weight from hist = " << weight << " diff = " << weight - wt << endl;
      }
      */

      /*
      if((name == "PiPlus" || name == "KPlus") && (binx == 2) && biny == nbinsy+1)
      {
         cout << "get weight called : name = " << name << " pz = " << pz << " pt = " << pT << " weight = " << weight << endl;
      }
      
      */

      
      return weight;
   }

   else
   {
      return CalcWeight(ptype,pT,pz);
   }
}

//-------------------------------------------------------------------------------------
double Lfluk::CalcWeight(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::CalcWeight (use SetParameters(vector<double>))"<<std::endl
        <<"Returning weight = "<<weight<<std::endl;
      return weight;
    }

  double xF=pz/120.;
  Lfluk::ParticleType_t eptype=GetParticleEnum(ptype);

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

     if(fParam.find("skzp") != string::npos) //default
     {
        weight = Lfluk::GetSKZPWeight(eptype, pT, xF);
     }
     if(fParam.find("skpc") != string::npos)
     {
        weight = Lfluk::GetSKPCWeight(eptype, pT, xF);
     }
     if(fParam.find("quad") != string::npos)
     {
        weight = Lfluk::GetQuadWeight(eptype, pT, xF);
     }
     if(fParam.find("meanpt") != string::npos)
     {
        weight = Lfluk::GetMeanPTWeight(eptype, pT, xF);
     }

     if(weight < fMinW) weight = fMinW;
     if(weight > fMaxW) weight = fMaxW;

  }
  else if (eptype==kKPlus)
  {
     weight=fKPiRatio[int(pz)+1];
  }
  else if (eptype==kPiMinus)
  {
     weight=CalcWeight(kPiPlus,pT,pz)*fPiRatio[int(pz)+1];  
  }
  else if (eptype==kKMinus)
  {
     weight=CalcWeight(kKPlus,pT,pz);
  }
  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 < -999.0) weight = fMinW;
  if(weight > 999.0) weight = fMaxW;


  return weight;
}


//------------------------------------------------------------------------------------------
double Lfluk::GetSKZPWeight(const Lfluk::ParticleType_t particle, const double pT, const double xF)
{
   double weight = 1.0;
   double A,B,C;
   long double Ap,Bp,Cp;
   
   if (particle==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 = Lfluk::GetSKZPA(particle, xF);
         B = Lfluk::GetSKZPB(particle, xF);
         C = Lfluk::GetSKZPC(particle, xF);
         
         Lfluk::SkewSKZP(xF, A, B, C, Ap, Bp, Cp);

         if(Cp < -100.0) Cp = -100.0;
         if(Cp > 100.0) Cp = 100.0;
         weight=(Ap+Bp*pT)/(A+B*pT)*exp(-(Cp-C)*pow(pT,3./2.));
      }
      else
      {
         weight = Lfluk::GetSKZPWeight(particle, 0.031, xF);  // for low pT
      }
   }
   else if (particle==kKPlus)
   {
      weight = 1.0;
   }
   else if (particle==kPiMinus)
   {
      weight = 1.0;
   }
   else if (particle==kKMinus)
   {
      weight = 1.0;
   }
   
   return weight;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKPCWeight(const Lfluk::ParticleType_t particle, const double pT, const double xF)
{
   double weight = 1.0;
   double A, B, C, D, E;
   double Ap, Bp, Cp, Dp, Ep;
   
   if (particle==kPiPlus)
   {
      if (pT>=0.03)
      {
         A = Lfluk::GetSKPCA(particle, xF);
         B = Lfluk::GetSKPCB(particle, xF);
         C = Lfluk::GetSKPCC(particle, xF);
         D = Lfluk::GetSKPCD(particle, xF);
         E = Lfluk::GetSKPCE(particle, xF);
         
         Lfluk::SkewSKPC(xF, A, B, C, D, E, Ap, Bp, Cp, Dp, Ep);

         if(Ep < -20.0) Ep = -20.0;
         double arg = -Dp*std::pow(pT,Ep)+D*std::pow(pT,E);
         if(arg > 100.0) arg = 100.0;
         if(arg < -100.0) arg = -100.0;
         weight=(Ap+Bp*pT+Cp*pT*pT)/(A+B*pT+C*pT*pT)*std::exp(arg);
      }
      else
      {
         weight = Lfluk::GetSKPCWeight(particle, 0.031, xF);  // for low pT
      }
   }
   else if (particle==kKPlus)
   {
      weight = 1.0;
   }
   else if (particle==kPiMinus)
   {
      weight = 1.0;
   }
   else if (particle==kKMinus)
   {
      weight = 1.0;
   }
   
   return weight;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetQuadWeight(const Lfluk::ParticleType_t particle, const double pT, const double xF)
{
   double weight = 1.0;
      
   if(particle == kPiPlus)
   {
      if(fParam == "quad0")
      {
         if(fPar.size() != 9) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 9: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    + fPar[5]*pT*xF*xF
                   + fPar[6]*pT*pT + fPar[7]*pT*pT*xF + fPar[8]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad1")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    + fPar[5]*pT*xF*xF
                   + fPar[6]*pT*pT + fPar[7]*pT*pT*xF ;
         }  
      }
      else if(fParam == "quad2")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    + fPar[5]*pT*xF*xF
                   + fPar[6]*pT*pT                    + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad3")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    + fPar[5]*pT*xF*xF
                                   + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad4")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad5")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT                       + fPar[4]*pT*xF*xF
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad6")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                                   + fPar[3]*pT*xF    + fPar[4]*pT*xF*xF
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad7")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       
                   + fPar[2]*pT    + fPar[3]*pT*xF    + fPar[4]*pT*xF*xF
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad8")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]                          + fPar[1]*xF*xF
                   + fPar[2]*pT    + fPar[3]*pT*xF    + fPar[4]*pT*xF*xF
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad9")
      {
         if(fPar.size() != 8) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 8: npars = " << fPar.size() << endl;
         }
         else
         {
            weight =                 fPar[0]*xF       + fPar[1]*xF*xF
                   + fPar[2]*pT    + fPar[3]*pT*xF    + fPar[4]*pT*xF*xF
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF + fPar[7]*pT*pT*xF*xF;
         }  
      }
      else if(fParam == "quad124")
      {
         if(fPar.size() != 6) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 7: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    
                   + fPar[5]*pT*pT  ;
         }  
      }
      else if(fParam == "quad134")
      {
         if(fPar.size() != 6) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 7: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    
                                   + fPar[5]*pT*pT*xF ;
         }  
      }
      else if(fParam == "quad1234")
      {
         if(fPar.size() != 5) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 7: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF ;
         }  
      }
      else if(fParam == "quad14")
      {
         if(fPar.size() != 7) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 7: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF       + fPar[2]*xF*xF
                   + fPar[3]*pT    + fPar[4]*pT*xF    
                   + fPar[5]*pT*pT + fPar[6]*pT*pT*xF ;
         }  
      }
      else if(fParam == "quad147")
      {
         if(fPar.size() != 6) 
         {
            cerr << "Problem in  Lfluk::GetQuadWeight - Number of pars is not 6: npars = " << fPar.size() << endl;
         }
         else
         {
            weight = fPar[0]       + fPar[1]*xF
                   + fPar[2]*pT    + fPar[3]*pT*xF    
                   + fPar[4]*pT*pT + fPar[5]*pT*pT*xF ;
         }  
      }
      else
      {
         cerr << "Problem in  Lfluk::GetQuadWeight - Not a valid Parameterization " << fParam << endl;
      }


   }

   //cout << "xF = " << xF << " pT = " << pT << " quad weight = " << weight << endl;
   return weight;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetMeanPTWeight(const Lfluk::ParticleType_t particle, const double pT, const double xF)
{
   double weight = 1.0;
      
   if(particle == kPiPlus)
   {
      if(fPar.size() != 5) 
      {
         cout << "Problem in  Lfluk::GetMeanPTWeight - Number of pars is not 4: npars = " << fPar.size() << endl;
      }
      else
      {
         double linxF = fPar[3]+fPar[4]*xF;
         if(linxF == 0.0)
         {
            linxF = 1e-20;
         }
         double exparg = fPar[2]*pT / linxF;
         if(exparg > 100.0) exparg = 100.0;
         if(exparg < -100.0) exparg = -100.0; 
         weight = (fPar[0] + fPar[1]*xF) * exp(exparg);
      }
   }

   
   //cout << "xF = " << xF << " pT = " << pT << " meanPT weight = " << weight << endl;
   return weight;
}

//-------------------------------------------------------------------------------------
void Lfluk::SkewSKPC(const double xF, const double A, const double B, const double C,
                 const double D, const double E, double &Ap, double &Bp, double &Cp,
                 double &Dp, double &Ep)
{
   
   if(fSkew.find("lin") != string::npos)
   {
      if(fPar.size() == 10) 
      {
         Ap=(fPar[0]+fPar[1]*xF)*A;
         Bp=(fPar[2]+fPar[3]*xF)*B;
         Cp=(fPar[4]+fPar[5]*xF)*C;
         Dp=(fPar[6]+fPar[7]*xF)*D;
         Ep=(fPar[8]+fPar[9]*xF)*E;

         return;
      }
   }

   if(fSkew.find("quad") != string::npos)
   {
      if(fPar.size() == 15) 
      {
         Ap=(fPar[0]+fPar[1]*xF+fPar[2]*xF*xF)*A;
         Bp=(fPar[3]+fPar[4]*xF+fPar[5]*xF*xF)*B;
         Cp=(fPar[6]+fPar[7]*xF+fPar[8]*xF*xF)*C;
         Dp=(fPar[9]+fPar[10]*xF+fPar[11]*xF*xF)*D;
         Ep=(fPar[12]+fPar[13]*xF+fPar[14]*xF*xF)*E;

         return;
      }
   }
   
   cout << "Lfluk::SkewSKPC() - No valid skew option or wrong number of parameters" << endl
        << "       fSkew  = " << fSkew << endl
        << "       NPars  = " << fPar.size() << endl;
   
   Ap = A;
   Bp = B;
   Cp = C;
   Dp = D;
   Ep = E;
   
}

//-------------------------------------------------------------------------------------
void Lfluk::SkewSKZP(const double xF, const double A, const double B, const double C, long double &Ap, long double &Bp, long double &Cp)
{
         // scale/skew A, B and C
   
   if(fSkew.find("fluka") != string::npos)
   {
      if(fPar.size() == 15) 
      {
         Ap= fPar[0]*pow((1.-xF),fPar[1])*(1.+fPar[2]*xF)*pow(xF,fPar[3]);
         Bp= fPar[4]*pow((1.-xF),fPar[5])*(1.+fPar[6]*xF)*pow(xF,fPar[7]);
         
         if (xF<0.22)
         {
            Cp= fPar[8]/pow(xF,fPar[9])+fPar[10];
         }
         else
         {
            double arg22 = (0.22-fPar[11])*fPar[12];
            if(arg22 < -100.0) arg22 = -100.0;
            if(arg22 > 100.0) arg22 = 100.0;
            double match = (fPar[8]/pow(0.22,fPar[9])+fPar[10]-
                            fPar[13]*0.22-fPar[14])*exp(arg22);
            double exparg = (xF-fPar[11])*fPar[12];
            if(exparg < -100.0) exparg = -100.0;
            if(exparg > 100.0) exparg = 100.0;
            Cp = ( match / exp(exparg)) + fPar[13]*xF+fPar[14];
            
         }

         return;
      }
   }

   if(fSkew.find("lin") != string::npos)
   {
      if(fPar.size() == 6) 
      {
         Ap=(fPar[0]+fPar[1]*xF)*A;
         Bp=(fPar[2]+fPar[3]*xF)*B;
         Cp=(fPar[4]+fPar[5]*xF)*C;

         return;
      }
   }

   if(fSkew.find("quad") != string::npos)
   {
      if(fPar.size() == 9) 
      {
         Ap=(fPar[0]+fPar[1]*xF+fPar[2]*xF*xF)*A;
         Bp=(fPar[3]+fPar[4]*xF+fPar[5]*xF*xF)*B;
         Cp=(fPar[6]+fPar[7]*xF+fPar[8]*xF*xF)*C;

         return;
      }
   }
   
   cout << "Lfluk::SkewSKZP() - No valid skew option or wrong number of parameters" << endl
        << "       fSkew  = " << fSkew << endl
        << "       NPars  = " << fPar.size() << endl;
   
   Ap = A;
   Bp = B;
   Cp = C;
   
}


//------------------------------------------------------------------------------------------
double Lfluk::GetSKZPA(const Lfluk::ParticleType_t particle, const double xf) const
{
   
   switch (particle)
   {
       case kPiPlus:  
          return -0.007607*std::pow(1.0-xf, 4.045)*(1.0+0.9620e+4*xf)*std::pow(xf, -2.975);
       case kPiMinus:
          return -0.006306*std::pow(1.0-xf, 5.730)*(1.0+0.1365e+5*xf)*std::pow(xf, -2.900);
       case kKPlus:
          return -0.005187*std::pow(1.0-xf, 4.119)*(1.0+0.2170e+4*xf)*std::pow(xf, -2.767);
       case kKMinus:
          return -8.854   *std::pow(1.0-xf, 6.778)*(1.0-0.6050   *xf)*std::pow(xf, -1.827);
       default:
          break;
   }

   return 1.0;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKZPB(const Lfluk::ParticleType_t particle, const double xf) const
{
   switch (particle)
   {
       case kPiPlus: 
          return 0.05465*std::pow(1.0-xf, 2.675)*(1.0+0.6959e+5*xf)*std::pow(xf, -3.144);
       case kPiMinus:
          return 0.04608*std::pow(1.0-xf, 3.291)*(1.0+0.5857e+5*xf)*std::pow(xf, -3.209);
       case kKPlus:
          return 0.4918 *std::pow(1.0-xf, 2.672)*(1.0+0.1373e+4*xf)*std::pow(xf, -2.927);
       case kKMinus:
          return 0.02857*std::pow(1.0-xf, 7.494)*(1.0+0.5879e+5*xf)*std::pow(xf, -2.577);
       default:
          break;
   }

   return 0.0;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKZPC(const Lfluk::ParticleType_t particle, const double xf) const
{
   if (xf<0.22)
   {
      switch (particle)
      {
          case kPiPlus: 
             return  -7.058/std::pow(xf, -0.1419 ) +  9.188;
          case kPiMinus:
             return -16.52 /std::pow(xf, -0.06204) + 18.12;
          case kKPlus:
             return -16.10 /std::pow(xf, -0.04582) + 17.92;
          case kKMinus:
             return -16.13 /std::pow(xf, -0.05678) + 17.39;
          default:
             break;
      }
   }
   else 
   {
      switch (particle)
      {
          case kPiPlus: 
             return 3.008/std::exp((xf-0.1948)*3.577)   + 2.616 *xf + 0.1225;
          case kPiMinus:
             return 2.972/std::exp((xf-0.1758)*2.266)   + 1.730 *xf + 0.04196;
          case kKPlus:
             return 6.905/std::exp((xf+0.1630)*6.718)   - 0.4257*xf + 2.486;
          case kKMinus:
             return 3.916/std::exp((xf-4.615) *0.03255) - 4.702 *xf + 0.4062;
          default:
             break;
      }
   }
   
   
   return 0.0;
}


//-------------------------------------------------------------------------------------
void Lfluk::RecalculateWeights()
{
   
   if (isFirst) {
    for (int i=0;i<121;i++) {
      fKPiRatio[i]=1.;
      fPiRatio[i]=1.;
    }
  }
  for (std::vector<ParticleType_t>::iterator itPlist=fPlist.begin();
       itPlist!=fPlist.end(); itPlist++)
  {
     double xf,pt;
     double meanpt(0.);
     double sum(0.);
     const int nbinsx = fPTXF[*itPlist]->GetNbinsX();
     const int nbinsy = fPTXF[*itPlist]->GetNbinsY();
     for (int i=0;i<=nbinsx+1;i++)
     {
        for (int j=0;j<=nbinsy+1;j++)
        {
           xf=fPTXF[*itPlist]->GetXaxis()->GetBinCenter(i);
           pt=fPTXF[*itPlist]->GetYaxis()->GetBinCenter(j);
           double weight = CalcWeight(*itPlist,pt,xf);
           
           fWeightedPTXF[*itPlist]
              ->SetBinContent(i,j,(fPTXF[*itPlist]->GetBinContent(i,j)
                                   *weight));

           Lfluk::ParticleType_t eptype= *itPlist;
           string name = GetParticleName(eptype);
           fWeightHist[*itPlist]->SetBinContent(i,j,weight);
           if(weight - fWeightHist[*itPlist]->GetBinContent(i,j) != 0)
           //if(name == "PiPlus" && weight !=1.0)
           {
              cout << "Lfluk::RecalculateWeights(): name = " << name << " pz = " << xf << " pt = " << pt
                   << " weight from calc = " << setprecision(10) << weight << " weight from hist = "
                   << fWeightHist[*itPlist]->GetBinContent(i,j)
                      << " diff = " << weight - fWeightHist[*itPlist]->GetBinContent(i,j) << endl;
           }
           
           //
           // Projection DOES include overflow bins
           // Get Mean DOES NOT
           //
           //need to include the overflow in y for sure in fWeightedPTXF
           //because it is used in the ratios below
           //but in meanpt we don't want to use this bin
           //I have checked that the fluka "signal" in i=nbinsx+1 is always 0
           //so it doesn't matter if this bin is included or not. I think the
           //proper thing is that is should be because the projection does
           //include it
           //
           if(j == nbinsy+1) continue;
           meanpt+=fWeightedPTXF[*itPlist]->GetBinContent(i,j)*pt;
           sum+=fWeightedPTXF[*itPlist]->GetBinContent(i,j);
           
           /*Lfluk::ParticleType_t eptype= *itPlist;
           string name = GetParticleName(eptype);
           if((name == "PiPlus" || name == "KPlus") && ((i == 2 && j == nbinsy+1) || (i==nbinsx+1 && j == 2)))
           {
              cout <<"name = " << name << " pz = " << xf << " pt = " << pt
                   << " weight = " << weight << " fluksig = " << fPTXF[*itPlist]->GetBinContent(i,j)
                   << " newsig = " << fWeightedPTXF[*itPlist]->GetBinContent(i,j) << endl;
                   }*/
           
                   
        }
     }
     
     meanpt/=sum;
     fWeightedMeanPT[*itPlist]=meanpt*1000.; //GeV to MeV
     fNWeighted[*itPlist]=sum;
     meanpt=0.;
     sum=0.;
     /* 
     if(!isFirst)
     {
        for (int i=0;i<=fWeightHist[*itPlist]->GetNbinsX()+1;i++)
        {
           for (int j=0;j<=fWeightHist[*itPlist]->GetNbinsY()+1;j++)
           {
              xf=fWeightHist[*itPlist]->GetXaxis()->GetBinCenter(i);
              pt=fWeightHist[*itPlist]->GetYaxis()->GetBinCenter(j);
              double weight = CalcWeight(*itPlist,pt,xf);
              
              Lfluk::ParticleType_t eptype= *itPlist;
              string name = GetParticleName(eptype);
              fWeightHist[*itPlist]->SetBinContent(i,j,weight);
             

              //if(name == "PiMinus" && weight - fWeightHist[*itPlist]->GetBinContent(i,j) != 0)
              //if(name == "PiPlus" && xf > 49 && xf < 50 && weight - fWeightHist[*itPlist]->GetBinContent(i,j) != 0)
              //if(name == "PiMinus" && xf > 117 && xf < 119 && weight - fWeightHist[*itPlist]->GetBinContent(i,j) != 0)
              if(weight - fWeightHist[*itPlist]->GetBinContent(i,j) != 0)
              {
                 cout << "Lfluk::RecalculateWeights(): name = " << name << " pz = " << xf << " pt = " << pt
                      << " weight from calc = " << setprecision(10) << weight << " weight from hist = "
                      << fWeightHist[*itPlist]->GetBinContent(i,j)
                      << " diff = " << weight - fWeightHist[*itPlist]->GetBinContent(i,j) << endl;
              }
           }
        }
     }
     */
  }
  
  
  if (isFirst)
  {
    isFirst=false;
    const TH1D* piplus=dynamic_cast<TH1D*> (GetReweightedPTXF(8)
                                      ->ProjectionX()->Clone("piplus"));
    const TH1D* piminus=dynamic_cast<TH1D*> (GetReweightedPTXF(9)
                                       ->ProjectionX()->Clone("piminus"));
    TH1D* kplus_f05=dynamic_cast<TH1D*> (GetReweightedPTXF(11)
                                         ->ProjectionX()->Clone("kplus_f05"));
    TH1D* piplus_rat=dynamic_cast<TH1D*> (GetReweightedPTXF(8)
                                      ->ProjectionX()->Clone("piplus_f05"));

    //fPiPlusReWProj = dynamic_cast<TH1D*>(piplus->Clone());
    
    piplus_rat -> Divide(piplus, piminus);
    
    //fKPlusf05Proj = dynamic_cast<TH1D*>(kplus_f05->Clone());   

    kplus_f05->Divide(piplus);

    //fKPlusf05PipReWRatio = dynamic_cast<TH1D*>(kplus_f05->Clone());
    
    for (int i=1;i<121;i++)
    {
       if (kplus_f05->GetBinContent(i)!=0)
       {
          fKPiRatio[i] = (fFluka05KPiRatio[i]*fKPiRatScale)/kplus_f05->GetBinContent(i);
       }

       if (piplus_rat->GetBinContent(i)!=0)
       {
          if(fUseFlukPi && fFluka05PiRatio[i] != 0) fPiRatio[i] = piplus_rat->GetBinContent(i)/fFluka05PiRatio[i];
          else fPiRatio[i] = piplus_rat->GetBinContent(i)/fNA49pirat->Eval(piplus_rat->GetBinCenter(i));
          
       }
       /*if (i<62)
       {
          fPiRatio[i] = piplus->GetBinContent(i)/fNA49pirat->Eval(piplus->GetBinCenter(i));
       }
       else
       {
          fPiRatio[i]=1.;
          }*/
    }
    delete piplus;
    piplus=0;
    delete piminus;
    piminus=0;
    delete kplus_f05;
    kplus_f05=0;
    delete  piplus_rat;
    piplus_rat = 0;
    
    RecalculateWeights();
  }
  
}

//-------------------------------------------------------------------------------------
double Lfluk::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

  Lfluk::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 Lfluk::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
  
  Lfluk::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 Lfluk::GetParticleName(Lfluk::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";
}

//-------------------------------------------------------------------------------------
Lfluk::ParticleType_t Lfluk::GetParticleEnum(const std::string &ptype)
{
   if     (ptype == "PiPlus")  return kPiPlus;
   else if(ptype == "PiMinus") return kPiMinus;
   else if(ptype == "KPlus")   return kKPlus;  
   else if(ptype == "KMinus")  return kKMinus; 
   else                        return kUnknown;
}

//-------------------------------------------------------------------------------------
Lfluk::ParticleType_t Lfluk::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;

}








//------------------------------------------------------------------------------------------
double Lfluk::GetSKPCA(const Lfluk::ParticleType_t particle, const double xf) const
{
   if (xf<0.25) 
   {
      //
      //form: Par[0]*std::pow(1.0-xf, Par[1])*(1.0+Par[2]*xf)*std::pow(xf, -Par[3]);
      //
      switch (particle)
      {
          case kPiPlus: 
             return 0.0462498*std::pow(1.0-xf, 20.8206)*(1.0+11911.3*xf)*std::pow(xf, -1.0*0.96808);
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   else 
   {
      //
      //form: Par[4]*(1.0-xf)*std::exp(Par[5]*std::pow(xf, Par[6]));
      //
      switch (particle)
      {
          case kPiPlus: 
             return 1250.0*(1.0-xf)*std::exp(-13.6097*std::pow(xf, 0.597619));      
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }

   return 1.0;
   
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKPCB(const Lfluk::ParticleType_t particle, const double xf) const
{
   
   //
   //form: Par[0]*std::pow(1.0-xf, Par[1])*(1.0+Par[2]*xf)*std::pow(xf, -Par[3]);
   //
   
   switch (particle)
   {
       case kPiPlus: 
          return 3.31629*std::pow(1.0-xf, 2.48421)*(1.0+22.7285*xf)*std::pow(xf, -1.0*2.38778);
       case kPiMinus:
          break;
       case kKPlus:
          break;
       case kKMinus:
          break;
       default:
          break;
   }

   return 0.0;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKPCC(const Lfluk::ParticleType_t particle, const double xf) const
{
   if (xf<0.35)
   {
      //
      //form: Par[0]*std::pow(1.0-xf, Par[1])*(1.0+Par[2]*xf)*std::pow(xf, -Par[3]);
      //
      switch (particle)
      {
          case kPiPlus: 
              return 0.074085*std::pow(1.0-xf, 11.2413)*(1.0+168369.0*xf)*std::pow(xf, -1.0*2.47553);
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   else 
   {
      //
      //form: Par[4]*(1.0-xf)*std::exp(Par[5]*xf);
      //
      switch (particle)
      {
          case kPiPlus: 
              return 49071.3*(1.0-xf)*std::exp(-12.1073*xf);  
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   
   
   return 0.0;
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKPCD(const Lfluk::ParticleType_t particle, const double xf) const
{
   if (xf<0.1)
   {
      //
      //form: Par[0]*std::exp(Par[1]*std::pow(xf, Par[2]));
      //
      switch (particle)
      {
          case kPiPlus: 
              return 35.8788*std::exp(-1.95477*std::pow(xf, 0.0638124));
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   else 
   {
      //
      //form: Par[3] + par[4]*xf + par[5]*xf*xf +par[6]*xf*xf*xf
      //
      switch (particle)
      {
          case kPiPlus: 
             return 7.47623 + -7.35376*xf + -14.6949*xf*xf + 22.7689*xf*xf*xf; 
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   
   
   return 0.0;
   
}

//------------------------------------------------------------------------------------------
double Lfluk::GetSKPCE(const Lfluk::ParticleType_t particle, const double xf) const
{
   if (xf<0.25)
   {
      //
      //form: Par[0]*std::exp(Par[1]*std::pow(xf, Par[2]))*(Par{3] + Par[4]*xf + Par[5]*xf*xf);
      //
      switch (particle)
      {
          case kPiPlus: 
              return 3.88469*std::exp(-6.20932*std::pow(xf, 0.26879))*(1.29961 + 53.0232*xf + 121.575*xf*xf);
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   else 
   {
      //
      //form: (Par[6]+Par[7]*xf)/std::pow(xf, Par[8]);
      //
      switch (particle)
      {
          case kPiPlus: 
             return (-0.26414 + 1.87321*xf)/std::pow(xf, 1.28207);
          case kPiMinus:
             break;
          case kKPlus:
             break;
          case kKMinus:
             break;
          default:
             break;
      }
   }
   
   
   return 0.0;
   
}

---