NCEventAdder.cxx

Go to the documentation of this file.

//$Id: NCEventAdder.cxx,v 1.3 2010/04/22 09:03:24 rodriges Exp $
//
//NCEventAdder.cxx
//
//C. Backhouse 7/2009


#include "NCEventAdder.h"

#include "NCExtrapolationModule.h"

#include "../NCUtility.h"

#include "../NCEventInfo.h"

#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "Conventions/SimFlag.h"
#include "MessageService/MsgService.h"

#include "TRandom3.h"
#include "TStopwatch.h"

#include <cassert>
#include <climits>

CVSID("");


//-----------------------------------------------------------------------
NC::EventAdderBase::EventAdderBase()
  : fRunLimitNearMC(-1),
    fRunLimitNearData(-1),
    fRunLimitFarMC(-1),
    fRunLimitFarData(-1)
{
}


//-----------------------------------------------------------------------
NC::EventAdderBase::~EventAdderBase()
{
}


//-----------------------------------------------------------------------
void NC::EventAdderBase::Config(const Registry& r)
{
  bool ok = true;

  ok = ok && r.Get("RunLimitNearMC",   fRunLimitNearMC);
  ok = ok && r.Get("RunLimitNearData", fRunLimitNearData);
  ok = ok && r.Get("RunLimitFarMC",    fRunLimitFarMC);
  ok = ok && r.Get("RunLimitFarData",  fRunLimitFarData);

  assert(ok);
}


//-----------------------------------------------------------------------
void NC::EventAdderBase::ReportOnEventInfo(const NCEventInfo* eventInfo) const
{
  using namespace Detector;
  using namespace SimFlag;
  using namespace BeamType;

  const Detector_t det = Detector_t(eventInfo->header->detector);
  const SimFlag_t sim = SimFlag_t(eventInfo->header->dataType);
  const BeamType_t beam = BeamType_t(eventInfo->beam->beamType);

  MSG("NC::EventAdderBase", Msg::kInfo)
    << "This is "
    << eventInfo->header->softVersion << " "
    << AsString(det) << " "
    << AsString(beam) << " "
    << AsString(sim)
    << endl;


  // Hackish check to refuse to run on old-style files
  if(beam == 2){
    MSG("NC::EventAdderBase", Msg::kFatal)
      << "This file has beam.beamType==2, which is L000z200i. "
      << "You are probably running on an old file with the Zarko "
      << "numbering scheme (where 2 is L010z185i), but if you "
      << "really do want L000z200i, remove this message and recompile."
      << endl;
    abort();
  }
}


//-----------------------------------------------------------------------
const Registry& NC::EventAdderDefaultConfig()
{
  static Registry r;

  r.UnLockValues();

  r.Set("RunLimitNearData", INT_MAX);
  r.Set("RunLimitFarData",  INT_MAX);
  r.Set("RunLimitNearMC",   INT_MAX);
  r.Set("RunLimitFarMC",    INT_MAX);

  r.Set("UseMCAsData",      false);
  r.Set("SplitMC",          false);
  r.Set("SplitMCSeed",      -1);
  r.Set("UseMockData",      false);
  r.Set("MockDataSet",      "F21910001");
  r.Set("MockDataSubRun",    -1);

  r.LockValues();
  return r;
}


//-----------------------------------------------------------------------
NC::IEventAdder* NC::EventAdderBuilder(const Registry& r)
{
  NC::IEventAdder* adder = 0;

  // These should be bools, but the registry is stupid
  int useMCAsData, splitMC, useMockData;

  bool ok = true;
  ok = ok && r.Get("UseMCAsData", useMCAsData);
  ok = ok && r.Get("SplitMC",     splitMC);
  ok = ok && r.Get("UseMockData", useMockData);
  assert(ok);

  if(useMockData){
    MSG("NCEventAdder", Msg::kInfo)
      << "Creating NC::MockDataAdder" << endl;
    adder = new NC::MockDataAdder;
  }
  else{
    if(useMCAsData){
      if(splitMC){
        MSG("NCEventAdder", Msg::kInfo)
          << "Creating NC::SplitFakeDataAdder" << endl;
        adder = new NC::SplitFakeDataAdder;
      }
      else{
        MSG("NCEventAdder", Msg::kInfo)
          << "Creating NC::FakeDataAdder" << endl;
        adder = new NC::FakeDataAdder;
      } // end if splitMC
    }
    else{
      MSG("NCEventAdder", Msg::kInfo) << "Creating NC::RealDataAdder" << endl;
      adder = new NC::RealDataAdder;
    } // end if useMCAsData
  } // end if useMockData

  assert(adder);
  adder->Config(r);
  return adder;
}


//-----------------------------------------------------------------------
class LoopRunLimitHelper
{
public:
  LoopRunLimitHelper(int runLimit, int eventTot, const NCEventInfo* eventInfo)
    : fEventTot(eventTot), fRunLimit(runLimit), fRunCounter(0)
  {
    fPrevSubRun = eventInfo->header->subRun;
    fPrevRun = eventInfo->header->run;

    if(fRunLimit != INT_MAX){
      MSG("RunLimitHelper", Msg::kInfo)
        << fEventTot << " total events "
        << "limited to " << fRunLimit  << " (sub)runs" << endl;
    }
    else{
      MSG("RunLimitHelper", Msg::kInfo)
        << fEventTot << " total events, no (sub)run limit" << endl;
    }
  }

  ~LoopRunLimitHelper()
  {
     NC::Utility::ReportProgress(1, fStopwatch);
     MSG("RunLimitHelper", Msg::kInfo) << endl;
  }

  bool AreDone(int i, const NCEventInfo* eventInfo)
  {
    const int thisSubRun = eventInfo->header->subRun;
    const int thisRun = eventInfo->header->run;

    if(thisRun != fPrevRun || thisSubRun != fPrevSubRun){
      fPrevRun = thisRun;
      fPrevSubRun = thisSubRun;
      ++fRunCounter;
    }

    // Don't do it on the very first few, other people want to print
    if(i % 1000 == 10)
      NC::Utility::ReportProgress(std::max(float(i)/fEventTot,
                                           float(fRunCounter)/fRunLimit),
                                  fStopwatch);

    return (fRunCounter > fRunLimit || i >= fEventTot);
  }

  int GetRunCounter(){return fRunCounter;}

protected:
  TStopwatch fStopwatch; 

  int fEventTot;   
  int fRunLimit;   
  int fPrevSubRun; 
  int fPrevRun;    
  int fRunCounter; 
};




//-----------------------------------------------------------------------
NC::FakeDataAdder::FakeDataAdder() : fSkipMC(false), fSkipData(false)
{
}


//-----------------------------------------------------------------------
NC::FakeDataAdder::~FakeDataAdder()
{
}


//-----------------------------------------------------------------------
void NC::FakeDataAdder::Config(const Registry& r)
{
  NC::EventAdderBase::Config(r);

  bool ok = true;

  const char* tmps;
  ok = ok && r.Get("LoadExtrapolationsFile", tmps);
  fSkipMC = TString(tmps) != "";
  int tmpi;
  ok = ok && r.Get("IncludeDataFromLoadedExtrapolations", tmpi);
  // This key only takes effect when we're loading anyway
  fSkipData = fSkipMC && tmpi;

  assert(ok);
}


//-----------------------------------------------------------------------
void NC::FakeDataAdder::AddEvents(NCExtrapolationModule* mod,
                                  NCEventInfo* eventInfo,
                                  TChain* /*nearData*/,
                                  TChain* nearMC,
                                  TChain* /*farData*/,
                                  TChain* farMC,
                                  TChain* farMCTau,
                                  TChain* farMCElectron) const
{
  if(fSkipMC && fSkipData) return;

  // Add near MC as both MC and fake data
  AddEventsToExtrapolations(mod, eventInfo, nearMC);

  // Add all the far MC and as data as well
  AddEventsToExtrapolations(mod, eventInfo, farMCTau);
  AddEventsToExtrapolations(mod, eventInfo, farMCElectron);
  AddEventsToExtrapolations(mod, eventInfo, farMC);
}


//-----------------------------------------------------------------------
void NC::FakeDataAdder::AddEventsToExtrapolations(NCExtrapolationModule* mod,
                                                  NCEventInfo* eventInfo,
                                                  TChain* chain) const
{
  if(chain->GetEntries() < 1){
    MSG("NC::FakeDataAdder", Msg::kError) << "No entries in chain!" << endl;
    return;
  }

  eventInfo->FillFromChain(chain, 0);

  ReportOnEventInfo(eventInfo);

  // The file limit
  Detector::Detector_t det = Detector::Detector_t(eventInfo->header->detector);
  const int limit = (det == Detector::kNear) ? fRunLimitNearMC : fRunLimitFarMC;

  const int eventTot = chain->GetEntries();

  LoopRunLimitHelper runLimitHelper(limit, eventTot, eventInfo);

  //loop over the entries in the tree
  for(int i = 0; i < eventTot; ++i){
    eventInfo->FillFromChain(chain, i);

    if(runLimitHelper.AreDone(i, eventInfo)) break;

    if(!fSkipMC) mod->AddEventToExtrapolations(false);
    //    chain->GetEntry(i);
    if(!fSkipData) mod->AddEventToExtrapolations(true);
  }//end loop over data chain
}




//-----------------------------------------------------------------------
NC::SplitFakeDataAdder::SplitFakeDataAdder() : fSplitMCSeed(-1)
{
}


//-----------------------------------------------------------------------
NC::SplitFakeDataAdder::~SplitFakeDataAdder()
{
}


//-----------------------------------------------------------------------
void NC::SplitFakeDataAdder::Config(const Registry& r)
{
  NC::EventAdderBase::Config(r);

  const bool ok = r.Get("SplitMCSeed", fSplitMCSeed);
  assert(ok);
}


//-----------------------------------------------------------------------
void NC::SplitFakeDataAdder::AddEvents(NCExtrapolationModule* mod,
                                       NCEventInfo* eventInfo,
                                       TChain* /*nearData*/,
                                       TChain* nearMC,
                                       TChain* /*farData*/,
                                       TChain* farMC,
                                       TChain* farMCTau,
                                       TChain* farMCElectron) const
{
  // Add near MC as both MC and fake data
  AddEventsToExtrapolations(mod, eventInfo, nearMC);

  // Add all the far MC and as data as well
  AddEventsToExtrapolations(mod, eventInfo, farMCTau);
  AddEventsToExtrapolations(mod, eventInfo, farMCElectron);
  AddEventsToExtrapolations(mod, eventInfo, farMC);
}


//-----------------------------------------------------------------------
void NC::SplitFakeDataAdder::
AddEventsToExtrapolations(NCExtrapolationModule* mod,
                          NCEventInfo* eventInfo,
                          TChain* chain) const
{
  if(chain->GetEntries() < 1){
    MSG("NC::SplitFakeDataAdder", Msg::kError) << "No entries in chain!" << endl;
    return;
  }

  eventInfo->FillFromChain(chain, 0);

  Detector::Detector_t det = Detector::Detector_t(eventInfo->header->detector);

  ReportOnEventInfo(eventInfo);

  // The file limit
  int limit = (det == Detector::kNear) ? fRunLimitNearMC : fRunLimitFarMC;
  // Need an even number of runs so we can split them equally
  if(TMath::Odd(limit)) --limit;

  const int eventTot = chain->GetEntries();

  LoopRunLimitHelper runLimitHelper(limit, eventTot, eventInfo);

  TRandom3 randomGen(fSplitMCSeed);

  //loop over the entries in the tree
  for(int i = 0; i < eventTot; ++i){
    eventInfo->FillFromChain(chain, i);

    if(runLimitHelper.AreDone(i, eventInfo)) break;

    bool asFakeData;
    if(fSplitMCSeed >= 0){
      // Put each event randomly into one spectrum or another
      asFakeData = randomGen.Uniform() > .5;
    }
    else{
      // Put alternate runs into data or MC
      asFakeData = (runLimitHelper.GetRunCounter()%2 == 0);
    }

    mod->AddEventToExtrapolations(asFakeData);

  }//end loop over data chain
}




//-----------------------------------------------------------------------
NC::RealDataAdder::RealDataAdder() : fSkipMC(false)
{
}


//-----------------------------------------------------------------------
NC::RealDataAdder::~RealDataAdder()
{
}


//-----------------------------------------------------------------------
void NC::RealDataAdder::Config(const Registry& r)
{
  NC::EventAdderBase::Config(r);

  const char* tmps;
  const bool ok = r.Get("LoadExtrapolationsFile", tmps);
  fSkipMC = TString(tmps) != "";

  assert(ok);
}


//-----------------------------------------------------------------------
void NC::RealDataAdder::AddEvents(NCExtrapolationModule* mod,
                                  NCEventInfo* eventInfo,
                                  TChain* nearData,
                                  TChain* nearMC,
                                  TChain* farData,
                                  TChain* farMC,
                                  TChain* farMCTau,
                                  TChain* farMCElectron) const
{
  if(!fSkipMC){
    AddEventsToExtrapolations(mod, eventInfo, nearMC);
    // Yes, really skip ND data if fSkipMC is set. fSkipMC is set if
    // we're reading the extrapolation from a file. The ND data is the
    // same every time, and so we can skip readding it
    AddEventsToExtrapolations(mod, eventInfo, nearData);
  }

  if(!fSkipMC){
    AddEventsToExtrapolations(mod, eventInfo, farMCTau);
    AddEventsToExtrapolations(mod, eventInfo, farMCElectron);
    AddEventsToExtrapolations(mod, eventInfo, farMC);
  }
  AddEventsToExtrapolations(mod, eventInfo, farData);
}


//-----------------------------------------------------------------------
void NC::RealDataAdder::
AddEventsToExtrapolations(NCExtrapolationModule* mod,
                          NCEventInfo* eventInfo,
                          TChain* chain) const
{
  if(chain->GetEntries() < 1){
    MSG("NC::RealDataAdder", Msg::kError) << "No entries in chain!" << endl;
    return;
  }

  eventInfo->FillFromChain(chain, 0);

  using namespace Detector;
  using namespace SimFlag;

  const Detector_t det = Detector_t(eventInfo->header->detector);
  const SimFlag_t sim = SimFlag_t(eventInfo->header->dataType);

  ReportOnEventInfo(eventInfo);

  // The file limit
  int limit = -1;
  if(det == kNear && sim == kMC)   limit = fRunLimitNearMC;
  if(det == kNear && sim == kData) limit = fRunLimitNearData;
  if(det == kFar  && sim == kMC)   limit = fRunLimitFarMC;
  if(det == kFar  && sim == kData) limit = fRunLimitFarData;
  assert(limit >= 0);

  const int eventTot = chain->GetEntries();

  LoopRunLimitHelper runLimitHelper(limit, eventTot, eventInfo);

  //loop over the entries in the tree
  for(int i = 0; i < eventTot; ++i){
    eventInfo->FillFromChain(chain, i);

    if(runLimitHelper.AreDone(i, eventInfo)) break;

    mod->AddEventToExtrapolations(false);

  }//end loop over data chain
}




//-----------------------------------------------------------------------
NC::MockDataAdder::MockDataAdder() : fMockDataSubRun(-1)
{
}


//-----------------------------------------------------------------------
NC::MockDataAdder::~MockDataAdder()
{
}


//-----------------------------------------------------------------------
void NC::MockDataAdder::Config(const Registry& r)
{
  NC::EventAdderBase::Config(r);

  bool ok = true;

  const char* tmps;
  ok = ok && r.Get("MockDataSet", tmps);
  fMockDataSet = tmps;

  ok = ok && r.Get("MockDataSubRun", fMockDataSubRun);

  assert(ok);
}


//-----------------------------------------------------------------------
void NC::MockDataAdder::AddEvents(NCExtrapolationModule* mod,
                                  NCEventInfo* eventInfo,
                                  TChain* /*nearData*/,
                                  TChain* nearMC,
                                  TChain* farData,
                                  TChain* farMC,
                                  TChain* farMCTau,
                                  TChain* farMCElectron) const
{
  AddEventsToExtrapolations(mod, eventInfo, nearMC);

  AddEventsToExtrapolations(mod, eventInfo, farMCTau);
  AddEventsToExtrapolations(mod, eventInfo, farMCElectron);
  AddEventsToExtrapolations(mod, eventInfo, farMC);
  AddEventsToExtrapolations(mod, eventInfo, farData);
}


//-----------------------------------------------------------------------
void NC::MockDataAdder::
AddEventsToExtrapolations(NCExtrapolationModule* mod,
                          NCEventInfo* eventInfo,
                          TChain* chain) const
{
  if(chain->GetEntries() < 1){
    MSG("NC::MockDataAdder", Msg::kError) << "No entries in chain!" << endl;
    return;
  }

  eventInfo->FillFromChain(chain, 0);

  using namespace Detector;
  using namespace SimFlag;

  const Detector_t det = Detector_t(eventInfo->header->detector);
  const SimFlag_t sim = SimFlag_t(eventInfo->header->dataType);
  const NCType::EFileType fileType = NCType::FindFileType(eventInfo->header);

  ReportOnEventInfo(eventInfo);

  // The file limit
  int limit = INT_MAX;
  if(det == kNear && sim == kMC) limit = fRunLimitNearMC;
  if(det == kFar  && sim == kMC) limit = fRunLimitFarMC;
  assert(limit >= 0);

  const int eventTot = chain->GetEntries();

  LoopRunLimitHelper runLimitHelper(limit, eventTot, eventInfo);

  //loop over the entries in the tree
  for(int i = 0; i < eventTot; ++i){
    eventInfo->FillFromChain(chain, i);

    const int thisSubRun = eventInfo->header->subRun;

    // For the case of just one mock data subrun, only use the subrun requested
    if(fileType == NCType::kMockFile &&
       fMockDataSubRun != -1 &&
       thisSubRun != fMockDataSubRun) continue;

    if(runLimitHelper.AreDone(i, eventInfo)) break;

    mod->AddEventToExtrapolations(false);
    if(det == kNear) mod->AddEventToExtrapolations(true);

  }//end loop over data chain
}

---