#!/bin/bash

#PBS -A PCON0003
#PBS -l walltime=01:00:00
#PBS -l nodes=1:ppn=20


touch commandfile
for value in {1..1000}
do
    line="/path/to/your_command_to_run $value (arg1) (arg2)..(argn)"
    echo ${line}>> commandfile
    
done


module load pcp
mpiexec parallel-command-processor commandfile

#!/bin/bash

#PBS -A PCON0003
#PBS -N effvol
#PBS -l walltime=01:00:00
#PBS -l nodes=1:ppn=20
#PBS -o ./log/out
#PBS -e ./log/err

source /users/PCON0003/osu10643/.bash_batch

cd $TMPDIR
touch effvolconf
for value in {1..1000}
do
    line="/users/PCON0003/osu10643/app/geant/app/nrt -m /users/PCON0003/osu10643/app/geant/app/nrt/effvol.mac -f \"$TMPDIR/effvol$value.root\""
    echo ${line}>> effvolconf
    
done


module load pcp
mpiexec parallel-command-processor effvolconf


cp $TMPDIR/*summary.root /users/PCON0003/osu10643/doc/root/summary/

Installing AnitaTools on OSC
Sam Stafford
04/27/2017

This document summarizes the issues I encountered installing AnitaTools on the OSC Oakley cluster.
I have indicated work-arounds I made for unexpected issues
  I do not know that this is the only valid process
  This process was developed by trial-and-error (mostly error) and may contain superfluous steps
    A person familiar with AnitaTools and cmake may be able to streamline it

Check out OSC's web site, particularly to find out about MODULES, which facilitate access to pre-installed software

export the following environment variables in your .bash_profile  (not .bashrc):

  ROOTSYS                             where you want ROOT to live
     install it in somewhere your user directory; at this time, ROOT is not pre-installed on Oakley as far as I can tell
  ANITA_UTIL_INSTALL_DIR              where you want anitaTools to live
  FFTWDIR                             where fftw is 
    look on OSC's website to find out where it is; you shouldn't have to install it locally

  PATH   should contain $FFTWDIR/bin  and $ROOTSYS/bin
  LD_LIBRARY_PATH should contain     $FFTWDIR/lib    $ROOTSYS/lib     $ANITA_UTIL_INSTALL_DIR/lib
  LD_INCLUDE_PATH should contain     $FFTWDIR/include    $ROOTSYS/include     $ANITA_UTIL_INSTALL_DIR/include

also put in your .bash_profile: (I put these after the exports)
  
    module load gnu/4.8.5     // loads g++ compiler
           (this should automatically load module fftw/3.3.4 also)

install ROOT  - follow ROOT's instructions to build from source.  It's a typical (configure / make / make install) sequence
  you probably need  ./configure --enable-Minuit2

get AnitaTools from github/anitaNeutrino "anitaBuildTool" (see Anita ELOG 672, by Cosmin Deaconu)

Change entry in which_event_reader to ANITA3, if you want to analyze ANITA-3 data
  (at least for now; I think they are developing smarts to make the SW adapt automatically to the anita data "version")

Do ./buildAnita.sh       //downloads the software and attempts a full build/install

    it may file on can't-find-fftw during configure:
      system fails to populate environment variable FFTW_ROOT, not sure why
      add the following line at beginning of anitaBuildTool/cmake/modules/FindFFTW.cmake:
        set( FFTW_ROOT /usr/local/fftw3/3.3.4-gnu)
          (this apparently tricks cmake into finding fftw)
  NOTE: ./buildAnita.sh always downloads the software from github.  IT WILL WIPE OUT ANY CHANGES YOU MADE TO AnitaTools!
      
Do "make" 
    May fail with   /usr/lib64/libstdc++.so.6: version `GLIBCXX_3.4.14' not found and/or a few other similar messages 
                   or may say c11 is not supported  
      need to change compiler setting for cmake:
        make the following change in anitaBuildTool/build/CMakeCache.txt   (points cmake to the g++ compiler instead of default intel/c++)
          #CMAKE_CXX_COMPILER:FILEPATH=/usr/bin/c++                  (comment this out)
           CMAKE_CXX_COMPILER:FILEPATH=/usr/local/gcc/4.8.5/bin/g++  (add this)  
       (you don't necessarfily have to use version, gcc/4.8.5, but it worked for me)

Then retry by doing "make"

Once make is completed, do "make install"

A couple of notes:
  Once AnitaTools is built, if you change source, just do make, then make install (from anitaBuildTool) (don't ./buildAnita.sh; see above)
      (actually make install will do the make step if source changes are detected)
  To start AnitaTools over from cmake, delete the anitaBuildTool/build directory and run make (not cmake: make will drive cmake for you)
      (don't do cmake directly unless you know what you're doing; it'll mess things up)
  

A few notes about installing icemc on OSC

Dependencies
  ROOT - download from CERN and install according to instructions
  FFTW - do "module load gnu/4.8.5"   (or put it in your .bash_profile)

The environment variable FFTWDIR must contain the directory where FFTW resides
  in my case this was /usr/local/fftw3/3.3.4-gnu 
  set this up in your .bash_profile (not .bashrc)

I copied my working instance of icemc from my laptop to a folder in my osc space
  Copy the whole icemc directory (maybe its icemc/trunk, depending on how you installed), EXCEPT for the "output" subdir because it's big and unnecessary
  in your icemc directory on OSC, do "mkdir output"

In icemc/Makefile
  find a statement like this:
    LIBS += -lMathMore $(FFTLIBS) -lAnitaEvent
  and modify it to include the directory where the FFTW library is:
    LIBS += -L$(FFTWDIR)/lib -lMathMore $(FFTLIBS) -lAnitaEvent

  note: FFTLIBS contains the list of libraries (e.g., -lfftw3), NOT the library search paths
  
Compile by doing "make"

Remember you should set up a batch job on OSC using PBS.

## annotated sample PBS batch job specification for OSC
## Sam Stafford 05/11/2017

#PBS -N j_ai06_${RUN_NUMBER}

##PBS -m abe       ##  request an email on job completion
#PBS -l mem=16GB   ##  request 16GB memory
##PBS -l walltime=06:00:00  ## set this in qsub
#PBS -j oe         ## merge stdout and stderr into a single output log file
#PBS -A PAS0174

echo "run number " $RUN_NUMBER
echo "cal pulser " $CAL_PULSER
echo "baseline file " $BASELINE_FILE
echo "temp dir is " $TMPDIR
echo "ANITA_DATA_REMOTE_DIR="$ANITA_DATA_REMOTE_DIR
set -x


## copy the files from kingbee to the temporary workspace  
##    (if you set up public key authentication between kingbee and OSC, you won't need a password; just google "public key authentication")
mkdir $TMPDIR/run${RUN_NUMBER}   ## make a directory for this run number
scp stafford.16@kingbee.mps.ohio-state.edu:/data/anita/anita3/flightData/copiedBySam/run${RUN_NUMBER}/calEventFile${RUN_NUMBER}.root $TMPDIR/run${RUN_NUMBER}/calEventFile${RUN_NUMBER}.root
scp stafford.16@kingbee.mps.ohio-state.edu:/data/anita/anita3/flightData/copiedBySam/newerData/run${RUN_NUMBER}/gpsEvent${RUN_NUMBER}.root $TMPDIR/run${RUN_NUMBER}/gpsEvent${RUN_NUMBER}.root
scp stafford.16@kingbee.mps.ohio-state.edu:/data/anita/anita3/flightData/copiedBySam/newerData/run${RUN_NUMBER}/timedHeadFile${RUN_NUMBER}.root $TMPDIR/run${RUN_NUMBER}/timedHeadFile${RUN_NUMBER}.root
scp stafford.16@kingbee.mps.ohio-state.edu:/data/anita/anita3/flightData/copiedBySam/newerData/run${RUN_NUMBER}/decBlindHeadFileV1_${RUN_NUMBER}.root $TMPDIR/run${RUN_NUMBER}/decBlindHeadFileV1_${RUN_NUMBER}.root

## set up the environment variables to point to the temporary work space
export ANITA_DATA_REMOTE_DIR=$TMPDIR
export ANITA_DATA_LOCAL_DIR=$TMPDIR
echo "ANITA_DATA_REMOTE_DIR="$ANITA_DATA_REMOTE_DIR

## run the analysis program
cd analysisSoftware
./analyzerIterator06 ${CAL_PULSER} -S1 -Noverlap --FILTER_OPTION=4 ${BASELINE_FILE} ${RUN_NUMBER} -O
echo "batch job ending"

#!/bin/bash
#Jude Rajasekera 3/20/2017
shelfmcDir=/users/PCON0003/cond0091/ShelfMC #put your shelfmc directory address here 
 

runName='TestRun' #name of run
NNU=500000 #total NNU per run
seed=42 #initial seed for every run,each processor will recieve a different seed (42,43,44,45...)
NNU="$(($NNU / 20))" #calculating NNU per processor, change 20 to however many processor your cluster has per node
ppn=20 #processors per node
########################### make changes for input.txt file here #####################################################
input1="#inputs for ARIANNA simulation, do not change order unless you change ReadInput()"
input2="$NNU #NNU, setting to 1 for unique neutrino"
input3="$seed   #seed Seed for Rand3"
input4="18.0    #EXPONENT, !should be exclusive with SPECTRUM"
input5="1000    #ATGap, m, distance between stations"
input6="4       #ST_TYPE, !restrict to 4 now!"
input7="4       #N_Ant_perST, not to be confused with ST_TYPE above"
input8="2       #N_Ant_Trigger, this is the minimum number of AT to trigger"
input9="30      #Z for ST_TYPE=2"
input10="575     #ICETHICK, thickness of ice including firn, 575m at Moore's Bay"
input11="1      #FIRN, KD: ensure DEPTH_DEPENDENT is off if FIRN is 0"
input12="1.30   #NFIRN 1.30"
input13="$122    #FIRNDEPTH in meters"
input14="1      #NROWS 12 initially, set to 3 for HEXAGONAL"
input15="1      #NCOLS 12 initially, set to 5 for HEXAGONAL"
input16="0      #SCATTER"
input17="1      #SCATTER_WIDTH,how many times wider after scattering"
input18="0      #SPECTRUM, use spectrum, ! was 1 initially!"
input19="0      #DIPOLE,  add a dipole to the station, useful for st_type=0 and 2"
input20="0      #CONST_ATTENLENGTH, use constant attenuation length if ==1"
input21="1000   #ATTEN_UP, this is the conjuction of the plot attenlength_up and attlength_down when setting REFLECT_RATE=0.5(3dB)"
input22="250    #ATTEN_DOWN, this is the average attenlength_down before Minna Bluff measurement(not used anymore except for CONST_ATTENLENGTH)"
input23="4      #NSIGMA, threshold of trigger"
input24="1      #ATTEN_FACTOR, change of the attenuation length"
input25="1      #REFLECT_RATE,power reflection rate at the ice bottom"
input26="0      #GZK, 1 means using GZK flux, 0 means E-2 flux"
input27="0      #FANFLUX, use fenfang's flux which only covers from 10^17 eV to 10^20 eV"
input28="0      #WIDESPECTRUM, use 10^16 eV to 10^21.5 eV as the energy spectrum, otherwise use 17-20"
input29="1      #SHADOWING"
input30="1      #DEPTH_DEPENDENT_N;0 means uniform firn, 1 means n_firn is a function of depth"
input31="0      #HEXAGONAL"
input32="1      #SIGNAL_FLUCT 1=add noise fluctuation to signal or 0=do not"
input33="4.0    #GAINV  gain dependency"
input34="1      #TAUREGENERATION if 1=tau regeneration effect, if 0=original"
input35="3.0    #ST4_R radius in meters between center of station and antenna"
input36="350    #TNOISE noise temperature in Kelvin"
input37="80     #FREQ_LOW low frequency of LPDA Response MHz #was 100"
input38="1000   #FREQ_HIGH high frequency of LPDA Response MHz"
input39="/users/PCON0003/cond0091/ShelfMC/temp/LP_gain_manual.txt     #GAINFILENAME"
###########################################################################################

cd $shelfmcDir #cd to dir containing shelfmc
mkdir $runName #make a folder for run
cd $runName    #cd into run folder
initSeed=$seed 

for (( i=1; i<=$ppn;i++)) #make 20 setup files for 20 processors
do
    mkdir Setup$i #make setup folder i
    cd Setup$i #go into setup folder i
    
    seed="$(($initSeed+$i-1))" #calculate seed for this iteration
    input3="$seed      #seed Seed for Rand3" #save new input line    
    for j in {1..40} #print all input.txt lines 
    do
	lineName=input$j
        echo "${!lineName}" >> input.txt #print line to input.txt file
    done
    cd ..
done



pwd=`pwd`

#create job file
echo '#!/bin/bash' >> run_shelfmc_multithread.sh
echo '#PBS -l nodes=1:ppn='$ppn >> run_shelfmc_multithread.sh #change depending on processors per node
echo '#PBS -l walltime=00:05:00' >> run_shelfmc_multithread.sh #change walltime depending on run size, will be 20x shorter than single processor run time
echo '#PBS -N shelfmc_'$runName'_job' >> run_shelfmc_multithread.sh
echo '#PBS -j oe'  >> run_shelfmc_multithread.sh
echo '#PBS -A PCON0003' >> run_shelfmc_multithread.sh #change to specify group
echo 'cd ' $shelfmcDir >> run_shelfmc_multithread.sh 
echo 'runName='$runName  >> run_shelfmc_multithread.sh
for (( k=1; k<=$ppn;k++))
do
    echo './shelfmc_stripped.exe $runName/Setup'$k' _'$k'$runName &' >> run_shelfmc_multithread.sh #execute commands for 20 setup files
done
echo 'wait' >> run_shelfmc_multithread.sh #wait until all runs are finished
echo 'cd $runName' >> run_shelfmc_multithread.sh #go into run folder 
echo 'for (( i=1; i<='$ppn';i++)) #20 iterations' >> run_shelfmc_multithread.sh 
echo 'do' >> run_shelfmc_multithread.sh
echo '  cd Setup$i #cd into setup dir' >> run_shelfmc_multithread.sh
echo '  mv *.root ..' >> run_shelfmc_multithread.sh #move root files to runDir
echo '  cd ..' >> run_shelfmc_multithread.sh
echo 'done' >> run_shelfmc_multithread.sh
echo 'hadd Result_'$runName'.root *.root' >> run_shelfmc_multithread.sh #add all root files
echo 'rm *ShelfMCTrees*' >> run_shelfmc_multithread.sh #delete all partial root files

chmod u+x run_shelfmc_multithread.sh

echo "Run files created"
echo "cd into run folder and do $ qsub run_shelfmc_multithread.sh"

This document will explain how to dowload, configure, and run multithread_shelfmc.sh in order to do large runs on computing clusters.

####DOWNLOAD####
1.Download multithread_shelfmc.sh
2.Move multithread_shelfmc.sh into ShelfMC directory
3.Do $chmod u+x multithread_shelfmc.sh

####CONFIGURE###
1.Open multithread_shelfmc.sh
2.On line 3, modify shelfmcDir to your ShelfMC dir
3.On line 6, add your run name
4.On line 7, add the total NNU
5.On line 8, add an intial seed
6.On line 10, specify number of processors per node for your cluster
7.On lines 12-49, edit the input.txt parameters
8.On line 50, add the location of your LP_gain_manual.txt
9.On line 80, specify a wall time for each run, remember this will be about 20x shorter than ShelfMC on a single processor
10.On line 83, Specify the group name for your cluster if needed
11.Save file

####RUN####
1.Do $./multithread_shelfmc.sh 
2.There should now be a new directory in the ShelfMC dir with 20 setup files and a run_shelfmc_multithread.sh script
3.Do $qsub run_shelfmc_multithread.sh

###RESULT####
1.After the run has completed, there will be a result .root file in the run directory
 

This document will explain how to dowload, configure, and run a parameter space search for ShelfMC on a computing cluster. 
These scripts explore the ShelfMC parameter space by varying ATTEN_UP, REFLECT_RATE, ICETHICK, FIRNDEPTH, and STATION_DEPTH for certain rages. 
The ranges and increments can be found in setup.sh. 

In order to vary STATION_DEPTH, some changes were made to the ShelfMC code. Follow these steps to allow STATION_DEPTH to be an input parameter.
1.cd to ShelfMC directory
2.Do $sed -i -e 's/ATDepth/STATION_DEPTH/g' *.cc
3.Open declaration.hh. Replace line 87 "const double ATDepth = 0.;" with "double STATION_DEPTH;"
4.In functions.cc go to line 1829. This is the ReadInput() method. Add the lines below to the end of this method. 
   GetNextNumber(inputfile, number); // new line for station Depth
   STATION_DEPTH  = (double) atof(number.c_str()); //new line
5.Do $make clean all

#######Script Descriptions########
setup.sh -> This script sets up the necessary directories and setup files for all the runs
scheduler.sh -> This script submits and monitors all jobs. 


#######DOWNLOAD########
1.Download setup.sh and scheduler.sh
2.Move both files into your ShelfMC directory
3.Do $chmod u+x setup.sh and $chmod u+x scheduler.sh

######CONFIGURE#######
1.Open setup.sh
2.On line 4, modify the job name
3.On line 6, modify group name
4.On line 10, specify your ShelfMC directory
5.On line 13, modify your run name
6.On line 14, specify the NNU per run
7.On line 15, specify the starting seed
8.On line 17, specify the number of processors per node on your cluster
9.On lines 19-56, edit the input.txt parameters that you want to keep constant for every run
10.On line 57, specify the location of the LP_gain_manual.txt
11.On line 126, change walltime depending on total NNU. Remember this wall time will be 20x shorter than a single processor run.
12.On line 127, change job prefix
13.On line 129, change the group name if needed 
14.Save file
15.Open scheduler.sh
16.On line 4, specify your ShelfMC directory
17.On line 5, modify run name. Make sure it is the same runName as you have in setup.sh
18.On lines 35 and 39, replace cond0091 with your username for the cluster
19.On line 42, you can pick how many nodes you want to use at any given time. It is set to 6 intially. 
20.Save file 

#######RUN#######
1.Do $qsub setup.sh
2.Wait for setup.sh to finish. This script is creating the setup files for all runs. This may take about an hour.
3.When setup.sh is done, there should be a new directory in your home directory. Move this directory to your ShelfMC directory.
4.Do $screen to start a new screen that the scheduler can run on. This is incase you lose connection to the cluster mid run. 
5.Do $./scheduler.sh to start script. This script automatically submits jobs and lets you see the status of the runs. This will run for several hours.
5.The scheduler makes a text file of all jobs called jobList.txt in the ShelfMC dir. Make sure to delete jobList.txt before starting a whole new run.


######RESULT#######
1.When Completed, there will be a great amount of data in the run files, about 460GB.  
2.The run directory is organized in tree, results for particular runs can be found by cd'ing deeper into the tree.
3.In each run directory, there will be a resulting root file, all the setup files, and a log file for the run.
 

#!/bin/bash
#PBS -l walltime=04:00:00
#PBS -l nodes=1:ppn=1,mem=4000mb
#PBS -N jude_SetupJob
#PBS -j oe
#PBS -A PCON0003
#Jude Rajasekera 3/20/17
#directories
WorkDir=$TMPDIR   
tmpShelfmc=$HOME/shelfmc/ShelfMC #set your ShelfMC directory here

#controlled variables for run
runName='ParamSpaceScanDir' #name of run
NNU=500000 #NNU per run
seed=42 #starting seed for every run, each processor will recieve a different seed (42,43,44,45...)
NNU="$(($NNU / 20))" #calculating processors per node, change 20 to however many processors your cluster has per node
ppn=5 #number of processors per node on cluster
########################### input.txt file ####################################################
input1="#inputs for ARIANNA simulation, do not change order unless you change ReadInput()"
input2="$NNU #NNU, setting to 1 for unique neutrino"
input3="$seed      #seed Seed for Rand3"
input4="18.0    #EXPONENT, !should be exclusive with SPECTRUM"
input5="1000    #ATGap, m, distance between stations"
input6="4       #ST_TYPE, !restrict to 4 now!"
input7="4 #N_Ant_perST, not to be confused with ST_TYPE above"
input8="2 #N_Ant_Trigger, this is the minimum number of AT to trigger"
input9="30      #Z for ST_TYPE=2"
input10="$T   #ICETHICK, thickness of ice including firn, 575m at Moore's Bay"
input11="1       #FIRN, KD: ensure DEPTH_DEPENDENT is off if FIRN is 0"
input12="1.30    #NFIRN 1.30"
input13="$FT      #FIRNDEPTH in meters"
input14="1 #NROWS 12 initially, set to 3 for HEXAGONAL"
input15="1 #NCOLS 12 initially, set to 5 for HEXAGONAL"
input16="0       #SCATTER"
input17="1       #SCATTER_WIDTH,how many times wider after scattering"
input18="0       #SPECTRUM, use spectrum, ! was 1 initially!"
input19="0       #DIPOLE,  add a dipole to the station, useful for st_type=0 and 2"
input20="0       #CONST_ATTENLENGTH, use constant attenuation length if ==1"
input21="$L     #ATTEN_UP, this is the conjuction of the plot attenlength_up and attlength_down when setting REFLECT_RATE=0.5(3dB)"
input22="250     #ATTEN_DOWN, this is the average attenlength_down before Minna Bluff measurement(not used anymore except for CONST_ATTENLENGTH)"
input23="4 #NSIGMA, threshold of trigger"
input24="1      #ATTEN_FACTOR, change of the attenuation length"
input25="$Rval    #REFLECT_RATE,power reflection rate at the ice bottom"
input26="0       #GZK, 1 means using GZK flux, 0 means E-2 flux"
input27="0       #FANFLUX, use fenfang's flux which only covers from 10^17 eV to 10^20 eV"
input28="0       #WIDESPECTRUM, use 10^16 eV to 10^21.5 eV as the energy spectrum, otherwise use 17-20"
input29="1       #SHADOWING"
input30="1       #DEPTH_DEPENDENT_N;0 means uniform firn, 1 means n_firn is a function of depth"
input31="0 #HEXAGONAL"
input32="1       #SIGNAL_FLUCT 1=add noise fluctuation to signal or 0=do not"
input33="4.0     #GAINV  gain dependency"
input34="1       #TAUREGENERATION if 1=tau regeneration effect, if 0=original"
input35="3.0     #ST4_R radius in meters between center of station and antenna"
input36="350     #TNOISE noise temperature in Kelvin"
input37="80      #FREQ_LOW low frequency of LPDA Response MHz #was 100"
input38="1000    #FREQ_HIGH high frequency of LPDA Response MHz"
input39="/home/rajasekera.3/shelfmc/ShelfMC/temp/LP_gain_manual.txt     #GAINFILENAME"
input40="$SD     #STATION_DEPTH"
#######################################################################################################

cd $TMPDIR   
mkdir $runName
cd $runName

initSeed=$seed
counter=0
for L in {500..1000..100} #attenuation length 500-1000
do
    mkdir Atten_Up$L
    cd Atten_Up$L

    for R in {0..100..25} #Reflection Rate 0-1
    do
        mkdir ReflectionRate$R
        cd ReflectionRate$R
        if [ "$R" = "100" ]; then #fixing reflection rate value
            Rval="1.0"
        else
            Rval="0.$R"
        fi

        for T in {500..2900..400} #Thickness of Ice 500-2900
        do
            mkdir IceThick$T
            cd IceThick$T
            for FT in {60..140..20} #Firn Thinckness 60-140
            do
                mkdir FirnThick$FT
                cd FirnThick$FT
                for SD in {0..200..50} #Station Depth
                do
                    mkdir StationDepth$SD
                    cd StationDepth$SD
                    #####Do file operations###########################################
                    counter=$((counter+1))
                    echo "Counter = $counter ; L = $L ; R = $Rval ; T = $T ; FT = $FT ; SD = $SD " #print variables

                    #define changing lines
                    input21="$L     #ATTEN_UP, this is the conjuction of the plot attenlength_up and attlength_down when setting REFLECT_RATE=0.5(3dB)"
                    input25="$Rval    #REFLECT_RATE,power reflection rate at the ice bottom"
                    input10="$T   #ICETHICK, thickness of ice including firn, 575m at Moore's Bay"
                    input13="$FT      #FIRNDEPTH in meters"
                    input40="$SD       #STATION_DEPTH"
		    
		    for (( i=1; i<=$ppn;i++)) #make 20 setup files for 20 processors
                    do

                        mkdir Setup$i #make setup folder
                        cd Setup$i #go into setup folder
                        seed="$(($initSeed + $i -1))" #calculate seed for this iteration
                        input3="$seed      #seed Seed for Rand3"

                        for j in {1..40} #print all input.txt lines
                        do
                            lineName=input$j
                            echo "${!lineName}" >> input.txt
                        done
			
                        cd ..
                    done
		    
		    pwd=`pwd`
                    #create job file
		    echo '#!/bin/bash' >> run_shelfmc_multithread.sh
		    echo '#PBS -l nodes=1:ppn='$ppn >> run_shelfmc_multithread.sh
		    echo '#PBS -l walltime=00:05:00' >> run_shelfmc_multithread.sh #change walltime as necessary
		    echo '#PBS -N jude_'$runName'_job' >> run_shelfmc_multithread.sh #change job name as necessary
		    echo '#PBS -j oe'  >> run_shelfmc_multithread.sh
		    echo '#PBS -A PCON0003' >> run_shelfmc_multithread.sh #change group if necessary
		    echo 'cd ' $tmpShelfmc >> run_shelfmc_multithread.sh
		    echo 'runName='$runName  >> run_shelfmc_multithread.sh
		    for (( i=1; i<=$ppn;i++))
		    do
			echo './shelfmc_stripped.exe $runName/'Atten_Up$L'/'ReflectionRate$R'/'IceThick$T'/'FirnThick$FT'/'StationDepth$SD'/Setup'$i' _'$i'$runName &' >> run_shelfmc_multithread.sh
		    done
		   # echo './shelfmc_stripped.exe $runName/'Atten_Up$L'/'ReflectionRate$R'/'IceThick$T'/'FirnThick$FT'/'StationDepth$SD'/Setup1 _01$runName &' >> run_shelfmc_multithread.sh
		    echo 'wait' >> run_shelfmc_multithread.sh
		    echo 'cd $runName/'Atten_Up$L'/'ReflectionRate$R'/'IceThick$T'/'FirnThick$FT'/'StationDepth$SD >> run_shelfmc_multithread.sh
		    echo 'for (( i=1; i<='$ppn';i++)) #20 iterations' >> run_shelfmc_multithread.sh
		    echo 'do' >> run_shelfmc_multithread.sh
		    echo '  cd Setup$i #cd into setup dir' >> run_shelfmc_multithread.sh
		    echo '  mv *.root ..' >> run_shelfmc_multithread.sh
		    echo '  cd ..' >> run_shelfmc_multithread.sh
		    echo 'done' >> run_shelfmc_multithread.sh
		    echo 'hadd Result_'$runName'.root *.root' >> run_shelfmc_multithread.sh
		    echo 'rm *ShelfMCTrees*' >> run_shelfmc_multithread.sh

		    chmod u+x run_shelfmc_multithread.sh # make executable

                    ##################################################################
                    cd ..
                done
                cd ..
            done
            cd ..
        done
        cd ..
    done
    cd ..
done
cd 

mv $WorkDir/$runName $HOME

#!/bin/bash
#Jude Rajasekera 3/20/17

tmpShelfmc=$HOME/shelfmc/ShelfMC #location of Shelfmc
runName=ParamSpaceScanDir #name of run

cd $tmpShelfmc #move to home directory

if [ ! -f ./jobList.txt ]; then #see if there is an existing job file
    echo "Creating new job List"
    for L in {500..1000..100} #attenuation length 500-1000
    do
	for R in {0..100..25} #Reflection Rate 0-1
	do
            for T in {500..2900..400} #Thickness of Ice 500-2900
            do
		for FT in {60..140..20} #Firn Thinckness 60-140
		do
                    for SD in {0..200..50} #Station Depth
                    do
		    echo "cd $runName/Atten_Up$L/ReflectionRate$R/IceThick$T/FirnThick$FT/StationDepth$SD" >> jobList.txt
                    done
		done
            done
	done
    done
else 
    echo "Picking up from last job"
fi


numbLeft=$(wc -l < ./jobList.txt)
while [ $numbLeft -gt 0 ];
do
    jobs=$(showq | grep "rajasekera.3") #change username here
    echo '__________Current Running Jobs__________'
    echo "$jobs"
    echo ''
    runningJobs=$(showq | grep "rajasekera.3" | wc -l) #change unsername here
    echo Number of Running Jobs = $runningJobs 
    echo Number of jobs left = $numbLeft
    if [ $runningJobs -le 6 ];then
	line=$(head -n 1 jobList.txt)
	$line
	echo Submit Job && pwd
	qsub run_shelfmc_multithread.sh
	cd $tmpShelfmc
	sed -i 1d jobList.txt
    else
	echo "Full Capacity"
    fi
    sleep 1
    numbLeft=$(wc -l < ./jobList.txt)
done

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* &

//////////////////////////////////////
//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;
}

#############################################################################
## 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)
#############################################################################

# -*- 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()
        
        

#include <iostream>
#include <thread>
#include <chrono>
#include <cmath>


using namespace std;

void nap(){
  // usleep(3000);
  // this_thread::sleep_for(30000ms);

  for (int i=0;i<1000000000;i++){
    double j = sqrt(i^2);
  }
}

int main(){

  cout << "taking a nap" << endl;

  nap();

  cout << "hello world" << endl;

}

SHELL=/bin/sh

# Note: Comments start with #.  $(FOOBAR) means: evaluate the variable 
#        defined by FOOBAR= (something).

# This file contains a set of rules used by the "make" command.
#   This makefile $(MAKEFILE) tells "make" how the executable $(COMMAND) 
#   should be create from the source files $(SRCS) and the header files 
#   $(HDRS) via the object files $(OBJS); type the command:
#        "make -f make_program"
#   where make_program should be replaced by the name of the makefile.
# 
# Programmer:  Dick Furnstahl (furnstahl.1@osu.edu)
# Latest revision: 12-Jan-2016 
# 
# Notes:
#  * If you are ok with the default options for compiling and linking, you
#     only need to change the entries in section 1.
#
#  * Defining BASE determines the name for the makefile (prepend "make_"), 
#     executable (append ".x"), zip archive (append ".zip") and gzipped 
#     tar file (append ".tar.gz"). 
#
#  * To remove the executable and object files, type the command:
#          "make -f $(MAKEFILE) clean"
#
#  * To create a zip archive with name $(BASE).zip containing this 
#     makefile and the SRCS and HDRS files, type the command:
#        "make -f $(MAKEFILE) zip"
#
#  * To create a gzipped tar file with name $(BASE).tar.gz containing this 
#     makefile and the source and header files, type the command:
#          "make -f $(MAKEFILE) tarz"
#
#  * Continuation lines are indicated by \ with no space after it.  
#     If you get a "missing separator" error, it is probably because there
#     is a space after a \ somewhere.
#

###########################################################################
# 1. Specify base name, source files, header files, input files
########################################################################### 

# The base for the names of the makefile, executable command, etc.
BASE= hello_world

# Put all C++ (or other) source files here.  NO SPACES after continuation \'s.
SRCS= \
hello_world.cpp

# Put all header files here.  NO SPACES after continuation \'s.
HDRS= \

# Put any input files you want to be saved in tarballs (e.g., sample files).
INPFILE= \

###########################################################################
# 2. Generate names for object files, makefile, command to execute, tar file
########################################################################### 

# *** YOU should not edit these lines unless to change naming conventions ***

OBJS= $(addsuffix .o, $(basename $(SRCS)))
MAKEFILE= make_$(BASE)
COMMAND=  $(BASE).x
TARFILE= $(BASE).tar.gz
ZIPFILE= $(BASE).zip

###########################################################################
# 3. Commands and options for different compilers
########################################################################### 

#
# Compiler parameters
#
# CXX           Name of the C++ compiler to use
# CFLAGS        Flags to the C++ compiler
# CWARNS        Warning options for C++ compiler
# F90           Name of the fortran compiler to use (if relevant) 
# FFLAGS        Flags to the fortran compiler 
# LDFLAGS       Flags to the loader
# LIBS          A list of libraries 
#

CXX= g++
CFLAGS=  -g
CWARNS= -Wall -W -Wshadow -fno-common 
MOREFLAGS= -Wpedantic -Wpointer-arith -Wcast-qual -Wcast-align \
           -Wwrite-strings -fshort-enums 

# add relevant libraries and link options
LIBS=           
# LDFLAGS= -lgsl -lgslcblas 
LDFLAGS=

###########################################################################
# 4. Instructions to compile and link, with dependencies
########################################################################### 
all:    $(COMMAND) 

.SUFFIXES:
.SUFFIXES: .o .mod .f90 .f .cpp

#%.o:   %.mod 

# This is the command to link all of the object files together. 
#  For fortran, replace CXX by F90.
$(COMMAND): $(OBJS) $(MAKEFILE) 
	$(CXX) -o $(COMMAND) $(OBJS) $(LDFLAGS) $(LIBS)

# Command to make object (.o) files from C++ source files (assumed to be .cpp).
#  Add $(MOREFLAGS) if you want additional warning options.
%.o: %.cpp $(HDRS) $(MAKEFILE)
	$(CXX) -c $(CFLAGS) $(CWARNS) -o $@ $<

# Commands to make object (.o) files from Fortran-90 (or beyond) and
#  Fortran-77 source files (.f90 and .f, respectively).
.f90.mod:
	$(F90) -c $(F90FLAGS) -o $@ $< 
 
.f90.o: 
	$(F90) -c $(F90FLAGS) -o $@ $<
 
.f.o:   
	$(F90) -c $(FFLAGS) -o $@ $<
      
##########################################################################
# 5. Additional tasks      
##########################################################################
      
# Delete the program and the object files (and any module files)
clean:
	/bin/rm -f $(COMMAND) $(OBJS)
	/bin/rm -f $(MODIR)/*.mod
 
# Pack up the code in a compressed gnu tar file 
tarz:
	tar cfvz $(TARFILE) $(MAKEFILE) $(SRCS) $(HDRS) $(MODIR) $(INPFILE) 

# Pack up the code in a zip archive
zip:
	zip -r $(ZIPFILE) $(MAKEFILE) $(SRCS) $(HDRS) $(MODIR) $(INPFILE) 

##########################################################################
# That's all, folks!     
##########################################################################

#include "TF1.h"

void fitSnr();

void fitSnr()

{
	gStyle->SetLineWidth(4); //set some style parameters
	TFile *stripe1file = new TFile("hstripe1snr.root"); //import the file containing the histogram to be fit
	TH1D *hstripe1 = (TH1D*)stripe1file->Get("stripe1snr"); //get the histogram to be fit
	TCanvas c1("c1","c1",1000,800); //make a canvas
	hstripe1->Draw(""); //draw it
	c1.SetLogy(); //set a log axis

	//need to declare an equation
	//I want to fit for two paramters, in the equation these are [0] and [1]
	//so, you'll need to re-write the equation to whatever you're trying to fit for
	//but ROOT wants the variables to find to be given as [0], [1], [2], etc.
	
	char equation1[150]; //declare a container for the equation
	sprintf(equation1,"([0]*(x^[1]))"); //declare the equation
	TF1 *fit1 = new TF1("PowerFit",equation1,20,50); //create a function to fit with, with the range being 20 to 50

	//now, we need to set the initial parameters of the fit
	//fit->SetParameter(0,H->GetRMS()); //this should be a good starting place for a standard deviation like variable
	//fit->SetParameter(1,H->GetMaximum()); //this should be a good starting place for amplitude like variable
	fit1->SetParameter(0,60000.); //for our example, we will manually choose this
	fit1->SetParameter(1,-3.);

	hstripe1->Fit("PowerFit","R"); //actually do the fit;
	fit1->Draw("same"); //draw the fit

	//now, we want to print out some parameters to see how good the fit was
	cout << "par0 " << fit1->GetParameter(0) << " par1 " << fit1->GetParameter(1) << endl;
	cout<<"chisquare "<<fit1->GetChisquare()<<endl;
	cout<<"Ndf "<<fit1->GetNDF()<<endl;
	cout<<"reduced chisquare "<<double(fit1->GetChisquare())/double(fit1->GetNDF())<<endl;
	cout<<"   "<<endl;
	c1.SaveAs("stripe1snrfit.png");

}

export ROOTSYS=/users/PAS0174/osu8620/root-6.08.06
eval 'source /users/PAS0174/osu8620/root-6.08.06/builddir/bin/thisroot.sh'
export LD_INCLUDE_PATH=/users/PAS0174/osu8620/cint/libcint/build/include:$LD_INCLUDE_PATH
module load fftw3/3.3.5
module load gnu/6.3.0
module load python/3.4.2
module load cmake/3.7.2

#might need this, but probably not
#export CC=/usr/local/gcc/6.3.0/bin/gcc

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.

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////		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;		
	}	
}

#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)
#############################################################################

///////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
///////		demo.cxx
//////		Nov 2018, Brian Clark 
//////		Demonstrate how to draw 2 TH2D's with two different color palettes (needs ROOT6)
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

//C++
#include <iostream>

//ROOT Includes
#include "TH2D.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TExec.h"
#include "TColor.h"
#include "TPaletteAxis.h"
#include "TRandom3.h"

using namespace std;

int main(int argc, char **argv)
{

	gStyle->SetOptStat(0);

	TH2D *h1 = new TH2D("h1","h1",100,-10,10,100,-10,10);
	TH2D *h2 = new TH2D("h2","h2",100,-10,10,100,-10,10);

	TRandom3 *R = new TRandom3(time(0));
	int count=0;
	while(count<10000000){
		h1->Fill(R->Gaus(-4.,1.),R->Gaus(-4.,1.));
		count++;
	}

	TRandom3 *R2 = new TRandom3(time(0)+100);
	count=0;
	while(count<10000000){
		h2->Fill(R2->Gaus(4.,1.),R2->Gaus(4.,1.));
		count++;
	}

	TCanvas *c = new TCanvas("c","c",1100,850);

	h1->Draw("colz");
		TExec *ex1 = new TExec("ex1","gStyle->SetPalette(kValentine);");
		ex1->Draw();
		h1->Draw("colz same");
	gPad->SetLogz();
	gPad->Modified();
	gPad->Update();

	TPaletteAxis *palette1 = (TPaletteAxis*)h1->GetListOfFunctions()->FindObject("palette");
	palette1->SetX1NDC(0.7);
	palette1->SetX2NDC(0.75);
	palette1->SetY1NDC(0.1);
	palette1->SetY2NDC(0.95);

	gPad->SetRightMargin(0.3);
	gPad->SetTopMargin(0.05);

	h2->Draw("same colz");
		TExec *ex2 = new TExec("ex2","gStyle->SetPalette(kCopper);"); //TColor::InvertPalette();
		ex2->Draw();
		h2->Draw("same colz");
	gPad->SetLogz();
	gPad->Modified();
	gPad->Update();

	TPaletteAxis *palette2 = (TPaletteAxis*)h2->GetListOfFunctions()->FindObject("palette");
	palette2->SetX1NDC(0.85);
	palette2->SetX2NDC(0.90);
	palette2->SetY1NDC(0.1);
	palette2->SetY2NDC(0.95);

	h1->SetTitle("");
	h1->GetXaxis()->SetTitle("X-Value");
	h1->GetYaxis()->SetTitle("Y-Value");

	h1->GetZaxis()->SetTitle("First Histogram Events");
	h2->GetZaxis()->SetTitle("Second Histogram Events");

	h1->GetYaxis()->SetTitleSize(0.05);
	h1->GetXaxis()->SetTitleSize(0.05);
	h1->GetYaxis()->SetLabelSize(0.045);
	h1->GetXaxis()->SetLabelSize(0.045);

	h1->GetZaxis()->SetTitleSize(0.045);
	h2->GetZaxis()->SetTitleSize(0.045);
	h1->GetZaxis()->SetLabelSize(0.04);
	h2->GetZaxis()->SetLabelSize(0.04);
	
	h1->GetYaxis()->SetTitleOffset(1);
	h1->GetZaxis()->SetTitleOffset(1);
	h2->GetZaxis()->SetTitleOffset(1);

	c->SetLogz();
	c->SaveAs("test.png");

}

LDFLAGS=-L${ARA_UTIL_INSTALL_DIR}/lib -L${shell root-config --libdir}
CXXFLAGS=-I${ARA_UTIL_INSTALL_DIR}/include -I${shell root-config --incdir}
LDLIBS += $(shell $(ROOTSYS)/bin/root-config --libs)

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////  demo.cxx 
////  demo
////
////  Nov 2018
////////////////////////////////////////////////////////////////////////////////

//Includes
#include <iostream>
#include <string>
#include <sstream>

//AraRoot Includes
#include "RawAtriStationEvent.h"
#include "UsefulAtriStationEvent.h"

//ROOT Includes
#include "TTree.h"
#include "TFile.h"
#include "TGraph.h"

using namespace std;

RawAtriStationEvent *rawAtriEvPtr;

int main(int argc, char **argv)
{

	if(argc<3) {
		std::cout << "Usage\n" << argv[0] << " <input_file> <output_location> "<<endl;
		return -1;
	}

	/*
	arguments
	0: exec
	1: input data file
	2: output location
	*/
	
	TFile *fpIn = TFile::Open(argv[1]);
	if(!fpIn) {
		std::cout << "Can't open file\n";
		return -1;
	}
	TTree *eventTree = (TTree*) fpIn->Get("eventTree");
	if(!eventTree) {
		std::cout << "Can't find eventTree\n";
		return -1;
	}
	eventTree->SetBranchAddress("event",&rawAtriEvPtr);
	int run;
	eventTree->SetBranchAddress("run",&run);
	eventTree->GetEntry(0);
	printf("Filter Run Number %d \n", run);

	char outfile_name[400];
	sprintf(outfile_name,"%s/outputs_run%d.root",argv[2],run);

	TFile *fpOut = TFile::Open(outfile_name, "RECREATE");
	TTree* outTree = new TTree("outTree", "outTree");
	int WaveformLength[16];
	outTree->Branch("WaveformLength", &WaveformLength, "WaveformLength[16]/D");
	
	Long64_t numEntries=eventTree->GetEntries();

	for(Long64_t event=0;event<numEntries;event++) {
		eventTree->GetEntry(event);
		UsefulAtriStationEvent *realAtriEvPtr = new UsefulAtriStationEvent(rawAtriEvPtr, AraCalType::kLatestCalib);
		for(int i=0; i<16; i++){
			TGraph *gr = realAtriEvPtr->getGraphFromRFChan(i);
			WaveformLength[i] = gr->GetN();
			delete gr;
		}		
		outTree->Fill();
		delete realAtriEvPtr;
	} //loop over events
	
	fpOut->Write();
	fpOut->Close();
	
	fpIn->Close();
	delete fpIn;
}

#/bin/bash
#PBS -l nodes=1:ppn=1
#PBS -l mem=4GB
#PBS -l walltime=00:05:00
#PBS -A PAS0654
#PBS -e /fs/scratch/PAS0654/shell_demo/err_out_logs
#PBS -o /fs/scratch/PAS0654/shell_demo/err_out_logs

# you should change the -e and -o to write your 
# log files to a location of your preference

# source your own shell script here
eval 'source /users/PAS0654/osu0673/A23_analysis/env.sh'

# $RUNDIR was defined in the submission script 
# along with $FILE and $OUTPUTDIR

cd $RUNDIR

# $TMPDIR is the local memory of this specific node
# it's the only variable we didn't have to define

./bin/demo $FILE $TMPDIR 

# after we're done
# we copy the results to the $OUTPUTDIR

pbsdcp $TMPDIR/'*.root' $OUTPUTDIR

/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1000.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1001.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1002.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1004.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1005.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1006.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1007.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1009.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1010.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1011.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1012.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1014.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1015.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1016.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1017.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1019.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1020.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1021.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1022.root
/fs/scratch/PAS0654/ara/10pct/RawData/A3/2013/sym_links/event1024.root

#!/bin/bash

#where should the outputs be stored?
OutputDir="/fs/scratch/PAS0654/shell_demo/outputs"
echo '[ Processed file output directory: ' $OutputDir ' ]'
export OutputDir

#where is your executable compiled?
RunDir="/users/PAS0654/osu0673/A23_analysis/araROOT"
export RunDir

#define the list of runs to execute on
readfile=run_list.txt

counter=0
while read line1
do
	qsub -v RUNDIR=$RunDir,OUTPUTDIR=$OutputDir,FILE=$line1 -N 'job_'$counter run.sh
	counter=$((counter+1))
done < $readfile

# .bash_profile

# Get the aliases and functions
if [ -f ~/.bashrc ]; then
        . ~/.bashrc
fi

# User specific environment and startup programs

PATH=$PATH:$HOME/.local/bin:$HOME/bin

export PATH

# .bashrc

source bashrc_anita.sh

# we also want to set two more environment variables that ANITA needs
# you should update ICEMC_SRC_DIR and ICEMC_BUILD_DIR to wherever you
# downloaded icemc too

export ICEMC_SRC_DIR=/path/to/icemc #change this line!
export ICEMC_BUILD_DIR=/path/to/icemc #change this line!
export DYLD_LIBRARY_PATH=${ICEMC_SRC_DIR}:${ICEMC_BUILD_DIR}:${DYLD_LIBRARY_PATH}

//C++ includes
#include <iostream>

//ROOT includes
#include "TCanvas.h"
#include "TStyle.h"
#include "TH1D.h"
#include "TFile.h"
#include "TTree.h"

using namespace std;

int main(int argc, char *argv[])
{

	if(argc<2)
	{
		cout << "Not enough arguments! Stop run. " << endl;
		return -1;
	}

	/*
	we're going to make a histogram, and set some parameters about it's X and Y axes
	*/
	TH1D *nuflavorint_hist = new TH1D("nuflavorint", "",3,1,4);  
	nuflavorint_hist->SetTitle("Neutrino Flavors");
	nuflavorint_hist->GetXaxis()->SetTitle("Neutrino Flavors (1=e, 2=muon, 3=tau)");
	nuflavorint_hist->GetYaxis()->SetTitle("Weigthed Fraction of Total Detected Events");
	nuflavorint_hist->GetXaxis()->SetTitleOffset(1.2);
	nuflavorint_hist->GetYaxis()->SetTitleOffset(1.2);
	nuflavorint_hist->GetXaxis()->CenterTitle();
	nuflavorint_hist->GetYaxis()->CenterTitle();

	for(int i=1; i < argc; i++)
		
	{  // loop over the input files

		//now we are going to load the icefinal.root file and draw in the "passing_events" tree, which stores info
		
		string readfile = string(argv[i]);
		TFile *AnitaFile = new TFile(( readfile ).c_str());
		cout << "AnitaFile" << endl;
		TTree *passing_events = (TTree*)AnitaFile->Get("passing_events");
		cout << "Reading AnitaFile..." << endl;

		//declare three variables we are going to use later
		
		int num_pass;               // number of entries (ultra-neutrinos);
		double weight;              // weight of neutrino counts;
		int nuflavorint;              // neutrino flavors;
		
		num_pass = passing_events->GetEntries();
		cout << "num_pass is " << num_pass << endl;
		
		/*PRIMARIES VARIABLES*/

		//set the "branch" of the tree which stores specific pieces of information
		
		passing_events->SetBranchAddress("weight", &weight);
		passing_events->SetBranchAddress("nuflavor", &nuflavorint);

		//loop over all the events in the tree
	
		for (int k=0; k <=num_pass; k++)
		{
			passing_events->GetEvent(k);
			nuflavorint_hist->Fill(nuflavorint, weight); //fill the histogram with this value and this weight
		
		} // CLOSE FOR LOOP OVER NUMBER OF EVENTS
		
	} // CLOSE FOR LOOP OVER NUMBER OF INPUT FILES

	//set up some parameters to make things loo pretty
	gStyle->SetHistFillColor(0);
	gStyle->SetHistFillStyle(1);
	gStyle->SetHistLineColor(1);
	gStyle->SetHistLineStyle(0);
	gStyle->SetHistLineWidth(2.5); //Setup plot Style

	//make a "canvas" to draw on

	TCanvas *c4 = new TCanvas("c4", "nuflavorint", 1100,850);
	gStyle->SetOptTitle(1);
	gStyle->SetStatX(0.33);
	gStyle->SetStatY(0.87);
	nuflavorint_hist->Draw("HIST"); //draw on it
			
	//Save Plots
	
	//make the line thicker and then save the result
	gStyle->SetHistLineWidth(9);
	c4->SaveAs("nuflavorint.png");
	gStyle->SetHistLineWidth(2);
	c4->SaveAs("nuflavorint.pdf");

	delete c4; //clean up
	return 0; //return successfully
	
}

# .bashrc

# Source global definitions
if [ -f /etc/bashrc ]; then
	. /etc/bashrc
fi

module load cmake/3.11.4
module load gnu/7.3.0
export CC=`which gcc`
export CXX=`which g++`

export BOOST_ROOT=/fs/project/PAS0654/shared_software/anita/owens_pitzer/build/boost_build
export LD_LIBRARY_PATH=${BOOST_ROOT}/stage/lib:$LD_LIBRARY_PATH
export BOOST_LIB=$BOOST_ROOT/stage/lib
export LD_LIBRARY_PATH=$BOOST_LIB:$LD_LIBRARY_PATH

export ROOTSYS=/fs/project/PAS0654/shared_software/anita/owens_pitzer/build/root
eval 'source /fs/project/PAS0654/shared_software/anita/owens_pitzer/build/root/bin/thisroot.sh'

# Makefile for the ROOT test programs.  # This Makefile shows nicely how to compile and link applications
# using the ROOT libraries on all supported platforms.
#
# Copyright (c) 2000 Rene Brun and Fons Rademakers
#
# Author: Fons Rademakers, 29/2/2000

include Makefile.arch

################################################################################
# Site specific flags
################################################################################
# Toggle these as needed to get things to install

#BOOSTFLAGS = -I boost_1_48_0
# commented out for kingbee and older versions of gcc
ANITA3_EVENTREADER=1

# Uncomment to enable healpix 
#USE_HEALPIX=1

# Uncomment to disable explicit vectorization (but will do nothing if ANITA_UTIL is not available) 
#VECTORIZE=1


# The ROOT flags are added to the CXXFLAGS in the .arch file
# so this should be simpler...
ifeq (,$(findstring -std=, $(CXXFLAGS)))
	ifeq ($(shell test $(GCC_MAJOR) -lt 5; echo $$?),0)
		ifeq ($(shell test $(GCC_MINOR) -lt 5; echo $$?),0)
			CXXFLAGS += -std=c++0x
		else
			CXXFLAGS += -std=c++11
		endif
	endif
endif

################################################################################

# If not compiling with C++11 (or later) support, all occurrences of "constexpr"
# must be replaced with "const", because "constexpr" is a keyword
# which pre-C++11 compilers do not support.
# ("constexpr" is needed in the code to perform in-class initialization
# of static non-integral member objects, i.e.:
#		static const double c_light = 2.99e8;
# which works in C++03 compilers, must be modified to:
#		static constexpr double c_light = 2.99e8;
# to work in C++11, but adding "constexpr" breaks C++03 compatibility.
# The following compiler flag defines a preprocessor macro which is
# simply:
#		#define constexpr const
# which replaces all instances of the text "constexpr" and replaces it
# with "const".
# This preserves functionality while only affecting very specific semantics.

ifeq (,$(findstring -std=c++1, $(CXXFLAGS)))
	CPPSTD_FLAGS = -Dconstexpr=const
endif



# Uses the standard ANITA environment variable to figure
# out if ANITA libs are installed
ifdef ANITA_UTIL_INSTALL_DIR
	ANITA_UTIL_EXISTS=1
	ANITA_UTIL_LIB_DIR=${ANITA_UTIL_INSTALL_DIR}/lib
	ANITA_UTIL_INC_DIR=${ANITA_UTIL_INSTALL_DIR}/include
	LD_ANITA_UTIL=-L$(ANITA_UTIL_LIB_DIR)
	LIBS_ANITA_UTIL=-lAnitaEvent -lRootFftwWrapper
	INC_ANITA_UTIL=-I$(ANITA_UTIL_INC_DIR)
	ANITA_UTIL_ETC_DIR=$(ANITA_UTIL_INSTALL_DIR)/etc
endif

ifdef ANITA_UTIL_EXISTS
	CXXFLAGS += -DANITA_UTIL_EXISTS
endif

ifdef VECTORIZE
	CXXFLAGS += -DVECTORIZE -march=native -fabi-version=0
endif

ifdef ANITA3_EVENTREADER
	CXXFLAGS += -DANITA3_EVENTREADER
endif

ifdef USE_HEALPIX
	CXXFLAGS += -DUSE_HEALPIX `pkg-config --cflags healpix_cxx`
	LIBS  += `pkg-config --libs healpix_cxx` 
endif


################################################################################

GENERAL_FLAGS = -g -O2 -pipe -m64 -pthread
WARN_FLAGS = -W -Wall -Wextra -Woverloaded-virtual
# -Wno-unused-variable -Wno-unused-parameter -Wno-unused-but-set-variable

CXXFLAGS += $(GENERAL_FLAGS) $(CPPSTD_FLAGS) $(WARN_FLAGS) $(ROOTCFLAGS) $(INC_ANITA_UTIL)

DBGFLAGS  = -pipe -Wall -W -Woverloaded-virtual -g -ggdb -O0 -fno-inline

DBGCXXFLAGS = $(DBGFLAGS) $(ROOTCFLAGS) $(BOOSTFLAGS)
LDFLAGS  += $(CPPSTD_FLAGS) $(LD_ANITA_UTIL) -I$(BOOST_ROOT) -L.
LIBS += $(LIBS_ANITA_UTIL)

# Mathmore not included in the standard ROOT libs
LIBS += -lMathMore

DICT = classdict

OBJS = vector.o position.o earthmodel.o balloon.o icemodel.o signal.o ray.o Spectra.o anita.o roughness.o secondaries.o Primaries.o Tools.o counting.o $(DICT).o Settings.o Taumodel.o screen.o GlobalTrigger.o ChanTrigger.o SimulatedSignal.o EnvironmentVariable.o source.o  random.o

BINARIES = test_plot$(ExeSuf)

################################################################################

.SUFFIXES: .$(SrcSuf) .$(ObjSuf) .$(DllSuf)

all:            $(BINARIES)

$(BINARIES): %: %.$(SrcSuf) $(OBJS)
		$(LD) $(CXXFLAGS) $(LDFLAGS) $(OBJS) $< $(LIBS) $(OutPutOpt) $@
		@echo "$@ done"


.PHONY: clean
clean:
		@rm -f $(BINARIES)


%.$(ObjSuf) : %.$(SrcSuf) %.h
	@echo "<**Compiling**> "$<
	$(LD) $(CXXFLAGS) -c $< -o $@

#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

# .bashrc

# Source global definitions
if [ -f /etc/bashrc ]; then
	. /etc/bashrc
fi

module load cmake/3.11.4
module load gnu/7.3.0
export CC=`which gcc`
export CXX=`which g++`

export BOOST_ROOT=/fs/project/PAS0654/shared_software/anita/owens_pitzer/build/boost_build
export LD_LIBRARY_PATH=${BOOST_ROOT}/stage/lib:$LD_LIBRARY_PATH
export BOOST_LIB=$BOOST_ROOT/stage/lib
export LD_LIBRARY_PATH=$BOOST_LIB:$LD_LIBRARY_PATH

export ROOTSYS=/fs/project/PAS0654/shared_software/anita/owens_pitzer/build/root
eval 'source /fs/project/PAS0654/shared_software/anita/owens_pitzer/build/root/bin/thisroot.sh'

# .bash_profile

# Get the aliases and functions
if [ -f ~/.bashrc ]; then
        . ~/.bashrc
fi

# User specific environment and startup programs

PATH=$PATH:$HOME/.local/bin:$HOME/bin

export PATH

# .bashrc

source bashrc_anita.sh

# we also want to set two more environment variables that ANITA needs
# you should update ICEMC_SRC_DIR and ICEMC_BUILD_DIR to wherever you
# downloaded icemc too

export ICEMC_SRC_DIR=/path/to/icemc #change this line!
export ICEMC_BUILD_DIR=/path/to/icemc #change this line!
export DYLD_LIBRARY_PATH=${ICEMC_SRC_DIR}:${ICEMC_BUILD_DIR}:${DYLD_LIBRARY_PATH}

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bcdfb138",
   "metadata": {},
   "source": [
    "# IC/ARA Coincident Simulation Events Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7ad7c80",
   "metadata": {},
   "source": [
    "### Settings and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b6915a86",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:75% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Makes this notebook maximally wide\n",
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:75% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6ef9e20b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Author: Alex Machtay (machtay.1@osu.edu)\n",
    "## Modified by: Alan Salcedo (salcedogomez.1@osu.edu)\n",
    "## Date: 4/26/23\n",
    "\n",
    "## Purpose:\n",
    "### This script will read the data files produced by AraRun_corrected_MultStat.py to make histograms relevant plots of\n",
    "### neutrino events passing through icecube and detected by ARA station (in AraSim)\n",
    "\n",
    "\n",
    "## Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "sys.path.append(\"/users/PAS0654/osu8354/root6_18_build/lib\") # go to parent dir\n",
    "sys.path.append(\"/users/PCON0003/cond0068/.local/lib/python3.6/site-packages\")\n",
    "import math\n",
    "import argparse\n",
    "import glob\n",
    "import pandas as pd\n",
    "pd.options.mode.chained_assignment = None  # default='warn'\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import jupyterthemes as jt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f24b8292",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/bin/bash: jt: command not found\r\n"
     ]
    }
   ],
   "source": [
    "## Set style for the jupyter notebook\n",
    "!jt -t grade3 -T -N -kl -lineh 160 -f code -fs 14 -ofs 14 -cursc o\n",
    "jt.jtplot.style('grade3', gridlines='')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5a812d6",
   "metadata": {},
   "source": [
    "### Set constants\n",
    "#### These are things like the position of ARA station holes, the South Pole, IceCube's position in ARA station's coordinates, and IceCube's radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cdb851de",
   "metadata": {},
   "outputs": [],
   "source": [
    "## What IceCube (IC) station are you analyzing\n",
    "station = 1\n",
    "\n",
    "## What's the radius around ARA where neutrinos were injected\n",
    "inj_rad = 5 #in km\n",
    "\n",
    "## IceCube's center relative to each ARA station\n",
    "IceCube = [[-1128.08, -2089.42, -1942.39], [-335.812, -3929.26, -1938.23],\n",
    "          [-2320.67, -3695.78, -1937.35], [-3153.04, -1856.05, -1942.81], [472.49, -5732.98, -1922.06]] #IceCube's position relative to A1, A2, or A3\n",
    "\n",
    "#To calculate this, we need to do some coordinate transformations. Refer to this notebook to see the calculations: \n",
    "# IceCube_Relative_to_ARA_Stations.ipynb (found here - http://radiorm.physics.ohio-state.edu/elog/How-To/48) \n",
    "\n",
    "## IceCube's radius\n",
    "IceCube_radius = 564.189583548 #Modelling IceCube as a cylinder, we find the radius with V = h * pi*r^2 with V = 1x10^9 m^3 and h = 1x10^3 m "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af5c6fc2",
   "metadata": {},
   "source": [
    "### Read the data\n",
    "\n",
    "#### Once we import the data, we'll make dataframes to concatenate it and make some calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "921f61f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import data files\n",
    "\n",
    "#Here, it's from OSC Connolly's group project space\n",
    "source = '/fs/project/PAS0654/IceCube_ARA_Coincident_Search/AraSim/outputs/Coincident_Search_Runs/20M_GZK_5km_S1_correct' \n",
    "num_files = 200  # Number of files to read in from the source directory\n",
    "\n",
    "## Make a list of all of the paths to check \n",
    "file_list = []\n",
    "for i in range(1, num_files + 1):\n",
    "        for name in glob.glob(source + \"/\" + str(i) + \"/*.csv\"):\n",
    "                file_list.append(str(name))\n",
    "                #file_list gets paths to .csv files\n",
    "                \n",
    "## Now read the csv files into a pandas dataframe\n",
    "dfs = []\n",
    "for filename in file_list:\n",
    "        df = pd.read_csv(filename, index_col=None, header=0) #Store each csv file into a pandas data frame\n",
    "        dfs.append(df) #Append the csv file to store all of them in one\n",
    "frame = pd.concat(dfs, axis=0, ignore_index = True) #Concatenate pandas dataframes "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "813fb3ee",
   "metadata": {},
   "source": [
    "### Work with the data\n",
    "\n",
    "#### All the data from our coincidence simulations (made by AraRun_MultStat.sh) is now stored in a pandas data frame that we can work with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8b449583",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19800000\n",
      "19799802\n"
     ]
    }
   ],
   "source": [
    "## Now let's clean up our data and calculate other relevant things\n",
    "print(len(frame))\n",
    "frame = frame[frame['weight'].between(0,1)] #Filter out events with ill-defined weight (should be between 0 and 1)\n",
    "print(len(frame))\n",
    "\n",
    "frame['x-positions'] = (frame['Radius (m)'] * np.cos(frame['Azimuth (rad)']) * np.sin(frame['Zenith (rad)']))\n",
    "frame['y-positions'] = (frame['Radius (m)'] * np.sin(frame['Azimuth (rad)']) * np.sin(frame['Zenith (rad)']))\n",
    "frame['z-positions'] = (frame['Radius (m)'] * np.cos(frame['Zenith (rad)']))\n",
    "\n",
    "## The energy in eV will be 10 raised to the number in the file, multiplied by 1-y (y is inelasticity)\n",
    "frame['Nu Energies (eV)'] = np.power(10, (frame['Energy (log10) (eV)']))\n",
    "frame['Mu Energies (eV)'] = ((1-frame['Inelasticity']) * frame['Nu Energies (eV)']) #Energy of the produced lepton\n",
    "#Here the lepton is not a muon necesarily, hence the label 'Mu Energies (eV)' may be misleading\n",
    "\n",
    "## Get a frame with only coincident events\n",
    "coincident_frame = frame[frame['Coincident'] == 1] \n",
    "\n",
    "## And a frame for strictly events *detected* by ARA\n",
    "detected_frame = frame[frame['Detected'] == 1]\n",
    "\n",
    "\n",
    "## Now let's calculate the energy of the lepton when reaching IceCube (IC)\n",
    "\n",
    "# To do this correctly, I need to find exactly the distance traveled by the muon and apply the equation\n",
    "# I need the trajectory of the muon to find the time it takes to reach IceCube, then I can find the distance it travels in that time\n",
    "# I should allow events that occur inside the icecube volume to have their full energy (but pretty much will happen anyway)\n",
    "## a = sin(Theta)*cos(Phi)\n",
    "## b = sin(Theta)*sin(Phi)\n",
    "## c = cos(Theta)\n",
    "## a_0 = x-position\n",
    "## b_0 = y-position\n",
    "## c_0 = z-position\n",
    "## x_0 = IceCube[0]\n",
    "## y_0 = IceCube[1]\n",
    "## z_0 = IceCube[2]\n",
    "## t = (-(a*(a_0-x_0) + b*(b_0-y_0))+D**0.5)/(a**2+b**2)\n",
    "## D = (a**2+b**2)*R_IC**2 - (a*(b_0-y_0)+b*(a_0-x_0))**2\n",
    "## d = ((a*t)**2 + (b*t)**2 + (c*t)**2)**0.5\n",
    "\n",
    "## Trajectories\n",
    "coincident_frame['a'] = (np.sin(coincident_frame['Theta (rad)'])*np.cos(coincident_frame['Phi (rad)']))\n",
    "coincident_frame['b'] = (np.sin(coincident_frame['Theta (rad)'])*np.sin(coincident_frame['Phi (rad)']))\n",
    "coincident_frame['c'] = (np.cos(coincident_frame['Theta (rad)']))\n",
    "\n",
    "## Discriminant\n",
    "coincident_frame['D'] = ((coincident_frame['a']**2 + coincident_frame['b']**2)*IceCube_radius**2 - \n",
    "                         (coincident_frame['a']*(coincident_frame['y-position (m)']-IceCube[station-1][1])- ## I think this might need to be a minus sign!\n",
    "                          coincident_frame['b']*(coincident_frame['x-position (m)']-IceCube[station-1][0]))**2)\n",
    "\n",
    "## Interaction time (this is actually the same as the distance traveled, at least for a straight line)\n",
    "coincident_frame['t_1'] = (-(coincident_frame['a']*(coincident_frame['x-position (m)']-IceCube[station-1][0])+\n",
    "                            coincident_frame['b']*(coincident_frame['y-position (m)']-IceCube[station-1][1]))+\n",
    "                          np.sqrt(coincident_frame['D']))/(coincident_frame['a']**2+coincident_frame['b']**2)\n",
    "coincident_frame['t_2'] = (-(coincident_frame['a']*(coincident_frame['x-position (m)']-IceCube[station-1][0])+\n",
    "                            coincident_frame['b']*(coincident_frame['y-position (m)']-IceCube[station-1][1]))-\n",
    "                          np.sqrt(coincident_frame['D']))/(coincident_frame['a']**2+coincident_frame['b']**2)\n",
    "\n",
    "## Intersection coordinates\n",
    "coincident_frame['x-intersect_1'] = (coincident_frame['a'] * coincident_frame['t_1'] + coincident_frame['x-position (m)'])\n",
    "coincident_frame['y-intersect_1'] = (coincident_frame['b'] * coincident_frame['t_1'] + coincident_frame['y-position (m)'])\n",
    "coincident_frame['z-intersect_1'] = (coincident_frame['c'] * coincident_frame['t_1'] + coincident_frame['z-position (m)'])\n",
    "\n",
    "coincident_frame['x-intersect_2'] = (coincident_frame['a'] * coincident_frame['t_2'] + coincident_frame['x-position (m)'])\n",
    "coincident_frame['y-intersect_2'] = (coincident_frame['b'] * coincident_frame['t_2'] + coincident_frame['y-position (m)'])\n",
    "coincident_frame['z-intersect_2'] = (coincident_frame['c'] * coincident_frame['t_2'] + coincident_frame['z-position (m)'])\n",
    "\n",
    "## Distance traveled (same as the parametric time, at least for a straight line)\n",
    "coincident_frame['d_1'] = (np.sqrt((coincident_frame['a']*coincident_frame['t_1'])**2+\n",
    "                          (coincident_frame['b']*coincident_frame['t_1'])**2+\n",
    "                          (coincident_frame['c']*coincident_frame['t_1'])**2))\n",
    "coincident_frame['d_2'] = (np.sqrt((coincident_frame['a']*coincident_frame['t_2'])**2+\n",
    "                          (coincident_frame['b']*coincident_frame['t_2'])**2+\n",
    "                          (coincident_frame['c']*coincident_frame['t_2'])**2))\n",
    "\n",
    "## Check if it started inside and set the distance based on if it needs to travel to reach icecube or not\n",
    "coincident_frame['Inside'] = (np.where((coincident_frame['t_1']/coincident_frame['t_2'] < 0) & (coincident_frame['z-position (m)'].between(-2450, -1450)), 1, 0))\n",
    "coincident_frame['preliminary d'] = (np.where(coincident_frame['d_1'] <= coincident_frame['d_2'], coincident_frame['d_1'], coincident_frame['d_2']))\n",
    "coincident_frame['d'] = (np.where(coincident_frame['Inside'] == 1, 0, coincident_frame['preliminary d']))\n",
    "\n",
    "## Check if the event lies in the cylinder\n",
    "coincident_frame['In IC'] = (np.where((np.sqrt((coincident_frame['x-position (m)']-IceCube[station-1][0])**2 + (coincident_frame['y-position (m)']-IceCube[station-1][1])**2) < IceCube_radius) &\n",
    "                                     ((coincident_frame['z-position (m)']).between(-2450, -1450)) , 1, 0))\n",
    "\n",
    "#Correct coincident_frame to only have electron neutrinos inside IC\n",
    "coincident_frame = coincident_frame[(((coincident_frame['In IC'] == 1) & (coincident_frame['flavor'] == 1)) | (coincident_frame['flavor'] == 2) | (coincident_frame['flavor'] == 3)) ]\n",
    "\n",
    "#Now calculate the lepton energies when they reach IC\n",
    "coincident_frame['IC Mu Energies (eV)'] = (coincident_frame['Mu Energies (eV)'] * np.exp(-10**-5 * coincident_frame['d']*100)) # convert d from meters to cm\n",
    "coincident_frame['weighted energies'] = (coincident_frame['weight'] * coincident_frame['Nu Energies (eV)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f1757103",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Add possible Tau decay to the frame\n",
    "coincident_frame['Tau decay'] = ''\n",
    "# Again, the label 'Tau Decay' may be misleading because not all leptons may be taus\n",
    "\n",
    "## Calculate distance from the interaction point to its walls and keep the shortest (the first interaction with the volume)\n",
    "\n",
    "coincident_frame['distance-to-IC_1'] = np.sqrt((coincident_frame['x-positions'] - coincident_frame['x-intersect_1'])**2 + \n",
    "                                        (coincident_frame['y-positions'] - coincident_frame['y-intersect_1'])**2)\n",
    "coincident_frame['distance-to-IC_2'] = np.sqrt((coincident_frame['x-positions'] - coincident_frame['x-intersect_2'])**2 + \n",

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a950b9af",
   "metadata": {},
   "source": [
    "**This script is simply for me to calculate the location of IceCube relative to the origin of any ARA station**\n",
    "\n",
    "The relevant documentation to understand the definitions after the imports can be found in https://elog.phys.hawaii.edu/elog/ARA/130712_170712/doc.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b926e2e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "901b442b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Definitions of translations in surveyor's coordinates:\n",
    "\n",
    "t_IC_to_ARAg = np.array([-24100, 1700, 6400])\n",
    "t_ARAg_to_A1 = np.array([16401.71, -2835.37, -25.67])\n",
    "t_ARAg_to_A2 = np.array([13126.7, -8519.62, -18.72])\n",
    "t_ARAg_to_A3 = np.array([9848.35, -2835.19, -12.7])\n",
    "\n",
    "#Definitions of rotations from surveyor's axes to the ARA Station's coordinate systems\n",
    "\n",
    "R1 = np.array([[-0.598647, 0.801013, -0.000332979], [-0.801013, -0.598647, -0.000401329], \\\n",
    "               [-0.000520806, 0.0000264661, 1]])\n",
    "R2 = np.array([[-0.598647, 0.801013, -0.000970507], [-0.801007, -0.598646,-0.00316072 ], \\\n",
    "               [-0.00311277, -0.00111477, 0.999995]])\n",
    "R3 = np.array([[-0.598646, 0.801011, -0.00198193],[-0.801008, -0.598649,-0.00247504], \\\n",
    "               [-0.00316902, 0.000105871, 0.999995]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab2d3206",
   "metadata": {},
   "source": [
    "**Using these definitions, I should be able to calculate the location of IceCube relative to each ARA station by:**\n",
    "\n",
    "$$\n",
    "\\vec{r}_{A 1}^{I C}=-R_1\\left(\\vec{t}_{I C}^{A R A}+\\vec{t}_{A R A}^{A 1}\\right)\n",
    "$$\n",
    "\n",
    "We have a write-up of how to get this. Contact salcedogomez.1@osu.edu if you need that.\n",
    "\n",
    "Alex had done this already, he got that \n",
    "\n",
    "$$\n",
    "\\vec{r}_{A 1}^{I C}=-3696.99^{\\prime} \\hat{x}-6843.56^{\\prime} \\hat{y}-6378.31^{\\prime} \\hat{z}\n",
    "$$\n",
    "\n",
    "Let me verify that I get the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "912163d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IC coordinates relative to A1 (in):  [-3696.98956579 -6843.55800868 -6378.30926681]\n",
      "IC coordinates relative to A1 (m):  [-1127.13096518 -2086.4506124  -1944.60648378]\n",
      "Distance of IC from A1 (m):  3066.788996234438\n"
     ]
    }
   ],
   "source": [
    "IC_A1 = -R1 @ np.add(t_ARAg_to_A1, t_IC_to_ARAg).T\n",
    "print(\"IC coordinates relative to A1 (in): \", IC_A1)\n",
    "print(\"IC coordinates relative to A1 (m): \", IC_A1/3.28)\n",
    "print(\"Distance of IC from A1 (m): \", np.sqrt((IC_A1[0]/3.28)**2 + (IC_A1[1]/3.28)**2 + (IC_A1[2]/3.28)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c9f252",
   "metadata": {},
   "source": [
    "Looks good!\n",
    "\n",
    "Now, I just get the other ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8afa27c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IC coordinates relative to A2 (in):  [ -1100.33577313 -12852.0589083   -6423.00776043]\n",
      "IC coordinates relative to A2 (m):  [ -335.46822352 -3918.31064277 -1958.2340733 ]\n",
      "Distance of IC from A2 (m):  4393.219537890439\n"
     ]
    }
   ],
   "source": [
    "IC_A2 = -R2 @ np.add(t_ARAg_to_A2, t_IC_to_ARAg).T\n",
    "print(\"IC coordinates relative to A2 (in): \", IC_A2)\n",
    "print(\"IC coordinates relative to A2 (m): \", IC_A2/3.28)\n",
    "print(\"Distance of IC from A2 (m): \", np.sqrt((IC_A2[0]/3.28)**2 + (IC_A2[1]/3.28)**2 + (IC_A2[2]/3.28)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9959d0a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IC coordinates relative to A3 (in):  [ -7609.73440732 -12079.45719852  -6432.31164368]\n",
      "IC coordinates relative to A3 (m):  [-2320.04097784 -3682.76134101 -1961.07062307]\n",
      "Distance of IC from A3 (m):  4774.00452685144\n"
     ]
    }
   ],
   "source": [
    "IC_A3 = -R3 @ np.add(t_ARAg_to_A3, t_IC_to_ARAg).T\n",
    "print(\"IC coordinates relative to A3 (in): \", IC_A3)\n",
    "print(\"IC coordinates relative to A3 (m): \", IC_A3/3.28)\n",
    "print(\"Distance of IC from A3 (m): \", np.sqrt((IC_A3[0]/3.28)**2 + (IC_A3[1]/3.28)**2 + (IC_A3[2]/3.28)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "093dff67",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}