ShowerLibraryProducer/test/AnalyzeTuples.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    AnalyzeTuples
0004 // Class:      AnalyzeTuples
0005 //
0006
0007 /**\class AnalyzeTuples AnalyzeTuples.cc Analysis/AnalyzeTuples/src/AnalyzeTuples.cc
0008
0009 Description: <one line class summary>
0010
0011 Implementation:
0012 <Notes on implementation>
0013 */
0014 //
0015 // Original Author: Taylan Yetkin
0016 // Created: Tue Feb 10 08:43:07 CST 2009
0017 //
0018 //
0019
0020 #include <memory>
0021 #include <vector>
0022 #include <string>
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/Utilities/interface/InputTag.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 #include "SimDataFormats/CaloHit/interface/HFShowerPhoton.h"
0034 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
0035
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0039 #include "TFile.h"
0040 #include "TTree.h"
0041 #include "TBranch.h"
0042 #include "TH1F.h"
0043 #include "TH1I.h"
0044
0045 class AnalyzeTuples : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0046 public:
0047   explicit AnalyzeTuples(const edm::ParameterSet&);
0048   ~AnalyzeTuples() override = default;
0049
0050 private:
0051   void beginJob() override;
0052   void analyze(const edm::Event&, const edm::EventSetup&) override;
0053   void endJob() override;
0054   void loadEventInfo(TBranch*);
0055   void getRecord(int type, int record);
0056   TFile* hf;
0057   TBranch *emBranch, *hadBranch;
0058
0059   int nMomBin, totEvents, evtPerBin;
0060   float libVers, listVersion;
0061   std::vector<double> pmom;
0062   std::vector<HFShowerPhoton> photon;
0063
0064   TH1I* hNPELongElec[12];
0065   TH1I* hNPEShortElec[12];
0066   TH1I* hNPELongPion[12];
0067   TH1I* hNPEShortPion[12];
0068 };
0069
0070 AnalyzeTuples::AnalyzeTuples(const edm::ParameterSet& iConfig) {
0071   usesResource(TFileService::kSharedResource);
0072
0073   edm::LogVerbatim("HcalSim") << "analyzetuples a buraya girdi";
0074   edm::ParameterSet m_HS = iConfig.getParameter<edm::ParameterSet>("HFShowerLibrary");
0075   edm::FileInPath fp = m_HS.getParameter<edm::FileInPath>("FileName");
0076   std::string pTreeName = fp.fullPath();
0077   std::string emName = m_HS.getParameter<std::string>("TreeEMID");
0078   std::string hadName = m_HS.getParameter<std::string>("TreeHadID");
0079   std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>(
0080       "BranchEvt", "HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
0081   std::string branchPre = m_HS.getUntrackedParameter<std::string>("BranchPre", "HFShowerPhotons_hfshowerlib_");
0082   std::string branchPost = m_HS.getUntrackedParameter<std::string>("BranchPost", "_R.obj");
0083
0084   if (pTreeName.find(".") == 0)
0085     pTreeName.erase(0, 2);
0086   const char* nTree = pTreeName.c_str();
0087   hf = TFile::Open(nTree);
0088
0089   if (!hf->IsOpen()) {
0090     edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
0091     throw cms::Exception("Unknown", "HFShowerLibrary") << "Opening of " << pTreeName << " fails\n";
0092   } else {
0093     edm::LogInfo("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
0094   }
0095
0096   TTree* event = (TTree*)hf->Get("Events");
0097   if (event) {
0098     std::string info = branchEvInfo + branchPost;
0099     TBranch* evtInfo = event->GetBranch(info.c_str());
0100     if (evtInfo) {
0101       loadEventInfo(evtInfo);
0102     } else {
0103       edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
0104                                 << " Branch does not exist in Event";
0105       throw cms::Exception("Unknown", "HFShowerLibrary") << "Event information absent\n";
0106     }
0107   } else {
0108     edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
0109                               << "exist";
0110     throw cms::Exception("Unknown", "HFShowerLibrary") << "Events tree absent\n";
0111   }
0112
0113   edm::LogInfo("HFShower") << "HFShowerLibrary: Library " << libVers << " ListVersion " << listVersion
0114                            << " Events Total " << totEvents << " and " << evtPerBin << " per bin";
0115   edm::LogInfo("HFShower") << "HFShowerLibrary: Energies (GeV) with " << nMomBin << " bins";
0116   for (int i = 0; i < nMomBin; i++)
0117     edm::LogInfo("HFShower") << "HFShowerLibrary: pmom[" << i << "] = " << pmom[i] << " GeV";
0118
0119   std::string nameBr = branchPre + emName + branchPost;
0120   emBranch = event->GetBranch(nameBr.c_str());
0121   nameBr = branchPre + hadName + branchPost;
0122   hadBranch = event->GetBranch(nameBr.c_str());
0123   edm::LogInfo("HFShower") << "HFShowerLibrary:Branch " << emName << " has " << emBranch->GetEntries()
0124                            << " entries and Branch " << hadName << " has " << hadBranch->GetEntries() << " entries";
0125   edm::LogInfo("HFShower") << "HFShowerLibrary::No packing information -"
0126                            << " Assume x, y, z are not in packed form";
0127 }
0128
0129 void AnalyzeTuples::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0130   for (int ibin = 0; ibin < 12; ++ibin) {
0131     int min = evtPerBin * (ibin);
0132     int max = evtPerBin * (ibin + 1);
0133     for (int i = min; i < max; ++i) {
0134       getRecord(0, i);
0135       int npe_long = 0;
0136       int npe_short = 0;
0137       edm::LogVerbatim("HcalSim") << "phptons size" << photon.size();
0138       for (int j = 0; j < int(photon.size()); ++j) {
0139         //int depth = 0;
0140         if (photon[j].z() < 0) {
0141           //depth = 2;
0142           ++npe_short;
0143         } else {
0144           //depth = 1;
0145           ++npe_long;
0146           edm::LogVerbatim("HcalSim") << photon[j].z();
0147         }
0148       }
0149       hNPELongElec[ibin]->Fill(npe_long);
0150       edm::LogVerbatim("HcalSim") << ibin << npe_long;
0151       hNPEShortElec[ibin]->Fill(npe_short);
0152     }
0153   }
0154   for (int ibin = 0; ibin < 12; ++ibin) {
0155     int min = evtPerBin * (ibin);
0156     int max = evtPerBin * (ibin + 1);
0157     for (int i = min; i < max; ++i) {
0158       getRecord(1, i);
0159       int npe_long = 0;
0160       int npe_short = 0;
0161       for (int j = 0; j < int(photon.size()); ++j) {
0162         //int depth = 0;
0163         if (photon[j].z() < 0) {
0164           //depth = 2;
0165           ++npe_short;
0166         } else {
0167           //depth = 1;
0168           ++npe_long;
0169         }
0170       }
0171       hNPELongPion[ibin]->Fill(npe_long);
0172       hNPEShortPion[ibin]->Fill(npe_short);
0173     }
0174   }
0175 }
0176
0177 void AnalyzeTuples::beginJob() {
0178   edm::Service<TFileService> fs;
0179   TFileDirectory HFDir = fs->mkdir("HF");
0180   char title[128];
0181   for (int i = 0; i < int(pmom.size()); ++i) {
0182     sprintf(title, "NPELongElec_Mom_%i", int(pmom[i]));
0183     int maxBin = int(pmom[i] + 50);
0184     hNPELongElec[i] = HFDir.make<TH1I>(title, "NPE Long", 140, 0, 140);
0185     sprintf(title, "NPEShortElec_Mom_%i", int(pmom[i]));
0186     hNPEShortElec[i] = HFDir.make<TH1I>(title, "NPE Short", maxBin, 0, maxBin);
0187     sprintf(title, "NPELongPion_Mom_%i", int(pmom[i]));
0188     hNPELongPion[i] = HFDir.make<TH1I>(title, "NPE Long", maxBin, 0, maxBin);
0189     sprintf(title, "NPEShortPion_Mom_%i", int(pmom[i]));
0190     hNPEShortPion[i] = HFDir.make<TH1I>(title, "NPE Short", maxBin, 0, maxBin);
0191   }
0192 }
0193
0194 void AnalyzeTuples::endJob() {}
0195
0196 void AnalyzeTuples::loadEventInfo(TBranch* branch) {
0197   std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
0198   branch->SetAddress(&eventInfoCollection);
0199   branch->GetEntry(0);
0200   edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
0201                            << " EventInfo Collection of size " << eventInfoCollection.size() << " records";
0202   totEvents = eventInfoCollection[0].totalEvents();
0203   nMomBin = eventInfoCollection[0].numberOfBins();
0204   evtPerBin = eventInfoCollection[0].eventsPerBin();
0205   libVers = eventInfoCollection[0].showerLibraryVersion();
0206   listVersion = eventInfoCollection[0].physListVersion();
0207   pmom = eventInfoCollection[0].energyBins();
0208 }
0209
0210 void AnalyzeTuples::getRecord(int type, int record) {
0211   int nrc = record - 1;
0212   photon.clear();
0213   if (type > 0) {
0214     hadBranch->SetAddress(&photon);
0215     hadBranch->GetEntry(nrc);
0216   } else {
0217     emBranch->SetAddress(&photon);
0218     emBranch->GetEntry(nrc);
0219   }
0220 #ifdef EDM_ML_DEBUG
0221   int nPhoton = photon.size();
0222   edm::LogVerbatim("HFShower") << "HFShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
0223                                << nPhoton << " photons";
0224   for (int j = 0; j < nPhoton; j++)
0225     edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon[j];
0226 #endif
0227 }
0228
0229 // define this as a plug-in
0230 DEFINE_FWK_MODULE(AnalyzeTuples);