EcalAlCaRecoProducers/test/recHitsValidation.py

0001 import math
0002 import ROOT
0003 from ROOT import gSystem, TH2F
0004 from DataFormats.FWLite import Events, Handle
0005 from PhysicsTools.PythonAnalysis import *
0006 #import print_options
0007
0008 #print_options.set_float_precision(4)
0009 gSystem.Load("libFWCoreFWLite.so")
0010 ROOT.FWLiteEnabler.enable()
0011
0012 #import EcalDetId
0013 #from DataFormats.EcalDetId import *
0014 #
0015 import sys,os
0016
0017 allRecHits=False
0018 lumi=393
0019 lumi=-1
0020 lumi=466
0021 eventNumber=-1
0022 eventNumber=685233340
0023 #lumi=376
0024 #eventNumber=344404776
0025 #eventNumber=552337102
0026 #eventNumber=577412878
0027
0028 #eventNumber=347006644
0029 #eventNumber=344911692
0030
0031 lumi=374
0032 eventNumber=548668449
0033 # kRecovered lumi=376 eventNumber=552337102
0034 # kRecovered lumi=393 eventNumber=576932539
0035 # kRecovered lumi=394 eventNumber=578483322
0036 lumi=394
0037 eventNumber=578490502
0038 lumi=394
0039 eventNumber=579700192
0040 lumi=395
0041 eventNumber=579843406
0042 lumi=401
0043 eventNumber=588810213
0044 lumi=402
0045 eventNumber=591275401
0046 # kRecovered lumi=403 eventNumber=593293410
0047 lumi=403
0048 eventNumber=591888388
0049 # kRecovered lumi=406 eventNumber=597401546
0050 # kRecovered lumi=407 eventNumber=598290564
0051 lumi=415
0052 eventNumber=610541757
0053 lumi=415
0054 eventNumber=610541436
0055 lumi=416
0056 eventNumber=612542602
0057 # kRecovered lumi=419 eventNumber=616682572
0058 # kRecovered lumi=419 eventNumber=615876590
0059 lumi=422
0060 eventNumber=620689835
0061
0062 lumi=466
0063 eventNumber=685900276
0064 lumi=467
0065 eventNumber=687572911
0066 lumi=472
0067 eventNumber=694966852
0068 # kRecovered
0069 lumi=466
0070 eventNumber=685233340
0071
0072 eventMin=-1
0073
0074 # for now look for events in two files with a given lumi section
0075 maxEvents=-1
0076 event_counter=0
0077
0078 for arg in sys.argv:
0079     if (arg=='testAlca1'):
0080         print("testAlca1")
0081         file="/tmp/"+os.environ["USER"]+"/testAlca1.root"
0082         file_format = "AlcaFromAOD"
0083         break
0084     elif(arg=='testAlca2'):
0085         print('testAlca2')
0086         file="/tmp/"+os.environ["USER"]+"/testAlca2.root"
0087         file_format = "AlcaFromAOD_Recalib"
0088         break
0089     elif(arg=='testAlca3'):
0090         print('testAlca3')
0091         file="/tmp/"+os.environ["USER"]+"/testAlca3.root"
0092         file_format = "AlcaFromAOD_Recalib"
0093         break
0094     elif(arg=='testAlca4'):
0095         print('testAlca4')
0096         file="/tmp/"+os.environ["USER"]+"/testAlca4.root"
0097         file_format = "AlcaFromAOD_Recalib"
0098         break
0099     elif(arg=='AOD'):
0100         print("AOD")
0101         file="/tmp/"+os.environ["USER"]+"/rereco30Nov-AOD.root"
0102         file_format="AOD"
0103         break
0104     elif(arg=='AlcaFromAOD'):
0105         print("AlcaFromAOD")
0106         file="/tmp/"+os.environ["USER"]+"/AlcarecoFromAOD.root"
0107         file_format="AlcaFromAOD"
0108         break
0109     elif(arg=='AlcaFromAOD-recalib'):
0110         print("AlcaFromAOD-recalib")
0111         file="/tmp/"+os.environ["USER"]+"/AlcarecoFromAOD-recalib.root"
0112         file_format="AlcaFromAOD_Recalib"
0113         break
0114     elif(arg=='sandbox'):
0115         print('sandbox')
0116         #        file="/tmp/"+os.environ["USER"]+"/sandbox.root"
0117         #        file="/tmp/"+os.environ["USER"]+"/alcaRecoSkim-2.root"
0118         file="/tmp/"+os.environ["USER"]+"/alcaSkimSandbox.root"
0119         #file="/tmp/"+os.environ["USER"]+"/alcaSkimSandbox-noADCtoGeV.root"
0120         file_format="sandbox"
0121 #sandbox"
0122         break
0123     elif(arg=='sandboxRecalib'):
0124         print('sandbox recalib')
0125 #        file="/tmp/"+os.environ["USER"]+"/Test-RecalibSandbox-GT_IC_LC.root"
0126 #        file="/tmp/"+os.environ["USER"]+"/SandboxReReco-noADCtoGeV.root"
0127         file="/tmp/"+os.environ["USER"]+"/SandboxReReco.root"
0128 #        file="/tmp/"+os.environ["USER"]+"/SandboxReReco-GTprompt-ALCARECO.root"
0129 #        file="/tmp/"+os.environ["USER"]+"/SANDBOX/SandboxReReco.root"
0130 #        file="/tmp/"+os.environ["USER"]+"/SANDBOX/SandboxReReco-noFranz2.root"
0131         file_format="sandboxRecalib"
0132
0133         break
0134     elif(arg=='RECO'):
0135         print('RECO')
0136         file="/tmp/"+os.environ["USER"]+"/SANDBOX/RAW-RECO.root"
0137         file_format="RECO"
0138         break
0139     else:
0140          continue
0141 #         exit(0)
0142
0143
0144
0145 events = Events (file)
0146 print(file)
0147 handleElectrons = Handle('std::vector<reco::GsfElectron>')
0148 handleRecHitsEB = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0149 handleRecHitsEE = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0150 handleRecHitsES = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0151
0152 handleRecHitsEB_RECO = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0153 handleRecHitsEE_RECO = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0154 handleRecHitsEB_ALCASKIM = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0155 handleRecHitsEE_ALCASKIM = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0156 handleRecHitsEB_ALCARECO = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0157 handleRecHitsEE_ALCARECO = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0158
0159 if (file_format == 'ALCARECO'):
0160     processName="ALCASKIM"
0161     electronTAG = 'electronRecalibSCAssociator'
0162 elif(file_format == 'sandboxRecalib'):
0163     processName = "ALCARERECO"
0164     electronTAG = 'electronRecalibSCAssociator'
0165     recHitsTAG = "alCaIsolatedElectrons"
0166 elif(file_format == 'sandbox'):
0167     processName = "ALCASKIM"
0168     electronTAG = 'electronRecalibSCAssociator'
0169 elif(file_format == "AOD"):
0170     processName = "RECO"
0171     electronTAG = 'gedGsfElectrons'
0172 elif(file_format == "AlcaFromAOD"):
0173     processName = "ALCASKIM"
0174     electronTAG = 'gedGsfElectrons'
0175 elif(file_format == "AlcaFromAOD_Recalib"):
0176     electronTAG = 'electronRecalibSCAssociator'
0177     processName = 'ALCASKIM'
0178 elif(file_format == "RECO"):
0179     electronTAG = "gedGsfElectrons"
0180     processName = "RECO"
0181
0182
0183
0184
0185 EErecHitmap_ele1 = TH2F("EErecHitmap_ele1", "EErecHitmap_ele1",
0186                    100,0,100,
0187                    100,0,100)
0188
0189 EBrecHitmap_ele1 = TH2F("EBrecHitmap_ele1", "EBrecHitmap_ele1",
0190                    171,-85,85,
0191                    360,0,360)
0192
0193 EErecHitmap_ele2 = TH2F("EErecHitmap_ele2", "EErecHitmap_ele2",
0194                    100,0,100,
0195                    100,0,100)
0196
0197 EBrecHitmap_ele2 = TH2F("EBrecHitmap_ele2", "EBrecHitmap_ele2",
0198                    171,-85,85,
0199                    360,0,360)
0200
0201 print(file_format, file, electronTAG, processName, maxEvents)
0202
0203 print("run\tlumi, event, energy, eSC, rawESC, e5x5, E_ES, etaEle, phiEle, etaSC, phiSC, clustersSize, nRecHits")
0204 for event in events:
0205
0206     if(maxEvents > 0 and event_counter > maxEvents):
0207         break
0208     #if(event.eventAuxiliary.run()== 145351895):
0209     if lumi > 0 and int(event.eventAuxiliary().luminosityBlock()) != lumi :
0210             continue
0211
0212     if(eventNumber > 0 and event.eventAuxiliary().event()!= eventNumber ):
0213         continue
0214
0215         #    event.getByLabel(electronTAG, "", processName, handleElectrons)
0216     event.getByLabel(electronTAG, handleElectrons)
0217     #    print(file_format, file, electronTAG        )
0218     electrons = handleElectrons.product()
0219
0220     #    event.getByLabel("reducedEcalRecHitsEB", "", processName, handleRecHitsEB)
0221     #    event.getByLabel("reducedEcalRecHitsEE", "", processName, handleRecHitsEE)
0222     if(file_format=="sandbox"):
0223         event.getByLabel("ecalRecHit",    "EcalRecHitsEB",   "RECO", handleRecHitsEB_RECO)
0224         event.getByLabel("ecalRecHit",    "EcalRecHitsEE",   "RECO", handleRecHitsEE_RECO)
0225         event.getByLabel("ecalRecHit",    "EcalRecHitsEB",   "ALCASKIM", handleRecHitsEB_ALCASKIM)
0226         event.getByLabel("ecalRecHit",    "EcalRecHitsEE",   "ALCASKIM", handleRecHitsEE_ALCASKIM)
0227     else:
0228         event.getByLabel("alCaIsolatedElectrons","alcaBarrelHits", "ALCARERECO", handleRecHitsEB_ALCARECO)
0229         event.getByLabel("alCaIsolatedElectrons","alcaEndcapHits", "ALCARERECO", handleRecHitsEE_ALCARECO)
0230
0231     if(file_format=="sandbox"):
0232
0233        for electron in electrons:
0234            if(abs(electron.eta()) < 1.4442):
0235                recHits_RECO = handleRecHitsEB_RECO.product()
0236                recHits_ALCASKIM = handleRecHitsEB_ALCASKIM.product()
0237            else:
0238                recHits_RECO = handleRecHitsEE_RECO.product()
0239                recHits_ALCASKIM = handleRecHitsEE_ALCASKIM.product()
0240
0241            nRecHits_RECO=0
0242            for recHit in recHits_RECO:
0243                nRecHits_RECO=nRecHits_RECO+1
0244 #                if(recHit.checkFlag(EcalRecHit.kTowerRecovered)):
0245 #                   print(recHit.id().rawId())
0246
0247            nRecHits_ALCASKIM=0
0248            for recHit in recHits_ALCASKIM:
0249                nRecHits_ALCASKIM=nRecHits_ALCASKIM+1
0250                #               if(recHit.checkFlag(EcalRecHit.kTowerRecovered)):
0251                print(recHit.id().rawId(), recHit.checkFlag(EcalRecHit.kTowerRecovered))
0252
0253            if(nRecHits_ALCASKIM != nRecHits_RECO):
0254                print(nRecHits_RECO, nRecHits_ALCASKIM)
0255                print(recHits_RECO)
0256                print("------------------------------")
0257                print(recHits_ALCASKIM)
0258     else:
0259        for electron in electrons:
0260            if(abs(electron.eta()) < 1.4442):
0261                recHits_ALCARECO = handleRecHitsEB_ALCARECO.product()
0262            else:
0263                recHits_ALCARECO = handleRecHitsEE_ALCARECO.product()
0264
0265            nRecHits_ALCARECO=0
0266            for recHit in recHits_ALCARECO:
0267                nRecHits_ALCARECO=nRecHits_ALCARECO+1
0268                #               if(recHit.checkFlag(EcalRecHit.kTowerRecovered)):
0269                print(recHit.id().rawId(), recHit.checkFlag(EcalRecHit.kTowerRecovered))
0270
0271 print(event_counter)
0272
0273
0274