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  {
126  reconfigure(pset);
127 
128  produces<std::vector<lem::PIDDetails>>();
129  produces<art::Assns<lem::PIDDetails, rb::Cluster>>();
130 
131  const std::string path = pset.get<std::string>("LibraryPath");
132  const std::string libPath = util::EnvExpansion(path) + "/";
133 
134  fReader = 0;
135  if(pset.get<bool>("FillTMVAs")){
136  fReader = new TMVA::Reader;
137 
138  fReader->AddVariable("pidExtrap", &gTMVAVar[0]); // vars.fPIDFit
139  fReader->AddVariable("yexp", &gTMVAVar[1]); // vars.fMeanYExp
140  fReader->AddVariable("qfracexp", &gTMVAVar[2]); // vars.fMeanQFracExp
141  fReader->AddVariable("ediff", &gTMVAVar[3]); // vars.fEnergyDiffExp
142  fReader->AddVariable("calE", &gTMVAVar[4]); // 2*totalRecoGeV
143  fReader->AddVariable("enrichExp", &gTMVAVar[5]); // vars.fEnrichFracExp
144 
145  fReader->AddSpectator("isSig", &gJunk);
146  fReader->AddSpectator("myW", &gJunk);
147 
148  const std::string tmvaPath = libPath+"weights/";
149 
150  fReader->BookMVA("BDTG", tmvaPath+"TMVAClassification_BDTG.weights.xml");
151  fReader->BookMVA("BDT", tmvaPath+"TMVAClassification_BDT.weights.xml");
152  fReader->BookMVA("BDTD", tmvaPath+"TMVAClassification_BDTD.weights.xml");
153  fReader->BookMVA("BDTMitFisher", tmvaPath+"TMVAClassification_BDTMitFisher.weights.xml");
154  fReader->BookMVA("MLP", tmvaPath+"TMVAClassification_MLP.weights.xml");
155 
156  fReader->BookMVA("BDTG_CC", tmvaPath+"TMVAClassification_BDTG_CC.weights.xml");
157  fReader->BookMVA("MLP_CC", tmvaPath+"TMVAClassification_MLP_CC.weights.xml");
158  }
159 
160  // fANN = fANNEnrich = 0;
161  // if(pset.get<bool>("FillANNs")){
162  // fANN = fann_create_from_file((libPath+"fann.net").c_str());
163  // assert(fANN);
164  // fann_print_connections(fANN);
165 
166  // fANNEnrich = fann_create_from_file((libPath+"fann_enrich.net").c_str());
167  // assert(fANNEnrich);
168  // fann_print_connections(fANNEnrich);
169  // }
170 
171  fForest = 0;
172  if(pset.get<bool>("FillDecTree")){
173  std::cout << "Loading LEM dectree from " << libPath << std::endl;
174  fForest = dec::Forest::FromFile(libPath+"dectree.root");
175  }
176 
177  fFlatHist = 0;
178  if(pset.get<bool>("Flatten")){
179  TFile* fFlat = new TFile((libPath+"flatten.root").c_str());
180  assert(!fFlat->IsZombie());
181  fFlatHist = (TH1*)fFlat->Get("conv");
182  assert(fFlatHist);
183  }
184  }
185 
186  //......................................................................
188  {
189  delete fReader;
190  delete fForest;
191  }
192 
193  //......................................................................
195  {
196  fScaleCalE = pset.get<bool>("ScaleCalE");
197 
198  fTrimMatchLists = pset.get<bool>("TrimMatchLists");
199  fUseFitInput = pset.get<bool>("UseFitInput");
200  fPositionalWeight = pset.get<bool>("PositionalWeight");
201  fMatchesLabel = pset.get<std::string>("MatchesLabel");
202  fSummarizeLabel = pset.get<std::string>("SummarizeLabel");
203  fExpLength = pset.get<double>("ExpLength");
204  fExpPow = pset.get<double>("ExpPow");
205  }
206 
207  //......................................................................
209  {
210  fExpTable.resize(kMaxNumMatches);
211  for(unsigned int i = 0; i < kMaxNumMatches; ++i){
213  }
214  }
215 
216  //......................................................................
218  {
220  run.getByLabel(fMatchesLabel, libsum);
221 
222  assert(!libsum.failedToGet());
223 
224  fNCScale = libsum->nTrueNC/double(libsum->nTrueNC+libsum->nSwapNC);
225  fSigScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totSig;
226  fEnrichScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totEnrich;
227 
228  std::cout << "Reweighted to sig, bkg, enrich = "
229  << libsum->totSig*fSigScale << ", "
230  << libsum->totBkg+(fNCScale-1)*libsum->nTrueNC-libsum->nSwapNC << ", "
231  << libsum->totEnrich*fEnrichScale;
232  }
233 
234  //......................................................................
235  void MakePID::AvgAndFit(const std::vector<MatchSummary>& matches,
236  WeightScheme weightScheme, MatchVar matchVar,
237  double& mean, double& fit) const
238  {
239  std::vector<double> xs, ys, ws;
240  double tot = 0;
241  mean = 0;
242 
243  const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
244 
245  for(unsigned int i = 0; i < matches.size(); ++i){
246  const MatchSummary& match = matches[i];
247 
248  double w = match.weight;
249 
250  if(weightScheme == kUnweighted) w = 1;
251 
252  if(weightScheme == kExp ||
253  weightScheme == kSigOnlyExp ||
254  weightScheme == kBkgOnlyExp){
255  if(fPositionalWeight){
256  w *= fExpTable[i];
257  }
258  else{
259  w *= exp(-pow(match.potential/maxPotential, fExpPow)/fExpLength);
260  }
261  }
262 
263  if(weightScheme == kSigOnly || weightScheme == kSigOnlyExp){
264  if(!match.isSig) w = 0;
265  }
266 
267  if(weightScheme == kBkgOnly || weightScheme == kBkgOnlyExp){
268  if(match.isSig) w = 0;
269  }
270 
271  if(weightScheme == kSig50Only){
272  if(!match.isSig || match.y > 0.50) w = 0;
273  }
274  if(weightScheme == kSig90Only){
275  if(!match.isSig || match.y > 0.90) w = 0;
276  }
277  if(weightScheme == kPi0Only) w = (match.photL0 > 0 || match.photL1 > 0 || match.photE0 > 0 || match.photE1 > 0);
278 
279  double var = 0;
280  switch(matchVar){
281  case kSigFrac: var = int(match.isSig); break;
282  case kBkgFrac: var = 1-int(match.isSig); break;
283  case kNumuFrac: var = int(match.ccnc == simb::kCC && abs(match.pdg) == 14); break;
284  case kY: var = match.y; break;
285  case kQFrac: var = match.qFrac; break;
286  case kSigFrac50: var = int(match.isSig && match.y < .50); break;
287  case kSigFrac90: var = int(match.isSig && match.y < .90); break;
288  case kL0: var = match.photL0; break;
289  case kL1: var = match.photL1; break;
290  case kTrueEVis: var = match.trueEVis; break;
291  case kQESigFrac: var = (match.isSig && match.mode == simb::kQE); break;
292  case kResSigFrac: var = (match.isSig && match.mode == simb::kRes); break;
293  case kDISSigFrac: var = (match.isSig && match.mode == simb::kDIS); break;
294  case kQENCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kQE); break;
295  case kResNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kRes); break;
296  case kDISNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kDIS); break;
297  case kQENumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kQE); break;
298  case kResNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kRes); break;
299  case kDISNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kDIS); break;
300  }
301 
302  xs.push_back(match.potential);
303  ys.push_back(var);
304  ws.push_back(w);
305 
306  tot += w;
307  mean += w*var;
308  } // end for i
309 
310  mean /= tot;
311 
312  if(ws.size() < 2){
313  fit = -1;
314  }
315  else{
316  double m, c;
317  if(weightScheme == kUnweighted)
318  util::LinFitUnweighted(xs, ys, m, c);
319  else
320  util::LinFit(xs, ys, ws, m, c);
321  fit = m*xs[0]+c;
322  }
323  }
324 
325  //......................................................................
327  {
328  std::unique_ptr<std::vector<lem::PIDDetails> > pidcol(new std::vector<lem::PIDDetails>);
329  std::unique_ptr<art::Assns<lem::PIDDetails, rb::Cluster>> assns(new art::Assns<lem::PIDDetails, rb::Cluster>);
330 
332  evt.getByLabel(fMatchesLabel, matchlists);
333 
335  evt.getByLabel(fSummarizeLabel, inputs);
336 
337  assert(inputs->size() == matchlists->size());
338 
340 
341  const unsigned int matchMax = matchlists->size();
342  for(unsigned int matchIdx = 0; matchIdx < matchMax; ++matchIdx){
343  MatchList matchlist = (*matchlists)[matchIdx];
344 
345  // Get all the slices associated with this match list
346  std::vector<art::Ptr<rb::Cluster> > slices = fmt.at(matchIdx);
347  // Associated slices can be empty when we're running over lemsum files,
348  // where the slices have been stripped out. In that case, don't make
349  // associations, LEMAssociator will sort it out later.
350  assert(slices.empty() || slices.size() == 1);
351 
352  const double totalRecoGeV = (*inputs)[matchIdx].totalGeV;
353 
354  std::vector<MatchSummary> matches, enrichedMatches;
355  for(unsigned int i = 0; i < matchlist.matches.size(); ++i){
356  if(matchlist.matches[i].enrich){
357  enrichedMatches.push_back(matchlist.matches[i]);
358  }
359  else{
360  matches.push_back(matchlist.matches[i]);
361  }
362  }
363 
364  if(matches.empty()){
365  lem::PIDDetails pid(-1);
366  pidcol->push_back(pid);
367  if(!slices.empty())
368  util::CreateAssn(*this, evt, *pidcol, slices[0], *assns);
369  continue;
370  }
371 
372  // Equalize stats
373  for(unsigned int n = 0; n < matches.size(); ++n){
374  if(matches[n].isSig) matches[n].weight *= fSigScale;
375  if(matches[n].ccnc == 1) matches[n].weight *= fNCScale;
376  }
377  for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
378  enrichedMatches[n].weight *= fEnrichScale;
379  }
380 
381 
383 
384  std::vector<double> xs, ys, ws, ws2;
385 
386  double meanESig = 0;
387  double meanEBkg = 0;
388  double meanESigExp = 0;
389  double meanEBkgExp = 0;
390  double totSig = 0;
391  double totBkg = 0;
392  double totSigExp = 0;
393  double totBkgExp = 0;
394  double totW = 0;
395 
396  double invE = 0, invEDenom = 0;
397  double expE = 0, expEDenom = 0;
398  double sigE = 0, sigEDenom = 0;
399 
400  const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
401  for(unsigned int i = 0; i < matches.size(); ++i){
402  const double w = matches[i].weight;
403 
404  const double invEVal = w/(matches[i].potential+1e-5);
405  if(matches[i].isSig) invE += invEVal; else invE -= invEVal;
406  invEDenom += invEVal;
407 
408  const double expEVal = w*exp(-matches[i].potential);
409  if(matches[i].isSig) expE += expEVal; else expE -= expEVal;
410  expEDenom += expEVal;
411 
412  const double sigEVal = w/(1+matches[i].potential);
413  if(matches[i].isSig) sigE += sigEVal; else sigE -= sigEVal;
414  sigEDenom += sigEVal;
415 
416  xs.push_back(matches[i].potential);
417  ys.push_back(matches[i].isSig);
418  ws.push_back(w);
419  double expVal = fExpTable[i];
420  if(!fPositionalWeight){
421  expVal = exp(-pow(matches[i].potential/maxPotential, fExpPow)/fExpLength);
422  }
423  ws2.push_back(w*expVal);
424 
425  if(matches[i].isSig){
426  totSig += w;
427  totSigExp += expVal*w;
428  meanESig += w*matches[i].potential;
429  meanESigExp += expVal*w*matches[i].potential;
430  }
431  else{
432  totBkg += w;
433  totBkgExp += expVal*w;
434  meanEBkg += w*matches[i].potential;
435  meanEBkgExp += expVal*w*matches[i].potential;
436  }
437 
438  totW += w;
439  } // end for i
440 
441  if(ws.size() < 2){
442  vars.fChiSig = 0;
443  vars.fChiBkg = 0;
444  vars.fChiSigExp = 0;
445  vars.fChiBkgExp = 0;
446  }
447  else{
448  for(int c = 0; c <= 1; ++c){
449  double num = 0, denom = 0;
450  for(unsigned int i = 0; i < xs.size(); ++i){
451  num += ws[i]*(xs[i]*ys[i]-c*xs[i]);
452  denom += xs[i]*xs[i]*ws[i];
453  }
454  const double m = num/denom;
455  double chi = 0;
456  double totW = 0;
457  for(unsigned int i = 0; i < xs.size(); ++i){
458  chi += ws[i]*util::sqr(ys[i]-m*xs[i]-c);
459  totW += ws[i];
460  }
461  chi /= totW;
462  if(c == 0) vars.fChiBkg = chi;
463  if(c == 1) vars.fChiSig = chi;
464  } // end for c
465 
466  for(int c = 0; c <= 1; ++c){
467  double num = 0, denom = 0;
468  for(unsigned int i = 0; i < xs.size(); ++i){
469  num += ws2[i]*(xs[i]*ys[i]-c*xs[i]);
470  denom += xs[i]*xs[i]*ws2[i];
471  }
472  const double m = num/denom;
473  double chi = 0;
474  double totW = 0;
475  for(unsigned int i = 0; i < xs.size(); ++i){
476  chi += ws2[i]*util::sqr(ys[i]-m*xs[i]-c);
477  totW += ws2[i];
478  }
479  chi /= totW;
480  if(c == 0) vars.fChiBkgExp = chi;
481  if(c == 1) vars.fChiSigExp = chi;
482  } // end for c
483  }
484 
485  if(totSig) meanESig /= totSig;
486  if(totBkg) meanEBkg /= totBkg;
487  if(totSigExp) meanESigExp /= totSigExp;
488  if(totBkgExp) meanEBkgExp /= totBkgExp;
489 
490 
491  double junk;
492 
493  AvgAndFit(matches, kSimple, kSigFrac, junk, vars.fPIDFit);
494  AvgAndFit(matches, kUnweighted, kSigFrac, junk, vars.fPIDFitUnweight);
495  AvgAndFit(matches, kExp, kSigFrac, vars.fPIDExp, vars.fPIDExpFit);
496 
497  AvgAndFit(matches, kSimple, kY, vars.fMeanY, vars.fMeanYFit);
498  AvgAndFit(matches, kExp, kY, vars.fMeanYExp, junk);
499 
500  AvgAndFit(matches, kSigOnly, kY, vars.fMeanYSig, vars.fMeanYSigFit);
501  AvgAndFit(matches, kSigOnlyExp, kY, vars.fMeanYSigExp, junk);
502 
503  AvgAndFit(matches, kSimple, kQFrac, vars.fMeanQFrac, vars.fMeanQFracFit);
504  AvgAndFit(matches, kExp, kQFrac, vars.fMeanQFracExp, junk);
505 
506  AvgAndFit(matches, kSigOnly, kQFrac, vars.fMeanQFracSig, vars.fMeanQFracSigFit);
507  AvgAndFit(matches, kSigOnlyExp, kQFrac, vars.fMeanQFracSigExp, junk);
508 
509  AvgAndFit(matches, kBkgOnly, kQFrac, vars.fMeanQFracBkg, vars.fMeanQFracBkgFit);
510  AvgAndFit(matches, kBkgOnlyExp, kQFrac, vars.fMeanQFracBkgExp, junk);
511 
512  AvgAndFit(matches, kExp, kTrueEVis, vars.fMeanMatchEExp, junk);
513 
514  AvgAndFit(matches, kExp, kQESigFrac, vars.fQESigFracExp, junk);
515  AvgAndFit(matches, kExp, kResSigFrac, vars.fResSigFracExp, junk);
516  AvgAndFit(matches, kExp, kDISSigFrac, vars.fDISSigFracExp, junk);
517  AvgAndFit(matches, kExp, kQENCFrac, vars.fQENCFracExp, junk);
518  AvgAndFit(matches, kExp, kResNCFrac, vars.fResNCFracExp, junk);
519  AvgAndFit(matches, kExp, kDISNCFrac, vars.fDISNCFracExp, junk);
520  AvgAndFit(matches, kExp, kQENumuFrac, vars.fQENumuFracExp, junk);
521  AvgAndFit(matches, kExp, kResNumuFrac, vars.fResNumuFracExp, junk);
522  AvgAndFit(matches, kExp, kDISNumuFrac, vars.fDISNumuFracExp, junk);
523 
524  // Mix in the usual signal matches
525  const int emOrig = enrichedMatches.size();
526  for(unsigned int n = 0; n < matches.size(); ++n){
527  if(matches[n].isSig) enrichedMatches.push_back(matches[n]);
528  }
529  std::sort(enrichedMatches.begin(), enrichedMatches.end());
530  enrichedMatches.resize(emOrig);
531 
532  AvgAndFit(enrichedMatches, kSimple, kBkgFrac, vars.fEnrichFrac, vars.fEnrichFracFit);
533  AvgAndFit(enrichedMatches, kExp, kBkgFrac, vars.fEnrichFracExp, junk);
534 
535  AvgAndFit(enrichedMatches, kBkgOnly, kQFrac, vars.fEnrichQFrac, vars.fEnrichQFracFit);
536  AvgAndFit(enrichedMatches, kBkgOnlyExp, kQFrac, vars.fEnrichQFracExp, junk);
537 
538  AvgAndFit(enrichedMatches, kSimple, kL0, junk, vars.fMeanL0Fit);
539  AvgAndFit(enrichedMatches, kSimple, kL1, junk, vars.fMeanL1Fit);
540  AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0Exp, junk);
541  AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1Exp, junk);
542 
543  AvgAndFit(enrichedMatches, kPi0Only, kL0, junk, vars.fMeanL0PiOnlyFit);
544  AvgAndFit(enrichedMatches, kPi0Only, kL1, junk, vars.fMeanL1PiOnlyFit);
545  AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0PiOnlyExp, junk);
546  AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1PiOnlyExp, junk);
547 
548  AvgAndFit(matches, kSimple, kSigFrac50, junk, vars.fPIDy50Fit);
549  AvgAndFit(matches, kSimple, kSigFrac90, junk, vars.fPIDy90Fit);
550 
551  AvgAndFit(matches, kSig50Only, kQFrac, junk, vars.fMeanQFracSigY50Fit);
552  AvgAndFit(matches, kSig90Only, kQFrac, junk, vars.fMeanQFracSigY90Fit);
553 
554  AvgAndFit(matches, kSig50Only, kY, junk, vars.fMeanYSigY50Fit);
555  AvgAndFit(matches, kSig90Only, kY, junk, vars.fMeanYSigY90Fit);
556 
557  AvgAndFit(matches, kSimple, kNumuFrac, junk, vars.fPIDFracNumuFit);
558 
559  vars.fEnergyDiff = meanESig-meanEBkg;
560  vars.fEnergyDiffExp = meanESigExp-meanEBkgExp;
561 
563  for(unsigned int n = 0; n < matches.size(); ++n){
564  if(matches[n].isSig){
565  vars.fBestSigPotential = matches[n].potential;
566  break;
567  }
568  }
569  for(unsigned int n = 0; n < matches.size(); ++n){
570  if(!matches[n].isSig){
571  vars.fBestBkgPotential = matches[n].potential;
572  break;
573  }
574  }
575  for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
576  if(!enrichedMatches[n].isSig){
577  vars.fBestEnrichPotential = enrichedMatches[n].potential;
578  break;
579  }
580  }
581 
582  vars.fAvgInvE = invE/invEDenom;
583  vars.fAvgExpE = expE/expEDenom;
584  vars.fAvgSigE = sigE/sigEDenom;
585 
586  // if(fANN || fANNEnrich){
587  // fann_type input[6];
588  // input[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
589  // input[1] = vars.fMeanYExp;
590  // input[2] = vars.fMeanQFracExp;
591  // input[3] = vars.fEnergyDiffExp;
592 
593  // if(fANN){
594  // fann_type* calc_out = fann_run(fANN, input);
595  // vars.fFannId = calc_out[0];
596  // }
597 
598  // input[4] = 2*totalRecoGeV;
599  // input[5] = vars.fEnrichFracExp;
600 
601  // if(fANNEnrich){
602  // fann_type* calc_out = fann_run(fANNEnrich, input);
603  // vars.fFannIdEnrich = calc_out[0];
604  // }
605  // }
606 
607  if(fReader){
608  gTMVAVar[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
609  gTMVAVar[1] = vars.fMeanYExp;
610  gTMVAVar[2] = vars.fMeanQFracExp;
611  gTMVAVar[3] = vars.fEnergyDiffExp;
612  if(fScaleCalE)
613  gTMVAVar[4] = 2*totalRecoGeV;
614  else
615  gTMVAVar[4] = totalRecoGeV;
616  gTMVAVar[5] = vars.fEnrichFracExp;
617 
618  vars.fTMVABDTG = fReader->EvaluateMVA("BDTG");
619  vars.fTMVABDT = fReader->EvaluateMVA("BDT");
620  vars.fTMVABDTD = fReader->EvaluateMVA("BDTD");
621  vars.fTMVABDTMF = fReader->EvaluateMVA("BDTMitFisher");
622  vars.fTMVAMLP = fReader->EvaluateMVA("MLP");
623 
624  vars.fTMVABDTG_CC = fReader->EvaluateMVA("BDTG_CC");
625  vars.fTMVAMLP_CC = fReader->EvaluateMVA("MLP_CC");
626  }
627 
628  if(fForest){
629  if(vars.fPIDExp == -1){
630  vars.fDecTreePID = -1;
631  }
632  else{
633  dec::Evt decevt;
634  decevt.vars[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
635  decevt.vars[1] = vars.fMeanYExp;
636  decevt.vars[2] = vars.fMeanQFracExp;
637  decevt.vars[3] = vars.fEnergyDiffExp;
638  if(fScaleCalE)
639  decevt.vars[4] = 2*totalRecoGeV;
640  else
641  decevt.vars[4] = totalRecoGeV;
642  decevt.vars[5] = vars.fEnrichFracExp;
643  vars.fDecTreePID = fForest->Classify(decevt);
644  }
645  }
646 
647  double defaultPID = vars.fDecTreePID;
648 
649  // Flattened version of the Decision Tree variable
650  if(fFlatHist) defaultPID = fFlatHist->GetBinContent(fFlatHist->FindBin(vars.fDecTreePID));
651 
652 
653  lem::PIDDetails pid(defaultPID);
654 
655  if(fTrimMatchLists) matchlist.matches.resize(1);
656 
657  // TODO: sample size equalization doesn't survive this
658  pid.matchList = matchlist;
659  pid.vars = vars;
660 
661  pidcol->push_back(pid);
662 
663  if(!slices.empty()){
664  util::CreateAssn(*this, evt, *pidcol, slices[0], *assns);
665  }
666  } // end for matchIdx
667 
668  evt.put(std::move(pidcol));
669  evt.put(std::move(assns));
670  }
671 
673 
674 } // end namespace lem
675 ////////////////////////////////////////////////////////////////////////
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
constexpr T pow(T x)
Definition: pow.h:75
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
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:31
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.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
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)
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
"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
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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
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:196
enum BeamMode string