Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-27 03:18:06

0001 #!/usr/bin/env python
0002 #____________________________________________________________
0003 #
0004 #  PlotLumiScan
0005 #
0006 # A very simple way to make lumiscan plots from beamfit.txt files
0007 # Needed files: A txt file specifying lumi section ranges eg. RunLumiScan.txt
0008 #   All the beam fit txt files in <data_dir> created after running AnalyzeLumiScan.py
0009 #
0010 # Geng-yuan Jeng
0011 # Geng-yuan.Jeng@cern.ch
0012 #
0013 # Fermilab, 2009
0014 #
0015 #____________________________________________________________
0016 
0017 from builtins import range
0018 import sys,os
0019 import string
0020 from array import array
0021 from math import sqrt
0022 from ROOT import gROOT,gPad,TH1F,TLegend,TFile,TStyle,TAxis
0023 
0024 def get_list_files(directory,pattern="txt"):
0025     dir = []
0026     dir = os.listdir(directory)
0027     dir.sort(cmp)
0028     lfiles = []
0029 
0030     for f in dir:
0031         if f.find(pattern) != -1:
0032         #print f
0033             lfiles.append(f)
0034     return lfiles;
0035 
0036 
0037 def plot_trending(type,label,x,xe):
0038     h = TH1F(type+'_lumi',type+'_lumi',len(x),0.5,0.5+len(x))
0039     h.SetStats(0)
0040     h.GetXaxis().SetTitle("Lumisection")
0041     h.GetYaxis().SetTitle(type[:len(type)-1]+"_{0} (cm)")
0042     h.GetYaxis().SetLabelSize(0.03)
0043     h.SetTitleOffset(1.1)
0044     h.SetOption("e1")
0045     
0046     for i in range(len(x)):
0047         h.SetBinContent(i+1, x[i])
0048         h.SetBinError(i+1, xe[i])
0049         h.GetXaxis().SetBinLabel(i+1,label[i])
0050     return h
0051 
0052 def main():
0053     
0054     if len(sys.argv) < 4:
0055         print("\n [Usage] python PlotLumiScan.py <LumiScanLists.txt> <data dir> <verbose:True/False>")
0056         sys.exit()
0057 
0058     lumilistfile = sys.argv[1]
0059     runinfofile = open(lumilistfile,"r")
0060     runinfolist = runinfofile.readlines()
0061     runsinfo = {}
0062 
0063     for line in runinfolist:
0064         npos=0
0065         for i in line.split():
0066             npos+=1
0067             if npos == 1:
0068                 run="Run"+str(i)+"/"
0069             else:
0070                 runsinfo.setdefault(run,[]).append(int(i))
0071 ##    print runsinfo
0072 
0073     infiledir = sys.argv[2]
0074     if infiledir.endswith("/") != 1:
0075         infiledir+="/"
0076     verbose = sys.argv[3]
0077     files = get_list_files(infiledir,"txt")
0078     nfiles = len(files)-1
0079     labels = []
0080     x0=[]
0081     y0=[]
0082     z0=[]
0083     sigZ=[]
0084     x0Err=[]
0085     y0Err=[]
0086     z0Err=[]
0087     sigZErr=[]
0088 
0089     ## Read files and put values into data containers
0090     ## Labels:
0091     
0092     lumilist = runsinfo.get(infiledir)
0093 
0094     for j in range((len(lumilist)+1)/2):
0095         labelName=str(lumilist[j*2])+"-"+str(lumilist[j*2+1])
0096         labels.append(labelName)
0097 ##    print labels
0098 
0099     for f in files:
0100         readfile = open(infiledir+f)
0101         for line in readfile:
0102             if line.find("X") != -1 and not "BeamWidth" in line and not "Emittance" in line:
0103                 count=0
0104                 for val in line.split():
0105                     count+=1
0106                     if count > 1:
0107                         x0.append(float(val))
0108             if line.find("Cov(0,j)") != -1:
0109                 count=0
0110                 for val in line.split():
0111                     count+=1
0112                     if count == 2:
0113                         valErr=sqrt(float(val))
0114                         x0Err.append(valErr)
0115             if line.find("Y") != -1 and not "BeamWidth" in line and not "Emittance" in line:
0116                 count=0
0117                 for val in line.split():
0118                     count+=1
0119                     if count > 1:
0120                         y0.append(float(val))
0121             if line.find("Cov(1,j)") != -1:
0122                 count=0
0123                 for val in line.split():
0124                     count+=1
0125                     if count == 3:
0126                         valErr=sqrt(float(val))
0127                         y0Err.append(valErr)
0128             if line.find("Z") != -1 and not "sigma" in line:
0129                 count=0
0130                 for val in line.split():
0131                     count+=1
0132                     if count > 1:
0133                         z0.append(float(val))
0134             if line.find("Cov(2,j)") != -1:
0135                 count=0
0136                 for val in line.split():
0137                     count+=1
0138                     if count == 4:
0139                         valErr=sqrt(float(val))
0140                         z0Err.append(valErr)
0141             if line.find("sigmaZ") != -1:
0142                 count=0
0143                 for val in line.split():
0144                     count+=1
0145                     if count > 1:
0146                         sigZ.append(float(val))
0147             if line.find("Cov(3,j)") != -1:
0148                 count=0
0149                 for val in line.split():
0150                     count+=1
0151                     if count == 5:
0152                         valErr=sqrt(float(val))
0153                         sigZErr.append(valErr)
0154 
0155     if verbose == "True":
0156         for i in range(len(x0)):
0157             print("     x0 = "+str(x0[i])+" +/- %1.8f (stats) [cm]" % (x0Err[i]))
0158             print("     y0 = "+str(y0[i])+" +/- %1.8f (stats) [cm]" % (y0Err[i]))
0159             print("     z0 = "+str(z0[i])+" +/- %1.6f (stats) [cm]" % (z0Err[i]))
0160             print("sigmaZ0 = "+str(sigZ[i])+" +/- %1.6f (stats) [cm]" % (sigZErr[i]))
0161 
0162     ## Make plots and save to root file
0163     rootFile = TFile("Summary.root","RECREATE");
0164     gROOT.SetStyle("Plain")
0165     
0166     hx0_lumi=plot_trending("x0",labels,x0,x0Err)
0167     hx0_lumi.SetTitle("x coordinate of beam spot vs. lumi")
0168     
0169     hy0_lumi=plot_trending("y0",labels,y0,y0Err)
0170     hy0_lumi.SetTitle("y coordinate of beam spot vs. lumi")
0171 
0172     hz0_lumi=plot_trending("z0",labels,z0,z0Err)
0173     hz0_lumi.SetTitle("z coordinate of beam spot vs. lumi")
0174 
0175     hsigZ_lumi=plot_trending("sigmaZ0",labels,sigZ,sigZErr)
0176     hsigZ_lumi.SetTitle("sigma z_{0} of beam spot vs. lumi")
0177     
0178     rootFile.Write();
0179     rootFile.Close();
0180     
0181 
0182 #_________________________________    
0183 if __name__ =='__main__':
0184     sys.exit(main())
0185 
0186