EcalElectronicsEmulation/src/EcalSimpleProducer.cc

0001 #include "SimCalorimetry/EcalElectronicsEmulation/interface/EcalSimpleProducer.h"
0002 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0003 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0004 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0009 #include "TFormula.h"
0010 #include <iostream>
0011 #include <memory>
0012
0013 #include <string>
0014
0015 using namespace std;
0016 using namespace edm;
0017
0018 void EcalSimpleProducer::produce(edm::Event &evt, const edm::EventSetup &) {
0019   const int ievt = evt.id().event();
0020   if (formula_.get() != nullptr) {
0021     unique_ptr<EBDigiCollection> digis(new EBDigiCollection);
0022
0023     digis->reserve(170 * 360);
0024
0025     const int nSamples = digis->stride();
0026     for (int iEta0 = 0; iEta0 < 170; ++iEta0) {
0027       for (int iPhi0 = 0; iPhi0 < 360; ++iPhi0) {
0028         int iEta1 = cIndex2iEta(iEta0);
0029         int iPhi = cIndex2iPhi(iPhi0);
0030         if (verbose_)
0031           cout << "(" << iEta0 << "," << iPhi0 << "): ";
0032         digis->push_back(EBDetId(iEta1, iPhi));
0033         DataFrame dframe(digis->back());
0034
0035         for (int t = 0; t < nSamples; ++t) {
0036           uint16_t encodedAdc = (uint16_t)formula_->Eval(iEta0, iPhi0, ievt - 1, t);
0037           if (verbose_)
0038             cout << encodedAdc << ((t < (nSamples - 1)) ? "\t" : "\n");
0039           dframe[t] = encodedAdc;
0040         }
0041       }
0042     }
0043     evt.put(std::move(digis));
0044     // puts an empty digi collecion for endcap:
0045     evt.put(std::make_unique<EEDigiCollection>());
0046   }
0047   if (tpFormula_.get() != nullptr) {
0048     unique_ptr<EcalTrigPrimDigiCollection> tps = std::make_unique<EcalTrigPrimDigiCollection>();
0049     tps->reserve(56 * 72);
0050     const int nSamples = 5;
0051     for (int iTtEta0 = 0; iTtEta0 < 56; ++iTtEta0) {
0052       for (int iTtPhi0 = 0; iTtPhi0 < 72; ++iTtPhi0) {
0053         int iTtEta1 = cIndex2iTtEta(iTtEta0);
0054         int iTtPhi = cIndex2iTtPhi(iTtPhi0);
0055
0056         if (verbose_)
0057           cout << "(" << iTtEta0 << "," << iTtPhi0 << "): ";
0058         int zside = iTtEta1 < 0 ? -1 : 1;
0059         EcalTriggerPrimitiveDigi tpframe(EcalTrigTowerDetId(zside, EcalTriggerTower, abs(iTtEta1), iTtPhi));
0060
0061         tpframe.setSize(nSamples);
0062
0063         if (verbose_)
0064           cout << "TP: ";
0065         for (int t = 0; t < nSamples; ++t) {
0066           uint16_t encodedTp = (uint16_t)tpFormula_->Eval(iTtEta0, iTtPhi0, ievt - 1, t);
0067
0068           if (verbose_)
0069             cout << "TP(" << iTtEta0 << "," << iTtPhi0 << ") = " << encodedTp << ((t < (nSamples - 1)) ? "\t" : "\n");
0070           tpframe.setSample(t, EcalTriggerPrimitiveSample(encodedTp));
0071         }
0072         tps->push_back(tpframe);
0073       }
0074     }
0075     evt.put(std::move(tps));
0076   }
0077   if (simHitFormula_.get() != nullptr) {  // generation of barrel sim hits
0078     unique_ptr<PCaloHitContainer> hits = std::make_unique<PCaloHitContainer>();
0079     for (int iEta0 = 0; iEta0 < 170; ++iEta0) {
0080       for (int iPhi0 = 0; iPhi0 < 360; ++iPhi0) {
0081         int iEta1 = cIndex2iEta(iEta0);
0082         int iPhi = cIndex2iPhi(iPhi0);
0083         if (verbose_)
0084           cout << "(" << iEta0 << "," << iPhi0 << "): ";
0085
0086         double em = simHitFormula_->Eval(iEta0, iPhi0, ievt - 1);
0087         double eh = 0.;
0088         double t = 0.;
0089         const PCaloHit hit(EBDetId(iEta1, iPhi).rawId(), em, eh, t, 0, 0);
0090         hits->push_back(hit);
0091       }
0092     }
0093     evt.put(std::move(hits), "EcalHitsEB");
0094     // puts an empty digi collecion for endcap:
0095     evt.put(std::make_unique<PCaloHitContainer>(), "EcalHitsEE");
0096   }
0097 }
0098
0099 EcalSimpleProducer::EcalSimpleProducer(const edm::ParameterSet &pset) : EDProducer() {
0100   string formula = pset.getParameter<string>("formula");
0101   string tpFormula = pset.getParameter<string>("tpFormula");
0102   string simHitFormula = pset.getParameter<string>("simHitFormula");
0103
0104   verbose_ = pset.getUntrackedParameter<bool>("verbose", false);
0105   //  replaceAll(formula, "itt0",
0106   //  "((((ieta0<85)*(84-ieta0)+(ieta0>=85)*(ieta0-85))/5-18)*4+((iphi0/5+2)%4))");
0107   replaceAll(formula, "ebm", "(ieta0<85)");
0108   replaceAll(formula, "ebp", "(ieta0>84)");
0109   replaceAll(formula, "ieta0", "x");
0110   replaceAll(formula, "iphi0", "y");
0111   replaceAll(formula, "ievt0", "z");
0112   replaceAll(formula, "isample0", "t");
0113   //  cout << "----------> " << formula << endl;
0114
0115   replaceAll(tpFormula, "itt0", "((ieta0<28)*(27-ieta0)+(ieta0>=28)*(ieta0-28))*4+(iphi0+2)%4");
0116   replaceAll(tpFormula, "eb", "(ieta0>10 && ieta0<45)");
0117   replaceAll(tpFormula, "ebm", "(ieta0>10 && ieta0<28)");
0118   replaceAll(tpFormula, "ebp", "(ieta0>27 && ieta0<45)");
0119   replaceAll(tpFormula, "ee", "(ieta0<11 || ieta0>44)");
0120   replaceAll(tpFormula, "eem", "(ieta0<11)");
0121   replaceAll(tpFormula, "eep", "(ieta0>44)");
0122   replaceAll(tpFormula, "ieta0", "x");
0123   replaceAll(tpFormula, "iphi0", "y");
0124   replaceAll(tpFormula, "ievt0", "z");
0125   replaceAll(tpFormula, "isample0", "t");
0126   // cout << "----------> " << tpFormula << endl;
0127
0128   // replaceAll(simHitormula, "itt0",
0129   // "((((ieta0<85)*(84-ieta0)+(ieta0>=85)*(ieta0-85))/5-18)*4+((iphi0/5+2)%4))");
0130   replaceAll(simHitFormula, "ebm", "(ieta0<85)");
0131   replaceAll(simHitFormula, "ebp", "(ieta0>84)");
0132   replaceAll(simHitFormula, "ieta0", "x");
0133   replaceAll(simHitFormula, "iphi0", "y");
0134   replaceAll(simHitFormula, "ievt0", "z");
0135
0136   if (!formula.empty()) {
0137     formula_ = std::make_unique<TFormula>("f", formula.c_str());
0138     Int_t err = formula_->Compile();
0139     if (err != 0) {
0140       throw cms::Exception("Error in EcalSimpleProducer 'formula' config.");
0141     }
0142     produces<EBDigiCollection>();
0143     produces<EEDigiCollection>();
0144   }
0145   if (!tpFormula.empty()) {
0146     tpFormula_ = std::make_unique<TFormula>("f", tpFormula.c_str());
0147     Int_t err = tpFormula_->Compile();
0148     if (err != 0) {
0149       throw cms::Exception("Error in EcalSimpleProducer 'tpFormula' config.");
0150     }
0151     produces<EcalTrigPrimDigiCollection>();
0152   }
0153   if (!simHitFormula.empty()) {
0154     simHitFormula_ = std::make_unique<TFormula>("f", simHitFormula.c_str());
0155     Int_t err = simHitFormula_->Compile();
0156     if (err != 0) {
0157       throw cms::Exception(
0158           "Error in EcalSimpleProducer "
0159           "'simHitFormula' config.");
0160     }
0161     produces<edm::PCaloHitContainer>("EcalHitsEB");
0162     produces<edm::PCaloHitContainer>("EcalHitsEE");
0163   }
0164 }
0165
0166 void EcalSimpleProducer::replaceAll(std::string &s, const std::string &from, const std::string &to) const {
0167   string::size_type pos = 0;
0168   //  cout << "replaceAll(" << s << "," << from << "," << to << ")\n";
0169   while ((pos = s.find(from, pos)) != string::npos) {
0170     // cout << "replace(" << pos << "," << from.size() << "," << to << ")\n";
0171     s.replace(pos, from.size(), to);
0172     //   cout << " -> " << s << "\n";
0173   }
0174 }