Zbeam.cxx

Go to the documentation of this file.

// $Author: med $ 
//
// $Revision: 1.29 $
// 
// $Name:  $
//
// $Id: Zbeam.cxx,v 1.29 2009/10/31 20:42:08 med Exp $
//
// $Log: Zbeam.cxx,v $
// Revision 1.29  2009/10/31 20:42:08  med
// Change kInfo to kFatal for non-existant beamsys_XYZ root file and change check to IsZombie()
//
// Revision 1.28  2008/02/22 16:31:08  rearmstr
// Added an option in BeamSys.h, Zbeam, and SKZPWeightCalculator.h
// to get the total error from focusing.
//
// Revision 1.27  2007/11/11 07:22:47  rhatcher
// Part of the ongoing purge of "DetectorType" in favor of "Detector"
// (the former is a synonym of the latter ...)
//
// Revision 1.26  2007/06/06 15:14:45  vahle
// Make Zbeam::SetReweightConfig return the user to the TDirectory they started in
//
// Revision 1.25  2007/05/13 13:47:21  vahle
// Adding new histogram files for Zbeam systematics.  Old hist file only worked with bleeding edge root (TDirectoryFiles are written in the new root instead of TDirectories).  Updated Zbeam to read in 4 separate files instead of 1 file with 4 directories.  Added different run periods to SKZPWeightCalculator to apply different weights depending on supplied run period.  Default is runI which does the same thing as it always used to.  If RunI is specified, weights are applied to ME and HE to correct for beam spot size effects.  If RunII is specified, LE010 beam is weighted to le9 beam.  Also supplied a function to supply a total weight which is the exposure weighted average of the two run periods.  Finally, the macro TestSKZPWeight has been updated to demonstrate the use of all these new features.
//
// Revision 1.24  2007/05/03 16:35:05  zarko
// Added nue hadprod error.
//
// Revision 1.23  2007/04/28 19:46:12  zarko
// *** empty log message ***
//
// Revision 1.22  2007/04/26 16:23:41  zarko
// *** empty log message ***
//
// Revision 1.21  2007/04/24 21:54:17  zarko
// Added BeamSys enumerator. Made Zbeam read from one root file
// instead of infinite amount of .vec files.
//
// Revision 1.20  2007/04/14 03:32:59  zarko
// Updated Zbeam to include target z correction for le10 (effect=15) and added the
// PiMinus_CedarDaikon hadron production errors.
// note: The target z correction takes -1 to give the weights for the le9 beam
//
// Revision 1.19  2007/03/16 17:05:15  vahle
// Added PiMinus weighting code to SKZPWeightCalculator
//
// Revision 1.18  2007/03/04 07:06:49  zarko
// Updated code to return latest hadron production error envelope.
//
// Revision 1.17  2007/01/10 16:02:15  zarko
// Renamed "PiMinus" config into "Boston".
//
// Revision 1.16  2007/01/05 23:21:49  zarko
// Added hadron production error envelope for Boston fit (this includes pi-).
// Use SetReweightConfig(config) to set either PRL or Boston fit
// hadron production error envelope (config="PRL" or "PiMinus")
//
// Revision 1.15  2006/10/03 16:00:56  zarko
// Updated code to reweight numus and numubars
//
// Revision 1.14  2006/09/20 21:18:51  zarko
// Added le150z200i beam.
//
// Revision 1.13  2006/06/03 03:42:33  zarko
// *** empty log message ***
//
// Revision 1.12  2006/05/27 15:40:10  zarko
// Added horn off beam.
//
// Revision 1.11  2006/04/15 04:28:18  zarko
//
// Last commit contained one unintendend change, so this is correcting for that.
//
// Revision 1.10  2006/04/14 02:36:24  zarko
// Zbeam now looks into SRT_PUBLIC_CONTEXT if the data cannot be found in SRT_PRIVATE_CONTEXT
//
// Revision 1.9  2006/02/28 05:50:58  zarko
// Calling GetWeight with nEffect=10 or nEffect=11 will now return the total
// error (all effects summed in quadrature). nEffect=10 includes hadron production
// error before hadron reweighting and nEffect=11 includes hadron production
// error after reweighting.
//
// Revision 1.8  2006/02/26 15:35:59  zarko
// Added hadron productin error after reweighting the MC as nEffect==8.
// Added le010z170i and le010z200i beams as iBeam==6 and iBeam==7.
//
// Revision 1.7  2006/02/22 01:22:46  zarko
// Added the reweighting for low intensity LE10 (Ibeam=5)
// During this run beam was narrower than in current MC,
// so as for ME and HE setting ns=0 for beamwidth will
// reweight to narrower beam.
//
// Revision 1.6  2006/01/30 20:53:52  zarko
// Updated beam width error vectors for ME and HE
// Put back hadron production error as nEffect=7
//
// Revision 1.5  2006/01/29 00:44:26  zarko
// Added beam width effect for ME and HE.
// Removed hadron production from Zbeam.
//
// Revision 1.4  2006/01/08 10:05:13  gmieg
// "#include <cassert>" added so it will compile.
//
// Revision 1.3  2005/12/28 19:34:33  brebel
// set the "Zbeam reading files from: " message to kDebug level and set the
// stream to Zbeam
//
// Revision 1.2  2005/12/18 17:45:40  kordosky
// reacting to instructions from zarko, modify weight function to return 1+deviation (e.g. a weight).
//
// Revision 1.1  2005/12/18 16:36:21  kordosky
// Z. Pavlovic zbeam code added . see macros/test_zbeam.C
//

#include "MCReweight/Zbeam.h"
#include "MessageService/MsgService.h"
#include "Conventions/Detector.h"
#include "Conventions/BeamType.h"
#include "MCReweight/BeamSys.h"

#include <fstream>
#include <string>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cassert>

#include "TSystem.h"
#include "TFile.h"
#include "TH1D.h"

ClassImp(Zbeam)
  
CVSID("$Id: Zbeam.cxx,v 1.29 2009/10/31 20:42:08 med Exp $");

Zbeam::Zbeam(const char* dir)
{
  // fTopDir is where the beamsys_hists.root is
  fTopDir=dir; // user may set location of input data
  if(fTopDir=="") { // by default, this code looks in a standard place
    fTopDir="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 + "/" + fTopDir;
        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);
    }
    fTopDir = base+ "/" + fTopDir;
  }
  MSG("Zbeam",Msg::kDebug) <<"Zbeam getting hists from: "<<fTopDir<<"/beamsys_hists.root"<<std::endl;
   
  SetReweightConfig("PRL"); //default config
}

Zbeam::~Zbeam() {

}

void Zbeam::SetReweightConfig(std::string cnf)
{
  if ( cnf == fCurrentConfig) return;
  TDirectory *save = gDirectory;

  fBeamSysFile = new TFile(Form("%s/beamsys_%s.root",fTopDir.c_str(),cnf.c_str()));
  if (fBeamSysFile->IsZombie()) {
    MSG("Zbeam",Msg::kFatal)<<"Don't recognize "
                            <<cnf
                            <<" configuration."
                            <<" Config set to :"
                            <<fCurrentConfig<<endl;
    
    return;
  }
  
  fBeamSysMap.clear();
  fCurrentConfig = cnf;
  std::string nus[]={"NuMu","NuMuBar","NuE","NuEBar"};
  int ntype[]    ={    56,       55  , 53  , 52};
  for (int inu=0;inu<4;inu++) {
    for (int eff=BeamSys::kFirst;eff<BeamSys::kEndOfList;eff++) {
      for (int beam=BeamType::kUnknown;beam<BeamType::kEndOfList;beam++) {
        for (int det=Detector::kNear;det<=Detector::kFar;det++) {
          std::string hname=nus[inu]+"_"+
            BeamSys::AsString(BeamSys::BeamSys_t(eff))+"_"+
            BeamType::AsString(BeamType::BeamType_t(beam))+"_"+
            Detector::AsString(Detector::Detector_t(det));
          TH1D* hist=dynamic_cast<TH1D*> (fBeamSysFile->Get(hname.c_str()));
          //      cout<<eff<<" hname is "<<hname<<endl;
          if (hist) FillVector(hist,ntype[inu],eff,beam,det);
        }
        //for far/near error bar internally use Detector::kUnknown
        std::string hname=nus[inu]+"_"+
          BeamSys::AsString(BeamSys::BeamSys_t(eff))+"_"+
          BeamType::AsString(BeamType::BeamType_t(beam))+"_FN";
        TH1D* hist=dynamic_cast<TH1D*> (fBeamSysFile->Get(hname.c_str()));
        if (hist) FillVector(hist,ntype[inu],eff,beam,Detector::kUnknown);
      }
    }
  }
  
  fBeamSysFile->Close();
  delete fBeamSysFile;
  fBeamSysFile = 0;
  save->cd();
}

double Zbeam::GetWeight(ZbeamData_t zbmdata, BeamSys::BeamSys_t sys, double nSigma)
{
  int ntype  = zbmdata.ntype;
  double Enu = zbmdata.true_enu;
  Detector::Detector_t det = zbmdata.detector;
  BeamType::BeamType_t beam = zbmdata.beam;

  int bin;
  bin=int(2*Enu);
  if (bin<0) bin=0;
  if (bin>99) bin=99;

  if (ntype == 14) ntype = 56;
  else if (ntype == -14) ntype = 55;
  else if (ntype == 12) ntype = 53;
  else if (ntype == -12) ntype = 52;

  double weight = 1. ;
  double totError=0.;
  BeamSys::BeamSys_t first=sys; 
  BeamSys::BeamSys_t last=sys;
  if (sys == BeamSys::kTotalBefore || sys == BeamSys::kTotalAfter) //Total Error
    {
      first = BeamSys::kFirst;
      last = BeamSys::kEndOfList; 
    }
  else if(sys == BeamSys::kTotalFocusing)
    {
      first = BeamSys::kFirst;
      last = BeamSys::kHornIDist;
    } 

  std::map<int, DetMap_t>::iterator bsmit = fBeamSysMap.find(ntype);   
  DetMap_t::iterator dmit;
  BeamMap_t::iterator bmit;
  bool beam_it_exists = true;
  if (bsmit == fBeamSysMap.end()) {
    weight = 0.; 
    beam_it_exists = false;
  } else {
    dmit = (bsmit->second).find(det);
    if (dmit == (bsmit->second).end() ) {
      weight = 0.; 
      beam_it_exists = false;
    } else {
      bmit = (dmit->second).find(beam);
      if (bmit == (dmit->second).end() ) {
        weight = 0.; 
        beam_it_exists = false;
      } 
    }
  }

  if (beam_it_exists) {
    for (int eff=first;eff<=last;eff++) {    
      if (BeamSys::EBeamSys(eff) != BeamSys::kTotalBefore && 
          BeamSys::EBeamSys(eff) != BeamSys::kTotalAfter &&
          BeamSys::EBeamSys(eff) != BeamSys::kHadProdBefore && 
          BeamSys::EBeamSys(eff) != BeamSys::kHadProdAfter &&
          BeamSys::EBeamSys(eff) != BeamSys::kTotalFocusing )
        {
          EffMap_t::iterator effit = (bmit->second).find(eff);
          weight =  ( effit == (bmit->second).end() ) ?  0. : ((*effit).second)[bin];
          totError += weight*weight;
        }
    } // loop over all effects 
    
    if ((sys == BeamSys::kTotalBefore || sys == BeamSys::kHadProdBefore )) {
      EffMap_t::iterator effit = (bmit->second).find(BeamSys::kHadProdBefore);
      weight =  ( effit == (bmit->second).end() ) ?  0. : ((*effit).second)[bin];
      totError += weight*weight;
    }
    if ((sys == BeamSys::kTotalAfter || sys == BeamSys::kHadProdAfter )) {
      EffMap_t::iterator effit = (bmit->second).find(BeamSys::kHadProdAfter);
      weight =  ( effit == (bmit->second).end() ) ?  0. : ((*effit).second)[bin];
      totError += weight*weight;
    }
  }

  totError=sqrt(totError);
  
  if (!(sys == BeamSys::kTotalBefore || sys == BeamSys::kTotalAfter 
        || sys == BeamSys::kTotalFocusing))  totError *= ( (weight>0.) - (weight<0.) );
  
  double nS = (sys == BeamSys::kBeamWidth) ? nSigma + 1. : nSigma; // by default reweights pME, pHE and lowint 

  //some of the effects are symmetric and only abs(nS) makes sense
  totError = ((sys == BeamSys::kHorn1Offset) ||
              (sys == BeamSys::kBaffleScraping)) ? totError*fabs(nS)+1.0 :  totError*nS+1.; 
 
  return totError;  
} 

double Zbeam::GetWeight(ZbeamData_t zbmdata,std::map<BeamSys::BeamSys_t, double> nSigma)
{ 
  double wgh=1.;
  std::map<BeamSys::BeamSys_t, double>::const_iterator bsit;

  for (bsit=nSigma.begin();bsit != nSigma.end() ;bsit++)
    {
      wgh*=GetWeight(zbmdata,bsit->first, bsit->second); 
    }
  return wgh; 
} 

double Zbeam::GetWeight(int iDet, int iBeam, int nEffect,  
                        double nSigma, int ntype, double Enu) 
{
  ZbeamData_t zbm;
  zbm.beam = BeamType::FromZarko(iBeam);
  zbm.detector = (iDet != 3) ?  Detector::EDetector(iDet) : Detector::kUnknown;
  zbm.ntype = ntype;
  zbm.true_enu = Enu; 
  BeamSys::BeamSys_t sys = BeamSys::EBeamSys(nEffect);

  return GetWeight(zbm,sys, nSigma);
}

double Zbeam::GetWeight(int iDet, int iBeam,  
                        std::vector<double> nSigma, int ntype, double Enu) 
{
  ZbeamData_t zbm;
  zbm.beam = BeamType::FromZarko(iBeam);
  zbm.detector = (iDet != 3) ?  Detector::EDetector(iDet) : Detector::kUnknown;
  zbm.ntype = ntype;
  zbm.true_enu = Enu;

  std::map<BeamSys::BeamSys_t, double> sgm;
  for (unsigned int nEffect=0; nEffect<nSigma.size();nEffect++)
    {
      BeamSys::BeamSys_t sys = BeamSys::EBeamSys(nEffect+1);
      const std::pair<BeamSys::BeamSys_t, double> sp(sys, nSigma[nEffect]);
      sgm.insert(sp);
    }
  return GetWeight(zbm, sgm);
}

void Zbeam::MakeHistogram(std::string name, std::vector<double> vec)
{
  //Function creates a histogram using a vector vec (assumes 0.5GeV bins)
  //Check if the name contains the detector, beam and effect
  //This doesn't guarantee that the name is correct. 
  //The name should look like this NuMu_Horn1Offset_L010z185i_Near ...
  bool validname=false;
  std::string numu="NuMu";
  std::string numubar="NuMuBar";
  std::string nue="NuE";
  std::string nuebar="NuEBar";
  
  for (int det=Detector::kNear;det<=Detector::kFar;det++) {
    for (int beam=BeamType::kUnknown;beam<BeamType::kEndOfList;beam++) {
      for (int eff=BeamSys::kFirst;eff<BeamSys::kEndOfList;eff++) {
        if (((name.find(Detector::AsString(Detector::Detector_t(det))) != string::npos)||name.find("FN") != string::npos) &&
            (name.find(BeamType::AsString(BeamType::BeamType_t(beam))) != string::npos) &&
            (name.find(BeamSys::AsString(BeamSys::BeamSys_t(eff))) != string::npos) &&
            ((name.find(numu) != string::npos) || (name.find(numubar) != string::npos) ||
             (name.find(nue) != string::npos) || (name.find(nuebar) != string::npos))) {
          validname=true;
        } 
      }
    }
  }
  if (!validname) {
    cout<<"Invalid histogram name! "<<name <<endl;
    cout<<"Name should include beam (as in BeamType), detector (as in Detector)" 
        <<" and beam systematics type"<<endl;
    return;
  }
  cout<<"Creating "<<name<<" histogram"<<endl;
  unsigned int nbins=vec.size()-1;
  double maxx=0.5*nbins;

  TH1D* hist=new TH1D(name.c_str(),name.c_str(),nbins,0.,maxx);
  for (unsigned int i=1;i<nbins+1;i++)
    hist->SetBinContent(i,vec[i-1]);
  
  fHistMap[name]=hist;
}

void Zbeam::SaveHistogram(std::string name, std::string config, std::string filename)
{
  map<std::string, TH1D* >::iterator it = fHistMap.find(name);
  if ( it == fHistMap.end() ) {
    cout << name <<" does not exist!"<<endl;
    return;
  }
  TFile* f=new TFile(filename.c_str(),"UPDATE");
  if (!f->GetListOfKeys()->Contains(config.c_str())) {
    cout << "Creating directory for " << config << " configuration" << endl;
    f->mkdir(config.c_str());
  }
  f->cd(config.c_str());
  fHistMap[name]->Write();
  f->Close();
  delete f;
  f=0;
}

void Zbeam::DeleteHistogram(std::string name, std::string config, std::string filename)
{
  TFile* f=new TFile(filename.c_str(),"UPDATE");
  f->cd(config.c_str());
  f->rmdir(name.c_str());
  f->Close();
  delete f;
  f=0;
}

void Zbeam::FillVector(TH1D* hist, int ntype, int eff, int beam, int det)
{
  EffMap_t map1;
  BeamMap_t map2;
  DetMap_t map3;

  const int NBINS=hist->GetNbinsX();
  double* vec=new double[NBINS];

  for (int i=1; i<hist->GetNbinsX()+1;i++) vec[i-1]=hist->GetBinContent(i);

  map1.insert(EffMap_t::value_type(eff, vec));
  map2.insert(BeamMap_t::value_type(beam, map1));
  map3.insert(DetMap_t::value_type(det, map2));

  std::map<int, DetMap_t>::iterator bsmit = fBeamSysMap.find(ntype);

  if (bsmit == fBeamSysMap.end()) {
    fBeamSysMap.insert(std::map<int, DetMap_t>::value_type(ntype, map3));
  } else {
    DetMap_t::iterator dmit = (bsmit->second).find(det);
    if (dmit == (bsmit->second).end() ) {
      ((*bsmit).second).insert(DetMap_t::value_type(det,map2));
    } else {
      BeamMap_t::iterator bmit = (dmit->second).find(beam);
      if (bmit == (dmit->second).end() ) {
        ((*dmit).second).insert(BeamMap_t::value_type(beam,map1));
      } else {
        EffMap_t::iterator effit = (bmit->second).find(eff);
        if (effit == (bmit->second).end() ) {
          ((*bmit).second).insert(EffMap_t::value_type(eff,vec));
        } else {
          cout << "Already loaded vector for: " 
               << "beam "<<beam<< " det "<<det<<" nu "<<ntype<<" eff "<<eff<<endl;
          return;
        }
      }
    }
  } 
}

---