BeamSpotProducer/scripts/plotBeamSpotDB.py

0001 #!/usr/bin/env python3
0002 #____________________________________________________________
0003 #
0004 #
0005 # A very simple way to make plots with ROOT via an XML file
0006 #
0007 # Francisco Yumiceva
0008 # yumiceva@fnal.gov
0009 #
0010 # Fermilab, 2010
0011 #
0012 #____________________________________________________________
0013
0014 """
0015    plotBeamSpotDB
0016
0017    A very simple script to plot the beam spot data stored in condDB
0018
0019    usage: %prog -t <tag name>
0020    -a, --auth    = AUTH: DB authorization path. online(/nfshome0/popcondev/conddb).
0021    -b, --batch : Run ROOT in batch mode.
0022    -c, --create  = CREATE: name for beam spot data file.
0023    -d, --data    = DATA: input beam spot data file.
0024    -D, --destDB  = DESTDB: destination DB string. online(oracle://cms_orcon_prod/CMS_COND_31X_BEAMSPOT).
0025    -i, --initial = INITIAL: First IOV. Options: run number, or run:lumi, eg. \"133200:21\"
0026    -f, --final   = FINAL: Last IOV. Options: run number, or run:lumi
0027    -g, --graph : create a TGraphError instead of a TH1 object
0028    -n, --noplot : Only extract beam spot data, plots are not created..
0029    -o, --output  = OUTPUT: filename of ROOT file with plots.
0030    -p, --payload = PAYLOAD: filename of output text file. Combine and splits lumi IOVs.
0031    -P, --Print : create PNG plots from canvas.
0032    -s, --suffix = SUFFIX: Suffix will be added to plots filename.
0033    -t, --tag     = TAG: Database tag name.
0034    -T, --Time : create plots with time axis.
0035    -I, --IOVbase = IOVBASE: options: runbase(default), lumibase, timebase
0036    -w, --wait : Pause script after plotting a new histograms.
0037    -W, --weighted : Create a weighted result for a range of lumi IOVs, skip lumi IOV combination and splitting.
0038    -x, --xcrossing = XCROSSING : Bunch crossing number.
0039
0040    Francisco Yumiceva (yumiceva@fnal.gov)
0041    Fermilab 2010
0042
0043 """
0044 from __future__ import print_function
0045
0046
0047 from builtins import range
0048 import os, string, re, sys, math
0049 import subprocess, time
0050 from BeamSpotObj import BeamSpot
0051 from IOVObj import IOV
0052 from CommonMethods import *
0053
0054 try:
0055     import ROOT
0056 except:
0057     print("\nCannot load PYROOT, make sure you have setup ROOT in the path")
0058     print("and pyroot library is also defined in the variable PYTHONPATH, try:\n")
0059     if (os.getenv("PYTHONPATH")):
0060         print(" setenv PYTHONPATH ${PYTHONPATH}:$ROOTSYS/lib\n")
0061     else:
0062         print(" setenv PYTHONPATH $ROOTSYS/lib\n")
0063         sys.exit()
0064
0065 from ROOT import TFile, TGraphErrors, TGaxis, TDatime
0066 from ROOT import TCanvas, TH1F
0067
0068 # ROOT STYLE
0069 #############################
0070 def SetStyle():
0071
0072     # canvas
0073     ROOT.gStyle.SetCanvasBorderMode(0)
0074     ROOT.gStyle.SetCanvasColor(0)
0075     ROOT.gStyle.SetCanvasDefH(600)
0076     ROOT.gStyle.SetCanvasDefW(600)
0077     ROOT.gStyle.SetCanvasDefX(0)
0078     ROOT.gStyle.SetCanvasDefY(0)
0079     # pad
0080     ROOT.gStyle.SetPadBorderMode(0)
0081     ROOT.gStyle.SetPadColor(0)
0082     ROOT.gStyle.SetPadGridX(False)
0083     ROOT.gStyle.SetPadGridY(False)
0084     ROOT.gStyle.SetGridColor(0)
0085     ROOT.gStyle.SetGridStyle(3)
0086     ROOT.gStyle.SetGridWidth(1)
0087
0088     ROOT.gStyle.SetFrameBorderMode(0)
0089     ROOT.gStyle.SetFrameFillColor(0)
0090     ROOT.gStyle.SetTitleColor(1)
0091     ROOT.gStyle.SetStatColor(0)
0092
0093     # set the paper & margin sizes
0094     ROOT.gStyle.SetPaperSize(20,26)
0095     ROOT.gStyle.SetPadTopMargin(0.04)
0096     ROOT.gStyle.SetPadRightMargin(0.04)
0097     ROOT.gStyle.SetPadBottomMargin(0.14)
0098     ROOT.gStyle.SetPadLeftMargin(0.11)
0099     ROOT.gStyle.SetPadTickX(1)
0100     ROOT.gStyle.SetPadTickY(1)
0101
0102     ROOT.gStyle.SetTextFont(42) #132
0103     ROOT.gStyle.SetTextSize(0.09)
0104     ROOT.gStyle.SetLabelFont(42,"xyz")
0105     ROOT.gStyle.SetTitleFont(42,"xyz")
0106     ROOT.gStyle.SetLabelSize(0.035,"xyz")
0107     ROOT.gStyle.SetTitleSize(0.045,"xyz")
0108     ROOT.gStyle.SetTitleOffset(1.1,"y")
0109
0110     # use bold lines and markers
0111     ROOT.gStyle.SetMarkerStyle(8)
0112     ROOT.gStyle.SetHistLineWidth(2)
0113     ROOT.gStyle.SetLineWidth(1)
0114     #ROOT.gStyle.SetLineStyleString(2,"[12 12]") // postscript dashes
0115
0116     ROOT.gStyle.SetMarkerSize(0.6)
0117
0118     # do not display any of the standard histogram decorations
0119     ROOT.gStyle.SetOptTitle(0)
0120     ROOT.gStyle.SetOptStat(0) #("m")
0121     ROOT.gStyle.SetOptFit(0)
0122
0123     #ROOT.gStyle.SetPalette(1,0)
0124     ROOT.gStyle.cd()
0125     ROOT.gROOT.ForceStyle()
0126 #########################################
0127
0128
0129 if __name__ == '__main__':
0130
0131     # style
0132     SetStyle()
0133
0134     printCanvas = False
0135     printFormat = "png"
0136     printBanner = False
0137     Banner = "CMS Preliminary"
0138
0139     # COMMAND LINE OPTIONS
0140     #################################
0141     option,args = parse(__doc__)
0142     if not args and not option: exit()
0143
0144     tag = ''
0145     if not option.tag and not option.data:
0146         print(" need to provide DB tag name or beam spot data file")
0147         exit()
0148     else:
0149         tag = option.tag
0150
0151     if option.batch:
0152         ROOT.gROOT.SetBatch()
0153
0154     datafilename = "tmp_beamspot.dat"
0155     if option.create:
0156         datafilename = option.create
0157
0158     getDBdata = True
0159     if option.data:
0160         getDBdata = False
0161
0162     IOVbase = 'runbase'
0163     if option.IOVbase:
0164         if option.IOVbase != "runbase" and option.IOVbase != "lumibase" and option.IOVbase != "timebase":
0165             print("\n\n unknown iov base option: "+ option.IOVbase +" \n\n\n")
0166             exit()
0167         IOVbase = option.IOVbase
0168
0169     firstRun = "0"
0170     lastRun  = "4999999999"
0171     if IOVbase == "lumibase":
0172         firstRun = "0:0"
0173         lastRun = "4999999999:4999999999"
0174
0175     if option.initial:
0176         firstRun = option.initial
0177     if option.final:
0178         lastRun = option.final
0179
0180     # GET IOVs
0181     ################################
0182
0183     if getDBdata:
0184
0185         print(" read DB to get list of IOVs for the given tag")
0186         mydestdb = 'frontier://PromptProd/CMS_COND_31X_BEAMSPOT'
0187         if option.destDB:
0188             mydestdb = option.destDB
0189         acommand = 'cmscond_list_iov -c '+mydestdb+' -P /afs/cern.ch/cms/DB/conddb -t '+ tag
0190         tmpstatus = subprocess.getstatusoutput( acommand )
0191         tmplistiov = tmpstatus[1].split('\n')
0192         #print tmplistiov
0193
0194         iovlist = []
0195         passline = False
0196         iline = jline = 0
0197         totlines = len(tmplistiov)
0198         for line in tmplistiov:
0199
0200             if line.find('since') != -1:
0201                 passline = True
0202                 jline = iline
0203             if passline and iline > jline and iline < totlines-1:
0204                 linedata = line.split()
0205                 #print linedata
0206                 aIOV = IOV()
0207                 aIOV.since = int(linedata[0])
0208                 aIOV.till = int(linedata[1])
0209                 iovlist.append( aIOV )
0210             iline += 1
0211
0212         print(" total number of IOVs = " + str(len(iovlist)))
0213
0214
0215         #  GET DATA
0216         ################################
0217
0218         otherArgs = ''
0219         if option.destDB:
0220             otherArgs = " -d " + option.destDB
0221             if option.auth:
0222                 otherArgs = otherArgs + " -a "+ option.auth
0223
0224         print(" get beam spot data from DB for IOVs. This can take a few minutes ...")
0225
0226         tmpfile = open(datafilename,'w')
0227
0228         for iIOV in iovlist:
0229             passiov = False
0230             tmprunfirst = firstRun
0231             tmprunlast = lastRun
0232             tmplumifirst = 1
0233             tmplumilast = 9999999
0234             if IOVbase=="lumibase":
0235                 #tmprunfirst = int(firstRun.split(":")[0])
0236                 #tmprunlast  = int(lastRun.split(":")[0])
0237                 #tmplumifirst = int(firstRun.split(":")[1])
0238                 #tmplumilast  = int(lastRun.split(":")[1])
0239                 tmprunfirst = pack( int(firstRun.split(":")[0]) , int(firstRun.split(":")[1]) )
0240                 tmprunlast  = pack( int(lastRun.split(":")[0]) , int(lastRun.split(":")[1]) )
0241             #print "since = " + str(iIOV.since) + " till = "+ str(iIOV.till)
0242             if iIOV.since >= int(tmprunfirst) and int(tmprunlast) < 0 and iIOV.since <= int(tmprunfirst):
0243                 print(" IOV: " + str(iIOV.since))
0244                 passiov = True
0245             if iIOV.since >= int(tmprunfirst) and int(tmprunlast) > 0 and iIOV.till <= int(tmprunlast):
0246                 print(" a IOV: " + str(iIOV.since) + " to " + str(iIOV.till))
0247                 passiov = True
0248             #if iIOV.since >= int(tmprunlast) and iIOV.till >= 4294967295:
0249             #    print " b IOV: " + str(iIOV.since) + " to " + str(iIOV.till)
0250             #    passiov = True
0251             if passiov:
0252                 acommand = 'getBeamSpotDB.py -t '+ tag + " -r " + str(iIOV.since) +otherArgs
0253                 if IOVbase=="lumibase":
0254                     tmprun = unpack(iIOV.since)[0]
0255                     tmplumi = unpack(iIOV.since)[1]
0256                     acommand = 'getBeamSpotDB.py -t '+ tag + " -r " + str(tmprun) +" -l "+str(tmplumi) +otherArgs
0257                     print(acommand)
0258                 status = subprocess.getstatusoutput( acommand )
0259                 tmpfile.write(status[1])
0260
0261         print(" beam spot data collected and stored in file " + datafilename)
0262
0263         tmpfile.close()
0264
0265
0266     # PROCESS DATA
0267     ###################################
0268
0269     # check if input data exists if given
0270     if option.data:
0271         if os.path.isdir(option.data):
0272             tmp = subprocess.getstatusoutput("ls "+option.data)
0273             files = tmp[1].split()
0274             datafilename = "combined_all.txt"
0275             output = open(datafilename,"w")
0276
0277             for f in files:
0278                 input = open(option.data +"/"+f)
0279                 output.writelines(input.readlines())
0280             output.close()
0281             print(" data files have been collected in "+datafilename)
0282
0283         elif os.path.exists(option.data):
0284             datafilename = option.data
0285         else:
0286             print(" input beam spot data DOES NOT exist, file " + option.data)
0287             exit()
0288
0289     listbeam = []
0290
0291     if option.xcrossing:
0292         listmap = readBeamSpotFile(datafilename,listbeam,IOVbase,firstRun,lastRun)
0293         # bx
0294         print("List of bunch crossings in the file:")
0295         print(listmap.keys())
0296         listbeam = listmap[option.Xrossing]
0297     else:
0298         readBeamSpotFile(datafilename,listbeam,IOVbase,firstRun,lastRun)
0299
0300     sortAndCleanBeamList(listbeam,IOVbase)
0301
0302     if IOVbase == "lumibase" and option.payload:
0303         weighted = True;
0304         if not option.weighted:
0305             weighted = False
0306         createWeightedPayloads(option.payload,listbeam,weighted)
0307     if option.noplot:
0308         print(" no plots requested, exit now.")
0309         sys.exit()
0310     # MAKE PLOTS
0311     ###################################
0312     TGaxis.SetMaxDigits(8)
0313
0314     graphlist = []
0315     graphnamelist = ['X','Y','Z','SigmaZ','dxdz','dydz','beamWidthX', 'beamWidthY']
0316     graphtitlelist = ['beam spot X','beam spot Y','beam spot Z','beam spot #sigma_Z','beam spot dX/dZ','beam spot dY/dZ','beam width X','beam width Y']
0317     graphXaxis = 'Run number'
0318     if IOVbase == 'runbase':
0319         graphXaxis = "Run number"
0320     if IOVbase == 'lumibase':
0321         graphXaxis = 'Lumi section'
0322     if IOVbase == 'timebase' or option.Time:
0323         graphXaxis = "Time"
0324         #dh = ROOT.TDatime(2010,06,01,00,00,00)
0325         ROOT.gStyle.SetTimeOffset(0) #dh.Convert())
0326
0327     graphYaxis = ['beam spot X [cm]','beam spot Y [cm]','beam spot Z [cm]', 'beam spot #sigma_{Z} [cm]', 'beam spot dX/dZ', 'beam spot dY/dZ','beam width X [cm]', 'beam width Y [cm]']
0328
0329     cvlist = []
0330
0331     for ig in range(0,8):
0332         cvlist.append( TCanvas(graphnamelist[ig],graphtitlelist[ig], 1200, 600) )
0333         if option.graph:
0334             graphlist.append( TGraphErrors( len(listbeam) ) )
0335         else:
0336             graphlist.append( TH1F("name","title",len(listbeam),0,len(listbeam)) )
0337
0338         graphlist[ig].SetName(graphnamelist[ig])
0339         graphlist[ig].SetTitle(graphtitlelist[ig])
0340         ipoint = 0
0341         for ii in range(0,len(listbeam)):
0342
0343             ibeam = listbeam[ii]
0344             datax = dataxerr = 0.
0345             datay = datayerr = 0.
0346             if graphnamelist[ig] == 'X':
0347                 datay = ibeam.X
0348                 datayerr = ibeam.Xerr
0349             if graphnamelist[ig] == 'Y':
0350                 datay = ibeam.Y
0351                 datayerr = ibeam.Yerr
0352             if graphnamelist[ig] == 'Z':
0353                 datay = ibeam.Z
0354                 datayerr = ibeam.Zerr
0355             if graphnamelist[ig] == 'SigmaZ':
0356                 datay = ibeam.sigmaZ
0357                 datayerr = ibeam.sigmaZerr
0358             if graphnamelist[ig] == 'dxdz':
0359                 datay = ibeam.dxdz
0360                 datayerr = ibeam.dxdzerr
0361             if graphnamelist[ig] == 'dydz':
0362                 datay = ibeam.dydz
0363                 datayerr = ibeam.dydzerr
0364             if graphnamelist[ig] == 'beamWidthX':
0365                 datay = ibeam.beamWidthX
0366                 datayerr = ibeam.beamWidthXerr
0367             if graphnamelist[ig] == 'beamWidthY':
0368                 datay = ibeam.beamWidthY
0369                 datayerr = ibeam.beamWidthYerr
0370
0371             datax = ibeam.IOVfirst
0372             if IOVbase=="lumibase":
0373                 datax = str(ibeam.Run) + ":" + str(ibeam.IOVfirst)
0374                 if ibeam.IOVfirst != ibeam.IOVlast:
0375                     datax = str(ibeam.Run) + ":" + str(ibeam.IOVfirst)+"-"+str(ibeam.IOVlast)
0376             #print datax
0377             if option.graph:
0378                 if IOVbase=="lumibase":
0379                     #first = int( pack( int(ibeam.Run) , int(ibeam.IOVfirst) ) )
0380                     #last = int( pack( int(ibeam.Run) , int(ibeam.IOVlast) ) )
0381                     first = ibeam.IOVfirst
0382                     last = ibeam.IOVlast
0383                     if option.Time:
0384                         atime = ibeam.IOVBeginTime
0385                         first = time.mktime( time.strptime(atime.split()[0] +  " " + atime.split()[1] + " " + atime.split()[2],"%Y.%m.%d %H:%M:%S %Z") )
0386                         atime = ibeam.IOVEndTime
0387                         last = time.mktime( time.strptime(atime.split()[0] +  " " + atime.split()[1] + " " + atime.split()[2],"%Y.%m.%d %H:%M:%S %Z") )
0388                         da_first = TDatime(time.strftime('%Y-%m-%d %H:%M:%S',time.localtime(first - time.timezone)))
0389                         da_last = TDatime(time.strftime('%Y-%m-%d %H:%M:%S',time.localtime(last - time.timezone)))
0390                         if ipoint == 0:
0391                             ## print local time
0392                             da_first.Print()
0393                             ## print gmt time
0394                             print("GMT = " + str(time.strftime('%Y-%m-%d %H:%M:%S',time.gmtime(first - time.timezone))))
0395                             reftime = first
0396                             ptm = time.localtime(reftime)
0397                             da = TDatime(time.strftime('%Y-%m-%d %H:%M:%S',ptm))
0398                             if time.daylight and ptm.tm_isdst:
0399                                 offset_daylight = time.timezone - time.altzone
0400                             ROOT.gStyle.SetTimeOffset(da.Convert(1) - 3600)
0401
0402                     datax = (float(last) - float(first))/2 + float(first) - da.Convert() + 3600
0403                     ## Comment out this block if running on Mac ##
0404                     if time.daylight and ptm.tm_isdst:
0405                         datax += offset_daylight
0406                     ##################################
0407
0408                     dataxerr =  (float(last) - float(first))/2
0409                 graphlist[ig].SetPoint(ipoint, float(datax), float(datay) )
0410                 graphlist[ig].SetPointError(ipoint, float(dataxerr), float(datayerr) )
0411             else:
0412                 graphlist[ig].GetXaxis().SetBinLabel(ipoint +1 , str(datax) )
0413                 graphlist[ig].SetBinContent(ipoint +1, float(datay) )
0414                 graphlist[ig].SetBinError(ipoint +1, float(datayerr) )
0415
0416             ipoint += 1
0417         if option.Time:
0418             ## print local time
0419             da_last.Print()
0420             print("GMT = " + str(time.strftime('%Y-%m-%d %H:%M:%S',time.gmtime(last - time.timezone))))
0421             graphlist[ig].GetXaxis().SetTimeDisplay(1);
0422             graphlist[ig].GetXaxis().SetTimeFormat("#splitline{%Y/%m/%d}{%H:%M}")
0423         if option.graph:
0424             graphlist[ig].Draw('AP')
0425         else:
0426             graphlist[ig].Draw('P E1 X0')
0427             graphlist[ig].GetXaxis().SetTitle(graphXaxis)
0428             graphlist[ig].GetYaxis().SetTitle(graphYaxis[ig])
0429             #graphlist[ig].Fit('pol1')
0430         cvlist[ig].Update()
0431         cvlist[ig].SetGrid()
0432         if option.Print:
0433             suffix = ''
0434             if option.suffix:
0435                 suffix = option.suffix
0436             cvlist[ig].Print(graphnamelist[ig]+"_"+suffix+".png")
0437         if option.wait:
0438             raw_input( 'Press ENTER to continue\n ' )
0439         #graphlist[0].Print('all')
0440
0441     if option.output:
0442         outroot = TFile(option.output,"RECREATE")
0443         for ig in graphlist:
0444             ig.Write()
0445
0446         outroot.Close()
0447         print(" plots have been written to "+option.output)
0448
0449
0450
0451     # CLEAN temporal files
0452     ###################################
0453     #os.system('rm tmp_beamspotdata.log')