Ecal/python/EcalCondTools.py

0001 from __future__ import print_function
0002 from __future__ import absolute_import
0003 #
0004 # Misc functions to manipulate Ecal records
0005 # author: Stefano Argiro
0006 # id: $Id: EcalCondTools.py,v 1.14 2013/05/06 08:29:44 argiro Exp $
0007 #
0008 #
0009 # WARNING: we assume that the list of iovs for a given tag
0010 #          contains one element only, the case of several elements
0011 #          will need to be addressed
0012
0013 #from pluginCondDBPyInterface import *
0014 from CondCore.Utilities import iovInspector as inspect
0015 from ROOT import TCanvas,TH1F, TH2F, gStyle, TChain, TTree, TLegend, TFile
0016 from . import EcalPyUtils
0017 import sys
0018 from math import sqrt
0019
0020 def listTags(db):
0021     '''List all available tags for a given db '''
0022     try:
0023         db.startReadOnlyTransaction()
0024         tags = db.allTags()
0025         db.commitTransaction()
0026         for tag in tags.split():
0027             print(tag)
0028     except:
0029         return ""
0030
0031 def listIovs(db,tag):
0032     '''List all available iovs for a given tag'''
0033
0034     try :
0035        iov = inspect.Iov(db,tag)
0036        iovlist = iov.list()
0037        print("Available iovs for tag: ",tag)
0038        for p in iovlist:
0039            print("  Since " , p[1], " Till " , p[2])
0040
0041     except Exception as er :
0042         print(" listIovs exception ",er)
0043
0044 def dumpXML(db,tag,since,filename='dump.xml'):
0045     '''Dump record in XML format for a given tag '''
0046     try :
0047        iov = inspect.Iov(db,tag)
0048        db.startReadOnlyTransaction()
0049        Plug = __import__(str(db.payloadModules(tag)[0]))
0050        payload = Plug.Object(db)
0051        listOfIovElem= [iovElem for iovElem in db.iov(tag).elements]
0052        inst = 0
0053        for elem in db.iov(tag).elements :
0054            inst = inst + 1
0055            if str(elem.since())==str(since):
0056                found = inst
0057                break
0058        db.commitTransaction()
0059        payload = inspect.PayLoad(db, tag, elem)
0060        out = open(filename,'w')
0061        print(payload, file=out)
0062
0063     except Exception as er :
0064         print(" dumpXML exception ",er)
0065
0066 def plot (db, tag,since,filename='plot.root'):
0067     '''Invoke the plot function from the wrapper and save to the specified \
0068        file. The file format will reflect the extension given.'''
0069
0070     try :
0071        iov = inspect.Iov(db,tag)
0072        db.startReadOnlyTransaction()
0073        Plug = __import__(str(db.payloadModules(tag)[0]))
0074        payload = Plug.Object(db)
0075        listOfIovElem= [iovElem for iovElem in db.iov(tag).elements]
0076        inst = 0
0077        for elem in db.iov(tag).elements :
0078            inst = inst + 1
0079            if str(elem.since())==str(since):
0080                found = inst
0081                break
0082        db.commitTransaction()
0083        payload = inspect.PayLoad(db, tag, elem)
0084        payload.plot(filename,"",[],[])
0085
0086     except Exception as er :
0087         print(" plot exception ",er)
0088
0089
0090 def compare(tag1,db1,since1,
0091             tag2,db2,since2,filename='compare.root'):
0092   '''Produce comparison plots for two records. Save plots to file \
0093      according to format. tag can be an xml file'''
0094   print("EcalCondTools.py compare tag1 ", tag1, "since1 : ", since1, " tag2 ", tag2," s2 ", since2)
0095
0096   coeff_1_b=[]
0097   coeff_2_b=[]
0098
0099   coeff_1_e=[]
0100   coeff_2_e=[]
0101
0102   if  tag1.find(".xml") < 0:
0103       found=0
0104       try:
0105         db1.startReadOnlyTransaction()
0106         Plug = __import__(str(db1.payloadModules(tag1)[0]))
0107         payload = Plug.Object(db1)
0108         listOfIovElem= [iovElem for iovElem in db1.iov(tag1).elements]
0109         what = {'how':'barrel'}
0110         w = inspect.setWhat(Plug.What(),what)
0111         exb = Plug.Extractor(w)
0112         what = {'how':'endcap'}
0113         w = inspect.setWhat(Plug.What(),what)
0114         exe = Plug.Extractor(w)
0115 #        p = getObject(db1,tag1,since1,ex)
0116 #        p.extract(ex)
0117         print(" before loop")
0118         for elem in db1.iov(tag1).elements :
0119           if str(elem.since())==str(since1):
0120             found=1
0121             payload.load(elem)
0122             payload.extract(exb)
0123             coeff_1_b = [i for i in exb.values()]# first set of coefficients
0124             payload.extract(exe)
0125             coeff_1_e = [i for i in exe.values()]
0126         db1.commitTransaction()
0127
0128       except Exception as er :
0129         print(" compare first set exception ",er)
0130       if not found :
0131         print("Could not retrieve payload for tag: " , tag1, " since: ", since1)
0132         sys.exit(0)
0133
0134   else:
0135       coeff_1_b,coeff_1_e = EcalPyUtils.fromXML(tag1)
0136
0137   if  tag2.find(".xml")<0:
0138       found=0
0139       try:
0140         db2.startReadOnlyTransaction()
0141         Plug = __import__(str(db2.payloadModules(tag2)[0]))
0142         what = {'how':'barrel'}
0143         w = inspect.setWhat(Plug.What(),what)
0144         exb = Plug.Extractor(w)
0145         what = {'how':'endcap'}
0146         w = inspect.setWhat(Plug.What(),what)
0147         exe = Plug.Extractor(w)
0148 #        p = getObject(db2,tag2,since2)
0149 #        p.extract(ex)
0150         for elem in db2.iov(tag2).elements :
0151           if str(elem.since())==str(since2):
0152             found=1
0153             payload.load(elem)
0154             payload.extract(exb)
0155             coeff_2_b = [i for i in exb.values()]# second set of coefficients
0156             payload.extract(exe)
0157             coeff_2_e = [i for i in exe.values()]
0158         db2.commitTransaction()
0159
0160       except Exception as er :
0161           print(" compare second set exception ",er)
0162       if not found :
0163         print("Could not retrieve payload for tag: " , tag2, " since: ", since2)
0164         sys.exit(0)
0165
0166   else:
0167       coeff_2_b,coeff_2_e = EcalPyUtils.fromXML(tag2)
0168
0169   gStyle.SetPalette(1)
0170
0171   savefile = TFile(filename,"RECREATE")
0172
0173   ebhisto,ebmap, profx, profy= compareBarrel(coeff_1_b,coeff_2_b)
0174   eephisto,eepmap,eemhisto,eemmap=compareEndcap(coeff_1_e,coeff_2_e)
0175
0176 #make more canvas (suppressed : use a root file)
0177
0178 #  cEBdiff = TCanvas("EBdiff","EBdiff")
0179 #  cEBdiff.Divide(2,2)
0180
0181 #  cEBdiff.cd(1)
0182 #  ebhisto.Draw()
0183 #  cEBdiff.cd(2)
0184 #  ebmap.Draw("colz")
0185 #  cEBdiff.cd(3)
0186 #  profx.Draw()
0187 #  cEBdiff.cd(4)
0188 #  profy.Draw()
0189
0190 #  EBfilename = "EB_"+filename
0191 #  cEBdiff.SaveAs(EBfilename)
0192
0193 #  cEEdiff = TCanvas("EEdiff","EEdiff")
0194 #  cEEdiff.Divide(2,2)
0195
0196 #  cEEdiff.cd(1)
0197 #  eephisto.Draw()
0198 #  cEEdiff.cd(2)
0199 #  eepmap.Draw("colz")
0200 #  cEEdiff.cd(3)
0201 #  eemhisto.Draw()
0202 #  cEEdiff.cd(4)
0203 #  eemmap.Draw("colz")
0204
0205 #  EEfilename = "EE_"+filename
0206 #  cEEdiff.SaveAs(EEfilename)
0207
0208
0209   eeborderphisto,eeborderpmap,eebordermhisto,eebordermmap=compareEndcapBorder(coeff_1_e,coeff_2_e)
0210   ebborderhisto,ebbordermap = compareBarrelBorder(coeff_1_b,coeff_2_b)
0211
0212   border_diff_distro_can = TCanvas("border_difference","borders difference")
0213   border_diff_distro_can.Divide(2,3)
0214
0215   border_diff_distro_can.cd(1)
0216   ebborderhisto.Draw()
0217   border_diff_distro_can.cd(2)
0218   ebbordermap.Draw("colz")
0219   border_diff_distro_can.cd(3)
0220   eeborderphisto.Draw()
0221   border_diff_distro_can.cd(4)
0222   eeborderpmap.Draw("colz")
0223   border_diff_distro_can.cd(5)
0224   eebordermhisto.Draw()
0225   border_diff_distro_can.cd(6)
0226   eebordermmap.Draw("colz")
0227
0228   bordersfilename = "borders_"+filename
0229   prof_filename = "profiles_"+filename
0230
0231 #  border_diff_distro_can.SaveAs(bordersfilename)
0232
0233   savefile.Write()
0234
0235
0236
0237 def histo (db, tag,since,filename='histo.root'):
0238     '''Make histograms and save to file. tag can be an xml file'''
0239
0240     coeff_barl=[]
0241     coeff_endc=[]
0242
0243     if  tag.find(".xml")< 0:
0244       found=0
0245       try:
0246 #          exec('import '+db.moduleName(tag)+' as Plug')
0247         db.startReadOnlyTransaction()
0248         Plug = __import__(str(db.payloadModules(tag)[0]))
0249         payload = Plug.Object(db)
0250         listOfIovElem= [iovElem for iovElem in db.iov(tag).elements]
0251         what = {'how':'barrel'}
0252         w = inspect.setWhat(Plug.What(),what)
0253         exb = Plug.Extractor(w)
0254         what = {'how':'endcap'}
0255         w = inspect.setWhat(Plug.What(),what)
0256         exe = Plug.Extractor(w)
0257         for elem in db.iov(tag).elements :
0258           if str(elem.since())==str(since):
0259             found=1
0260             payload.load(elem)
0261             payload.extract(exb)
0262             coeff_barl = [i for i in exb.values()]
0263             payload.extract(exe)
0264             coeff_endc = [i for i in exe.values()]
0265         db.commitTransaction()
0266
0267       except Exception as er :
0268           print(" histo exception ",er)
0269       if not found :
0270         print("Could not retrieve payload for tag: " , tag, " since: ", since)
0271         sys.exit(0)
0272
0273     else :
0274       coeff_barl,coeff_endc=EcalPyUtils.fromXML(tag)
0275
0276     savefile = TFile(filename,"RECREATE")
0277
0278     ebmap, ebeta, ebphi, eePmap, ebdist, eeMmap, prof_eePL, prof_eePR, prof_eeML, prof_eeMR, ebBorderdist = makedist(coeff_barl, coeff_endc)
0279
0280     gStyle.SetPalette(1)
0281
0282     c =  TCanvas("CCdist")
0283     c.Divide(2,2)
0284
0285     c.cd(1)
0286     ebmap.Draw("colz")
0287     c.cd(2)
0288     eePmap.Draw("colz")
0289     c.cd(3)
0290     ebdist.Draw()
0291     c.cd(4)
0292     eeMmap.Draw("colz")
0293
0294 #    c.SaveAs(filename)
0295
0296     prof_eb_eta = ebeta.ProfileX()
0297     prof_eb_phi = ebphi.ProfileX()
0298
0299     c2 = TCanvas("CCprofiles")
0300     c2.Divide(2,2)
0301
0302     leg = TLegend(0.1,0.7,0.48,0.9)
0303
0304     c2.cd(1)
0305     prof_eb_eta.Draw()
0306     c2.cd(2)
0307     prof_eb_phi.Draw()
0308     c2.cd(3)
0309     prof_eePL.SetMarkerStyle(8)
0310     prof_eePL.Draw()
0311     prof_eePR.SetMarkerStyle(8)
0312     prof_eePR.SetMarkerColor(2)
0313     prof_eePR.Draw("same")
0314     prof_eeML.SetMarkerStyle(8)
0315     prof_eeML.SetMarkerColor(3)
0316     prof_eeML.Draw("same")
0317     prof_eeMR.SetMarkerStyle(8)
0318     prof_eeMR.SetMarkerColor(5)
0319     prof_eeMR.Draw("same")
0320     leg.AddEntry(prof_eePL,"Dee1","lp")
0321     leg.AddEntry(prof_eePR,"Dee2","lp")
0322     leg.AddEntry(prof_eeMR,"Dee3","lp")
0323     leg.AddEntry(prof_eeML,"Dee4","lp")
0324     leg.Draw()
0325     c2.cd(4)
0326     ebBorderdist.Draw()
0327
0328     extrafilename = "profiles_"+filename
0329  #   c2.SaveAs(extrafilename)
0330
0331     savefile.Write()
0332
0333
0334 def getToken(db,tag,since):
0335     ''' Return payload token for a given iov, tag, db'''
0336     try :
0337        iov = inspect.Iov(db,tag)
0338        iovlist = iov.list()
0339        for p in iovlist:
0340            tmpsince=p[1]
0341            if str(tmpsince)==str(since) :
0342                return p[0]
0343        print("Could not retrieve token for tag: " , tag, " since: ", since)
0344        sys.exit(0)
0345
0346     except Exception as er :
0347        print(er)
0348
0349
0350 def getObject(db,tag,since):
0351     ''' Return payload object for a given iov, tag, db'''
0352     found=0
0353     try:
0354 #       exec('import '+db.moduleName(tag)+' as Plug')
0355        db.startReadOnlyTransaction()
0356        Plug = __import__(str(db.payloadModules(tag)[0]))
0357        print(" getObject Plug")
0358        payload = Plug.Object(db)
0359        db.commitTransaction()
0360        listOfIovElem= [iovElem for iovElem in db.iov(tag).elements]
0361        print(" getObject before loop")
0362        for elem in db.iov(tag).elements :
0363            if str(elem.since())==str(since):
0364                found=1
0365                print(" getObject found ", elem.since())
0366 #               return Plug.Object(elem)
0367                return elem
0368
0369     except Exception as er :
0370         print(" getObject exception ",er)
0371
0372     if not found :
0373         print("Could not retrieve payload for tag: " , tag, " since: ", since)
0374         sys.exit(0)
0375
0376
0377 def makedist(coeff_barl, coeff_endc) :
0378
0379     ebmap = TH2F("EB","EB",360,1,261,171, -85,86)
0380     eePmap = TH2F("EE","EE",100, 1,101,100,1,101)
0381     eeMmap = TH2F("EEminus","EEminus",100,1,101,100,1,101)
0382     ebdist = TH1F("EBdist","EBdist",100,-2,2)
0383     ebBorderdist = TH1F("EBBorderdist","EBBorderdist",100,-2,2)
0384
0385     ebeta = TH2F("ebeta","ebeta",171,-85,86,100,-2,2)
0386     ebphi = TH2F("ebphi","ebphi",360,1,361,100,-2,2)
0387
0388     eePL = TH2F("EEPL","EEPlus Left",50,10,55,100,-2,2)
0389     eePR = TH2F("EEPR","EEPlus Right",50,10,55,100,-2,2)
0390     eeML = TH2F("EEML","EEMinus Left",50,10,55,100,-2,2)
0391     eeMR = TH2F("EEMR","EEMinus Right",50,10,55,100,-2,2)
0392
0393     for i,c in enumerate(coeff_barl):
0394         ieta,iphi = EcalPyUtils.unhashEBIndex(i)
0395         ebmap.Fill(iphi,ieta,c)
0396         ebdist.Fill(c)
0397         ebeta.Fill(ieta,c)
0398         ebphi.Fill(iphi,c)
0399
0400         if (abs(ieta)==85 or abs(ieta)==65 or abs(ieta)==64 or abs(ieta)==45 or abs(ieta)==44 or abs(ieta)==25 or abs(ieta)==24 or abs(ieta)==1 or iphi%20==1 or iphi%20==0):
0401             ebBorderdist.Fill(c)
0402
0403
0404     for i,c in enumerate(coeff_endc):
0405         ix,iy,iz = EcalPyUtils.unhashEEIndex(i)
0406         R = sqrt((ix-50)*(ix-50)+(iy-50)*(iy-50))
0407
0408         if  iz>0:
0409             eePmap.Fill(ix,iy,c)
0410             if ix<50:
0411                 eePL.Fill(R,c,1)
0412             if ix>50:
0413                 eePR.Fill(R,c,1)
0414
0415         if iz<0:
0416             eeMmap.Fill(ix,iy,c)
0417             if ix<50:
0418                 eeML.Fill(R,c,1)
0419             if ix>50:
0420                 eeMR.Fill(R,c,1)
0421
0422     prof_eePL = eePL.ProfileX()
0423     prof_eePR = eePR.ProfileX()
0424     prof_eeML = eeML.ProfileX()
0425     prof_eeMR = eeMR.ProfileX()
0426
0427     return ebmap, ebeta, ebphi, eePmap, ebdist, eeMmap, prof_eePL, prof_eePR, prof_eeML, prof_eeMR, ebBorderdist
0428
0429 def compareBarrel(coeff_barl_1,coeff_barl_2) :
0430   '''Return an histogram and a map of the differences '''
0431
0432   diff_distro_h   = TH1F("diffh","diffh",100,-2,2)
0433   diff_distro_m   = TH2F("diffm","diffm",360,1,361,171,-85,86)
0434   diff_distro_m.SetStats(0)
0435   ebeta = TH2F("ebeta","ebeta",171,-85,86,100,-2,2)
0436   ebphi = TH2F("ebphi","ebphi",360,1,361,100,-2,2)
0437
0438
0439   for i,c in enumerate(coeff_barl_1):
0440       diff = c - coeff_barl_2[i]
0441       ieta,iphi= EcalPyUtils.unhashEBIndex(i)
0442       diff_distro_h.Fill(diff)
0443       diff_distro_m.Fill(iphi,ieta,diff)
0444       ebeta.Fill(ieta,diff)
0445       ebphi.Fill(iphi,diff)
0446
0447   prof_x_h = ebeta.ProfileX()
0448   prof_y_h = ebphi.ProfileX()
0449
0450   return diff_distro_h, diff_distro_m, prof_x_h, prof_y_h
0451
0452
0453
0454 def compareBarrelBorder(coeff_barl_1,coeff_barl_2) :
0455   '''Return an histogram and a map of the differences '''
0456
0457   diff_distro_border_h   = TH1F("diffborderh","diffh",100,-2,2)
0458   diff_distro_border_m   = TH2F("diffborderm","diffm",360,1,361,171,-85,86)
0459   diff_distro_border_m.SetStats(0)
0460
0461   for i,c in enumerate(coeff_barl_1):
0462       diff = c - coeff_barl_2[i]
0463       ieta,iphi= EcalPyUtils.unhashEBIndex(i)
0464       if (abs(ieta)==85 or abs(ieta)==65 or abs(ieta)==64 or abs(ieta)==45 or abs(ieta)==44 or abs(ieta)==25 or abs(ieta)==24 or abs(ieta)==1 or iphi%20==1 or iphi%20==0):
0465           diff_distro_border_h.Fill(diff)
0466       if (abs(ieta)==85 or abs(ieta)==65 or abs(ieta)==64 or abs(ieta)==45 or abs(ieta)==44 or abs(ieta)==25 or abs(ieta)==24 or abs(ieta)==1 or iphi%20==0 or iphi%20==1):
0467           diff_distro_border_m.Fill(iphi,ieta,diff)
0468
0469   return diff_distro_border_h, diff_distro_border_m
0470
0471
0472
0473
0474
0475 def compareEndcap(coeff_endc_1, coeff_endc_2) :
0476     ''' Return an histogram and a map of the differences for each endcap'''
0477
0478     diff_distro_h_eep   = TH1F("diff EE+","diff EE+",100,-2,2)
0479     diff_distro_h_eem   = TH1F("diff EE-","diff EE-",100,-2,2)
0480
0481
0482     diff_distro_m_eep   = TH2F("map EE+","map EE+",100,1,101,100,1,101)
0483     diff_distro_m_eem   = TH2F("map EE-","map EE-",100,1,101,100,1,101)
0484
0485     temp_h = TH1F("tempR","tempR",50,0,50)
0486
0487     diff_distro_m_eep.SetStats(0)
0488     diff_distro_m_eem.SetStats(0)
0489
0490
0491     for i,c in enumerate(coeff_endc_1):
0492       diff = c - coeff_endc_2[i]
0493       ix,iy,iz = EcalPyUtils.unhashEEIndex(i)
0494       R = sqrt((ix-50)*(ix-50)+(iy-50)*(iy-50))
0495
0496       if iz >0:
0497           diff_distro_h_eep.Fill(diff)
0498           diff_distro_m_eep.Fill(ix,iy,diff)
0499
0500       else:
0501           diff_distro_h_eem.Fill(diff)
0502           diff_distro_m_eem.Fill(ix,iy,diff)
0503
0504     return diff_distro_h_eep, \
0505            diff_distro_m_eep, \
0506            diff_distro_h_eem, \
0507            diff_distro_m_eem
0508
0509
0510
0511 def compareEndcapBorder(coeff_endc_1, coeff_endc_2) :
0512     ''' Return an histogram and a map of the differences for each endcap'''
0513
0514     border_diff_distro_h_eep   = TH1F("borderdiff EE+","diff EE+",100,-2,2)
0515     border_diff_distro_h_eem   = TH1F("borderdiff EE-","diff EE-",100,-2,2)
0516
0517
0518     border_diff_distro_m_eep   = TH2F("bordermap EE+","map EE+",100,1,101,100,1,101)
0519     border_diff_distro_m_eem   = TH2F("bordermap EE-","map EE-",100,1,101,100,1,101)
0520
0521     border_diff_distro_m_eep.SetStats(0)
0522     border_diff_distro_m_eem.SetStats(0)
0523
0524
0525     for i,c in enumerate(coeff_endc_1):
0526       diff = c - coeff_endc_2[i]
0527       ix,iy,iz = EcalPyUtils.unhashEEIndex(i)
0528       Rsq = ((ix-50.0)**2+(iy-50.0)**2)
0529
0530       if (iz >0 and (Rsq<144.0 or Rsq>2500.0)):
0531           border_diff_distro_h_eep.Fill(diff)
0532           border_diff_distro_m_eep.Fill(ix,iy,diff)
0533       elif (iz<0 and (Rsq<144.0 or Rsq>2500.0)):
0534           border_diff_distro_h_eem.Fill(diff)
0535           border_diff_distro_m_eem.Fill(ix,iy,diff)
0536
0537
0538     return border_diff_distro_h_eep, \
0539            border_diff_distro_m_eep, \
0540            border_diff_distro_h_eem, \
0541            border_diff_distro_m_eem