Updates and Results Talks and Posters Advice Ideas Important Figures Write-Ups Outreach How-To Funding Opportunities GENETIS
  Important Plots, Tables, and Measurements  ELOG logo
Entry  Sun Apr 23 14:54:50 2017, Hannah Hasan, Other, Simulation, Plotting ShelfMC Parameter Space, Other Veff_A_I.pdfVeff_R_S.pdf
    Reply  Fri Oct 6 15:15:53 2017, Hannah Hasan, Other, Simulation, Plotting ShelfMC Parameter Space, Other 6x
Message ID: 27     Entry time: Fri Oct 6 15:15:53 2017     In reply to: 16
Author: Hannah Hasan 
Type: Other 
Category: Simulation 
Subject: Plotting ShelfMC Parameter Space 
Project: Other 

Attached are instructions and scripts for carrying out a parameter space scan with ShelfMC, the simulation package for the ARIANNA detector.

Because some of the plotted outputs looked like colored stripes and did not offer any insight into how effective volume changed with some variables, I made some changes to the simulation and plotting scripts so that different maximum, minimum, and increment values can be chosen for each variable. Now rather than having fixed, hard-coded values for all variables, the parameter space scan and plotting is more flexible for use with variable inputs.

Quote:

I am trying to write a script that will plot a 2d histogram of effective volume versus two of ShelfMC's parameters.

The script prompts the user for which two parameters (out of five that we vary in our parameter space scan) to plot along the x- and y-axes, as well as what values to hold the other 3 parameters constant at. It then collects the necessary root files from simulation results, generates a plotting script, and runs the plotting script to produce a plot in pdf form.

After many struggles I have the script written to the point where it functions, but the plots don't look right. Some plots look like they could be actual data (like Veff_A_I), and others just look flat-out wrong (like Veff_R_S).

I have yet to pin down the cause of this, but hopefully will be able to sometime in the near future.

 

Attachment 1: ParameterSpaceScan_instructions.txt  3 kB  | Show | Hide all | Show all
Attachment 2: setup.sh  8 kB  | Hide | Hide all | Show all
#!/bin/bash
#PBS -l walltime=04:00:00
#PBS -l nodes=1:ppn=1,mem=4000mb
#PBS -N hannah_SetupJob
#PBS -j oe
#PBS -A PCON0003
#Jude Rajasekera 3/20/17
#Modified by Hannah Hasan on 09/17/2017

#directories
WorkDir=$TMPDIR   
tmpShelfmc=$HOME/ShelfMC/git_shelfmc #set your ShelfMC directory here


AttenMin=600                      # MINIMUM attenuation length for the simulation      *****ALL UNITS ARE IN METERS*****
AttenMax=1000                     # MAXIMUM attenuation length
AttenInc=400                      # INCREMENT for attenuation length

RadiusMin=3                       # MINIMUM radius
RadiusMax=31                      # MAXIMUM radius
RadiusInc=7                       # INCREMENT

IceMin=500                        # etc...
IceMax=2900
IceInc=400

FirnMin=60
FirnMax=120
FirnInc=60

StDepthMin=0
StDepthMax=200
StDepthInc=50

#controlled variables for run
runName='e18_ParamSpaceScan' #name of run
NNU=10000000 #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=20 #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="0.85    #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="$AR     #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/cond0092/ShelfMC/git_shelfmc/GainFiles/LP_gain_manual.txt     #GAINFILENAME"
input40="$SD     #STATION_DEPTH"
#######################################################################################################

cd $TMPDIR   
mkdir $runName
cd $runName

initSeed=$seed
counter=0

for ((L=$AttenMin;L<=$AttenMax;L+=$AttenInc)) #Attenuation length
do
    mkdir Atten_Up$L
    cd Atten_Up$L
    for ((AR=$RadiusMin;AR<=$RadiusMax;AR+=$RadiusInc)) #Station radius (measured from center of radius to antenna
    do
	mkdir AntennaRadius$AR
        cd AntennaRadius$AR
	for ((T=$IceMin;T<=$IceMax;T+=$IceInc)) #Ice thickness
        do
            mkdir IceThick$T
            cd IceThick$T
            for ((FT=$FirnMin;FT<=$FirnMax;FT+=$FirnInc)) #Firn thickness
            do
		mkdir FirnThick$FT
                cd FirnThick$FT
		for ((SD=$StDepthMin;SD<=$StDepthMax;SD+=$StDepthInc)) #Station depth
                do
                    mkdir StationDepth$SD
                    cd StationDepth$SD
                    #####Do file operations###########################################
                    counter=$((counter+1))
                    echo "Counter = $counter ; L = $L ; AR = $AR ; 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)"
		    input35="$AR     #ST4_R radius in meters between center of station and antenna"
                    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:10:00' >> run_shelfmc_multithread.sh #change walltime as necessary
		    echo '#PBS -N hannah_'$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'/'AntennaRadius$AR'/'IceThick$T'/'FirnThick$FT'/'StationDepth$SD'/Setup'$i' _'$i'$runName &' >> run_shelfmc_multithread.sh
		    done
		   # echo './shelfmc_stripped.exe $runName/'Atten_Up$L'/'AntennaGap$ATGap'/'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'/'AntennaRadius$AR'/'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
		    echo 'rm -rf Setup*' >> 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
Attachment 3: scheduler.sh  2 kB  | Show | Hide all | Show all
Attachment 4: colorplot_instructions.txt  8 kB  | Show | Hide all | Show all
Attachment 5: colorplot.sh  14 kB  | Hide | Hide all | Show all
#!/bin/bash
#04/18/2017
# script will collect root files and plot as dictated by user
#Hannah Hasan

shelfmc=$HOME/ShelfMC/git_shelfmc ##set to your ShelfMC directory - should have the parameter space scan directory and a directory called "outputs". outputs may be empty

AttenMin=500
AttenMax=1000
AttenInc=100

RadiusMin=3
RadiusMax=31
RadiusInc=7

IceMin=500
IceMax=2900
IceInc=400

FirnMin=60
FirnMax=140
FirnInc=20

StDepthMin=0
StDepthMax=200
StDepthInc=50

line1='{'
line2='  TFile *f = new TFile("'"Veff_${xvar}_${yvar}.root"'");'   #set later
line3='  TTree *SimParam = (TTree*)f->Get("sim_parameters");'
line4='  TCanvas *c = new TCanvas("c", "V_{eff}", 2000, 1000);'
line5='  TPad* p = new TPad("p","p",0.0,0.0,1.0,1.0);'
line6='  p->Draw();'
line7='  p->SetLeftMargin(0.15);'
line8='  p->SetRightMargin(0.15);'
line9='  p->cd();'
line10='  TH2F *h = new TH2F("''h''", "''V_{eff} vs '"$var1"' vs '"$var2"'"'",$xbins,$xmin,$xmax"', '"$ybins,$ymin,$ymax);"   set later
line11='  double x;'
line12='  double y;'
line13='  double z;'
line14='  double veff=0;'
line15='  double energy;'
line16="  double avgs[$SIZE];"   #set later; SIZE varies depending on which variables we are plotting
line17='  SimParam->SetBranchAddress("'"$xbranch"'", &x );'   #set later
line18='  SimParam->SetBranchAddress("'"$ybranch"'", &y );'   #set later
line19='  SimParam->SetBranchAddress("Veff", &z);'
line20='  int v=0;'
line21="  for(int k=1; k<=$SIZE; k++){"   #set later for SIZE
line22='    veff = 0;'
line23='    for(v; v<(k*20.0); v++){'       #change 20 to number of processors ShelfMC runs were split across
line24='      SimParam->GetEntry(v);'
line25='      veff += z;'
line26='    }'
line27='    avgs[(k-1)] = (veff/20.0);'      #change 20 to number of processors
line28='  }'
line29="  for(int i=0; i<$SIZE; i++){"   #set later for SIZE
line30='    SimParam->GetEntry((i*20));'       #change 20 to number of processors
line31='    h->Fill(x,y,avgs[i]);'
line32='  }'
line33='  gStyle->SetPalette(1,0);'
line34='  h->GetXaxis()->SetTitle("'"$var1"' (m)");'   #set later
line35='  h->GetYaxis()->SetTitle("'"$var2"' (m)");'   #set later
line36='  h->GetZaxis()->SetTitle("V_{eff} (km^{3} str)");'
line37='  h->GetXaxis()->SetTitleOffset(1.2);'
line38='  h->GetYaxis()->SetTitleOffset(1.0);'
line39='  h->GetZaxis()->SetTitleOffset(1.3);'
line40='  gStyle->SetOptStat(0);'
line41='  h->Draw("COLZ");'
line42='  double const1;'
line43='  double const2;'
line44='  double const3;'
line45='  SimParam->SetBranchAddress("'"$const1br"'", &const1);'   #set later
line46='  SimParam->SetBranchAddress("'"$const2br"'", &const2);'   #set later
line47='  SimParam->SetBranchAddress("'"$const3br"'", &const3);'   #set later
line48='  SimParam->SetBranchAddress("exponent", &energy);'
line49='  SimParam->GetEntry(0);'
line50='  char const_1[30];'
line51='  char const_2[30];'
line52='  char const_3[30];'
line53='  char enrgy[30];'
line54='  sprintf(enrgy, "Energy (exponent): %d", energy);'
line55='  sprintf(const_1, "'"$const1br"' %d m", const1);'   
line56='  sprintf(const_2, "'"$const2br"' %d m", const2);'   
line57='  sprintf(const_3, "'"$const3br"' %d m", const3);'    #set later
line58='  TPaveText *t = new TPaveText(0,0,0.07,0.3,"brNDC");'
line59='  t->SetLabel("Simulation Parameters:");'
line60='  t->AddText(enrgy);'
line61='  t->AddText(const_1);'
line62='  t->AddText(const_2);'
line63='  t->AddText(const_3);'
line64='  t->Draw();'
line65='  c->Modified();'
line66='  c->Update();'
line67='  c->SaveAs("'"Veff_${xvar}_${yvar}.pdf"'");'   #set later
line68='}'


echo
echo
echo 'This script will allow you to choose two variables to plot effective volume against.'
echo
echo 'Enter the name of the run you wish to plot data for.'
echo
read -p 'Name of run: ' runName
echo
echo 'Thank you!'
echo
echo 'Please make your selections from the following list:'
echo '(A) - Attenuation Length'
echo '(R) - Antenna Radius'
echo '(I) - Ice Thickness'
echo '(F) - Firn Thickness'
echo '(S) - Station Depth'
echo
echo '** please use capital letters in your selection **'
echo
read -p 'Select the variable you wish to plot along the X-axis: ' xvar
read -p 'Now select the variable you wish to plot along the Y-axis: ' yvar
echo


cd $shelfmc/outputs
mkdir ${runName}_${xvar}_${yvar}

cd $shelfmc/$runName

attenprompt='Select what value at which to hold constant Attenuation Length'
atteninc=' > Choose only values between 500-1000, at increments of 100: '

radiusprompt='Select what value at which to hold constant Antenna Radius'
radiusinc=' > Choose only values between 3-31, at increments of 7: '

iceprompt='Select what value at which to hold constant Ice Thickness'
iceinc=' > Choose only values between 500-2900, at increments of 400: '

firnprompt='Select what value at which to hold constant Firn Thickness'
firninc=' > Choose only values between 60-140, at increments of 20: '

stprompt='Select what value at which to hold constant Station Depth'
stinc=' > Choose only values between 0-200, at increments of 50: '

if ([ "$xvar" = "A" ] && [ "$yvar" = "R" ]) || ([ "$yvar" = "A" ] && [ "$xvar" = "R" ]); then
     echo $iceprompt
     echo $iceinc
     read const1
    const1br='icethick'
     echo $firnprompt
     echo $firninc
     read const2
    const2br='firndepth'
     echo $stprompt
     echo $stinc
     read const3
    const3br='station_depth'
    echo
    echo 'Adding files...'
    echo
    for ((L=$AttenMin;L<=$AttenMax;L+=$AttenInc))
    do
	cd Atten_Up$L
	for ((AR=$RadiusMin;AR<=$RadiusMax;AR+=$RadiusInc))
	do
	    cd AntennaRadius$AR/IceThick$const1/FirnThick$const2/StationDepth$const3
            cp *.root $shelfmc/outputs/${runName}_${xvar}_${yvar}/Result_${L}_${AR}.root
	    cd ../../../../
	done
	cd ..
    done

elif ([ "$xvar" = "A" ] && [ "$yvar" = "I" ]) || ([ "$yvar" = "A" ] && [ "$xvar" = "I" ]); then
     echo $radiusprompt
     echo $radiusinc
     read const1
    const1br='st4_r'
     echo $firnprompt
     echo $firninc
     read const2
    const2br='firndepth'
     echo $stprompt
     echo $stinc
     read const3
    const3br='station_depth'
    echo
    for ((L=$AttenMin;L<=$AttenMax;L+=$AttenInc))
    do
        cd Atten_Up$L/AntennaRadius$const1
        for ((Ice=$IceMin;Ice<=$IceMax;Ice+=$IceInc))
        do
            cd IceThick$Ice/FirnThick$const2/StationDepth$const3
            cp *.root $shelfmc/outputs/${runName}_${xvar}_${yvar}/Result_${L}_${Ice}.root
            cd ../../../
        done
        cd ../../
    done

elif ([ "$xvar" = "A" ] && [ "$yvar" = "F" ]) || ([ "$yvar" = "A" ] && [ "$xvar" = "F" ]); then
     echo $radiusprompt
     echo $radiusinc
     read const1
    const1br='st4_r'
     echo $iceprompt
     echo $iceinc
     read const2
    const2br='icethick'
     echo $stprompt
     echo $stinc
     read const3
    const3br='station_depth'
    echo
    for ((L=$AttenMin;L<=$AttenMax;L+=$AttenInc))
    do
	cd Atten_Up$L/AntennaRadius$const1/IceThick$const2
	for ((FT=$FirnMin;FT<=$FirnMax;FT+=$FirnInc))
	do
	    cd FirnThick$FT/StationDepth$const3
	    cp *.root $shelfmc/outputs/${runName}_${xvar}_${yvar}/Result_${L}_${FT}.root
	    cd ../..
	done
	cd ../../..
    done

elif ([ "$xvar" = "A" ] && [ "$yvar" = "S" ]) || ([ "$yvar" = "A" ] && [ "$xvar" = "S" ]); then
     echo $radiusprompt
     echo $radiusinc
     read const1
    const1br='st4_r'
     echo $iceprompt
     echo $iceinc
     read const2
    const2br='icethick'
     echo $firnprompt
     echo $firninc
     read const3
    const3br='firndepth'
    echo
    for ((L=$AttenMin;L<=$AttenMax;L+=$AttenInc))
    do
	cd Atten_Up$L/AntennaRadius$const1/IceThick$const2/FirnThick$const3
	for ((SD=$StDepthMin;SD<=$StDepthMax;SD+=$StDepthInc))
	do
	    cd StationDepth$SD
	    cp *.root $shelfmc/outputs/${runName}_${xvar}_${yvar}/Result_${L}_${SD}.root
	    cd ../
	done
	cd ../../../..
    done

elif ([ "$xvar" = "R" ] && [ "$yvar" = "I" ]) || ([ "$yvar" = "R" ] && [ "$xvar" = "I" ]); then
     echo $attenprompt
     echo $atteninc
     read const1
    const1br='atten_up'
     echo $firnprompt
     echo $firninc
     read const2
    const2br='firndepth'
     echo $stprompt
     echo $stinc
     read const3
    const3br='station_depth'
    echo
    for ((AR=$RadiusMin;AR<=$RadiusMax;AR+=$RadiusInc))
    do
	cd Atten_Up$const1/AntennaRadius$AR
	for ((Ice=$IceMin;Ice<=$IceMax;Ice+=$IceInc))
	do
	    cd IceThick$Ice/FirnThick$const2/StationDepth$const3
            cp *.root $shelfmc/outputs/${runName}_${xvar}_${yvar}/Result_${AR}_${Ice}.root
	    cd ../../..
	done
	cd ../..
    done

elif ([ "$xvar" = "R" ] && [ "$yvar" = "F" ]) || ([ "$yvar" = "R" ] && [ "$xvar" = "F" ]); then
     echo $attenprompt
     echo $atteninc
     read const1
    const1br='atten_up'
     echo $iceprompt
     echo $iceinc
     read const2
    const2br='icethick'
     echo $stprompt
     echo $stinc
     read const3
    const3br='station_depth'
    echo
    for ((AR=$RadiusMin;AR<=$RadiusMax;AR+=$RadiusInc))
    do
	cd Atten_Up$const1/AntennaRadius$AR/IceThick$const2
	for ((FT=$FirnMin;FT<=$FirnMax;FT+=$FirnInc))
	do
	    cd FirnThick$FT/StationDepth$const3
            cp *.root $shelfmc/outputs/${runName}_${xvar}_${yvar}/Result_${AR}_${FT}.root
	    cd ../..
	done
	cd ../../..
    done

elif ([ "$xvar" = "R" ] && [ "$yvar" = "S" ]) || ([ "$yvar" = "R" ] && [ "$xvar" = "S" ]); then
... 215 more lines ...
Attachment 6: colorplot_loop.sh  32 kB  | Show | Hide all | Show all
ELOG V3.1.5-fc6679b