MakePID_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file MakePID_module.cc
3 // \brief Calculate final PID variables from match details
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "LEM/PIDDetails.h"
8 
9 #include "LEM/LEMInput.h"
10 #include "LEM/func/FindMatches.h" // For num matches constant
12 #include "LEM/func/MatchList.h"
13 
14 #include "RecoBase/CellHit.h"
15 #include "RecoBase/Cluster.h"
16 #include "nusimdata/SimulationBase/MCTruth.h" // mode enum
17 #include "Utilities/AssociationUtil.h"
19 #include "Utilities/func/MathUtil.h"
20 
21 //#include "doublefann.h" // from the FANN external
22 
23 #include "TMVA/Reader.h"
24 
25 #include "LEM/func/DecisionTree.h"
26 
27 // Framework includes
36 #include "fhiclcpp/ParameterSet.h"
38 
39 float gTMVAVar[6];
40 float gJunk;
41 
42 #include "TFile.h"
43 #include "TH2.h"
44 
45 namespace lem
46 {
47  /// Calculate final %PID variables from match details
48  class MakePID: public art::EDProducer
49  {
50  public:
51  explicit MakePID(const fhicl::ParameterSet& pset);
52  ~MakePID();
53 
54  virtual void beginJob();
55  virtual void beginRun(art::Run& run);
56  virtual void produce(art::Event& evt);
57  virtual void reconfigure(const fhicl::ParameterSet& pset);
58 
59  protected:
60  bool fScaleCalE;
61 
73  };
74 
75  enum MatchVar{
79  kY,
83  kL0,
84  kL1,
95  };
96 
97  void AvgAndFit(const std::vector<MatchSummary>& matches,
98  WeightScheme weightScheme, MatchVar matchVar,
99  double& mean, double& fit) const;
100 
101  std::vector<double> fExpTable;
102 
103  double fNCScale;
104  double fSigScale;
105  double fEnrichScale;
106 
108 
109  // fann *fANN, *fANNEnrich;
110 
116 
117  TMVA::Reader* fReader;
119 
120  TH1* fFlatHist;
121  };
122 
123  //.......................................................................
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  }
186 
187  //......................................................................
189  {
190  delete fReader;
191  delete fForest;
192  }
193 
194  //......................................................................
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  }
207 
208  //......................................................................
210  {
211  fExpTable.resize(kMaxNumMatches);
212  for(unsigned int i = 0; i < kMaxNumMatches; ++i){
213  fExpTable[i] = exp(-pow(i/fExpLength, fExpPow));
214  }
215  }
216 
217  //......................................................................
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  }
234 
235  //......................................................................
236  void MakePID::AvgAndFit(const std::vector<MatchSummary>& matches,
237  WeightScheme weightScheme, MatchVar matchVar,
238  double& mean, double& fit) const
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  }
325 
326  //......................................................................
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 
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 
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  }
672 
674 
675 } // end namespace lem
676 ////////////////////////////////////////////////////////////////////////
Decision Tree PID.
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.
std::vector< MatchSummary > matches
Definition: MatchList.h:24
double fBestBkgPotential
Definition: PIDExtraVars.h:69
double fMeanL1PiOnlyExp
Definition: PIDExtraVars.h:97
virtual void reconfigure(const fhicl::ParameterSet &pset)
double fMeanQFracSigY90Fit
Definition: PIDExtraVars.h:90
const unsigned int kMaxNumMatches
Definition: FindMatches.h:29
double fBestSigPotential
Definition: PIDExtraVars.h:68
static Forest * FromFile(const std::string &fname)
Load PID from a file.
double fPIDFitUnweight
Definition: PIDExtraVars.h:26
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
constexpr T pow(T x)
Definition: pow.h:72
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
var_value< double > var
Definition: StanTypedefs.h:14
void abs(TH1 *hist)
double fMeanQFracBkgExp
Definition: PIDExtraVars.h:47
double fMeanL1PiOnlyFit
Definition: PIDExtraVars.h:95
std::vector< double > fExpTable
double fMeanYSigY90Fit
Definition: PIDExtraVars.h:92
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Definition: Run.h:21
double fMeanQFracSigFit
Definition: PIDExtraVars.h:43
A training or trial event for the decision tree.
Definition: DecisionTree.h:28
Collection of MatchSummary objects.
PID
Definition: FillPIDs.h:14
Core of the LEM match-finding algorithm.
dec::Forest * fForest
double Classify(const Evt &evt) const
Calculate the PID value for evt.
Calculate final PID variables from match details.
double fMeanQFracSigY50Fit
Definition: PIDExtraVars.h:89
Details of the library LEM matches were made against.
Simplified Match information, suitable for serialization.
Definition: MatchSummary.h:16
virtual void produce(art::Event &evt)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
Collection of MatchSummary objects.
Definition: MatchList.h:17
float gTMVAVar[6]
double fMeanYSigY50Fit
Definition: PIDExtraVars.h:91
double fMeanL0PiOnlyFit
Definition: PIDExtraVars.h:95
Outputs of the MakePID module.
Definition: PIDExtraVars.h:15
std::void_t< T > n
"Random forest" of decision trees
Definition: DecisionTree.h:42
MatchList matchList
Definition: PIDDetails.h:31
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
double fMeanQFracSigExp
Definition: PIDExtraVars.h:44
const std::map< std::pair< std::string, std::string >, Variable > vars
double fMeanQFracBkgFit
Definition: PIDExtraVars.h:46
double fEnrichQFracFit
Definition: PIDExtraVars.h:57
MakePID(const fhicl::ParameterSet &pset)
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
PIDExtraVars vars
Definition: PIDDetails.h:32
double fPIDFracNumuFit
Definition: PIDExtraVars.h:85
void LinFitUnweighted(const std::vector< double > &x, const std::vector< double > &y, double &m, double &c)
Simplified version of LinFit.
Definition: MathUtil.cxx:8
std::string fSummarizeLabel
std::string fMatchesLabel
double fMeanL0PiOnlyExp
Definition: PIDExtraVars.h:97
int num
Definition: f2_nu.C:119
double fEnrichQFracExp
Definition: PIDExtraVars.h:58
double vars[kNumPIDVars]
Definition: DecisionTree.h:30
virtual void beginRun(art::Run &run)
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
assert(nhit_max >=nhit_nbins)
Attach LEM-specific info to the base PID object.
Definition: PIDDetails.h:20
float gJunk
Float_t e
Definition: plot.C:35
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
virtual void beginJob()
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
double fBestEnrichPotential
Definition: PIDExtraVars.h:70
bool failedToGet() const
Definition: Handle.h:190
enum BeamMode string