Public Types | Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
lem::MakePID Class Reference

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

Inheritance diagram for lem::MakePID:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

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< TconsumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< TmayConsumeView (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< Tconsumes (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< TmayConsume (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::ForestfForest
 
TH1 * fFlatHist
 

Detailed Description

Calculate final PID variables from match details.

Definition at line 48 of file MakePID_module.cc.

Member Typedef Documentation

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.

Enumerator
kUnweighted 
kSimple 
kExp 
kSigOnly 
kSigOnlyExp 
kBkgOnly 
kBkgOnlyExp 
kSig50Only 
kSig90Only 
kPi0Only 

Definition at line 62 of file MakePID_module.cc.

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.

125  : EDProducer(pset)
126  {
127  reconfigure(pset);
128 
129  produces<std::vector<lem::PIDDetails>>();
130  produces<art::Assns<lem::PIDDetails, rb::Cluster>>();
131 
132  const std::string path = pset.get<std::string>("LibraryPath");
133  const std::string libPath = util::EnvExpansion(path) + "/";
134 
135  fReader = 0;
136  if(pset.get<bool>("FillTMVAs")){
137  fReader = new TMVA::Reader;
138 
139  fReader->AddVariable("pidExtrap", &gTMVAVar[0]); // vars.fPIDFit
140  fReader->AddVariable("yexp", &gTMVAVar[1]); // vars.fMeanYExp
141  fReader->AddVariable("qfracexp", &gTMVAVar[2]); // vars.fMeanQFracExp
142  fReader->AddVariable("ediff", &gTMVAVar[3]); // vars.fEnergyDiffExp
143  fReader->AddVariable("calE", &gTMVAVar[4]); // 2*totalRecoGeV
144  fReader->AddVariable("enrichExp", &gTMVAVar[5]); // vars.fEnrichFracExp
145 
146  fReader->AddSpectator("isSig", &gJunk);
147  fReader->AddSpectator("myW", &gJunk);
148 
149  const std::string tmvaPath = libPath+"weights/";
150 
151  fReader->BookMVA("BDTG", tmvaPath+"TMVAClassification_BDTG.weights.xml");
152  fReader->BookMVA("BDT", tmvaPath+"TMVAClassification_BDT.weights.xml");
153  fReader->BookMVA("BDTD", tmvaPath+"TMVAClassification_BDTD.weights.xml");
154  fReader->BookMVA("BDTMitFisher", tmvaPath+"TMVAClassification_BDTMitFisher.weights.xml");
155  fReader->BookMVA("MLP", tmvaPath+"TMVAClassification_MLP.weights.xml");
156 
157  fReader->BookMVA("BDTG_CC", tmvaPath+"TMVAClassification_BDTG_CC.weights.xml");
158  fReader->BookMVA("MLP_CC", tmvaPath+"TMVAClassification_MLP_CC.weights.xml");
159  }
160 
161  // fANN = fANNEnrich = 0;
162  // if(pset.get<bool>("FillANNs")){
163  // fANN = fann_create_from_file((libPath+"fann.net").c_str());
164  // assert(fANN);
165  // fann_print_connections(fANN);
166 
167  // fANNEnrich = fann_create_from_file((libPath+"fann_enrich.net").c_str());
168  // assert(fANNEnrich);
169  // fann_print_connections(fANNEnrich);
170  // }
171 
172  fForest = 0;
173  if(pset.get<bool>("FillDecTree")){
174  std::cout << "Loading LEM dectree from " << libPath << std::endl;
175  fForest = dec::Forest::FromFile(libPath+"dectree.root");
176  }
177 
178  fFlatHist = 0;
179  if(pset.get<bool>("Flatten")){
180  TFile* fFlat = new TFile((libPath+"flatten.root").c_str());
181  assert(!fFlat->IsZombie());
182  fFlatHist = (TH1*)fFlat->Get("conv");
183  assert(fFlatHist);
184  }
185  }
virtual void reconfigure(const fhicl::ParameterSet &pset)
static Forest * FromFile(const std::string &fname)
Load PID from a file.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
dec::Forest * fForest
float gTMVAVar[6]
T get(std::string const &key) const
Definition: ParameterSet.h:231
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
assert(nhit_max >=nhit_nbins)
float gJunk
TMVA::Reader * fReader
enum BeamMode string
lem::MakePID::~MakePID ( )

Definition at line 188 of file MakePID_module.cc.

References fForest, and fReader.

189  {
190  delete fReader;
191  delete fForest;
192  }
dec::Forest * fForest
TMVA::Reader * fReader

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().

239  {
240  std::vector<double> xs, ys, ws;
241  double tot = 0;
242  mean = 0;
243 
244  const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
245 
246  for(unsigned int i = 0; i < matches.size(); ++i){
247  const MatchSummary& match = matches[i];
248 
249  double w = match.weight;
250 
251  if(weightScheme == kUnweighted) w = 1;
252 
253  if(weightScheme == kExp ||
254  weightScheme == kSigOnlyExp ||
255  weightScheme == kBkgOnlyExp){
256  if(fPositionalWeight){
257  w *= fExpTable[i];
258  }
259  else{
260  w *= exp(-pow(match.potential/maxPotential, fExpPow)/fExpLength);
261  }
262  }
263 
264  if(weightScheme == kSigOnly || weightScheme == kSigOnlyExp){
265  if(!match.isSig) w = 0;
266  }
267 
268  if(weightScheme == kBkgOnly || weightScheme == kBkgOnlyExp){
269  if(match.isSig) w = 0;
270  }
271 
272  if(weightScheme == kSig50Only){
273  if(!match.isSig || match.y > 0.50) w = 0;
274  }
275  if(weightScheme == kSig90Only){
276  if(!match.isSig || match.y > 0.90) w = 0;
277  }
278  if(weightScheme == kPi0Only) w = (match.photL0 > 0 || match.photL1 > 0 || match.photE0 > 0 || match.photE1 > 0);
279 
280  double var = 0;
281  switch(matchVar){
282  case kSigFrac: var = int(match.isSig); break;
283  case kBkgFrac: var = 1-int(match.isSig); break;
284  case kNumuFrac: var = int(match.ccnc == simb::kCC && abs(match.pdg) == 14); break;
285  case kY: var = match.y; break;
286  case kQFrac: var = match.qFrac; break;
287  case kSigFrac50: var = int(match.isSig && match.y < .50); break;
288  case kSigFrac90: var = int(match.isSig && match.y < .90); break;
289  case kL0: var = match.photL0; break;
290  case kL1: var = match.photL1; break;
291  case kTrueEVis: var = match.trueEVis; break;
292  case kQESigFrac: var = (match.isSig && match.mode == simb::kQE); break;
293  case kResSigFrac: var = (match.isSig && match.mode == simb::kRes); break;
294  case kDISSigFrac: var = (match.isSig && match.mode == simb::kDIS); break;
295  case kQENCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kQE); break;
296  case kResNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kRes); break;
297  case kDISNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kDIS); break;
298  case kQENumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kQE); break;
299  case kResNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kRes); break;
300  case kDISNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kDIS); break;
301  }
302 
303  xs.push_back(match.potential);
304  ys.push_back(var);
305  ws.push_back(w);
306 
307  tot += w;
308  mean += w*var;
309  } // end for i
310 
311  mean /= tot;
312 
313  if(ws.size() < 2){
314  fit = -1;
315  }
316  else{
317  double m, c;
318  if(weightScheme == kUnweighted)
319  util::LinFitUnweighted(xs, ys, m, c);
320  else
321  util::LinFit(xs, ys, ws, m, c);
322  fit = m*xs[0]+c;
323  }
324  }
constexpr T pow(T x)
Definition: pow.h:72
var_value< double > var
Definition: StanTypedefs.h:14
void abs(TH1 *hist)
std::vector< double > fExpTable
void LinFitUnweighted(const std::vector< double > &x, const std::vector< double > &y, double &m, double &c)
Simplified version of LinFit.
Definition: MathUtil.cxx:8
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
Float_t w
Definition: plot.C:20
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().

210  {
211  fExpTable.resize(kMaxNumMatches);
212  for(unsigned int i = 0; i < kMaxNumMatches; ++i){
213  fExpTable[i] = exp(-pow(i/fExpLength, fExpPow));
214  }
215  }
const unsigned int kMaxNumMatches
Definition: FindMatches.h:29
constexpr T pow(T x)
Definition: pow.h:72
std::vector< double > fExpTable
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.

219  {
221  run.getByLabel(fMatchesLabel, libsum);
222 
223  assert(!libsum.failedToGet());
224 
225  fNCScale = libsum->nTrueNC/double(libsum->nTrueNC+libsum->nSwapNC);
226  fSigScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totSig;
227  fEnrichScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totEnrich;
228 
229  std::cout << "Reweighted to sig, bkg, enrich = "
230  << libsum->totSig*fSigScale << ", "
231  << libsum->totBkg+(fNCScale-1)*libsum->nTrueNC-libsum->nSwapNC << ", "
232  << libsum->totEnrich*fEnrichScale;
233  }
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
OStream cout
Definition: OStream.cxx:6
std::string fMatchesLabel
assert(nhit_max >=nhit_nbins)
bool failedToGet() const
Definition: Handle.h:190
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.

56  {
57  return collector_.consumes<T, BT>(tag);
58  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
ProductToken< T > consumes(InputTag const &)
double T
Definition: Xdiff_gwt.C:5
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.

70  {
71  collector_.consumesMany<T, BT>();
72  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
double T
Definition: Xdiff_gwt.C:5
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.

63  {
64  return collector_.consumesView<T, BT>(tag);
65  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
ViewToken< Element > consumesView(InputTag const &)
double T
Definition: Xdiff_gwt.C:5
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.

77  {
78  return collector_.mayConsume<T, BT>(tag);
79  }
ProductToken< T > mayConsume(InputTag const &)
ConsumesCollector collector_
Definition: ModuleBase.h:50
double T
Definition: Xdiff_gwt.C:5
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.

91  {
93  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
double T
Definition: Xdiff_gwt.C:5
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.

84  {
85  return collector_.mayConsumeView<T, BT>(tag);
86  }
ConsumesCollector collector_
Definition: ModuleBase.h:50
ViewToken< Element > mayConsumeView(InputTag const &)
double T
Definition: Xdiff_gwt.C:5
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.

328  {
329  std::unique_ptr<std::vector<lem::PIDDetails> > pidcol(new std::vector<lem::PIDDetails>);
330  std::unique_ptr<art::Assns<lem::PIDDetails, rb::Cluster>> assns(new art::Assns<lem::PIDDetails, rb::Cluster>);
331 
333  evt.getByLabel(fMatchesLabel, matchlists);
334 
336  evt.getByLabel(fSummarizeLabel, inputs);
337 
338  assert(inputs->size() == matchlists->size());
339 
341 
342  const unsigned int matchMax = matchlists->size();
343  for(unsigned int matchIdx = 0; matchIdx < matchMax; ++matchIdx){
344  MatchList matchlist = (*matchlists)[matchIdx];
345 
346  // Get all the slices associated with this match list
347  std::vector<art::Ptr<rb::Cluster> > slices = fmt.at(matchIdx);
348  // Associated slices can be empty when we're running over lemsum files,
349  // where the slices have been stripped out. In that case, don't make
350  // associations, LEMAssociator will sort it out later.
351  assert(slices.empty() || slices.size() == 1);
352 
353  const double totalRecoGeV = (*inputs)[matchIdx].totalGeV;
354 
355  std::vector<MatchSummary> matches, enrichedMatches;
356  for(unsigned int i = 0; i < matchlist.matches.size(); ++i){
357  if(matchlist.matches[i].enrich){
358  enrichedMatches.push_back(matchlist.matches[i]);
359  }
360  else{
361  matches.push_back(matchlist.matches[i]);
362  }
363  }
364 
365  if(matches.empty()){
366  lem::PIDDetails pid(-1);
367  pidcol->push_back(pid);
368  if(!slices.empty())
369  util::CreateAssn(evt, *pidcol, slices[0], *assns);
370  continue;
371  }
372 
373  // Equalize stats
374  for(unsigned int n = 0; n < matches.size(); ++n){
375  if(matches[n].isSig) matches[n].weight *= fSigScale;
376  if(matches[n].ccnc == 1) matches[n].weight *= fNCScale;
377  }
378  for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
379  enrichedMatches[n].weight *= fEnrichScale;
380  }
381 
382 
383  PIDExtraVars vars;
384 
385  std::vector<double> xs, ys, ws, ws2;
386 
387  double meanESig = 0;
388  double meanEBkg = 0;
389  double meanESigExp = 0;
390  double meanEBkgExp = 0;
391  double totSig = 0;
392  double totBkg = 0;
393  double totSigExp = 0;
394  double totBkgExp = 0;
395  double totW = 0;
396 
397  double invE = 0, invEDenom = 0;
398  double expE = 0, expEDenom = 0;
399  double sigE = 0, sigEDenom = 0;
400 
401  const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
402  for(unsigned int i = 0; i < matches.size(); ++i){
403  const double w = matches[i].weight;
404 
405  const double invEVal = w/(matches[i].potential+1e-5);
406  if(matches[i].isSig) invE += invEVal; else invE -= invEVal;
407  invEDenom += invEVal;
408 
409  const double expEVal = w*exp(-matches[i].potential);
410  if(matches[i].isSig) expE += expEVal; else expE -= expEVal;
411  expEDenom += expEVal;
412 
413  const double sigEVal = w/(1+matches[i].potential);
414  if(matches[i].isSig) sigE += sigEVal; else sigE -= sigEVal;
415  sigEDenom += sigEVal;
416 
417  xs.push_back(matches[i].potential);
418  ys.push_back(matches[i].isSig);
419  ws.push_back(w);
420  double expVal = fExpTable[i];
421  if(!fPositionalWeight){
422  expVal = exp(-pow(matches[i].potential/maxPotential, fExpPow)/fExpLength);
423  }
424  ws2.push_back(w*expVal);
425 
426  if(matches[i].isSig){
427  totSig += w;
428  totSigExp += expVal*w;
429  meanESig += w*matches[i].potential;
430  meanESigExp += expVal*w*matches[i].potential;
431  }
432  else{
433  totBkg += w;
434  totBkgExp += expVal*w;
435  meanEBkg += w*matches[i].potential;
436  meanEBkgExp += expVal*w*matches[i].potential;
437  }
438 
439  totW += w;
440  } // end for i
441 
442  if(ws.size() < 2){
443  vars.fChiSig = 0;
444  vars.fChiBkg = 0;
445  vars.fChiSigExp = 0;
446  vars.fChiBkgExp = 0;
447  }
448  else{
449  for(int c = 0; c <= 1; ++c){
450  double num = 0, denom = 0;
451  for(unsigned int i = 0; i < xs.size(); ++i){
452  num += ws[i]*(xs[i]*ys[i]-c*xs[i]);
453  denom += xs[i]*xs[i]*ws[i];
454  }
455  const double m = num/denom;
456  double chi = 0;
457  double totW = 0;
458  for(unsigned int i = 0; i < xs.size(); ++i){
459  chi += ws[i]*util::sqr(ys[i]-m*xs[i]-c);
460  totW += ws[i];
461  }
462  chi /= totW;
463  if(c == 0) vars.fChiBkg = chi;
464  if(c == 1) vars.fChiSig = chi;
465  } // end for c
466 
467  for(int c = 0; c <= 1; ++c){
468  double num = 0, denom = 0;
469  for(unsigned int i = 0; i < xs.size(); ++i){
470  num += ws2[i]*(xs[i]*ys[i]-c*xs[i]);
471  denom += xs[i]*xs[i]*ws2[i];
472  }
473  const double m = num/denom;
474  double chi = 0;
475  double totW = 0;
476  for(unsigned int i = 0; i < xs.size(); ++i){
477  chi += ws2[i]*util::sqr(ys[i]-m*xs[i]-c);
478  totW += ws2[i];
479  }
480  chi /= totW;
481  if(c == 0) vars.fChiBkgExp = chi;
482  if(c == 1) vars.fChiSigExp = chi;
483  } // end for c
484  }
485 
486  if(totSig) meanESig /= totSig;
487  if(totBkg) meanEBkg /= totBkg;
488  if(totSigExp) meanESigExp /= totSigExp;
489  if(totBkgExp) meanEBkgExp /= totBkgExp;
490 
491 
492  double junk;
493 
494  AvgAndFit(matches, kSimple, kSigFrac, junk, vars.fPIDFit);
495  AvgAndFit(matches, kUnweighted, kSigFrac, junk, vars.fPIDFitUnweight);
496  AvgAndFit(matches, kExp, kSigFrac, vars.fPIDExp, vars.fPIDExpFit);
497 
498  AvgAndFit(matches, kSimple, kY, vars.fMeanY, vars.fMeanYFit);
499  AvgAndFit(matches, kExp, kY, vars.fMeanYExp, junk);
500 
501  AvgAndFit(matches, kSigOnly, kY, vars.fMeanYSig, vars.fMeanYSigFit);
502  AvgAndFit(matches, kSigOnlyExp, kY, vars.fMeanYSigExp, junk);
503 
504  AvgAndFit(matches, kSimple, kQFrac, vars.fMeanQFrac, vars.fMeanQFracFit);
505  AvgAndFit(matches, kExp, kQFrac, vars.fMeanQFracExp, junk);
506 
507  AvgAndFit(matches, kSigOnly, kQFrac, vars.fMeanQFracSig, vars.fMeanQFracSigFit);
508  AvgAndFit(matches, kSigOnlyExp, kQFrac, vars.fMeanQFracSigExp, junk);
509 
510  AvgAndFit(matches, kBkgOnly, kQFrac, vars.fMeanQFracBkg, vars.fMeanQFracBkgFit);
511  AvgAndFit(matches, kBkgOnlyExp, kQFrac, vars.fMeanQFracBkgExp, junk);
512 
513  AvgAndFit(matches, kExp, kTrueEVis, vars.fMeanMatchEExp, junk);
514 
515  AvgAndFit(matches, kExp, kQESigFrac, vars.fQESigFracExp, junk);
516  AvgAndFit(matches, kExp, kResSigFrac, vars.fResSigFracExp, junk);
517  AvgAndFit(matches, kExp, kDISSigFrac, vars.fDISSigFracExp, junk);
518  AvgAndFit(matches, kExp, kQENCFrac, vars.fQENCFracExp, junk);
519  AvgAndFit(matches, kExp, kResNCFrac, vars.fResNCFracExp, junk);
520  AvgAndFit(matches, kExp, kDISNCFrac, vars.fDISNCFracExp, junk);
521  AvgAndFit(matches, kExp, kQENumuFrac, vars.fQENumuFracExp, junk);
522  AvgAndFit(matches, kExp, kResNumuFrac, vars.fResNumuFracExp, junk);
523  AvgAndFit(matches, kExp, kDISNumuFrac, vars.fDISNumuFracExp, junk);
524 
525  // Mix in the usual signal matches
526  const int emOrig = enrichedMatches.size();
527  for(unsigned int n = 0; n < matches.size(); ++n){
528  if(matches[n].isSig) enrichedMatches.push_back(matches[n]);
529  }
530  std::sort(enrichedMatches.begin(), enrichedMatches.end());
531  enrichedMatches.resize(emOrig);
532 
533  AvgAndFit(enrichedMatches, kSimple, kBkgFrac, vars.fEnrichFrac, vars.fEnrichFracFit);
534  AvgAndFit(enrichedMatches, kExp, kBkgFrac, vars.fEnrichFracExp, junk);
535 
536  AvgAndFit(enrichedMatches, kBkgOnly, kQFrac, vars.fEnrichQFrac, vars.fEnrichQFracFit);
537  AvgAndFit(enrichedMatches, kBkgOnlyExp, kQFrac, vars.fEnrichQFracExp, junk);
538 
539  AvgAndFit(enrichedMatches, kSimple, kL0, junk, vars.fMeanL0Fit);
540  AvgAndFit(enrichedMatches, kSimple, kL1, junk, vars.fMeanL1Fit);
541  AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0Exp, junk);
542  AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1Exp, junk);
543 
544  AvgAndFit(enrichedMatches, kPi0Only, kL0, junk, vars.fMeanL0PiOnlyFit);
545  AvgAndFit(enrichedMatches, kPi0Only, kL1, junk, vars.fMeanL1PiOnlyFit);
546  AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0PiOnlyExp, junk);
547  AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1PiOnlyExp, junk);
548 
549  AvgAndFit(matches, kSimple, kSigFrac50, junk, vars.fPIDy50Fit);
550  AvgAndFit(matches, kSimple, kSigFrac90, junk, vars.fPIDy90Fit);
551 
552  AvgAndFit(matches, kSig50Only, kQFrac, junk, vars.fMeanQFracSigY50Fit);
553  AvgAndFit(matches, kSig90Only, kQFrac, junk, vars.fMeanQFracSigY90Fit);
554 
555  AvgAndFit(matches, kSig50Only, kY, junk, vars.fMeanYSigY50Fit);
556  AvgAndFit(matches, kSig90Only, kY, junk, vars.fMeanYSigY90Fit);
557 
558  AvgAndFit(matches, kSimple, kNumuFrac, junk, vars.fPIDFracNumuFit);
559 
560  vars.fEnergyDiff = meanESig-meanEBkg;
561  vars.fEnergyDiffExp = meanESigExp-meanEBkgExp;
562 
563  vars.fBestSigPotential = vars.fBestBkgPotential = vars.fBestEnrichPotential = 2;
564  for(unsigned int n = 0; n < matches.size(); ++n){
565  if(matches[n].isSig){
566  vars.fBestSigPotential = matches[n].potential;
567  break;
568  }
569  }
570  for(unsigned int n = 0; n < matches.size(); ++n){
571  if(!matches[n].isSig){
572  vars.fBestBkgPotential = matches[n].potential;
573  break;
574  }
575  }
576  for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
577  if(!enrichedMatches[n].isSig){
578  vars.fBestEnrichPotential = enrichedMatches[n].potential;
579  break;
580  }
581  }
582 
583  vars.fAvgInvE = invE/invEDenom;
584  vars.fAvgExpE = expE/expEDenom;
585  vars.fAvgSigE = sigE/sigEDenom;
586 
587  // if(fANN || fANNEnrich){
588  // fann_type input[6];
589  // input[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
590  // input[1] = vars.fMeanYExp;
591  // input[2] = vars.fMeanQFracExp;
592  // input[3] = vars.fEnergyDiffExp;
593 
594  // if(fANN){
595  // fann_type* calc_out = fann_run(fANN, input);
596  // vars.fFannId = calc_out[0];
597  // }
598 
599  // input[4] = 2*totalRecoGeV;
600  // input[5] = vars.fEnrichFracExp;
601 
602  // if(fANNEnrich){
603  // fann_type* calc_out = fann_run(fANNEnrich, input);
604  // vars.fFannIdEnrich = calc_out[0];
605  // }
606  // }
607 
608  if(fReader){
609  gTMVAVar[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
610  gTMVAVar[1] = vars.fMeanYExp;
611  gTMVAVar[2] = vars.fMeanQFracExp;
612  gTMVAVar[3] = vars.fEnergyDiffExp;
613  if(fScaleCalE)
614  gTMVAVar[4] = 2*totalRecoGeV;
615  else
616  gTMVAVar[4] = totalRecoGeV;
617  gTMVAVar[5] = vars.fEnrichFracExp;
618 
619  vars.fTMVABDTG = fReader->EvaluateMVA("BDTG");
620  vars.fTMVABDT = fReader->EvaluateMVA("BDT");
621  vars.fTMVABDTD = fReader->EvaluateMVA("BDTD");
622  vars.fTMVABDTMF = fReader->EvaluateMVA("BDTMitFisher");
623  vars.fTMVAMLP = fReader->EvaluateMVA("MLP");
624 
625  vars.fTMVABDTG_CC = fReader->EvaluateMVA("BDTG_CC");
626  vars.fTMVAMLP_CC = fReader->EvaluateMVA("MLP_CC");
627  }
628 
629  if(fForest){
630  if(vars.fPIDExp == -1){
631  vars.fDecTreePID = -1;
632  }
633  else{
634  dec::Evt decevt;
635  decevt.vars[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
636  decevt.vars[1] = vars.fMeanYExp;
637  decevt.vars[2] = vars.fMeanQFracExp;
638  decevt.vars[3] = vars.fEnergyDiffExp;
639  if(fScaleCalE)
640  decevt.vars[4] = 2*totalRecoGeV;
641  else
642  decevt.vars[4] = totalRecoGeV;
643  decevt.vars[5] = vars.fEnrichFracExp;
644  vars.fDecTreePID = fForest->Classify(decevt);
645  }
646  }
647 
648  double defaultPID = vars.fDecTreePID;
649 
650  // Flattened version of the Decision Tree variable
651  if(fFlatHist) defaultPID = fFlatHist->GetBinContent(fFlatHist->FindBin(vars.fDecTreePID));
652 
653 
654  lem::PIDDetails pid(defaultPID);
655 
656  if(fTrimMatchLists) matchlist.matches.resize(1);
657 
658  // TODO: sample size equalization doesn't survive this
659  pid.matchList = matchlist;
660  pid.vars = vars;
661 
662  pidcol->push_back(pid);
663 
664  if(!slices.empty()){
665  util::CreateAssn(evt, *pidcol, slices[0], *assns);
666  }
667  } // end for matchIdx
668 
669  evt.put(std::move(pidcol));
670  evt.put(std::move(assns));
671  }
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
constexpr T pow(T x)
Definition: pow.h:72
std::vector< double > fExpTable
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
dec::Forest * fForest
double Classify(const Evt &evt) const
Calculate the PID value for evt.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
float gTMVAVar[6]
std::void_t< T > n
const std::map< std::pair< std::string, std::string >, Variable > vars
std::string fSummarizeLabel
std::string fMatchesLabel
int num
Definition: f2_nu.C:119
assert(nhit_max >=nhit_nbins)
Attach LEM-specific info to the base PID object.
Definition: PIDDetails.h:20
Float_t e
Definition: plot.C:35
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
void AvgAndFit(const std::vector< MatchSummary > &matches, WeightScheme weightScheme, MatchVar matchVar, double &mean, double &fit) const
Float_t w
Definition: plot.C:20
TMVA::Reader * fReader
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().

196  {
197  fScaleCalE = pset.get<bool>("ScaleCalE");
198 
199  fTrimMatchLists = pset.get<bool>("TrimMatchLists");
200  fUseFitInput = pset.get<bool>("UseFitInput");
201  fPositionalWeight = pset.get<bool>("PositionalWeight");
202  fMatchesLabel = pset.get<std::string>("MatchesLabel");
203  fSummarizeLabel = pset.get<std::string>("SummarizeLabel");
204  fExpLength = pset.get<double>("ExpLength");
205  fExpPow = pset.get<double>("ExpPow");
206  }
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fSummarizeLabel
std::string fMatchesLabel
enum BeamMode string
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

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: