lem::MakePID Class Reference

Calculate final PID variables from match details. More...

Inheritance diagram for lem::MakePID:

Public Types
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Member Functions
	MakePID (const fhicl::ParameterSet &pset)
	~MakePID ()
virtual void	beginJob ()
virtual void	beginRun (art::Run &run)
virtual void	produce (art::Event &evt)
virtual void	reconfigure (const fhicl::ParameterSet &pset)
std::string	workerType () const
void	doBeginJob ()
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Protected Types
enum	WeightScheme {   kUnweighted, kSimple, kExp, kSigOnly,   kSigOnlyExp, kBkgOnly, kBkgOnlyExp, kSig50Only,   kSig90Only, kPi0Only }
enum	MatchVar {   kSigFrac, kBkgFrac, kNumuFrac, kY,   kQFrac, kSigFrac50, kSigFrac90, kL0,   kL1, kTrueEVis, kQESigFrac, kResSigFrac,   kDISSigFrac, kQENCFrac, kResNCFrac, kDISNCFrac,   kQENumuFrac, kResNumuFrac, kDISNumuFrac }

Protected Member Functions
void	AvgAndFit (const std::vector< MatchSummary > &matches, WeightScheme weightScheme, MatchVar matchVar, double &mean, double &fit) const
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Protected Attributes
bool	fScaleCalE
std::vector< double >	fExpTable
double	fNCScale
double	fSigScale
double	fEnrichScale
double	fExpLength
double	fExpPow
bool	fTrimMatchLists
bool	fUseFitInput
bool	fPositionalWeight
std::string	fMatchesLabel
std::string	fSummarizeLabel
TMVA::Reader *	fReader
dec::Forest *	fForest
TH1 *	fFlatHist

Detailed Description

Calculate final PID variables from match details.

Definition at line 48 of file MakePID_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 25 of file Producer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>

inherited

Definition at line 18 of file EDProducer.h.

Member Enumeration Documentation

enum lem::MakePID::MatchVar

protected

Enumerator
kSigFrac
kBkgFrac
kNumuFrac
kY
kQFrac
kSigFrac50
kSigFrac90
kL0
kL1
kTrueEVis
kQESigFrac
kResSigFrac
kDISSigFrac
kQENCFrac
kResNCFrac
kDISNCFrac
kQENumuFrac
kResNumuFrac
kDISNumuFrac

Definition at line 75 of file MakePID_module.cc.

                 {
      kSigFrac,
      kBkgFrac,
      kNumuFrac,
      kY,
      kQFrac,
      kSigFrac50,
      kSigFrac90,
      kL0,
      kL1,
      kTrueEVis,
      kQESigFrac,
      kResSigFrac,
      kDISSigFrac,
      kQENCFrac,
      kResNCFrac,
      kDISNCFrac,
      kQENumuFrac,
      kResNumuFrac,
      kDISNumuFrac
    };

enum lem::MakePID::WeightScheme

protected

Enumerator
kUnweighted
kSimple
kExp
kSigOnly
kSigOnlyExp
kBkgOnly
kBkgOnlyExp
kSig50Only
kSig90Only
kPi0Only

Definition at line 62 of file MakePID_module.cc.

                     {
      kUnweighted,
      kSimple,
      kExp,
      kSigOnly,
      kSigOnlyExp,
      kBkgOnly,
      kBkgOnlyExp,
      kSig50Only,
      kSig90Only,
      kPi0Only
    };

Constructor & Destructor Documentation

lem::MakePID::MakePID ( const fhicl::ParameterSet &  pset )

explicit

Definition at line 124 of file MakePID_module.cc.

References ana::assert(), om::cout, allTimeWatchdog::endl, util::EnvExpansion(), fFlatHist, fForest, fReader, lem::dec::Forest::FromFile(), fhicl::ParameterSet::get(), gJunk, gTMVAVar, path, reconfigure(), and string.

  : EDProducer(pset)
  {
    reconfigure(pset);

    produces<std::vector<lem::PIDDetails>>();
    produces<art::Assns<lem::PIDDetails, rb::Cluster>>();

    const std::string path     = pset.get<std::string>("LibraryPath");
    const std::string libPath  = util::EnvExpansion(path) + "/";

    fReader = 0;
    if(pset.get<bool>("FillTMVAs")){
      fReader = new TMVA::Reader;

      fReader->AddVariable("pidExtrap", &gTMVAVar[0]); // vars.fPIDFit
      fReader->AddVariable("yexp",      &gTMVAVar[1]); // vars.fMeanYExp
      fReader->AddVariable("qfracexp",  &gTMVAVar[2]); // vars.fMeanQFracExp
      fReader->AddVariable("ediff",     &gTMVAVar[3]); // vars.fEnergyDiffExp
      fReader->AddVariable("calE",      &gTMVAVar[4]); // 2*totalRecoGeV
      fReader->AddVariable("enrichExp", &gTMVAVar[5]); // vars.fEnrichFracExp

      fReader->AddSpectator("isSig", &gJunk);
      fReader->AddSpectator("myW",   &gJunk);

      const std::string tmvaPath = libPath+"weights/";

      fReader->BookMVA("BDTG", tmvaPath+"TMVAClassification_BDTG.weights.xml");
      fReader->BookMVA("BDT",  tmvaPath+"TMVAClassification_BDT.weights.xml");
      fReader->BookMVA("BDTD", tmvaPath+"TMVAClassification_BDTD.weights.xml");
      fReader->BookMVA("BDTMitFisher", tmvaPath+"TMVAClassification_BDTMitFisher.weights.xml");
      fReader->BookMVA("MLP",  tmvaPath+"TMVAClassification_MLP.weights.xml");

      fReader->BookMVA("BDTG_CC", tmvaPath+"TMVAClassification_BDTG_CC.weights.xml");
      fReader->BookMVA("MLP_CC",  tmvaPath+"TMVAClassification_MLP_CC.weights.xml");
    }

    // fANN = fANNEnrich = 0;
    // if(pset.get<bool>("FillANNs")){
    //   fANN = fann_create_from_file((libPath+"fann.net").c_str());
    //   assert(fANN);
    //   fann_print_connections(fANN);

    //   fANNEnrich = fann_create_from_file((libPath+"fann_enrich.net").c_str());
    //   assert(fANNEnrich);
    //   fann_print_connections(fANNEnrich);
    // }

    fForest = 0;
    if(pset.get<bool>("FillDecTree")){
      std::cout << "Loading LEM dectree from " << libPath << std::endl;
      fForest = dec::Forest::FromFile(libPath+"dectree.root");
    }

    fFlatHist = 0;
    if(pset.get<bool>("Flatten")){
      TFile* fFlat = new TFile((libPath+"flatten.root").c_str());
      assert(!fFlat->IsZombie());
      fFlatHist = (TH1*)fFlat->Get("conv");
      assert(fFlatHist);
    }
  }

lem::MakePID::~MakePID

(

)

Definition at line 188 of file MakePID_module.cc.

References fForest, and fReader.

  {
    delete fReader;
    delete fForest;
  }

Member Function Documentation

void lem::MakePID::AvgAndFit ( const std::vector< MatchSummary > &  matches, WeightScheme  weightScheme, MatchVar  matchVar, double &  mean, double &  fit  ) const

protected

Definition at line 236 of file MakePID_module.cc.

References abs(), plot_validation_datamc::c, lem::MatchSummary::ccnc, fExpLength, fExpPow, fExpTable, fPositionalWeight, MECModelEnuComparisons::i, makeTrainCVSamples::int, lem::MatchSummary::isSig, kBkgFrac, kBkgOnly, kBkgOnlyExp, simb::kCC, simb::kDIS, kDISNCFrac, kDISNumuFrac, kDISSigFrac, kExp, kL0, kL1, simb::kNC, kNumuFrac, kPi0Only, simb::kQE, kQENCFrac, kQENumuFrac, kQESigFrac, kQFrac, simb::kRes, kResNCFrac, kResNumuFrac, kResSigFrac, kSig50Only, kSig90Only, kSigFrac, kSigFrac50, kSigFrac90, kSigOnly, kSigOnlyExp, kTrueEVis, kUnweighted, kY, util::LinFit(), util::LinFitUnweighted(), m, cafExposure::match, lem::MatchSummary::mode, lem::MatchSummary::pdg, lem::MatchSummary::photE0, lem::MatchSummary::photE1, lem::MatchSummary::photL0, lem::MatchSummary::photL1, lem::MatchSummary::potential, cet::pow(), lem::MatchSummary::qFrac, PandAna.Demos.pi0_spectra::tot, lem::MatchSummary::trueEVis, w, lem::MatchSummary::weight, and lem::MatchSummary::y.

Referenced by produce().

  {
    std::vector<double> xs, ys, ws;
    double tot = 0;
    mean = 0;

    const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;

    for(unsigned int i = 0; i < matches.size(); ++i){
      const MatchSummary& match = matches[i];

      double w = match.weight;

      if(weightScheme == kUnweighted) w = 1;

      if(weightScheme == kExp ||
         weightScheme == kSigOnlyExp ||
         weightScheme == kBkgOnlyExp){
        if(fPositionalWeight){
          w *= fExpTable[i];
        }
        else{
          w *= exp(-pow(match.potential/maxPotential, fExpPow)/fExpLength);
        }
      }

      if(weightScheme == kSigOnly || weightScheme == kSigOnlyExp){
        if(!match.isSig) w = 0;
      }

      if(weightScheme == kBkgOnly || weightScheme == kBkgOnlyExp){
        if(match.isSig) w = 0;
      }

      if(weightScheme == kSig50Only){
        if(!match.isSig || match.y > 0.50) w = 0;
      }
      if(weightScheme == kSig90Only){
        if(!match.isSig || match.y > 0.90) w = 0;
      }
      if(weightScheme == kPi0Only) w = (match.photL0 > 0 || match.photL1 > 0 || match.photE0 > 0 || match.photE1 > 0);

      double var = 0;
      switch(matchVar){
      case kSigFrac: var = int(match.isSig); break;
      case kBkgFrac: var = 1-int(match.isSig); break;
      case kNumuFrac: var = int(match.ccnc == simb::kCC && abs(match.pdg) == 14); break;
      case kY: var = match.y; break;
      case kQFrac: var = match.qFrac; break;
      case kSigFrac50: var = int(match.isSig && match.y < .50); break;
      case kSigFrac90: var = int(match.isSig && match.y < .90); break;
      case kL0: var = match.photL0; break;
      case kL1: var = match.photL1; break;
      case kTrueEVis: var = match.trueEVis; break;
      case kQESigFrac:  var = (match.isSig && match.mode == simb::kQE);  break;
      case kResSigFrac: var = (match.isSig && match.mode == simb::kRes); break;
      case kDISSigFrac: var = (match.isSig && match.mode == simb::kDIS); break;
      case kQENCFrac:  var = (match.ccnc == simb::kNC && match.mode == simb::kQE);  break;
      case kResNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kRes); break;
      case kDISNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kDIS); break;
      case kQENumuFrac:  var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kQE);  break;
      case kResNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kRes); break;
      case kDISNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kDIS); break;
      }

      xs.push_back(match.potential);
      ys.push_back(var);
      ws.push_back(w);

      tot += w;
      mean += w*var;
    } // end for i

    mean /= tot;

    if(ws.size() < 2){
      fit = -1;
    }
    else{
      double m, c;
      if(weightScheme == kUnweighted)
        util::LinFitUnweighted(xs, ys, m, c);
      else
        util::LinFit(xs, ys, ws, m, c);
      fit = m*xs[0]+c;
    }
  }

void lem::MakePID::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 209 of file MakePID_module.cc.

References fExpLength, fExpPow, fExpTable, MECModelEnuComparisons::i, lem::kMaxNumMatches, and cet::pow().

  {
    fExpTable.resize(kMaxNumMatches);
    for(unsigned int i = 0; i < kMaxNumMatches; ++i){
      fExpTable[i] = exp(-pow(i/fExpLength, fExpPow));
    }
  }

void lem::MakePID::beginRun ( art::Run &  run )

virtual

Reimplemented from art::EDProducer.

Definition at line 218 of file MakePID_module.cc.

References ana::assert(), om::cout, art::Handle< T >::failedToGet(), fEnrichScale, fMatchesLabel, fNCScale, fSigScale, art::DataViewImpl::getByLabel(), lem::LibrarySummary::nSwapNC, lem::LibrarySummary::nTrueNC, lem::LibrarySummary::totBkg, lem::LibrarySummary::totEnrich, and lem::LibrarySummary::totSig.

  {
    art::Handle<lem::LibrarySummary> libsum;
    run.getByLabel(fMatchesLabel, libsum);

    assert(!libsum.failedToGet());

    fNCScale = libsum->nTrueNC/double(libsum->nTrueNC+libsum->nSwapNC);
    fSigScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totSig;
    fEnrichScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totEnrich;

    std::cout << "Reweighted to sig, bkg, enrich = "
              << libsum->totSig*fSigScale << ", "
              << libsum->totBkg+(fNCScale-1)*libsum->nTrueNC-libsum->nSwapNC << ", "
              << libsum->totEnrich*fEnrichScale;
  }

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag )

protectedinherited

Definition at line 55 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumes(), T, and getGoodRuns4SAM::tag.

  {
    return collector_.consumes<T, BT>(tag);
  }

ConsumesCollector& art::ModuleBase::consumesCollector ( )

protectedinherited

template<typename T , BranchType BT>

void art::ModuleBase::consumesMany ( )

protectedinherited

Definition at line 69 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumesMany(), and T.

  {
    collector_.consumesMany<T, BT>();
  }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag )

inherited

Definition at line 62 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumesView(), T, and getGoodRuns4SAM::tag.

  {
    return collector_.consumesView<T, BT>(tag);
  }

void art::detail::Producer::doBeginJob ( )

inherited

bool art::detail::Producer::doBeginRun ( RunPrincipal &  rp, ModuleContext const &  mc  )

inherited

bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal &  srp, ModuleContext const &  mc  )

inherited

void art::detail::Producer::doEndJob ( )

inherited

bool art::detail::Producer::doEndRun ( RunPrincipal &  rp, ModuleContext const &  mc  )

inherited

bool art::detail::Producer::doEndSubRun ( SubRunPrincipal &  srp, ModuleContext const &  mc  )

inherited

bool art::detail::Producer::doEvent ( EventPrincipal &  ep, ModuleContext const &  mc, std::atomic< std::size_t > &  counts_run, std::atomic< std::size_t > &  counts_passed, std::atomic< std::size_t > &  counts_failed  )

inherited

void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb )

inherited

void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb )

inherited

void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb )

inherited

void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb )

inherited

std::array<std::vector<ProductInfo>, NumBranchTypes> const& art::ModuleBase::getConsumables ( ) const

inherited

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag )

protectedinherited

Definition at line 76 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsume(), T, and getGoodRuns4SAM::tag.

  {
    return collector_.mayConsume<T, BT>(tag);
  }

template<typename T , BranchType BT>

void art::ModuleBase::mayConsumeMany ( )

protectedinherited

Definition at line 90 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsumeMany(), and T.

  {
    collector_.mayConsumeMany<T, BT>();
  }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag )

inherited

Definition at line 83 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsumeView(), T, and getGoodRuns4SAM::tag.

  {
    return collector_.mayConsumeView<T, BT>(tag);
  }

ModuleDescription const& art::ModuleBase::moduleDescription ( ) const

inherited

void lem::MakePID::produce ( art::Event &  evt )

virtual

Implements art::EDProducer.

Definition at line 327 of file MakePID_module.cc.

References ana::assert(), AvgAndFit(), plot_validation_datamc::c, lem::dec::Forest::Classify(), util::CreateAssn(), DEFINE_ART_MODULE(), e, lem::PIDExtraVars::fAvgExpE, lem::PIDExtraVars::fAvgInvE, lem::PIDExtraVars::fAvgSigE, lem::PIDExtraVars::fBestBkgPotential, lem::PIDExtraVars::fBestEnrichPotential, lem::PIDExtraVars::fBestSigPotential, lem::PIDExtraVars::fChiBkg, lem::PIDExtraVars::fChiBkgExp, lem::PIDExtraVars::fChiSig, lem::PIDExtraVars::fChiSigExp, lem::PIDExtraVars::fDecTreePID, lem::PIDExtraVars::fDISNCFracExp, lem::PIDExtraVars::fDISNumuFracExp, lem::PIDExtraVars::fDISSigFracExp, lem::PIDExtraVars::fEnergyDiff, lem::PIDExtraVars::fEnergyDiffExp, lem::PIDExtraVars::fEnrichFrac, lem::PIDExtraVars::fEnrichFracExp, lem::PIDExtraVars::fEnrichFracFit, lem::PIDExtraVars::fEnrichQFrac, lem::PIDExtraVars::fEnrichQFracExp, lem::PIDExtraVars::fEnrichQFracFit, fEnrichScale, fExpLength, fExpPow, fExpTable, fFlatHist, fForest, fMatchesLabel, lem::PIDExtraVars::fMeanL0Exp, lem::PIDExtraVars::fMeanL0Fit, lem::PIDExtraVars::fMeanL0PiOnlyExp, lem::PIDExtraVars::fMeanL0PiOnlyFit, lem::PIDExtraVars::fMeanL1Exp, lem::PIDExtraVars::fMeanL1Fit, lem::PIDExtraVars::fMeanL1PiOnlyExp, lem::PIDExtraVars::fMeanL1PiOnlyFit, lem::PIDExtraVars::fMeanMatchEExp, lem::PIDExtraVars::fMeanQFrac, lem::PIDExtraVars::fMeanQFracBkg, lem::PIDExtraVars::fMeanQFracBkgExp, lem::PIDExtraVars::fMeanQFracBkgFit, lem::PIDExtraVars::fMeanQFracExp, lem::PIDExtraVars::fMeanQFracFit, lem::PIDExtraVars::fMeanQFracSig, lem::PIDExtraVars::fMeanQFracSigExp, lem::PIDExtraVars::fMeanQFracSigFit, lem::PIDExtraVars::fMeanQFracSigY50Fit, lem::PIDExtraVars::fMeanQFracSigY90Fit, lem::PIDExtraVars::fMeanY, lem::PIDExtraVars::fMeanYExp, lem::PIDExtraVars::fMeanYFit, lem::PIDExtraVars::fMeanYSig, lem::PIDExtraVars::fMeanYSigExp, lem::PIDExtraVars::fMeanYSigFit, lem::PIDExtraVars::fMeanYSigY50Fit, lem::PIDExtraVars::fMeanYSigY90Fit, PandAna.Demos.pi0_spectra::fmt, fNCScale, lem::PIDExtraVars::fPIDExp, lem::PIDExtraVars::fPIDExpFit, lem::PIDExtraVars::fPIDFit, lem::PIDExtraVars::fPIDFitUnweight, lem::PIDExtraVars::fPIDFracNumuFit, lem::PIDExtraVars::fPIDy50Fit, lem::PIDExtraVars::fPIDy90Fit, fPositionalWeight, lem::PIDExtraVars::fQENCFracExp, lem::PIDExtraVars::fQENumuFracExp, lem::PIDExtraVars::fQESigFracExp, fReader, lem::PIDExtraVars::fResNCFracExp, lem::PIDExtraVars::fResNumuFracExp, lem::PIDExtraVars::fResSigFracExp, fScaleCalE, fSigScale, fSummarizeLabel, lem::PIDExtraVars::fTMVABDT, lem::PIDExtraVars::fTMVABDTD, lem::PIDExtraVars::fTMVABDTG, lem::PIDExtraVars::fTMVABDTG_CC, lem::PIDExtraVars::fTMVABDTMF, lem::PIDExtraVars::fTMVAMLP, lem::PIDExtraVars::fTMVAMLP_CC, fTrimMatchLists, fUseFitInput, art::DataViewImpl::getByLabel(), gTMVAVar, MECModelEnuComparisons::i, dumpEventsToText::inputs, kBkgFrac, kBkgOnly, kBkgOnlyExp, kDISNCFrac, kDISNumuFrac, kDISSigFrac, kExp, kL0, kL1, kNumuFrac, kPi0Only, kQENCFrac, kQENumuFrac, kQESigFrac, kQFrac, kResNCFrac, kResNumuFrac, kResSigFrac, kSig50Only, kSig90Only, kSigFrac, kSigFrac50, kSigFrac90, kSigOnly, kSigOnlyExp, kSimple, kTrueEVis, kUnweighted, kY, m, lem::MatchList::matches, overlay_prestage_def::matches, lem::PIDDetails::matchList, num, BlessedPlots::pid, cet::pow(), art::DataViewImpl::put(), util::sqr(), lem::dec::Evt::vars, lem::PIDDetails::vars, vars, and w.

  {
    std::unique_ptr<std::vector<lem::PIDDetails> > pidcol(new std::vector<lem::PIDDetails>);
    std::unique_ptr<art::Assns<lem::PIDDetails, rb::Cluster>> assns(new art::Assns<lem::PIDDetails, rb::Cluster>);

    art::Handle<std::vector<lem::MatchList> > matchlists;
    evt.getByLabel(fMatchesLabel, matchlists);

    art::Handle<std::vector<LEMInput> > inputs;
    evt.getByLabel(fSummarizeLabel, inputs);

    assert(inputs->size() == matchlists->size());

    art::FindManyP<rb::Cluster> fmt(matchlists, evt, fMatchesLabel);

    const unsigned int matchMax = matchlists->size();
    for(unsigned int matchIdx = 0; matchIdx < matchMax; ++matchIdx){
      MatchList matchlist = (*matchlists)[matchIdx];

      // Get all the slices associated with this match list
      std::vector<art::Ptr<rb::Cluster> > slices = fmt.at(matchIdx);
      // Associated slices can be empty when we're running over lemsum files,
      // where the slices have been stripped out. In that case, don't make
      // associations, LEMAssociator will sort it out later.
      assert(slices.empty() || slices.size() == 1);

      const double totalRecoGeV = (*inputs)[matchIdx].totalGeV;

      std::vector<MatchSummary> matches, enrichedMatches;
      for(unsigned int i = 0; i < matchlist.matches.size(); ++i){
        if(matchlist.matches[i].enrich){
          enrichedMatches.push_back(matchlist.matches[i]);
        }
        else{
          matches.push_back(matchlist.matches[i]);
        }
      }

      if(matches.empty()){
        lem::PIDDetails pid(-1);
        pidcol->push_back(pid);
        if(!slices.empty())
          util::CreateAssn(evt, *pidcol, slices[0], *assns);
        continue;
      }

      // Equalize stats
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(matches[n].isSig) matches[n].weight *= fSigScale;
        if(matches[n].ccnc == 1) matches[n].weight *= fNCScale;
      }
      for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
        enrichedMatches[n].weight *= fEnrichScale;
      }


      PIDExtraVars vars;

      std::vector<double> xs, ys, ws, ws2;

      double meanESig = 0;
      double meanEBkg = 0;
      double meanESigExp = 0;
      double meanEBkgExp = 0;
      double totSig = 0;
      double totBkg = 0;
      double totSigExp = 0;
      double totBkgExp = 0;
      double totW = 0;

      double invE = 0, invEDenom = 0;
      double expE = 0, expEDenom = 0;
      double sigE = 0, sigEDenom = 0;

      const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
      for(unsigned int i = 0; i < matches.size(); ++i){
        const double w = matches[i].weight;

        const double invEVal = w/(matches[i].potential+1e-5);
        if(matches[i].isSig) invE += invEVal; else invE -= invEVal;
        invEDenom += invEVal;

        const double expEVal = w*exp(-matches[i].potential);
        if(matches[i].isSig) expE += expEVal; else expE -= expEVal;
        expEDenom += expEVal;

        const double sigEVal = w/(1+matches[i].potential);
        if(matches[i].isSig) sigE += sigEVal; else sigE -= sigEVal;
        sigEDenom += sigEVal;

        xs.push_back(matches[i].potential);
        ys.push_back(matches[i].isSig);
        ws.push_back(w);
        double expVal = fExpTable[i];
        if(!fPositionalWeight){
          expVal = exp(-pow(matches[i].potential/maxPotential, fExpPow)/fExpLength);
        }
        ws2.push_back(w*expVal);

        if(matches[i].isSig){
          totSig += w;
          totSigExp += expVal*w;
          meanESig += w*matches[i].potential;
          meanESigExp += expVal*w*matches[i].potential;
        }
        else{
          totBkg += w;
          totBkgExp += expVal*w;
          meanEBkg += w*matches[i].potential;
          meanEBkgExp += expVal*w*matches[i].potential;
        }

        totW += w;
      } // end for i

      if(ws.size() < 2){
        vars.fChiSig = 0;
        vars.fChiBkg = 0;
        vars.fChiSigExp = 0;
        vars.fChiBkgExp = 0;
      }
      else{
        for(int c = 0; c <= 1; ++c){
          double num = 0, denom = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            num += ws[i]*(xs[i]*ys[i]-c*xs[i]);
            denom += xs[i]*xs[i]*ws[i];
          }
          const double m = num/denom;
          double chi = 0;
          double totW = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            chi += ws[i]*util::sqr(ys[i]-m*xs[i]-c);
            totW += ws[i];
          }
          chi /= totW;
          if(c == 0) vars.fChiBkg = chi;
          if(c == 1) vars.fChiSig = chi;
        } // end for c

        for(int c = 0; c <= 1; ++c){
          double num = 0, denom = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            num += ws2[i]*(xs[i]*ys[i]-c*xs[i]);
            denom += xs[i]*xs[i]*ws2[i];
          }
          const double m = num/denom;
          double chi = 0;
          double totW = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            chi += ws2[i]*util::sqr(ys[i]-m*xs[i]-c);
            totW += ws2[i];
          }
          chi /= totW;
          if(c == 0) vars.fChiBkgExp = chi;
          if(c == 1) vars.fChiSigExp = chi;
        } // end for c
      }

      if(totSig) meanESig /= totSig;
      if(totBkg) meanEBkg /= totBkg;
      if(totSigExp) meanESigExp /= totSigExp;
      if(totBkgExp) meanEBkgExp /= totBkgExp;


      double junk;

      AvgAndFit(matches, kSimple,     kSigFrac, junk, vars.fPIDFit);
      AvgAndFit(matches, kUnweighted, kSigFrac, junk, vars.fPIDFitUnweight);
      AvgAndFit(matches, kExp,        kSigFrac, vars.fPIDExp, vars.fPIDExpFit);

      AvgAndFit(matches, kSimple, kY, vars.fMeanY, vars.fMeanYFit);
      AvgAndFit(matches, kExp,    kY, vars.fMeanYExp, junk);

      AvgAndFit(matches, kSigOnly,    kY, vars.fMeanYSig, vars.fMeanYSigFit);
      AvgAndFit(matches, kSigOnlyExp, kY, vars.fMeanYSigExp, junk);

      AvgAndFit(matches, kSimple, kQFrac, vars.fMeanQFrac, vars.fMeanQFracFit);
      AvgAndFit(matches, kExp,    kQFrac, vars.fMeanQFracExp, junk);

      AvgAndFit(matches, kSigOnly,    kQFrac, vars.fMeanQFracSig, vars.fMeanQFracSigFit);
      AvgAndFit(matches, kSigOnlyExp, kQFrac, vars.fMeanQFracSigExp, junk);

      AvgAndFit(matches, kBkgOnly,    kQFrac, vars.fMeanQFracBkg, vars.fMeanQFracBkgFit);
      AvgAndFit(matches, kBkgOnlyExp, kQFrac, vars.fMeanQFracBkgExp, junk);

      AvgAndFit(matches, kExp, kTrueEVis, vars.fMeanMatchEExp, junk);

      AvgAndFit(matches, kExp, kQESigFrac,   vars.fQESigFracExp,   junk);
      AvgAndFit(matches, kExp, kResSigFrac,  vars.fResSigFracExp,  junk);
      AvgAndFit(matches, kExp, kDISSigFrac,  vars.fDISSigFracExp,  junk);
      AvgAndFit(matches, kExp, kQENCFrac,    vars.fQENCFracExp,    junk);
      AvgAndFit(matches, kExp, kResNCFrac,   vars.fResNCFracExp,   junk);
      AvgAndFit(matches, kExp, kDISNCFrac,   vars.fDISNCFracExp,   junk);
      AvgAndFit(matches, kExp, kQENumuFrac,  vars.fQENumuFracExp,  junk);
      AvgAndFit(matches, kExp, kResNumuFrac, vars.fResNumuFracExp, junk);
      AvgAndFit(matches, kExp, kDISNumuFrac, vars.fDISNumuFracExp, junk);

      // Mix in the usual signal matches
      const int emOrig = enrichedMatches.size();
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(matches[n].isSig) enrichedMatches.push_back(matches[n]);
      }
      std::sort(enrichedMatches.begin(), enrichedMatches.end());
      enrichedMatches.resize(emOrig);

      AvgAndFit(enrichedMatches, kSimple, kBkgFrac, vars.fEnrichFrac, vars.fEnrichFracFit);
      AvgAndFit(enrichedMatches, kExp,    kBkgFrac, vars.fEnrichFracExp, junk);

      AvgAndFit(enrichedMatches, kBkgOnly,    kQFrac, vars.fEnrichQFrac, vars.fEnrichQFracFit);
      AvgAndFit(enrichedMatches, kBkgOnlyExp, kQFrac, vars.fEnrichQFracExp, junk);

      AvgAndFit(enrichedMatches, kSimple, kL0, junk, vars.fMeanL0Fit);
      AvgAndFit(enrichedMatches, kSimple, kL1, junk, vars.fMeanL1Fit);
      AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0Exp, junk);
      AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1Exp, junk);

      AvgAndFit(enrichedMatches, kPi0Only, kL0, junk, vars.fMeanL0PiOnlyFit);
      AvgAndFit(enrichedMatches, kPi0Only, kL1, junk, vars.fMeanL1PiOnlyFit);
      AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0PiOnlyExp, junk);
      AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1PiOnlyExp, junk);

      AvgAndFit(matches, kSimple, kSigFrac50, junk, vars.fPIDy50Fit);
      AvgAndFit(matches, kSimple, kSigFrac90, junk, vars.fPIDy90Fit);

      AvgAndFit(matches, kSig50Only, kQFrac, junk, vars.fMeanQFracSigY50Fit);
      AvgAndFit(matches, kSig90Only, kQFrac, junk, vars.fMeanQFracSigY90Fit);

      AvgAndFit(matches, kSig50Only, kY, junk, vars.fMeanYSigY50Fit);
      AvgAndFit(matches, kSig90Only, kY, junk, vars.fMeanYSigY90Fit);

      AvgAndFit(matches, kSimple, kNumuFrac, junk, vars.fPIDFracNumuFit);

      vars.fEnergyDiff = meanESig-meanEBkg;
      vars.fEnergyDiffExp = meanESigExp-meanEBkgExp;

      vars.fBestSigPotential = vars.fBestBkgPotential = vars.fBestEnrichPotential = 2;
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(matches[n].isSig){
          vars.fBestSigPotential = matches[n].potential;
          break;
        }
      }
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(!matches[n].isSig){
          vars.fBestBkgPotential = matches[n].potential;
          break;
        }
      }
      for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
        if(!enrichedMatches[n].isSig){
          vars.fBestEnrichPotential = enrichedMatches[n].potential;
          break;
        }
      }

      vars.fAvgInvE = invE/invEDenom;
      vars.fAvgExpE = expE/expEDenom;
      vars.fAvgSigE = sigE/sigEDenom;

      // if(fANN || fANNEnrich){
      //   fann_type input[6];
      //   input[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
      //   input[1] = vars.fMeanYExp;
      //   input[2] = vars.fMeanQFracExp;
      //   input[3] = vars.fEnergyDiffExp;

      //   if(fANN){
      //     fann_type* calc_out = fann_run(fANN, input);
      //     vars.fFannId = calc_out[0];
      //   }

      //   input[4] = 2*totalRecoGeV;
      //   input[5] = vars.fEnrichFracExp;

      //   if(fANNEnrich){
      //     fann_type* calc_out = fann_run(fANNEnrich, input);
      //     vars.fFannIdEnrich = calc_out[0];
      //   }
      // }

      if(fReader){
        gTMVAVar[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
        gTMVAVar[1] = vars.fMeanYExp;
        gTMVAVar[2] = vars.fMeanQFracExp;
        gTMVAVar[3] = vars.fEnergyDiffExp;
        if(fScaleCalE)
          gTMVAVar[4] = 2*totalRecoGeV;
        else
          gTMVAVar[4] = totalRecoGeV;
        gTMVAVar[5] = vars.fEnrichFracExp;

        vars.fTMVABDTG  = fReader->EvaluateMVA("BDTG");
        vars.fTMVABDT   = fReader->EvaluateMVA("BDT");
        vars.fTMVABDTD  = fReader->EvaluateMVA("BDTD");
        vars.fTMVABDTMF = fReader->EvaluateMVA("BDTMitFisher");
        vars.fTMVAMLP   = fReader->EvaluateMVA("MLP");

        vars.fTMVABDTG_CC = fReader->EvaluateMVA("BDTG_CC");
        vars.fTMVAMLP_CC  = fReader->EvaluateMVA("MLP_CC");
      }

      if(fForest){
        if(vars.fPIDExp == -1){
          vars.fDecTreePID = -1;
        }
        else{
          dec::Evt decevt;
          decevt.vars[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
          decevt.vars[1] = vars.fMeanYExp;
          decevt.vars[2] = vars.fMeanQFracExp;
          decevt.vars[3] = vars.fEnergyDiffExp;
          if(fScaleCalE)
            decevt.vars[4] = 2*totalRecoGeV;
          else
            decevt.vars[4] = totalRecoGeV;
          decevt.vars[5] = vars.fEnrichFracExp;
          vars.fDecTreePID = fForest->Classify(decevt);
        }
      }

      double defaultPID = vars.fDecTreePID;

      // Flattened version of the Decision Tree variable
      if(fFlatHist) defaultPID = fFlatHist->GetBinContent(fFlatHist->FindBin(vars.fDecTreePID));


      lem::PIDDetails pid(defaultPID);

      if(fTrimMatchLists) matchlist.matches.resize(1);

      // TODO: sample size equalization doesn't survive this
      pid.matchList = matchlist;
      pid.vars = vars;

      pidcol->push_back(pid);

      if(!slices.empty()){
        util::CreateAssn(evt, *pidcol, slices[0], *assns);
      }
    } // end for matchIdx

    evt.put(std::move(pidcol));
    evt.put(std::move(assns));
  }

void lem::MakePID::reconfigure ( const fhicl::ParameterSet &  pset )

virtual

Definition at line 195 of file MakePID_module.cc.

References fExpLength, fExpPow, fMatchesLabel, fPositionalWeight, fScaleCalE, fSummarizeLabel, fTrimMatchLists, fUseFitInput, fhicl::ParameterSet::get(), and string.

Referenced by MakePID().

  {
    fScaleCalE = pset.get<bool>("ScaleCalE");

    fTrimMatchLists = pset.get<bool>("TrimMatchLists");
    fUseFitInput = pset.get<bool>("UseFitInput");
    fPositionalWeight = pset.get<bool>("PositionalWeight");
    fMatchesLabel = pset.get<std::string>("MatchesLabel");
    fSummarizeLabel = pset.get<std::string>("SummarizeLabel");
    fExpLength = pset.get<double>("ExpLength");
    fExpPow = pset.get<double>("ExpPow");
  }

void art::ModuleBase::setModuleDescription ( ModuleDescription const &  )

inherited

void art::ModuleBase::sortConsumables ( std::string const &  current_process_name )

inherited

std::string art::EDProducer::workerType ( ) const

inherited

Referenced by art::EDProducer::EDProducer().

Member Data Documentation

double lem::MakePID::fEnrichScale

protected

Definition at line 105 of file MakePID_module.cc.

Referenced by beginRun(), and produce().

double lem::MakePID::fExpLength

protected

Definition at line 107 of file MakePID_module.cc.

Referenced by AvgAndFit(), beginJob(), produce(), and reconfigure().

double lem::MakePID::fExpPow

protected

Definition at line 107 of file MakePID_module.cc.

Referenced by AvgAndFit(), beginJob(), produce(), and reconfigure().

std::vector<double> lem::MakePID::fExpTable

protected

Definition at line 101 of file MakePID_module.cc.

Referenced by AvgAndFit(), beginJob(), and produce().

TH1* lem::MakePID::fFlatHist

protected

Definition at line 120 of file MakePID_module.cc.

Referenced by MakePID(), and produce().

dec::Forest* lem::MakePID::fForest

protected

Definition at line 118 of file MakePID_module.cc.

Referenced by MakePID(), produce(), and ~MakePID().

std::string lem::MakePID::fMatchesLabel

protected

Definition at line 114 of file MakePID_module.cc.

Referenced by beginRun(), produce(), and reconfigure().

double lem::MakePID::fNCScale

protected

Definition at line 103 of file MakePID_module.cc.

Referenced by beginRun(), and produce().

bool lem::MakePID::fPositionalWeight

protected

Definition at line 113 of file MakePID_module.cc.

Referenced by AvgAndFit(), produce(), and reconfigure().

TMVA::Reader* lem::MakePID::fReader

protected

Definition at line 117 of file MakePID_module.cc.

Referenced by MakePID(), produce(), and ~MakePID().

bool lem::MakePID::fScaleCalE

protected

Definition at line 60 of file MakePID_module.cc.

Referenced by produce(), and reconfigure().

double lem::MakePID::fSigScale

protected

Definition at line 104 of file MakePID_module.cc.

Referenced by beginRun(), and produce().

std::string lem::MakePID::fSummarizeLabel

protected

Definition at line 115 of file MakePID_module.cc.

Referenced by produce(), and reconfigure().

bool lem::MakePID::fTrimMatchLists

protected

Definition at line 111 of file MakePID_module.cc.

Referenced by produce(), and reconfigure().

bool lem::MakePID::fUseFitInput

protected

Definition at line 112 of file MakePID_module.cc.

Referenced by produce(), and reconfigure().

---

The documentation for this class was generated from the following file:

• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-04-21/LEM/MakePID_module.cc