BeamEnergyCalculator.cxx

Go to the documentation of this file.

// BeamEnergyCalculator
// 
// Calculates weights for the different beam energy configurations
// as well as a 1/E flux.  
//
//
// Created: M. Kordosky  -- June 2, 2005
//
// $Author: rhatcher $ 
//
// $Revision: 1.6 $
// 
// $Name:  $
//
// $Id: BeamEnergyCalculator.cxx,v 1.6 2007/01/30 03:13:14 rhatcher Exp $
//

#include "MCReweight/BeamEnergyCalculator.h"
#include <cmath>
//#include <string>
#include <ctype.h>
#include <cassert>
#include <utility>
#include <iterator>
#include <iostream>

#include "TH1.h"
#include "TFile.h"
#include "TSystem.h"

#include "Registry/Registry.h"
#include "MCReweight/MCEventInfo.h"
#include "MCReweight/NuParent.h"
#include "MessageService/MsgService.h"


CVSID("$Id: BeamEnergyCalculator.cxx,v 1.6 2007/01/30 03:13:14 rhatcher Exp $");
const double BeamEnergyCalculator::low_energy_cut=1e-6;

BeamEnergyCalculator::BeamEnergyCalculator(Registry* stdconfig) 
   : 
  WeightCalculator(stdconfig),
  fFluxDir(0),fFluxFileName(""),
  fFluxFileChanged(true),
  fDetector(Detector::kUnknown),
  fBeam(BeamType::kUnknown),
  fLowRange(low_energy_cut), fHighRange(120.0),fCache()
{
  fWCname="BeamEnergyCalculator";
  fStandardConfigChanged=true;
  Config();
}


void BeamEnergyCalculator::Config() {
  MSG("BEC", Msg::kInfo)<<"configuring..."<<endl;

  int tempi=-1;
  double tempd=0.0;
  std::string ff;
  if(fStandardConfigChanged){
    // this is a little clumsy
    const char* c=fStandardConfig->GetCharString("beam:flux_file");    
    if(c){
      ff=c;
      if(ff!=fFluxFileName){
        fFluxFileName=ff;
        fFluxFileChanged=true;
      }
    }
    /*    
    if(fStandardConfig->Get("beam:nu_id",tempi)){
      fNuID=tempi;
    }
    if(fStandardConfig->Get("beam:nu_energy",tempd)){
      fNuEnergy=tempd;
    }
    */

    if(fStandardConfig->Get("beam:beam_type",tempi)){
      fBeam=static_cast<BeamType::BeamType_t>(tempi);
    }
    if(fStandardConfig->Get("beam:detector",tempi)){
      fDetector=static_cast<Detector::Detector_t>(tempi);
    }    
    if(fStandardConfig->Get("beam:low_energy_limit",tempd)){
      fLowRange=tempd;
    }
    if(fStandardConfig->Get("beam:high_energy_limit",tempd)){
      fHighRange=tempd;
    }
    
    fStandardConfigChanged=0;
  }

}

void BeamEnergyCalculator::ReweightConfigReset(){
  if(fStandardConfig) SetReweightConfig(fStandardConfig);
}


double BeamEnergyCalculator::GetWeight(Registry* r){
  
  BeamType::BeamType_t bt=BeamType::kUnknown;
  Detector::Detector_t det = Detector::kUnknown;
  int nu_idhep=0;
  double nu_energy=0;
  double high=0;
  double low=0;
  double w=1.0;

  // these must be passed in
  bool bad=false;
  if(!r->Get("beam:nu_id",nu_idhep)) bad=true;
  if(!r->Get("beam:nu_energy",nu_energy)) bad=true;

  if(bad){
    MSG("BEC", Msg::kFatal)<<"neutrino id or energy not set"<<endl;
    assert(0);
  }
  // look for certain values set in the registry
  // if not found, use defaults
  int tmpi=0;
  if(!r->Get("beam:beam_type",tmpi)){
    //    if(!fStandardConfig->Get("beam:beam_type",bt)) return w;
    bt=fBeam;
  }
  else{
    bt=static_cast<BeamType::BeamType_t>(tmpi);
  }
  
  if(!r->Get("beam:detector",tmpi)){
    //    if(!fStandardConfig->Get("beam:detector",det)) return w;
    det=fDetector;
  }
  else{
    det=static_cast<Detector::Detector_t>(tmpi);
  }
  if(!r->Get("beam:high_energy_limit",high)){
    //    if(!fStandardConfig->Get("beam:high_energy_limit",high)) return w;
    high=fHighRange;
  }
  if(!r->Get("beam:low_energy_limit",low)){
    //    if(!fStandardConfig->Get("beam:low_energy_limit",low)) return w;
    low=fLowRange;
  }

  w=GetWeight(bt,det,nu_idhep,nu_energy,high,low);
  
  return w;
}
double BeamEnergyCalculator::GetWeight(MCEventInfo* /* event */,
                                       NuParent* /* parent */){
  MSG("BEC", Msg::kFatal)<<"Function Not Implemented Yet!"<<endl;
  assert(0);
  return 1;
}


double BeamEnergyCalculator::GetWeight(BeamType::BeamType_t bt, 
                                       Detector::Detector_t det,
                                       int nu_idhep, double nu_energy, 
                                       double high_energy_limit,
                                       double low_energy_limit){
  double w=0.0;
  if(bt==BeamType::kInverseE){
    // special case of 1/E flux
    w=DeriveInverseEWeight(nu_energy,high_energy_limit, low_energy_limit);
  }
  else if(bt==BeamType::kUnknown){
    MSG("BEC", Msg::kFatal)<<"Requested weight for unknown beam type"<<endl;
    assert(0);
  }
  else{
    TH1* h = GetHist(bt,det,nu_idhep);
    
    if(!h) {w=1.0; return w;}
    
    w=DeriveWeight(h,nu_energy,high_energy_limit,low_energy_limit);
  }

  return w;
}
  
double BeamEnergyCalculator::DeriveWeight(TH1* h,double energy,
                                          double high, double low){
  // simply uses the binned histogram to estimate the weights
  // certainly, it's reliant on the quality of the histogram...

  if(energy>high || energy<low) return 0;
  if(!h) return 0;
  if(low<low_energy_cut) low=low_energy_cut;
  Double_t* integral=h->GetIntegral();
  if(!integral){ h->ComputeIntegral(); integral=h->GetIntegral();}
  Int_t binlow=h->FindBin(low);
  Int_t binhigh=h->FindBin(high);
  if(binlow<1 || binhigh>h->GetNbinsX()) return 0;
  Int_t bin = h->FindBin(energy);
  double weight = h->GetBinContent(bin);
  double norm = integral[binhigh]-integral[binlow-1];
  Stat_t stats[4];
  h->GetStats(stats);
  Stat_t sumw=stats[0];
  /*
  static bool first=true;
  if(first){
    std::cout<<"Histogram: "<<std::endl;
    h->Print("all");
    for(int i=1; i<=h->GetNbinsX(); i++){
      std::cout<<"integral["<<i<<"] = "<<integral[i]<<std::endl;
    }
    std::cout<<"sumw: "<<sumw<<std::endl;
    first=false;
  }
  */
  if(norm<=0) {
    MSG("BEC",Msg::kFatal)<<"normalization=0"<<endl;
    assert(0);
  }

  weight/=(norm*sumw);

  return weight;
}


double BeamEnergyCalculator::DeriveInverseEWeight(double energy,
                                                  double high,double low){
  // 1/E flux over range [low,high]
  // range is necessary for proper normalization of weights!
  if(low<low_energy_cut) low = low_energy_cut;
  if(energy<low || energy>high) return 0;
  return (1/energy)/std::log(high/low);

}

void BeamEnergyCalculator::ConstructName(Detector::Detector_t det,
                                         int nu_idhep, std::string& s){
  // given the detector and type of neutrino
  // construct (part of) a histogram name
  // note, adds to 's'
  
  using std::string;
  string nu;
  switch(nu_idhep){
  case 12: nu="nue"; break;
  case -12: nu="nuebar"; break;
  case 14: nu="numu"; break;
  case -14: nu="numubar"; break;
  case 16: nu="nutau"; break;
  case -16: nu="nutaubar"; break;
  default: nu="???"; break;
  }
  string dets=Detector::AsString(det);
  for(string::size_type i=0; i<dets.size(); i++){
    dets[i]=std::tolower(dets[i]);
  }

  s+=dets+'_'+nu;

}

TH1* BeamEnergyCalculator::GetHist(BeamType::BeamType_t bt, 
                                   Detector::Detector_t det, 
                                   int nu_idhep){
  // get a histogram from fFluxDir
  // histograms are taken from the file according to a naming convention
  // that encompasses the type of neutrino, Detector, and BeamType

  // open new file if needed
  if(fFluxFileChanged){
    
    if(fFluxDir){
      delete fFluxDir; fFluxDir=0;
    }
    fCache.clear();
    MSG("BEC",Msg::kInfo)<<"opening new flux file: "
                         <<fFluxFileName<<endl;
    // check that file exists
    if(gSystem->AccessPathName(fFluxFileName.c_str()) == kTRUE){
      // file *doesn't* exist... who established that convention!!
      MSG("BEC",Msg::kFatal)<<"cannot find flux file: "
                            <<fFluxFileName<<endl;
      assert(0); //die,die,die
    }
    else{
      fFluxDir = new TFile(fFluxFileName.c_str());
      MSG("BEC",Msg::kInfo)<<"successfully opened: "
                           <<fFluxFileName<<endl;
    }
    fFluxFileChanged=false;
  }

  // look for histogram in our cache
  // initially I just wanted to use the TFile for this purpose
  // however, if you do TH1* h = (TH1*) file->Get("some_dir/some_hist")
  // multiple times, the TFile code looks in the list of objects
  // for "some_dir/some_hist" after the second time,
  // but, of course, the actual thing in the list is called "some_hist"
  // so, it rereads the histogram from the file, which is real slow
  //  I couldn't find a workaround, so, hash_map..
  TH1* h = 0;

  // make the histogram name
  // beamtype tag is the directory name
  std::string s=BeamType::AsTag(bt);
  // pick a flux histogram
  // there are also 'evt' histograms but
  // they don't play well with 1/E spectrum.. must know xsection
  s+="/flux_";
  ConstructName(det,nu_idhep,s);
  // so, for LE beam, muon-neutrinos in the near detector
  // s=z_000/flux_near_numu


  CacheMapIter it = fCache.find(s);
  if(it!=fCache.end()){
    h=it->second;
  }
  else{
    // if not found, look in the file    
    //TH1* h = FindHist(fFluxDir,bt,det,nu_idhep);
    h = dynamic_cast<TH1*>(fFluxDir->Get(s.c_str()));
    if(h){
      //      fCache[s]=h; // add it!
      fCache.insert(std::pair<std::string,TH1*>(s,h));
    }
    else{
      MAXMSG("BEC",Msg::kWarning,10)
        <<"Could not find a flux histogram with the tag: "
        <<s<<"\n It's possible to continue but"
        <<" these events will not be weighted properly."<<endl;      
    }
  }
  return h;
}

---