diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 392ecd56..ddafd052 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -130,21 +130,22 @@ get_data: include: # - local: 'benchmarks/backgrounds/config.yml' - - local: 'benchmarks/ecal_gaps/config.yml' - - local: 'benchmarks/tracking_detectors/config.yml' - - local: 'benchmarks/barrel_ecal/config.yml' - - local: 'benchmarks/barrel_hcal/config.yml' - - local: 'benchmarks/zdc/config.yml' - - local: 'benchmarks/material_maps/config.yml' - - local: 'benchmarks/material_scan/config.yml' - - local: 'benchmarks/pid/config.yml' - - local: 'benchmarks/timing/config.yml' - - local: 'benchmarks/b0_tracker/config.yml' +# - local: 'benchmarks/ecal_gaps/config.yml' +# - local: 'benchmarks/tracking_detectors/config.yml' +# - local: 'benchmarks/barrel_ecal/config.yml' +# - local: 'benchmarks/barrel_hcal/config.yml' +# - local: 'benchmarks/zdc/config.yml' +# - local: 'benchmarks/material_maps/config.yml' +# - local: 'benchmarks/material_scan/config.yml' +# - local: 'benchmarks/pid/config.yml' +# - local: 'benchmarks/timing/config.yml' +# - local: 'benchmarks/b0_tracker/config.yml' + - local: 'benchmarks/LOWQ2/analysis/config.yml' deploy_results: stage: deploy needs: - - ["collect_results:zdc","collect_results:barrel_ecal","collect_results:barrel_hcal","collect_results:material_scan"] +# - ["collect_results:zdc","collect_results:barrel_ecal","collect_results:barrel_hcal","collect_results:material_scan"] script: - echo "deploy results!" - find results -print | sort | tee summary.txt @@ -222,5 +223,3 @@ pages: artifacts: paths: - public - - diff --git a/benchmarks/LOWQ2/analysis/.gitignore b/benchmarks/LOWQ2/analysis/.gitignore new file mode 100644 index 00000000..c96503c7 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/.gitignore @@ -0,0 +1,5 @@ +plots* +calibrations* +fieldmaps* +gdml* +.snakemake* \ No newline at end of file diff --git a/benchmarks/LOWQ2/analysis/LOWQ2Benchmarks.h b/benchmarks/LOWQ2/analysis/LOWQ2Benchmarks.h new file mode 100755 index 00000000..13ac9f4d --- /dev/null +++ b/benchmarks/LOWQ2/analysis/LOWQ2Benchmarks.h @@ -0,0 +1,289 @@ +#include "TString.h" +#include "TCanvas.h" +#include "ROOT/RDataFrame.hxx" +#include +#include "TH1.h" +#include "TFile.h" +#include "LOWQ2hits.h" +#include "LOWQ2acceptance.h" +#include "LOWQ2clusters.h" +#include "LOWQ2reconstruction.h" + +// Define alias +using RNode = ROOT::RDF::RNode; +using H1ResultPtr = ROOT::RDF::RResultPtr; +using H2ResultPtr = ROOT::RDF::RResultPtr; +using H3ResultPtr = ROOT::RDF::RResultPtr; + +using RVecS = ROOT::VecOps::RVec; +using RVecI = ROOT::VecOps::RVec; + +std::map,std::map,std::map>> histMap; + +//--------------------------------------------------------------------------------------------- +// Create dataframe from input file(s) +//--------------------------------------------------------------------------------------------- +RNode initialise( string fileName ){ + + ROOT::RDataFrame d0("events",fileName); + return d0; + +} + +//--------------------------------------------------------------------------------------------- +// Create histograms showing the variation of flux in various system components +//--------------------------------------------------------------------------------------------- +void WriteFluxPlots(H2ResultPtr hist, std::string parts){ + // Make projection of 2D pixel binned histogram + auto nBins = hist->GetNcells(); + TString pixelFluxName = TString(parts)+"Flux"; + TString logPixelFluxName = TString(parts)+"LogFlux"; + + //Get maximum bin content to set range + double maxBinContent = hist->GetMaximum(); + double logMax = log10(maxBinContent); + + //Create pixel flux histogram + TH1D* pixelFlux = new TH1D(pixelFluxName,pixelFluxName,100,0,maxBinContent); + + //Create log pixel flux histogram + TH1D* logPixelFlux = new TH1D(logPixelFluxName,logPixelFluxName,100,0,logMax); + for(int i=0; iFill(hist->GetBinContent(i)); + logPixelFlux->Fill(log10(hist->GetBinContent(i))); + } + pixelFlux->GetXaxis()->SetTitle("N_{hits}"); + logPixelFlux->GetXaxis()->SetTitle("log(N_{hits})"); + pixelFlux->GetYaxis()->SetTitle("N_{pixels}"); + logPixelFlux->GetYaxis()->SetTitle("N_{pixels}"); + + pixelFlux->Write(); + logPixelFlux->Write(); + +} + +//--------------------------------------------------------------------------------------------- +//Create histograms showing the variation of resolution accross the kinematic phase space +//--------------------------------------------------------------------------------------------- +void WriteResolutionPlots(H2ResultPtr hist, std::string parts){ + + hist->FitSlicesX(nullptr,0,-1,5); + TH1D* h1 = (TH1D*)gDirectory->Get(TString(hist->GetName())+"_2"); + h1->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + h1->GetYaxis()->SetTitle(TString(hist->GetYaxis()->GetTitle())+"Resolution"); + + h1->Write(); + +} + + +//--------------------------------------------------------------------------------------------- +// Format and write plots +//--------------------------------------------------------------------------------------------- +void writePlots( TString outName ){ + + TFile* rootfile = new TFile(outName,"RECREATE"); + auto benchmarkDir = rootfile->mkdir("LOWQ2"); + + //Loop over 1D histograms saving + for(auto &[bmkey,hists]: histMap){ + + //Split mapped distributions name into parts + std::stringstream ss(bmkey.Data()); + std::string part; + std::vector parts; + while (std::getline(ss, part, '/')) { + parts.push_back(part); + } + + TDirectory* dir = benchmarkDir; + for (size_t i = 0; i < parts.size(); ++i) { + // This is not the last part, create or get the directory + if (!dir->GetDirectory(parts[i].c_str())) { + dir = dir->mkdir(parts[i].c_str()); + } else { + dir = dir->GetDirectory(parts[i].c_str()); + } + dir->cd(); + } + + TDirectory* currentDir = gDirectory; + + for(auto &[key,hist]: get<0>(hists)){ + + //Split histogram name into parts + std::stringstream ss2(key.Data()); + std::string part2; + std::vector parts2; + while (std::getline(ss2, part2, '/')) { + parts2.push_back(part2); + } + + TDirectory* dir = currentDir; + for (size_t i = 0; i < parts2.size(); ++i) { + if (i == parts2.size() - 1) { + // This is the last part, write the histogram + hist->SetMinimum(0); + hist->Write(parts2[i].c_str()); + } else { + // This is not the last part, create or get the directory + if (!dir->GetDirectory(parts2[i].c_str())) { + dir = dir->mkdir(parts2[i].c_str()); + } else { + dir = dir->GetDirectory(parts2[i].c_str()); + } + dir->cd(); + } + } + currentDir->cd(); + // hist->Write(key); + } + + //Loop over 2D histograms saving and calculating fluxes as needed + for(auto &[key,hist]: get<1>(hists)){ + + TDirectory* currentDir = gDirectory; + std::stringstream ss(key.Data()); + std::string part; + TDirectory* dir = currentDir; + std::vector parts; + while (std::getline(ss, part, '/')) { + parts.push_back(part); + } + for (size_t i = 0; i < parts.size(); ++i) { + if (i == parts.size() - 1) { + // If histogram name contains rate or hits then calculate flux + if(parts[i].find("Rate") != std::string::npos || parts[i].find("Hits") != std::string::npos){ + WriteFluxPlots(hist,parts[i].c_str()); + } + + // If histogram name contains Res then create a profile + if(parts[i].find("Res") != std::string::npos){ + WriteResolutionPlots(hist,parts[i].c_str()); + } + + hist->Sumw2(false); + hist->Write(); + } else { + // This is not the last part, create or get the directory + if (!dir->GetDirectory(parts[i].c_str())) { + dir = dir->mkdir(parts[i].c_str()); + } else { + dir = dir->GetDirectory(parts[i].c_str()); + } + dir->cd(); + } + } + currentDir->cd(); + } + + //Loop over 3D histograms saving + for(auto &[key,hist]: get<2>(hists)){ + TDirectory* currentDir = gDirectory; + std::stringstream ss(key.Data()); + std::string part; + TDirectory* dir = currentDir; + std::vector parts; + while (std::getline(ss, part, '/')) { + parts.push_back(part); + } + for (size_t i = 0; i < parts.size(); ++i) { + if (i == parts.size() - 1) { + // This is the last part, write the histogram + hist->Write(parts[i].c_str()); + + //Fit histogram z slices and save 2d histogram + hist->FitSlicesZ(nullptr,0,-1,5); + TH2D* h2 = (TH2D*)gDirectory->Get(TString(hist->GetName())+"_2"); + h2->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + h2->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle()); + h2->GetZaxis()->SetTitle(TString(hist->GetZaxis()->GetTitle())+"Resolution"); + h2->SetTitle(hist->GetTitle()); + h2->Write(); + + } else { + // This is not the last part, create or get the directory + if (!dir->GetDirectory(parts[i].c_str())) { + dir = dir->mkdir(parts[i].c_str()); + } else { + dir = dir->GetDirectory(parts[i].c_str()); + } + dir->cd(); + } + } + currentDir->cd(); + // hist->Write(key); + } + + benchmarkDir->cd(); + + } + + rootfile->Close(); + +} + +//--------------------------------------------------------------------------------------------- +// Create the benchmark plots +//--------------------------------------------------------------------------------------------- +void LOWQ2Benchmarks( string inName, TString outName, dd4hep::Detector& detector, double eventRate ){ + + auto node = initialise( inName ); + + int events = *node.Count(); + double eventWeight = eventRate / events; + + RVecS colNames = node.GetColumnNames(); + + std::string readoutName = "TaggerTrackerHits"; + + if(Any(colNames==readoutName)){ + //----------------------------------------- + // Hit detector IDs + //----------------------------------------- + auto ids = detector.readout(readoutName).idSpec().fields(); + + node = node.Define("hitParticleIndex","_TaggerTrackerHits_MCParticle.index") + .Define("generatorStat","MCParticles.generatorStatus") + .Define("primary",[](RVecI index, RVecI status){ + RVecI result(index.size()); + for (size_t i = 0; i < index.size(); ++i) { + result[i] = status[index[i]] == 1; + } + return result; + },{"hitParticleIndex","generatorStat"}); + + std::string filterName = "TaggerTrackerHitsFilter"; + auto allnode = node.Define(filterName,"TaggerTrackerHits"); + auto primarynode = node.Define(filterName,"TaggerTrackerHits[primary==1]"); + auto secondarynode = node.Define(filterName,"TaggerTrackerHits[primary==0]"); + + for(auto &[key,id]: ids){ + TString colName = key+"ID"; + node = node.Define(colName,getSubID(key,detector),{readoutName}); + primarynode = primarynode.Define(colName,getSubID(key,detector),{filterName}); + secondarynode = secondarynode.Define(colName,getSubID(key,detector),{filterName}); + } + + //Create Rate Plots + histMap["Rates/AllHits"] = createHitPlots(node,eventWeight); + histMap["Rates/PrimaryHits"] = createHitPlots(primarynode,eventWeight); + histMap["Rates/SecondaryHits"] = createHitPlots(secondarynode,eventWeight); + + } + + if((Any(colNames==readoutName) || Any(colNames=="InclusiveKinematicsElectron")) && Any(colNames=="MCParticles")){ + histMap["Acceptance"] = createAcceptancePlots(node); + } + + if(Any(colNames=="TaggerTrackerClusterPositions")){ + histMap["Clusters"] = createClusterPlots(node); + } + + if(Any(colNames=="LowQ2TrackParameters") && Any(colNames=="MCParticles")){ + histMap["Reconstruction"] = createReconstructionPlots(node); + } + + writePlots( outName ); + +} \ No newline at end of file diff --git a/benchmarks/LOWQ2/analysis/LOWQ2acceptance.h b/benchmarks/LOWQ2/analysis/LOWQ2acceptance.h new file mode 100644 index 00000000..d7d1037b --- /dev/null +++ b/benchmarks/LOWQ2/analysis/LOWQ2acceptance.h @@ -0,0 +1,186 @@ +#pragma once + +// #include "functors.h" + +#include "ROOT/RDF/RInterface.hxx" +#include +#include +#include "ROOT/RVec.hxx" + +// Define alias +using RNode = ROOT::RDF::RNode; +using H1ResultPtr = ROOT::RDF::RResultPtr; +using H2ResultPtr = ROOT::RDF::RResultPtr; +using H3ResultPtr = ROOT::RDF::RResultPtr; +using RVecI = ROOT::VecOps::RVec; +using RVecD = ROOT::VecOps::RVec; +using RVecS = ROOT::VecOps::RVec; + +// Lazy execution methods +std::tuple,std::map,std::map> createAcceptancePlots( RNode d1 ){ + + std::map hHists1D; + std::map hHists2D; + std::map hHists3D; + + int ePrimID = 2; + int pBeamID = 1; + int eBeamID = 0; + + // d1 = d1.Define("primMom","MCParticles[2].momentum.z"); + d1 = d1.Define("eBeam", [eBeamID] + (std::vector h) + {return ROOT::Math::PxPyPzMVector (h[eBeamID].momentum.x,h[eBeamID].momentum.y,h[eBeamID].momentum.z,h[eBeamID].mass);}, + {"MCParticles"}) + .Define("pBeam", [pBeamID] + (std::vector h) + {return ROOT::Math::PxPyPzMVector (h[pBeamID].momentum.x,h[pBeamID].momentum.y,h[pBeamID].momentum.z,h[pBeamID].mass);}, + {"MCParticles"}) + .Define("primMom", [ePrimID] + (std::vector h) + {return ROOT::Math::PxPyPzMVector (h[ePrimID].momentum.x,h[ePrimID].momentum.y,h[ePrimID].momentum.z,h[ePrimID].mass);}, + {"MCParticles"}) + .Define("primE","primMom.E()") + .Define("primTheta","1000*(TMath::Pi()-primMom.Theta())") + .Define("primEta","primMom.Eta()") + .Define("scatteredV","eBeam-primMom") + .Define("Q2","-scatteredV.M2()") + .Define("logQ2","log10(Q2)") + .Define("x","Q2/(2*scatteredV.Dot(pBeam))") + .Define("logx","log10(x)") + .Define("primPhi","primMom.Phi()") + .Define("W2","(pBeam+eBeam-primMom).M2()"); + + RVecS colNames = d1.GetColumnNames(); + std::vector> filters; + filters.push_back(std::make_pair("Generated",d1)); + filters.push_back(std::make_pair("Theta<10mrad",d1.Filter("primTheta<10"))); + if(Any(colNames=="TaggerTrackerHits")){ + filters.push_back(std::make_pair("Module1-Hit",d1.Filter("moduleID[moduleID==1].size()>0"))); + filters.push_back(std::make_pair("Module2-Hit",d1.Filter("moduleID[moduleID==2].size()>0"))); + filters.push_back(std::make_pair("Any-Hit",d1.Filter("TaggerTrackerHits.size() > 0"))); + } + + if(Any(colNames=="LowQ2Tracks")){ + filters.push_back(std::make_pair("Reconstructed-Track",d1.Filter("LowQ2Tracks.size()>0"))); + } + + + for (auto& [rawString,filterNode] : filters) + { + + //Histogram names + TString energyHistName = "hEnergy"; + TString thetaHistName = "hTheta"; + TString etaHistName = "hEta"; + TString logQ2HistName = "hlogQ2"; + TString logxHistName = "hlogx"; + TString Q2HistName = "hQ2"; + TString xHistName = "hx"; + TString W2HistName = "hW2"; + TString Q2xHistName = "hQ2x"; + TString logQ2logxHistName = "hlogQ2logx"; + TString WlogQ2HistName = "hWlogQ2"; + TString WlogxHistName = "hWlogx"; + TString primElogQ2HistName = "hprimElogQ2"; + TString primEthetaHistName = "hprimEtheta"; + + //Histogram Ranges + double energyMin = 0; + double energyMax = 18; + double thetaMin = 0; + double thetaMax = 12; + double etaMin = -15; + double etaMax = 5; + double logQ2Min = -8.5; + double logQ2Max = 0; + double logxMin = -12; + double logxMax = 0; + double Q2Min = 0; + double Q2Max = 500; + double xMin = 0; + double xMax = 1; + double W2Min = 0; + double W2Max = 500; + + int nBins1D = 400; + int nBins2D = 200; + + //Histogram axis names + TString energyAxisName = "Energy [GeV]"; + TString thetaAxisName = "Theta [mrad]"; + TString etaAxisName = "Eta"; + TString logQ2AxisName = "log(Q2) [GeV]"; + TString logxAxisName = "log(x) [GeV]"; + TString Q2AxisName = "Q2 [GeV]"; + TString xAxisName = "x [GeV]"; + TString W2AxisName = "W2 [GeV]"; + TString Q2xAxisName = Q2AxisName+";"+xAxisName; + TString logQ2logxAxisName = logQ2AxisName+";"+logxAxisName; + TString WlogQ2AxisName = W2AxisName+";"+logQ2AxisName; + TString WlogxAxisName = W2AxisName+";"+logxAxisName; + TString primElogQ2AxisName = energyAxisName+";"+logQ2AxisName; + TString primEthetaAxisName = energyAxisName+";"+thetaAxisName; + + + hHists1D[rawString+"/"+energyHistName] = filterNode.Histo1D({energyHistName, energyHistName + ";Energy [GeV]", nBins1D, energyMin, energyMax}, "primE"); + hHists1D[rawString+"/"+thetaHistName] = filterNode.Histo1D({thetaHistName, thetaHistName + ";Theta [mrad]", nBins1D, thetaMin, thetaMax}, "primTheta"); + hHists1D[rawString+"/"+etaHistName] = filterNode.Histo1D({etaHistName, etaHistName + ";Eta", nBins1D, etaMin, etaMax}, "primEta"); + hHists1D[rawString+"/"+logQ2HistName] = filterNode.Histo1D({logQ2HistName, logQ2HistName + ";log(Q2) [GeV]", nBins1D, logQ2Min, logQ2Max}, "logQ2"); + hHists1D[rawString+"/"+logxHistName] = filterNode.Histo1D({logxHistName, logxHistName + ";log(x) [GeV]", nBins1D, logxMin, logxMax}, "logx"); + hHists1D[rawString+"/"+Q2HistName] = filterNode.Histo1D({Q2HistName, Q2HistName + ";Q2 [GeV]", nBins1D, Q2Min, Q2Max}, "Q2"); + hHists1D[rawString+"/"+xHistName] = filterNode.Histo1D({xHistName, xHistName + ";x [GeV]", nBins1D, xMin, xMax}, "x"); + hHists1D[rawString+"/"+W2HistName] = filterNode.Histo1D({W2HistName, W2HistName + ";W2 [GeV]", nBins1D, W2Min, W2Max}, "W2"); + + hHists2D[rawString+"/"+Q2xHistName] = filterNode.Histo2D({Q2xHistName, Q2xHistName + ";Q2 [GeV];x [GeV]", nBins2D, Q2Min, Q2Max, nBins2D, xMin, xMax}, "Q2", "x"); + hHists2D[rawString+"/"+logQ2logxHistName] = filterNode.Histo2D({logQ2logxHistName, logQ2logxHistName + ";log(Q2) [GeV];log(x) [GeV]", nBins2D, logQ2Min, logQ2Max, nBins2D, logxMin, logxMax}, "logQ2", "logx"); + hHists2D[rawString+"/"+WlogQ2HistName] = filterNode.Histo2D({WlogQ2HistName, WlogQ2HistName + ";W2 [GeV];log(Q2) [GeV]", nBins2D, W2Min, W2Max, nBins2D, logQ2Min, logQ2Max}, "W2", "logQ2"); + hHists2D[rawString+"/"+WlogxHistName] = filterNode.Histo2D({WlogxHistName, WlogxHistName + ";W2 [GeV];log(x) [GeV]", nBins2D, W2Min, W2Max, nBins2D, logxMin, logxMax}, "W2", "logx"); + hHists2D[rawString+"/"+primElogQ2HistName] = filterNode.Histo2D({primElogQ2HistName, primElogQ2HistName + ";E [GeV];log(Q2) [GeV]", nBins2D, energyMin, energyMax, nBins2D, logQ2Min, logQ2Max}, "primE", "logQ2"); + hHists2D[rawString+"/"+primEthetaHistName] = filterNode.Histo2D({primEthetaHistName, primEthetaHistName + ";E [GeV];Theta [mrad]", nBins2D, energyMin, energyMax, nBins2D, thetaMin, thetaMax}, "primE", "primTheta"); + + } + + //---------------------------------------------------------------------------------------------------- + // Creates integrated acceptance plots for various filters + //---------------------------------------------------------------------------------------------------- + + //Store the counts for each filter + RVecI filterCounts; + int rawCount = 0; + TString filtersName = ""; + + //Calculate the counts for each filter + for (auto& [rawString,filterNode] : filters) + { + filtersName += rawString + "_"; + int counts = *filterNode.Count(); + std::cout << rawString << " " << counts << std::endl; + filterCounts.push_back(counts); + if (rawString == "Generated") { + rawCount = counts; + } + } + + filtersName = filtersName(0,filtersName.Length()-1); + + // Number of filters + int Nfilters = filterCounts.size(); + RVecI entryInt = ROOT::VecOps::Range(Nfilters); + + // Fraction of events passing each filter + RVecD fractions = filterCounts/static_cast(rawCount); + + ROOT::RDataFrame df(Nfilters); + auto dfWithCounts = df.Define("entry", [entryInt](ULong64_t entry) { return (double)entryInt[(int)entry]; }, {"rdfentry_"}) + .Define("count", [filterCounts](ULong64_t entry) { return (double)filterCounts[(int)entry]; }, {"rdfentry_"}) + .Define("fraction", [fractions](ULong64_t entry) { return (double)fractions[(int)entry]; }, {"rdfentry_"}); + + // Use the new column in the histograms + hHists1D["TotalCounts"] = dfWithCounts.Histo1D({"hTotalCounts", "Total Counts;"+filtersName+";Count", Nfilters, 0, static_cast(Nfilters)}, "entry","count"); + hHists1D["IntegratedAcceptance"] = dfWithCounts.Histo1D({"hFractions", "Fractions;"+filtersName+";Fraction", Nfilters, 0, static_cast(Nfilters)}, "entry","fraction"); + hHists1D["TotalCounts"]->Sumw2(false); + hHists1D["IntegratedAcceptance"]->Sumw2(false); + + return {hHists1D,hHists2D,hHists3D}; +} \ No newline at end of file diff --git a/benchmarks/LOWQ2/analysis/LOWQ2clusters.h b/benchmarks/LOWQ2/analysis/LOWQ2clusters.h new file mode 100644 index 00000000..2022cc40 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/LOWQ2clusters.h @@ -0,0 +1,21 @@ +// Define alias +using RNode = ROOT::RDF::RNode; +using H1ResultPtr = ROOT::RDF::RResultPtr; +using H2ResultPtr = ROOT::RDF::RResultPtr; +using H3ResultPtr = ROOT::RDF::RResultPtr; +using RVecI = ROOT::VecOps::RVec; + +// Lazy execution methods +std::tuple,std::map,std::map> createClusterPlots( RNode d1 ){ + + std::map hHists1D; + std::map hHists2D; + std::map hHists3D; + + // auto d2 = d1.Define("ClusterSize",[](RVecI clu){return clu.size();},{"TaggerTrackerClusterPositions_0.index"}); + auto d2 = d1.Define("ClusterSize","TaggerTrackerClusterPositions.rawHits_end-TaggerTrackerClusterPositions.rawHits_begin"); + hHists1D["hClusterSize"] = d2.Histo1D({"hClusterSize","hClusterSize",100,0,100}, "ClusterSize"); + + return {hHists1D,hHists2D,hHists3D}; + +} diff --git a/benchmarks/LOWQ2/analysis/LOWQ2hits.h b/benchmarks/LOWQ2/analysis/LOWQ2hits.h new file mode 100644 index 00000000..ac3a26df --- /dev/null +++ b/benchmarks/LOWQ2/analysis/LOWQ2hits.h @@ -0,0 +1,127 @@ +#pragma once + +#include "functors.h" + +// Define alias +using RNode = ROOT::RDF::RNode; +using H1ResultPtr = ROOT::RDF::RResultPtr; +using H2ResultPtr = ROOT::RDF::RResultPtr; +using H3ResultPtr = ROOT::RDF::RResultPtr; +using RVecI = ROOT::VecOps::RVec; +using RVecD = ROOT::VecOps::RVec; + +// Lazy execution methods +std::tuple,std::map,std::map> createHitPlots( RNode d2, double eventWeight ){ + + int xChip = 448; + int yChip = 512; + int xChips = 6; + int yChips = 6; + std::map xPixelMin = {{1,-xChip*xChips/2},{2,-xChip*xChips/2}}; + std::map xPixelMax = {{1, xChip*xChips/2},{2, xChip*xChips/2}}; + + std::map yPixelMin = {{1,-yChip*yChips/2},{2,-yChip*yChips/2}}; + std::map yPixelMax = {{1, yChip*yChips/2},{2, yChip*yChips/2}}; + + + std::map hHists1D; + std::map hHists2D; + std::map hHists3D; + + d2 = d2.Define("boardID", [xChip](RVecI xID){ return (xID + 3*xChip) / (2 * xChip); },{"xID"}) + .Define("xChipID", [xChip](RVecI xID){ return (xID + 3*xChip) / (xChip); },{"xID"}) + .Define("xColumnID", [xChip](RVecI xID){ return (xID + 3*xChip) / 2; },{"xID"}) + .Define("yChipID", [yChip](RVecI yID){ return (yID + 2.75*yChip) / (yChip); },{"yID"}) + .Define("yHalfChipID", [yChip](RVecI yID){ return (yID + 2.75*yChip) / (0.5*yChip); },{"yID"}) + .Define("eventWeight", [eventWeight](){return eventWeight;},{}); + + + // Plot number of hits per moduleID + hHists1D["hmoduleIDRate"] = d2.Histo1D({"hmoduleIDRate", "hmoduleID;Module ID;Rate [Hz]", 2, 1, 3 }, "moduleID", "eventWeight"); + + // Module Loop + for(int i=1; i<=2; i++){ + + double xMin = xPixelMin[i]; + double xMax = xPixelMax[i]; + double yMin = yPixelMin[i]; + double yMax = yPixelMax[i]; + int xRange = xMax-xMin; + int yRange = yMax-yMin; + + TString modNum = std::to_string(i); + TString moduleTag = "module"+modNum; + TString layerName = "layerID"+moduleTag; + TString layerHistName = "h"+layerName+"Rate"; + + + auto d3 = d2.Define("ModuleFilter",[i](RVecI modID){return modID==i;},{"moduleID"}) + .Define(layerName,"layerID[ModuleFilter]"); + + + hHists1D[moduleTag+"/"+layerHistName] = d3.Histo1D({layerHistName, layerHistName+";Layer ID;Rate [Hz]", 4, 0, 4 }, layerName, "eventWeight"); + + // Layer Loop + for(int j=0; j<=3; j++){ + + TString layerNum = std::to_string(j); + TString layerTag = "layer"+layerNum; + + TString xName = "xPixel"; + TString xHistName = "h"+xName+"Rate"; + + TString yName = "yPixel"; + TString yHistName = "h"+yName+"Rate"; + + TString xChipName = "xChip"; + TString yChipName = "yChip"; + + TString yHalfChipName = "yHalfChip"; + + TString xyHistName = "h"+xName+yName+"Rate"; + + TString xColumnName = "xColumn"; + + TString boardName = "board"; + + std::vector layerSizeInput = {xName.Data()}; + TString layerSizeName = "HitsPerEvent"; + TString sizeHistName = "h"+layerSizeName; + + auto d4 = d3.Define("LayerFilter",[j](RVecI layID){return layID==j;},{"layerID"}) + .Define(xName,"xID[LayerFilter&&ModuleFilter]") + .Define(yName,"yID[LayerFilter&&ModuleFilter]") + .Define(boardName,"boardID[LayerFilter&&ModuleFilter]") + .Define(xChipName,"xChipID[LayerFilter&&ModuleFilter]") + .Define(yChipName,"yChipID[LayerFilter&&ModuleFilter]") + .Define(yHalfChipName,"yHalfChipID[LayerFilter&&ModuleFilter]") + .Define(xColumnName,"xColumnID[LayerFilter&&ModuleFilter]") + .Define(layerSizeName,[](RVecI lay){return lay.size();},layerSizeInput); + + + hHists1D[moduleTag+"/"+layerTag+"/"+xHistName] = d4.Histo1D({xHistName, xHistName+";x pixel column;Rate [Hz]", xRange, xMin, xMax }, xName, "eventWeight"); + + hHists1D[moduleTag+"/"+layerTag+"/"+yHistName] = d4.Histo1D({yHistName, yHistName+";y pixel column;Rate [Hz]", yRange, yMin, yMax }, yName, "eventWeight"); + + hHists2D[moduleTag+"/"+layerTag+"/"+xyHistName] = d4.Histo2D({xyHistName, xyHistName+";x pixel;y pixel;Rate [Hz]", xRange, xMin, xMax, yRange, yMin, yMax}, xName, yName, "eventWeight"); + + hHists1D[moduleTag+"/"+layerTag+"/"+sizeHistName] = d4.Histo1D({sizeHistName, sizeHistName+";hits per event", 100, 0, 100 }, layerSizeName ); + + // Plot number of hits per boardID for each layer + TString boardIDHistName = "hBoardIDRate"; + hHists1D[moduleTag+"/"+layerTag+"/"+boardIDHistName] = d4.Histo1D({boardIDHistName, boardIDHistName+";board ID;Rate [Hz]", 3, 0, 3}, boardName, "eventWeight"); + + // Plot number of hits per chipID for each layer + TString chipIDHistName = "hChipIDRate"; + hHists2D[moduleTag+"/"+layerTag+"/"+chipIDHistName] = d4.Histo2D({chipIDHistName, chipIDHistName+";x Chip;y Chip;Rate [Hz]", 6, 0, 6, 6, 0, 6}, xChipName, yChipName, "eventWeight"); + + // Plot number of hits per chipID for each layer + TString xColumnIDHistName = "hxColumnIDRate"; + hHists2D[moduleTag+"/"+layerTag+"/"+xColumnIDHistName] = d4.Histo2D({xColumnIDHistName, xColumnIDHistName+";x Column;y Half Chip;Rate [Hz]", 3*xChip, 0, 3.0*xChip, 12, 0, 12}, xColumnName, yHalfChipName, "eventWeight"); + + } + } + + return {hHists1D,hHists2D,hHists3D}; + +} diff --git a/benchmarks/LOWQ2/analysis/LOWQ2reconstruction.h b/benchmarks/LOWQ2/analysis/LOWQ2reconstruction.h new file mode 100644 index 00000000..01b0dfa2 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/LOWQ2reconstruction.h @@ -0,0 +1,150 @@ +#pragma once + +// Define alias +using RNode = ROOT::RDF::RNode; +using H1ResultPtr = ROOT::RDF::RResultPtr; +using H2ResultPtr = ROOT::RDF::RResultPtr; +using H3ResultPtr = ROOT::RDF::RResultPtr; +using RVecI = ROOT::VecOps::RVec; + +// Lazy execution methods +std::tuple,std::map,std::map> createReconstructionPlots( RNode d1 ){ + + int ePrimID = 2; + int pBeamID = 1; + int eBeamID = 0; + + std::map hHists1D; + std::map hHists2D; + std::map hHists3D; + + d1 = d1.Filter("LowQ2TrackParameters.size()>0") + .Define("eBeam", [eBeamID] + (std::vector h) + {return ROOT::Math::PxPyPzMVector (h[eBeamID].momentum.x,h[eBeamID].momentum.y,h[eBeamID].momentum.z,h[eBeamID].mass);}, + {"MCParticles"}) + .Define("pBeam", [pBeamID] + (std::vector h) + {return ROOT::Math::PxPyPzMVector (h[pBeamID].momentum.x,h[pBeamID].momentum.y,h[pBeamID].momentum.z,h[pBeamID].mass);}, + {"MCParticles"}) + .Define("primMom", [ePrimID] + (std::vector h) + {return ROOT::Math::PxPyPzMVector (h[ePrimID].momentum.x,h[ePrimID].momentum.y,h[ePrimID].momentum.z,h[ePrimID].mass);}, + {"MCParticles"}) + .Define("primE","primMom.E()") + .Define("primPx","primMom.Px()") + .Define("primPy","primMom.Py()") + .Define("primPz","primMom.Pz()") + .Define("primTheta","(TMath::Pi()-primMom.Theta())") + .Define("primTheta_mrad","1000*primTheta") + .Define("primEta","primMom.Eta()") + .Define("scatteredV","eBeam-primMom") + .Define("Q2","-scatteredV.M2()") + .Define("logQ2","log10(Q2)") + .Define("x","Q2/(2*scatteredV.Dot(pBeam))") + .Define("logx","log10(x)") + .Define("primPhi","primMom.Phi()") + .Define("primPhi_deg","primPhi * TMath::RadToDeg()") + .Define("W2","(pBeam+eBeam-primMom).M2()") + .Define("reconTheta","(double)(TMath::Pi()-LowQ2TrackParameters[0].theta)") + .Define("reconTheta_mrad","1000.0*reconTheta") + .Define("reconPhi","(double)LowQ2TrackParameters[0].phi") + .Define("reconPhi_deg","reconPhi * TMath::RadToDeg()") + .Define("reconP","(18./10.)*(double)(-1/(LowQ2TrackParameters[0].qOverP))") + .Define("reconMom",[](double p, double theta, double phi){return ROOT::Math::PxPyPzMVector(p*sin(theta)*cos(phi),p*sin(theta)*sin(phi),p*cos(theta),0.000511);},{"reconP","reconTheta","reconPhi"}) + .Define("reconE","reconMom.E()") + .Define("reconScatteredV","eBeam-reconMom") + .Define("reconQ2","-reconScatteredV.M2()") + .Define("reconlogQ2","log10(reconQ2)") + .Define("reconX","reconQ2/(2*reconScatteredV.Dot(pBeam))") + .Define("reconlogx","log10(reconX)") + .Define("reconW2","(pBeam+eBeam-reconMom).M2()") + .Define("reconPx","(double)reconMom.Px()") + .Define("reconPy","(double)reconMom.Py()") + .Define("reconPz","(double)reconMom.Pz()") + .Define("thetaRes","(reconTheta-primTheta)") + .Define("thetaRes_mrad","1000.0*thetaRes") + .Define("phiRes","reconPhi_deg-primPhi_deg") + .Define("ERes","(reconMom.E()-primE)/primE") + .Define("logQ2Res","reconlogQ2-logQ2") + .Define("XRes","reconX-x") + .Define("W2Res","reconW2-W2"); + + + + double thetaMin = 0; // mrad + double thetaMax = 12; // mrad + double phiMin = -180; // deg + double phiMax = 180; // deg + double pxMin = -0.2; // GeV + double pxMax = 0.2; // GeV + double pyMin = -0.2; // GeV + double pyMax = 0.2; // GeV + double pzMin = -18; // GeV + double pzMax = 0; // GeV + double eMin = 0; // GeV + double eMax = 18; // GeV + double logQ2Min = -8.5; // GeV^2 + double logQ2Max = 0; // GeV^2 + double logxMin = -11; + double logxMax = 0; + double w2Min = 0; // GeV^2 + double w2Max = 0.5; // GeV^2 + double thetaResMin = -5; // mrad + double thetaResMax = 5; // mrad + double phiResMin = -180; // deg + double phiResMax = 180; // deg + double eResMin = -0.1; // GeV + double eResMax = 0.1; // GeV + double q2ResMin = -0.01; // GeV^2 + double q2ResMax = 0.01; // GeV^2 + double xResMin = -0.1; + double xResMax = 0.1; + double w2ResMin = -0.1; // GeV^2 + double w2ResMax = 0.1; // GeV^2 + + int bins1D = 400; + int bins2D = 200; + int bins3D = 50; + int bins3DRes = 100; + + hHists2D["reconThetaVsPrimTheta"] = d1.Histo2D({"reconThetaVsPrimTheta","reconThetaVsPrimTheta;#theta_{prim} [mrad];#theta_{recon} [mrad]",bins2D,thetaMin,thetaMax,bins2D,thetaMin,thetaMax},"primTheta_mrad","reconTheta_mrad"); + hHists2D["reconPhiVsPrimPhi"] = d1.Histo2D({"reconPhiVsPrimPhi","reconPhiVsPrimPhi;#phi_{prim} [deg];#phi_{recon} [deg]",bins2D,phiMin,phiMax,bins2D,phiMin,phiMax},"primPhi_deg","reconPhi_deg"); + + hHists2D["reconPxVsPrimPx"] = d1.Histo2D({"reconPxVsPrimPx","reconPxVsPrimPx;p_{x,prim} [GeV];p_{x,recon} [GeV]",bins2D,pxMin,pxMax,bins2D,pxMin,pxMax},"primPx","reconPx"); + hHists2D["reconPyVsPrimPy"] = d1.Histo2D({"reconPyVsPrimPy","reconPyVsPrimPy;p_{y,prim} [GeV];p_{y,recon} [GeV]",bins2D,pyMin,pyMax,bins2D,pyMin,pyMax},"primPy","reconPy"); + hHists2D["reconPzVsPrimPz"] = d1.Histo2D({"reconPzVsPrimPz","reconPzVsPrimPz;p_{z,prim} [GeV];p_{z,recon} [GeV]",bins2D,pzMin,pzMax,bins2D,pzMin,pzMax},"primPz","reconPz"); + hHists2D["reconEVsPrimE"] = d1.Histo2D({"reconEVsPrimE","reconEVsPrimE;E_{prim} [GeV];E_{recon} [GeV]",bins2D,eMin,eMax,bins2D,eMin,eMax},"primE","reconE"); + hHists2D["reconQ2VsPrimQ2"] = d1.Histo2D({"reconQ2VsPrimQ2","reconQ2VsPrimQ2;log(Q^{2}_{prim}) [GeV^{2}];log(Q^{2}_{recon}) [GeV^{2}]",bins2D,logQ2Min,logQ2Max,bins2D,logQ2Min,logQ2Max},"logQ2","reconlogQ2"); + hHists2D["reconXVsPrimX"] = d1.Histo2D({"reconXVsPrimX","reconXVsPrimX;log(x_{prim});log(x_{recon})",bins2D,logxMin,logxMax,bins2D,logxMin,logxMax},"logx","reconlogx"); + hHists2D["reconW2VsPrimW2"] = d1.Histo2D({"reconW2VsPrimW2","reconW2VsPrimW2;W^{2}_{prim} [GeV^{2}];W^{2}_{recon} [GeV^{2}]",bins2D,w2Min,w2Max,bins2D,w2Min,w2Max},"W2","reconW2"); + + hHists1D["thetaRes"] = d1.Histo1D({"thetaRes","thetaRes;#theta_{recon}-#theta_{prim} [mrad]",bins1D,thetaResMin,thetaResMax},"thetaRes_mrad"); + hHists1D["phiRes"] = d1.Histo1D({"phiRes","phiRes;#phi_{recon}-#phi_{prim} [rad]",bins1D,phiResMin,phiResMax},"phiRes"); + hHists1D["ERes"] = d1.Histo1D({"ERes","ERes;E_{recon}-E_{prim} [GeV]",bins1D,eResMin,eResMax},"ERes"); + hHists1D["Q2Res"] = d1.Histo1D({"Q2Res","Q2Res;log(Q^{2}_{recon})-log(Q^{2}_{prim}) [GeV^{2}]",bins1D,q2ResMin,q2ResMax},"logQ2Res"); + hHists1D["XRes"] = d1.Histo1D({"XRes","XRes;log(x_{recon})-log(x_{prim})",bins1D,xResMin,xResMax},"XRes"); + hHists1D["W2Res"] = d1.Histo1D({"W2Res","W2Res;W^{2}_{recon}-W^{2}_{prim} [GeV^{2}]",bins1D,w2ResMin,w2ResMax},"W2Res"); + + hHists2D["thetaResVsE"] = d1.Histo2D({"thetaResVsE", "thetaResVsE;#theta_{recon}-#theta_{prim} [mrad];E_{prim} [GeV]", bins2D, thetaResMin, thetaResMax, bins2D, eMin, eMax}, "thetaRes_mrad", "primE"); + hHists2D["phiResVsE"] = d1.Histo2D({"phiResVsE", "phiResVsE;#phi_{recon}-#phi_{prim} [rad];E_{prim} [GeV]", bins2D, phiResMin, phiResMax, bins2D, eMin, eMax}, "phiRes", "primE"); + hHists2D["EResVsE"] = d1.Histo2D({"EResVsE", "EResVsE;E_{recon}-E_{prim} [GeV];E_{prim} [GeV]", bins2D, eResMin, eResMax, bins2D, eMin, eMax}, "ERes", "primE"); + hHists2D["Q2ResVsE"] = d1.Histo2D({"Q2ResVsE", "Q2ResVsE;log(Q^{2}_{recon})-log(Q^{2}_{prim}) [GeV^{2}];E_{prim} [GeV]", bins2D, q2ResMin, q2ResMax, bins2D, eMin, eMax}, "logQ2Res", "primE"); + + hHists2D["phiResVsTheta"] = d1.Histo2D({"phiResVsTheta", "phiResVsTheta;#phi_{recon}-#phi_{prim} [rad];#theta_{prim} [mrad]", bins2D, phiResMin, phiResMax, bins2D, thetaMin, thetaMax}, "phiRes", "primTheta_mrad"); + + //3D histograms to extract the resolution as a function of E and theta + hHists3D["thetaResVsETheta"] = d1.Histo3D({"thetaResVsETheta","thetaResVsETheta;E_{prim} [GeV];#theta_{prim} [mrad];#theta_{recon}-#theta_{prim} [mrad]",bins3D,eMin,eMax,bins3D,thetaMin,thetaMax,bins3DRes,thetaResMin,thetaResMax},"primE","primTheta_mrad","thetaRes_mrad"); + hHists3D["phiResVsETheta"] = d1.Histo3D({"phiResVsETheta","phiResVsETheta;E_{prim} [GeV];#theta_{prim} [mrad];#phi_{recon}-#phi_{prim} [rad]",bins3D,eMin,eMax,bins3D,thetaMin,thetaMax,bins3DRes,phiResMin,phiResMax},"primE","primTheta_mrad","phiRes"); + hHists3D["EResVsETheta"] = d1.Histo3D({"EResVsETheta","EResVsETheta;E_{prim} [GeV];#theta_{prim} [mrad];E_{recon}-E_{prim} [GeV]",bins3D,eMin,eMax,bins3D,thetaMin,thetaMax,bins3DRes,eResMin,eResMax},"primE","primTheta_mrad","ERes"); + + auto d2 = d1.Filter("reconTheta_mrad>1"); + + hHists1D["thetacut/phiRes"] = d2.Histo1D({"phiRes","phiRes;#phi_{recon}-#phi_{prim} [rad]",bins1D,phiResMin,phiResMax},"phiRes"); + hHists2D["thetacut/phiResVsE"] = d2.Histo2D({"phiResVsE", "phiResVsE;#phi_{recon}-#phi_{prim} [rad];E_{prim} [GeV]", bins2D, phiResMin, phiResMax, bins2D, eMin, eMax}, "phiRes", "primE"); + hHists2D["thetacut/reconPhiVsPrimPhi"] = d2.Histo2D({"reconPhiVsPrimPhi","reconPhiVsPrimPhi;#phi_{prim} [deg];#phi_{recon} [deg]",bins2D,phiMin,phiMax,bins2D,phiMin,phiMax},"primPhi_deg","reconPhi_deg"); + + return {hHists1D, hHists2D, hHists3D}; + +} + diff --git a/benchmarks/LOWQ2/analysis/PostprocessLOWQ2.C b/benchmarks/LOWQ2/analysis/PostprocessLOWQ2.C new file mode 100644 index 00000000..3bc5b380 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/PostprocessLOWQ2.C @@ -0,0 +1,533 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +//---------------------------------------------------------------------- +// Set global plot format variables +//---------------------------------------------------------------------- +void SetStyle() { + // Set global plot format variables + gStyle->SetOptStat(0); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetPadRightMargin(0.18); + gStyle->SetPadBottomMargin(0.12); + gStyle->SetTitleSize(0.04, "XYZ"); + gStyle->SetTitleOffset(4.0, "Z"); + +} + +//---------------------------------------------------------------------- +// Create acceptance plots +//---------------------------------------------------------------------- +TH1* AcceptancePlot(TDirectory* inputDir, TString ReconHistName, TString AllHistName, TString Tag="Quasi-Real") { + + // Read in the plots from the input file + TH1* ReconPlot = (TH1*)inputDir->Get(ReconHistName); + TH1* AllPlot = (TH1*)inputDir->Get(AllHistName); + + // Check plots exist + if (!ReconPlot || !AllPlot) { + std::cout << "Error: plots "<< ReconHistName <<" and/or "<< AllHistName <<" not found in input file" << std::endl; + return nullptr; + } + + //Divide the reconstructed plot by the all plot + ReconPlot->Divide(AllPlot); + + return ReconPlot; + +} + +//---------------------------------------------------------------------- +// Create rate plots +//---------------------------------------------------------------------- +TH2* RatePlot(TDirectory* inputDir, int Module, int Layer, TString Tag="Quasi-Real", TString inTag="AllHits") { + + TString histName = inTag+"/module"+std::to_string(Module)+"/layer"+std::to_string(Layer)+"/hxPixelyPixelRate"; + + // Read in the plots from the input file + TH2* RatePlot = (TH2*)inputDir->Get(histName); + // Check plots exist + if (!RatePlot) { + std::cout << "Error: plot "<< histName <<" not found in input file" << std::endl; + return nullptr; + } + + // Format the plot + int rebin = 32; + RatePlot->Rebin2D(rebin,rebin); + RatePlot->Scale(1.0/(rebin*rebin)); + TString title = "Tagger module "+std::to_string(Module)+", layer "+std::to_string(Layer)+" - Mean "+Tag+" rate per "+std::to_string(rebin)+"x"+std::to_string(rebin)+" pixel group"; + RatePlot->SetTitle(title); + RatePlot->GetXaxis()->SetTitle("x pixel"); + RatePlot->GetYaxis()->SetTitle("y pixel"); + RatePlot->GetZaxis()->SetTitle("Average Rate [Hz/55 #mum pixel]"); + + + RatePlot->SetStats(0); + + return RatePlot; + +} + +//---------------------------------------------------------------------- +// Create formatted acceptance plots on a canvas +//---------------------------------------------------------------------- +void FormatAcceptancePlots(TDirectory* inputDir, TFile* outputFile, TString Tag="Quasi-Real") { + + //---------------------------------------------------------------------- + // E-Q2 acceptance plot + //---------------------------------------------------------------------- + TString ReconEQ2Name = "Reconstructed-Track/hprimElogQ2"; + TString AllEQ2Name = "Generated/hprimElogQ2"; + + TH1* AcceptEQ2 = AcceptancePlot(inputDir,ReconEQ2Name,AllEQ2Name); + + TCanvas* canvasEQ2 = new TCanvas("AcceptanceEQ2Canvas", "AcceptanceEQ2Canvas", 1600, 1200); + + // Draw the plots on the canvas + canvasEQ2->cd(); + + TString title = "Acceptance as a function of scattered electron energy and reaction log_{10}(Q^{2})"; + AcceptEQ2->SetTitle(title); + AcceptEQ2->GetXaxis()->SetTitle("E_{e'} [GeV]"); + AcceptEQ2->GetYaxis()->SetTitle("log_{10}(Q^{2}) [GeV^{2}]"); + AcceptEQ2->GetZaxis()->SetTitle("Acceptance"); + + AcceptEQ2->SetStats(0); + AcceptEQ2->Draw("colz"); + + // Save the canvas to output file + outputFile->WriteTObject(canvasEQ2); + + // Clean up + delete canvasEQ2; + + //---------------------------------------------------------------------- + // E-Theta acceptance plot + //---------------------------------------------------------------------- + TString ReconEThetaName = "Reconstructed-Track/hprimEtheta"; + TString AllEThetaName = "Generated/hprimEtheta"; + + TH1* AcceptETheta = AcceptancePlot(inputDir,ReconEThetaName,AllEThetaName); + + TCanvas* canvasETheta = new TCanvas("AcceptanceEThetaCanvas", "AcceptanceEThetaCanvas", 1600, 1200); + + // Draw the plots on the canvas + canvasETheta->cd(); + + title = "Acceptance as a function of scattered electron energy and theta"; + AcceptETheta->SetTitle(title); + AcceptETheta->GetXaxis()->SetTitle("E_{e'} [GeV]"); + AcceptETheta->GetYaxis()->SetTitle("#theta_{e'} [mrad]"); + AcceptETheta->GetZaxis()->SetTitle("Acceptance"); + + AcceptETheta->SetStats(0); + AcceptETheta->Draw("colz"); + + // Save the canvas to output file + outputFile->WriteTObject(canvasETheta); + + // Clean up + delete canvasETheta; + + //---------------------------------------------------------------------- + // Integrated acceptance plot + //---------------------------------------------------------------------- + TString ReconIntegratedName = "IntegratedAcceptance"; + + // Read in the plots from the input file + TH1* IntegratedAcceptancePlot = (TH1*)inputDir->Get(ReconIntegratedName); + // Check plots exist + if (!IntegratedAcceptancePlot) { + std::cout << "Error: plot "<< ReconIntegratedName <<" not found in input file" << std::endl; + return; + } + + TCanvas* canvasIntegratedAcceptance = new TCanvas("IntegratedAcceptance", "IntegratedAcceptance", 1600, 1200); + + //Get xAxis title + TString xAxisTitle = IntegratedAcceptancePlot->GetXaxis()->GetTitle(); + //Break up xAxis title into words by "_" + TObjArray* xAxisTitleWords = xAxisTitle.Tokenize("_"); + //SetBinLabel for each bin + for (int i=1; i<=xAxisTitleWords->GetEntries(); i++) { + IntegratedAcceptancePlot->GetXaxis()->SetBinLabel(i, xAxisTitleWords->At(i-1)->GetName()); + } + + // Format the plot + IntegratedAcceptancePlot->SetTitle("Integrated acceptance"); + IntegratedAcceptancePlot->GetXaxis()->SetTitle(""); + IntegratedAcceptancePlot->GetYaxis()->SetTitle("Acceptance"); + + IntegratedAcceptancePlot->SetStats(0); + IntegratedAcceptancePlot->Draw(); + + + // Get the number of bins in the histogram + int nBins = IntegratedAcceptancePlot->GetNbinsX(); + + // Create a TText object to draw the text + TText t; + t.SetTextAlign(22); // Center the text + t.SetTextSize(0.02); // Set the text size + + // Loop over the bins + for (int i = 1; i <= nBins; ++i) { + // Get the bin content + double binContent = IntegratedAcceptancePlot->GetBinContent(i); + + // Get the bin center + double binCenter = IntegratedAcceptancePlot->GetBinCenter(i); + + // Draw the bin content at the bin center + t.DrawText(binCenter, binContent+0.02, Form("%.1f %s", binContent*100,"%")); + } + + // Update the canvas to show the changes + gPad->Update(); + + // Save the canvas to output file + outputFile->WriteTObject(canvasIntegratedAcceptance); + + // Clean up + delete canvasIntegratedAcceptance; + +} + +//---------------------------------------------------------------------- +// Create formatted rate plots on a canvas +//---------------------------------------------------------------------- +void FormatRatePlots(TDirectory* inputDir, TFile* outputFile, TString Tag="Quasi-Real") { + + + TCanvas* canvas = new TCanvas("RateCanvas", "RateCanvas", 3200, 1200); + canvas->Divide(2,1); + + // Read in the plots from the input file + TH2* RatePlot1_0 = RatePlot(inputDir,1,0,Tag); + TH2* RatePlot2_0 = RatePlot(inputDir,2,0,Tag); + + // Draw the plots on the canvas + canvas->cd(1); + gPad->SetLogz(); + RatePlot1_0->Draw("colz"); + + canvas->cd(2); + gPad->SetLogz(); + RatePlot2_0->Draw("colz"); + + // Save the canvas to output file + outputFile->WriteTObject(canvas); + + // Clean up + delete canvas; + + TCanvas* canvas2 = new TCanvas("RateCanvasOverlay", "RateCanvasOverlay", 3200, 1200); + canvas2->Divide(2,1); + + // Draw the plots on the canvas + canvas2->cd(1); + gPad->SetLogz(); + RatePlot1_0->Draw("colz"); + + //Add a grid ontop of the histogram showing 448x512 pixel chip boundaries + double xMin = RatePlot1_0->GetXaxis()->GetXmin(); + double xMax = RatePlot1_0->GetXaxis()->GetXmax(); + double yMin = RatePlot1_0->GetYaxis()->GetXmin()+512/4; + double yMax = RatePlot1_0->GetYaxis()->GetXmax()+512/4; + std::vector verticalLineWidths = {3,1,3,1,3,1,3}; + std::vector horizontalLineWidths = {3,1,2,1,2,1,2,1,2,1,2,1,3}; + std::vector horizontalLineStyles = {1,7,1,7,1,7,1,7,1,7,1,7,1}; + std::vector verticalLineColors = {kRed,kBlue,kRed,kBlue,kRed,kBlue,kRed}; + std::vector horizontalLineColors = {kRed,kBlue,kBlue,kBlue,kBlue,kBlue,kBlue,kBlue,kBlue,kBlue,kBlue,kBlue,kRed}; + //Vertical lines + for (int i=0; i<=6; i++) { + TLine* line = new TLine(xMin+i*448, yMin, xMin+i*448, yMax); + line->SetLineColor(verticalLineColors[i]); + line->SetLineWidth(verticalLineWidths[i]); + line->Draw(); + } + //Horizontal lines + for (int i=0; i<=12; i++) { + TLine* line = new TLine(xMin, yMin+i*256, xMax, yMin+i*256); + line->SetLineColor(horizontalLineColors[i]); + line->SetLineWidth(horizontalLineWidths[i]); + line->SetLineStyle(horizontalLineStyles[i]); + line->Draw(); + } + + gPad->Update(); + + canvas2->cd(2); + gPad->SetLogz(); + RatePlot2_0->Draw("colz"); + + //Add a grid on top of the histogram showing 448x512 pixel chip boundaries + xMin = RatePlot2_0->GetXaxis()->GetXmin(); + xMax = RatePlot2_0->GetXaxis()->GetXmax(); + yMin = RatePlot2_0->GetYaxis()->GetXmin()+512/4; + yMax = RatePlot2_0->GetYaxis()->GetXmax()+512/4; + //Vertical lines + for (int i = 0; i <= 6; i++) { + TLine* line = new TLine(xMin+i*448, yMin, xMin+i*448, yMax); + line->SetLineColor(verticalLineColors[i]); + line->SetLineWidth(verticalLineWidths[i]); + line->Draw(); + } + //Horizontal lines + for (int i = 0; i <= 12; i++) { + TLine* line = new TLine(xMin, yMin+i*256, xMax, yMin+i*256); + line->SetLineColor(horizontalLineColors[i]); + line->SetLineWidth(horizontalLineWidths[i]); + line->SetLineStyle(horizontalLineStyles[i]); + line->Draw(); + } + + gPad->Update(); + + // Save the canvas to output file + outputFile->WriteTObject(canvas); + + // Clean up + delete canvas; + + // Canvas showing primary and secondary hits + // Todo: Neaten up + TCanvas* canvas3 = new TCanvas("PrimarySecondary-RateCanvas", "PrimarySecondary-RateCanvas", 2400, 2400); + canvas3->Divide(2,2); + + TH2* RatePlotPrimary1_0 = RatePlot(inputDir,1,0,Tag,"PrimaryHits"); + TH2* RatePlotPrimary2_0 = RatePlot(inputDir,2,0,Tag,"PrimaryHits"); + + TH2* RatePlotSecondary1_0 = RatePlot(inputDir,1,0,Tag,"SecondaryHits"); + TH2* RatePlotSecondary2_0 = RatePlot(inputDir,2,0,Tag,"SecondaryHits"); + + // Draw the plots on the canvas + canvas3->cd(1); + gPad->SetLogz(); + RatePlotPrimary1_0->Draw("colz"); + + canvas3->cd(2); + gPad->SetLogz(); + RatePlotPrimary2_0->Draw("colz"); + + canvas3->cd(3); + gPad->SetLogz(); + RatePlotSecondary1_0->Draw("colz"); + + canvas3->cd(4); + gPad->SetLogz(); + RatePlotSecondary2_0->Draw("colz"); + + // Save the canvas to output file + outputFile->WriteTObject(canvas3); + + // Clean up + delete canvas3; + + +} + +//---------------------------------------------------------------------- +// Create formatted reconstruction plots on a canvas +//---------------------------------------------------------------------- +void FormatReconstructionPlots(TDirectory* inputDir, TFile* outputFile, TString Tag="Quasi-Real") { + + TString EHistName = "reconEVsPrimE"; + TString ThetaHistName = "reconThetaVsPrimTheta"; + TString PhiHistName = "thetacut/reconPhiVsPrimPhi"; + + TString EResName = "ERes"; + TString ThetaResName = "thetaRes"; + TString PhiResName = "thetacut/phiRes"; + + TCanvas* canvas = new TCanvas("ReconCanvas", "ReconCanvas", 2400, 1200); + + // Read in the plots from the input file + TH2* EPlot = (TH2*)inputDir->Get(EHistName); + TH2* ThetaPlot = (TH2*)inputDir->Get(ThetaHistName); + TH2* PhiPlot = (TH2*)inputDir->Get(PhiHistName); + + TH1* EResPlot = (TH1*)inputDir->Get(EResName); + TH1* ThetaResPlot = (TH1*)inputDir->Get(ThetaResName); + TH1* PhiResPlot = (TH1*)inputDir->Get(PhiResName); + + // Check plots exist + if (!ThetaPlot || !EPlot || !PhiPlot || !ThetaResPlot || !EResPlot || !PhiResPlot) { + std::cout << "Error: plots "<< ThetaHistName <<", "<< EHistName <<", "<< PhiHistName <<", "<< ThetaResName <<", "<< EResName <<", "<< PhiResName <<" not found in input file" << std::endl; + return; + } + + // Draw the plots on the canvas + canvas->Divide(3,2); + canvas->cd(1); + gPad->SetLogz(); + + EPlot->SetTitle("Reconstructed electron energy vs. Primary electron energy"); + EPlot->GetXaxis()->SetTitle("E_{prim} [GeV]"); + EPlot->GetYaxis()->SetTitle("E_{recon} [GeV]"); + EPlot->GetZaxis()->SetTitle("Counts"); + EPlot->Rebin2D(2,2); + + EPlot->SetStats(0); + EPlot->Draw("colz"); + + canvas->cd(2); + gPad->SetLogz(); + + ThetaPlot->SetTitle("Reconstructed #theta vs. Primary #theta"); + ThetaPlot->GetXaxis()->SetTitle("#theta_{prim} [mrad]"); + ThetaPlot->GetYaxis()->SetTitle("#theta_{recon} [mrad]"); + ThetaPlot->GetZaxis()->SetTitle("Counts"); + ThetaPlot->Rebin2D(2,2); + + ThetaPlot->SetStats(0); + ThetaPlot->Draw("colz"); + + canvas->cd(3); + gPad->SetLogz(); + + PhiPlot->SetTitle("Reconstructed #varphi vs. Primary #varphi (#theta>1mrad)"); + PhiPlot->GetXaxis()->SetTitle("#phi_{prim} [deg]"); + PhiPlot->GetYaxis()->SetTitle("#phi_{recon} [deg]"); + PhiPlot->GetZaxis()->SetTitle("Counts"); + PhiPlot->Rebin2D(2,2); + + PhiPlot->SetStats(0); + PhiPlot->Draw("colz"); + + canvas->cd(4); + + EResPlot->SetTitle("Electron energy resolution"); + EResPlot->GetXaxis()->SetTitle("(E_{recon} - E_{prim}) / E_{prim}"); + EResPlot->GetYaxis()->SetTitle("Counts"); + + EResPlot->SetStats(0); + EResPlot->Draw(); + + // Write fitted Gaussian standard deviation in pad + //Fit Gaussian to histogram maximum bin +- 10 bins + int maxBin = EResPlot->GetMaximumBin(); + int fitMin = maxBin-10; + int fitMax = maxBin+10; + EResPlot->Fit("gaus", "Q", "", EResPlot->GetBinCenter(fitMin), EResPlot->GetBinCenter(fitMax)); + double EResStdDev = EResPlot->GetFunction("gaus")->GetParameter(2); + TLatex* latex = new TLatex(); + latex->SetTextSize(0.03); + latex->DrawLatexNDC(0.2, 0.8, Form("#sigma_{E} = %.2f %s", EResStdDev*100, "%")); + + // Remove fitted Gaussian from histogram + EResPlot->GetListOfFunctions()->Remove(EResPlot->GetFunction("gaus")); + + canvas->cd(5); + + ThetaResPlot->SetTitle("Theta resolution"); + ThetaResPlot->GetXaxis()->SetTitle("#theta_{recon} - #theta_{prim} [mrad]"); + ThetaResPlot->GetYaxis()->SetTitle("Counts"); + + ThetaResPlot->SetStats(0); + ThetaResPlot->Draw(); + + // Write fitted Gaussian standard deviation in pad + //Fit Gaussian to histogram maximum bin +- 10 bins + maxBin = ThetaResPlot->GetMaximumBin(); + fitMin = maxBin-10; + fitMax = maxBin+10; + ThetaResPlot->Fit("gaus", "Q", "", ThetaResPlot->GetBinCenter(fitMin), ThetaResPlot->GetBinCenter(fitMax)); + double ThetaResStdDev = ThetaResPlot->GetFunction("gaus")->GetParameter(2); + latex->DrawLatexNDC(0.2, 0.8, Form("#sigma_{#theta} = %.2f mrad", ThetaResStdDev)); + + // Remove fitted Gaussian from histogram + ThetaResPlot->GetListOfFunctions()->Remove(ThetaResPlot->GetFunction("gaus")); + + canvas->cd(6); + + PhiResPlot->SetTitle("Phi resolution (#theta>1mrad)"); + PhiResPlot->GetXaxis()->SetTitle("#phi_{recon} - #phi_{prim} [deg]"); + PhiResPlot->GetYaxis()->SetTitle("Counts"); + + PhiResPlot->SetStats(0); + PhiResPlot->Draw(); + + // Write fitted Gaussian standard deviation in pad + //Fit Gaussian to histogram maximum bin +- 10 bins + maxBin = PhiResPlot->GetMaximumBin(); + fitMin = maxBin-10; + fitMax = maxBin+10; + PhiResPlot->Fit("gaus", "Q", "", PhiResPlot->GetBinCenter(fitMin), PhiResPlot->GetBinCenter(fitMax)); + double PhiResStdDev = PhiResPlot->GetFunction("gaus")->GetParameter(2); + latex->DrawLatexNDC(0.2, 0.8, Form("#sigma_{#phi} = %.1f deg", PhiResStdDev)); + + // Remove fitted Gaussian from histogram + PhiResPlot->GetListOfFunctions()->Remove(PhiResPlot->GetFunction("gaus")); + + // Save the canvas to output file + outputFile->WriteTObject(canvas); + + // Clean up + delete canvas; + +} + +//---------------------------------------------------------------------- +// This function is called by the benchmarking script +//---------------------------------------------------------------------- +void PostprocessLOWQ2(TString inName, TString outName, TString Tag) { + + SetStyle(); + + // Open the input file + TFile* inputFile = TFile::Open(inName); + if (!inputFile) { + std::cout << "Error opening input file:" << inName << std::endl; + return; + } + + // Check the directory LOWQ2 exists and cd into it + TDirectory* dir = inputFile->GetDirectory("LOWQ2"); + if (!dir) { + std::cout << "Error: directory LOWQ2 not found in input file" << std::endl; + return; + } + + // Open the output file + TFile* outputFile = TFile::Open(outName, "RECREATE"); + + // Format the plots + + //Check if AcceptanceDistributions directory exists + TDirectory* dirA = dir->GetDirectory("Acceptance"); + if (dirA) { + FormatAcceptancePlots(dirA, outputFile,Tag); + } + + //Check if SimDistributions directory exists + TDirectory* dirS = dir->GetDirectory("Rates"); + if (dirS) { + FormatRatePlots(dirS, outputFile,Tag); + } + + //Check if ReconstructedDistributions directory exists + TDirectory* dirR = dir->GetDirectory("Reconstruction"); + if (dirR) { + FormatReconstructionPlots(dirR, outputFile,Tag); + } + + inputFile->Close(); + outputFile->Close(); + +} + +//---------------------------------------------------------------------- +// Main function to create canvases +//---------------------------------------------------------------------- +// void PostprocessLOWQ2() { +// // Postprocess("plots/LOWQ2QRRecon2.root", "plots/LOWQ2QR_FormattedPlots.root", "Quasi-Real"); +// Postprocess("plots/LOWQ2BremsRecon3.root", "plots/LOWQ2Brems_FormattedPlots3.root", "Bremsstrahlung"); +// } \ No newline at end of file diff --git a/benchmarks/LOWQ2/analysis/RunLOWQ2.C b/benchmarks/LOWQ2/analysis/RunLOWQ2.C new file mode 100644 index 00000000..e73e9849 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/RunLOWQ2.C @@ -0,0 +1,40 @@ +#include "LOWQ2Benchmarks.h" +#include + +void RunLOWQ2( std::string inputFileName = "Brems_input.root", + std::string outputFileName = "plots/LOWQ2QRRecon3.root", + std::string compactName = "/opt/detector/epic-nightly/share/epic/epic.xml", + bool inputIsTimeBased = false, // true if the event sample is time-based, false if it is event-based + double timeWindow = 10.15*1e-9, //[s] + double eventCrossSection = 0.0551, // [mb] + double luminosity = 1e34 // [cm^-2 s^-1] + ) { + + //Set implicit multi-threading + ROOT::EnableImplicitMT(); + + // Output script running conditions + std::cout << "Running LOWQ2 benchmarks with the following parameters:" << std::endl; + std::cout << " - input file: " << inputFileName << std::endl; + std::cout << " - output file: " << outputFileName << std::endl; + std::cout << " - xml file: " << compactName << std::endl; + std::cout << " - input is time-based: " << inputIsTimeBased << std::endl; + std::cout << " - time window: " << timeWindow << " s" << std::endl; + std::cout << " - event cross section: " << eventCrossSection << " mb" << std::endl; + std::cout << " - luminosity: " << luminosity << " cm^-2 s^-1" << std::endl; + + dd4hep::Detector& detector = dd4hep::Detector::getInstance(); + detector.fromCompact(compactName); + + // eventCrossSectionBrems = 171.3; [mb] + // eventCrossSectionQR = 0.0551; [mb] + + double eventRate = luminosity * eventCrossSection * 1e-27; // [Hz] + + if(inputIsTimeBased){ + eventRate = 1.0 / timeWindow; // [Hz] + } + + LOWQ2Benchmarks(inputFileName,outputFileName,detector,eventRate); + +} \ No newline at end of file diff --git a/benchmarks/LOWQ2/analysis/Snakefile b/benchmarks/LOWQ2/analysis/Snakefile new file mode 100644 index 00000000..27a45002 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/Snakefile @@ -0,0 +1,93 @@ +# Possible tags for local/remote running +tag_params = { + "QR": {"timeBased": 0, "timeWindow": 0, "eventCrossSection": 0.0551, "tagName": "Quasi-Real"}, + "Brems": {"timeBased": 1, "timeWindow": 10.15*1e-9, "eventCrossSection": 0, "tagName": "Bremsstrahlung"}, + "pythia6": {"timeBased": 0, "timeWindow": 0, "eventCrossSection": 0.054, "tagName": "Deep Inelastic Scattering"}, + # Add more mappings here if needed +} + +configfile: "../local_config.yml" + +EVENT_EXTENSION = ".ab.hepmc3.tree.root" +SIM_EXTENSION = ".edm4hep.root" +RECO_EXTENSION = ".eicrecon.tree.edm4eic.root" + + +REMOTE_EVENTS_SERVER = "root://dtn-eic.jlab.org/" +REMOTE_EVENTS_DIRECTORY = "/work/eic2/EPIC/EVGEN/SIDIS/pythia6-eic/1.0.0/10x100/q2_0to1/" + +FILE_BASE = "pythia_ep_noradcor_10x100_q2_0.000000001_1.0_run" + +BEAM_ENERGY = "10" +XML_FILE = "/opt/detector/epic-24.05.0/share/epic/epic.xml", + +# Function to check if the remote file exists +def remote_file_exists(server,url): + try: + subprocess.check_output(['xrdfs', server, 'stat', url]) + return url + except subprocess.CalledProcessError: + return None + +rule run_simulation_remote: + params: + XML=XML_FILE, + input=lambda wildcards: remote_file_exists(REMOTE_EVENTS_SERVER,REMOTE_EVENTS_DIRECTORY+FILE_BASE+wildcards.index+EVENT_EXTENSION), + output: + config["SIM_DIRECTORY"]+FILE_BASE+"{index}_tagger"+SIM_EXTENSION, + shell: """ + npsim \ + --inputFiles {params.input} \ + --outputFile {output[0]} \ + --compactFile {params.XML} \ + --runType run \ + --numberOfEvents 100 \ + --physics.list FTFP_BERT \ + --field.eps_min 5e-06 \ + --field.eps_max 1e-04 \ + --physics.rangecut 50 \ + """ + +rule run_reconstruction: + params: + XML=XML_FILE, + beam_energy=BEAM_ENERGY, + collections="TaggerTrackerProjectedTracks,MCScatteredElectrons,MCParticles,EventHeader", + input: + expand( + config["SIM_DIRECTORY"]+FILE_BASE+"{index}_tagger"+SIM_EXTENSION, + index=range(1,4), + ), + output: + config["RECO_IN_DIRECTORY"]+FILE_BASE+".tagger_recon"+RECO_EXTENSION, + shell: """ + eicrecon {input} -Pdd4hep:xml_files={params.XML} -Ppodio:output_include_collections={params.collections} -Ppodio:output_file={output} -PLOWQ2:LowQ2Trajectories:electron_beamE={params.beam_energy} + """ + +rule run_benchmarks: + input: + config["RECO_IN_DIRECTORY"]+FILE_BASE+".tagger_recon"+RECO_EXTENSION + params: + xmlName=XML_FILE, + timeBased = lambda wildcards: tag_params[wildcards.tag]["timeBased"], + timeWindow = lambda wildcards: tag_params[wildcards.tag]["timeWindow"], + eventCrossSection = lambda wildcards: tag_params[wildcards.tag]["eventCrossSection"] + + output: + "{dir}/LOWQ2{tag}_Plots.root" + shell: + """ + root -l -b -q 'RunLOWQ2.C++("{input}", "{output}", "{params.xmlName}", {params.timeBased}, {params.timeWindow}, {params.eventCrossSection})' + """ + +rule run_plots: + input: + "{dir}/LOWQ2{tag}_Plots.root" + params: + tagName = lambda wildcards: tag_params[wildcards.tag]["tagName"], + output: + "{dir}/LOWQ2{tag}_FormattedPlots.root" + shell: + """ + root -l -b -q 'PostprocessLOWQ2.C("{input}", "{output}","{params.tagName}")' + """ diff --git a/benchmarks/LOWQ2/analysis/config.yml b/benchmarks/LOWQ2/analysis/config.yml new file mode 100644 index 00000000..4d45650a --- /dev/null +++ b/benchmarks/LOWQ2/analysis/config.yml @@ -0,0 +1,5 @@ +bench:LOWQ2_events: + extends: .det_benchmark + stage: benchmarks + script: + - snakemake --cores 8 ${RESULTS_DIRECTORY}/LOWQ2_FormattedPlots.edm4hep.root --configfile ../remote_config.yml diff --git a/benchmarks/LOWQ2/analysis/functors.h b/benchmarks/LOWQ2/analysis/functors.h new file mode 100644 index 00000000..e978ff89 --- /dev/null +++ b/benchmarks/LOWQ2/analysis/functors.h @@ -0,0 +1,36 @@ +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHit.h" +#include "edm4hep/SimCalorimeterHit.h" + +using RVecHits = ROOT::VecOps::RVec; + +//----------------------------------------------------------------------------------------- +// Grab Component functor +//----------------------------------------------------------------------------------------- +struct getSubID{ + getSubID(std::string cname, dd4hep::Detector& det, std::string rname = "TaggerTrackerHits") : componentName(cname), detector(det), readoutName(rname){} + + ROOT::VecOps::RVec operator()(const RVecHits& evt) { + auto decoder = detector.readout(readoutName).idSpec().decoder(); + auto indexID = decoder->index(componentName); + ROOT::VecOps::RVec result; + for(const auto& h: evt) { + result.push_back(decoder->get(h.cellID,indexID)); + } + return result; + }; + + void SetComponent(std::string cname){ + componentName = cname; + } + void SetReadout(std::string rname){ + readoutName = rname; + } + + private: + std::string componentName; + dd4hep::Detector& detector; + std::string readoutName; +}; diff --git a/benchmarks/LOWQ2/local_config.yml b/benchmarks/LOWQ2/local_config.yml new file mode 100644 index 00000000..437975b8 --- /dev/null +++ b/benchmarks/LOWQ2/local_config.yml @@ -0,0 +1,3 @@ +SIM_DIRECTORY: "/scratch/EIC/SimOut/S3in/" +RECO_IN_DIRECTORY: "/scratch/EIC/ReconOut/S3in/" +MODEL_DIRECTORY: "/scratch/EIC/LowQ2Model/" \ No newline at end of file diff --git a/benchmarks/LOWQ2/reconstruction_training/.gitignore b/benchmarks/LOWQ2/reconstruction_training/.gitignore new file mode 100644 index 00000000..c96503c7 --- /dev/null +++ b/benchmarks/LOWQ2/reconstruction_training/.gitignore @@ -0,0 +1,5 @@ +plots* +calibrations* +fieldmaps* +gdml* +.snakemake* \ No newline at end of file diff --git a/benchmarks/LOWQ2/remote_config.yml b/benchmarks/LOWQ2/remote_config.yml new file mode 100644 index 00000000..b793dbac --- /dev/null +++ b/benchmarks/LOWQ2/remote_config.yml @@ -0,0 +1,3 @@ +SIM_DIRECTORY: "LowQ2_Sim/" +RECO_IN_DIRECTORY: "LowQ2_ReconIn/" +MODEL_DIRECTORY: "LowQ2_Model/" \ No newline at end of file diff --git a/benchmarks/LOWQ2/signal_training/.gitignore b/benchmarks/LOWQ2/signal_training/.gitignore new file mode 100644 index 00000000..c96503c7 --- /dev/null +++ b/benchmarks/LOWQ2/signal_training/.gitignore @@ -0,0 +1,5 @@ +plots* +calibrations* +fieldmaps* +gdml* +.snakemake* \ No newline at end of file