EcalSimpleProducer.cc

Go to the documentation of this file.

#include "SimCalorimetry/EcalElectronicsEmulation/interface/EcalSimpleProducer.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDigi/interface/EBDataFrame.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
#include "TFormula.h"
#include <iostream>
#include <memory>

#include <string>

using namespace std;
using namespace edm;

void EcalSimpleProducer::produce(edm::Event &evt, const edm::EventSetup &) {
  const int ievt = evt.id().event();
  if (formula_.get() != nullptr) {
    unique_ptr<EBDigiCollection> digis(new EBDigiCollection);

    digis->reserve(170 * 360);

    const int nSamples = digis->stride();
    for (int iEta0 = 0; iEta0 < 170; ++iEta0) {
      for (int iPhi0 = 0; iPhi0 < 360; ++iPhi0) {
        int iEta1 = cIndex2iEta(iEta0);
        int iPhi = cIndex2iPhi(iPhi0);
        if (verbose_)
          cout << "(" << iEta0 << "," << iPhi0 << "): ";
        digis->push_back(EBDetId(iEta1, iPhi));
        DataFrame dframe(digis->back());

        for (int t = 0; t < nSamples; ++t) {
          uint16_t encodedAdc = (uint16_t)formula_->Eval(iEta0, iPhi0, ievt - 1, t);
          if (verbose_)
            cout << encodedAdc << ((t < (nSamples - 1)) ? "\t" : "\n");
          dframe[t] = encodedAdc;
        }
      }
    }
    evt.put(std::move(digis));
    // puts an empty digi collecion for endcap:
    evt.put(std::make_unique<EEDigiCollection>());
  }
  if (tpFormula_.get() != nullptr) {
    unique_ptr<EcalTrigPrimDigiCollection> tps = std::make_unique<EcalTrigPrimDigiCollection>();
    tps->reserve(56 * 72);
    const int nSamples = 5;
    for (int iTtEta0 = 0; iTtEta0 < 56; ++iTtEta0) {
      for (int iTtPhi0 = 0; iTtPhi0 < 72; ++iTtPhi0) {
        int iTtEta1 = cIndex2iTtEta(iTtEta0);
        int iTtPhi = cIndex2iTtPhi(iTtPhi0);

        if (verbose_)
          cout << "(" << iTtEta0 << "," << iTtPhi0 << "): ";
        int zside = iTtEta1 < 0 ? -1 : 1;
        EcalTriggerPrimitiveDigi tpframe(EcalTrigTowerDetId(zside, EcalTriggerTower, abs(iTtEta1), iTtPhi));

        tpframe.setSize(nSamples);

        if (verbose_)
          cout << "TP: ";
        for (int t = 0; t < nSamples; ++t) {
          uint16_t encodedTp = (uint16_t)tpFormula_->Eval(iTtEta0, iTtPhi0, ievt - 1, t);

          if (verbose_)
            cout << "TP(" << iTtEta0 << "," << iTtPhi0 << ") = " << encodedTp << ((t < (nSamples - 1)) ? "\t" : "\n");
          tpframe.setSample(t, EcalTriggerPrimitiveSample(encodedTp));
        }
        tps->push_back(tpframe);
      }
    }
    evt.put(std::move(tps));
  }
  if (simHitFormula_.get() != nullptr) {  // generation of barrel sim hits
    unique_ptr<PCaloHitContainer> hits = std::make_unique<PCaloHitContainer>();
    for (int iEta0 = 0; iEta0 < 170; ++iEta0) {
      for (int iPhi0 = 0; iPhi0 < 360; ++iPhi0) {
        int iEta1 = cIndex2iEta(iEta0);
        int iPhi = cIndex2iPhi(iPhi0);
        if (verbose_)
          cout << "(" << iEta0 << "," << iPhi0 << "): ";

        double em = simHitFormula_->Eval(iEta0, iPhi0, ievt - 1);
        double eh = 0.;
        double t = 0.;
        const PCaloHit hit(EBDetId(iEta1, iPhi).rawId(), em, eh, t, 0, 0);
        hits->push_back(hit);
      }
    }
    evt.put(std::move(hits), "EcalHitsEB");
    // puts an empty digi collecion for endcap:
    evt.put(std::make_unique<PCaloHitContainer>(), "EcalHitsEE");
  }
}

EcalSimpleProducer::EcalSimpleProducer(const edm::ParameterSet &pset) : EDProducer() {
  string formula = pset.getParameter<string>("formula");
  string tpFormula = pset.getParameter<string>("tpFormula");
  string simHitFormula = pset.getParameter<string>("simHitFormula");

  verbose_ = pset.getUntrackedParameter<bool>("verbose", false);
  //  replaceAll(formula, "itt0",
  //  "((((ieta0<85)*(84-ieta0)+(ieta0>=85)*(ieta0-85))/5-18)*4+((iphi0/5+2)%4))");
  replaceAll(formula, "ebm", "(ieta0<85)");
  replaceAll(formula, "ebp", "(ieta0>84)");
  replaceAll(formula, "ieta0", "x");
  replaceAll(formula, "iphi0", "y");
  replaceAll(formula, "ievt0", "z");
  replaceAll(formula, "isample0", "t");
  //  cout << "----------> " << formula << endl;

  replaceAll(tpFormula, "itt0", "((ieta0<28)*(27-ieta0)+(ieta0>=28)*(ieta0-28))*4+(iphi0+2)%4");
  replaceAll(tpFormula, "eb", "(ieta0>10 && ieta0<45)");
  replaceAll(tpFormula, "ebm", "(ieta0>10 && ieta0<28)");
  replaceAll(tpFormula, "ebp", "(ieta0>27 && ieta0<45)");
  replaceAll(tpFormula, "ee", "(ieta0<11 || ieta0>44)");
  replaceAll(tpFormula, "eem", "(ieta0<11)");
  replaceAll(tpFormula, "eep", "(ieta0>44)");
  replaceAll(tpFormula, "ieta0", "x");
  replaceAll(tpFormula, "iphi0", "y");
  replaceAll(tpFormula, "ievt0", "z");
  replaceAll(tpFormula, "isample0", "t");
  // cout << "----------> " << tpFormula << endl;

  // replaceAll(simHitormula, "itt0",
  // "((((ieta0<85)*(84-ieta0)+(ieta0>=85)*(ieta0-85))/5-18)*4+((iphi0/5+2)%4))");
  replaceAll(simHitFormula, "ebm", "(ieta0<85)");
  replaceAll(simHitFormula, "ebp", "(ieta0>84)");
  replaceAll(simHitFormula, "ieta0", "x");
  replaceAll(simHitFormula, "iphi0", "y");
  replaceAll(simHitFormula, "ievt0", "z");

  if (!formula.empty()) {
    formula_ = std::make_unique<TFormula>("f", formula.c_str());
    Int_t err = formula_->Compile();
    if (err != 0) {
      throw cms::Exception("Error in EcalSimpleProducer 'formula' config.");
    }
    produces<EBDigiCollection>();
    produces<EEDigiCollection>();
  }
  if (!tpFormula.empty()) {
    tpFormula_ = std::make_unique<TFormula>("f", tpFormula.c_str());
    Int_t err = tpFormula_->Compile();
    if (err != 0) {
      throw cms::Exception("Error in EcalSimpleProducer 'tpFormula' config.");
    }
    produces<EcalTrigPrimDigiCollection>();
  }
  if (!simHitFormula.empty()) {
    simHitFormula_ = std::make_unique<TFormula>("f", simHitFormula.c_str());
    Int_t err = simHitFormula_->Compile();
    if (err != 0) {
      throw cms::Exception(
          "Error in EcalSimpleProducer "
          "'simHitFormula' config.");
    }
    produces<edm::PCaloHitContainer>("EcalHitsEB");
    produces<edm::PCaloHitContainer>("EcalHitsEE");
  }
}

void EcalSimpleProducer::replaceAll(std::string &s, const std::string &from, const std::string &to) const {
  string::size_type pos = 0;
  //  cout << "replaceAll(" << s << "," << from << "," << to << ")\n";
  while ((pos = s.find(from, pos)) != string::npos) {
    // cout << "replace(" << pos << "," << from.size() << "," << to << ")\n";
    s.replace(pos, from.size(), to);
    //   cout << " -> " << s << "\n";
  }
}