Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-17 01:48:02

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015 
0016 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
0017 #include "SimDataFormats/CaloHit/interface/HFShowerPhoton.h"
0018 
0019 #include "CLHEP/Units/SystemOfUnits.h"
0020 #include "TFile.h"
0021 #include "TTree.h"
0022 #include "TBranch.h"
0023 #include "TH1F.h"
0024 
0025 #include <string>
0026 #include <vector>
0027 
0028 class HFShowerLibraryAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0029 public:
0030   HFShowerLibraryAnalyzer(const edm::ParameterSet& ps);
0031   ~HFShowerLibraryAnalyzer() override;
0032   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0033 
0034 protected:
0035   void beginJob() override {}
0036   void analyze(edm::Event const&, edm::EventSetup const&) override {}
0037   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0038   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0039 
0040 private:
0041   void bookHistos();
0042   void getRecord(int, int);
0043   void loadEventInfo(TBranch*);
0044 
0045   TFile* hf_;
0046   TBranch *emBranch_, *hadBranch_;
0047   bool verbose_, newForm_, v3version_;
0048   int nMomBin_, totEvents_, evtPerBin_;
0049   float libVers_, listVersion_;
0050   std::vector<double> pmom_;
0051   HFShowerPhotonCollection* photo_;
0052   HFShowerPhotonCollection photon_;
0053   std::vector<TH1F*> h_x_[2], h_y_[2], h_z_[2], h_t_[2], h_l_[2], h_nl_[2], h_ns_[2];
0054 };
0055 
0056 HFShowerLibraryAnalyzer::HFShowerLibraryAnalyzer(edm::ParameterSet const& ps)
0057     : hf_(nullptr), emBranch_(nullptr), hadBranch_(nullptr) {
0058   usesResource(TFileService::kSharedResource);
0059 
0060   edm::FileInPath fp("SimG4CMS/Calo/data/" + ps.getParameter<std::string>("FileName"));
0061   std::string pTreeName = fp.fullPath();
0062   std::string emName = ps.getParameter<std::string>("TreeEMID");
0063   std::string hadName = ps.getParameter<std::string>("TreeHadID");
0064   std::string branchEvInfo = ps.getParameter<std::string>("BranchEvt");
0065   std::string branchPre = ps.getParameter<std::string>("BranchPre");
0066   std::string branchPost = ps.getParameter<std::string>("BranchPost");
0067   verbose_ = ps.getParameter<bool>("Verbosity");
0068   evtPerBin_ = ps.getParameter<bool>("EventPerBin");
0069 
0070   if (pTreeName.find('.') == 0)
0071     pTreeName.erase(0, 2);
0072   const char* nTree = pTreeName.c_str();
0073   hf_ = TFile::Open(nTree);
0074 
0075   if (hf_->IsOpen()) {
0076     edm::LogVerbatim("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
0077     newForm_ = (branchEvInfo.empty());
0078     TTree* event = ((newForm_) ? static_cast<TTree*>(hf_->Get("HFSimHits")) : static_cast<TTree*>(hf_->Get("Events")));
0079     if (event) {
0080       TBranch* evtInfo(nullptr);
0081       if (!newForm_) {
0082         std::string info = branchEvInfo + branchPost;
0083         evtInfo = event->GetBranch(info.c_str());
0084       }
0085       if (evtInfo || newForm_) {
0086         loadEventInfo(evtInfo);
0087       } else {
0088         edm::LogError("HFShower") << "HFShowerLibrayEventInfo Branch does not exist in Event";
0089         throw cms::Exception("Unknown", "HFShowerLibraryAnalyzer") << "Event information absent\n";
0090       }
0091     } else {
0092       edm::LogError("HFShower") << "Events Tree does not exist";
0093       throw cms::Exception("Unknown", "HFShowerLibraryAnalyzer") << "Events tree absent\n";
0094     }
0095 
0096     std::stringstream ss;
0097     ss << "HFShowerLibraryAnalyzer: Library " << libVers_ << " ListVersion " << listVersion_ << " Events Total "
0098        << totEvents_ << " and " << evtPerBin_ << " per bin\n";
0099     ss << "HFShowerLibraryAnalyzer: Energies (GeV) with " << nMomBin_ << " bins\n";
0100     for (int i = 0; i < nMomBin_; ++i) {
0101       if (i / 10 * 10 == i && i > 0) {
0102         ss << "\n";
0103       }
0104       ss << "  " << pmom_[i] / CLHEP::GeV;
0105     }
0106     edm::LogVerbatim("HFShower") << ss.str();
0107 
0108     std::string nameBr = branchPre + emName + branchPost;
0109     emBranch_ = event->GetBranch(nameBr.c_str());
0110     if (verbose_)
0111       emBranch_->Print();
0112     nameBr = branchPre + hadName + branchPost;
0113     hadBranch_ = event->GetBranch(nameBr.c_str());
0114     if (verbose_)
0115       hadBranch_->Print();
0116 
0117     v3version_ = (emBranch_->GetClassName() == std::string("vector<float>")) ? true : false;
0118     edm::LogVerbatim("HFShower")
0119         << " HFShowerLibraryAnalyzer:Branch " << emName << " has " << emBranch_->GetEntries() << " entries and Branch "
0120         << hadName << " has " << hadBranch_->GetEntries()
0121         << " entries\n HFShowerLibraryAnalyzer::No packing information - Assume x, y, z are not in packed form";
0122     photo_ = new HFShowerPhotonCollection;
0123 
0124     bookHistos();
0125   } else {
0126     edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
0127   }
0128 }
0129 
0130 HFShowerLibraryAnalyzer::~HFShowerLibraryAnalyzer() {
0131   if (hf_)
0132     hf_->Close();
0133   delete photo_;
0134 }
0135 
0136 void HFShowerLibraryAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0137   edm::ParameterSetDescription desc;
0138   desc.add<std::string>("FileName", "HFShowerLibrary_10000.root");
0139   desc.add<std::string>("TreeEMID", "emParticles");
0140   desc.add<std::string>("TreeHadID", "hadParticles");
0141   desc.add<std::string>("BranchEvt", "");
0142   desc.add<std::string>("BranchPre", "");
0143   desc.add<std::string>("BranchPost", "");
0144   desc.add<bool>("Verbosity", false);
0145   desc.add<int>("EventPerBin", 10000);
0146   descriptions.add("hfShowerLibaryAnalysis", desc);
0147 }
0148 
0149 void HFShowerLibraryAnalyzer::bookHistos() {
0150   edm::Service<TFileService> tfile;
0151   if (!tfile.isAvailable())
0152     throw cms::Exception("Unknown", "HFShowerLibraryAnalyzer")
0153         << "TFileService unavailable: please add it to config file";
0154   char name[20], title[40], titlx[120];
0155   TH1F* hist;
0156   for (int i = 0; i < 2; ++i) {
0157     std::string type = (i == 0) ? "EM" : "Hadron";
0158     for (int k = 0; k < nMomBin_; ++k) {
0159       sprintf(title, "Showers for p = %6.1f GeV", pmom_[k] / CLHEP::GeV);
0160       sprintf(name, "X%d%d", k, i);
0161       sprintf(titlx, "x coordinate of %s shower (mm)", type.c_str());
0162       hist = tfile->make<TH1F>(name, title, 200, -500.0, 500.0);
0163       hist->GetXaxis()->SetTitle(titlx);
0164       hist->GetYaxis()->SetTitle("Shower");
0165       hist->GetYaxis()->SetTitleOffset(1.2);
0166       hist->Sumw2();
0167       h_x_[i].emplace_back(hist);
0168       sprintf(name, "Y%d%d", k, i);
0169       sprintf(titlx, "y coordinate of %s shower (mm)", type.c_str());
0170       hist = tfile->make<TH1F>(name, title, 200, -500.0, 500.0);
0171       hist->GetXaxis()->SetTitle(titlx);
0172       hist->GetYaxis()->SetTitle("Shower");
0173       hist->GetYaxis()->SetTitleOffset(1.2);
0174       hist->Sumw2();
0175       h_y_[i].emplace_back(hist);
0176       sprintf(name, "Z%d%d", k, i);
0177       sprintf(titlx, "z coordinate of %s shower (mm)", type.c_str());
0178       hist = tfile->make<TH1F>(name, title, 200, -2000.0, 2000.0);
0179       hist->GetXaxis()->SetTitle(titlx);
0180       hist->GetYaxis()->SetTitle("Shower");
0181       hist->GetYaxis()->SetTitleOffset(1.2);
0182       hist->Sumw2();
0183       h_z_[i].emplace_back(hist);
0184       sprintf(name, "T%d%d", k, i);
0185       sprintf(titlx, "Time of %s shower (ns)", type.c_str());
0186       hist = tfile->make<TH1F>(name, title, 200, 0.0, 50.0);
0187       hist->GetXaxis()->SetTitle(titlx);
0188       hist->GetYaxis()->SetTitle("Shower");
0189       hist->GetYaxis()->SetTitleOffset(1.2);
0190       hist->Sumw2();
0191       h_t_[i].emplace_back(hist);
0192       sprintf(name, "L%d%d", k, i);
0193       sprintf(titlx, "Lambda of %s shower photon (nm)", type.c_str());
0194       hist = tfile->make<TH1F>(name, title, 200, 300.0, 800.0);
0195       hist->GetXaxis()->SetTitle(titlx);
0196       hist->GetYaxis()->SetTitle("Shower");
0197       hist->GetYaxis()->SetTitleOffset(1.2);
0198       h_l_[i].emplace_back(hist);
0199       hist->Sumw2();
0200       sprintf(name, "NL%d%d", k, i);
0201       sprintf(titlx, "Number of %s shower photons in long fibers", type.c_str());
0202       hist = tfile->make<TH1F>(name, title, 3000, 0.0, 3000.0);
0203       hist->GetXaxis()->SetTitle(titlx);
0204       hist->GetYaxis()->SetTitle("Shower");
0205       hist->GetYaxis()->SetTitleOffset(1.2);
0206       h_nl_[i].emplace_back(hist);
0207       sprintf(name, "NS%d%d", k, i);
0208       sprintf(titlx, "Number of %s shower photons in short fibes", type.c_str());
0209       hist = tfile->make<TH1F>(name, title, 3000, 0.0, 3000.0);
0210       hist->GetXaxis()->SetTitle(titlx);
0211       hist->GetYaxis()->SetTitle("Shower");
0212       hist->GetYaxis()->SetTitleOffset(1.2);
0213       h_ns_[i].emplace_back(hist);
0214     }
0215   }
0216 
0217   // Now fill them up
0218   for (int type = 0; type < 2; ++type) {
0219     for (int k = 0; k < nMomBin_; ++k) {
0220       for (int j = 0; j < evtPerBin_; ++j) {
0221         int irc = k * evtPerBin_ + j + 1;
0222         getRecord(type, irc);
0223         int nPhoton = (newForm_) ? photo_->size() : photon_.size();
0224         int nlong = 0, nshort = 0;
0225         for (int i = 0; i < nPhoton; i++) {
0226           if (newForm_) {
0227             if (photo_->at(i).z() > 0) {
0228               nlong++;
0229             } else {
0230               nshort++;
0231             }
0232             h_x_[type][k]->Fill((photo_->at(i)).x());
0233             h_y_[type][k]->Fill((photo_->at(i)).y());
0234             h_z_[type][k]->Fill((photo_->at(i)).z());
0235             h_t_[type][k]->Fill((photo_->at(i)).t());
0236             h_l_[type][k]->Fill((photo_->at(i)).lambda());
0237           } else {
0238             if (photon_[i].z() > 0) {
0239               nlong++;
0240             } else {
0241               nshort++;
0242             }
0243             h_x_[type][k]->Fill(photon_[i].x());
0244             h_y_[type][k]->Fill(photon_[i].y());
0245             h_z_[type][k]->Fill(photon_[i].z());
0246             h_t_[type][k]->Fill(photon_[i].t());
0247             h_l_[type][k]->Fill(photon_[i].lambda());
0248           }
0249         }
0250         h_nl_[type][k]->Fill(nlong);
0251         h_ns_[type][k]->Fill(nshort);
0252       }
0253     }
0254   }
0255 }
0256 
0257 void HFShowerLibraryAnalyzer::getRecord(int type, int record) {
0258   int nrc = record - 1;
0259   photon_.clear();
0260   photo_->clear();
0261   if (type > 0) {
0262     if (newForm_) {
0263       if (!v3version_) {
0264         hadBranch_->SetAddress(&photo_);
0265         hadBranch_->GetEntry(nrc + totEvents_);
0266       } else {
0267         std::vector<float> t;
0268         std::vector<float>* tp = &t;
0269         hadBranch_->SetAddress(&tp);
0270         hadBranch_->GetEntry(nrc + totEvents_);
0271         unsigned int tSize = t.size() / 5;
0272         photo_->reserve(tSize);
0273         for (unsigned int i = 0; i < tSize; i++) {
0274           photo_->push_back(
0275               HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
0276         }
0277       }
0278     } else {
0279       hadBranch_->SetAddress(&photon_);
0280       hadBranch_->GetEntry(nrc);
0281     }
0282   } else {
0283     if (newForm_) {
0284       if (!v3version_) {
0285         emBranch_->SetAddress(&photo_);
0286         emBranch_->GetEntry(nrc);
0287       } else {
0288         std::vector<float> t;
0289         std::vector<float>* tp = &t;
0290         emBranch_->SetAddress(&tp);
0291         emBranch_->GetEntry(nrc);
0292         unsigned int tSize = t.size() / 5;
0293         photo_->reserve(tSize);
0294         for (unsigned int i = 0; i < tSize; i++) {
0295           photo_->push_back(
0296               HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
0297         }
0298       }
0299     } else {
0300       emBranch_->SetAddress(&photon_);
0301       emBranch_->GetEntry(nrc);
0302     }
0303   }
0304   if (verbose_) {
0305     int nPhoton = (newForm_) ? photo_->size() : photon_.size();
0306     edm::LogVerbatim("HFShower") << "getRecord: Record " << record << " of type " << type << " with " << nPhoton
0307                                  << " photons";
0308     for (int j = 0; j < nPhoton; j++)
0309       if (newForm_)
0310         edm::LogVerbatim("HFShower") << "Photon " << j << " " << photo_->at(j);
0311       else
0312         edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon_[j];
0313   }
0314 }
0315 
0316 void HFShowerLibraryAnalyzer::loadEventInfo(TBranch* branch) {
0317   if (branch) {
0318     std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
0319     branch->SetAddress(&eventInfoCollection);
0320     branch->GetEntry(0);
0321     edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
0322                                  << eventInfoCollection.size() << " records";
0323     totEvents_ = eventInfoCollection[0].totalEvents();
0324     nMomBin_ = eventInfoCollection[0].numberOfBins();
0325     evtPerBin_ = eventInfoCollection[0].eventsPerBin();
0326     libVers_ = eventInfoCollection[0].showerLibraryVersion();
0327     listVersion_ = eventInfoCollection[0].physListVersion();
0328     pmom_ = eventInfoCollection[0].energyBins();
0329   } else {
0330     edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired numbers";
0331     nMomBin_ = 16;
0332     totEvents_ = nMomBin_ * evtPerBin_;
0333     libVers_ = 1.1;
0334     listVersion_ = 3.6;
0335     pmom_ = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
0336   }
0337   for (int i = 0; i < nMomBin_; i++)
0338     pmom_[i] *= CLHEP::GeV;
0339 }
0340 
0341 //define this as a plug-in
0342 DEFINE_FWK_MODULE(HFShowerLibraryAnalyzer);