PhotonIdentification/test/testPhotonMVA_cfg.py

0001 import FWCore.ParameterSet.Config as cms
0002 from PhysicsTools.SelectorUtils.tools.vid_id_tools import *
0003 from Configuration.AlCa.GlobalTag import GlobalTag
0004
0005 process = cms.Process("PhotonMVANtuplizer")
0006
0007 process.load("Configuration.StandardSequences.GeometryDB_cff")
0008 process.load("Configuration.StandardSequences.MagneticField_cff")
0009 process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
0010
0011 process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc', '')
0012
0013 # File with the ID variables form the text file to include in the Ntuplizer
0014 mvaVariablesFile = "RecoEgamma/PhotonIdentification/data/PhotonMVAEstimatorRun2VariablesFall17V1p1.txt"
0015
0016 outputFile = "photon_ntuple.root"
0017
0018 process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000) )
0019
0020 process.source = cms.Source("PoolSource",
0021     fileNames = cms.untracked.vstring(
0022         '/store/mc/RunIIFall17MiniAODv2/GJet_Pt-20to40_DoubleEMEnriched_MGG-80toInf_TuneCP5_13TeV_Pythia8/MINIAODSIM/PU2017_12Apr2018_94X_mc2017_realistic_v14-v1/00000/00AE0E2A-6F42-E811-8EA2-0025905B85AA.root'
0023     )
0024 )
0025
0026 useAOD = False
0027
0028 from PhysicsTools.SelectorUtils.tools.vid_id_tools import *
0029 # turn on VID producer, indicate data format  to be
0030 # DataFormat.AOD or DataFormat.MiniAOD, as appropriate
0031 if useAOD == True :
0032     dataFormat = DataFormat.AOD
0033     input_tags = dict(
0034         src = cms.InputTag("gedPhotons"),
0035         vertices = cms.InputTag("offlinePrimaryVertices"),
0036         pileup = cms.InputTag("addPileupInfo"),
0037         genParticles = cms.InputTag("genParticles"),
0038         ebReducedRecHitCollection = cms.InputTag("reducedEcalRecHitsEB"),
0039         eeReducedRecHitCollection = cms.InputTag("reducedEcalRecHitsEE"),
0040     )
0041 else :
0042     dataFormat = DataFormat.MiniAOD
0043     input_tags = dict()
0044
0045 switchOnVIDPhotonIdProducer(process, dataFormat)
0046
0047 # define which IDs we want to produce
0048 my_id_modules = [
0049         'RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Spring16_nonTrig_V1_cff',
0050         'RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Fall17_94X_V1_cff',
0051         'RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Fall17_94X_V1p1_cff',
0052         'RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Fall17_94X_V2_cff',
0053                  ]
0054
0055 #add them to the VID producer
0056 for idmod in my_id_modules:
0057     setupAllVIDIdsInModule(process,idmod,setupVIDPhotonSelection)
0058
0059 process.ntuplizer = cms.EDAnalyzer('PhotonMVANtuplizer',
0060         phoMVAs              = cms.vstring(
0061                                           ),
0062         phoMVALabels         = cms.vstring(
0063                                           ),
0064         phoMVAValMaps        = cms.vstring(
0065                                            "photonMVAValueMapProducer:PhotonMVAEstimatorRun2Spring16NonTrigV1Values",
0066                                            "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v1Values",
0067                                            "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v1p1Values",
0068                                            "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v2Values",
0069                                            ),
0070         phoMVAValMapLabels   = cms.vstring(
0071                                            "Spring16NonTrigV1",
0072                                            "Fall17v1",
0073                                            "Fall17v1p1",
0074                                            "Fall17v2",
0075                                            ),
0076         phoMVACats           = cms.vstring(
0077                                            "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v1Categories",
0078                                            ),
0079         phoMVACatLabels      = cms.vstring(
0080                                            "PhoMVACats",
0081                                            ),
0082         variableDefinition = cms.string(mvaVariablesFile),
0083         #
0084         doEnergyMatrix = cms.bool(False), # disabled by default due to large size
0085         energyMatrixSize = cms.int32(2), # corresponding to 5x5
0086         #
0087         **input_tags
0088         )
0089 """
0090 The energy matrix is the n x n of raw rec-hit energies around the seed
0091 crystal.
0092
0093 The size of the energy matrix is controlled with the parameter
0094 "energyMatrixSize", which controlls the extension of crystals in each
0095 direction away from the seed, in other words n = 2 * energyMatrixSize + 1.
0096
0097 The energy matrix gets saved as a vector but you can easily unroll it
0098 to a two dimensional numpy array later, for example like that:
0099
0100 >>> import uproot
0101 >>> import numpy as np
0102 >>> import matplotlib.pyplot as plt
0103
0104 >>> tree = uproot.open("photon_ntuple.root")["ntuplizer/tree"]
0105 >>> n = 5
0106
0107 >>> for a in tree.array("ele_energyMatrix"):
0108 >>>     a = a.reshape((n,n))
0109 >>>     plt.imshow(np.log10(a))
0110 >>>     plt.colorbar()
0111 >>>     plt.show()
0112 """
0113
0114 process.TFileService = cms.Service("TFileService", fileName = cms.string(outputFile))
0115
0116 process.p = cms.Path(process.egmPhotonIDSequence * process.ntuplizer)