GFlash/TB/TreeProducerCalibSimul.cc

0001 // framework
0002 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011
0012 // for reconstruction
0013 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0014 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0015 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0016 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0017 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRecInfo.h"
0018 #include "TBDataFormats/EcalTBObjects/interface/EcalTBTDCRecInfo.h"
0019 #include "TBDataFormats/EcalTBObjects/interface/EcalTBEventHeader.h"
0020
0021 // geometry
0022 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0023 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0024 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0025
0026 // my include files
0027 #include "SimG4Core/GFlash/TB/TreeMatrixCalib.h"
0028
0029 // root includes
0030 #include "TROOT.h"
0031 #include "TSystem.h"
0032 #include "TFile.h"
0033 #include "TTree.h"
0034 #include "TF1.h"
0035 #include "TH1.h"
0036 #include "TH2.h"
0037 #include "TGraph.h"
0038 #include "TStyle.h"
0039 #include "TCanvas.h"
0040 #include "TSelector.h"
0041 #include "TApplication.h"
0042
0043 // c++ includes
0044 #include <string>
0045 #include <stdio.h>
0046 #include <sstream>
0047 #include <iostream>
0048 #include <unistd.h>
0049 #include <fstream>
0050 #include <math.h>
0051 #include <memory>
0052 #include <stdexcept>
0053
0054 class TreeProducerCalibSimul : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0055 public:
0056   explicit TreeProducerCalibSimul(const edm::ParameterSet&);
0057   ~TreeProducerCalibSimul() override = default;
0058
0059   void analyze(const edm::Event&, const edm::EventSetup&) override;
0060   void beginJob() override;
0061   void endJob() override;
0062
0063 private:
0064   std::string rootfile_;
0065   std::string txtfile_;
0066   std::string EBRecHitCollection_;
0067   std::string RecHitProducer_;
0068   std::string hodoRecInfoCollection_;
0069   std::string hodoRecInfoProducer_;
0070   std::string tdcRecInfoCollection_;
0071   std::string tdcRecInfoProducer_;
0072   std::string eventHeaderCollection_;
0073   std::string eventHeaderProducer_;
0074   double posCluster_;
0075
0076   std::unique_ptr<TreeMatrixCalib> myTree_;
0077
0078   edm::EDGetTokenT<EBRecHitCollection> tokEBRecHit_;
0079   edm::EDGetTokenT<EcalTBHodoscopeRecInfo> tokEcalHodo_;
0080   edm::EDGetTokenT<EcalTBTDCRecInfo> tokEcalTDC_;
0081   edm::EDGetTokenT<EcalTBEventHeader> tokEventHeader_;
0082
0083   int xtalInBeam_;
0084   int tot_events_;
0085   int tot_events_ok_;
0086   int noHits_;
0087   int noHodo_;
0088   int noTdc_;
0089   int noHeader_;
0090 };
0091
0092 // -------------------------------------------------
0093 // contructor
0094 TreeProducerCalibSimul::TreeProducerCalibSimul(const edm::ParameterSet& iConfig) {
0095   usesResource(TFileService::kSharedResource);
0096
0097   // now do what ever initialization is needed
0098   xtalInBeam_ = iConfig.getUntrackedParameter<int>("xtalInBeam", -1000);
0099   rootfile_ = iConfig.getUntrackedParameter<std::string>("rootfile", "mySimMatrixTree.root");
0100   txtfile_ = iConfig.getUntrackedParameter<std::string>("txtfile", "mySimMatrixTree.txt");
0101   EBRecHitCollection_ = iConfig.getParameter<std::string>("EBRecHitCollection");
0102   RecHitProducer_ = iConfig.getParameter<std::string>("RecHitProducer");
0103   hodoRecInfoCollection_ = iConfig.getParameter<std::string>("hodoRecInfoCollection");
0104   hodoRecInfoProducer_ = iConfig.getParameter<std::string>("hodoRecInfoProducer");
0105   tdcRecInfoCollection_ = iConfig.getParameter<std::string>("tdcRecInfoCollection");
0106   tdcRecInfoProducer_ = iConfig.getParameter<std::string>("tdcRecInfoProducer");
0107   eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
0108   eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
0109
0110   // summary
0111   edm::LogVerbatim("GFlash") << "\nConstructor\n\nTreeProducerCalibSimul\nxtal in beam = " << xtalInBeam_;
0112   edm::LogVerbatim("GFlash") << "Fetching hitCollection: " << EBRecHitCollection_.c_str() << " prod by "
0113                              << RecHitProducer_.c_str();
0114   edm::LogVerbatim("GFlash") << "Fetching hodoCollection: " << hodoRecInfoCollection_.c_str() << " prod by "
0115                              << hodoRecInfoProducer_.c_str();
0116   edm::LogVerbatim("GFlash") << "Fetching tdcCollection: " << tdcRecInfoCollection_.c_str() << " prod by "
0117                              << tdcRecInfoProducer_.c_str();
0118   edm::LogVerbatim("GFlash") << "Fetching evHeaCollection: " << eventHeaderCollection_.c_str() << " prod by "
0119                              << eventHeaderProducer_.c_str() << "\n";
0120
0121   tokEBRecHit_ = consumes<EBRecHitCollection>(edm::InputTag(RecHitProducer_, EBRecHitCollection_));
0122   tokEcalHodo_ = consumes<EcalTBHodoscopeRecInfo>(edm::InputTag(hodoRecInfoProducer_, hodoRecInfoCollection_));
0123   tokEcalTDC_ = consumes<EcalTBTDCRecInfo>(edm::InputTag(tdcRecInfoProducer_, tdcRecInfoCollection_));
0124   tokEventHeader_ = consumes<EcalTBEventHeader>(edm::InputTag(eventHeaderProducer_));
0125 }
0126
0127 // ------------------------------------------------------
0128 // initializations
0129 void TreeProducerCalibSimul::beginJob() {
0130   edm::LogVerbatim("GFlash") << "\nBeginJob\n";
0131
0132   // tree
0133   myTree_ = std::make_unique<TreeMatrixCalib>(rootfile_.c_str());
0134
0135   // counters
0136   tot_events_ = 0;
0137   tot_events_ok_ = 0;
0138   noHits_ = 0;
0139   noHodo_ = 0;
0140   noTdc_ = 0;
0141   noHeader_ = 0;
0142 }
0143
0144 // -------------------------------------------
0145 // finalizing
0146 void TreeProducerCalibSimul::endJob() {
0147   edm::LogVerbatim("GFlash") << "\nEndJob\n";
0148
0149   std::ofstream MyOut(txtfile_.c_str());
0150   MyOut << "total events: " << tot_events_ << std::endl;
0151   MyOut << "events skipped because of no hits: " << noHits_ << std::endl;
0152   MyOut << "events skipped because of no hodos: " << noHodo_ << std::endl;
0153   MyOut << "events skipped because of no tdc: " << noTdc_ << std::endl;
0154   MyOut << "events skipped because of no header: " << noHeader_ << std::endl;
0155   MyOut << "total OK events (passing the basic selection): " << tot_events_ok_ << std::endl;
0156   MyOut.close();
0157 }
0158
0159 // -----------------------------------------------
0160 // my analysis
0161 void TreeProducerCalibSimul::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0162   // counting events
0163   ++tot_events_;
0164
0165   if (tot_events_ % 5000 == 0)
0166     edm::LogVerbatim("GFlash") << "event " << tot_events_;
0167
0168   // ---------------------------------------------------------------------
0169   // taking what I need: hits
0170   const EBRecHitCollection* EBRecHits = &iEvent.get(tokEBRecHit_);
0171
0172   // taking what I need: hodoscopes
0173   const EcalTBHodoscopeRecInfo* recHodo = &iEvent.get(tokEcalHodo_);
0174
0175   // taking what I need: tdc
0176   const EcalTBTDCRecInfo* recTDC = &iEvent.get(tokEcalTDC_);
0177
0178   // taking what I need: event header
0179   const EcalTBEventHeader* evtHeader = &iEvent.get(tokEventHeader_);
0180
0181   // checking everything is there and fine
0182   if ((!EBRecHits) || (EBRecHits->size() == 0)) {
0183     ++noHits_;
0184     return;
0185   }
0186   if (!recTDC) {
0187     ++noTdc_;
0188     return;
0189   }
0190   if (!recHodo) {
0191     ++noHodo_;
0192     return;
0193   }
0194   if (!evtHeader) {
0195     ++noHeader_;
0196     return;
0197   }
0198   ++tot_events_ok_;
0199
0200   // ---------------------------------------------------------------------
0201   // info on the event
0202   int run = -999;
0203   int tbm = -999;
0204   int event = evtHeader->eventNumber();
0205
0206   // ---------------------------------------------------------------------
0207   // xtal-in-beam
0208   int nomXtalInBeam = -999;
0209   int nextXtalInBeam = -999;
0210
0211   EBDetId xtalInBeamId(1, xtalInBeam_, EBDetId::SMCRYSTALMODE);
0212   if (xtalInBeamId == EBDetId(0)) {
0213     return;
0214   }
0215   int mySupCry = xtalInBeamId.ic();
0216   int mySupEta = xtalInBeamId.ieta();
0217   int mySupPhi = xtalInBeamId.iphi();
0218
0219   // ---------------------------------------------------------------------
0220   // hodoscope information
0221   double x = recHodo->posX();
0222   double y = recHodo->posY();
0223   double sx = recHodo->slopeX();
0224   double sy = recHodo->slopeY();
0225   double qx = recHodo->qualX();
0226   double qy = recHodo->qualY();
0227
0228   // ---------------------------------------------------------------------
0229   // tdc information
0230   double tdcOffset = recTDC->offset();
0231
0232   // ---------------------------------------------------------------------
0233   // Find EBDetId in a 7x7 Matrix
0234   EBDetId Xtals7x7[49];
0235   double energy[49];
0236   int crystal[49];
0237   int allMatrix = 1;
0238   for (unsigned int icry = 0; icry < 49; icry++) {
0239     unsigned int row = icry / 7;
0240     unsigned int column = icry % 7;
0241     Xtals7x7[icry] = EBDetId(xtalInBeamId.ieta() + column - 3, xtalInBeamId.iphi() + row - 3, EBDetId::ETAPHIMODE);
0242
0243     if (Xtals7x7[icry].ism() == 1) {
0244       energy[icry] = EBRecHits->find(Xtals7x7[icry])->energy();
0245       crystal[icry] = Xtals7x7[icry].ic();
0246     } else {
0247       energy[icry] = -100.;
0248       crystal[icry] = -100;
0249       allMatrix = 0;
0250     }
0251   }
0252
0253   // ---------------------------------------------------------------------
0254   // Looking for the max energy crystal
0255   double maxEne = -999.;
0256   int maxEneCry = 9999;
0257   for (int ii = 0; ii < 49; ii++) {
0258     if (energy[ii] > maxEne) {
0259       maxEne = energy[ii];
0260       maxEneCry = crystal[ii];
0261     }
0262   }
0263
0264   // Position reconstruction - skipped here
0265   double Xcal = -999.;
0266   double Ycal = -999.;
0267
0268   // filling the tree
0269   myTree_->fillInfo(run,
0270                     event,
0271                     mySupCry,
0272                     maxEneCry,
0273                     nomXtalInBeam,
0274                     nextXtalInBeam,
0275                     mySupEta,
0276                     mySupPhi,
0277                     tbm,
0278                     x,
0279                     y,
0280                     Xcal,
0281                     Ycal,
0282                     sx,
0283                     sy,
0284                     qx,
0285                     qy,
0286                     tdcOffset,
0287                     allMatrix,
0288                     energy,
0289                     crystal);
0290   myTree_->store();
0291 }
0292
0293 #include "FWCore/PluginManager/interface/ModuleDef.h"
0294 #include "FWCore/Framework/interface/MakerMacros.h"
0295
0296 //define this as a plug-in
0297
0298 DEFINE_FWK_MODULE(TreeProducerCalibSimul);