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.33 2009/10/30 17:52:04 rodriges Exp $");

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

#include <cassert>

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


NCPOTCounter::Cache NCPOTCounter::fgCache;

NCPOTCounter::NCPOTCounter()
{
}


void NCPOTCounter::Config(const Registry& r,
                          vector<BeamType::BeamType_t> beamIndex,
                          vector<NC::RunUtil::ERunType> runsToUse)
{
  const char* tmps;
  int tmpb, tmpi;

  // TODO - this code is copy-pasted from NCExtrapolationModule
  if(r.Get("DataMCPath", tmps)){
    //for some reason the TSystemDirectory constructor can't resolve
    //environmental variables, so just take care of it here
    TString temp = tmps;
    // Assume the form $ENVVAR/path
    if(temp.Contains("$")){
      int pos = temp.Index("/"); // find the first /
      temp.Resize(pos); // and truncate at it
      temp.Remove(TString::kLeading, '$'); // drop the $ off the front
      fDataMCPath = tmps; // Copy the original string
      // and then replace the part that we think is the environment variable
      fDataMCPath.Replace(0, pos, gSystem->Getenv(temp));
    }
    else fDataMCPath = temp;
    // Add a trailing slash so that concatenating the path later will work
    if(!fDataMCPath.EndsWith("/")) fDataMCPath += "/";
  }

  if(r.Get("MockDataSet", tmps)) fMockDataSet = tmps;
  if(r.Get("UseMCAsData", tmpb)) fMCAsData = tmpb;
  if(r.Get("UseMockData", tmpb)) fMockData = tmpb;

  if(r.Get("RunLimitNearData", tmpi))
    fRunLimits[Detector::kNear][NCType::kData] = tmpi;
  if(r.Get("RunLimitFarData", tmpi))
    fRunLimits[Detector::kFar][NCType::kData] = tmpi;
  if(r.Get("RunLimitNearMC", tmpi))
    fRunLimits[Detector::kNear][NCType::kMC] = tmpi;
  if(r.Get("RunLimitFarMC", tmpi))
    fRunLimits[Detector::kFar][NCType::kMC] = tmpi;

  if(r.Get("MockDataSubRun", tmpi)) fMockDataSubRun = tmpi;

  fBeamIndex = beamIndex;
  fRunsToUse = runsToUse;
}


void NCPOTCounter::SetPOTValues(map<NCBeam::Info, double>& pot,
                                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;
    return;
  }


  pot.clear();

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


  // Make POT map with an entry for each run/beam type
  // combination. If we actually have nonzero POTs for a particular
  // combination, it'll get added to the 'pot' map at the end.
  map<NCBeam::Info, double> POTtmp;

  for(int r = NC::RunUtil::kRunI; 
      r <= NC::RunUtil::kMaxRun; ++r){
    for(unsigned int j = 0; j < fBeamIndex.size(); ++j){
      NCBeam::Info info(fBeamIndex[j], NC::RunUtil::ERunType(r));
      POTtmp[info] = 0.;
    }
  }

  //loop over the beams/runperiods we're using and fill the pot values
  for(unsigned int bt = 0; bt < fBeamIndex.size(); ++bt){
    for(unsigned int runIdx = 0; runIdx<fRunsToUse.size(); ++runIdx){
      NCBeam::Info info(fBeamIndex[bt], fRunsToUse[runIdx]);

      // The total number of runs we've read from files,
      // so we can apply the file limits
      int nRunsRead = 0;
      
      const vector<TString> fileList = GetListOfFiles(info, 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, nRunsRead);
        assert(potToAdd >= 0);
        
        //const int runType = NC::RunUtil::FindRunType(fileName);

        assert(POTtmp.find(info)!=POTtmp.end());
        POTtmp[info] += potToAdd;
        
        MSG("NCPOTCounter", Msg::kInfo) << "  Adding " << potToAdd << " POT to "
                                        << FileDesignator(fileType, det, dataMC)
                                        << endl;

        nRunsRead += GetRunsInFile(fileName);
        // If we've already hit the file limit, stop
        if(nRunsRead > fRunLimits[det][dataMC]) break;

      } // 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 << "] = "
        << POTtmp[0].at(bt)
        << " run II POT[" << beamType << "] = "
        << POTtmp[1].at(bt)
        << endl;
      */

      if(POTtmp[info] > 0){
        if(dataMC == NCType::kData && fMCAsData){
          assert(fFarMCExposures.find(info) != fFarMCExposures.end() &&
                 "Did you call set_mc_exposure_for_beam for this beam?");
          pot[info] = fFarMCExposures[info];
        }
        else{
          pot[info] = POTtmp[info];
        }
      }//end if pot > 0 for this combination of run type and beam      
     
    }//end loop over run types 

    // Now we've looped over all the run periods, add their total to
    // the RunAll NCBeam for this beamtype

    double totalPOT=0;

    for(unsigned int n = 0; n < fRunsToUse.size(); ++n){
      NCBeam::Info info(fBeamIndex[bt], fRunsToUse[n]);
      totalPOT+=pot[info];
    }
      
    const NCBeam::Info infoAll(fBeamIndex[bt], NC::RunUtil::kRunAll);
    pot[infoAll] = totalPOT;
      
  }//end loop over beam types
    
  for(map<NCBeam::Info, double>::const_iterator it = pot.begin();
      it != pot.end(); ++it){
    MSG("NCPOTCounter", Msg::kInfo) << "POT INFO: "
                                    << (det == Detector::kFar ? "Far" : "Near")
                                    << " "
                                    << (dataMC == NCType::kMC ? "MC" : "data")
                                    << " from \""
                                    << FileDesignator(fileType, det, dataMC)
                                    << "\" files " << it->first <<  " has "
                                    << it->second << " POT" << endl;
  }
    
  fgCache[key] = pot;
}


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

  SetPOTValues(pot, det, fileType, dataMC);

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

  return pot[beamInfo];
}


vector<TString> NCPOTCounter::GetListOfFiles(NCBeam::Info info,
                                             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(info.GetBeamType());
  const TString runPeriod= Form("run%d", (int)info.GetRunType());
  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)){

      // Check it's the right run type - this should work for both data and mc
      const NC::RunUtil::ERunType fileRunPeriod = NC::RunUtil::FindRunType(fileName);
      if(fileRunPeriod != info.GetRunType()) 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;
}

int NCPOTCounter::GetRunsInFile(TString fileName)
{
  TFile tf(fileName);
  assert(!tf.IsZombie());
  TH1F* potHist = (TH1F*)(tf.Get("mcPOT"));
  if(!potHist) potHist=(TH1F*)(tf.Get("dataPOT"));
  assert(potHist);
  // Hopefully this cast does what I want. Why does GetEntries() not return an int? Gah
  return int(potHist->GetEntries());
}


void NCPOTCounter::SetMCExposureForBeam(NCBeam::Info beam, double pot)
{
  fFarMCExposures[beam] = pot;
}


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

  TFile tf(fileName);
  assert(!tf.IsZombie());
  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 POTPerRun
    if (MockFileSubRun(fileName) != -1)
      return potHist->Integral(0, potHist->GetNbinsX()+1)*1e12;
    else return POTPerRun(bt, det, fileType);
  }

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


  //use all pot's for data
  if(fileDesignator == "data" || fileDesignator == "blind"){
    return potHist->Integral(0, potHist->GetNbinsX()+1)*1e12;
  }

  const double potPerRun = POTPerRun(bt, det, fileType);

  const int runLimit = fRunLimits[det][dataMC]-nRunsOffset;
  if(runLimit <= potHist->GetEntries()){
    return potPerRun*runLimit;
  }
  else{
    return potHist->Integral(0, potHist->GetNbinsX()+1)*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 && !(fMCAsData || fMockData)){
    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 ||
     (dataMC == NCType::kData && !(fMCAsData || fMockData))) return "dataPOT";

  return "mcPOT";
}


double NCPOTCounter::POTPerRun(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;
  }
}

---