EcalTrigPrimFunctionalAlgo.h

Go to the documentation of this file.

#ifndef EcalTrigPrimFunctionalAlgo_h
#define EcalTrigPrimFunctionalAlgo_h

#include <iostream>
#include <sys/time.h>
#include <vector>

#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"

#include "SimCalorimetry/EcalTrigPrimAlgos/interface/EcalFenixStrip.h"
#include "SimCalorimetry/EcalTrigPrimAlgos/interface/EcalFenixTcp.h"

#include "DataFormats/Common/interface/SortedCollection.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <map>
#include <utility>

class EcalTrigTowerDetId;
class ETPCoherenceTest;
class EcalTriggerPrimitiveSample;
class CaloSubdetectorGeometry;
class EBDataFrame;
class EEDataFrame;
class EcalElectronicsMapping;

class EcalTrigPrimFunctionalAlgo {
public:
  explicit EcalTrigPrimFunctionalAlgo(const edm::EventSetup &setup,
                                      int binofmax, bool tcpFormat,
                                      bool barrelOnly, bool debug, bool famos);

  virtual ~EcalTrigPrimFunctionalAlgo();

  void run(const edm::EventSetup &, const EBDigiCollection *col,
           EcalTrigPrimDigiCollection &result,
           EcalTrigPrimDigiCollection &resultTcp);
  void run(const edm::EventSetup &, const EEDigiCollection *col,
           EcalTrigPrimDigiCollection &result,
           EcalTrigPrimDigiCollection &resultTcp);
  void run_part1_EB(EBDigiCollection const *col);
  void run_part1_EE(EEDigiCollection const *col);
  template <class Coll>
  void run_part2(
      const edm::EventSetup &, Coll const *col,
      std::vector<std::vector<std::pair<int, std::vector<typename Coll::Digi>>>>
          &towerMap,
      EcalTrigPrimDigiCollection &result,
      EcalTrigPrimDigiCollection &resultTcp);

  void setPointers(const EcalTPGLinearizationConst *ecaltpLin,
                   const EcalTPGPedestals *ecaltpPed,
                   const EcalTPGSlidingWindow *ecaltpgSlidW,
                   const EcalTPGWeightIdMap *ecaltpgWeightMap,
                   const EcalTPGWeightGroup *ecaltpgWeightGroup,
                   const EcalTPGFineGrainStripEE *ecaltpgFgStripEE,
                   const EcalTPGCrystalStatus *ecaltpgBadX,
                   const EcalTPGStripStatus *ecaltpgStripStatus) {
    estrip_->setPointers(ecaltpPed, ecaltpLin, ecaltpgWeightMap,
                         ecaltpgWeightGroup, ecaltpgSlidW, ecaltpgFgStripEE,
                         ecaltpgBadX, ecaltpgStripStatus);
  }
  void setPointers2(const EcalTPGFineGrainEBGroup *ecaltpgFgEBGroup,
                    const EcalTPGLutGroup *ecaltpgLutGroup,
                    const EcalTPGLutIdMap *ecaltpgLut,
                    const EcalTPGFineGrainEBIdMap *ecaltpgFineGrainEB,
                    const EcalTPGFineGrainTowerEE *ecaltpgFineGrainTowerEE,
                    const EcalTPGTowerStatus *ecaltpgBadTT,
                    const EcalTPGSpike *ecaltpgSpike) {

    etcp_->setPointers(ecaltpgFgEBGroup, ecaltpgLutGroup, ecaltpgLut,
                       ecaltpgFineGrainEB, ecaltpgFineGrainTowerEE,
                       ecaltpgBadTT, ecaltpgSpike);
  }

private:
  void init(const edm::EventSetup &);
  template <class T>
  void initStructures(
      std::vector<std::vector<std::pair<int, std::vector<T>>>> &towMap);
  template <class T>
  void
  clean(std::vector<std::vector<std::pair<int, std::vector<T>>>> &towerMap);
  template <class Coll>
  void fillMap(
      Coll const *col,
      std::vector<std::vector<std::pair<int, std::vector<typename Coll::Digi>>>>
          &towerMap);
  int findStripNr(const EBDetId &id);
  int findStripNr(const EEDetId &id);

  // FIXME: temporary until hashedIndex works alsom for endcap
  int getIndex(const EBDigiCollection *, EcalTrigTowerDetId &id) {
    return id.hashedIndex();
  }
  // mind that eta is continuous between barrel+endcap
  int getIndex(const EEDigiCollection *, EcalTrigTowerDetId &id) {
    int ind = (id.ietaAbs() - 18) * 72 + id.iphi();
    if (id.zside() < 0)
      ind += 792;
    return ind;
  }

  EcalFenixStrip *estrip_;
  EcalFenixTcp *etcp_;

  edm::ESHandle<EcalTrigTowerConstituentsMap> eTTmap_;
  const CaloSubdetectorGeometry *theEndcapGeometry;
  const EcalElectronicsMapping *theMapping_;

  float threshold;

  int binOfMaximum_;
  int maxNrSamples_;

  bool tcpFormat_;
  bool barrelOnly_;
  bool debug_;
  bool famos_;

  static const unsigned int
      nrSamples_; // nr samples to write, should not be changed since by
                  // convention the size means that it is coming from simulation
  static const unsigned int
      maxNrSamplesOut_; // to be placed in the intermediate samples
  static const unsigned int
      maxNrTowers_; // would be better to get from somewhere..
  static const unsigned int maxNrTPs_; // would be better to get from
                                       // somewhere..

  int nrTowers_; // nr of towers found by fillmap method

  // data structures kept during the whole run
  std::vector<std::vector<int>> striptp_;
  std::vector<std::vector<std::pair<int, std::vector<EBDataFrame>>>>
      towerMapEB_;
  std::vector<std::vector<std::pair<int, std::vector<EEDataFrame>>>>
      towerMapEE_;
  std::vector<std::pair<int, EcalTrigTowerDetId>> hitTowers_;
  std::vector<EcalTriggerPrimitiveSample> towtp_;
  std::vector<EcalTriggerPrimitiveSample> towtp2_;

  enum { nbMaxStrips_ = 5 };
  enum { nbMaxXtals_ = 5 };
};

//=================================== implementations
//=============================================

template <class Coll>
void EcalTrigPrimFunctionalAlgo::run_part2(
    const edm::EventSetup &setup, Coll const *col,
    std::vector<std::vector<std::pair<int, std::vector<typename Coll::Digi>>>>
        &towerMap,
    EcalTrigPrimDigiCollection &result, EcalTrigPrimDigiCollection &resultTcp) {
  typedef typename Coll::Digi Digi;

  // prepare writing of TP-s

  int firstSample = binOfMaximum_ - 1 - nrSamples_ / 2;
  int lastSample = binOfMaximum_ - 1 + nrSamples_ / 2;
  int nrTP = 0;
  std::vector<typename Coll::Digi> dummy;
  EcalTriggerPrimitiveDigi tptow[2];
  EcalTriggerPrimitiveDigi tptowTcp[2];

  estrip_->getFGVB()->setbadStripMissing(false);

  for (int itow = 0; itow < nrTowers_; ++itow) {

    int index = hitTowers_[itow].first;
    const EcalTrigTowerDetId &thisTower = hitTowers_[itow].second;

    // loop over all strips assigned to this trigger tower
    int nstr = 0;
    for (unsigned int i = 0; i < towerMap[itow].size(); ++i) {
      std::vector<Digi> &df =
          (towerMap[index])[i].second; // vector of dataframes for this strip,
                                       // size; nr of crystals/strip

      if ((towerMap[index])[i].first > 0) {
        estrip_->process(setup, df, (towerMap[index])[i].first,
                         striptp_[nstr++]);
      }
    } // loop over strips in one tower

    bool isInInnerRings = false;
    if (thisTower.subDet() == EcalEndcap &&
        (thisTower.ietaAbs() == 27 || thisTower.ietaAbs() == 28))
      isInInnerRings = true;
    etcp_->process(setup, dummy, striptp_, nstr, towtp_, towtp2_,
                   isInInnerRings, thisTower);

    // prepare TP-s
    // special treatment for 2 inner endcap rings
    int nrTowers;
    if (isInInnerRings) {
      nrTowers = 2;
      int phi = 2 * ((thisTower.iphi() - 1) / 2);
      tptow[0] = EcalTriggerPrimitiveDigi(EcalTrigTowerDetId(
          thisTower.zside(), thisTower.subDet(), thisTower.ietaAbs(), phi + 1));
      tptow[1] = EcalTriggerPrimitiveDigi(EcalTrigTowerDetId(
          thisTower.zside(), thisTower.subDet(), thisTower.ietaAbs(), phi + 2));

      if (tcpFormat_) {
        tptowTcp[0] = EcalTriggerPrimitiveDigi(
            EcalTrigTowerDetId(thisTower.zside(), thisTower.subDet(),
                               thisTower.ietaAbs(), phi + 1));
        tptowTcp[1] = EcalTriggerPrimitiveDigi(
            EcalTrigTowerDetId(thisTower.zside(), thisTower.subDet(),
                               thisTower.ietaAbs(), phi + 2));
      }
    } else {
      nrTowers = 1;
      tptow[0] = EcalTriggerPrimitiveDigi(thisTower);
      if (tcpFormat_)
        tptowTcp[0] = EcalTriggerPrimitiveDigi(thisTower);
    }

    // now fill in
    for (int nrt = 0; nrt < nrTowers; nrt++) {
      (tptow[nrt]).setSize(nrSamples_);
      if (towtp_.size() < nrSamples_) { // FIXME: only once
        edm::LogWarning("")
            << "Too few samples produced, nr is " << towtp_.size();
        break;
      }
      int isam = 0;
      for (int i = firstSample; i <= lastSample; ++i) {
        tptow[nrt].setSample(isam++, EcalTriggerPrimitiveSample(towtp_[i]));
      }
      nrTP++;
      LogDebug("EcalTPG") << " For tower " << itow << " created TP nr " << nrTP
                          << " with Et " << tptow[nrt].compressedEt();
      result.push_back(tptow[nrt]);
    }

    if (tcpFormat_) {

      for (int nrt = 0; nrt < nrTowers; nrt++) {
        tptowTcp[nrt].setSize(nrSamples_);
        if (towtp2_.size() < nrSamples_) { // FIXME: only once
          edm::LogWarning("")
              << "Too few samples produced, nr is " << towtp2_.size();
          break;
        }
        int isam = 0;
        for (int i = firstSample; i <= lastSample; ++i) {
          if (nrTowers <= 1)
            tptowTcp[nrt].setSample(isam++,
                                    EcalTriggerPrimitiveSample(towtp2_[i]));
          else {
            int et = towtp2_[i].compressedEt() / 2;
            tptowTcp[nrt].setSample(
                isam++, EcalTriggerPrimitiveSample(et, towtp2_[i].fineGrain(),
                                                   towtp2_[i].ttFlag()));
          }
        }
        resultTcp.push_back(tptowTcp[nrt]);
      }
    }
  }
  return;
}

template <class Coll>
void EcalTrigPrimFunctionalAlgo::fillMap(
    Coll const *col,
    std::vector<std::vector<std::pair<int, std::vector<typename Coll::Digi>>>>
        &towerMap) {
  typedef typename Coll::Digi Digi;

  // implementation for Barrel and Endcap

  if (col) {
    nrTowers_ = 0;
    LogDebug("EcalTPG") << "Fill mapping, Collection size = " << col->size();
    for (unsigned int i = 0; i < col->size(); ++i) {
      Digi samples((*col)[i]);
      EcalTrigTowerDetId coarser = (*eTTmap_).towerOf(samples.id());
      int index = getIndex(col, coarser);
      int stripnr = findStripNr(samples.id());

      int filled = 0;
      for (unsigned int ij = 0; ij < towerMap[index].size(); ++ij)
        filled += towerMap[index][ij].first;
      if (!filled) {
        hitTowers_[nrTowers_++] =
            std::pair<int, EcalTrigTowerDetId>(index, coarser);
      }

      // FIXME: temporary protection
      int ncryst = towerMap[index][stripnr - 1].first;
      if (ncryst >= nbMaxXtals_) {
        edm::LogError("EcalTrigPrimFunctionAlgo")
            << "! Too many xtals for TT " << coarser << " stripnr " << stripnr
            << " xtalid " << samples.id();
        continue;
      }
      ((towerMap[index])[stripnr - 1].second)[ncryst] = samples;
      (towerMap[index])[stripnr - 1].first++;
    }

    LogDebug("EcalTPG") << "fillMap"
                        << "[EcalTrigPrimFunctionalAlgo] (found " << col->size()
                        << " frames in " << towerMap.size() << " towers) ";
  } else {
    LogDebug("EcalTPG") << "FillMap - FillMap Collection size=0 !!!!";
  }
}

template <class T>
void EcalTrigPrimFunctionalAlgo::clean(
    std::vector<std::vector<std::pair<int, std::vector<T>>>> &towMap) {
  // clean internal data structures
  for (unsigned int i = 0; i < maxNrTowers_; ++i)
    for (int j = 0; j < nbMaxStrips_; ++j)
      (towMap[i])[j].first = 0;
  return;
}

template <class T>
void EcalTrigPrimFunctionalAlgo::initStructures(
    std::vector<std::vector<std::pair<int, std::vector<T>>>> &towMap) {
  // initialise internal data structures

  std::vector<T> vec0(nbMaxXtals_);
  std::vector<std::pair<int, std::vector<T>>> vec1(nbMaxStrips_);
  for (int i = 0; i < nbMaxStrips_; ++i)
    vec1[i] = std::pair<int, std::vector<T>>(0, vec0);
  towMap.resize(maxNrTowers_);
  for (unsigned int i = 0; i < maxNrTowers_; ++i)
    towMap[i] = vec1;

  std::vector<int> vecint(maxNrSamples_);
  striptp_.resize(nbMaxStrips_);
  for (int i = 0; i < nbMaxStrips_; ++i)
    striptp_[i] = vecint;
}

#endif