EcalAlCaRecoProducers/test/reRecoValidation.py

0001 import math
0002 import ROOT
0003 from ROOT import gSystem, TH2F, gStyle
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 FWLiteEnabler::enable()
0011
0012 #import EcalDetId
0013 #from DataFormats.EcalDetId import *
0014 #
0015 import sys,os
0016
0017 allRecHits=False
0018 eventListPrint=False
0019 lumi=-1
0020 #lumi=247
0021 eventNumbers=[
0022 343500081,
0023 343652254,
0024 344326842,
0025 344404776,
0026 345023043,
0027 344758536,
0028 344533037,
0029 343780994,
0030 344736267,
0031 344828101,
0032 344846924,
0033 344263362,
0034 344466033,
0035 343687292,
0036 344102889,
0037 344150484,
0038 343978421,
0039 344657593,
0040 344966978,
0041 345140483,
0042 343618665,
0043 344354767,
0044 344911692]
0045 #eventNumbers=[64944437] # evento con trackerDrivenEle
0046 eventNumbers=[]
0047 eventMin=-1
0048
0049 # for now look for events in two files with a given lumi section
0050 maxEvents=-1
0051 event_counter=0
0052
0053 for arg in sys.argv:
0054     if(arg=='AOD'):
0055         print "AOD"
0056         file="/tmp/"+os.environ["USER"]+"/rereco30Nov-AOD.root"
0057         file_format="AOD"
0058         break
0059     elif(arg=='AlcaFromAOD'):
0060         print "AlcaFromAOD"
0061         file="/tmp/"+os.environ["USER"]+"/AlcarecoFromAOD.root"
0062         file_format="AlcaFromAOD"
0063         break
0064     elif(arg=='AlcaFromAOD-recalib'):
0065         print "AlcaFromAOD-recalib"
0066         file="/tmp/"+os.environ["USER"]+"/AlcarecoFromAOD-recalib.root"
0067         file_format="AlcaFromAOD_Recalib"
0068         break
0069     elif(arg=='sandbox'):
0070         print 'sandbox'
0071         #        file="/tmp/"+os.environ["USER"]+"/sandbox.root"
0072         #        file="/tmp/"+os.environ["USER"]+"/alcaRecoSkim-2.root"
0073         file="./scratch/sandbox.root"
0074         #file="/tmp/"+os.environ["USER"]+"/alcaSkimSandbox-noADCtoGeV.root"
0075         file_format="sandbox"
0076 #sandbox"
0077         break
0078     elif(arg=='sandboxRereco'):
0079         print 'sandbox recalib'
0080
0081
0082         #file="/tmp/"+os.environ["USER"]+"/SANDBOX/sandboxRereco.root"
0083         file="./alcaSkimSandbox1.root"
0084 #        file="root://eoscms//eos/cms/store/group/alca_ecalcalib/sandboxRecalib/7TeV/fromUncalib/rereco30Nov/DoubleElectron-RUN2011B-v1/sandboxRereco-007.root"
0085         file_format="sandboxRecalib"
0086         break
0087     elif(arg=='RECO'):
0088         print 'RECO'
0089         file="root://eoscms//eos/cms/store/data/Run2012A/DoubleElectron/RECO/PromptReco-v1/000/193/336/BC442450-7997-E111-8177-003048D3C90E.root"
0090         #        file="/tmp/"+os.environ["USER"]+"/SANDBOX/RAW-RECO.root"
0091         file_format="RECO"
0092         break
0093     elif(arg=='MC'):
0094         print 'MC'
0095         file="/tmp/"+os.environ["USER"]+"/MC-AODSIM.root"
0096         file_format="AOD"
0097         break
0098     elif(arg=='ALCARECO'):
0099         print 'ALCARECO'
0100         file="alcaSkimSandbox_numEvent10000.root"
0101 #        file="alcaSkimSandbox.root"
0102         file_format="ALCARECO"
0103         break
0104
0105     else:
0106          continue
0107 #         exit(0)
0108
0109
0110 # AOD
0111 # sandbox
0112 #file_format = "ALCARECO"
0113 #file_format = 'sandbox'
0114 #file_format = 'AOD'
0115 #file_format = "AlcaFromAOD"
0116 #file_format = "AlcaFromAOD_Recalib"
0117
0118
0119
0120 #file='./alcaRecoSkim.root'
0121 # if(file_format == 'AOD'):
0122 #     file='/tmp/shervin/AOD.root'
0123 # elif(file_format == 'AlcaFromAOD'):
0124 #     file='/tmp/shervin/alcaRecoSkimFromAOD.root'
0125 # elif(file_format == 'AlcaFromAOD_Recalib'):
0126 #     file='/tmp/shervin/alcaRecoSkimFromAOD_Recalib.root'
0127 #     file='/tmp/shervin/alcaRecoSkimFromAOD_Recalib_IC_LC-2.root'
0128
0129 #file='./root/alcaRecoSkimFromAOD_LC_IC.root'
0130 #file='./root/alcaRecoSkimFromAOD_Recalib_IC_LC.root'
0131 #file='./root/alcaRecoSkimFromSandBox_Recalib.root'
0132 #file='/tmp/shervin/alcaRecoSkim-1.root'
0133 #file='/tmp/shervin/myAlcaRecoSkim.root'
0134
0135 print file
0136 events = Events (file)
0137
0138 handleElectrons = Handle('std::vector<reco::GsfElectron>')
0139 handleRecHitsEB = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0140 handleRecHitsEE = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0141 handleRecHitsES = Handle('edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> >')
0142 handleRhoFastJet = Handle('double')
0143
0144 if (file_format == 'ALCARECO'):
0145     processName="ALCARECO"
0146     electronTAG = 'gedGsfElectrons'
0147 elif(file_format == 'sandboxRecalib'):
0148     processName = "ALCARERECO"
0149     electronTAG = 'electronRecalibSCAssociator'
0150 #    recHitsTAG = "alCaIsolatedElectrons"
0151 elif(file_format == 'sandbox'):
0152     processName = "ALCASKIM"
0153 #    electronTAG = 'electronRecalibSCAssociator'
0154     electronTAG = 'gedGsfElectrons'
0155 elif(file_format == "AOD"):
0156     processName = "RECO"
0157     electronTAG = 'gedGsfElectrons'
0158 elif(file_format == "AlcaFromAOD"):
0159     processName = "ALCASKIM"
0160     electronTAG = 'gedGsfElectrons'
0161 elif(file_format == "AlcaFromAOD_Recalib"):
0162     electronTAG = 'electronRecalibSCAssociator'
0163     processName = 'ALCASKIM'
0164 elif(file_format == "RECO"):
0165     electronTAG = "gedGsfElectrons"
0166     processName = "RECO"
0167
0168
0169
0170
0171 EErecHitmap_ele1 = TH2F("EErecHitmap_ele1", "EErecHitmap_ele1",
0172                    100,0,100,
0173                    100,0,100)
0174
0175 EBrecHitmap_ele1 = TH2F("EBrecHitmap_ele1", "EBrecHitmap_ele1",
0176                    171,-85,85,
0177                    360,0,360)
0178
0179 EErecHitmap_ele2 = TH2F("EErecHitmap_ele2", "EErecHitmap_ele2",
0180                    100,0,100,
0181                    100,0,100)
0182
0183 EBrecHitmap_ele2 = TH2F("EBrecHitmap_ele2", "EBrecHitmap_ele2",
0184                    171,-85,85,
0185                    360,0,360)
0186
0187 print file_format, file, electronTAG, processName, maxEvents
0188
0189 print "run    lumi  event     isEB    energy     eSC  rawESC    e5x5    E_ES, etaEle, phiEle, etaSC, phiSC, clustersSize, nRecHits"
0190 for event in events:
0191
0192 #    if(event_counter % 100):
0193 #        print "[STATUS] event ", event_counter
0194
0195     if(maxEvents > 0 and event_counter > maxEvents):
0196         break
0197     if(eventListPrint==True):
0198         print event.eventAuxiliary().run(), event.eventAuxiliary().luminosityBlock(), event.eventAuxiliary().event()
0199         continue
0200     #if(event.eventAuxiliary.run()== 145351895):
0201     if lumi > 0 and int(event.eventAuxiliary().luminosityBlock()) != lumi :
0202             continue
0203
0204     if(len(eventNumbers)>0 and not (event.eventAuxiliary().event() in eventNumbers)):
0205        continue
0206
0207
0208         #    event.getByLabel(electronTAG, "", processName, handleElectrons)
0209     event.getByLabel(electronTAG, handleElectrons)
0210     #    print file_format, file, electronTAG
0211     electrons = handleElectrons.product()
0212
0213
0214     if(file_format=="AOD"):
0215         event.getByLabel("reducedEcalRecHitsEB", "", processName, handleRecHitsEB)
0216         event.getByLabel("reducedEcalRecHitsEE", "", processName, handleRecHitsEE)
0217         event.getByLabel("reducedEcalRecHitsES", "", processName, handleRecHitsES)
0218 #        print "##############",
0219         #        rhoTAG=edm.InputTag()
0220         #        rhoTAG=("kt6PFJets","rho","RECO")
0221 #        event.getByLabel("kt6PFJets","rho","RECO",handleRhoFastJet)
0222         # elif(file_format=="sandboxRecalib" or file_format=="RECO"):
0223     elif(file_format=="RECO"):
0224         event.getByLabel("ecalRecHit", "EcalRecHitsEB", processName, handleRecHitsEB)
0225         event.getByLabel("ecalRecHit", "EcalRecHitsEE", processName, handleRecHitsEE)
0226     elif(file_format=="ALCARECO" or file_format=="sandboxRecalib"):
0227         event.getByLabel("alCaIsolatedElectrons", "alCaRecHitsEB", processName, handleRecHitsEB)
0228         event.getByLabel("alCaIsolatedElectrons", "alCaRecHitsEE", processName, handleRecHitsEE)
0229 #        event.getByLabel("kt6PFJets","rho","RECO",handleRhoFastJet)
0230 #        event.getByLabel("reducedEcalRecHitsES", "", processName, handleRecHitsES)
0231
0232         #        event.getByLabel("ecalRecHit", "EcalRecHitsES", processName, handleRecHitsES)
0233
0234 #    elif(file_format=="sandbox"):
0235 #        event.getByLabel("ecalRecHit", "EcalRecHitsEB", processName, handleRecHitsEB)
0236 #        event.getByLabel("ecalRecHit", "EcalRecHitsEE", processName, handleRecHitsEE)
0237 #    else:
0238
0239
0240 #    print "Num of electrons: ",len(electrons)
0241     if(len(electrons)>=2):
0242      ele_counter=0
0243      for electron in electrons:
0244         if(not electron.ecalDrivenSeed()):
0245             print "trackerDriven",
0246 #            sys.exit(0)
0247         electron.superCluster().energy()
0248
0249
0250         #        ESrecHits = handleRecHitsES.product()
0251         #        if(abs(electron.eta()) > 1.566):
0252         #            for ESrecHit in ESrecHits:
0253         #                if(eventNumber >0):
0254         #                    esrecHit = ESDetId(ESrecHit.id().rawId())
0255         #                    print ESrecHit.id()(), esrecHit.strip(), esrecHit.six(), esrecHit.siy(), esrecHit.plane()
0256         print "------------------------------"
0257         if(not file_format=="sandbox"):
0258          if(electron.isEB()):
0259              recHits = handleRecHitsEB.product()
0260          else:
0261              recHits = handleRecHitsEE.product()
0262          nRecHits=0
0263          for recHit in recHits:
0264              nRecHits=nRecHits+1
0265              if(len(eventNumbers)==1):
0266                  if(electron.isEB()):
0267                      EBrecHit = EBDetId(recHit.id().rawId())
0268                      if(allRecHits):
0269                          if(ele_counter==0):
0270                              EBrecHitmap_ele1.Fill(EBrecHit.ieta(), EBrecHit.iphi(), recHit.energy());
0271                          elif(ele_counter==1):
0272                              EBrecHitmap_ele2.Fill(EBrecHit.ieta(), EBrecHit.iphi(), recHit.energy());
0273
0274                      print recHit.id()(), EBrecHit.ieta(), EBrecHit.iphi(), recHit.energy(), recHit.checkFlag(0)
0275                  else:
0276                      EErecHit = EEDetId(recHit.id().rawId())
0277                      if(allRecHits):
0278                          if(ele_counter==0):
0279                              EErecHitmap_ele1.Fill(EErecHit.ix(), EErecHit.iy(), recHit.energy());
0280                          elif(ele_counter==1):
0281                              EErecHitmap_ele2.Fill(EErecHit.ix(), EErecHit.iy(), recHit.energy());
0282                      print recHit.id()(), EErecHit.ix(), EErecHit.iy(), recHit.energy()
0283
0284          hits = electron.superCluster().hitsAndFractions()
0285          nRecHitsSC=0
0286          for hit in hits:
0287              nRecHitsSC=nRecHitsSC+1
0288              if(len(eventNumbers)==1):
0289                  if(electron.isEB()):
0290                      EBrecHit = EBDetId(hit.first.rawId())
0291                      if(not allRecHits):
0292                          if(ele_counter==0):
0293                              EBrecHitmap_ele1.Fill(EBrecHit.ieta(), EBrecHit.iphi(), hit.second*electron.superCluster().energy());
0294                          elif(ele_counter==1):
0295                              EBrecHitmap_ele2.Fill(EBrecHit.ieta(), EBrecHit.iphi(), hit.second*electron.superCluster().energy());
0296                      print "SC", (hit.first).rawId(), (hit.first)(), EBrecHit.ieta(), EBrecHit.iphi(),  EBrecHit.tower().iTT(), EBrecHit.ism(), EBrecHit.im()
0297                  else:
0298                      EErecHit = EEDetId(hit.first.rawId())
0299                      if(not allRecHits):
0300                          if(ele_counter==0):
0301                              EErecHitmap_ele1.Fill(EErecHit.ix(), EErecHit.iy(), recHit.energy());
0302                          elif(ele_counter==1):
0303                              EErecHitmap_ele2.Fill(EErecHit.ix(), EErecHit.iy(), recHit.energy());
0304
0305  #                print "SC ", (hit.first)()
0306
0307         print event.eventAuxiliary().run(), event.eventAuxiliary().luminosityBlock(), event.eventAuxiliary().event(),
0308         print "isEB=",electron.isEB(),
0309         print '{0:7.3f} {1:7.3f} {2:7.3f} {3:7.3f} {4:7.3f}'.format(electron.energy(), electron.superCluster().energy(), electron.superCluster().rawEnergy(), electron.e5x5(), electron.superCluster().preshowerEnergy()),
0310         print '{0:6.3f} {1:6.3f} {2:6.3f} {3:6.3f}'.format(electron.eta(), electron.phi(), electron.superCluster().eta(), electron.superCluster().phi()),
0311         print electron.superCluster().clustersSize(), nRecHits, nRecHitsSC
0312         ele_counter+=1
0313
0314
0315     event_counter+=1
0316
0317 print event_counter
0318
0319 # setting maps to -999
0320 gStyle.SetPaintTextFormat("1.1f")
0321 EBrecHitmap_ele1.SaveAs("EBrecHitmap_ele1.root")
0322 EErecHitmap_ele1.SaveAs("EErecHitmap_ele1.root")
0323 EBrecHitmap_ele2.SaveAs("EBrecHitmap_ele2.root")
0324 EErecHitmap_ele2.SaveAs("EErecHitmap_ele2.root")
0325
0326
0327