Pax v640 newtruth qing

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("<detector: XENON100, XENON1T>")
+    print("<pkl file for Log(S2/S1) verse S1>")
+    print("<nominal g1 value>")
+    print("<nominal g2 value>")
+    print("<number of events>")
+    print("<recoil type: ER, NR>")
+    print("<output file (abs. path)>")
+    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<R2upper)
+        if(Zcut & Rcut):
+            return True
+
+    return False
+
+
+def RandomizeFV():
+
+    # randomize the X, Y, Z according to X48kg FV
+    if Detector == "XENON100":
+        Zlower, Zupper = -14.6-15.0, -14.6+15.0
+        Rlower, Rupper = -np.sqrt(200.), np.sqrt(200.)
+
+    elif Detector == "XENON1T": # NEED TO UPDATE THIS
+        Zlower, Zupper = -90*scalecmtomm, -15*scalecmtomm
+        Rlower, Rupper = -46*scalecmtomm, 46*scalecmtomm
+
+    for i in range(100000):
+        x = np.random.uniform(Rlower,Rupper)
+        y = np.random.uniform(Rlower,Rupper)
+        z = np.random.uniform(Zlower,Zupper)
+        if IfPassFV(x,y,z):
+            return (x,y,z)
+    return (0,0,0)
+
+####################################
+## S1&S2 generator loading
+####################################
+import PhotonChargeGenerator
+from PhotonChargeGenerator import MyPhotonChargeGenerator
+
+pGen = MyPhotonChargeGenerator(
+                                                          Band2DPDFFilename,
+                                                          NominalG1Value,
+                                                          NominalG2Value,
+                                                         )
+
+
+####################################
+## Starts to create
+####################################
+# Some default
+DefaultEventTime = MaxDriftTime*1000.
+##########
+fout = open(OutputFilename, 'w')
+# headers
+fout.write("instruction,recoil_type,x,y,depth,s1_photons,s2_electrons,t\n")
+# events loop
+for i in range(NumEvents):
+    fout.write(str(i)+",")
+    fout.write(DefaultType+",")
+    X, Y, Z = RandomizeFV()
+    fout.write(str(X)+",")
+    fout.write(str(Y)+",")
+    fout.write(str(-Z)+",")
+    NumPhoton, NumElectron = pGen.GetPhotonChargeNum()
+    fout.write(str(int(NumPhoton))+",")
+    fout.write(str(int(NumElectron))+",")
+    fout.write(str(DefaultEventTime)+"\n")
+

montecarlo/fax_waveform/MakePickleFromHistogram.py

@@ -0,0 +1,56 @@
+####################################
+## This code translate the histogram to pickle format
+####################################
+
+import numpy as np
+import pickle
+
+import ROOT
+from ROOT import TFile
+from ROOT import TH2D
+
+import sys, os
+
+if len(sys.argv)<2:
+    print("========== Syntax ===========")
+    print("python MakePickleFromHistogram.py .....")
+    print("<input root file>")
+    print("<histogram name>")
+    print("<output pickle file>")
+    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') )

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("<Output path (abs.)>")
+    print("<number of jobs>")
+    print("<number of events in each job>")
+    print("<enable PMT after pulses ?(0 for disable)>")
+    print("<enable S2 after pulses ?(0 for disable)>")
+    print("<2D band for S1-S2 area correlation (abs.)>")
+    print("<g1 value>")
+    print("<g2 value>")
+    print("<If use Public node (0 for no(xenon1t nodes); 1 for yes; 2 for kicp nodes)>")
+    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<MaxNumJob:
+            #sbatch it 
+            subp.call("cd "+SubmitPath+";sbatch "+SubmitFile+";cd -", shell=True)
+            IfSubmitted=1
+            time.sleep(2.0)
+        else:
+            time.sleep(30) 

montecarlo/fax_waveform/PhotonChargeGenerator.py

@@ -0,0 +1,101 @@
+####################################
+## This code generate number of photons and charges
+## according to a 2-D pdf input
+## by Qing Lin
+## Created @ 2017-02-21 
+####################################
+
+import sys, os
+import pickle
+import numpy as np
+import scipy as sp
+from scipy.interpolate import interp1d
+
+class MyPhotonChargeGenerator:
+    def __init__(self, Input2DFilename, g1, g2):
+        self.m_pG1Value = g1
+        self.m_pG2Value = g2
+        self.Load(Input2DFilename)
+        return
+
+    def __eq__(self, other):
+        self.m_pG1Value = other.m_pG1Value
+        self.m_pG2Value = other.m_pG2Value
+        # range of S1 and log
+        self.m_pS1Nbins = other.m_pS1Nbins
+        self.m_pS1Lower = other.m_pS1Lower
+        self.m_pS1Upper = other.m_pS1Upper
+        self.m_pS1Step = other.m_pS1Step
+        self.m_pLogNbins = other.m_pLogNbins
+        self.m_pLogLower = other.m_pLogLower
+        self.m_pLogUpper = other.m_pLogUpper
+        self.m_pLogStep = other.m_pLogStep
+        # interpolators
+        self.m_pS1Interpolator = other.m_pS1Interpolator
+        self.m_pLogInterpolators = other.m_pLogInterpolators
+        return
+
+
+    def Load(self, Input2DFilename):
+        data = pickle.load( open(Input2DFilename, 'rb') )
+        if not self.CheckIfComplete(data):
+            raise ValueError("Pickle file not complete")
+        self.m_pS1Nbins = data['s1nbins']
+        self.m_pS1Lower = data['s1lower']
+        self.m_pS1Upper = data['s1upper']
+        self.m_pS1Step = (self.m_pS1Upper - self.m_pS1Lower) / float( self.m_pS1Nbins)
+        self.m_pLogNbins = data['lognbins']
+        self.m_pLogLower = data['loglower']
+        self.m_pLogUpper = data['logupper']
+        self.m_pLogStep = (self.m_pLogUpper - self.m_pLogLower) / float( self.m_pLogNbins)
+        # load S1 pdf and calculate cdf
+        S1s = np.linspace(self.m_pS1Lower, self.m_pS1Upper, self.m_pS1Nbins+1)
+        S1s = S1s[1:]
+        S1PDFs = [ np.sum(A) for A in data['map'] ]
+        self.m_pS1Interpolator = interp1d(np.cumsum(S1PDFs)/np.sum(S1PDFs), S1s, bounds_error=False, fill_value=(0, 1))
+        # load S2 pdfs and calculate cdf
+        Logs = np.linspace(self.m_pLogLower, self.m_pLogUpper, self.m_pLogNbins+1)
+        Logs = Logs[1:]
+        self.m_pLogInterpolators = []
+        for LogPDFs in data['map']:
+            self.m_pLogInterpolators.append(
+                interp1d(np.cumsum(LogPDFs)/np.sum(LogPDFs), Logs, bounds_error=False, fill_value=(0,1) )
+                )
+        return
+
+
+    def CheckIfComplete(self, Data):
+        Items = [
+                      's1nbins',
+                      's1lower',
+                      's1upper',
+                      'lognbins',
+                      'loglower',
+                      'logupper',
+                      'map',
+                     ]
+        for item in Items:
+            if item not in Data:
+                return False
+        return True
+
+    def GetPhotonChargeNum(self):
+        s1_cdf = np.random.uniform(0, 1)
+        s1 = self.m_pS1Interpolator(s1_cdf)
+        s1_bin = int( (s1 - self.m_pS1Lower) / self.m_pS1Step )
+        s1_bin_upper = s1_bin + 1
+        if s1_bin_upper>=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)
+
+
+

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

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