Kfluk.cxx

Go to the documentation of this file.

#include <cassert>
#include<iostream>
#include<fstream>
#include "TSystem.h"
#include "MCReweight/Kfluk.h"
#include "MCNtuple/NtpMCFluxInfo.h"
#include "MessageService/MsgService.h"

CVSID("$Id: Kfluk.cxx,v 1.3 2006/12/03 02:31:45 gmieg Exp $");


using namespace std;

Kfluk::Kfluk(const char* pathname):
  ELO(),
  EHI(),
  CK0E3(),
  CK0M3(),
  CKE3(),
  CKM3()
{

  //initialize vectors by reading in vector files
  //copied code from Zbeam
  std::string topDir=pathname; // 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 after all those shenanigans, it still can't find the file, give up
    if(base=="") {
      MSG("MCReweight",Msg::kFatal)<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
      assert(false);
    }
    topDir = base+ "/" + topDir;
  }

  MSG("Kfluk",Msg::kDebug) <<"Kfluk reading files from: "<<topDir<<std::endl;

  float low, high, val;
  ifstream in1((topDir+"/k0m3.vec").c_str());
  if(in1.is_open()){
    in1>>low;
    while(!in1.eof()){
      in1>>high>>val;
      ELO.push_back(low);
      EHI.push_back(high);
      CK0M3.push_back(val);
      //      cout<<"high "<<high<<" low "<<low<<" val "<<val<<endl;
      in1>>low;
    }
    in1.close();
  }
  else{
    MSG("Kfluk",Msg::kError)<<"Could not open "<<topDir+"/k0m3.vec"<<endl;
  }

  ifstream in2((topDir+"/k0e3.vec").c_str());
  if(in2.is_open()){
    in2>>low;
    while(!in2.eof()){
      in2>>high>>val;
      CK0E3.push_back(val);
      in2>>low;
    }
    in2.close();
  }
  else{
    MSG("Kfluk",Msg::kError)<<"Could not open "<<topDir+"/k0e3.vec"<<endl;
  }

  ifstream in3((topDir+"/km3.vec").c_str());
  if(in3.is_open()){
    in3>>low;
    while(!in3.eof()){
      in3>>high>>val;
      CKM3.push_back(val);
      in3>>low;
    }
    in3.close();
  }
  else{
    MSG("Kfluk",Msg::kError)<<"Could not open "<<topDir+"/km3.vec"<<endl;
  }


  ifstream in4((topDir+"/ke3.vec").c_str());
  if(in4.is_open()){
    in4>>low;
    while(!in4.eof()){
      in4>>high>>val;
      CKE3.push_back(val);
      in4>>low;
    }
    in4.close();
  }
  else{
    MSG("Kfluk",Msg::kError)<<"Could not open "<<topDir+"/ke3.vec"<<endl;
  }

  bool sizematch=true;
  if(ELO.size()!=EHI.size()) sizematch=false;
  if(ELO.size()!=CK0E3.size()) sizematch=false;
  if(ELO.size()!=CK0M3.size()) sizematch=false;
  if(ELO.size()!=CKE3.size()) sizematch=false;
  if(ELO.size()!=CKM3.size()) sizematch=false;

  if(!sizematch){
    MSG("Kfluk",Msg::kError) <<"Sizes of vectors do not match!"
                             <<"Can not continue"<<std::endl;
    assert(false);
  }

}

Kfluk::~Kfluk()
{}

double Kfluk::GetWeight(NtpMCFluxInfo *fluxinfo)
{
  return GetWeight(fluxinfo->ndecay,fluxinfo->necm);
}

double Kfluk::GetWeight(int ndecay, float ecm)
{
  //C implementation of D. Jaffe's correction to kl3 decays
  //an approximation of the real fix that is going into gnumi
  //as of April, 2006.  This function is meant to be a stop gap measure.

  // correction factor for kl3 decays to take
  // matrix element into account (approximately)
  // see relflux#k3corr 27mar06
  //
  // 28mar06 add correction factors for BR(recent)/BR(gnumi)
  //        for the 2 Ke3 decays
  // KL: PRL93 (2004) 181802
  // K+: PRL91 (2003) 261802

  double corrk3 = 1.;
  if(ndecay==5) return corrk3;  //Km2 +
  if(ndecay==8) return corrk3;  //Km2 -
  if(ndecay>=11) return corrk3; //mu,pi decays
                                                        
  unsigned int i = 0;
  while(ecm>EHI[i]&&i<EHI.size()){
    i = i+1;
  }
  if(i>=EHI.size()){
    i = EHI.size()-1;
  }
  //  cout<<"Will look in "<<i<<"bin"<<endl;
  //  cout<<" factors are: "<<CK0E3[i]<<" "<<CK0M3[i]<<" "<<CKE3[i]<<" "<<CKM3[i]<<endl;

  if(ndecay==1||ndecay==2){
    corrk3 = CK0E3[i]*(40.67/38.70);
  }
  else if(ndecay==3||ndecay==4){
    corrk3 = CK0M3[i];
  }
  else if(ndecay==6||ndecay==9){
    corrk3 = CKE3[i]*(5.13/4.82);
  }
  else if(ndecay==7||ndecay==10){
    corrk3 = CKM3[i];
  }

  return corrk3;

}

---