root/BxAnalysis/lumiregperbunch.py

0001 #!/usr/bin/env python3
0002 from __future__ import print_function
0003 from builtins import range
0004 import sys,os,calendar
0005 from ROOT import gDirectory,TFile
0006
0007 from online_beamspot_reader import BeamspotMeasurement
0008 from online_beamspot_reader import paragraphs, \
0009      start_of_new_beamspot_measurement
0010
0011 FILL=''
0012 OUTDIR="./perbunchfiles/"
0013
0014 runlstime={}
0015
0016
0017 class bsmeas(object):
0018     def __init__(self, x=None,y=None,z=None,ex=None,ey=None,ez=None,
0019                  wx=None,wy=None,wz=None,ewx=None,ewy=None,ewz=None,
0020                  dxdz=None,dydz=None,edxdz=None,edydz=None):
0021         self.x = x
0022         self.y = y
0023         self.z = z
0024         self.ex = ex
0025         self.ey = ey
0026         self.ez = ez
0027         self.wx = ex
0028         self.wy = ey
0029         self.wz = ez
0030         self.ewx = ex
0031         self.ewy = ey
0032         self.ewz = ez
0033         self.dxdz = dxdz
0034         self.dydz = dydz
0035         self.edxdz = edxdz
0036         self.edydz = edydz
0037
0038
0039
0040
0041 def timeof(run,lumisection):
0042     # first check if this run is already in the list, otherwise read it
0043     if run not in runlstime.keys():
0044         print("Reading lumi time from lumireg localcopy files")
0045         filename="localcopy/BeamFitResults_Run"+run+".txt"
0046         if not os.path.exists(filename):
0047             print("WARNING: file ",filename," does not exist. Returning null.")
0048             return -1
0049
0050         # reading file
0051         lstime={}
0052         in_file = open(filename)
0053         pieces = paragraphs(in_file,start_of_new_beamspot_measurement,True)
0054         for piece in pieces:
0055             if len(piece) < 1:
0056                 continue
0057             try:
0058                 tmp = BeamspotMeasurement(piece)
0059             except Exception as err:
0060                 print("    ERROR Found corrupt " \
0061                       "beamspot measurement entry!", file=sys.stderr)
0062                 print("    !!! %s !!!" % str(err), file=sys.stderr)
0063                 continue
0064             # Argh!
0065             runfromfile=tmp.run_number
0066             (lumimin,lumimax)=tmp.lumi_range
0067             time_begin=tmp.time_begin
0068             time_end=tmp.time_end
0069             time_begin=calendar.timegm(time_begin.timetuple())
0070             time_end=calendar.timegm(time_end.timetuple())-23 # assume end of lumisection
0071             lstime[lumimin]=time_begin
0072             lstime[lumimax]=time_end
0073
0074         # order lumisections and make a list
0075         lslist=sorted(lstime.keys())
0076         lstimesorted=[]
0077         for ls in lslist:
0078             lstimesorted.append((ls,lstime[ls]))
0079         runlstime[run]=lstimesorted
0080
0081         #            print runfromfile
0082         #            print lumirange
0083         #            print time_begin, calendar.timegm(time_begin.timetuple())
0084         #            print time_end, calendar.timegm(time_end.timetuple())
0085
0086         in_file.close()
0087     # now give a time
0088     dcloselumi=999999
0089     closelumi=-1
0090     closetime=-1
0091     lstimesorted=runlstime[run]
0092
0093     for pair in lstimesorted:
0094         (lumi,time)=pair
0095         if abs(lumisection-lumi)<dcloselumi:
0096             dcloselumi=abs(lumisection-lumi)
0097             closelumi=lumi
0098             closetime=time
0099     if closelumi!=-1:
0100         finaltime=closetime+(lumisection-closelumi)*23
0101     else:
0102         finaltime=-1
0103
0104     return finaltime
0105
0106
0107 def readroot():
0108     rls=[]
0109     bxlist=[]
0110     allmeas={}
0111
0112     DIRES=['X0','Y0','Z0','width_X0','Width_Y0','Sigma_Z0','dxdz','dydz']
0113     # DIRES=['X0']
0114     rootfile="BxAnalysis_Fill_"+FILL+".root"
0115     filein=TFile(rootfile)
0116     for dire in DIRES:
0117         filein.cd(dire)
0118         # get list of histograms
0119         histolist=gDirectory.GetListOfKeys()
0120         iter = histolist.MakeIterator()
0121         key = iter.Next()
0122         while key:
0123             if key.GetClassName() == 'TH1F':
0124                 td = key.ReadObj()
0125                 histoname = td.GetName()
0126                 if "bx" in histoname:
0127 #                    print histoname
0128                     bx=histoname.split('_')[-1]
0129                     if bx not in bxlist:
0130 # this is to be removed
0131 #                        if len(bxlist)>=2:
0132 #                            key = iter.Next()
0133 #                            continue
0134 # end to be removed
0135                         bxlist.append(bx)
0136                         allmeas[bx]={}
0137 #                    print bx,histoname
0138                     histo=gDirectory.Get(histoname)
0139                     nbin=histo.GetNbinsX()
0140
0141                     thisbx=allmeas[bx]
0142
0143                     for bin in range(1,nbin+1):
0144                         label=histo.GetXaxis().GetBinLabel(bin)
0145                         label=label.strip()
0146                         if ":" not in label:
0147                             # not a valid label of type run:lumi-lumi, skip it
0148                             continue
0149
0150                         cont=histo.GetBinContent(bin)
0151                         if cont!=cont:
0152                             # it's a nan
0153                             cont=-999.0
0154                         err=histo.GetBinError(bin)
0155                         if err!=err:
0156                             err=-999.0
0157                         #                        if len(bxlist)==1:
0158                         #                            rls.append(label)
0159                         #                            print label
0160                         #                        else:
0161                         if label not in rls:
0162                             print("New range:",label," found in ",histoname)
0163                             rls.append(label)
0164
0165                         if label in thisbx.keys():
0166                             thismeas=thisbx[label]
0167                         else:
0168                             thisbx[label]=bsmeas()
0169                             thismeas=thisbx[label]
0170                         #  now filling up
0171                         if dire=='X0':
0172                             thismeas.x=cont
0173                             thismeas.ex=err
0174                         if dire=='Y0':
0175                             thismeas.y=cont
0176                             thismeas.ey=cont
0177                         if dire=='Z0':
0178                             thismeas.z=cont
0179                             thismeas.ez=err
0180                         if dire=='width_X0':
0181                             thismeas.wx=cont
0182                             thismeas.ewx=err
0183                         if dire=='Width_Y0':
0184                             thismeas.wy=cont
0185                             thismeas.ewy=err
0186                         if dire=='Sigma_Z0':
0187                             thismeas.wz=cont
0188                             thismeas.ewz=err
0189                         if dire=='dxdz':
0190                             thismeas.dxdz=cont
0191                             thismeas.edxdz=err
0192                         if dire=='dydz':
0193                             thismeas.dydz=cont
0194                             thismeas.edydz=err
0195
0196
0197             key = iter.Next()
0198
0199
0200     #    for name in pippo:
0201     #        print name
0202
0203     filein.Close()
0204
0205     # let's try to show it
0206 #    for bx in allmeas.keys():
0207 #        print "bx=",bx
0208 #        bxmeas=allmeas[bx]
0209 #        for meas in bxmeas.keys():
0210 #            print "meas time=",meas
0211 #            thismeas=bxmeas[meas]
0212 #            print thismeas.x,thismeas.ex,thismeas.y,thismeas.ey,thismeas.z,thismeas.ez
0213 #            print thismeas.wx,thismeas.ewx,thismeas.wy,thismeas.ewy,thismeas.wz,thismeas.ewz
0214 #            print thismeas.dxdz,thismeas.edxdz,thismeas.dydz,thismeas.edydz
0215     return allmeas
0216
0217 if __name__ == '__main__':
0218     if len(sys.argv)!=2:
0219         print("Usage: :",sys.argv[0]," <fillnr>")
0220         sys.exit(1)
0221     FILL=sys.argv[1]
0222
0223     allmeas=readroot()
0224     # now write it
0225
0226     for bx in allmeas.keys():
0227         print("writing bx=",bx)
0228         bxmeas=allmeas[bx]
0229         lines={}
0230         for meas in bxmeas.keys():
0231             # first derive time in unix time
0232             runno=meas.split(':')[0]
0233             runno=runno.strip()
0234             lumirange=meas.split(':')[1]
0235             lumimin=lumirange.split('-')[0]
0236             lumimin=lumimin.strip()
0237             lumimax=lumirange.split('-')[1]
0238             lumimax=lumimax.strip()
0239             lumimid=int((int(lumimin)+int(lumimax))/2.)
0240             meastime=timeof(runno,lumimid)
0241             print(runno,str(lumimid),meastime)
0242
0243             thismeas=bxmeas[meas]
0244 #            print thismeas.x,thismeas.ex,thismeas.y,thismeas.ey,thismeas.z,thismeas.ez
0245 #            print thismeas.wx,thismeas.ewx,thismeas.wy,thismeas.ewy,thismeas.wz,thismeas.ewz
0246 #            print thismeas.dxdz,thismeas.edxdz,thismeas.dydz,thismeas.edydz
0247             line=str(meastime)+" "
0248             line+="1 "
0249             line+="%11.7f %11.7f %11.7f %11.7f %11.7f %11.7f " % (-thismeas.x*10,thismeas.ex*10,
0250                                                             thismeas.y*10,thismeas.ey*10,
0251                                                             -thismeas.z*10,thismeas.ez*10)
0252             line+="%11.7f %11.7f %11.7f %11.7f %11.7f %11.7f " % (thismeas.wx*10,thismeas.ewx*10,
0253                                                             thismeas.wy*10,thismeas.ewy*10,
0254                                                             thismeas.wz*10,thismeas.ewz*10)
0255             line+="%11.7f %11.7f %11.7f %11.7f" % (thismeas.dxdz,thismeas.edxdz,-thismeas.dydz,thismeas.edydz)
0256             line+="\n"
0257
0258             # validate it
0259             if (thismeas.x != 0.0 and thismeas.y != 0.0 and thismeas.z != 0.0 and
0260                 thismeas.wx != 0.0 and thismeas.wy != 0.0 and thismeas.wz != 0.0 and
0261                 thismeas.dxdz != 0.0 and thismeas.dydz != 0.0 ):
0262                 lines[meastime]=line
0263
0264
0265         # now write it
0266         WORKDIR=OUTDIR+FILL
0267         os.system("mkdir -p "+WORKDIR)
0268         rfbucket=(int(bx)-1)*10+1
0269         filename=WORKDIR+"/"+FILL+"_lumireg_"+str(rfbucket)+"_CMS.txt"
0270         file=open(filename,'w')
0271         sortedtimes=sorted(lines.keys())
0272         for meastime in sortedtimes:
0273             file.write(lines[meastime])
0274         file.close()
0275
0276     #    print timeof("168264",25)