|
Updates and Results
Talks and Posters
Advice
Ideas
Important Figures
Write-Ups
Outreach
How-To
Funding Opportunities
GENETIS
|
| Place to document instructions for how to do things |
 |
|
|
Message ID: 27
Entry time: Mon Oct 1 19:06:59 2018
|
| Author: |
Brian Clark |
| Subject: |
Code to Compute Effective Volumes in AraSim |
| Project: |
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...." |
|
#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;
}
|
|
#############################################################################
##
##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)
#############################################################################
|