LZfluk.cxx

Go to the documentation of this file.

// $Author: gmieg $ 
//
// $Revision: 1.3 $
// 
// $Name:  $
//
// $Id: LZfluk.cxx,v 1.3 2009/03/17 18:10:06 gmieg Exp $
//
// $Log: LZfluk.cxx,v $
// Revision 1.3  2009/03/17 18:10:06  gmieg
// #include <cstring> for strcmp.
//
// Revision 1.2  2009/03/08 19:30:07  loiacono
// commit changes for Muon Fitting
//
// Revision 1.1  2009/02/15 14:23:27  loiacono
// commit new classes for muon monitor fitting
//
// Revision 1.12  2008/03/15 17:23:31  zarko
//
// Renamed the old GetWeight functions in Zfluk to
// GetWeightPRL and GetWeightBoston to avoid confusion between
// different functional forms used for those fits.
// Updated SKZPWeightCalculator to reflect the changes.
//
// Revision 1.11  2007/12/20 21:29:14  zarko
// *** empty log message ***
//
// 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/LZfluk.h"
#include "MessageService/MsgService.h"
#include "Registry/Registry.h"

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

CVSID("$Id: LZfluk.cxx,v 1.3 2009/03/17 18:10:06 gmieg Exp $");

LZfluk::LZfluk(/*const char* dir*/) 
   :fPar(0),
    fWghtFromHist(true)
{

}

//------------------------------------------------------------------------------------------
void LZfluk::Init()
{
   cout << "LZfluk::Init" << endl;

    std::string dir = "";
    TH1::AddDirectory(kFALSE);
   
  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
        gSystem->FreeDirectory(dir_ptr);
      }
    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("LZfluk",Msg::kDebug) <<"LZfluk reading data from: "<<topDir<<std::endl;
  std::string fileName =topDir+"/fluka05ptxf.root";

  TFile* hFile= 0;

  cout << "Opening ptxf file" << endl;

  hFile = new TFile(fileName.c_str());

  if(!hFile || !(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;
     }
  }
  
  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 = hist -> GetNbinsX();
      const int nbinsx = 120;
      const int nbinsy = hist -> GetNbinsY();
      //cout << "nbinsx = " << nbinsx << "nbinsy = " << nbinsy << endl;
      //
      //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);
         }
      }
      fWeightHist[*itPlist]->SetDirectory(0);

      std::pair<LZfluk::ParticleType_t, TH2F* > p(*itPlist, hist);
      std::pair<LZfluk::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;
}

LZfluk::~LZfluk() {
 for (std::vector<ParticleType_t>::iterator itPlist=fPlist.begin();itPlist!=fPlist.end(); itPlist++)
    {
      delete fWeightHist[*itPlist];
      delete fPTXF[*itPlist];
      delete fWeightedPTXF[*itPlist];
    }
 fWeightHist.clear();
 fPTXF.clear();
 fWeightedPTXF.clear();
 fMeanPT.clear();
 fWeightedMeanPT.clear();
 fN.clear();
 fNWeighted.clear();
 fPlist.clear();
 fPar.clear();
 fNBinsX.clear();
 fNBinsY.clear();
}

//------------------------------------------------------------------------------------------
void LZfluk::Config(const Registry &reg)
{

   const char *value_char = 0;

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

        
   cout << "LZfluk::Config" << endl;
   cout << "   Get Weight From Hist = " <<  fWghtFromHist << endl;
   
   
}

//-------------------------------------------------------------------------------------
double LZfluk::GetWeight(int ptype, double pT, double pz)
{
    
   if(fWghtFromHist)
   {
      
      double weight = 1.0;
      
      LZfluk::ParticleType_t eptype=GetParticleEnum(ptype);
      
      if(eptype != kPiPlus && eptype != kKPlus && eptype != kPiMinus && eptype != kPiMinus)
      {
         return 1.0;
      }
      
      if (pz > 120.0 || pz < 0.0 || pT < 0.0 || pT > 1.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(!(binx < 1) && !(biny < 1) && !(binx > nbinsx) && !(biny > nbinsy))
      {
         weight = fWeightHist[eptype] -> GetBinContent(binx, biny);
      }
      
      return weight;
   }

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

//-------------------------------------------------------------------------------------
double LZfluk::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 "
        <<"LZfluk::CalcWeight (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;
  LZfluk::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=CalcWeight(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=CalcWeight(ptype,0.051,pz);  // for low pT
    }
  else if (eptype==kPiMinus)
    {
      weight=CalcWeight(kPiPlus,pT,pz);
      if ( pz > 4. ) weight *= ( fPar[12] + fPar[13] * xF );
    }
  else if (eptype==kKMinus)
    {
      weight=CalcWeight(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 LZfluk::RecalculateWeights()
{
  
  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));

           
           fWeightHist[*itPlist]->SetBinContent(i,j,weight);
           if(weight - fWeightHist[*itPlist]->GetBinContent(i,j) != 0)
           {
              LZfluk::ParticleType_t eptype= *itPlist;
              string name = GetParticleName(eptype);
              cout << "LZfluk::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);
           

        }
     }
     meanpt/=sum;
     fWeightedMeanPT[*itPlist]=meanpt*1000.; //GeV to MeV
     fNWeighted[*itPlist]=sum;
     meanpt=0.;
     sum=0.;
     
     
  }
}

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

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

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

}

---