#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 "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 "TVector3.h"
#include "TRotation.h"
#include "TSpline.h"
#include "Math/InterpolationTypes.h"
#include "Math/Interpolator.h"
#include "Math/Integrator.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TLatex.h"
#include "Constants.h"
#include "Settings.h"
#include "Position.h"
#include "EarthModel.h"
#include "Tools.h"
#include "Vector.h"
#include "IceModel.h"
#include "Trigger.h"
#include "Spectra.h"
#include "signal.hh"
#include "secondaries.hh"
#include "Ray.h"
#include "counting.hh"
#include "Primaries.h"
#include "Efficiencies.h"
#include "Event.h"
#include "Detector.h"
#include "Report.h"
using namespace std;
/////////////////Plotting Script for nuflavorint for AraSIMQC
/////////////////Created by Kaeli Hughes, modified by Brian Clark
/////////////////Prepared on April 28 2016 as an "introductor script" to plotting with AraSim; for new users of queenbee
class EarthModel;
int main(int argc, char **argv)
{
gStyle->SetOptStat(111111); //this is a selection of statistics settings; you should do some googling and figure out exactly what this particular comibation does
gStyle->SetOptDate(1); //this tells root to put a date and timestamp on whatever plot we output
if(argc<2){ //if you don't have at least 1 file to run over, then you haven't given it a file to analyze; this checks this
cout<<"Not enough arguments! Abort run."<<endl;
}
//Create the histogram
TCanvas *c2 = new TCanvas("c2", "nuflavorint", 1100,850); //make a canvas on which to plot the data
TH1F *nuflavorint_hist = new TH1F("nuflavorint_hist", "nuflavorint histogram", 3, 0.5, 3.5); //create a histogram with three bins
nuflavorint_hist->GetXaxis()->SetNdivisions(3); //set the number of divisions
int total_thrown=0; // a variable to hold the grant total number of events thrown for all input files
for(int i=1; i<argc;i++){ //loop over the input files
string readfile; //create a variable called readfile that will hold the title of the simulation file
readfile = string(argv[i]); //set the readfile variable equal to the filename
Event *event = 0; //create a Event class pointer called event; note that it is set equal to zero to avoid creating a bald pointer
Report *report=0; //create a Event class pointer called event; note that it is set equal to zero to avoid creating a bald pointer
Settings *settings = 0; //create a Event class pointer called event; note that it is set equal to zero to avoid creating a bald pointer
TFile *AraFile=new TFile(( readfile ).c_str()); //make a new file called "AraFile" that will be simulation file we are reading in
if(!(AraFile->IsOpen())) return 0; //checking to see if the file we're trying to read in opened correctly; if not, bounce out of the program
int num_pass;//number of passing events
TTree *AraTree=(TTree*)AraFile->Get("AraTree"); //get the AraTree
TTree *AraTree2=(TTree*)AraFile->Get("AraTree2"); //get the AraTree2
AraTree2->SetBranchAddress("event",&event); //get the event branch
AraTree2->SetBranchAddress("report",&report); //get the report branch
AraTree->SetBranchAddress("settings",&settings); //get the settings branch
num_pass=AraTree2->GetEntries(); //get the number of passed events in the data file
AraTree->GetEvent(0); //get the first event; sometimes the tree does not instantiate properly if you'd explicitly "active" the first event
total_thrown+=(settings->NNU); //add the number of events from this file (The NNU variable) to the grand total; NNU is the number of THROWN neutrinos
for (int k=0; k<num_pass; k++){ //going to fill the histograms for as many events as were in this input file
AraTree2->GetEvent(k); //get the even from the tree
int nuflavorint; //make the container variable
double weight; //the weight of the event
int trigger; //the global trigger value for the event
nuflavorint=event->nuflavorint; //draw the event out; one of the objects in the event class is the nuflavorint, and this is the syntax for accessing it
weight=event->Nu_Interaction[0].weight; //draw out the weight of the event
trigger=report->stations[0].Global_Pass; //find out if the event was a triggered event or not
/*
if(trigger!=0){ //check if the event triggered
nuflavorint_hist->Fill(nuflavorint,weight); //fill the event into the histogram; the first argument of the fill (which is mandatory) is the value you're putting into the histogram; the second value is optional, and is the weight of the event in the histogram
//in this particular version of the code then, we are only plotting the TRIGGERED events; if you wanted to plot all of the events, you could instead remove this "if" condition and just Fill everything
}
*/
nuflavorint_hist->Fill(nuflavorint,weight); //fill the event into the histogram; the first argument of the fill (which is mandatory)
}
}
//After looping over all the files, make the plots and save them
//do some stuff to get the "total events thrown" text box ready; this is useful because on the plot itself you can then see how many events you THREW into the simulation; this is especially useful if you're only plotting passed events, but want to know what fraction of your total thrown that is
char *buffer= (char*) malloc (250); //declare a buffer pointer, and allocate it some chunk of memory
int a = snprintf(buffer, 250,"Total Events Thrown: %d",total_thrown); //print the words "Total Events Thrown: %d" to the variable "buffer" and tell the system how long that phrase is; the %d sign tells C++ to replace that "%d" with the next argument, or in this case, the number "total_thrown"
if(a>=250){ //if the phrase is longer than the pre-allocated space, increase the size of the buffer until it's big enough
buffer=(char*) realloc(buffer, a+1);
snprintf(buffer, 250,"Total Events Thrown: %d",total_thrown);
}
TLatex *u = new TLatex(.3,.01,buffer); //create a latex tex object which we can draw on the cavas
u->SetNDC(kTRUE); //changes the coordinate system for the tex object plotting
u->SetIndiceSize(.1); //set the size of the latex index
u->SetTextSize(.025); //set the size of the latex text
nuflavorint_hist->Draw(); //draw the histogram
nuflavorint_hist->GetXaxis()->SetTitle("Neutrino Flavor"); //set the x-axis label
nuflavorint_hist->GetYaxis()->SetTitle("Number of Events (weighted)"); //set the y-axis label
nuflavorint_hist->GetYaxis()->SetTitleOffset(1.5); //set the separation between the y-axis and it's label; root natively makes this smaller than is ideal
nuflavorint_hist->SetLineColor(kBlack); //set the color of the histogram to black, instead of root's default navy blue
u->Draw(); //draw the statistics box information
c2->SaveAs("outputs/plotting_example.png"); //save the canvas as a JPG file for viewing
c2->SaveAs("outputs/plotting_example.pdf"); //save the canvas as a PDF file for viewing
c2->SaveAs("outputs/plotting_example.root"); //save the canvas as a ROOT file for viewing or editing later
} //end main; this is the end of the script
|
#############################################################################
## 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 plotting_example on 28 April 2016
##
##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 plotting_example.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 plotting_example.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 = plotting_example
all : $(PROGRAMS)
plotting_example : $(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)
#############################################################################
|
NFOUR=1024
EXPONENT=21
NNU=300 // number of neutrino events
NNU_PASSED=10 // number of neutrino events that are allowed to pass the trigger
ONLY_PASSED_EVENTS=0 // 0 (default): AraSim throws NNU events whether or not they pass; 1: AraSim throws events until the number of events that pass the trigger is equal to NNU_PASSED (WARNING: may cause long run times if reasonable values are not chosen)
NOISE_WAVEFORM_GENERATE_MODE=0 // generate new noise waveforms for each events
NOISE_EVENTS=16 // number of pure noise waveforms
TRIG_ANALYSIS_MODE=0 // 0 = signal + noise, 1 = signal only, 2 = noise only
DETECTOR=1 // ARA stations 1 to 7
NOFZ=1
core_x=10000
core_y=10000
TIMESTEP=5.E-10 // value for 2GHz actual station value
TRIG_WINDOW=1.E-7 // 100ns which is actual testbed trig window
POWERTHRESHOLD=-6.06 // 100Hz global trig rate for 3 out of 16 ARA stations
POSNU_RADIUS=3000
V_MIMIC_MODE=0 // 0 : global trig is located center of readout windows
DATA_SAVE_MODE=0 // 2 : don't save any waveform informations at all
DATA_LIKE_OUTPUT=0 // 0 : don't save any waveform information to eventTree
BORE_HOLE_ANTENNA_LAYOUT=0
SECONDARIES=0
TRIG_ONLY_BH_ON=0
CALPULSER_ON=0
USE_MANUAL_GAINOFFSET=0
USE_TESTBED_RFCM_ON=0
NOISE_TEMP_MODE=0
TRIG_THRES_MODE=0
READGEOM=0 // reads geometry information from the sqlite file or not (0 : don't read)
TRIG_MODE=0 // use vpol, hpol separated trigger mode. by default N_TRIG_V=3, N_TRIG_H=3. You can change this values
number_of_stations=1
core_x=10000
core_y=10000
DETECTOR=1
DETECTOR_STATION=2
DATA_LIKE_OUTPUT=0
NOISE_WAVEFORM_GENERATE_MODE=0 // generates new waveforms for every event
NOISE=0 //flat thermal noise
NOISE_CHANNEL_MODE=0 //using different noise temperature for each channel
NOISE_EVENTS=16 // number of noise events which will be store in the trigger class for later use
ANTENNA_MODE=1
APPLY_NOISE_FIGURE=0
|