| ID |
Date |
Author |
Subject |
Project |
|
18
|
Tue Dec 12 17:38:36 2017 |
Brian Clark | Data Analysis in R from day w/ Brian Connolly | Analysis | On Nov 28 2017, Brian Connolly came and visited and taught us how to do basic data analysis in R.
He in particular showed us how to do a Linear Discriminant Analysis (LDA) and Principal Component Analysis (PCA).
Attached are three of the data files Carl Pfendner prepared for us to analyze (ARA data, including simulation, force triggers, and RF triggers).
Also attached is some R code that shows how to set up the LDA and the PCA and how to plot their result. You are meant to run each line of the code in r by hand (this is not a functioning R script I don't think).
Go here (https://www.r-project.org/) to learn how to download R. You will also probably have to download GGFortify. To do that, open an r session, and type "install.packages('ggfortify')". |
| Attachment 1: pca_and_lda.R
|
Notes
#First we need to read in the data
dataSim <- read.table("output.txt.Simulation",sep=" ",header=TRUE)
dataRF <- read.table("output.txt.RFTriggers",sep=" ",header=TRUE)
dataPulser <- read.table("output.txt.CalPulsers",sep=" ",header=TRUE)
#Now we combine them into one data object
data <- rbind(dataRF,dataSim)
#Now we actually have to specify that we
#want to use the first 10 columns
data <- data[,1:10]
#make a principal component analysis object
pca <- prcomp(data,center=TRUE,scale=TRUE,retx=TRUE)
#Load a plotting library
library(ggfortify)
#then can plot
autplot(prcomp(data))
#or, we can plot like this
#to get some color data points
labels<-c(rep(0,nrow(dataSim)),rep(1,nrow(dataRF)))
plot(pca$x[,1],pca$x[,2],col=labels+1)
#we can also do an LDA analysis
#we need to import the MASS library
library(MASS)
#and now we make the lda object
lda(data,grouping=labels)
#didn't write down any plotting unfortunately.
|
| Attachment 2: data.zip
|
|
27
|
Mon Oct 1 19:06:59 2018 |
Brian Clark | Code to Compute Effective Volumes in AraSim | Analysis | Here is some C++ code and an associated makefile to find effective volumes from AraSim output files.
It computes error bars on the effective volumes using the relevant AraSim function.
Compile like "make -f veff.mk"
Run like "./veff thrown_radius thrown_depth AraOut.1.root AraOut.2.root...." |
| Attachment 1: veff.cc
|
#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <time.h>
#include <cstdio>
#include <cstdlib>
#include <unistd.h>
#include <cstring>
#include <unistd.h>
#include "TTreeIndex.h"
#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "Event.h"
#include "Detector.h"
#include "Report.h"
#include "Vector.h"
#include "counting.hh"
using namespace std;
//Get effective volumes with error bars
class EarthModel;
int main(int argc, char **argv)
{
gStyle->SetOptStat(111111);
gStyle->SetOptDate(1);
if(argc<4){
cout<<"Not enough arguments! Abort run."<<endl;
cout<<"Run like: ./veff volradius voldepth AraOut1.root AraOut2.rot ..."<<endl;
return -1;
}
double volradius = atof(argv[1]);
double voldepth = atof(argv[2]);
TChain *AraTree = new TChain("AraTree");
TChain *AraTree2 = new TChain("AraTree2");
for(int i=3; i<argc;i++){
AraTree->Add(string(argv[i]).c_str());
AraTree2->Add(string(argv[i]).c_str());
}
Report *report = 0;
Event *event=0;
AraTree2->SetBranchAddress("report",&report);
AraTree2->SetBranchAddress("event",&event);
int totnthrown = AraTree2->GetEntries();
cout << "Total number of events: " << totnthrown << endl;
int NBINS=10;
double eventsfound_binned[NBINS];
for(int i=0; i<NBINS; i++) eventsfound_binned[i]=0.;
double totweight=0;
for(int iEvt2=0; iEvt2<totnthrown; iEvt2++){
AraTree2->GetEntry(iEvt2);
if(report->stations[0].Global_Pass<=0) continue;
double weight = event->Nu_Interaction[0].weight;
if(weight > 1.0) continue;
totweight += weight;
int index_weights = Counting::findWeightBin(log10(weight));
if(index_weights<NBINS) eventsfound_binned[index_weights]++;
}
double error_minus=0.;
double error_plus=0.;
Counting::findErrorOnSumWeights(eventsfound_binned,error_plus,error_minus);
double vtot = TMath::Pi() * double(volradius) * double(volradius) * double(voldepth) / 1.e9; //answer in km^3
double veff = vtot * totweight / double(totnthrown) * 4. * TMath::Pi(); //answer in km^3 sr
double veff_p = vtot * (error_plus) / double(totnthrown) * 4. * TMath::Pi(); //answer in km^3 sr
double veff_m = vtot * (error_minus) / double(totnthrown) * 4. * TMath::Pi(); //answer in km^3 sr
printf("volthrown: %.6f \n totweight: %.6f + %.6f - %.6f \n Veff: %.6f + %.6f - %.6f \n",
vtot,
totweight, error_plus, error_minus,
veff, veff_p, veff_m
);
return 0;
}
|
| Attachment 2: veff.mk
|
#############################################################################
##
##Changes:
##line 54 - OBJS = .... add filename.o .... del oldfilename.o
##line 55 - CCFILE = .... add filename.cc .... del oldfilename.cc
##line 58 - PROGRAMS = filename
##line 62 - filename : $(OBJS)
##
##############################################################################
include StandardDefinitions.mk
#Site Specific Flags
ifeq ($(strip $(BOOST_ROOT)),)
BOOST_ROOT = /usr/local/include
endif
SYSINCLUDES = -I/usr/include -I$(BOOST_ROOT)
SYSLIBS = -L/usr/lib
DLLSUF = ${DllSuf}
OBJSUF = ${ObjSuf}
SRCSUF = ${SrcSuf}
CXX = g++
#Generic and Site Specific Flags
CXXFLAGS += $(INC_ARA_UTIL) $(SYSINCLUDES)
LDFLAGS += -g $(LD_ARA_UTIL) -I$(BOOST_ROOT) $(ROOTLDFLAGS) -L.
# copy from ray_solver_makefile (removed -lAra part)
# added for Fortran to C++
LIBS = $(ROOTLIBS) -lMinuit $(SYSLIBS)
GLIBS = $(ROOTGLIBS) $(SYSLIBS)
LIB_DIR = ./lib
INC_DIR = ./include
#ROOT_LIBRARY = libAra.${DLLSUF}
OBJS = Vector.o EarthModel.o IceModel.o Trigger.o Ray.o Tools.o Efficiencies.o Event.o Detector.o Position.o Spectra.o RayTrace.o RayTrace_IceModels.o signal.o secondaries.o Settings.o Primaries.o counting.o RaySolver.o Report.o eventSimDict.o veff.o
CCFILE = Vector.cc EarthModel.cc IceModel.cc Trigger.cc Ray.cc Tools.cc Efficiencies.cc Event.cc Detector.cc Spectra.cc Position.cc RayTrace.cc signal.cc secondaries.cc RayTrace_IceModels.cc Settings.cc Primaries.cc counting.cc RaySolver.cc Report.cc veff.cc
CLASS_HEADERS = Trigger.h Detector.h Settings.h Spectra.h IceModel.h Primaries.h Report.h Event.h secondaries.hh #need to add headers which added to Tree Branch
PROGRAMS = veff
all : $(PROGRAMS)
veff : $(OBJS)
$(LD) $(OBJS) $(LDFLAGS) $(LIBS) -o $(PROGRAMS)
@echo "done."
#The library
$(ROOT_LIBRARY) : $(LIB_OBJS)
@echo "Linking $@ ..."
ifeq ($(PLATFORM),macosx)
# We need to make both the .dylib and the .so
$(LD) $(SOFLAGS)$@ $(LDFLAGS) $(G77LDFLAGS) $^ $(OutPutOpt) $@
ifneq ($(subst $(MACOSX_MINOR),,1234),1234)
ifeq ($(MACOSX_MINOR),4)
ln -sf $@ $(subst .$(DllSuf),.so,$@)
else
$(LD) -dynamiclib -undefined $(UNDEFOPT) $(LDFLAGS) $(G77LDFLAGS) $^ \
$(OutPutOpt) $(subst .$(DllSuf),.so,$@)
endif
endif
else
$(LD) $(SOFLAGS) $(LDFLAGS) $(G77LDFLAGS) $(LIBS) $(LIB_OBJS) -o $@
endif
##-bundle
#%.$(OBJSUF) : %.$(SRCSUF)
# @echo "<**Compiling**> "$<
# $(CXX) $(CXXFLAGS) -c $< -o $@
%.$(OBJSUF) : %.C
@echo "<**Compiling**> "$<
$(CXX) $(CXXFLAGS) $ -c $< -o $@
%.$(OBJSUF) : %.cc
@echo "<**Compiling**> "$<
$(CXX) $(CXXFLAGS) $ -c $< -o $@
# added for fortran code compiling
%.$(OBJSUF) : %.f
@echo "<**Compiling**> "$<
$(G77) -c $<
eventSimDict.C: $(CLASS_HEADERS)
@echo "Generating dictionary ..."
@ rm -f *Dict*
rootcint $@ -c ${INC_ARA_UTIL} $(CLASS_HEADERS) ${ARA_ROOT_HEADERS} LinkDef.h
clean:
@rm -f *Dict*
@rm -f *.${OBJSUF}
@rm -f $(LIBRARY)
@rm -f $(ROOT_LIBRARY)
@rm -f $(subst .$(DLLSUF),.so,$(ROOT_LIBRARY))
@rm -f $(TEST)
#############################################################################
|
|
28
|
Fri Oct 26 18:08:43 2018 |
Jorge Torres | Analyzing effective volumes | Analysis | Attaching some scripts that help processing the effective volumes. This is an extension of what Brian Clark did in a previous post (http://radiorm.physics.ohio-state.edu/elog/How-To/27)
There are 4 files attached:
- veff_aeff2.C and veff_aeff2.mk. veff_aeff2.C produces Veff_des$1.txt ($1 can be A or B or C). This file contains the following columns: energy, veff, veff_error, veff1 (PA), veff2 (LPDA), veff3 (bicone), respectively. However, the energies are not sorted.
-veff.sh: this bash executable runs veff_aeff2.C for all (that's what the "*" in the executable is for) the root output files, for a given design (A, B, C). You need to modify the location of your output files, though. Run like "./veff.sh A", which will execute veff_aeff2.C and produce the veff text files. Do the same for B or C.
-make_plot.py: takes Veff_des$1.txt, sorts energies out, plots the effective volumes vs. energies, and produces a csv file containing the veffs (just for the sake of copying and pastting on the spreadsheets). Run like "pyhton make_plot.py".
|
| Attachment 1: veff.sh
|
nohup ./veff_aeff2 3000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_16.5.txt.run* &
nohup ./veff_aeff2 3000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_17.txt.run* &
nohup ./veff_aeff2 3000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_17.5.txt.run* &
nohup ./veff_aeff2 5000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_18.txt.run* &
nohup ./veff_aeff2 5000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_18.5.txt.run* &
nohup ./veff_aeff2 7000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_19.txt.run* &
nohup ./veff_aeff2 7000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_19.5.txt.run* &
nohup ./veff_aeff2 7000 3000 $1 /users/PAS0654/osu8354/outputs_signal_noise/des"$1"/AraOut.des"$1"_20.txt.run* &
|
| Attachment 2: veff_aeff2.C
|
//////////////////////////////////////
//To calculate various veff and aeff
//Run as ./veff_aeff2 $RADIUS $DEPTH $FILE
//For example ./veff_aeff2 8000 3000 /data/user/ypan/bin/AraSim/trunk/outputs/AraOut.setup_single_station_2_energy_20.run0.root
//
/////////////////////////////////////
#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <time.h>
#include "TTreeIndex.h"
#include "TChain.h"
#include "TH1.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TRandom.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TTree.h"
#include "TLegend.h"
#include "TLine.h"
#include "TROOT.h"
#include "TPostScript.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TText.h"
#include "TProfile.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TMath.h"
#include <unistd.h>
#include "TVector3.h"
#include "TRotation.h"
#include "TSpline.h"
//#include "TObject.h"
#include "Tools.h"
#include "Constants.h"
#include "Vector.h"
#include "Position.h"
#include "EarthModel.h"
#include "IceModel.h"
#include "Efficiencies.h"
#include "Spectra.h"
#include "Event.h"
#include "Trigger.h"
#include "Detector.h"
#include "Settings.h"
#include "counting.hh"
#include "Primaries.h"
#include "signal.hh"
#include "secondaries.hh"
#include "Ray.h"
#include "RaySolver.h"
#include "Report.h"
using namespace std;
class EarthModel;
int main(int argc, char **argv){
//string readfile;
//readfile = string(argv[1]);
//readfile = "/data/user/ypan/bin/AraSim/branches/new_geom/outputs/AraOut20.0.root";
int ifile = 0;
double totweightsq;
double totweight;
int totnthrown;
int typetotnthrown[12];
double tottrigeff;
double sigma[12];
double typetotweight[12];
double typetotweightsq[12];
double totsigmaweight;
double totsigmaweightsq;
double volradius;
double voldepth;
const int nstrings = 9;
const int nantennas = 1;
double veff1 = 0.0;
double veff2 = 0.0;
double veff3 = 0.0;
double weight1 = 0.0;
double weight2 = 0.0;
double weight3 = 0.0;
double veffT[6], vefferrT[6], aeffT[6], aefferrT[6];
double veffF[3], vefferrF[3], aeffF[3], aefferrF[3];
double veffNu[2], vefferrNu[2], aeffNu[2], aefferrNu[2];
double veff, vefferr, aeff, aefferr, aeff2;
double pnu;
Detector *detector = 0;
//Settings *settings = 0;
//IceModel *icemodel = 0;
Event *event = 0;
Report *report = 0;
cout<<"construct detector"<<endl;
//TFile *AraFile=new TFile(readfile.c_str());
//TFile *AraFile=new TFile((outputdir+"/AraOut.root").c_str());
TChain *AraTree = new TChain("AraTree");
TChain *AraTree2 = new TChain("AraTree2");
TChain *eventTree = new TChain("eventTree");
//AraTree->SetBranchAddress("detector",&detector);
//AraTree->SetBranchAddress("settings",&settings);
//AraTree->SetBranchAddress("icemodel",&icemodel);
cout << "trees set" << endl;
for(ifile = 3; ifile < (argc - 1); ifile++){
AraTree->Add(string(argv[ifile + 1]).c_str());
AraTree2->Add(string(argv[ifile + 1]).c_str());
eventTree->Add(string(argv[ifile + 1]).c_str());
}
AraTree2->SetBranchAddress("event",&event);
AraTree2->SetBranchAddress("report",&report);
cout<<"branch detector"<<endl;
for(int i=0; i<12; i++) {
typetotweight[i] = 0.0;
typetotweightsq[i] = 0.0;
typetotnthrown[i] = 0;
}
totweightsq = 0.0;
totweight = 0.0;
totsigmaweight = 0.0;
totsigmaweightsq = 0.0;
totnthrown = AraTree2->GetEntries();
cout << "Total number of events: " << totnthrown << endl;
//totnthrown = settings->NNU;
//volradius = settings->POSNU_RADIUS;
volradius = atof(argv[1]);
voldepth = atof(argv[2]);
AraTree2->GetEntry(0);
pnu = event->pnu;
cout << "Energy " << pnu << endl;
for(int iEvt2=0; iEvt2<totnthrown; iEvt2++) {
AraTree2->GetEntry(iEvt2);
double sigm = event->Nu_Interaction[0].sigma;
int iflavor = (event->nuflavorint)-1;
int inu = event->nu_nubar;
int icurr = event->Nu_Interaction[0].currentint;
sigma[inu+2*icurr+4*iflavor] = sigm;
typetotnthrown[inu+2*icurr+4*iflavor]++;
if( (iEvt2 % 10000 ) == 0 ) cout << "*";
if(report->stations[0].Global_Pass<=0) continue;
double weight = event->Nu_Interaction[0].weight;
if(weight > 1.0){
cout << weight << "; " << iEvt2 << endl;
continue;
}
// cout << weight << endl;
totweightsq += pow(weight,2);
totweight += weight;
typetotweight[inu+2*icurr+4*iflavor] += weight;
typetotweightsq[inu+2*icurr+4*iflavor] += pow(weight,2);
totsigmaweight += weight*sigm;
totsigmaweightsq += pow(weight*sigm,2);
int trig1 = 0;
int trig2 = 0;
int trig3 = 0;
for (int i = 0; i < nstrings; i++){
if (i == 0 && report->stations[0].strings[i].antennas[0].Trig_Pass > 0) trig1++;
if (i > 0 && i < 5 && report->stations[0].strings[i].antennas[0].Trig_Pass > 0) trig2++;
if (i > 4 && i < 9 && report->stations[0].strings[i].antennas[0].Trig_Pass > 0) trig3++;
}
if ( trig1 > 0)//phase array
weight1 += event->Nu_Interaction[0].weight;
if ( trig2 > 1)//lpda
weight2 += event->Nu_Interaction[0].weight;
if (trig3 > 3)//bicone
weight3 += event->Nu_Interaction[0].weight;
}
tottrigeff = totweight / double(totnthrown);
double nnucleon = 5.54e29;
double vtot = PI * double(volradius) * double(volradius) * double(voldepth) / 1e9;
veff = vtot * totweight / double(totnthrown) * 4.0 * PI;
//vefferr = sqrt(SQ(sqrt(double(totnthrown))/double(totnthrown))+SQ(sqrt(totweightsq)/totweight));
vefferr = sqrt(totweightsq) / totweight * veff;
aeff = vtot * (1e3) * nnucleon * totsigmaweight / double(totnthrown);
//aefferr = sqrt(SQ(sqrt(double(totnthrown))/double(totnthrown))+SQ(sqrt(totsigmaweightsq)/totsigmaweight));
//aefferr = sqrt(SQ(sqrt(double(totnthrown))/double(totnthrown))+SQ(sqrt(totweightsq)/totweight));
aefferr = sqrt(totweightsq) / totweight * aeff;
double sigmaave = 0.0;
for(int iflavor=0; iflavor<3; iflavor++) {
double flavorweight = 0.0;
double flavorweightsq = 0.0;
double flavorsigmaave = 0.0;
int flavortotthrown = 0;
double temptotweightnu[2] = {0};
double tempsignu[2] = {0};
double temptotweight = 0.0;
for(int inu=0; inu<2; inu++) {
double tempsig = 0.0;
double tempweight = 0.0;
for(int icurr=0; icurr<2; icurr++) {
tempsig += sigma[inu+2*icurr+4*iflavor];
tempsignu[inu] += sigma[inu+2*icurr+4*iflavor];
tempweight += typetotweight[inu+2*icurr+4*iflavor];
flavorweight += typetotweight[inu+2*icurr+4*iflavor];
flavorweightsq += typetotweightsq[inu+2*icurr+4*iflavor];
temptotweight += typetotweight[inu+2*icurr+4*iflavor];
temptotweightnu[inu] += typetotweight[inu+2*icurr+4*iflavor];
flavortotthrown += typetotnthrown[inu+2*icurr+4*iflavor];
}
//printf("Temp Sigma: "); cout << tempsig << "\n";
sigmaave += tempsig*(tempweight/totweight);
}
flavorsigmaave += tempsignu[0]*(temptotweightnu[0]/temptotweight)+tempsignu[1]*(temptotweightnu[1]/temptotweight);
veffF[iflavor] = vtot*flavorweight/double(flavortotthrown);
vefferrF[iflavor] = sqrt(flavorweightsq)/flavorweight;
//printf("Volume: %.9f*%.9f/%.9f \n",vtot,flavorweight,double(totnthrown));
aeffF[iflavor] = veffF[iflavor]*(1e3)*nnucleon*flavorsigmaave;
aefferrF[iflavor] = sqrt(flavorweightsq)/flavorweight;
}
for(int inu=0; inu<2; inu++) {
double tempsig = 0.0;
double tempweight = 0.0;
double tempweightsq = 0.0;
int nutotthrown = 0;
for(int iflavor=0; iflavor<3; iflavor++) {
for(int icurr=0; icurr<2; icurr++) {
tempweight += typetotweight[inu+2*icurr+4*iflavor];
tempweightsq += typetotweightsq[inu+2*icurr+4*iflavor];
nutotthrown += typetotnthrown[inu+2*icurr+4*iflavor];
}
}
tempsig += sigma[inu+2*0+4*0];
tempsig += sigma[inu+2*1+4*0];
veffNu[inu] = vtot*tempweight/double(nutotthrown);
vefferrNu[inu] = sqrt(tempweightsq)/tempweight;
aeffNu[inu] = veffNu[inu]*(1e3)*nnucleon*tempsig;
aefferrNu[inu] = sqrt(tempweightsq)/tempweight;
}
double totalveff = 0.0;
double totalaeff = 0.0;
for(int inu=0; inu<2; inu++) {
for(int iflavor=0; iflavor<3; iflavor++) {
int typetotthrown = 0;
for(int icurr=0; icurr<2; icurr++) {
typetotthrown += typetotnthrown[inu+2*icurr+4*iflavor];
}
totalveff += veffT[iflavor+3*inu]*(double(typetotthrown)/double(totnthrown));
totalaeff += aeffT[iflavor+3*inu]*(double(typetotthrown)/double(totnthrown));
}
}
aeff2 = veff*(1e3)*nnucleon*sigmaave;
aeff = aeff2;
veff1 = weight1 / totnthrown * vtot * 4.0 * PI;
veff2 = weight2 / totnthrown * vtot * 4.0 * PI;
veff3 = weight3 / totnthrown * vtot * 4.0 * PI;
printf("\nvolthrown: %.6f; totweight: %.6f; Veff: %.6f +- %.6f\n", vtot, totweight, veff, vefferr);
printf("veff1: %.3f; veff2: %.3f; veff3: %.3f\n", veff1, veff2, veff3);
//string des = string(argv[4]);
char buf[100];
std::ostringstream stringStream;
stringStream << string(argv[3]);
std::string copyOfStr = stringStream.str();
snprintf(buf, sizeof(buf), "Veff_des%s.txt", copyOfStr.c_str());
FILE *fout = fopen(buf, "a+");
fprintf(fout, "%e, %.6f, %.6f, %.3f, %.3f, %.3f \n", pnu, veff, vefferr, veff1, veff2, veff3);
fclose(fout);
return 0;
}
|
| Attachment 3: veff_aeff2.mk
|
#############################################################################
## Makefile -- New Version of my Makefile that works on both linux
## and mac os x
## Ryan Nichol <rjn@hep.ucl.ac.uk>
##############################################################################
##############################################################################
##############################################################################
##
##This file was copied from M.readGeom and altered for my use 14 May 2014
##Khalida Hendricks.
##
##Modified by Brian Clark for use on CosTheta_NuTraject on 20 February 2015
##
##Changes:
##line 54 - OBJS = .... add filename.o .... del oldfilename.o
##line 55 - CCFILE = .... add filename.cc .... del oldfilename.cc
##line 58 - PROGRAMS = filename
##line 62 - filename : $(OBJS)
##
##
##############################################################################
##############################################################################
##############################################################################
include StandardDefinitions.mk
#Site Specific Flags
ifeq ($(strip $(BOOST_ROOT)),)
BOOST_ROOT = /usr/local/include
endif
SYSINCLUDES = -I/usr/include -I$(BOOST_ROOT)
SYSLIBS = -L/usr/lib
DLLSUF = ${DllSuf}
OBJSUF = ${ObjSuf}
SRCSUF = ${SrcSuf}
CXX = g++
#Generic and Site Specific Flags
CXXFLAGS += $(INC_ARA_UTIL) $(SYSINCLUDES)
LDFLAGS += -g $(LD_ARA_UTIL) -I$(BOOST_ROOT) $(ROOTLDFLAGS) -L.
# copy from ray_solver_makefile (removed -lAra part)
# added for Fortran to C++
LIBS = $(ROOTLIBS) -lMinuit $(SYSLIBS)
GLIBS = $(ROOTGLIBS) $(SYSLIBS)
LIB_DIR = ./lib
INC_DIR = ./include
#ROOT_LIBRARY = libAra.${DLLSUF}
OBJS = Vector.o EarthModel.o IceModel.o Trigger.o Ray.o Tools.o Efficiencies.o Event.o Detector.o Position.o Spectra.o RayTrace.o RayTrace_IceModels.o signal.o secondaries.o Settings.o Primaries.o counting.o RaySolver.o Report.o eventSimDict.o veff_aeff2.o
CCFILE = Vector.cc EarthModel.cc IceModel.cc Trigger.cc Ray.cc Tools.cc Efficiencies.cc Event.cc Detector.cc Spectra.cc Position.cc RayTrace.cc signal.cc secondaries.cc RayTrace_IceModels.cc Settings.cc Primaries.cc counting.cc RaySolver.cc Report.cc veff_aeff2.cc
CLASS_HEADERS = Trigger.h Detector.h Settings.h Spectra.h IceModel.h Primaries.h Report.h Event.h secondaries.hh #need to add headers which added to Tree Branch
PROGRAMS = veff_aeff2
all : $(PROGRAMS)
veff_aeff2 : $(OBJS)
$(LD) $(OBJS) $(LDFLAGS) $(LIBS) -o $(PROGRAMS)
@echo "done."
#The library
$(ROOT_LIBRARY) : $(LIB_OBJS)
@echo "Linking $@ ..."
ifeq ($(PLATFORM),macosx)
# We need to make both the .dylib and the .so
$(LD) $(SOFLAGS)$@ $(LDFLAGS) $(G77LDFLAGS) $^ $(OutPutOpt) $@
ifneq ($(subst $(MACOSX_MINOR),,1234),1234)
ifeq ($(MACOSX_MINOR),4)
ln -sf $@ $(subst .$(DllSuf),.so,$@)
else
$(LD) -dynamiclib -undefined $(UNDEFOPT) $(LDFLAGS) $(G77LDFLAGS) $^ \
$(OutPutOpt) $(subst .$(DllSuf),.so,$@)
endif
endif
else
$(LD) $(SOFLAGS) $(LDFLAGS) $(G77LDFLAGS) $(LIBS) $(LIB_OBJS) -o $@
endif
##-bundle
#%.$(OBJSUF) : %.$(SRCSUF)
# @echo "<**Compiling**> "$<
# $(CXX) $(CXXFLAGS) -c $< -o $@
%.$(OBJSUF) : %.C
@echo "<**Compiling**> "$<
$(CXX) $(CXXFLAGS) $ -c $< -o $@
%.$(OBJSUF) : %.cc
@echo "<**Compiling**> "$<
$(CXX) $(CXXFLAGS) $ -c $< -o $@
# added for fortran code compiling
%.$(OBJSUF) : %.f
@echo "<**Compiling**> "$<
$(G77) -c $<
eventSimDict.C: $(CLASS_HEADERS)
@echo "Generating dictionary ..."
@ rm -f *Dict*
rootcint $@ -c ${INC_ARA_UTIL} $(CLASS_HEADERS) ${ARA_ROOT_HEADERS} LinkDef.h
clean:
@rm -f *Dict*
@rm -f *.${OBJSUF}
@rm -f $(LIBRARY)
@rm -f $(ROOT_LIBRARY)
@rm -f $(subst .$(DLLSUF),.so,$(ROOT_LIBRARY))
@rm -f $(TEST)
#############################################################################
|
| Attachment 4: make_plot.py
|
# -*- coding: utf-8 -*-
import numpy as np
import sys
import matplotlib.pyplot as plt
from pylab import setp
from matplotlib.pyplot import rcParams
import csv
import pandas as pd
rcParams['mathtext.default'] = 'regular'
def read_file(finame):
fi = open(finame, 'r')
rdr = csv.reader(fi, delimiter=',', skipinitialspace=True)
table = []
for row in rdr:
# print(row)
energy = float(row[0])
veff = float(row[1])
veff_err = float(row[2])
veff1 = float(row[3])
veff2 = float(row[4])
veff3 = float(row[5])
row = {'energy':energy, 'veff':veff, 'veff_err':veff_err, 'veff1':veff1, 'veff2':veff2, 'veff3':veff3}
table.append(row)
df=pd.DataFrame(table)
df_ordered=df.sort_values('energy',ascending=True)
# print(df_ordered)
return df_ordered
def beautify_veff(this_ax):
sizer=20
xlow = 1.e16 #the lower x limit
xup = 2.e20 #the uppper x limit
ylow =1e-3 #the lower x limit
yup = 6.e1 #the uppper x limit
this_ax.set_xlabel('Energy [eV]',size=sizer) #give it a title
this_ax.set_ylabel('[V$\Omega]_{eff}$ [km$^3$sr]',size=sizer)
this_ax.set_yscale('log')
this_ax.set_xscale('log')
this_ax.tick_params(labelsize=sizer)
this_ax.set_xlim([xlow,xup]) #set the x limits of the plot
this_ax.set_ylim([ylow,yup]) #set the y limits of the plot
this_ax.grid()
this_legend = this_ax.legend(loc='upper left')
setp(this_legend.get_texts(), fontsize=17)
setp(this_legend.get_title(), fontsize=17)
def main():
"""
arasim_energies = np.array([3.16e+16, 1e+17, 3.16e+17, 1e+18, 3.16e+18, 1e+19, 3.16e+19, 1e+20])
arasim_energies2 = np.array([3.16e+16, 1e+17, 3.16e+17, 1e+18, 1e+19, 3.16e+19, 1e+20])
#arasim_desA_veff = np.array([0.080894,0.290695,0.943223,2.388708,4.070498,6.824112,10.506490,13.969418])
arasim_desA_veff_s = np.array([0.067384,0.289591,0.996509,2.464464,4.945600,8.735506,13.357300,18.751915])
# arasim_desA_veff_ice = np.array([0.066291,0.303620,0.927647,2.427554,4.962093,8.465895,13.425852,18.706528])
arasim_desA_error = np.array([0.008367,0.017401,0.032222,0.084223,0.119240,0.221266,0.272460,0.320456])
# arasim_desA_error_ice = np.array([0.008345,0.017748,0.031295,0.083747,0.119370,0.217717,0.273075,0.319872])
# arasim_desB_veff = np.array([0.124111,0.417102,1.310555,3.648494,4.070E+0,11.675189,18.393961,13.688909])
arasim_desB_veff_s = np.array([0.064937,0.355747,1.289947,3.821705,8.002805,15.352981,25.391282,18.009545])
# arasim_desB_veff_ice = np.array([0.080070,0.340927,1.369081,3.938550,8.211407,15.190858,25.066541,18.021033])
arasim_desB_error = np.array([0.008268,0.019317,0.036926,0.105445,0.152488,0.296617,0.377997,0.314314])
# arasim_desB_error_ice = np.array([0.009220,0.019144,0.037966,0.107189,0.154430,0.293728,0.375827,0.314244])
arasim_desA_veff1_sm = np.array([0.054,0.231,0.868,2.239,4.586,8.339,12.910,18.301])
arasim_desA_veff2_sm = np.array([0.009,0.036,0.098,0.323,0.705,1.167,2.169,3.130])
arasim_desA_veff3_sm = np.array([0.011,0.074,0.193,0.664,1.208,2.347,3.689,4.814])
arasim_desA_veff1 = np.array([0.053,0.282,1.108,3.462,7.486,14.613,24.682,17.557])
arasim_desA_veff2 = np.array([0.006,0.040,0.121,0.322,0.636,1.308,1.961,2.708])
arasim_desA_veff3 = np.array([0.006,0.065,0.188,0.638,1.187,2.270,3.301,4.322])
"""
veff_A=read_file("Veff_desA.txt")
veff_B=read_file("Veff_desB.txt")
print("desA is \n", veff_A)
print("desB is \n", veff_B)
dfA = veff_A[['veff','veff_err']]
dfB = veff_B[['veff','veff_err']]
dfA.to_csv('veffA.csv', sep='\t',index=False)
dfB.to_csv('veffB.csv', sep='\t',index=False)
fig = plt.figure(figsize=(11,8.5))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(veff_A['energy'], veff_A['veff'],'bs-',label='Strawman',markersize=8,linewidth=2)
ax1.plot(veff_B['energy'], veff_B['veff'],'gs-',label='Punch @100 m',markersize=8,linewidth=2)
ax1.fill_between(veff_A['energy'], veff_A['veff']-veff_A['veff_err'], veff_A['veff']+veff_A['veff_err'], alpha=0.2, color='red')
ax1.fill_between(veff_B['energy'], veff_B['veff']-veff_B['veff_err'], veff_B['veff']+veff_B['veff_err'], alpha=0.2, color='red')
beautify_veff(ax1)
ax1.set_title("Punch vs Strawman, noise + signal",fontsize=20)
fig.savefig("desAB.png",edgecolor='none',bbox_inches="tight") #save the figure
fig2 = plt.figure(figsize=(11,8.5))
ax2 = fig2.add_subplot(1,1,1)
#Triggers plot
ax2.plot(veff_B['energy'], veff_B['veff1'],'g^-',label='Phased array (Punch)',markersize=8,linewidth=2)
ax2.plot(veff_A['energy'], veff_A['veff1'],'gs-',label='Phased array (Strawman)',markersize=8,linewidth=2)
ax2.plot(veff_B['energy'], veff_B['veff2'],'b^-',label='LPDAs (Punch)',markersize=8,linewidth=2)
ax2.plot(veff_A['energy'], veff_A['veff2'],'bs-',label='LPDAs (Strawman)',markersize=8,linewidth=2)
ax2.plot(veff_B['energy'], veff_B['veff3'],'y^-',label='Bicones (Punch)',markersize=8,linewidth=2)
ax2.plot(veff_A['energy'], veff_A['veff3'],'ys-',label='Bicones (Strawman)',markersize=8,linewidth=2)
ax2.set_title("Triggers contribution, noise + signal",fontsize=20)
beautify_veff(ax2)
fig2.savefig("desAB_triggers.png",edgecolor='none',bbox_inches="tight") #save the figure
main()
|
|
34
|
Tue Feb 26 16:19:20 2019 |
Julie Rolla | All of the group GitHub account links | Software | ANITA Binned Analysis: https://github.com/osu-particle-astrophysics/BinnedAnalysis
GENETIS Bicone: https://github.com/mclowdus/BiconeEvolution
GENETIS Dipole: https://github.com/hchasan/XF-Scripts
ANITA Build tool: https://github.com/anitaNeutrino/anitaBuildTool
ANITA Hackathon: https://github.com/anitaNeutrino/hackathon2017
ICEMC: https://github.com/anitaNeutrino/icemc
Brian's Github: https://github.com/clark2668?tab=repositories
Note that you *may* need permissions for some of these. Please email Lauren (ennesser.1@buckeyemail.osu.edu ), Julie (JulieRolla@gmail.com), AND Brian (clark.2668@buckeyemail.osu.edu ) if you have any issues with permssions. Please state which GitHub links you are looking to view. |
|
3
|
Wed Mar 22 18:01:23 2017 |
Brian Clark | Advice for Using the Ray Trace Correlator | Analysis | If you are trying to use the Ray Trace Correlator with AraRoot, you will probably encounter some issues as you go. Here is some advice that Carl Pfendner found, and Brian Clark compiled.
Please note that it is extremely important that your AntennaInfo.sqlite table in araROOT contain the ICRR versions of both Testbed and Station1. Testbed seems to have fallen out of practice of being included in the SQL table. Also, Station1 is the ICRR (earliest) version of A1, unlike the ATRI version which is logged as ARA01. This will cause seg faults in the intial setup of the timing and geometry arrays that seem unrelated to pure geometry files. If you get a seg-fault in the "setupSizes" function or the Detector call of the "setupPairs" function, checking your SQL file is a good idea. araROOT branch 3.13 has such a source table with Testbed and Station1 included.
Which combination of Makefile/Makefile.arch/StandardDefinitions.mk works can be machine specific (frustratingly). Sometimes the best StandardDefinitions.mk is found int he make_timing_arrays example.
Common Things to Check
1: Did you "make install" the Ray Trace Correlator after you made it?
2: Do you have the setup.txt file?
3: Do you have the "data" directory?
Common Errors
1: If the Ray Trace Correlator compiles, and you execute a binary, and get the following:
******** Begin Correlator ********, this one!
Pre-icemodel test
terminate called after throwing an instance of 'std::out_of_range'
what(): basic_string::substr
Aborted
Check to make sure have the "data" directory.
|
|
40
|
Thu Jul 25 16:50:43 2019 |
Dustin Nguyen | Advice (not mine) on writing HEP stuff | Other | PDF of advice by Andy Buckley (U Glasgow) on writing a HEP thesis (and presumably HEP papers too) that was forwarded by John Beacom to the CCAPP mailing list a few months back. |
| Attachment 1: thesis-writing-gotchas.pdf
|
|
22
|
Sun Apr 29 21:44:15 2018 |
Brian Clark | Access Deep ARA Station Data | Analysis | Quick C++ program for pulling waveforms out of deep ARA station data. If you are using AraRoot, you would put this inside your "analysis" directory and add it to your CMakeLists.txt. |
| Attachment 1: plot_deep_station_event.cxx
|
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//// plot_event.cxx
//// plot deep station event
////
//// Apr 2018, clark.2668@osu.edu
////////////////////////////////////////////////////////////////////////////////
//Includes
#include <iostream>
//AraRoot Includes
#include "RawAtriStationEvent.h"
#include "UsefulAtriStationEvent.h"
#include "AraEventCalibrator.h"
//ROOT Includes
#include "TTree.h"
#include "TFile.h"
#include "TGraph.h"
#include "TCanvas.h"
using namespace std;
int main(int argc, char **argv)
{
//check to make sure they've given me a run and a pedestal
if(argc<2) {
std::cout << "Usage\n" << argv[0] << " <run file>\n";
std::cout << "e.g.\n" << argv[0] << " event1841.root\n";
return 0;
}
char pedFileName[200];
sprintf(pedFileName, "%s", argv[1]);
printf("------------------------------------------------------------------------\n");
printf("%s\n", argv[0]);
printf("runFileName %s\n", argv[1]);
printf("------------------------------------------------------------------------\n");
TFile *fp = TFile::Open(argv[2]);
if(!fp) {
std::cerr << "Can't open file\n";
return -1;
}
TTree *eventTree = (TTree*) fp->Get("eventTree");
if(!eventTree) {
std::cerr << "Can't find eventTree\n";
return -1;
}
RawAtriStationEvent *rawAtriEvPtr = 0; //empty pointer
eventTree->SetBranchAddress("event",&rawAtriEvPtr); //set the branch address
Int_t run_num; //run number of event
eventTree->SetBranchAddress("run", &run_num); //set the branch address
int numEntries=eventTree->GetEntries(); //get the number of events
int stationId=0;
eventTree->GetEntry(0);
stationId = rawAtriEvPtr->stationId; //assign the statio id number
AraEventCalibrator *calib = AraEventCalibrator::Instance(); //make a calibrator
for(int event=0;event<numEntries;event++) {
//for(int event=0;event<700;event++) {
eventTree->GetEntry(event); //get the event
int evt_num = rawAtriEvPtr->eventNumber; //check the event umber
if(rawAtriEvPtr->isCalpulserEvent()==0) continue; //bounce out if it's not a cal pulser
UsefulAtriStationEvent *realAtriEvPtr_fullcalib = new UsefulAtriStationEvent(rawAtriEvPtr, AraCalType::kLatestCalib); //make the event
TGraph *waveforms[16]={0};
for(int i=0; i<16; i++){
waveforms[i]=realAtriEvPtr_fullcalib->getGraphFromRFChan(i);
}
TCanvas *canvas = new TCanvas("","",1000,1000);
canvas->Divide(4,4);
for(int i=0; i<16; i++){
canvas->cd(i+1);
waveforms[i]->Draw("alp");
}
char title[200];
sprintf(title,"waveforms_run%d_event%d.pdf",stationId,run_num,evt_num);
canvas->SaveAs(title);
delete canvas;
}
}
|
|
45
|
Fri Feb 4 13:06:25 2022 |
William Luszczak | "Help! AnitaBuildTools/PueoBuilder can't seem to find FFTW!" | Software | Disclaimer: This might not be the best solution to this problem. I arrived here after a lot of googling and stumbling across this thread with a similar problem for an unrelated project: https://github.com/xtensor-stack/xtensor-fftw/issues/52. If you're someone who actually knows cmake, maybe you have a better solution.
When compiling both pueoBuilder and anitaBuildTools, I have run into a cmake error that looks like:
CMake Error at /apps/cmake/3.17.2/share/cmake-3.17/Modules/FindPackageHandleStandardArgs.cmake:164 (message):
Could NOT find FFTW (missing: FFTW_LIBRARIES)
(potentially also missing FFTW_INCLUDES). Directing CMake to the pre-existing FFTW installations on OSC does not seem to do anything to resolve this error. From what I can tell, this might be related to how FFTW is built, so to get around this we need to build our own installation of FFTW using cmake instead of the recommended build process. To do this, grab the whatever version of FFTW you need from here: http://www.fftw.org/download.html (for example, I needed 3.3.9). Untar the source file into whatever directory you're working in:
tar -xzvf fftw-3.3.9.tar.gz
Then make a build directory and cd into it:
mkdir install
cd install
Now build using cmake, using the flags shown below.
cmake -DCMAKE_INSTALL_PREFIX=$(path_to_install_loc) -DBUILD_SHARED_LIBS=ON -DENABLE_OPENMP=ON -DENABLE_THREADS=ON ../fftw-3.3.9
For example, I downloaded and untarred the source file in `/scratch/wluszczak/fftw/`, and my install prefix was `/scratch/wluszczak/fftw/install/`. In principle this installation prefix can be anywhere you have write access, but for the sake of organization I usually try to keep everything in one place.
Once you have configured cmake, go ahead and install:
make install -j $(nproc)
Where $(nproc) is the number of threads you want to use. On OSC I used $(nproc)=4 for compiling the ANITA tools and it finished in a reasonable amount of time.
Once this has finished, cd to your install directory and remove everything except the `include` and `lib64` folders:
cd $(path_to_install_dir) #You might already be here if you never left
rm *
rm -r CMakeFiles
Now we need to rebuild with slightly different flags:
cmake -DCMAKE_INSTALL_PREFIX=$(path_to_install_loc) -DBUILD_SHARED_LIBS=ON -DENABLE_OPENMP=ON -DENABLE_THREADS=ON -DENABLE_FLOAT=ON ../fftw-3.3.9
make install -j $(nproc)
At the end of the day, your fftw install directory should have the following files:
include/fftw3.f
include/fftw3.f03
include/fftw3.h
include/fftw3l.f03
include/fftw3q.f03
lib64/libfftw3f.so
lib64/libfftw3f_threads.so.3
lib64/libfftw3_omp.so.3.6.9
lib64/libfftw3_threads.so
lib64/libfftw3f_omp.so
lib64/libfftw3f.so.3
lib64/libfftw3f_threads.so.3.6.9
lib64/libfftw3.so
lib64/libfftw3_threads.so.3
lib64/libfftw3f_omp.so.3
lib64/libfftw3f.so.3.6.9
lib64/libfftw3_omp.so
lib64/libfftw3.so.3
lib64/libfftw3_threads.so.3.6.9
lib64/libfftw3f_omp.so.3.6.9
lib64/libfftw3f_threads.so
lib64/libfftw3_omp.so.3
lib64/libfftw3.so.3.6.9
Once fftw has been installed, export your install directory (the one with the include and lib64 folders) to the following environment variable:
export FFTWDIR=$(path_to_install_loc)
Now you should be able to cd to your anitaBuildTools directory (or pueoBuilder directory) and run their associated build scripts:
./buildAnita.sh
or:
./pueoBuilder.sh
And hopefully your tools will magically compile (or at least, you'll get a new set of errors that are no longer related to this problem).
If you're running into an error that looks like:
CMake Error: The following variables are used in this project, but they are set to NOTFOUND.
Please set them or make sure they are set and tested correctly in the CMake files:
FFTWF_LIB (ADVANCED)
then pueoBuilder/anitaBuildTools can't seem to find your fftw installation (or files that are supposed to be included in that installation), try rebuilding FFTW with different flags according to which files it seems to think are missing.
If it seems like pueoBuilder can't seem to find your FFTW installation at all (i.e. you're getting some error that looks like missing: FFTW_LIBRARIES or missing: FFTW_INCLUDES, check the environment variables that are supposed to point to your local FFTW installation (`$FFTWDIR`) and make sure there are the correct files in the `lib` and `include` subdirectories. |
|
4
|
Fri Mar 31 11:36:54 2017 |
Brian Clark, Hannah Hasan, Jude Rajasekera, and Carl Pfendner | Installing Software Pre-Requisites for Simulation and Analysis | Software | Instructions on installing simulation software prerequisites (ROOT, Boost, etc) on Linux computers. |
| Attachment 1: installation.pdf
|
| Attachment 2: Installation-Instructions.tar.gz
|
|
Draft
|
Fri Jul 28 17:57:06 2017 |
Brian Clark and Ian Best | | | |
|