EcalEBTrigPrimTestAlgo.cc

Go to the documentation of this file.

#include <string>
#include <algorithm>
#include <numeric>
#include <functional>

#include "FWCore/MessageLogger/interface/MessageLogger.h"

//#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
//#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
#include "Geometry/EcalMapping/interface/EcalMappingRcd.h"

#include "SimCalorimetry/EcalEBTrigPrimAlgos/interface/EcalEBTrigPrimTestAlgo.h"

#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EcalDigi/interface/EBDataFrame.h"
#include "DataFormats/EcalDigi/interface/EEDataFrame.h"
#include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EcalTriggerElectronicsId.h"

#include "CondFormats/EcalObjects/interface/EcalTPGPedestals.h"
#include "CondFormats/DataRecord/interface/EcalTPGPedestalsRcd.h"

#include <TTree.h>
#include <TMath.h>

//----------------------------------------------------------------------

const unsigned int EcalEBTrigPrimTestAlgo::nrSamples_ = 5;
const unsigned int EcalEBTrigPrimTestAlgo::maxNrTowers_ = 2448;
const unsigned int EcalEBTrigPrimTestAlgo::maxNrSamplesOut_ = 10;

// not BarrelOnly
EcalEBTrigPrimTestAlgo::EcalEBTrigPrimTestAlgo(const EcalTrigTowerConstituentsMap *eTTmap,
                                               const CaloGeometry *theGeometry,
                                               int nSam,
                                               int binofmax,
                                               bool tcpFormat,
                                               bool debug,
                                               bool famos)
    : eTTmap_(eTTmap),
      theGeometry_(theGeometry),
      nSamples_(nSam),
      binOfMaximum_(binofmax),
      tcpFormat_(tcpFormat),
      barrelOnly_(false),
      debug_(debug),
      famos_(famos) {
  maxNrSamples_ = 10;
  this->init();
}

//barrel only
EcalEBTrigPrimTestAlgo::EcalEBTrigPrimTestAlgo(int nSam, int binofmax, bool tcpFormat, bool debug, bool famos)
    : nSamples_(nSam),
      binOfMaximum_(binofmax),
      tcpFormat_(tcpFormat),
      barrelOnly_(true),
      debug_(debug),
      famos_(famos)

{
  maxNrSamples_ = 10;
  this->init();
}

//----------------------------------------------------------------------
void EcalEBTrigPrimTestAlgo::init() {
  // initialise data structures
  theMapping_ = new EcalElectronicsMapping();
  initStructures(towerMapEB_);
  hitTowers_.resize(maxNrTowers_);

  linearizer_.resize(nbMaxXtals_);
  for (int i = 0; i < nbMaxXtals_; i++)
    linearizer_[i] = new EcalEBFenixLinearizer(famos_);

  //
  std::vector<int> v;
  v.resize(maxNrSamples_);
  lin_out_.resize(nbMaxXtals_);
  for (int i = 0; i < 5; i++)
    lin_out_[i] = v;
  //
  amplitude_filter_ = new EcalEBFenixAmplitudeFilter();
  filt_out_.resize(maxNrSamples_);
  peak_out_.resize(maxNrSamples_);
  // these two are dummy
  fgvb_out_.resize(maxNrSamples_);
  fgvb_out_temp_.resize(maxNrSamples_);
  //
  peak_finder_ = new EcalEBFenixPeakFinder();
  fenixFormatterEB_ = new EcalEBFenixStripFormatEB();
  format_out_.resize(maxNrSamples_);
  //
  fenixTcpFormat_ = new EcalEBFenixTcpFormat(tcpFormat_, debug_, famos_, binOfMaximum_);
  tcpformat_out_.resize(maxNrSamples_);
}
//----------------------------------------------------------------------

EcalEBTrigPrimTestAlgo::~EcalEBTrigPrimTestAlgo() {
  for (int i = 0; i < nbMaxXtals_; i++)
    delete linearizer_[i];
  delete amplitude_filter_;
  delete peak_finder_;
  delete fenixFormatterEB_;
  delete fenixTcpFormat_;
  delete theMapping_;
}

void EcalEBTrigPrimTestAlgo::run(EBDigiCollection const *digi,
                                 EcalEBTrigPrimDigiCollection &result,
                                 EcalEBTrigPrimDigiCollection &resultTcp) {
  //typedef typename Coll::Digi Digi;
  if (debug_) {
    std::cout << "  EcalEBTrigPrimTestAlgo: Testing that the algorythm with digis is well plugged " << std::endl;
    std::cout << "  EcalEBTrigPrimTestAlgo: digi size " << digi->size() << std::endl;
  }

  uint16_t etInADC;
  EcalEBTriggerPrimitiveDigi tp;
  int firstSample = binOfMaximum_ - 1 - nrSamples_ / 2;
  int lastSample = binOfMaximum_ - 1 + nrSamples_ / 2;

  if (debug_) {
    std::cout << "  binOfMaximum_ " << binOfMaximum_ << " nrSamples_" << nrSamples_ << std::endl;
    std::cout << " first sample " << firstSample << " last " << lastSample << std::endl;
  }

  clean(towerMapEB_);
  fillMap(digi, towerMapEB_);

  for (int itow = 0; itow < nrTowers_; ++itow) {
    int index = hitTowers_[itow].first;
    const EcalTrigTowerDetId &thisTower = hitTowers_[itow].second;
    if (debug_)
      std::cout << " Data for TOWER num " << itow << " index " << index << " TowerId " << thisTower << " size "
                << towerMapEB_[itow].size() << std::endl;
    // loop over all strips assigned to this trigger tower
    int nxstals = 0;
    for (unsigned int iStrip = 0; iStrip < towerMapEB_[itow].size(); ++iStrip) {
      if (debug_)
        std::cout << " Data for STRIP num " << iStrip << std::endl;
      std::vector<EBDataFrame> &dataFrames =
          (towerMapEB_[index])[iStrip].second;  //vector of dataframes for this strip, size; nr of crystals/strip

      nxstals = (towerMapEB_[index])[iStrip].first;
      if (nxstals <= 0)
        continue;
      if (debug_)
        std::cout << " Number of xTals " << nxstals << std::endl;

      const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(dataFrames[0].id());
      uint32_t stripid = elId.rawId() & 0xfffffff8;

      // loop over the xstals in a strip
      for (int iXstal = 0; iXstal < nxstals; iXstal++) {
        const EBDetId &myid = dataFrames[iXstal].id();
        tp = EcalEBTriggerPrimitiveDigi(myid);
        tp.setSize(nrSamples_);

        if (debug_) {
          std::cout << std::endl;
          std::cout << "iXstal= " << iXstal << " id " << dataFrames[iXstal].id() << " EBDataFrame is: " << std::endl;
          for (int i = 0; i < dataFrames[iXstal].size(); i++) {
            std::cout << " " << std::dec << dataFrames[iXstal][i].adc();
          }
          std::cout << std::endl;
        }
        //   Call the linearizer
        this->getLinearizer(iXstal)->setParameters(
            dataFrames[iXstal].id().rawId(), ecaltpPed_, ecaltpLin_, ecaltpgBadX_);
        this->getLinearizer(iXstal)->process(dataFrames[iXstal], lin_out_[iXstal]);

        for (unsigned int i = 0; i < lin_out_[iXstal].size(); i++) {
          if ((lin_out_[iXstal])[i] > 0X3FFFF)
            (lin_out_[iXstal])[i] = 0X3FFFF;
        }

        if (debug_) {
          std::cout << "output of linearizer for channel " << iXstal << std::endl;
          std::cout << " lin_out[iXstal].size()= " << std::dec << lin_out_[iXstal].size() << std::endl;
          for (unsigned int i = 0; i < lin_out_[iXstal].size(); i++) {
            std::cout << " " << std::dec << (lin_out_[iXstal])[i];
          }
          std::cout << std::endl;
        }

        // Call the amplitude filter
        this->getFilter()->setParameters(stripid, ecaltpgWeightMap_, ecaltpgWeightGroup_);
        this->getFilter()->process(lin_out_[iXstal], filt_out_, fgvb_out_temp_, fgvb_out_);

        if (debug_) {
          std::cout << "output of filter is a vector of size: " << std::dec << filt_out_.size() << std::endl;
          for (unsigned int ix = 0; ix < filt_out_.size(); ix++) {
            std::cout << std::dec << filt_out_[ix] << " ";
          }
          std::cout << std::endl;
        }

        // call peakfinder
        this->getPeakFinder()->process(filt_out_, peak_out_);

        if (debug_) {
          std::cout << "output of peakfinder is a vector of size: " << std::dec << peak_out_.size() << std::endl;
          for (unsigned int ix = 0; ix < peak_out_.size(); ix++) {
            std::cout << std::dec << peak_out_[ix] << " ";
          }
          std::cout << std::endl;
        }

        // call formatter
        this->getFormatterEB()->setParameters(stripid, ecaltpgSlidW_);
        this->getFormatterEB()->process(fgvb_out_, peak_out_, filt_out_, format_out_);

        if (debug_) {
          std::cout << "output of formatter is a vector of size: " << format_out_.size() << std::endl;
          for (unsigned int i = 0; i < format_out_.size(); i++) {
            std::cout << " " << std::dec << format_out_[i] << " ";
          }
          std::cout << std::endl;
        }

        // call final tcp formatter
        this->getFormatter()->setParameters(
            thisTower.rawId(), ecaltpgLutGroup_, ecaltpgLut_, ecaltpgBadTT_, ecaltpgSpike_);
        this->getFormatter()->process(format_out_, tcpformat_out_);
        // loop over the time samples and fill the TP
        int nSam = 0;
        for (int iSample = firstSample; iSample <= lastSample; ++iSample) {
          etInADC = tcpformat_out_[iSample];
          if (debug_)
            std::cout << " format_out " << tcpformat_out_[iSample] << " etInADC " << etInADC << std::endl;

          bool isASpike = false;  // no spikes for now
          int timing = 0;         //  set to 0  value for now
          tp.setSample(nSam, EcalEBTriggerPrimitiveSample(etInADC, isASpike, timing));

          nSam++;
          if (debug_)
            std::cout << "in TestAlgo"
                      << " tp size " << tp.size() << std::endl;
        }

        if (!tcpFormat_)
          result.push_back(tp);
        else
          resultTcp.push_back(tp);

      }  // Loop over the xStals

    }  //loop over strips in one tower
  }

  /*
  for (unsigned int i=0;i<digi->size();i++) {
    EBDataFrame myFrame((*digi)[i]);  
    const EBDetId & myid1 = myFrame.id();
    tp=  EcalTriggerPrimitiveDigi(  myid1);   
    tp.setSize( myFrame.size());
    int nSam=0;

    if (debug_) {
      std::cout << " data frame size " << myFrame.size() << " Id " <<  myFrame.id()  << std::endl;
      std::cout << " Sample data ADC: " << std::endl;
      for (int iSample=0; iSample<myFrame.size(); iSample++) {
    std::cout << " " << std::dec<< myFrame.sample(iSample).adc() ;
      }
      std::cout<<std::endl;
    }

    
    this->getLinearizer(i)->setParameters( myFrame.id().rawId(),ecaltpPed_,ecaltpLin_,ecaltpgBadX_) ; 
    //this->getLinearizer(i)->process( myFrame,lin_out_[i]);

    if (debug_) {
      std::cout<< "cryst: "<< i <<"  value : "<<std::dec<<std::endl;
      std::cout<<" lin_out[i].size()= "<<std::dec<<lin_out_[i].size()<<std::endl;
      for (unsigned int j =0; j<lin_out_[i].size();j++){
    std::cout <<" "<<std::dec<<(lin_out_[i])[j];
      }
      std::cout<<std::endl;
    }


    for (int iSample=0; iSample<myFrame.size(); iSample++) {
      etInADC= myFrame.sample(iSample).adc();
      EcalEBTriggerPrimitiveSample mysam(etInADC);
      tp.setSample(nSam, mysam );
      nSam++;
      if (debug_) std::cout << "in TestAlgo" <<" tp size "<<tp.size() << std::endl;
    }

    if (!tcpFormat_)
      result.push_back(tp);
    else 
      resultTcp.push_back(tp);
    
   
    if (debug_) std::cout << " result size " << result.size() << std::endl;
    
    
    
  }
  */
}

//----------------------------------------------------------------------

int EcalEBTrigPrimTestAlgo::findStripNr(const EBDetId &id) {
  int stripnr;
  int n = ((id.ic() - 1) % 100) / 20;  //20 corresponds to 4 * ecal_barrel_crystals_per_strip FIXME!!
  if (id.ieta() < 0)
    stripnr = n + 1;
  else
    stripnr = nbMaxStrips_ - n;
  return stripnr;
}