Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
  • Loading branch information
zihlmann committed Feb 18, 2022
2 parents a59147b + f232b36 commit 704f7ea
Show file tree
Hide file tree
Showing 47 changed files with 6,442 additions and 517 deletions.
2 changes: 1 addition & 1 deletion AnalysisHowTo/ThrownTopology/plotTopology.C
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ void plotTopology() {
TH1F *hThrownTopologies = (TH1F*)f->Get("hThrownTopologies");
hThrownTopologies->GetXaxis()->LabelsOption(">"); // order by most common topology
hThrownTopologies->GetXaxis()->SetRangeUser(0, 10); // only plot first 20 topologies
hThrownTopologies->Scale(100./hThrownTopologies->GetEntries()); // turn histogram into percentage
hThrownTopologies->Scale(100./hThrownTopologies->Integral()); // turn histogram into percentage
hThrownTopologies->GetYaxis()->SetTitle("Thrown Topology %");
hThrownTopologies->Draw("htext");

Expand Down
Binary file modified FlattenForFSRoot/Documentation/GlueXFSRootFormat.pdf
Binary file not shown.
10 changes: 10 additions & 0 deletions FlattenForFSRoot/Documentation/GlueXFSRootFormat.tex
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,16 @@ \section{Track Information}
TkNDFP(n): number of degrees of freedom for the track fit
\end{verbatim}

If PID information is requested, the following variables are added to the tree:
\begin{verbatim}
TkTOFBetaP(n): beta from the TOF (using kinematic fit if available)
TkTOFChi2P(n): chi2 from the TOF (using kinematic fit if available)
TkDEDXCDCP(n): DEDX from the CDC
TkDEDXFDCP(n): DEDX from the FDC
TkDEDXChi2P(n): chi2 for the DEDX measurement
TkDEDXNDFP(n): NDF for the DEDX measurement
\end{verbatim}


\section{Shower Information}
\label{sec:nt:shower}
Expand Down
46 changes: 44 additions & 2 deletions FlattenForFSRoot/flatten.cc
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ int main(int argc, char** argv){
cout << " -numUnusedNeutrals [optional cut (<= cut)] (default: -1 (no cut))" << endl;
cout << " -numNeutralHypos [optional cut (<= cut)] (default: -1 (no cut))" << endl;
cout << " -usePolarization [get polarization angle from RCDB? 0 or 1] (default: 0)" << endl;
cout << " -addPID [include PID info in the output tree? 0 or 1] (default: 0)" << endl;
cout << " -mcChecks [check for baryon number violation, etc.," << endl;
cout << " when parsing truth information? 0 or 1] (default: 1)" << endl;
cout << " -safe [check array sizes? 0 or 1] (default: 1)" << endl;
Expand Down Expand Up @@ -121,6 +122,7 @@ int main(int argc, char** argv){
int gNumNeutralHyposCut = -1;
bool gSafe = true;
bool gUsePolarization = false;
bool gAddPID = false;
bool gMCChecks = true;
int gPrint = 0;
{
Expand All @@ -130,7 +132,7 @@ int main(int argc, char** argv){
if ((argi == "-in")||(argi == "-out")||(argi == "-mc")||(argi == "-mctag")
||(argi == "-chi2")||(argi == "-shQuality")||(argi == "-massWindows")
||(argi == "-numUnusedTracks")||(argi == "-usePolarization")||(argi == "-numUnusedNeutrals")
||(argi == "-mcChecks")
||(argi == "-mcChecks")||(argi == "-addPID")
||(argi == "-numNeutralHypos")||(argi == "-safe")||(argi == "-print")){
flag = argi;
continue;
Expand All @@ -147,6 +149,7 @@ int main(int argc, char** argv){
if (flag == "-numUnusedNeutrals"){ gNumUnusedNeutralsCut = atoi(argi); }
if (flag == "-numNeutralHypos"){ gNumNeutralHyposCut = atoi(argi); }
if (flag == "-usePolarization"){ if (argi == "1") gUsePolarization = true; }
if (flag == "-addPID"){ if (argi == "1") gAddPID = true; }
if (flag == "-mcChecks"){ if (argi == "0") gMCChecks = false; }
if (flag == "-safe"){ if (argi == "0") gSafe = false; }
if (flag == "-print"){ gPrint = atoi(argi); }
Expand Down Expand Up @@ -182,6 +185,7 @@ int main(int argc, char** argv){
cout << " numUnusedNeutrals cut: " << gNumUnusedNeutralsCut << endl;
cout << " numNeutralHypos cut: " << gNumNeutralHyposCut << endl;
cout << " use polarization? " << gUsePolarization << endl;
cout << " use PID? " << gAddPID << endl;
cout << " MC checks? " << gMCChecks << endl;
cout << " safe mode? " << gSafe << endl;
cout << endl;
Expand Down Expand Up @@ -642,6 +646,14 @@ int main(int argc, char** argv){
if (gUseParticles) gInTree->SetBranchAddress("ChargedHypo__ChiSq_Tracking", inChargedHypo__ChiSq_Tracking);
UInt_t inChargedHypo__NDF_Tracking[MAXTRACKS] = {};
if (gUseParticles) gInTree->SetBranchAddress("ChargedHypo__NDF_Tracking", inChargedHypo__NDF_Tracking);
UInt_t inChargedHypo__NDF_DCdEdx[MAXTRACKS] = {};
if (gUseParticles && gAddPID) gInTree->SetBranchAddress("ChargedHypo__NDF_DCdEdx", inChargedHypo__NDF_DCdEdx);
Float_t inChargedHypo__ChiSq_DCdEdx[MAXTRACKS] = {};
if (gUseParticles && gAddPID) gInTree->SetBranchAddress("ChargedHypo__ChiSq_DCdEdx", inChargedHypo__ChiSq_DCdEdx);
Float_t inChargedHypo__dEdx_CDC[MAXTRACKS] = {};
if (gUseParticles && gAddPID) gInTree->SetBranchAddress("ChargedHypo__dEdx_CDC", inChargedHypo__dEdx_CDC);
Float_t inChargedHypo__dEdx_FDC[MAXTRACKS] = {};
if (gUseParticles && gAddPID) gInTree->SetBranchAddress("ChargedHypo__dEdx_FDC", inChargedHypo__dEdx_FDC);


// *** Neutral Particle Hypotheses (indexed using <particleName>__NeutralIndex) ***
Expand Down Expand Up @@ -693,6 +705,8 @@ int main(int argc, char** argv){

TClonesArray *inP4_KinFit[MAXPARTICLES] = {};
Int_t inChargedIndex[MAXPARTICLES][MAXCOMBOS] = {};
Float_t inBeta_Timing[MAXPARTICLES][MAXCOMBOS] = {};
Float_t inChiSq_Timing[MAXPARTICLES][MAXCOMBOS] = {};
Int_t inNeutralIndex[MAXPARTICLES][MAXCOMBOS] = {};
TClonesArray *inX4[MAXPARTICLES] = {};
Float_t inPathLengthSigma[MAXPARTICLES][MAXCOMBOS] = {};
Expand All @@ -711,6 +725,14 @@ int main(int argc, char** argv){
if (gUseParticles && gUseKinFit) gInTree->SetBranchAddress(var_P4_KinFit,&(inP4_KinFit[pIndex]));
TString var_ChargedIndex(name); var_ChargedIndex += "__ChargedIndex";
if (gUseParticles) gInTree->SetBranchAddress(var_ChargedIndex,inChargedIndex[pIndex]);
TString var_Beta_Timing(name);
if (gUseParticles && !gUseKinFit) var_Beta_Timing += "__Beta_Timing_Measured";
if (gUseParticles && gUseKinFit) var_Beta_Timing += "__Beta_Timing_KinFit";
if (gUseParticles && gAddPID) gInTree->SetBranchAddress(var_Beta_Timing,inBeta_Timing[pIndex]);
TString var_ChiSq_Timing(name);
if (gUseParticles && !gUseKinFit) var_ChiSq_Timing += "__ChiSq_Timing_Measured";
if (gUseParticles && gUseKinFit) var_ChiSq_Timing += "__ChiSq_Timing_KinFit";
if (gUseParticles && gAddPID) gInTree->SetBranchAddress(var_ChiSq_Timing,inChiSq_Timing[pIndex]);
}

// *** Combo Neutrals ***
Expand Down Expand Up @@ -813,6 +835,9 @@ int main(int argc, char** argv){
double outMCPx[MAXPARTICLES]={}, outMCPy[MAXPARTICLES]={}, outMCPz[MAXPARTICLES]={}, outMCEn[MAXPARTICLES]={};
double outVeeL[MAXPARTICLES]={}, outVeeLSigma[MAXPARTICLES]={};
double outTkChi2[MAXPARTICLES]={}, outTkNDF[MAXPARTICLES]={};
double outTkDEDXChi2[MAXPARTICLES]={}, outTkDEDXNDF[MAXPARTICLES]={};
double outTkDEDXCDC[MAXPARTICLES]={}, outTkDEDXFDC[MAXPARTICLES]={};
double outTkTOFBeta[MAXPARTICLES]={}, outTkTOFChi2[MAXPARTICLES]={};
double outShQuality[MAXPARTICLES]={};
{
for (unsigned int im = 0; im < gOrderedParticleNames.size(); im++){
Expand All @@ -835,6 +860,14 @@ int main(int argc, char** argv){
TString vTkNDF("TkNDFP"); vTkNDF += fsIndex; gOutTree->Branch(vTkNDF, &outTkNDF [pIndex]);
TString vTkChi2("TkChi2P"); vTkChi2 += fsIndex; gOutTree->Branch(vTkChi2,&outTkChi2[pIndex]);
}
if (gAddPID && GlueXParticleClass(name) == "Charged"){
TString vTkTOFBeta ("TkTOFBetaP"); vTkTOFBeta += fsIndex; gOutTree->Branch(vTkTOFBeta, &outTkTOFBeta[pIndex]);
TString vTkTOFChi2 ("TkTOFChi2P"); vTkTOFChi2 += fsIndex; gOutTree->Branch(vTkTOFChi2, &outTkTOFChi2[pIndex]);
TString vTkDEDXChi2("TkDEDXChi2P"); vTkDEDXChi2 += fsIndex; gOutTree->Branch(vTkDEDXChi2,&outTkDEDXChi2[pIndex]);
TString vTkDEDXNDF("TkDEDXNDFP"); vTkDEDXNDF += fsIndex; gOutTree->Branch(vTkDEDXNDF, &outTkDEDXNDF [pIndex]);
TString vTkDEDXCDC("TkDEDXCDCP"); vTkDEDXCDC += fsIndex; gOutTree->Branch(vTkDEDXCDC, &outTkDEDXCDC [pIndex]);
TString vTkDEDXFDC("TkDEDXFDCP"); vTkDEDXFDC += fsIndex; gOutTree->Branch(vTkDEDXFDC, &outTkDEDXFDC [pIndex]);
}
if (GlueXParticleClass(name) == "Neutral"){
TString vQual("ShQualityP"); vQual += fsIndex; gOutTree->Branch(vQual, &outShQuality[pIndex]);
}
Expand Down Expand Up @@ -1134,6 +1167,14 @@ int main(int argc, char** argv){
outREn[pIndex] = p4->E();
outTkNDF [pIndex] = inChargedHypo__NDF_Tracking [(inChargedIndex[pIndex][ic])];
outTkChi2[pIndex] = inChargedHypo__ChiSq_Tracking[(inChargedIndex[pIndex][ic])];
if (gAddPID){
outTkTOFBeta[pIndex] = inBeta_Timing[pIndex][ic];
outTkTOFChi2[pIndex] = inChiSq_Timing[pIndex][ic];
outTkDEDXChi2[pIndex] = inChargedHypo__ChiSq_DCdEdx[(inChargedIndex[pIndex][ic])];
outTkDEDXNDF [pIndex] = inChargedHypo__NDF_DCdEdx [(inChargedIndex[pIndex][ic])];
outTkDEDXCDC [pIndex] = inChargedHypo__dEdx_CDC [(inChargedIndex[pIndex][ic])];
outTkDEDXFDC [pIndex] = inChargedHypo__dEdx_FDC [(inChargedIndex[pIndex][ic])];
}
}
if (gUseMCParticles && outMCSignal > 0.1){
p4 = (TLorentzVector*)inThrown__P4->At(tIndex);
Expand Down Expand Up @@ -1897,7 +1938,8 @@ map<TString, vector<TString> > GlueXDecayProductMap(int fsCode1, int fsCode2){
tmp = "Photon"; tmp += (pNumber++); names.push_back(tmp); }
if (name.Contains("KShort")){ tmp = "PiPlus"; tmp += (pNumber++); names.push_back(tmp);
tmp = "PiMinus"; tmp += (pNumber++); names.push_back(tmp); }
if (name.Contains("Lambda")){ tmp = "Proton"; tmp += (pNumber++); names.push_back(tmp);
if (name.Contains("Lambda")&&!name.Contains("AntiLambda"))
{ tmp = "Proton"; tmp += (pNumber++); names.push_back(tmp);
tmp = "PiMinus"; tmp += (pNumber++); names.push_back(tmp); }
if (name.Contains("AntiLambda")){ tmp = "AntiProton"; tmp += (pNumber++); names.push_back(tmp);
tmp = "PiPlus"; tmp += (pNumber++); names.push_back(tmp); }
Expand Down
14 changes: 14 additions & 0 deletions PWA_scripts/benchmark/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
# Introduction
This directory contains scripts for submitting batch jobs with MPI (and GPU). They are meant to determine how your fit speed improves (or not) by adding additioinal resources
* submit.py -- submits MPI jobs with various # of cores
* submitGPU.py -- submits MPI+GPU jobs with various # of GPUs
* plotBenchmark.C -- plots fit speed for the benchmark jobs

# Required user modifications
* In submit.py and submitGPU.py you should replace the MyEnv, MyConfig and MyOutDir variables with your own environment setup script, AmpTools fit configuration and output directory location
* You can change the MyCPU or MyGPU list to contain different amounts of cores for your benchmark if you want to test with more or fewer cores/GPUs

# Notes:
* These fits require using the MPI compiled version of AmpTools and halld_sim, see https://halldweb.jlab.org/wiki/index.php/HOWTO_use_AmpTools_on_the_JLab_farm_with_MPI for more details
* Only run the GPU version of the fitter if your fit utilizes a GPU accelerated amplitude and you've compiled that amplitude with the GPU (CUDA) libraries on one of the sciml nodes, see https://halldweb.jlab.org/wiki/index.php/HOWTO_use_AmpTools_on_the_JLab_farm_GPUs for more details
* Some of these default benchmarks require many CPUs or GPUs and may take some time for those nodes to become available on the ifarm/sciml nodes, so be patient.
123 changes: 123 additions & 0 deletions PWA_scripts/benchmark/plot_benchmark.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <bits/stdc++.h>
#include <string>
#include <sstream>

void plot_benchmark(TString dir = "./") {

gStyle->SetOptStat(0);

// initialize list of nCores to plot
vector<int> numThreadsCPU = {1,2,4,8,16,32,64,96,128};
int numTestCPU = numThreadsCPU.size();

// for GPU fits, only add if desired
vector<int> numThreadsGPUT4 = {1,2,3,4,6,8,10,12};
vector<int> numThreadsGPURTX = {};

// names of directories containing benchmark results
vector<TString> types = {"cpu"};
vector<TGraphErrors*> grBenchmarkScan;
if(numThreadsGPUT4.size() > 0) types.push_back("gpuT4");
if(numThreadsGPURTX.size() > 0) types.push_back("gpuTitanRTX");

TH1F *hBenchmarkScan = new TH1F("hBenchmarkScan","; Number of GPUs or CPUs; Fit speed (Likelihood function call rate [Hz])", 200, 0, 200);
double maxRate = 0;

// loop over different achitecture types to plot results
for(int itype=0; itype<types.size(); itype++) {
vector<int> numThreads = numThreadsCPU;
if(types[itype] == "gpuT4") numThreads = numThreadsGPUT4;
if(types[itype] == "gpuTitanRTX") numThreads = numThreadsGPURTX;
grBenchmarkScan.push_back(new TGraphErrors(numThreads.size()));

// loop over number of threads in test
for(int ithread=0; ithread<numThreads.size(); ithread++) {

int nThreads = numThreads[ithread];
string spath = Form("%s/%s%03d/log/fit.out", dir.Data(), types[itype].Data(), nThreads);
//cout << spath << endl;

std::string read_line;
ifstream file(spath);
double parValue = 0;
double parAvg = 0;
vector<double> parSq;
int nValues = 0;
while (std::getline(file, read_line)) {

TString line = read_line;
if(line.Contains("time ")) {
line.ReplaceAll("average time per function call: ","");
line.ReplaceAll(" ms.","");
parValue = 1./(atof(line)/1000);
parAvg += parValue;
parSq.push_back(parValue*parValue);
nValues++;
}
else continue;

}

if(nValues > 0) {
parAvg /= float(nValues);
double parRms = 0;
for(uint ip=0; ip<parSq.size(); ip++)
parRms += (parSq.at(ip) + parAvg*parAvg - 2*sqrt(parSq.at(ip))*parAvg);
parRms /= float(nValues);
parRms = sqrt(parRms);
if(parAvg > maxRate) maxRate = parAvg;
//cout<<parAvg<<" "<<parRms<<endl;
if(parRms < 1e-9) parRms = 0.01;
grBenchmarkScan[itype]->SetPoint(ithread, nThreads, parAvg);
grBenchmarkScan[itype]->SetPointError(ithread, 0, parRms);
}
}
}

TCanvas *cc = new TCanvas("cc","cc",800,400);
auto legend = new TLegend(0.47,0.17,0.9,0.42);

hBenchmarkScan->SetMaximum(maxRate*2.5);
hBenchmarkScan->SetMinimum(0.1);
hBenchmarkScan->Draw();
vector<TF1*> fit;
for(int itype=0; itype<types.size(); itype++) {
grBenchmarkScan[itype]->SetMarkerStyle(20);
grBenchmarkScan[itype]->SetMarkerColor(kBlack+itype);
grBenchmarkScan[itype]->Draw("same pl");

fit.push_back(new TF1(types[itype],"pol1",1,200));
fit[itype]->FixParameter(0,0);
grBenchmarkScan[itype]->Fit(fit[itype],"N","",0.5,2);
fit[itype]->SetLineColor(kBlack+itype); fit[itype]->SetLineStyle(kDashed);
fit[itype]->Draw("same");

if(itype==0)
legend->AddEntry(grBenchmarkScan[0],"ifarm19 CPU (2 thread/core)","pl");
if(types[itype] == "gpuT4")
legend->AddEntry(grBenchmarkScan[itype],"sciml21 T4 GPU","pl");
if(types[itype] == "gpuTitanRTX")
legend->AddEntry(grBenchmarkScan[itype],"sciml19 Titan RTX GPU","pl");
}

gPad->SetLeftMargin(0.09);
gPad->SetBottomMargin(0.15);
gPad->SetTopMargin(0.05);
gPad->SetRightMargin(0.05);
gPad->SetLogx(); gPad->SetLogy();
gPad->SetGridy(); gPad->SetGridx();
hBenchmarkScan->GetXaxis()->SetTitleSize(0.05);
hBenchmarkScan->GetYaxis()->SetTitleSize(0.05);
hBenchmarkScan->GetXaxis()->SetTitleOffset(1.3);
hBenchmarkScan->GetYaxis()->SetTitleOffset(0.8);

legend->SetFillColor(0);
legend->Draw();

cc->Print("benchmark.png");

return;
}
72 changes: 72 additions & 0 deletions PWA_scripts/benchmark/submit.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
#!/usr/bin/env python

import sys
import os
import subprocess
import math
import pwd
from optparse import OptionParser

########################################################## MAIN ##########################################################
def main(argv):

# SLURM INFO (see options at https://scicomp.jlab.org/scicomp/slurmJob/slurmInfo)
PARTITION = "ifarm"
CONSTRAINT = "farm19"
TIMELIMIT = "24:00:00" # Max walltime
MyCPUs = [1,2,4,8,16,32,64,96,128,192] # List of CPU cores to use in benchmark fits

# User provided environment, fit configuration and options
MyEnv = "/work/halld2/home/jrsteven/analysisGluexI/builds/setup_gluex_scanParam.csh"
MyConfig = "/work/halld2/home/jrsteven/forBenchmark/benchmark.cfg"
MyMPIOpt = "--mca btl_openib_allow_ib 1"
MyFitOpt = "-m 100000 -r 5"
MyOutDir = "/volatile/halld/home/" + pwd.getpwuid( os.getuid() )[0] + "/benchmark/"

# LOOP OVER # OF CORES FOR BENCHMARK
for nCores in MyCPUs:
# nodes used in fit (for every 64 CPUs allow an additional node)
nNodes = nCores/64 + 1

# create output directories
MyRunningDir = MyOutDir + "cpu%03d" % nCores
MyLogOutDir = MyRunningDir + "/log"
if not os.path.exists(MyOutDir):
os.makedirs(MyOutDir)
if not os.path.exists(MyRunningDir):
os.makedirs(MyRunningDir)
if not os.path.exists(MyLogOutDir):
os.makedirs(MyLogOutDir)

# create slurm submission script
slurmOut = open("tempSlurm.txt",'w')
slurmOut.write("#!/bin/csh \n")
slurmOut.write("#SBATCH --nodes=%d \n" % nNodes)
slurmOut.write("#SBATCH --partition=%s \n" % PARTITION)
slurmOut.write("#SBATCH --constraint=%s \n" % CONSTRAINT)
slurmOut.write("#SBATCH --cpus-per-task=1 \n")
slurmOut.write("#SBATCH --ntasks-per-core=1 \n")
slurmOut.write("#SBATCH --threads-per-core=1 \n")
slurmOut.write("#SBATCH --mem=%dGB \n" % nCores) # 1 GB per core
slurmOut.write("#SBATCH --time=%s \n" % TIMELIMIT)
slurmOut.write("#SBATCH --ntasks=%d \n" % (nCores+1))

slurmOut.write("#SBATCH --chdir=%s \n" % MyRunningDir)
slurmOut.write("#SBATCH --error=%s/fit.err \n" % (MyLogOutDir))
slurmOut.write("#SBATCH --output=%s/fit.out \n" % (MyLogOutDir))
slurmOut.write("#SBATCH --job-name=benchfit_%03d \n\n\n" % nCores)

# commands to execute during job
slurmOut.write("pwd \n")
slurmOut.write("source %s \n" % MyEnv)
slurmOut.write("mpirun %s fitMPI -c %s %s \n" % (MyMPIOpt, MyConfig, MyFitOpt))
slurmOut.close()

# submit individual job
print("Submitting %d core job on %d nodes" % (nCores, nNodes))
subprocess.call(["sbatch", "tempSlurm.txt"])
os.remove("tempSlurm.txt")


if __name__ == "__main__":
main(sys.argv[1:])
Loading

0 comments on commit 704f7ea

Please sign in to comment.