From 5a3d7572f982758c6be7b95bc4111841d5588ec7 Mon Sep 17 00:00:00 2001 From: Connor Pecar <39525107+cpecar@users.noreply.github.com> Date: Tue, 27 Dec 2022 12:39:26 -0500 Subject: [PATCH] HFS simple tree + neural network predictions using pybind (#214) --- .gitignore | 1 + config.mk | 10 ++- macro/analysis_testml.C | 48 ++++++++++++++ src/Analysis.cxx | 28 +++++---- src/Analysis.h | 5 +- src/AnalysisEcce.cxx | 69 ++++++++++++++++----- src/HFSTree.cxx | 31 ++++++++++ src/HFSTree.h | 47 ++++++++++++++ src/Kinematics.cxx | 134 ++++++++++++++++++++++++++++++++++++++-- src/Kinematics.h | 41 ++++++++++-- src/LinkDef.h | 1 + testEFlowimport.py | 20 ++++++ 12 files changed, 392 insertions(+), 43 deletions(-) create mode 100644 macro/analysis_testml.C create mode 100644 src/HFSTree.cxx create mode 100644 src/HFSTree.h create mode 100644 testEFlowimport.py diff --git a/.gitignore b/.gitignore index 6955a74e..c5a37c64 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ # main ignores .*.sw* +*~ *.hepmc* *.root *.pcm diff --git a/config.mk b/config.mk index 80847587..40f9f6b3 100644 --- a/config.mk +++ b/config.mk @@ -1,9 +1,10 @@ # compiler and flags CXX = g++ -FLAGS = -g -Wno-deprecated -fPIC -fno-inline -Wno-write-strings +FLAGS = -Wno-deprecated -fPIC FLAGS += -fmax-errors=3 -# extra flags -#FLAGS += -O0 +# FLAGS += -fvisibility=hidden # FIXME: required by pybind, but causes unresolved symbols in cling... +# FLAGS += -g # build with debugging symbols +# FLAGS += -O0 # extra flags for Mac OS UNAME := $(shell uname) @@ -19,6 +20,9 @@ LIBS = $(shell root-config --glibs) # Data Model (PODIO + EDM4hep + EDM4eic) LIBS += -L/usr/local/lib -lpodio -lpodioRootIO -ledm4hep -ledm4eic +# PYTHON (for pybind) +LIBS += -lpython3.10 + # Miscellaneous LIBS += -lfmt diff --git a/macro/analysis_testml.C b/macro/analysis_testml.C new file mode 100644 index 00000000..09ca14c0 --- /dev/null +++ b/macro/analysis_testml.C @@ -0,0 +1,48 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Connor Pecar + +R__LOAD_LIBRARY(EpicAnalysis) +#include +namespace py = pybind11; +// using ML prediction for vecQ, and writing out tree of HFS four-vectors +// requires pybind includes above +void analysis_testml( + TString configFile="datarec/in.config", /* delphes tree(s) */ + TString outfilePrefix="resolutions" /* output filename prefix*/ +) { + // object needed at start of script using pybind11 + py::scoped_interpreter guard{}; + //outfilePrefix+="_DA"; + // setup analysis ======================================== + AnalysisEcce *A = new AnalysisEcce( + configFile, + outfilePrefix + ); + + A->maxEvents = 10000; // use this to limit the number of events + A->writeSimpleTree = true; // write SimpleTree (for one bin) + A->writeHFSTree = true; // write HFSTree (for one bin) + A->SetReconMethod("ML"); // set reconstruction method + A->AddFinalState("pipTrack"); // pion final state + //A->AddFinalState("KpTrack"); // kaon final state + + + // define cuts ==================================== + A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV + //A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0 + //A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit) + + + // set binning scheme ==================================== + // z ranges + //A->AddBinScheme("z"); + //A->BinScheme("z")->BuildBin("Min",0.2); // needed? + + // y minima + A->AddBinScheme("y"); + A->BinScheme("y")->BuildBin("Full"); // a bin with no y-cut + + // perform the analysis ================================== + A->Execute(); + +}; diff --git a/src/Analysis.cxx b/src/Analysis.cxx index c296d04f..fa2ee695 100644 --- a/src/Analysis.cxx +++ b/src/Analysis.cxx @@ -92,16 +92,18 @@ Analysis::Analysis( // common settings defaults // - these settings can be set at the macro level - verbose = false; - writeSimpleTree = false; + verbose = false; + writeSimpleTree = false; + writeHFSTree = false; writeParticleTree = false; - maxEvents = 0; - useBreitJets = false; - errorCntMax = 1000; - jetAlg = 0; // Default to kT Algorithm - jetRad = 0.8; - jetMin = 1.0; // Minimum Jet pT - jetMatchDR = 0.5; // Delta R between true and reco jet to be considered matched + + maxEvents = 0; + useBreitJets = false; + errorCntMax = 1000; + jetAlg = 0; // Default to kT Algorithm + jetRad = 0.8; + jetMin = 1.0; // Minimum Jet pT + jetMatchDR = 0.5; // Delta R between true and reco jet to be considered matched weightInclusive = new WeightsUniform(); weightTrack = new WeightsUniform(); @@ -310,8 +312,9 @@ void Analysis::Prepare() { // instantiate shared objects kin = new Kinematics(eleBeamEn,ionBeamEn,crossingAngle); kinTrue = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle); - ST = new SimpleTree("tree",kin,kinTrue); - PT = new ParticleTree("ptree"); + ST = new SimpleTree("tree",kin,kinTrue); + HFST = new HFSTree("hfstree",kin,kinTrue); + PT = new ParticleTree("ptree"); // if including jets, define a `jet` final state #ifndef EXCLUDE_DELPHES @@ -660,7 +663,8 @@ void Analysis::Finish() { cout << sep << endl; cout << "writing ROOT file..." << endl; outFile->cd(); - if(writeSimpleTree) ST->WriteTree(); + if(writeSimpleTree) ST->WriteTree(); + if(writeHFSTree) HFST->WriteTree(); if(writeParticleTree) PT->WriteTree(); HD->Payload([this](Histos *H){ H->WriteHists(outFile); }); HD->ExecuteAndClearOps(); HD->Payload([this](Histos *H){ H->Write(); }); HD->ExecuteAndClearOps(); diff --git a/src/Analysis.h b/src/Analysis.h index 0d59a1c2..e98b1832 100644 --- a/src/Analysis.h +++ b/src/Analysis.h @@ -33,6 +33,7 @@ #include "HistosDAG.h" #include "Kinematics.h" #include "SimpleTree.h" +#include "HFSTree.h" #include "ParticleTree.h" #include "Weights.h" #include "CommonConstants.h" @@ -60,7 +61,8 @@ class Analysis : public TNamed // common settings Bool_t verbose; // if true, print a lot more information - Bool_t writeSimpleTree; // if true, write SimpleTree (not binned) + Bool_t writeSimpleTree; // if true, write SimpleTree (not binned) + Bool_t writeHFSTree; // if true, write HFSTree (not binned) Bool_t writeParticleTree; // if true, write ParticleTree (not binned) Long64_t maxEvents; /* default=0, which runs all events; * if > 0, run a maximum number of `maxEvents` events (useful for quick tests) @@ -124,6 +126,7 @@ class Analysis : public TNamed // shared objects SimpleTree *ST; + HFSTree *HFST; ParticleTree *PT; Kinematics *kin, *kinTrue; HistosDAG *HD; diff --git a/src/AnalysisEcce.cxx b/src/AnalysisEcce.cxx index 71adbaa3..ae78b325 100644 --- a/src/AnalysisEcce.cxx +++ b/src/AnalysisEcce.cxx @@ -101,14 +101,11 @@ void AnalysisEcce::Execute() tr.SetEntriesRange(1,maxEvents); do { if(tr.GetCurrentEntry()%10000==0) cout << tr.GetCurrentEntry() << " events..." << endl; - + // resets kin->ResetHFS(); kinTrue->ResetHFS(); - - - // a few maps needed to get the associated info between tracks, true particles, etec std::map mcidmap; @@ -140,16 +137,16 @@ void AnalysisEcce::Execute() * - find beam particles */ std::vector mcpart; - double maxP = 0; + double maxPt = 0; int genEleID = -1; bool foundBeamElectron = false; bool foundBeamIon = false; + int genEleBCID = -1; for(int imc=0; imcsecond; } - + + if (imcpart >-1){ px_ = mcpart_psx[imcpart]; py_ = mcpart_psy[imcpart]; pz_ = mcpart_psz[imcpart]; e_ = mcpart_E[imcpart]; p_ = sqrt(pow(mcpart_psx[imcpart],2) + pow(mcpart_psy[imcpart],2) + pow(mcpart_psz[imcpart],2)); + pt_ = sqrt(pow(mcpart_psx[imcpart],2) + pow(mcpart_psy[imcpart],2)); part.mcID = mcpart_ID[imcpart]; } else part.mcID = -1; - + //cout << genStatus_ << " " << pid_ << " " << part.mcID << " " << hepmcp_BCID[imc] << " " << hepmcp_m1[imc] << " " << hepmcp_m2[imc] << " " + // << px_ << " " << py_ << " " << pz_ << " " << e_ << endl; + part.pid = pid_; part.vecPart.SetPxPyPzE(px_, py_, pz_, e_); @@ -193,15 +198,16 @@ void AnalysisEcce::Execute() // identify scattered electron by max momentum if(pid_ == 11) { - if(p_ > maxP) { - maxP = p_; + if(pt_ > maxPt) { + maxPt = pt_; kinTrue->vecElectron.SetPxPyPzE(px_, py_, pz_, e_); genEleID = part.mcID; //mcpart_ID[imc]; + genEleBCID = hepmcp_BCID[imc]; // cout << "\t\t\t found scattered electron " << Form(" %6.2f %6.2f %6.2f %6.2f %6.2f %5.3f %6.2f %6.2f id %3d\n",px_,py_,pz_, sqrt(p_*p_ + mass_*mass_),p_,mass_,hepmcp_E[imc],mcpart_E[imcpart],genEleID); } } } - + else if(genStatus_ == 4) { // beam particles if(pid_ == 11) { // electron beam if(!foundBeamElectron) { @@ -221,11 +227,37 @@ void AnalysisEcce::Execute() } } } // end truth loop + + // looking for radiated photons (mother particle == scattered electron) + // requires that we already know scattered electron index + for(int imc=0; imcvecElectron + FSRmom; + // correcting true scattered electron with FSR + kinTrue->vecElectron.SetPxPyPzE(eleCorr.Px(), eleCorr.Py(), eleCorr.Pz(), eleCorr.E()); + } + } + } + + + } // check beam finding if(!foundBeamElectron || !foundBeamIon) { numNoBeam++; continue; }; - // reconstructed particles loop /* - add reconstructed particle to `recopart` * - find the scattered electron @@ -268,14 +300,15 @@ void AnalysisEcce::Execute() part.vecPart.SetPxPyPzE(reco_px, reco_py, reco_pz, sqrt(reco_p*reco_p + reco_mass*reco_mass)); - // cout << "\t\t\t track " << Form(" %4.2f %4.2f %4.2f true id %4d imc %3d mcid %3d \n",reco_px,reco_py,reco_pz,tracks_trueID[ireco],imc,part.mcID); + //cout << "\t\t\t track " << Form(" %4.2f %4.2f %4.2f true id %4d imc %3d mcid %3d \n",reco_px,reco_py,reco_pz,tracks_trueID[ireco],imc,part.mcID); // add to `recopart` and hadronic final state sums only if there is a matching truth particle - if(part.mcID > 0) { + if(part.mcID > 0 && part.mcID != genEleID) { if(imc>-1) { // cout << "\t\t\t add to hadfs \n" ; recopart.push_back(part); kin->AddToHFS(part.vecPart); + kin->hfspid[kin->nHFS - 1] = pid_; } } @@ -288,6 +321,7 @@ void AnalysisEcce::Execute() } // end reco loop + //add all clusters for hadronic final state int jentryt = tr.GetCurrentEntry(); Long64_t localEntry = chain->LoadTree(jentryt); @@ -386,10 +420,13 @@ void AnalysisEcce::Execute() // calculate DIS kinematics if(!(kin->CalculateDIS(reconMethod))) continue; // reconstructed if(!(kinTrue->CalculateDIS(reconMethod))) continue; // generated (truth) - + // Get the weight for this event's Q2 auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]); + //fill HFS tree for event + if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor); + // fill inclusive histograms, if only `inclusive` is included in output // (otherwise they will be filled in track and jet loops) if(includeOutputSet["inclusive_only"]) { diff --git a/src/HFSTree.cxx b/src/HFSTree.cxx new file mode 100644 index 00000000..d18ab2a3 --- /dev/null +++ b/src/HFSTree.cxx @@ -0,0 +1,31 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Connor Pecar + +#include "HFSTree.h" + +ClassImp(HFSTree) + +// constructor +HFSTree::HFSTree(TString treeName_, Kinematics *K_, Kinematics *Ktrue_) + : treeName(treeName_) + , K(K_) + , Ktrue(Ktrue_) +{ + T = new TTree(treeName,treeName); + T->Branch("vecElectron", &(K->vecElectron)); + T->Branch("vecElectronTrue", &(Ktrue->vecElectron)); + T->Branch("vecEleBeamTrue", &(Ktrue->vecEleBeam)); + T->Branch("vecIonBeamTrue", &(Ktrue->vecIonBeam)); + T->Branch("nHFS", &(K->nHFS), "nHFS/I"); + T->Branch("hfsp4", &(K->hfsp4)); + //T->Branch("hfspy", &(K->hfspy), "hfspy[nHFS]/F"); + //T->Branch("hfspz", &(K->hfspz), "hfspz[nHFS]/F"); + //T->Branch("hfsE", &(K->hfsE), "hfsE[nHFS]/F"); + T->Branch("hfspid", &(K->hfspid), "hfspid[nHFS]/F"); + T->Branch("weight", &(weight), "weight/D"); +}; + +// destructor +HFSTree::~HFSTree() { +}; + diff --git a/src/HFSTree.h b/src/HFSTree.h new file mode 100644 index 00000000..3a71d204 --- /dev/null +++ b/src/HFSTree.h @@ -0,0 +1,47 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Connor Pecar + +/* HFSTree + - Produces a tree containing information needed for kinematic + reconstruction studies: hadronic final state four-momenta + and PID, scattered electron information, and true scattered + electron and beam information. + */ +#ifndef HFSTree_ +#define HFSTree_ + +#include +#include +#include +#include + +// epic-analysis +#include "Kinematics.h" + +// ROOT +#include "TTree.h" + +class HFSTree : public TObject +{ + public: + HFSTree(TString treeName_, Kinematics *K_, Kinematics *Ktrue_); + ~HFSTree(); + + TTree *GetTree() { return T; }; + Kinematics *GetKinematics() { return K; }; + Kinematics *GetKinematicsTrue() { return Ktrue; }; + void FillTree(Double_t w) { weight = w; + T->Fill(); }; + void WriteTree() { T->Write(); }; + + private: + Double_t weight; + TTree *T; + Kinematics *K; + Kinematics *Ktrue; + TString treeName; + + ClassDef(HFSTree,1); +}; + +#endif diff --git a/src/Kinematics.cxx b/src/Kinematics.cxx index 8f540e0f..b85ac0c8 100644 --- a/src/Kinematics.cxx +++ b/src/Kinematics.cxx @@ -6,19 +6,24 @@ */ #include "Kinematics.h" - ClassImp(Kinematics) Kinematics::Kinematics( Double_t enEleBeam, /*GeV*/ Double_t enIonBeam, /*GeV*/ - Double_t crossAng /*mrad*/ + Double_t crossAng /*mrad*/ ) { - + srand(time(NULL)); + // importing from local python script for ML predictions + // requires tensorflow, energyflow packages installed +#ifdef SIDIS_MLPRED + efnpackage = py::module_::import("testEFlowimport"); + pfnimport = efnpackage.attr("eflowPredict"); +#endif // set ion mass IonMass = ProtonMass(); - + // revise crossing angle crossAng *= 1e-3; // mrad -> rad crossAng = -1*TMath::Abs(crossAng); // take -1*abs(crossAng) to enforce the correct sign @@ -196,6 +201,97 @@ void Kinematics::GetQWNu_electronic(){ Nu = vecIonBeam.Dot(vecQ) / IonMass; }; +void Kinematics::GetQWNu_ML(){ + hfsinfo.clear(); + float pidadj = 0; + if(nHFS >= 2){ + std::vector partHold; + for(int i = 0; i < nHFS; i++){ + double pidsgn=(hfspid[i]/abs(hfspid[i])); + if(abs(hfspid[i])==211) pidadj = 0.4*pidsgn; + if(abs(hfspid[i])==22) pidadj = 0.2*pidsgn; + if(abs(hfspid[i])==321) pidadj = 0.6*pidsgn; + if(abs(hfspid[i])==2212) pidadj = 0.8*pidsgn; + if(abs(hfspid[i])==11) pidadj = 1.0*pidsgn; + partHold.push_back(hfseta[i]); + partHold.push_back(hfsphi[i]); + partHold.push_back(hfspx[i]); + partHold.push_back(hfspy[i]); + partHold.push_back(hfspz[i]); + partHold.push_back(hfsE[i]); + partHold.push_back(pidadj); + hfsinfo.push_back(partHold); + partHold.clear(); + } + double Q2ele, Q2DA, Q2JB; + double xele, xDA, xJB; + TLorentzVector vecQEle; + globalinfo.clear(); + this->CalculateDISbyElectron(); + vecQEle.SetPxPyPzE(vecQ.Px(), vecQ.Py(), vecQ.Pz(), vecQ.E()); + Q2ele = Q2; + xele = x; + this->CalculateDISbyDA(); + Q2DA = Q2; + xDA = x; + this->CalculateDISbyJB(); + Q2JB = Q2; + xJB = x; + if( Q2DA > 0 && Q2DA < 1e4){ + globalinfo.push_back(log10(Q2DA)); + } + else{ + globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0))); + } + if( Q2ele > 0 && Q2ele < 1e4){ + globalinfo.push_back(log10(Q2ele)); + } + else{ + globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0))); + } + if( Q2JB > 0 && Q2JB < 1e4){ + globalinfo.push_back(log10(Q2JB)); + } + else{ + globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0))); + } + if(xDA>0 && xDA < 1){ + globalinfo.push_back(-1*log10(xDA)); + } + else{ + globalinfo.push_back(-1*log10( (float) (rand()) / (float) (RAND_MAX/1.0) )); + } + if(xele>0 && xele < 1){ + globalinfo.push_back(-1*log10(xele)); + } + else{ + globalinfo.push_back(-1*log10( (float) (rand()) / (float) (RAND_MAX/1.0) )); + } + if(xJB>0 && xJB < 1){ + globalinfo.push_back(-1*log10(xJB)); + } + else{ + globalinfo.push_back( -1*log10((float) (rand()) / (float) (RAND_MAX/1.0) ) ); + } + globalinfo.push_back(vecQEle.Px()); + globalinfo.push_back(vecQEle.Py()); + globalinfo.push_back(vecQEle.Pz()); + globalinfo.push_back(vecQEle.E()); +#ifdef SIDIS_MLPRED + py::object nnoutput = pfnimport(hfsinfo, globalinfo); + std::vector nnvecq = nnoutput.cast>(); + vecQ.SetPxPyPzE(nnvecq[0],nnvecq[1],nnvecq[2],nnvecq[3]); +#endif + } + else{ + this->CalculateDISbyElectron(); + } + vecW = vecEleBeam + vecIonBeam - vecElectron; + W = vecW.M(); + Nu = vecIonBeam.Dot(vecQ) / IonMass; + +} + // ------------------------------------------------------ @@ -209,7 +305,7 @@ Bool_t Kinematics::CalculateDIS(TString recmethod){ this->TransformToHeadOnFrame(vecEleBeam,HvecEleBeam); this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam); this->TransformToHeadOnFrame(vecElectron,HvecElectron); - + // calculate primary DIS variables, including Q2,x,y,W,nu if (recmethod.CompareTo( "Ele", TString::kIgnoreCase)==0) { this->CalculateDISbyElectron(); } else if(recmethod.CompareTo( "DA", TString::kIgnoreCase)==0) { this->CalculateDISbyDA(); } @@ -217,6 +313,7 @@ Bool_t Kinematics::CalculateDIS(TString recmethod){ else if(recmethod.CompareTo( "Mixed", TString::kIgnoreCase)==0) { this->CalculateDISbyMixed(); } else if(recmethod.CompareTo( "Sigma", TString::kIgnoreCase)==0) { this->CalculateDISbySigma(); } else if(recmethod.CompareTo( "eSigma", TString::kIgnoreCase)==0) { this->CalculateDISbyeSigma(); } + else if(recmethod.CompareTo( "ML", TString::kIgnoreCase)==0) { this->CalculateDISbyML(); } else { cerr << "ERROR: unknown reconstruction method" << endl; return false; @@ -352,6 +449,12 @@ void Kinematics::CalculateDISbyeSigma(){ } }; +void Kinematics::CalculateDISbyML() { + this->GetQWNu_ML(); // set `vecQ`, `vecW`, `W`, `Nu` + Q2 = -1 * vecQ.M2(); + x = Q2 / ( 2 * vecQ.Dot(vecIonBeam) ); + y = vecIonBeam.Dot(vecQ) / vecIonBeam.Dot(vecEleBeam); +}; // calculate hadron kinematics @@ -429,6 +532,13 @@ void Kinematics::ValidateHeadOnFrame() { // add a 4-momentum to the hadronic final state void Kinematics::AddToHFS(TLorentzVector p4_) { TLorentzVector p4 = p4_; + hfspx[nHFS] = p4.Px(); + hfspy[nHFS] = p4.Py(); + hfspz[nHFS] = p4.Pz(); + hfsE[nHFS] = p4.E(); + new(ar[nHFS]) TLorentzVector(p4); + nHFS++; + if(mainFrame==fHeadOn) this->TransformToHeadOnFrame(p4,p4); sigmah += (p4.E() - p4.Pz()); Pxh += p4.Px(); @@ -437,6 +547,12 @@ void Kinematics::AddToHFS(TLorentzVector p4_) { countHadrons++; }; +void Kinematics::AddPion(TLorentzVector p4_){ + TLorentzVector p4 = p4_; + if(mainFrame==fHeadOn) this->TransformToHeadOnFrame(p4,p4); + new(arpi[nPi]) TLorentzVector(p4); + nPi++; +}; // subtract electron from hadronic final state variables void Kinematics::SubtractElectronFromHFS() { @@ -457,6 +573,8 @@ void Kinematics::SubtractElectronFromHFS() { break; } countHadrons--; + ar.Remove( &vecElectron ); + nHFS--; } else { cerr << "ERROR: electron energy is NaN" << endl; // TODO: kill event @@ -465,10 +583,14 @@ void Kinematics::SubtractElectronFromHFS() { // reset some variables for the hadronic final state -void Kinematics::ResetHFS() { +void Kinematics::ResetHFS() { sigmah = Pxh = Pyh = 0; hadronSumVec.SetPxPyPzE(0,0,0,0); countHadrons = 0; + nHFS = 0; + nPi = 0; + ar.Clear(); + arpi.Clear(); }; diff --git a/src/Kinematics.h b/src/Kinematics.h index 33b8b31b..aee5e953 100644 --- a/src/Kinematics.h +++ b/src/Kinematics.h @@ -24,6 +24,7 @@ #include "TLorentzVector.h" #include "TRandom.h" #include "TRandomGen.h" +#include "TClonesArray.h" // Delphes #ifndef EXCLUDE_DELPHES @@ -33,6 +34,12 @@ #include "fastjet/plugins/Centauro/Centauro.hh" #endif #endif +// pybind (for ML models using python packages) +#include +#include +#include +#include +namespace py = pybind11; using std::map; using std::cout; @@ -51,6 +58,7 @@ class Kinematics : public TObject // hadronic final state (HFS) void AddToHFS(TLorentzVector p4_); + void AddPion(TLorentzVector p4_); void SubtractElectronFromHFS(); void ResetHFS(); @@ -104,7 +112,7 @@ class Kinematics : public TObject Double_t pLab,pTlab,phiLab,etaLab,z,pT,qT,mX,xF,phiH,phiS; // hadron Double_t sigmah, Pxh, Pyh; // hadronic final state variables TLorentzVector hadronSumVec; - + // nucleon transverse spin; if you set this externally, // it must be done before calculating `phiS` (before // `CalculateHadronKinematics`) @@ -153,7 +161,23 @@ class Kinematics : public TObject // struck quark information Double_t quarkpT; #endif - + // HFS tree objects + Int_t nHFS; + Int_t nPi; + Double_t hfspx[100]; + Double_t hfspy[100]; + Double_t hfspz[100]; + Double_t hfsE[100]; + Double_t hfseta[100]; + Double_t hfsphi[100]; + Int_t hfspid[100]; + TClonesArray *hfsp4 = new TClonesArray("TLorentzVector"); + TClonesArray &ar = *hfsp4; + TClonesArray *pip4 = new TClonesArray("TLorentzVector"); + TClonesArray &arpi = *pip4; + // TMVA for ML sidis reconstruction + std::vector> hfsinfo; + std::vector globalinfo; // particle masses static Double_t ElectronMass() { return 0.000511; }; @@ -269,13 +293,13 @@ class Kinematics : public TObject void CalculateDISbyMixed(); void CalculateDISbySigma(); void CalculateDISbyeSigma(); - + void CalculateDISbyML(); // calculate 4-momenta components of q and W (`vecQ` and `vecW`) as well as // derived invariants `W` and `nu` void GetQWNu_electronic(); void GetQWNu_hadronic(); void GetQWNu_quadratic(); - + void GetQWNu_ML(); private: static const Int_t asymInjectN = 2; Double_t moduVal[asymInjectN]; @@ -308,7 +332,14 @@ class Kinematics : public TObject Double_t rotAboutX, rotAboutY; // other TLorentzVector vecSpin, IvecSpin; - +#ifdef SIDIS_MLPRED + py::object keras, tensorflow; + py::object efnpackage; + py::function pfnimport; + py::object model; + py::object modelload; + std::string modelname = "pfn_testEpic_000-2_vecQele_nHFS2_500_bs10k_bestValLoss"; +#endif ClassDef(Kinematics,1); }; diff --git a/src/LinkDef.h b/src/LinkDef.h index bdb67e51..e9c023ff 100644 --- a/src/LinkDef.h +++ b/src/LinkDef.h @@ -17,6 +17,7 @@ // analysis objects #pragma link C++ class Kinematics+; #pragma link C++ class SimpleTree+; +#pragma link C++ class HFSTree+; #pragma link C++ class ParticleTree+; #pragma link C++ class Weights+; #pragma link C++ class WeightsUniform+; diff --git a/testEFlowimport.py b/testEFlowimport.py new file mode 100644 index 00000000..4b74c3a3 --- /dev/null +++ b/testEFlowimport.py @@ -0,0 +1,20 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Connor Pecar + +import tensorflow +from tensorflow import keras +import energyflow +from energyflow.archs import PFN +import numpy as np + +modelname = "pfn_testEpic_000-2_vecQele_nHFS2_500_bs10k_bestValLoss" +model = keras.models.load_model(modelname) + +def eflowPredict(feat, globalfeat): + feat = np.asarray(feat) + feat = np.reshape(feat, (1,len(feat),7)) + globalfeat = np.asarray(globalfeat) + globalfeat = np.reshape(globalfeat, (1,10)) + pred = model([feat,globalfeat]) + pred = np.reshape(pred,(4)) + return pred