Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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