Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-03 05:26:36

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 ROOT.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