Reference Manual

Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 ########################################################################
00003 # Command to produce perbunch and specific lumi                        #
00004 #                                                                      #
00005 # Author:      Zhen Xie                                                #
00006 ########################################################################
00007 #
00008 # dump all fills into files.
00009 # allfills.txt all the existing fills.
00010 # fill_num.txt all the runs in the fill
00011 # dumpFill -o outputdir
00012 # dumpFill -f fillnum generate runlist for the given fill
00013 #
00014 import os,os.path,sys,math,array,datetime,time,calendar,re
00015 import coral
00016 
00017 from RecoLuminosity.LumiDB import argparse,sessionManager,lumiTime,CommonUtil,lumiCalcAPI,lumiParameters,revisionDML,normDML
00018 MINFILL=1800
00019 MAXFILL=9999
00020 allfillname='allfills.txt'
00021 
00022 def getFillFromDB(schema,fillnum):
00023     '''
00024     output: {run:starttime}
00025     '''
00026     runtimesInFill={}
00027     fillrundict=lumiCalcAPI.fillrunMap(schema,fillnum)
00028     if len(fillrundict)>0:
00029         runs=fillrundict.values()[0]
00030         runlsdict=dict(zip(runs,[None]*len(runs)))
00031         runresult=lumiCalcAPI.runsummary(schema,runlsdict)    
00032         for perrundata in runresult:
00033             runtimesInFill[perrundata[0]]=perrundata[7]
00034     return runtimesInFill
00035 
00036 def listfilldir(indir):
00037     '''
00038     list all fills contained in the given dir
00039     input: indir
00040     output: [fill]
00041     '''
00042     fillnamepat=r'^[0-9]{4}$'
00043     p=re.compile(fillnamepat)
00044     processedfills=[]
00045     dirList=os.listdir(indir)
00046     for fname in dirList:
00047         if p.match(fname) and os.path.isdir(os.path.join(indir,fname)):#found fill dir
00048             allfs=os.listdir(os.path.join(indir,fname))
00049             for myfile in allfs:
00050                 sumfilenamepat=r'^[0-9]{4}_bxsum_CMS.txt$'
00051                 s=re.compile(sumfilenamepat)
00052                 if s.match(myfile):
00053                     #only if fill_summary_CMS.txt file exists
00054                     processedfills.append(int(fname))
00055     return processedfills
00056 
00057 def lastcompleteFill(infile):
00058     '''
00059     parse infile to find LASTCOMPLETEFILL
00060     input: input file name
00061     output: last completed fill number
00062     '''
00063     lastfill=None
00064     hlinepat=r'(LASTCOMPLETEFILL )([0-9]{4})'
00065     h=re.compile(hlinepat)
00066     dqmfile=open(infile,'r')
00067     for line in dqmfile:
00068         result=h.match(line)
00069         if result:
00070             lastfill=result.group(2)
00071             break
00072     return int(lastfill)
00073 
00074 def calculateSpecificLumi(lumi,lumierr,beam1intensity,beam1intensityerr,beam2intensity,beam2intensityerr):
00075     '''
00076     calculate specific lumi
00077     input: instlumi, instlumierror,beam1intensity,beam1intensityerror,beam2intensity,beam2intensityerror
00078     output (specific lumi value,specific lumi error)
00079     '''
00080     specificlumi=0.0
00081     specificlumierr=0.0
00082     if beam1intensity<0: beam1intensity=0
00083     if beam2intensity<0: beam2intensity=0
00084     if beam1intensity>0.0 and  beam2intensity>0.0:
00085         specificlumi=float(lumi)/(float(beam1intensity)*float(beam2intensity))
00086         specificlumierr=specificlumi*math.sqrt(lumierr**2/lumi**2+beam1intensityerr**2/beam1intensity**2+beam2intensityerr**2/beam2intensity**2)
00087     return (specificlumi,specificlumierr)
00088 
00089 def getFillFromFile(fillnum,inputdir):
00090     '''
00091     parse fill_xxx.txt files in the input directory for runs, starttime in the fill
00092     input: fillnumber, input dir
00093     output: {run:tarttime}
00094     '''
00095     runtimesInFill={}
00096     #look for files 'fill_num.txt' in inputdir
00097     for filename in os.listdir(inputdir):
00098         mpat=r'^fill_[0-9]{4}.txt$'
00099         m=re.compile(mpat)
00100         if m.match(filename) is None:
00101             continue
00102         filename=filename.strip()
00103         if filename.find('.')==-1: continue            
00104         basename,extension=filename.split('.')        
00105         if not extension or extension!='txt':
00106             continue
00107         if basename.find('_')==-1: continue
00108         prefix,number=basename.split('_')
00109         if not number : continue
00110         if fillnum!=int(number):continue
00111         f=open(os.path.join(inputdir,'fill_'+number+'.txt'),'r')
00112         for line in f:
00113             l=line.strip()
00114             fields=l.split(',')
00115             if len(fields)<2 : continue
00116             runtimesInFill[int(fields[0])]=fields[1]
00117         f.close()
00118     return runtimesInFill
00119 
00120 #####output methods####
00121 def filltofiles(allfills,runsperfill,runtimes,dirname):
00122     '''
00123     write runnumber:starttime map per fill to files
00124     '''
00125     f=open(os.path.join(dirname,allfillname),'w')
00126     for fill in allfills:
00127         print >>f,'%d'%(fill)
00128     f.close()
00129     for fill,runs in runsperfill.items():
00130         filename='fill_'+str(fill)+'.txt'
00131         if len(runs)!=0:
00132             f=open(os.path.join(dirname,filename),'w')
00133             for run in runs:
00134                 print >>f,'%d,%s'%(run,runtimes[run])
00135             f.close()
00136             
00137 def specificlumiTofile(fillnum,filldata,outdir):
00138     #
00139     #input : fillnum
00140     #        filldata: {bxidx:[[lstime,beamstatusfrac,lumivalue,lumierror,speclumi,speclumierr]],[]}
00141     #sorted by bxidx, sorted by lstime inside list
00142     #check outdir/fillnum subdir exists; if not, create it; else outdir=outdir/fillnum
00143     #
00144     if not filldata:
00145         print 'empty input data, do nothing for fill ',fillnum
00146         return
00147     timedict={}#{lstime:[[stablebeamfrac,lumi,lumierr,speclumi,speclumierr]]}
00148     filloutdir=os.path.join(outdir,str(fillnum))
00149     if not os.path.exists(filloutdir):
00150         os.mkdir(filloutdir)
00151     for cmsbxidx,perbxdata in filldata.items():
00152         lhcbucket=0
00153         if cmsbxidx!=0:
00154             lhcbucket=(cmsbxidx-1)*10+1
00155         a=sorted(perbxdata,key=lambda x:x[0])
00156         filename=str(fillnum)+'_lumi_'+str(lhcbucket)+'_CMS.txt'
00157         linedata=[]
00158         for perlsdata in a:
00159             ts=int(perlsdata[0])
00160             beamstatusfrac=perlsdata[1]
00161             lumi=perlsdata[2]
00162             lumierror=perlsdata[3]
00163             #beam1intensity=perlsdata[4]
00164             #beam2intensity=perlsdata[5]
00165             speclumi=perlsdata[4]
00166             speclumierror= perlsdata[5]
00167             if lumi>0:
00168                 linedata.append([ts,beamstatusfrac,lumi,lumierror,speclumi,speclumierror])
00169             if not timedict.has_key(ts):
00170                 timedict[ts]=[]
00171             timedict[ts].append([beamstatusfrac,lumi,lumierror,speclumi,speclumierror])
00172         if len(linedata)>10:#at least 10 good ls
00173             f=open(os.path.join(filloutdir,filename),'w')
00174             for line in linedata:
00175                 print >>f, '%d\t%e\t%e\t%e\t%e\t%e'%(line[0],line[1],line[2],line[3],line[4],line[5])
00176             f.close()
00177     #print 'writing avg file'
00178     summaryfilename=str(fillnum)+'_lumi_CMS.txt'
00179     f=None
00180     lstimes=timedict.keys()
00181     lstimes.sort()
00182     fillseg=[]
00183     lscounter=0
00184     for lstime in lstimes:
00185         allvalues=timedict[lstime]
00186         transposedvalues=CommonUtil.transposed(allvalues,0.0)
00187         bstatfrac=transposedvalues[0][0]#beamstatus does not change with bx position
00188         lumivals=transposedvalues[1]
00189         lumitot=sum(lumivals)
00190         if bstatfrac==1.0 :
00191             fillseg.append([lstime,lumitot])
00192         lumierrs=transposedvalues[2]
00193         lumierrortot=math.sqrt(sum(map(lambda x:x**2,lumierrs)))
00194         specificvals=transposedvalues[3]
00195         specificavg=sum(specificvals)/float(len(specificvals))#avg spec lumi
00196         specificerrs=transposedvalues[4]
00197         specifictoterr=math.sqrt(sum(map(lambda x:x**2,specificerrs)))
00198         specificerravg=specifictoterr/float(len(specificvals))
00199         if lscounter==0:
00200             f=open(os.path.join(filloutdir,summaryfilename),'w')
00201         lscounter+=1
00202         print >>f,'%d\t%e\t%e\t%e\t%e\t%e'%(lstime,bstatfrac,lumitot,lumierrortot,specificavg,specificerravg)
00203     if f is not None:
00204         f.close()
00205     #print 'writing summary file'
00206     fillsummaryfilename=str(fillnum)+'_bxsum_CMS.txt'
00207     f=open(os.path.join(filloutdir,fillsummaryfilename),'w')    
00208     if len(fillseg)==0:
00209         print >>f,'%s'%('#no stable beams')
00210         f.close()
00211         return
00212     previoustime=fillseg[0][0]
00213     boundarytime=fillseg[0][0]
00214     #print 'boundary time ',boundarytime
00215     summaryls={}
00216     summaryls[boundarytime]=[]
00217     for [lstime,lumitot] in fillseg:#fillseg is everything with stable beam flag
00218         if lstime-previoustime>50.0:
00219             boundarytime=lstime
00220             #print 'found new boundary ',boundarytime
00221             summaryls[boundarytime]=[]
00222      #   print 'appending ',boundarytime,lstime,lumitot
00223         summaryls[boundarytime].append([lstime,lumitot])
00224         previoustime=lstime
00225     #print summaryls
00226    
00227     summarylstimes=summaryls.keys()
00228     summarylstimes.sort()
00229     lumip=lumiParameters.ParametersObject()
00230     for bts in summarylstimes:
00231         startts=bts
00232         tsdatainseg=summaryls[bts]
00233         #print 'tsdatainseg ',tsdatainseg
00234         stopts=tsdatainseg[-1][0]
00235         plu=max(CommonUtil.transposed(tsdatainseg,0.0)[1])
00236         lui=sum(CommonUtil.transposed(tsdatainseg,0.0)[1])*lumip.lslengthsec()
00237         print >>f,'%d\t%d\t%e\t%e'%(startts,stopts,plu,lui)
00238     f.close()
00239 
00240 def getSpecificLumi(schema,fillnum,inputdir,dataidmap,normmap,xingMinLum=0.0,amodetag='PROTPHYS',bxAlgo='OCC1'):
00241     '''
00242     specific lumi in 1e-30 (ub-1s-1) unit
00243     lumidetail occlumi in 1e-27
00244     1309_lumi_401_CMS.txt
00245     time(in seconds since January 1,2011,00:00:00 UTC) stab(fraction of time spent in stable beams for this time bin) l(lumi in Hz/ub) dl(point-to-point error on lumi in Hz/ub) sl(specific lumi in Hz/ub) dsl(error on specific lumi)
00246     20800119.0 1 -0.889948 0.00475996848729 0.249009 0.005583287562 -0.68359 6.24140208607 0.0 0.0 0.0 0.0 0.0 0.0 0.0383576 0.00430892097862 0.0479095 0.00430892097862 66.6447 4.41269758764 0.0 0.0 0.0
00247     result [(time,beamstatusfrac,lumi,lumierror,speclumi,speclumierror)]
00248     '''
00249     t=lumiTime.lumiTime()
00250     fillbypos={}#{bxidx:[[ts,beamstatusfrac,lumi,lumierror,spec1,specerror],[]]}
00251     runtimesInFill=getFillFromDB(schema,fillnum)#{runnum:starttimestr}
00252     runlist=runtimesInFill.keys()
00253     if not runlist: return fillbypos
00254     irunlsdict=dict(zip(runlist,[None]*len(runlist)))
00255     #prirunlsdict
00256     GrunsummaryData=lumiCalcAPI.runsummaryMap(session.nominalSchema(),irunlsdict)
00257     lumidetails=lumiCalcAPI.deliveredLumiForIds(schema,irunlsdict,dataidmap,GrunsummaryData,beamstatusfilter=None,normmap=normmap,withBXInfo=True,bxAlgo=bxAlgo,xingMinLum=xingMinLum,withBeamIntensity=True,lumitype='HF')
00258     #
00259     #output: {run:[lumilsnum(0),cmslsnum(1),timestamp(2),beamstatus(3),beamenergy(4),deliveredlumi(5),calibratedlumierr(6),(bxvalues,bxerrs)(7),(bxidx,b1intensities,b2intensities)(8),fillnum(9)]}
00260     #
00261     totalstablebeamls=0
00262     orderedrunlist=sorted(lumidetails)
00263     for run in orderedrunlist:
00264         perrundata=lumidetails[run]
00265         for perlsdata in perrundata:
00266             beamstatus=perlsdata[3]
00267             if beamstatus=='STABLE BEAMS':
00268                 totalstablebeamls+=1
00269     #print 'totalstablebeamls in fill ',totalstablebeamls
00270     if totalstablebeamls<10:#less than 10 LS in a fill has 'stable beam', it's no a good fill
00271         print 'fill ',fillnum,' , having less than 10 stable beam lS, is not good, skip'
00272         return fillbypos
00273     lumiparam=lumiParameters.ParametersObject()
00274     for run in orderedrunlist:
00275         perrundata=lumidetails[run]
00276         for perlsdata in perrundata:
00277             beamstatusfrac=0.0
00278             tsdatetime=perlsdata[2]
00279             ts=calendar.timegm(tsdatetime.utctimetuple())
00280             beamstatus=perlsdata[3]
00281             if beamstatus=='STABLE BEAMS':
00282                 beamstatusfrac=1.0
00283             (bxidxlist,bxvaluelist,bxerrolist)=perlsdata[7]
00284             #instbxvaluelist=[x/lumiparam.lslengthsec() for x in bxvaluelist if x]
00285             instbxvaluelist=[x for x in bxvaluelist if x]
00286             maxlumi=0.0
00287             if len(instbxvaluelist)!=0:
00288                 maxlumi=max(instbxvaluelist)
00289             avginstlumi=0.0
00290             if len(instbxvaluelist)!=0:
00291                 avginstlumi=sum(instbxvaluelist)
00292             (intbxidxlist,b1intensities,b2intensities)=perlsdata[8]#contains only non-zero bx
00293             for bxidx in bxidxlist:
00294                 idx=bxidxlist.index(bxidx)
00295                 instbxvalue=bxvaluelist[idx]
00296                 bxerror=bxerrolist[idx]
00297                 if instbxvalue<max(xingMinLum,maxlumi*0.2):
00298                     continue
00299                 bintensityPos=-1
00300                 try:
00301                     bintensityPos=intbxidxlist.index(bxidx)
00302                 except ValueError:
00303                     pass
00304                 if bintensityPos<=0:
00305                     fillbypos.setdefault(bxidx,[]).append([ts,beamstatusfrac,instbxvalue,bxerror,0.0,0.0])
00306                     continue
00307                 b1intensity=b1intensities[bintensityPos]
00308                 b2intensity=b2intensities[bintensityPos]
00309                 speclumi=calculateSpecificLumi(instbxvalue,bxerror,b1intensity,0.0,b2intensity,0.0)
00310                 fillbypos.setdefault(bxidx,[]).append([ts,beamstatusfrac,instbxvalue,bxerror,speclumi[0],speclumi[1]])
00311     return fillbypos
00312 
00313 
00314 ##############################
00315 ## ######################## ##
00316 ## ## ################## ## ##
00317 ## ## ## Main Program ## ## ##
00318 ## ## ################## ## ##
00319 ## ######################## ##
00320 ##############################
00321         
00322 if __name__ == '__main__':
00323     parser = argparse.ArgumentParser(prog=os.path.basename(sys.argv[0]),description = "specific lumi",formatter_class=argparse.ArgumentDefaultsHelpFormatter)
00324     amodetagChoices = [ "PROTPHYS","IONPHYS",'PAPHYS' ]
00325     xingAlgoChoices =[ "OCC1","OCC2","ET"]
00326     # parse arguments
00327     parser.add_argument('-c',dest='connect',
00328                         action='store',
00329                         required=False,
00330                         help='connect string to lumiDB,optional',
00331                         default='frontier://LumiCalc/CMS_LUMI_PROD')
00332     parser.add_argument('-P',dest='authpath',
00333                         action='store',
00334                         help='path to authentication file,optional')
00335     parser.add_argument('-i',dest='inputdir',
00336                         action='store',
00337                         required=False,
00338                         help='output dir',
00339                         default='.')
00340     parser.add_argument('-o',dest='outputdir',
00341                         action='store',
00342                         required=False,
00343                         help='output dir',
00344                         default='.')
00345     parser.add_argument('-f','--fill',dest='fillnum',
00346                         action='store',
00347                         required=False,
00348                         help='specific fill',
00349                         default=None)
00350     parser.add_argument('--minfill',dest='minfill',
00351                         type=int,
00352                         action='store',
00353                         required=False,
00354                         default=MINFILL,
00355                         help='min fill')
00356     parser.add_argument('--maxfill',dest='maxfill',
00357                         type=int,
00358                         action='store',
00359                         required=False,
00360                         default=MAXFILL,
00361                         help='maximum fillnumber '
00362                         )
00363     parser.add_argument('--amodetag',dest='amodetag',
00364                         action='store',
00365                         choices=amodetagChoices,
00366                         required=False,
00367                         help='specific accelerator mode choices [PROTOPHYS,IONPHYS,PAPHYS] (optional)')
00368     parser.add_argument('--xingMinLum', dest = 'xingMinLum',
00369                         type=float,
00370                         default=1e-03,
00371                         required=False,
00372                         help='Minimum luminosity considered for lumibylsXing action')
00373     parser.add_argument('--xingAlgo', dest = 'bxAlgo',
00374                         default='OCC1',
00375                         required=False,
00376                         help='algorithm name for per-bunch lumi ')
00377     parser.add_argument('--normtag',dest='normtag',action='store',
00378                         required=False,
00379                         help='norm tag',
00380                         default=None)
00381     parser.add_argument('--datatag',dest='datatag',action='store',
00382                         required=False,
00383                         help='data tag',
00384                         default=None)
00385     #
00386     #command configuration 
00387     #
00388     parser.add_argument('--siteconfpath',dest='siteconfpath',action='store',
00389                         help='specific path to site-local-config.xml file, optional. If path undefined, fallback to cern proxy&server')
00390     #
00391     #switches
00392     #
00393     parser.add_argument('--without-correction',dest='withoutNorm',action='store_true',
00394                         help='without any correction/calibration' )
00395     parser.add_argument('--debug',dest='debug',action='store_true',
00396                         help='debug')
00397     options=parser.parse_args()
00398     if options.authpath:
00399         os.environ['CORAL_AUTH_PATH'] = options.authpath
00400     ##
00401     #query DB for all fills and compare with allfills.txt
00402     #if found newer fills, store  in mem fill number
00403     #reprocess anyway the last 1 fill in the dir
00404     #redo specific lumi for all marked fills
00405     ##    
00406     svc=sessionManager.sessionManager(options.connect,authpath=options.authpath,debugON=options.debug)
00407     session=svc.openSession(isReadOnly=True,cpp2sqltype=[('unsigned int','NUMBER(10)'),('unsigned long long','NUMBER(20)')])
00408 
00409     fillstoprocess=[]
00410     maxfillnum=options.maxfill
00411     minfillnum=options.minfill
00412     if options.fillnum is not None: #if process a specific single fill
00413         fillstoprocess.append(int(options.fillnum))
00414     else:
00415         session.transaction().start(True)
00416         schema=session.nominalSchema()
00417         allfillsFromDB=lumiCalcAPI.fillInRange(schema,fillmin=minfillnum,fillmax=maxfillnum,amodetag=options.amodetag)
00418         processedfills=listfilldir(options.outputdir)
00419         lastcompletedFill=lastcompleteFill(os.path.join(options.inputdir,'runtofill_dqm.txt'))
00420         for pf in processedfills:
00421             if pf>lastcompletedFill:
00422                 print '\tremove unfinished fill from processed list ',pf
00423                 processedfills.remove(pf)
00424         for fill in allfillsFromDB:
00425             if fill not in processedfills :
00426                 if int(fill)<=lastcompletedFill:
00427                     if int(fill)>minfillnum and int(fill)<maxfillnum:
00428                         fillstoprocess.append(fill)
00429                 else:
00430                     print 'ongoing fill...',fill
00431         session.transaction().commit()
00432     print 'fills to process : ',fillstoprocess
00433     if len(fillstoprocess)==0:
00434         print 'no fill to process, exit '
00435         exit(0)
00436 
00437     
00438     print '===== Start Processing Fills',fillstoprocess
00439     print '====='
00440     filldata={}
00441     #
00442     # check datatag
00443     #
00444     reqfillmin=min(fillstoprocess)
00445     reqfillmax=max(fillstoprocess)
00446     session.transaction().start(True)
00447     runlist=lumiCalcAPI.runList(session.nominalSchema(),options.fillnum,runmin=None,runmax=None,fillmin=reqfillmin,fillmax=reqfillmax,startT=None,stopT=None,l1keyPattern=None,hltkeyPattern=None,amodetag=options.amodetag,nominalEnergy=None,energyFlut=None,requiretrg=False,requirehlt=False)
00448     
00449     datatagname=options.datatag
00450     if not datatagname:
00451         (datatagid,datatagname)=revisionDML.currentDataTag(session.nominalSchema())
00452         dataidmap=revisionDML.dataIdsByTagId(session.nominalSchema(),datatagid,runlist=runlist,withcomment=False)
00453         #{run:(lumidataid,trgdataid,hltdataid,())}
00454     else:
00455         dataidmap=revisionDML.dataIdsByTagName(session.nominalSchema(),datatagname,runlist=runlist,withcomment=False)
00456 
00457     #
00458     # check normtag and get norm values if required
00459     #
00460     normname='NONE'
00461     normid=0
00462     normvalueDict={}
00463     if not options.withoutNorm:
00464         normname=options.normtag
00465         if not normname:
00466             normmap=normDML.normIdByType(session.nominalSchema(),lumitype='HF',defaultonly=True)
00467             if len(normmap):
00468                 normname=normmap.keys()[0]
00469                 normid=normmap[normname]
00470         else:
00471             normid=normDML.normIdByName(session.nominalSchema(),normname)
00472         if not normid:
00473             raise RuntimeError('[ERROR] cannot resolve norm/correction')
00474             sys.exit(-1)
00475         normvalueDict=normDML.normValueById(session.nominalSchema(),normid) #{since:[corrector(0),{paramname:paramvalue}(1),amodetag(2),egev(3),comment(4)]}
00476     session.transaction().commit()
00477     for fillnum in fillstoprocess:# process per fill
00478         session.transaction().start(True)
00479         filldata=getSpecificLumi(session.nominalSchema(),fillnum,options.inputdir,dataidmap,normvalueDict,xingMinLum=options.xingMinLum,amodetag=options.amodetag,bxAlgo=options.bxAlgo)
00480         specificlumiTofile(fillnum,filldata,options.outputdir)
00481         session.transaction().commit()
00482 
00483
/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/RecoLuminosity/LumiDB/scripts/specificLumi.py
