BMPTHistoWeightCalculator.cxx

Go to the documentation of this file.

#ifndef bmpthistoweightcalculator_cxx
#define bmpthistoweightcalculator_cxx
#include "MessageService/MsgService.h"
#include "MCReweight/BMPTHistoWeightCalculator.h"
#include "TMath.h"

CVSID("$Id: BMPTHistoWeightCalculator.cxx,v 1.7 2007/01/15 20:20:26 rhatcher Exp $");

//*********************************************
BMPTHistoWeightCalculator::BMPTHistoWeightCalculator(Registry *config)
  : WeightCalculator(config)
{

  fWCname = "BMPTHistoWeightCalculator";

  fStandardConfigChanged = false;
  fReweightConfigChanged = false;

  if(config) SetStandardConfig(config);
  
  // 1 = std -> random, 0 = off, -1 = random -> std
  fUseBMPTHistoReweighting = 0; 
  fBMPTHistofile = 0;  
  fPidHist.clear();

  OpenReweightFile();

}

//*********************************************
void BMPTHistoWeightCalculator::OpenReweightFile()
{

  //get directory from env variable
  string directory = "";
  char *dir = getenv("GMINOSAUX");
  if(dir!=NULL) directory = dir;
  else directory = ".";

  string file = "/bmpt-wght.root";
  string fullpath = directory+file;

  //open file:  
  fBMPTHistofile = new TFile(fullpath.c_str(),"READ");

  //check file is opened:  
  if(!fBMPTHistofile->IsOpen() || fBMPTHistofile->IsZombie()) {
    MSG("BMPTHistoWC",Msg::kError) 
      << "Could not open BMPTHisto reweight file!\n" 
      << "Check that $GMINOSAUX points to the correct directory.\n" 
      << "All weights will be 1" << endl;
    delete fBMPTHistofile;
    fBMPTHistofile = 0;
    return;
  }

  //get histograms from file:
  TH2F *piplushist  = (TH2F*) fBMPTHistofile->Get("h9901008");
  TH2F *piminushist = (TH2F*) fBMPTHistofile->Get("h9901009");
  TH2F *kplushist   = (TH2F*) fBMPTHistofile->Get("h9901011");
  TH2F *kminushist  = (TH2F*) fBMPTHistofile->Get("h9901012");
  TH2F *k0longhist  = (TH2F*) fBMPTHistofile->Get("h9901010");
  TH2F *k0shorthist = (TH2F*) fBMPTHistofile->Get("h9901016");
  TH2F *phist       = (TH2F*) fBMPTHistofile->Get("h9901014");
  TH2F *pbarhist    = (TH2F*) fBMPTHistofile->Get("h9901015");

  //add histograms to map with appropriate PID:
  fPidHist[8]  = piplushist;
  fPidHist[9]  = piminushist;
  fPidHist[11] = kplushist;
  fPidHist[12] = kminushist;
  fPidHist[10] = k0longhist;
  fPidHist[16] = k0shorthist;
  fPidHist[14] = phist;
  fPidHist[15] = pbarhist;

}

//*********************************************
BMPTHistoWeightCalculator::~BMPTHistoWeightCalculator()
{
  
  if(fBMPTHistofile) {
    fBMPTHistofile->Close();
    delete fBMPTHistofile; 
  }
  
}

//*********************************************
void BMPTHistoWeightCalculator::ReweightConfigReset()
{
  if(fStandardConfig) SetReweightConfig(fStandardConfig);
}

//*********************************************
void BMPTHistoWeightCalculator::Config()
{

  Int_t tempi = 0;
  if(fStandardConfigChanged) {
    if(fStandardConfig->Get("bmpt_histo:use",tempi)) 
      fUseBMPTHistoReweighting = tempi;
    fStandardConfigChanged = false;
  }
  if(fReweightConfigChanged) {
    if(fReweightConfig->Get("bmpt_histo:use",tempi)) 
      fUseBMPTHistoReweighting = tempi;
    fReweightConfigChanged = false;
  }
}

//*********************************************
double BMPTHistoWeightCalculator::GetWeight(MCEventInfo *event,NuParent *parent)
{
  if(!event) return 1;
  if(!parent) return 1;
  return 1;
}

//*********************************************
double BMPTHistoWeightCalculator::GetWeight(Registry *event)
{

  if(!fBMPTHistofile) return 1;
  if(!abs(fUseBMPTHistoReweighting)==1) return 1;
  if(!event) return 1;

  //neutrino parent characteristics:
  double px = 0;
  double py = 0;
  double pz = 0;
  int pid = 0;
  if(!event->Get("event:nuparent_px",px)) return 1;
  if(!event->Get("event:nuparent_py",py)) return 1;
  if(!event->Get("event:nuparent_pz",pz)) return 1;
  if(!event->Get("event:nuparent_pid",pid)) return 1;

  //get the appropriate histogram  
  std::map<int,TH2F*>::iterator pidHistIter = fPidHist.find(pid);
  if( pidHistIter == fPidHist.end() ) {
    MSG("BMPTHistoWC",Msg::kWarning) << "Unknown PID: " << pid << endl;
    return 1;
  }

  TH2F *hist = pidHistIter->second;
  if(!hist) return 1;

  Double_t log_pz = TMath::Log10(pz);
  Double_t log_pt = 0.5*TMath::Log10(px*px + py*py);

  Int_t binx = hist->GetXaxis()->FindBin(log_pz);
  Int_t biny = hist->GetYaxis()->FindBin(log_pt);

  if(binx>hist->GetNbinsX()) binx = hist->GetNbinsX();
  if(binx<1) binx = 1;
  if(biny>hist->GetNbinsY()) biny = hist->GetNbinsY();
  if(biny<1) biny = 1;

  double weight = hist->GetBinContent(binx,biny);
  if(fUseBMPTHistoReweighting==-1) return 1./weight; 
  return weight;

}

//*********************************************
void BMPTHistoWeightCalculator::PrintReweightConfig(ostream & stream)
{
  if(!fUseBMPTHistoReweighting) stream << "Not ";
  stream << "Using BMPT Histo Reweighting" << endl;
}

#endif

---