NCPOTCounter.cxx

Go to the documentation of this file.

#include "NCPOTCounter.h"

#include "NCUtils/NCType.h"
using namespace NCType;

#include "MessageService/MsgService.h"

CVSID("$Id: NCPOTCounter.cxx,v 1.14 2009/02/13 18:21:37 bckhouse Exp $");

#include "TFile.h"
#include "TH2F.h"
#include "TList.h"
#include "TString.h"
#include "TSystemDirectory.h"
#include "TRegexp.h"

#include <cassert>

using std::map;
using std::vector;


NCPOTCounter::Cache NCPOTCounter::fgCache;

NCPOTCounter::NCPOTCounter(TString dataMCPath,
                           TString mockDataSet,
                           bool useMCAsData,
                           bool useMockData,
                           int fileLimitNearData,
                           int fileLimitFarData,
                           int fileLimitNearMC,
                           int fileLimitFarMC,
                           vector<BeamType::BeamType_t> beamIndex,
                           int runToUse,
                           int mockDataSubRun,
                           bool scaleNearExposures,
                           double farExposureScaleVal)
  :fDataMCPath(dataMCPath),
   fMockDataSet(mockDataSet),
   fBeamIndex(beamIndex),
   fRunToUse(runToUse),
   fMockDataSubRun(mockDataSubRun),
   fScaleNearExposures(scaleNearExposures),
   fFarExposureScaleVal(farExposureScaleVal)
{
  fFileLimits[Detector::kNear][NCType::kData] = fileLimitNearData;
  fFileLimits[Detector::kFar] [NCType::kData] = fileLimitFarData;
  fFileLimits[Detector::kNear][NCType::kMC]   = fileLimitNearMC;
  fFileLimits[Detector::kFar] [NCType::kMC]   = fileLimitFarMC;

  fMockOrMCAsData = useMCAsData || useMockData;
}


void NCPOTCounter::SetPOTValues(map<NCBeam::Info, double>& pot,
                                map<NCBeam::Info, double>& scale,
                                Detector::Detector_t det,
                                NCType::EFileType fileType,
                                NCType::EDataMC dataMC)
{
  POTKey key(3);
  key[0] = det;
  key[1] = fileType;
  key[2] = dataMC;

  Cache::iterator it = fgCache.find(key);

  if(it != fgCache.end()){
    pot = it->second.first;
    scale = it->second.second;
    return;
  }


  pot.clear();
  scale.clear();

  MSG("NCPOTCounter", Msg::kInfo)
    << "NCPOTCounter::SetPOTValues:"
    << " fDataMCPath = " << fDataMCPath << endl;

  //pot vector needs to have entry for each run type
  vector<vector<double> > POT;
  POT.resize(NCType::kMaxRun);
  for(int i = 0; i < NCType::kMaxRun; ++i){
    POT[i].resize(fBeamIndex.size());
    for(unsigned int j = 0; j < fBeamIndex.size(); ++j){
      POT[i][j] = 0.;
    }
  }

  //loop over the BeamType::AsString values and fill the
  //data and mc POT numbers
  for(UInt_t bt = 0; bt < fBeamIndex.size(); ++bt){

    const vector<TString> fileList = GetListOfFiles(bt, fileType, det, dataMC);
    for(unsigned int nf = 0; nf < fileList.size(); ++nf){
      const TString fileName = fileList[nf];

      const double potToAdd = GetPOTForSingleFile(fileName, fileType,
                                                  det, dataMC, bt);
      assert(potToAdd >= 0);

      const int runType = NCType::FindRunType(fileName);
      POT[runType-1].at(bt) += potToAdd;

      MSG("NCPOTCounter", Msg::kInfo) << "  Adding " << potToAdd << " POT to "
                                      << FileDesignator(fileType, det, dataMC)
                                      << endl;
    } // end for nf


    const TString beamType = BeamType::AsString(fBeamIndex[bt]);
    MSG("NCPOTCounter", Msg::kInfo) << "POT INFORMATION("
                                    << (det == Detector::kFar ? "far" : "near")
                                    << "): " << endl
                                    << FileDesignator(fileType, det, dataMC)
                                    << " run I POT[" << beamType << "] = "
                                    << POT[0].at(bt)
                                    << " run II POT[" << beamType << "] = "
                                    << POT[1].at(bt)
                                    << endl;

    for(int i = 0; i < NCType::kMaxRun; ++i){
      const NCBeam::Info beamInfo(fBeamIndex[bt], NCType::ERunType(i+1));
      if(POT[i].at(bt) > 0.){
        pot[beamInfo] = POT[i][bt];
        scale[beamInfo] = 1.;
      }//end if pot > 0 for this combination of run type and beam
    }//end loop over run types
  }//end loop to fill pot map

  DoScaling(pot, scale, det, dataMC);

  fgCache[key] = POTVal(pot, scale);
}


double NCPOTCounter::GetPOTValue(Detector::Detector_t det,
                                 NCType::EFileType fileType,
                                 NCType::EDataMC dataMC,
                                 BeamType::BeamType_t beamType,
                                 NCType::ERunType run)
{
  map<NCBeam::Info, double> pot;
  map<NCBeam::Info, double> scale;
  SetPOTValues(pot, scale, det, fileType, dataMC);

  const NCBeam::Info beamInfo(beamType, run);

  return pot[beamInfo];
}


double NCPOTCounter::GetScaleValue(Detector::Detector_t det,
                                   NCType::EFileType fileType,
                                   NCType::EDataMC dataMC,
                                   BeamType::BeamType_t beamType,
                                   NCType::ERunType run)
{
  POTMap pot;
  ScaleMap scale;
  SetPOTValues(pot, scale, det, fileType, dataMC);

  const NCBeam::Info beamInfo(beamType, run);

  return scale[beamInfo];
}


void NCPOTCounter::DoScaling(POTMap& pot,
                             ScaleMap& scale,
                             Detector::Detector_t det,
                             NCType::EDataMC dataMC) const
{
  if(fScaleNearExposures &&
     det == Detector::kNear &&
     dataMC == NCType::kData){
    // find smallest exposure
    double smallest = 1e25;
    for(POTMap::iterator it = pot.begin(); it != pot.end(); ++it){
      if(it->second < smallest)
        smallest = it->second;
    }

    ToFixedScale(pot, scale, smallest);
  }


  if(fFarExposureScaleVal != 1 &&
     det == Detector::kFar &&
     dataMC == NCType::kData){
    //scale the far detector data pot to whatever the scale factor is
    ToFixedScale(pot, scale, fFarExposureScaleVal);
  }
}


void NCPOTCounter::ToFixedScale(POTMap& pot, ScaleMap& scale,
                                double fixedScale) const
{
  for(POTMap::iterator itr = pot.begin(); itr != pot.end(); ++itr){
    scale[itr->first] = fixedScale/itr->second;
    itr->second = fixedScale;

    // unfortunately don't know what detector we are by now...
    MSG("NCPOTCounter", Msg::kInfo) << "set "//far "
                                    << itr->first.GetDescription()
                                    << " scale to "
                                    << scale[itr->first]
                                    << " - " << itr->second << " POT"
                                    << endl;
  }
}


vector<TString> NCPOTCounter::GetListOfFiles(int bt,
                                             NCType::EFileType fileType,
                                             Detector::Detector_t det,
                                             NCType::EDataMC dataMC) const
{
//   MSG("NCPOTCounter", Msg::kInfo) << "GetListOfFiles(" << bt << ", " << fileType
//                                   << ", " << det << ", " << dataMC << endl;
  TDirectory* saveDir = gDirectory;

  vector<TString> ret;

  const TString fileDesignator = FileDesignator(fileType, det, dataMC);

  const TString beamType = BeamType::AsString(fBeamIndex[bt]);
  const TString detName = (det == Detector::kFar) ? "far" : "near";

  //have a path to the relevant uDST's.  each uDST has the
  //BeamType::AsString() for the beam in its name, also has near
  //or far in the name.
  //make a TSystemDirectory object to get a list of files in the
  //path for the data and mc

  //the files for different beam configurations are in uDST/LXXXzYYYi
  //directories.  that is directories named for the beam type
  TSystemDirectory dataDir(beamType, fDataMCPath+beamType+"/");
  TList* files = dataDir.GetListOfFiles();
  assert(files);

  for(int de = 0; de < files->GetEntries(); ++de){
    TString fileName = files->At(de)->GetName();

    //Do nothing for "." and ".." entries
    if(fileName == "." || fileName == "..") continue;

    //check that this file is the right kind
    if(fileName.Contains(beamType) &&
       fileName.Contains("strip") &&
       fileName.Contains(detName) &&
       fileName.Contains(fileDesignator)){

      //run I = 1, runII = 2, etc
      const int runType = NCType::FindRunType(fileName);
      if(dataMC == NCType::kData
         && fRunToUse != NCType::kRunAll
         && runType != fRunToUse) continue;

      if(fileType == NCType::kMockFile && fMockDataSubRun != -1){
        // For mock data, only add files with subrun requested
        int subRun = MockFileSubRun(fileName);
        // If it's not the one we want, don't add it to the chain.
        //
        // subRun==-1 if the file is the one-uDST-per-run type. In
        // that case, the fileDesignator check above is enough to tell
        // us that we want to use this file, so don't skip it
        if(subRun != -1 && subRun != fMockDataSubRun) continue;
      }
      ret.push_back(fileName);
    }//end if the right detector and beam and designator
  }//end loop over data files

  delete files;

  saveDir->cd();

  return ret;
}



double NCPOTCounter::GetPOTForSingleFile(TString fileName,
                                         NCType::EFileType fileType,
                                         Detector::Detector_t det,
                                         NCType::EDataMC dataMC,
                                         int bt) const
{
  //  TDirectory* saveDir = gDirectory;

  TFile tf(fileName);

  const TString histName = HistName(fileType, dataMC);
  const TString fileDesignator = FileDesignator(fileType, det, dataMC);

  TH1F* potHist = (TH1F*)(tf.Get(histName));

  if(fileType==NCType::kMockFile) {
    // The one-subrun-per-uDST files have the correct POT stored in
    // the dataPOT histogram. The one-run-per-uDST files have to use
    // the hardcoded value from POTPerFile
    if (MockFileSubRun(fileName) != -1) return potHist->Integral()*1e12;
    else return POTPerFile(bt, det, fileType);
  }

  MSG("NCPOTCounter", Msg::kInfo) << "  Running on "
                                  << fileName << " "
                                  << potHist->GetEntries()
                                  << " files" << endl;


  //use all pot's for data
  if(fileDesignator == "data" || fileDesignator == "blind"){
    double potToAdd = potHist->Integral()*1e12;

    //2008/03/06 BJR: lots of network outages at soudan in june 2005
    //means that spill triggers were not produced during those times
    //and the pot counting done in the NCAnalysisModule is off by
    //2.56e18 POT.  Add it back in explicitly here.
    if(fileDesignator == "blind" && fileName.Contains("2005-06"))
      potToAdd += 2.56e18;

    return potToAdd;
  }
  
  const double potPerFile = POTPerFile(bt, det, fileType);

  const int fileLimit = fFileLimits[det][dataMC];
  if(fileLimit <= potHist->GetEntries()){
    return potPerFile*fileLimit;
  }
  else{
    return potHist->Integral()*1e12;
  }
}


TString NCPOTCounter::FileDesignator(NCType::EFileType fileType,
                                     Detector::Detector_t det,
                                     NCType::EDataMC dataMC) const
{
  if(fileType == NCType::kMockFile) return fMockDataSet;

  if(dataMC == NCType::kData && !fMockOrMCAsData){
    if(det == Detector::kNear) return "data";
    else if(det == Detector::kFar) return "blind";
  }

  if(fileType == NCType::kTauFile) return "tau";

  if(fileType == NCType::kElectronFile) return "electron";

  return "mc";
}


TString NCPOTCounter::HistName(NCType::EFileType fileType,
                               NCType::EDataMC dataMC) const
{
  if(fileType == NCType::kMockFile) return "dataPOT";

  if(dataMC == NCType::kData && !fMockOrMCAsData) return "dataPOT";

  return "mcPOT";
}


double NCPOTCounter::POTPerFile(int bt, Detector::Detector_t det,
                                NCType::EFileType fileType) const
{
  //really need to do something better than the next few lines
  //to account for changes in MC version

  if(det == Detector::kFar){
    if(fileType == NCType::kMockFile) {
      const double kPotFarMock = 2.5e20;
      return kPotFarMock;
    }

    // TODO: Is this a good name for it? If so should go in NCType
    const double kPotFarMC = 6.5e20;
    return kPotFarMC;
  }

  const double ret = kPOTNear[fBeamIndex[bt]];
  assert(ret >= 0);
  return ret;
}

int NCPOTCounter::MockFileSubRun(TString fileName) const
{
  // Filenames look like:
  // "far_mock_L010z185i_F21910001_0000.uDST_strip.root"
  //
  // ROOT is useless and doesn't support character classes, or
  // specifying a number of times to match, so we get this RE
  // to find the subrun from the filename.

  TRegexp re("_[0-9][0-9][0-9][0-9]\\.uDST");
  int stupidROOT;
  int index = re.Index(fileName, &stupidROOT);
  if(index != kNPOS) {
    // We found a subrun string in the filename, so use that.
    TString subRunStr(fileName(index+1, 4));
    return subRunStr.Atoi();
  } else {
    // This will happen for one-uDST-per-run files
    return -1;
  }
}

---