17 #include "Utilities/AssociationUtil.h" 19 #include "Utilities/func/MathUtil.h" 23 #include "TMVA/Reader.h" 99 double&
mean,
double& fit)
const;
128 produces<std::vector<lem::PIDDetails>>();
129 produces<art::Assns<lem::PIDDetails, rb::Cluster>>();
135 if(pset.
get<
bool>(
"FillTMVAs")){
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");
156 fReader->BookMVA(
"BDTG_CC", tmvaPath+
"TMVAClassification_BDTG_CC.weights.xml");
157 fReader->BookMVA(
"MLP_CC", tmvaPath+
"TMVAClassification_MLP_CC.weights.xml");
172 if(pset.
get<
bool>(
"FillDecTree")){
178 if(pset.
get<
bool>(
"Flatten")){
179 TFile* fFlat =
new TFile((libPath+
"flatten.root").c_str());
180 assert(!fFlat->IsZombie());
228 std::cout <<
"Reweighted to sig, bkg, enrich = " 237 double&
mean,
double& fit)
const 239 std::vector<double> xs, ys, ws;
243 const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
245 for(
unsigned int i = 0;
i < matches.size(); ++
i){
252 if(weightScheme ==
kExp ||
264 if(!match.
isSig) w = 0;
268 if(match.
isSig) w = 0;
272 if(!match.
isSig || match.
y > 0.50) w = 0;
275 if(!match.
isSig || match.
y > 0.90) w = 0;
284 case kY: var = match.
y;
break;
328 std::unique_ptr<std::vector<lem::PIDDetails> > pidcol(
new std::vector<lem::PIDDetails>);
337 assert(inputs->size() == matchlists->size());
341 const unsigned int matchMax = matchlists->size();
342 for(
unsigned int matchIdx = 0; matchIdx < matchMax; ++matchIdx){
343 MatchList matchlist = (*matchlists)[matchIdx];
346 std::vector<art::Ptr<rb::Cluster> > slices = fmt.at(matchIdx);
350 assert(slices.empty() || slices.size() == 1);
352 const double totalRecoGeV = (*inputs)[matchIdx].totalGeV;
354 std::vector<MatchSummary>
matches, enrichedMatches;
355 for(
unsigned int i = 0;
i < matchlist.
matches.size(); ++
i){
357 enrichedMatches.push_back(matchlist.
matches[
i]);
360 matches.push_back(matchlist.
matches[
i]);
366 pidcol->push_back(pid);
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;
377 for(
unsigned int n = 0;
n < enrichedMatches.size(); ++
n){
384 std::vector<double> xs, ys, ws, ws2;
388 double meanESigExp = 0;
389 double meanEBkgExp = 0;
392 double totSigExp = 0;
393 double totBkgExp = 0;
396 double invE = 0, invEDenom = 0;
397 double expE = 0, expEDenom = 0;
398 double sigE = 0, sigEDenom = 0;
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;
404 const double invEVal = w/(matches[
i].potential+1
e-5);
405 if(matches[
i].isSig) invE += invEVal;
else invE -= invEVal;
406 invEDenom += invEVal;
408 const double expEVal = w*
exp(-matches[
i].potential);
409 if(matches[
i].isSig) expE += expEVal;
else expE -= expEVal;
410 expEDenom += expEVal;
412 const double sigEVal = w/(1+matches[
i].potential);
413 if(matches[
i].isSig) sigE += sigEVal;
else sigE -= sigEVal;
414 sigEDenom += sigEVal;
416 xs.push_back(matches[
i].potential);
417 ys.push_back(matches[
i].isSig);
423 ws2.push_back(w*expVal);
425 if(matches[
i].isSig){
427 totSigExp += expVal*
w;
428 meanESig += w*matches[
i].potential;
429 meanESigExp += expVal*w*matches[
i].potential;
433 totBkgExp += expVal*
w;
434 meanEBkg += w*matches[
i].potential;
435 meanEBkgExp += expVal*w*matches[
i].potential;
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];
454 const double m = num/denom;
457 for(
unsigned int i = 0;
i < xs.size(); ++
i){
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];
472 const double m = num/denom;
475 for(
unsigned int i = 0;
i < xs.size(); ++
i){
485 if(totSig) meanESig /= totSig;
486 if(totBkg) meanEBkg /= totBkg;
487 if(totSigExp) meanESigExp /= totSigExp;
488 if(totBkgExp) meanEBkgExp /= totBkgExp;
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]);
529 std::sort(enrichedMatches.begin(), enrichedMatches.end());
530 enrichedMatches.resize(emOrig);
563 for(
unsigned int n = 0;
n < matches.size(); ++
n){
564 if(matches[
n].isSig){
569 for(
unsigned int n = 0;
n < matches.size(); ++
n){
570 if(!matches[
n].isSig){
575 for(
unsigned int n = 0;
n < enrichedMatches.size(); ++
n){
576 if(!enrichedMatches[
n].isSig){
639 decevt.
vars[4] = 2*totalRecoGeV;
641 decevt.
vars[4] = totalRecoGeV;
661 pidcol->push_back(pid);
668 evt.
put(std::move(pidcol));
669 evt.
put(std::move(assns));
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
virtual void reconfigure(const fhicl::ParameterSet &pset)
const unsigned int kMaxNumMatches
static Forest * FromFile(const std::string &fname)
Load PID from a file.
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
std::vector< double > fExpTable
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
A training or trial event for the decision tree.
Collection of MatchSummary objects.
Core of the LEM match-finding algorithm.
ProductID put(std::unique_ptr< PROD > &&product)
double Classify(const Evt &evt) const
Calculate the PID value for evt.
Calculate final PID variables from match details.
Details of the library LEM matches were made against.
Simplified Match information, suitable for serialization.
virtual void produce(art::Event &evt)
Collection of MatchSummary objects.
"Random forest" of decision trees
T get(std::string const &key) const
fvar< T > exp(const fvar< T > &x)
MakePID(const fhicl::ParameterSet &pset)
void LinFitUnweighted(const std::vector< double > &x, const std::vector< double > &y, double &m, double &c)
Simplified version of LinFit.
std::string fSummarizeLabel
std::string fMatchesLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
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...
assert(nhit_max >=nhit_nbins)
Attach LEM-specific info to the base PID object.
void AvgAndFit(const std::vector< MatchSummary > &matches, WeightScheme weightScheme, MatchVar matchVar, double &mean, double &fit) const