diff --git a/.github/workflows/amptools.yml b/.github/workflows/amptools.yml index b902fee..9eb3d24 100644 --- a/.github/workflows/amptools.yml +++ b/.github/workflows/amptools.yml @@ -76,7 +76,20 @@ jobs: fi pushd extern/AmpTools; make -j$(nproc); popd make -j$(nproc) + - uses: actions/setup-python@v4 + with: + python-version: "3.10" + - run: pip install . - name: Run executables run: | - cd script - ./print_amplitudes -c zlm.cfg + ./print_amplitudes -c zlm.cfg > .print_amplitudes.log + python ./makeComparison.py + ./draw_zlm + working-directory: script + - uses: actions/upload-artifact@v3 + if: always() + with: + path: | + script/**.log + script/**.root + !script/data.root diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 00fdf9c..a676240 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -17,6 +17,11 @@ repos: - id: mixed-line-ending - id: trailing-whitespace + - repo: https://github.com/psf/black + rev: 23.7.0 + hooks: + - id: black + - repo: https://github.com/pre-commit/mirrors-clang-format rev: v16.0.6 hooks: @@ -26,6 +31,11 @@ repos: - c - cuda + - repo: https://github.com/PyCQA/isort + rev: 5.12.0 + hooks: + - id: isort + - repo: https://github.com/nbQA-dev/nbQA rev: 1.7.0 hooks: @@ -57,6 +67,6 @@ repos: metadata.vscode - repo: https://github.com/pre-commit/mirrors-prettier - rev: v3.0.0 + rev: v3.0.1 hooks: - id: prettier diff --git a/.vscode/extensions.json b/.vscode/extensions.json index 82a0632..8c616cb 100644 --- a/.vscode/extensions.json +++ b/.vscode/extensions.json @@ -1,5 +1,6 @@ { "recommendations": [ + "albertopdrf.root-file-viewer", "eamodio.gitlens", "editorconfig.editorconfig", "esbenp.prettier-vscode", diff --git a/pyproject.toml b/pyproject.toml index 076aae2..a7527c1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,9 @@ dependencies = [ "ipywidgets", "ipywidgets", "matplotlib", + "pandas", "sympy", + "uproot", ] name = "gluex-amplitude" requires-python = ">=3.7" diff --git a/script/.gitignore b/script/.gitignore index 9e30a54..a182956 100644 --- a/script/.gitignore +++ b/script/.gitignore @@ -1,5 +1,8 @@ * +!.gitignore !*.cc !*.cfg -!.gitignore +!*.md +!*.py +!data.root !Makefile diff --git a/script/Makefile b/script/Makefile index 533d132..c591899 100644 --- a/script/Makefile +++ b/script/Makefile @@ -4,7 +4,7 @@ LIB_DIR := -L$(GLUEXAMPLITUDE)/lib # add extra library linking flags here if they are needed -- be sure the locations # are specified in the LIB_DIR arugments above -EXTRA_LIBS := -lGlueXAmplitude +EXTRA_LIBS := -lGlueXAmplitude -lAmpTools$(ATSUFFIX) ###### most things below probably don't need adjustment @@ -25,7 +25,7 @@ CXX_FLAGS += $(shell root-config --cflags) # need some linking information available LIB_DIR += -L$(shell root-config --libdir) -L$(AMPTOOLS)/lib -L$(AMPPLOTTER)/lib -LIBS := -lAmpTools$(ATSUFFIX) $(EXTRA_LIBS) $(shell root-config --glibs) -lstdc++ +LIBS := $(EXTRA_LIBS) $(shell root-config --glibs) -lstdc++ # any source file falls into one of these two categories -- we'll # assume that anything that has MPI in the name is an MPI executable diff --git a/script/README.md b/script/README.md new file mode 100644 index 0000000..5b2fec0 --- /dev/null +++ b/script/README.md @@ -0,0 +1,8 @@ +# Usage + +The data comes from the spring 2017 run, thrown monte-carlo simulations used in the etapi analysis + +```shell +./print_amplitudes -c zlm.cfg > .print_amplitudes.log +python ./makeComparison.py +``` diff --git a/script/data.root b/script/data.root new file mode 100644 index 0000000..aa0f672 Binary files /dev/null and b/script/data.root differ diff --git a/script/draw_zlm.cc b/script/draw_zlm.cc new file mode 100644 index 0000000..1832abb --- /dev/null +++ b/script/draw_zlm.cc @@ -0,0 +1,94 @@ +#include "AMPTOOLS_AMPS/wignerD.h" +#include "TAxis.h" +#include "TCanvas.h" +#include "TF3.h" +#include "TH3F.h" +#include "TLorentzRotation.h" +#include "TLorentzVector.h" +#include "TMath.h" +#include "TVector3.h" + +// complex< GDouble > zlm(Double_t *xs, Int_t *pars) +Double_t zlm(Double_t *xs, Double_t *pars) { + // Parameters + Int_t m_j = (Int_t)pars[0]; // spin + Int_t m_m = (Int_t)pars[1]; // projection + Int_t m_r = (Int_t)pars[2]; // reflectivity + Int_t m_s = (Int_t)pars[3]; // sign + Int_t pAngle = (Int_t)pars[4]; // polarization angle + GDouble pGamma = (GDouble)pars[5]; // polarization fraction + GDouble factor = sqrt(1 + m_s * pGamma); + // cout << "j: " << m_j << endl; + // cout << "m: " << m_m << endl; + // cout << "reflectivity: " << m_r << endl; + // cout << "sign: " << m_s << endl; + // cout << "pAngle: " << pAngle << endl; + // cout << "pGamma: " << pGamma << endl; + // cout << "factor: " << factor << endl; + + // Variables + GDouble cosTheta = xs[0]; + GDouble phi = xs[1]; + GDouble bigPhi = xs[2]; + + // Calculate Intensity + GDouble zlm = 0; + complex rotateY = polar((GDouble)1., (GDouble)(-1. * bigPhi)); + if (m_r == 1) + zlm = real(Y(m_j, m_m, cosTheta, phi) * rotateY); + if (m_r == -1) + zlm = imag(Y(m_j, m_m, cosTheta, phi) * rotateY); + + complex amplitude(factor * zlm); + Double_t intensity = std::norm(amplitude); + return intensity; +} + +int main(int argc, char *argv[]) { + cout << "Drawing Zlm Amplitude" << endl; + // Double_t xs[] = {0,0,0}; + // Int_t pars[] = {2,0,1,1,0}; + // cout << zlm(xs, pars) << endl; // Test + + // Define Function + TF3 *function = new TF3("zlm", zlm, -1, 1, -3.14, 3.14, -3.14, 3.14, 6); + Int_t m_j = 2; + Int_t m_m = 0; + Int_t m_r = 1; + Int_t m_s = 1; + Int_t pAngle = 0; + GDouble pGamma = 1; + function->SetParameters(m_j, m_m, m_r, m_s, pAngle, pGamma); + + // Draw Function + TCanvas canvas = new TCanvas("canvas", "canvas", 600, 600); + int nBins = 50; + string par_string( + Form("l:%i m:%i e:%i s:%i a:%i g:%f;cos#theta;#phi [deg];#Phi [deg]", m_j, + m_m, m_r, m_s, pAngle, pGamma)); + TH3F *fcnHist = new TH3F("zlm", par_string.c_str(), nBins, -1, +1, nBins, + -3.14, 3.14, nBins, -3.14, 3.14); + TAxis *xAxis = fcnHist->GetXaxis(); + TAxis *yAxis = fcnHist->GetYaxis(); + TAxis *zAxis = fcnHist->GetZaxis(); + double x, y, z; + for (int xBin = 1; xBin < nBins + 1; ++xBin) { + for (int yBin = 1; yBin < nBins + 1; ++yBin) { + for (int zBin = 1; zBin < nBins + 1; ++zBin) { + x = xAxis->GetBinCenter(xBin); + y = yAxis->GetBinCenter(yBin); + z = zAxis->GetBinCenter(zBin); + fcnHist->SetBinContent(xBin, yBin, zBin, function->Eval(x, y, z)); + // cout << function->Eval(x, y, z) << endl; + } + } + } + fcnHist->SetMinimum(0); + fcnHist->GetXaxis()->SetTitleOffset(1.5); + fcnHist->GetYaxis()->SetTitleOffset(2); + fcnHist->GetZaxis()->SetTitleOffset(1.5); + fcnHist->Draw("BOX2"); + // canvas.SaveAs("zlm_3d.pdf"); + canvas.SaveAs("zlm_3d.root"); + return 0; +} diff --git a/script/makeComparison.py b/script/makeComparison.py new file mode 100644 index 0000000..20dd578 --- /dev/null +++ b/script/makeComparison.py @@ -0,0 +1,38 @@ +import re + +import numpy as np +import uproot + +# Load the data from the root file +data = uproot.open("data.root")["kin"].arrays( + ["cosTheta_eta_hel_thrown", "phi_eta_hel_thrown", "Phi_thrown"], library="pd" +) + +# Load the intensities extracted by print_amplitudes.cc +intensities = [] +with open(".print_amplitudes.log", "r") as file: + for line in file: + if "INTENSITY =" in line: + line = re.sub(" +", " ", line.rstrip().lstrip()) + intensities.append(float(line.split()[2])) +intensities = np.array(intensities) + +assert len(intensities) == len(data["cosTheta_eta_hel_thrown"]) + +# Rename dictionary keys +data["intensity"] = intensities +data = data.rename( + { + "cosTheta_eta_hel_thrown": "theta", + "phi_eta_hel_thrown": "phi", + "Phi_thrown": "Phi", + }, + axis=1, +) + +# Standardize the data to be angles in radians +DEGREE_TO_RADIAN = np.pi / 180.0 +data["theta"] = np.arccos(data["theta"]) +data["Phi"] = data["Phi"] * DEGREE_TO_RADIAN + +data.to_csv("comparison.csv", index=False, sep=",") diff --git a/script/print_amplitudes.cc b/script/print_amplitudes.cc index 69a564b..0937042 100644 --- a/script/print_amplitudes.cc +++ b/script/print_amplitudes.cc @@ -9,6 +9,7 @@ #include "TSystem.h" +#include "AMPTOOLS_AMPS/ROOTDataReader.h" #include "AMPTOOLS_AMPS/Zlm.h" #include "IUAmpTools/AmpToolsInterface.h" @@ -23,7 +24,6 @@ int main(int argc, char *argv[]) { for (int i = 1; i < argc; i++) { string arg(argv[i]); - if (arg == "-c") { configfile = argv[++i]; } @@ -36,25 +36,45 @@ int main(int argc, char *argv[]) { ConfigFileParser parser(configfile); ConfigurationInfo *cfgInfo = parser.getConfigurationInfo(); - cfgInfo->display(); + string reactionName = cfgInfo->reactionList()[0]->reactionName(); + + // cfgInfo->display(); AmpToolsInterface::registerAmplitude(Zlm()); + AmpToolsInterface::registerDataReader(ROOTDataReader()); AmpToolsInterface ati(cfgInfo); - vector p4List; - // list lab frame p4s, in the order specified in the AmpTools config file - p4List.push_back(TLorentzVector(0, 0, 8.209470, 8.209470)); // beam photon - p4List.push_back(TLorentzVector(-0.291810, -0.315238, 0.695296, - 1.244320)); // recoil proton - p4List.push_back(TLorentzVector(0.183147, 0.002392, 0.384435, - 0.448137)); // resonance daughter 1 - p4List.push_back(TLorentzVector(0.151104, 0.333728, 1.914300, - 1.953710)); // resonance daughter 2 + int iDataSet = 0; // if use used {data, accmc, genmc datasets}, we can choose + // one of them here - Kinematics kin(p4List); + ////////////////////////////////////////////////////////////////////////// + /// TESTING PURPOSES + // vector< TLorentzVector > p4List; + ////// list lab frame p4s, in the order specified in the AmpTools config file + // p4List.push_back( TLorentzVector(0.000000,0.000000,8.419198,8.419198) ); + // p4List.push_back( TLorentzVector(-0.105035,0.014094,1.106842,1.454883) ); + // p4List.push_back( TLorentzVector(-0.330669,-0.777849,0.458704,0.971092) ); + // p4List.push_back( TLorentzVector(0.435704,0.763756,6.855184,6.933026) ); + // Kinematics kin_test( p4List ); + // ati.printIntensity( reactionName, &kin_test ); + ////////////////////////////////////////////////////////////////////////// - ati.printEventDetails(cfgInfo->reactionList()[0]->reactionName(), &kin); + // Have to loadEvents before using kinematics (and numEvents) but cannot run + // processEvents. + // To use the intensity method, we need to processEvents first + Kinematics kin; + ati.loadEvents(ati.dataReader(reactionName), iDataSet); + int nentries = ati.numEvents(iDataSet); + printf("Number of entries: %i\n", nentries); + // cout << ati.processEvents(reactionName,iDataSet) << endl; + for (int i = 0; i < nentries; ++i) { + kin = *ati.kinematics(i, iDataSet); + // for (auto four_vector: kin.particleList()){ + // four_vector.Print(); + // } + ati.printIntensity(reactionName, &kin); + } return 0; } diff --git a/script/zlm.cfg b/script/zlm.cfg index 86c7189..8ca8a0b 100644 --- a/script/zlm.cfg +++ b/script/zlm.cfg @@ -1,91 +1,95 @@ -# -##################################### -#### THIS IS A CONFIG FILE #### -##################################### -## -## Blank lines or lines beginning with a "#" are ignored. -## -## Double colons (::) are treated like a space. -## This is sometimes useful for grouping (for example, -## grouping strings like "reaction::sum::amplitudeName") -## -## All non-comment lines must begin with one of the following keywords. -## -## (note: means necessary -## (word) means optional) -## -## include -## define (defn1) (defn2) (defn3) ... -## fit -## keyword -## reaction (particle3) ... -## data (arg1) (arg2) (arg3) ... -## genmc (arg1) (arg2) (arg3) ... -## accmc (arg1) (arg2) (arg3) ... -## normintfile -## sum (sum2) (sum3) ... -## amplitude (arg1) (arg2) ([par]) ... -## initialize <"events"/"polar"/"cartesian"> -## ("fixed"/"real") -## scale -## constrain ... -## permute ... -## parameter ("fixed"/"bounded"/"gaussian") -## (lower/central) (upper/error) -## DEPRECATED: -## datafile (file2) (file3) ... -## genmcfile (file2) (file3) ... -## accmcfile (file2) (file3) ... -## -##################################### - -#define beamConfig beamconfig.cfg -define polVal 1.0 -define polAngle 0.0 - -reaction EtaPi0_000 gamma Proton Pi0 Eta - -sum EtaPi0_000 PositiveRe -sum EtaPi0_000 PositiveIm - -# ------------------------------------------------ -# Define the amplitudes, both positive and negative reflectivity -# ------------------------------------------------ -# Zlm as suggested in GlueX doc-4094 (M. Shepherd) -# argument 1 : j -# argument 2 : m -# argument 3 : real (+1) or imaginary (-1) part -# argument 4 : 1 + (+1/-1) * P_gamma -# argument 5 : polarization angle -# argument 6 : beam properties config file or fixed polarization -#parameter polAngle 1.77 fixed - -# S-waves -amplitude EtaPi0_000::PositiveIm::S0+ Zlm 0 0 -1 -1 polAngle polVal -amplitude EtaPi0_000::PositiveRe::S0+ Zlm 0 0 +1 +1 polAngle polVal -constrain EtaPi0_000::PositiveIm::S0+ EtaPi0_000::PositiveRe::S0+ -initialize EtaPi0_000::PositiveIm::S0+ cartesian 1.0 0.0 real - -## Some D-waves -#amplitude EtaPi0_000::PositiveIm::D0+ Zlm 2 0 -1 -1 polAngle polVal -#amplitude EtaPi0_000::PositiveRe::D0+ Zlm 2 0 +1 +1 polAngle polVal -#amplitude EtaPi0_000::PositiveIm::D1+ Zlm 2 1 -1 -1 polAngle polVal -#amplitude EtaPi0_000::PositiveRe::D1+ Zlm 2 1 +1 +1 polAngle polVal -#amplitude EtaPi0_000::PositiveIm::D2+ Zlm 2 2 -1 -1 polAngle polVal -#amplitude EtaPi0_000::PositiveRe::D2+ Zlm 2 2 +1 +1 polAngle polVal -#initialize EtaPi0_000::PositiveIm::D0+ cartesian 1.0 0.0 -#initialize EtaPi0_000::PositiveIm::D1+ cartesian 1.0 0.0 -#initialize EtaPi0_000::PositiveIm::D2+ cartesian 1.0 0.0 -#constrain EtaPi0_000::PositiveIm::D0+ EtaPi0_000::PositiveRe::D0+ -#constrain EtaPi0_000::PositiveIm::D1+ EtaPi0_000::PositiveRe::D1+ -#constrain EtaPi0_000::PositiveIm::D2+ EtaPi0_000::PositiveRe::D2+ -# -## P-waves -#amplitude EtaPi0_000::PositiveIm::P0+ Zlm 1 0 -1 -1 polAngle polVal -#amplitude EtaPi0_000::PositiveRe::P0+ Zlm 1 0 +1 +1 polAngle polVal -#amplitude EtaPi0_000::PositiveIm::P1+ Zlm 1 1 -1 -1 polAngle polVal -#amplitude EtaPi0_000::PositiveRe::P1+ Zlm 1 1 +1 +1 polAngle polVal -#initialize EtaPi0_000::PositiveIm::P0+ cartesian 1.0 0.0 -#initialize EtaPi0_000::PositiveIm::P1+ cartesian 1.0 0.0 -#constrain EtaPi0_000::PositiveIm::P0+ EtaPi0_000::PositiveRe::P0+ -#constrain EtaPi0_000::PositiveIm::P1+ EtaPi0_000::PositiveRe::P1+ +# +##################################### +#### THIS IS A CONFIG FILE #### +##################################### +## +## Blank lines or lines beginning with a "#" are ignored. +## +## Double colons (::) are treated like a space. +## This is sometimes useful for grouping (for example, +## grouping strings like "reaction::sum::amplitudeName") +## +## All non-comment lines must begin with one of the following keywords. +## +## (note: means necessary +## (word) means optional) +## +## include +## define (defn1) (defn2) (defn3) ... +## fit +## keyword +## reaction (particle3) ... +## data (arg1) (arg2) (arg3) ... +## genmc (arg1) (arg2) (arg3) ... +## accmc (arg1) (arg2) (arg3) ... +## normintfile +## sum (sum2) (sum3) ... +## amplitude (arg1) (arg2) ([par]) ... +## initialize <"events"/"polar"/"cartesian"> +## ("fixed"/"real") +## scale +## constrain ... +## permute ... +## parameter ("fixed"/"bounded"/"gaussian") +## (lower/central) (upper/error) +## DEPRECATED: +## datafile (file2) (file3) ... +## genmcfile (file2) (file3) ... +## accmcfile (file2) (file3) ... +## +##################################### + +#define beamConfig beamconfig.cfg +define polVal 1.0 +define polAngle 0.0 + +#genmc EtaPi0 ROOTDataReader data.root +#accmc EtaPi0 ROOTDataReader data.root +data EtaPi0 ROOTDataReader data.root + +reaction EtaPi0 gamma Proton Pi0 Eta + +sum EtaPi0 PositiveRe +sum EtaPi0 PositiveIm + +# ------------------------------------------------ +# Define the amplitudes, both positive and negative reflectivity +# ------------------------------------------------ +# Zlm as suggested in GlueX doc-4094 (M. Shepherd) +# argument 1 : j +# argument 2 : m +# argument 3 : real (+1) or imaginary (-1) part +# argument 4 : 1 + (+1/-1) * P_gamma +# argument 5 : polarization angle +# argument 6 : beam properties config file or fixed polarization +#parameter polAngle 1.77 fixed + +# S-waves +amplitude EtaPi0::PositiveIm::S0+ Zlm 0 0 -1 -1 polAngle polVal +amplitude EtaPi0::PositiveRe::S0+ Zlm 0 0 +1 +1 polAngle polVal +constrain EtaPi0::PositiveIm::S0+ EtaPi0::PositiveRe::S0+ +initialize EtaPi0::PositiveIm::S0+ cartesian 1.0 0.0 real + +## Some D-waves +#amplitude EtaPi0::PositiveIm::D0+ Zlm 2 0 -1 -1 polAngle polVal +#amplitude EtaPi0::PositiveRe::D0+ Zlm 2 0 +1 +1 polAngle polVal +#amplitude EtaPi0::PositiveIm::D1+ Zlm 2 1 -1 -1 polAngle polVal +#amplitude EtaPi0::PositiveRe::D1+ Zlm 2 1 +1 +1 polAngle polVal +#amplitude EtaPi0::PositiveIm::D2+ Zlm 2 2 -1 -1 polAngle polVal +#amplitude EtaPi0::PositiveRe::D2+ Zlm 2 2 +1 +1 polAngle polVal +#initialize EtaPi0::PositiveIm::D0+ cartesian 1.0 0.0 +#initialize EtaPi0::PositiveIm::D1+ cartesian 1.0 0.0 +#initialize EtaPi0::PositiveIm::D2+ cartesian 1.0 0.0 +#constrain EtaPi0::PositiveIm::D0+ EtaPi0::PositiveRe::D0+ +#constrain EtaPi0::PositiveIm::D1+ EtaPi0::PositiveRe::D1+ +#constrain EtaPi0::PositiveIm::D2+ EtaPi0::PositiveRe::D2+ +# +## P-waves +#amplitude EtaPi0::PositiveIm::P0+ Zlm 1 0 -1 -1 polAngle polVal +#amplitude EtaPi0::PositiveRe::P0+ Zlm 1 0 +1 +1 polAngle polVal +#amplitude EtaPi0::PositiveIm::P1+ Zlm 1 1 -1 -1 polAngle polVal +#amplitude EtaPi0::PositiveRe::P1+ Zlm 1 1 +1 +1 polAngle polVal +#initialize EtaPi0::PositiveIm::P0+ cartesian 1.0 0.0 +#initialize EtaPi0::PositiveIm::P1+ cartesian 1.0 0.0 +#constrain EtaPi0::PositiveIm::P0+ EtaPi0::PositiveRe::P0+ +#constrain EtaPi0::PositiveIm::P1+ EtaPi0::PositiveRe::P1+ diff --git a/src/GlueXAmplitude/AMPTOOLS_AMPS/ROOTDataReader.cc b/src/GlueXAmplitude/AMPTOOLS_AMPS/ROOTDataReader.cc new file mode 100644 index 0000000..2ae5e42 --- /dev/null +++ b/src/GlueXAmplitude/AMPTOOLS_AMPS/ROOTDataReader.cc @@ -0,0 +1,97 @@ + +#include +#include +#include + +#include "TLorentzVector.h" + +#include "AMPTOOLS_AMPS/ROOTDataReader.h" +#include "IUAmpTools/Kinematics.h" + +#include "TFile.h" +#include "TH1.h" +#include "TTree.h" + +using namespace std; + +ROOTDataReader::ROOTDataReader(const vector &args) + : UserDataReader(args), m_eventCounter(0), + m_useWeight(false) { + assert(args.size() == 2 || args.size() == 1); + + TH1::AddDirectory(kFALSE); + + // this way of opening files works with URLs of the form + // root://xrootdserver/path/to/myfile.root + m_inFile = TFile::Open(args[0].c_str()); + + // default to tree name of "kin" if none is provided + if (args.size() == 1) { + + m_inTree = dynamic_cast(m_inFile->Get("kin")); + } else { + + m_inTree = dynamic_cast(m_inFile->Get(args[1].c_str())); + } + + m_inTree->SetBranchAddress("NumFinalState", &m_nPart); + m_inTree->SetBranchAddress("E_FinalState", m_e); + m_inTree->SetBranchAddress("Px_FinalState", m_px); + m_inTree->SetBranchAddress("Py_FinalState", m_py); + m_inTree->SetBranchAddress("Pz_FinalState", m_pz); + m_inTree->SetBranchAddress("E_Beam", &m_eBeam); + m_inTree->SetBranchAddress("Px_Beam", &m_pxBeam); + m_inTree->SetBranchAddress("Py_Beam", &m_pyBeam); + m_inTree->SetBranchAddress("Pz_Beam", &m_pzBeam); + + if (m_inTree->GetBranch("Weight") != NULL) { + + m_useWeight = true; + m_inTree->SetBranchAddress("Weight", &m_weight); + } else { + + m_useWeight = false; + } +} + +ROOTDataReader::~ROOTDataReader() { + if (m_inFile != NULL) + m_inFile->Close(); +} + +void ROOTDataReader::resetSource() { + + cout << "Resetting source " << m_inTree->GetName() << " in " + << m_inFile->GetName() << endl; + + // this will cause the read to start back at event 0 + m_eventCounter = 0; +} + +Kinematics *ROOTDataReader::getEvent() { + if (m_eventCounter < static_cast(m_inTree->GetEntries())) { + // if( m_eventCounter < 10 ){ + + m_inTree->GetEntry(m_eventCounter++); + assert(m_nPart < Kinematics::kMaxParticles); + + vector particleList; + + particleList.push_back( + TLorentzVector(m_pxBeam, m_pyBeam, m_pzBeam, m_eBeam)); + + for (int i = 0; i < m_nPart; ++i) { + + particleList.push_back(TLorentzVector(m_px[i], m_py[i], m_pz[i], m_e[i])); + } + + return new Kinematics(particleList, m_useWeight ? m_weight : 1.0); + } else { + + return NULL; + } +} + +unsigned int ROOTDataReader::numEvents() const { + return static_cast(m_inTree->GetEntries()); +} diff --git a/src/GlueXAmplitude/AMPTOOLS_AMPS/ROOTDataReader.h b/src/GlueXAmplitude/AMPTOOLS_AMPS/ROOTDataReader.h new file mode 100644 index 0000000..96935df --- /dev/null +++ b/src/GlueXAmplitude/AMPTOOLS_AMPS/ROOTDataReader.h @@ -0,0 +1,62 @@ +#if !defined(ROOTDATAREADER) +#define ROOTDATAREADER + +#include "IUAmpTools/Kinematics.h" +#include "IUAmpTools/UserDataReader.h" + +#include "TFile.h" +#include "TString.h" +#include "TTree.h" + +#include + +using namespace std; + +class ROOTDataReader : public UserDataReader { + +public: + /** + * Default constructor for ROOTDataReader + */ + ROOTDataReader() : UserDataReader(), m_inFile(NULL) {} + + ~ROOTDataReader(); + + /** + * Constructor for ROOTDataReader + * \param[in] args vector of string arguments + */ + ROOTDataReader(const vector &args); + + string name() const { return "ROOTDataReader"; } + + virtual Kinematics *getEvent(); + virtual void resetSource(); + + /** + * This function returns a true if the file was open + * with weight-reading enabled and had this tree branch, + * false, if these criteria are not met. + */ + virtual bool hasWeight() { return m_useWeight; }; + virtual unsigned int numEvents() const; + +private: + TFile *m_inFile; + TTree *m_inTree; + unsigned int m_eventCounter; + bool m_useWeight; + + int m_nPart; + float m_e[Kinematics::kMaxParticles]; + float m_px[Kinematics::kMaxParticles]; + float m_py[Kinematics::kMaxParticles]; + float m_pz[Kinematics::kMaxParticles]; + float m_eBeam; + float m_pxBeam; + float m_pyBeam; + float m_pzBeam; + float m_weight; +}; + +#endif