diff --git a/montecarlo/fax_waveform/CreateFakeCSV_CorrelatedS1S2.py b/montecarlo/fax_waveform/CreateFakeCSV_CorrelatedS1S2.py new file mode 100644 index 0000000..2a476f9 --- /dev/null +++ b/montecarlo/fax_waveform/CreateFakeCSV_CorrelatedS1S2.py @@ -0,0 +1,122 @@ +################################# +## Sub-code used in WF simulation +## It creates a csv file for the input of fax +## using the 2d-pdf from an input histogram +## by Qing Lin +## @ 2016-09-12 +## +## HARDCODE WARNING: The FV dimensions below need to be modified +## according to the detector you wish to simulate +################################# +import sys +import numpy as np +import scipy as sp + +if len(sys.argv)<2: + print("========= Syntax ==========") + print("python CreateFakeCSV.py ..... ") + print("") + print("") + print("") + print("") + print("") + print("") + print("") + exit() + +Detector = sys.argv[1] +Band2DPDFFilename = sys.argv[2] +NominalG1Value = float(sys.argv[3]) +NominalG2Value = float(sys.argv[4]) +NumEvents = int(sys.argv[5]) +DefaultType = sys.argv[6] +OutputFilename = sys.argv[7] + +#################################### +## Some nuisance parameters (HARDCODE WARNING): +#################################### +MaxDriftTime = 650. # us + + +#################################### +## Some functions (HARDCODE WARNING): +#################################### + +# Current FV cut for Xe1T +scalecmtomm=1 +def radius2_cut(zpos): + return 1400*scalecmtomm**2+(zpos+100*scalecmtomm)*(2250-1900)*scalecmtomm/100 + +def IfPassFV(x,y,z): + + if Detector == "XENON100": + # check if the x,y,z passing X48kg0 + I = np.power( (z+15.)/14.6, 4.) + I += np.power( (x**2+y**2)/20000., 4.) + if I<1: + return True + elif Detector == "XENON1T": # NEED TO UPDATE THIS + Zlower, Zupper = -90*scalecmtomm, -15*scalecmtomm + Zcut = ((z>=Zlower) & (z<=Zupper)) + R2upper=radius2_cut(z) + Rcut = (x**2+y**2") + print("") + print("") + exit() + +InputROOTFilename = sys.argv[1] +HistogramName = sys.argv[2] +OutputPickleFilename = sys.argv[3] + + +#################################### +# Open the root file and load the TH2D +#################################### +pfile = TFile(InputROOTFilename) +hist = pfile.Get(HistogramName) + +#################################### +# translate into a dictionary +#################################### +data = {} +data['s1nbins'] = int( hist.GetXaxis().GetNbins() ) +data['s1lower'] = float( hist.GetXaxis().GetXmin() ) +data['s1upper'] = float( hist.GetXaxis().GetXmax() ) +data['lognbins'] = int( hist.GetYaxis().GetNbins() ) +data['loglower'] = float( hist.GetYaxis().GetXmin() ) +data['logupper'] = float( hist.GetYaxis().GetXmax() ) +data['map'] = [] + +for i in range(data['s1nbins']): + TmpList = [] + for j in range(data['lognbins']): + cont = float(hist.GetBinContent(i+1, j+1) ) + TmpList.append(cont) + data['map'].append( list(TmpList) ) + +#################################### +## output to pickle +#################################### + +pickle.dump( data, open(OutputPickleFilename, 'wb') ) diff --git a/montecarlo/fax_waveform/MidwayBatch_AreaCorrelatedS1S2.py b/montecarlo/fax_waveform/MidwayBatch_AreaCorrelatedS1S2.py new file mode 100644 index 0000000..144ddb2 --- /dev/null +++ b/montecarlo/fax_waveform/MidwayBatch_AreaCorrelatedS1S2.py @@ -0,0 +1,101 @@ +#################################### +## batch code for WF simulation +#################################### +import sys, array, os, getpass +from subprocess import call +import subprocess as subp +import time +import math as math +from subprocess import Popen, PIPE + +if len(sys.argv)<2: + print("========= Syntax ========") + print("python MidwayBatch_AreaCorrelatedS1S2.py ....") + print("") + print("") + print("") + print("") + print("") + print("<2D band for S1-S2 area correlation (abs.)>") + print("") + print("") + print("") + exit() + +OutputGeneralPath = sys.argv[1] +NumJobs = int(sys.argv[2]) +NumEvents = int(sys.argv[3]) +PMTAfterpulseFlag = int(sys.argv[4]) +S2AfterpulseFlag = int(sys.argv[5]) +Input2DBandFile = sys.argv[6] +Nomial_g1 = float(sys.argv[7]) +Nomial_g2 = float(sys.argv[8]) +IfUsePublicNodes = int(sys.argv[9]) + +MaxNumJob = 64 +if not IfUsePublicNodes: + MaxNumJob=200 + + +##### Start batching ######### +CurrentPath = os.getcwd() +print (CurrentPath) +CurrentUser = getpass.getuser() +for i in range(NumJobs): + + RunString = "%06d" % i + + # create folder + OutputPath = OutputGeneralPath + "/" + RunString + if os.path.exists(OutputPath): + subp.call("rm -r "+OutputPath, shell=True) + subp.call("mkdir -p "+OutputPath, shell=True) + + # define filenames + SubmitFile = OutputPath+"/submit_"+ RunString + ".sh" + SubmitOutputFilename = OutputPath+"/submit_"+ RunString + ".log" + SubmitErrorFilename = OutputPath+"/submit_"+ RunString + ".log" + + # create the basic submit + subp.call("echo '#!/bin/bash\n' >> "+SubmitFile, shell=True) + subp.call("echo '#SBATCH --output="+SubmitOutputFilename+"' >> "+SubmitFile, shell=True) + subp.call("echo '#SBATCH --error="+SubmitErrorFilename+"' >> "+SubmitFile, shell=True) + subp.call("echo '#SBATCH --time=03:59:00' >> "+SubmitFile, shell=True) + subp.call("echo '#SBATCH --account=pi-lgrandi' >> "+SubmitFile, shell=True) + if IfUsePublicNodes==0: + subp.call("echo '#SBATCH --qos=xenon1t' >> "+SubmitFile, shell=True) + subp.call("echo '#SBATCH --partition=xenon1t\n' >> "+SubmitFile, shell=True) + elif IfUsePublicNodes==2: + subp.call("echo '#SBATCH --qos=xenon1t-kicp' >> "+SubmitFile, shell=True) + subp.call("echo '#SBATCH --partition=kicp\n' >> "+SubmitFile, shell=True) + + Command = CurrentPath+"/./run_fax_AreaCorrelatedS1S2.sh "+Input2DBandFile+" "+str(Nomial_g1)+" "+str(Nomial_g2)+" "+str(PMTAfterpulseFlag)+" "+str(S2AfterpulseFlag)+" "+str(NumEvents)+" "+OutputGeneralPath+" "+RunString + subp.call("echo '"+Command+"\n' >> "+SubmitFile, shell=True) + + SubmitPath = OutputPath + + #submit + IfSubmitted=0 + while IfSubmitted==0: + Partition = "sandyb" # public + if not IfUsePublicNodes: + Partition = "xenon1t" + elif IfUsePublicNodes==2: + Partition = "kicp" + p1 = Popen(["squeue","--partition="+Partition, "--user="+CurrentUser], stdout=PIPE) + p2 = Popen(["wc", "-l"], stdin=p1.stdout, stdout=PIPE) + p1.stdout.close() # Allow p1 to receive a SIGPIPE if p2 exits. + output = p2.communicate()[0] + Status=subp.call("squeue --partition="+Partition+" --user="+CurrentUser +" | wc -l", shell=True) + Output=int(output) + #print(Status) + + print("Current job running number "+str(Output)) + + if Status==0 and Output=len(self.m_pLogInterpolators): + s1_bin_upper = s1_bin + s1lower = float(s1_bin)*self.m_pS1Step + self.m_pS1Lower + s1upper = float(s1_bin_upper)*self.m_pS1Step + self.m_pS1Lower + log_cdf = np.random.uniform(0, 1) + loglower = self.m_pLogInterpolators[s1_bin](log_cdf) + logupper = self.m_pLogInterpolators[s1_bin_upper](log_cdf) + log = loglower + if s1lower!=s1upper: + log = (s1 - s1lower) / (s1upper - s1lower) * (logupper - loglower) + loglower + s2 = np.power(10., log)*s1 + return (s1/self.m_pG1Value, s2/self.m_pG2Value) + + + diff --git a/montecarlo/fax_waveform/README.md b/montecarlo/fax_waveform/README.md index dd0c4ce..8ca4dcb 100644 --- a/montecarlo/fax_waveform/README.md +++ b/montecarlo/fax_waveform/README.md @@ -28,4 +28,4 @@ You can use this code to produce fax-only data with Qing's framework in one term - `nodetype` : 0 for xenon1t, 1 for public, 2 for kicp - `use_array_truth` : set to 1 to put truth information in arrays - `save_ap_truth` : if `use_array_truth=1`, set this to 1 to save afterpulse truth info -- `minitree_type` : 0 for basics, 1 for S1S2Properties minitrees, 2 for PeakEfficiency minitrees +- `minitree_type` : 0 for basics, 1 for S1S2Properties minitrees, 2 for PeakEfficiency minitrees \ No newline at end of file diff --git a/montecarlo/fax_waveform/ReduceDataNormal.py b/montecarlo/fax_waveform/ReduceDataNormal.py index 47f87a9..c98b382 100644 --- a/montecarlo/fax_waveform/ReduceDataNormal.py +++ b/montecarlo/fax_waveform/ReduceDataNormal.py @@ -52,7 +52,7 @@ # hax.init(experiment='XENON1T') # my own builder -from collections import defaultdict +#~ from collections import defaultdict class S1S2Properties(hax.minitrees.TreeMaker): """Computing properties of the S1 diff --git a/montecarlo/fax_waveform/TruthSorting.py b/montecarlo/fax_waveform/TruthSorting.py index 40a9cb6..8cb76fb 100644 --- a/montecarlo/fax_waveform/TruthSorting.py +++ b/montecarlo/fax_waveform/TruthSorting.py @@ -96,6 +96,11 @@ tag = 0 if truth_data['peak_type'][iteration_id] == 's2': tag = 1 + elif ifcounteds1==0: + tag=0 + ifcounteds1=1 + else: + tag=2 if tag==0: #print("Iterator: "+str(iteration_id)+" -> S1") s1_time_truth = truth_data['t_mean_photons'][iteration_id] diff --git a/montecarlo/fax_waveform/TruthSorting_arrays.py b/montecarlo/fax_waveform/TruthSorting_arrays.py index da7610e..53abd5f 100644 --- a/montecarlo/fax_waveform/TruthSorting_arrays.py +++ b/montecarlo/fax_waveform/TruthSorting_arrays.py @@ -76,7 +76,7 @@ result[field] = [] ifcounteds1 = 0 - + while truth_data['event'][iteration_id]==event_id: tag = 2 # 0 for s1, 1 for s2, 2 for photoionization if truth_data['peak_type'][iteration_id] == 's1': diff --git a/montecarlo/fax_waveform/logs/desc.txt b/montecarlo/fax_waveform/logs/desc.txt new file mode 100644 index 0000000..535f0df --- /dev/null +++ b/montecarlo/fax_waveform/logs/desc.txt @@ -0,0 +1 @@ +placeholder, this is where log files will be diff --git a/montecarlo/fax_waveform/run_fax.sh b/montecarlo/fax_waveform/run_fax.sh index 18d6420..b3b9a61 100755 --- a/montecarlo/fax_waveform/run_fax.sh +++ b/montecarlo/fax_waveform/run_fax.sh @@ -30,7 +30,8 @@ PMTAfterpulseEnableFlag=$5 S2AfterpulseEnableFlag=$6 # Select fax+pax version -PAXVERSION=head +#~ PAXVERSION=head +PAXVERSION=v6.5.1 # temporary change # Specify number of events NumEvents=$7 @@ -78,7 +79,11 @@ FAX_FILENAME=${FILENAME}_truth # fax truth info PKL_FILENAME=${FILENAME}_truth.pkl # converted fax truth info RAW_FILENAME=${FILENAME}_raw # fax simulated raw data PAX_FILENAME=${FILENAME}_pax # pax processed data +PAX_FILENAME_WOPATH=${FILEROOT}_pax # pax processed data without path +MERGEDTRUTH_FILENAME=${FILENAME}_merged_truth +MERGED_FILENAME=${FILENAME}_merged HAX_FILENAME=${FILENAME}_hax # hax reduced data +MINITREE_FILENAME=${FILENAME}_pax_S1S2Properties.root CustomIniFilename=${RELEASEDIR}/NoS2Afterpulses.ini NoPMTAfterpulseIniFilename=${RELEASEDIR}/NoPMTAfterpulses.ini echo ${CustomIniFilename} @@ -117,12 +122,21 @@ python ${RELEASEDIR}/ConvertFaxTruthToPickle.py ${FAX_FILENAME} ${PKL_FILENAME} # pax stage (time paxer --ignore_rundb --input ${RAW_FILENAME} --config ${Detector} --output ${PAX_FILENAME};) &> ${PAX_FILENAME}.log -# hax stage -HAXPYTHON="import hax; " -HAXPYTHON+="hax.init(main_data_paths=['${OUTDIR}'], minitree_paths=['${OUTDIR}'], pax_version_policy = 'loose'); " -HAXPYTHON+="hax.minitrees.load('${PAX_FILENAME##*/}', ['Basics', 'Fundamentals']);" +#~ # hax stage +#~ HAXPYTHON="import hax; " +#~ HAXPYTHON+="hax.init(main_data_paths=['${OUTDIR}'], minitree_paths=['${OUTDIR}'], pax_version_policy = 'loose'); " +#~ HAXPYTHON+="hax.minitrees.load('${PAX_FILENAME##*/}', ['Basics', 'Fundamentals']);" + +#~ (time python -c "${HAXPYTHON}";) &> ${HAX_FILENAME}.log + + +# custom minitree +(time python ${RELEASEDIR}/ReduceDataNormal.py ${PAX_FILENAME_WOPATH} ${OUTDIR};) &> ${HAX_FILENAME}.log + +# merge +(time python ${RELEASEDIR}/TruthSorting.py ${FAX_FILENAME}.root ${MERGEDTRUTH_FILENAME}.pkl;) &> ${MERGEDTRUTH_FILENAME}.log +(time python ${RELEASEDIR}/MergeTruthAndProcessed.py ${RELEASEDIR}/Configs/QingConfig ${MERGEDTRUTH_FILENAME}.pkl ${MINITREE_FILENAME} ${MERGED_FILENAME}.pkl;) &> ${MERGED_FILENAME}.log -(time python -c "${HAXPYTHON}";) &> ${HAX_FILENAME}.log # Cleanup rm -f pax* diff --git a/montecarlo/fax_waveform/run_fax_AreaCorrelatedS1S2.sh b/montecarlo/fax_waveform/run_fax_AreaCorrelatedS1S2.sh new file mode 100644 index 0000000..8e92f27 --- /dev/null +++ b/montecarlo/fax_waveform/run_fax_AreaCorrelatedS1S2.sh @@ -0,0 +1,141 @@ +#!/usr/bin/env bash +############################################ +# +# Usage: Need to modify all parameters and paths below, then: +# ./run_fax.sh +# +############################################ + +echo "Start time: " `/bin/date` +echo "Job is running on node: " `/bin/hostname` +echo "Job running as user: " `/usr/bin/id` +echo "Job is running in directory: $PWD" + +###### General parameters ##### +Detector=XENON1T + +###### Simulation parameters ##### +Input2DBandFile=$1 +Nomial_g1=$2 +Nomial_g2=$3 + +RecoilType=ER + +PMTAfterpulseEnableFlag=$4 +S2AfterpulseEnableFlag=$5 + +# Select fax+pax version +#~ PAXVERSION=head +PAXVERSION=v6.5.1 # temporary change + +# Specify number of events +NumEvents=$6 + +# This run number (from command line argument) +SUBRUN=$8 + +######################################## + +# Setup the software +CVMFSDIR=/cvmfs/xenon.opensciencegrid.org +#export PATH="${CVMFSDIR}/releases/anaconda/2.4/bin:$PATH" +export PATH="/project/lgrandi/anaconda3/bin:$PATH" +source activate pax_${PAXVERSION} &> /dev/null + +# Use path of this script for Python scripts below +# (In case user modified them) +#MY_PATH=`dirname \"$0\"` +MY_PATH=$(cd `dirname $0`; pwd) +echo $MY_PATH +RELEASEDIR=`( cd "$MY_PATH" && pwd )` + +# Setting up directories +#start_dir=$PWD + + +OUTDIR=$7/${SUBRUN} +mkdir -p ${OUTDIR} +cd ${OUTDIR} + +#if [ "$OSG_WN_TMP" == "" ]; +#then +# OSG_WN_TMP=$PWD +#fi +#cd $OSG_WN_TMP +# +#work_dir=`mktemp -d --tmpdir=$OSG_WN_TMP` +#cd $work_dir + +# Filenaming +FILEROOT=FakeWaveform_${Detector}_${SUBRUN} +FILENAME=${OUTDIR}/${FILEROOT} +CSV_FILENAME=${FILENAME}.csv # Fake input data +FAX_FILENAME=${FILENAME}_truth # fax truth info +PKL_FILENAME=${FILENAME}_truth.pkl # converted fax truth info +RAW_FILENAME=${FILENAME}_raw # fax simulated raw data +PAX_FILENAME=${FILENAME}_pax # pax processed data +PAX_FILENAME_WOPATH=${FILEROOT}_pax # pax processed data without path +MERGEDTRUTH_FILENAME=${FILENAME}_merged_truth +MERGED_FILENAME=${FILENAME}_merged +HAX_FILENAME=${FILENAME}_hax # hax reduced data +MINITREE_FILENAME=${FILENAME}_pax_S1S2Properties.root +CustomIniFilename=${RELEASEDIR}/NoS2Afterpulses.ini +NoPMTAfterpulseIniFilename=${RELEASEDIR}/NoPMTAfterpulses.ini +echo ${CustomIniFilename} + + +# Create the fake input data +python ${RELEASEDIR}/CreateFakeCSV_CorrelatedS1S2.py ${Detector} ${Input2DBandFile} ${Nomial_g1} ${Nomial_g2} ${NumEvents} ${RecoilType} ${CSV_FILENAME} + +# Start of simulations # + +# fax stage +if (($S2AfterpulseEnableFlag==0)); then + if (($PMTAfterpulseEnableFlag==0)); then + echo 'Both S2 and PMT afterpulse disabled' + (time paxer --input ${CSV_FILENAME} --config ${Detector} reduce_raw_data Simulation --config_path ${NoPMTAfterpulseIniFilename} ${CustomIniFilename} --config_string "[WaveformSimulator]truth_file_name=\"${FAX_FILENAME}\"" --output ${RAW_FILENAME};) &> ${RAW_FILENAME}.log + else + echo 'Only S2 afterpulse disabled' + (time paxer --input ${CSV_FILENAME} --config ${Detector} reduce_raw_data Simulation --config_path ${CustomIniFilename} --config_string "[WaveformSimulator]truth_file_name=\"${FAX_FILENAME}\"" --output ${RAW_FILENAME};) &> ${RAW_FILENAME}.log + fi +else + if (($PMTAfterpulseEnableFlag==0)); then + echo 'Only PMT afterpulse disabled' + (time paxer --input ${CSV_FILENAME} --config ${Detector} reduce_raw_data Simulation --config_path ${NoPMTAfterpulseIniFilename} --config_string "[WaveformSimulator]truth_file_name=\"${FAX_FILENAME}\"" --output ${RAW_FILENAME};) &> ${RAW_FILENAME}.log + else + echo 'Both S2 and PMT afterpulse enabled' + (time paxer --input ${CSV_FILENAME} --config ${Detector} reduce_raw_data Simulation --config_string "[WaveformSimulator]truth_file_name=\"${FAX_FILENAME}\"" --output ${RAW_FILENAME};) &> ${RAW_FILENAME}.log + fi +fi + +# (time paxer --input ${CSV_FILENAME} --config ${Detector} reduce_raw_data Simulation --config_string "[WaveformSimulator]truth_file_name=\"${FAX_FILENAME}\"" --output ${RAW_FILENAME};) &> ${RAW_FILENAME}.log + + +# convert fax truth to pickle +python ${RELEASEDIR}/ConvertFaxTruthToPickle.py ${FAX_FILENAME} ${PKL_FILENAME} + +# pax stage +(time paxer --ignore_rundb --input ${RAW_FILENAME} --config ${Detector} --output ${PAX_FILENAME};) &> ${PAX_FILENAME}.log + +#~ # hax stage +#~ HAXPYTHON="import hax; " +#~ HAXPYTHON+="hax.init(main_data_paths=['${OUTDIR}'], minitree_paths=['${OUTDIR}'], pax_version_policy = 'loose'); " +#~ HAXPYTHON+="hax.minitrees.load('${PAX_FILENAME##*/}', ['Basics', 'Fundamentals']);" + +#~ (time python -c "${HAXPYTHON}";) &> ${HAX_FILENAME}.log + + +# custom minitree +(time python ${RELEASEDIR}/ReduceDataNormal.py ${PAX_FILENAME_WOPATH} ${OUTDIR};) &> ${HAX_FILENAME}.log + +# merge +(time python ${RELEASEDIR}/TruthSorting.py ${FAX_FILENAME}.root ${MERGEDTRUTH_FILENAME}.pkl;) &> ${MERGEDTRUTH_FILENAME}.log +(time python ${RELEASEDIR}/MergeTruthAndProcessed.py ${RELEASEDIR}/Configs/QingConfig ${MERGEDTRUTH_FILENAME}.pkl ${MINITREE_FILENAME} ${MERGED_FILENAME}.pkl;) &> ${MERGED_FILENAME}.log + + +# Cleanup +rm -f pax* + + +#cd $start_dir +#rm -fr $work_dir