Pax v640 newtruth qing

montecarlo/fax_waveform/BatchMergeTruthAndProcessed.py

@@ -14,7 +14,7 @@
 
 if len(sys.argv)<=1:
     print("======== Usage =========")
-    print("python BatchMergeTruthAndProcessed.py <config file> <truth root path> <processed root path> <output path> <(opt)top-to-total fraction in truth(default 0.68)> <(opt)relative path for submission> <(opt) if use public node (1) optional (2 for use kicp nodes)> <(opt)Submit ID> <(opt) if use arrays in output (1) (default 0)>")
+    print("python BatchMergeTruthAndProcessed.py <config file> <truth root path> <processed root path> <output path> <(opt)relative path for submission> <(opt) if use public node (1) optional (2 for use kicp nodes)> <(opt)Submit ID> <(opt) if use arrays in output (1) (default 0)> <(opt) minitree type; 0(default): Basics, 1: S1S2Properties, 2: PeakEfficiency")
     exit()
 
 CurrentEXE = sys.argv[0]
@@ -23,20 +23,20 @@
 ProcessedRootPath = sys.argv[3]
 OutputPath = sys.argv[4]
 IfPublicNode = 1
-TopFraction = 0.68
-if len(sys.argv)>5:
-    TopFraction = float(sys.argv[5])
 RelativeSubmitPath = "Submit"
+if len(sys.argv)>5:
+    RelativeSubmitPath = sys.argv[5]
 if len(sys.argv)>6:
-    RelativeSubmitPath = sys.argv[6]
-if len(sys.argv)>7:
-    IfPublicNode = int(sys.argv[7])
+    IfPublicNode = int(sys.argv[6])
 SubmitID = 0
-if len(sys.argv)>8:
-    SubmitID = int(sys.argv[8])
+if len(sys.argv)>7:
+    SubmitID = int(sys.argv[7])
 ArrayOutput = 0
+if len(sys.argv)>8:
+    ArrayOutput = int(sys.argv[8])
+MinitreeType=0
 if len(sys.argv)>9:
-    ArrayOutput = int(sys.argv[9])
+    MinitreeType = int(sys.argv[9])
 
 
 #######################
@@ -52,7 +52,10 @@
     EXE1 = CurrentPath+"/"+EXE_Path+"TruthSorting_arrays.py"
 else:
     EXE1 = CurrentPath+"/"+EXE_Path+"TruthSorting.py"
-EXE2 = CurrentPath+"/"+EXE_Path+"MergeTruthAndProcessed.py"
+if MinitreeType==2:
+    EXE2 = CurrentPath+"/"+EXE_Path+"MergeTruthAndProcessed_peaks.py"
+else:
+    EXE2 = CurrentPath+"/"+EXE_Path+"MergeTruthAndProcessed.py"
 
 
 #######################
@@ -102,7 +105,7 @@
     subp.call("echo '#!/bin/bash\n' >> "+SubmitFile, shell=True)
     subp.call("echo '#SBATCH --output="+SubmitPath+"/myout_"+str(SubmitID)+"_"+str(j)+".txt \n' >> "+SubmitFile, shell=True)
     subp.call("echo '#SBATCH --error="+SubmitPath+"/myerr_"+str(SubmitID)+"_"+str(j)+".txt\n' >> "+SubmitFile, shell=True)
-    subp.call("echo '#SBATCH --time=01:19:00\n' >> "+SubmitFile, shell=True)
+    subp.call("echo '#SBATCH --time=00:05:00\n' >> "+SubmitFile, shell=True)
     if not IfPublicNode:
         subp.call("echo '#SBATCH --account=pi-lgrandi\n' >> "+SubmitFile, shell=True)
         subp.call("echo '#SBATCH --qos=xenon1t\n' >> "+SubmitFile, shell=True)
@@ -112,7 +115,7 @@
         subp.call("echo '#SBATCH --qos=xenon1t-kicp\n' >> "+SubmitFile, shell=True)
         subp.call("echo '#SBATCH --partition=kicp\n' >> "+SubmitFile, shell=True)
     subp.call("echo '. /home/mcfate/Env/GlobalPAXEnv.sh\n\n' >> "+SubmitFile, shell=True)
-    subp.call("echo 'python "+EXE1+" "+TruthRootFilename+" "+TmpOutputFilename+" "+str(TopFraction)+"' >> "+SubmitFile, shell=True)
+    subp.call("echo 'python "+EXE1+" "+TruthRootFilename+" "+TmpOutputFilename+"' >> "+SubmitFile, shell=True)
     subp.call("echo 'python "+EXE2+" "+AbsoluteConfigFile+" "+TmpOutputFilename+" "+ProcessedRootFilename+" "+OutputFilename+"' >> "+SubmitFile, shell=True)
 
     #submit

montecarlo/fax_waveform/BatchReduceDataSubmission.py

@@ -13,7 +13,7 @@
 
 if len(sys.argv)<=1:
     print("======== Usage =========")
-    print("python ReduceDataNormal.py <filelist> <data path> <output path> <absolute path for submission> <if use public node (1) optional (2 for use kicp nodes)> <Submit ID>")
+    print("python BatchReduceDataSubmission.py <filelist> <data path> <output path> <absolute path for submission> <if use public node (1) optional (2 for use kicp nodes)> <Submit ID> <(opt) minitree type; 1: S1S2Properties, 2: PeakEfficiency>")
     print("======== List file format: ==========")
     print("ex.:")
     print("FakeWaveform_XENON1T_000000_pax")
@@ -35,13 +35,20 @@
 SubmitID = 0
 if len(sys.argv)>6:
     SubmitID = int(sys.argv[6])
+minitree_type = '1'
+if len(sys.argv)>7:
+    minitree_type = sys.argv[7]
+
 
 ##########################
 ## Some nuisance settings
 ##########################
 CurrentPath = os.getcwd()
 CurrentUser = getpass.getuser()
-EXE = CurrentPath+"/"+EXE_Path+"/ReduceDataNormal.py"
+if minitree_type == '1':
+    EXE = CurrentPath+"/"+EXE_Path+"/ReduceDataNormal.py"
+elif minitree_type == '2':
+    EXE = CurrentPath+"/"+EXE_Path+"/reduce_peak_level.py"
 MaxNumJob = 64
 if not IfPublicNode:
     MaxNumJob = 200
@@ -87,7 +94,10 @@
         subp.call("echo '#SBATCH --partition=kicp\n' >> "+SubmitFile, shell=True)
     subp.call("echo '. /home/mcfate/Env/GlobalPAXEnv.sh\n\n' >> "+SubmitFile, shell=True)
     subp.call("echo 'python "+EXE+" "+filename+" "+DataPath+"' >> "+SubmitFile, shell=True)
-    subp.call("echo 'mv "+SubmitPath+"/"+filename+"_S1S2Properties.root  "+OutputPath+"' >> "+SubmitFile, shell=True)
+    if minitree_type=='1':
+        subp.call("echo 'mv "+SubmitPath+"/"+filename+"_S1S2Properties.root  "+OutputPath+"' >> "+SubmitFile, shell=True)
+    elif minitree_type=='2':
+        subp.call("echo 'mv "+SubmitPath+"/"+filename+"_PeakEfficiency.root  "+OutputPath+"' >> "+SubmitFile, shell=True)
 
     #submit
     IfSubmitted=0
@@ -114,6 +124,3 @@
         else:
             time.sleep(30) 
 
-
-
-

montecarlo/fax_waveform/Configs/PeakEfficiency

@@ -0,0 +1,13 @@
+Name: PeakEfficiency
+event_number event_number_processed
+run_number run_number
+area area
+hit_time_mean hit_time_mean
+hit_time_std hit_time_std
+type type
+area_fraction_top area_fraction_top
+height height
+n_contributing_channels n_contributing_channels
+range_50p_area range_50p_area
+range_70p_area range_70p_area
+range_90p_area range_90p_area

montecarlo/fax_waveform/Configs/basics_config

@@ -0,0 +1,18 @@
+Name: Basics
+cs1 cs1_processed
+cs2 cs2_processed
+drift_time drift_time_processed
+event_number event_number_processed
+largest_coincidence largest_coincidence_processed
+largest_other_s1 largest_other_s1_processed
+largest_other_s2 largest_other_s2_processed
+run_number run_number_processed
+s1 s1_processed
+s1_area_fraction_top s1_area_fraction_top_processed
+s1_range_50p_area s1_range_50p_area_processed
+s2 s2_processed
+s2_area_fraction_top s2_area_fraction_top_processed
+s2_range_50p_area s2_range_50p_area_processed
+x x_processed
+y y_processed
+z z_processed

montecarlo/fax_waveform/Configs/s1s2_preserve_all

@@ -0,0 +1,57 @@
+Name: S1S2Properties
+index index
+NbPeaks NbPeaks
+NbS1s NbS1s
+NbS2s NbS2s
+S1PatternLnL S1PatternLnL
+S1S2InMainInteraction S1S2InMainInteraction
+S1TopFraction S1TopFraction
+S1_X_PosRecMaxPMT S1_X_PosRecMaxPMT
+S1_X_PosRecNeuralNet S1_X_PosRecNeuralNet
+S1_X_PosRecRobustWeightedMean S1_X_PosRecRobustWeightedMean
+S1_X_PosRecTopPatternFit S1_X_PosRecTopPatternFit
+S1_X_PosRecWeightedSum S1_X_PosRecWeightedSum
+S1_Y_PosRecMaxPMT S1_Y_PosRecMaxPMT
+S1_Y_PosRecNeuralNet S1_Y_PosRecNeuralNet
+S1_Y_PosRecRobustWeightedMean S1_Y_PosRecRobustWeightedMean
+S1_Y_PosRecTopPatternFit S1_Y_PosRecTopPatternFit
+S1_Y_PosRecWeightedSum S1_Y_PosRecWeightedSum
+S1sCoin S1sCoin
+S1sCorrection S1sCorrection
+S1sHeight S1sHeight
+S1sLowWidth S1sLowWidth
+S1sNbSaturationChannels S1sNbSaturationChannels
+S1sPeakTime S1sPeakTime
+S1sPeakTimeStd S1sPeakTimeStd
+S1sTot S1sTot
+S1sTotGained S1sTotGained
+S1sTotSecond S1sTotSecond
+S1sTotSecondGained S1sTotSecondGained
+S1sWidth S1sWidth
+S2PosGoodnessOfFit S2PosGoodnessOfFit
+S2TopFraction S2TopFraction
+S2_X_PosRecMaxPMT S2_X_PosRecMaxPMT
+S2_X_PosRecNeuralNet S2_X_PosRecNeuralNet
+S2_X_PosRecRobustWeightedMean S2_X_PosRecRobustWeightedMean
+S2_X_PosRecTopPatternFit S2_X_PosRecTopPatternFit
+S2_X_PosRecWeightedSum S2_X_PosRecWeightedSum
+S2_Y_PosRecMaxPMT S2_Y_PosRecMaxPMT
+S2_Y_PosRecNeuralNet S2_Y_PosRecNeuralNet
+S2_Y_PosRecRobustWeightedMean S2_Y_PosRecRobustWeightedMean
+S2_Y_PosRecTopPatternFit S2_Y_PosRecTopPatternFit
+S2_Y_PosRecWeightedSum S2_Y_PosRecWeightedSum
+S2sCoin S2sCoin
+S2sCorrection S2sCorrection
+S2sHeight S2sHeight
+S2sLowWidth S2sLowWidth
+S2sNbSaturationChannels S2sNbSaturationChannels
+S2sPeakTime S2sPeakTime
+S2sPeakTimeStd S2sPeakTimeStd
+S2sTot S2sTot
+S2sTotGained S2sTotGained
+S2sTotSecond S2sTotSecond
+S2sTotSecondGained S2sTotSecondGained
+S2sWidth S2sWidth
+event_number event_number
+run_number run_number
+time time

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