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;
 
EcalEBTrigPrimTestAlgo::EcalEBTrigPrimTestAlgo(
    const edm::EventSetup &setup, int nSam, int binofmax, bool tcpFormat, bool barrelOnly, bool debug, bool famos)
    : nSamples_(nSam),
      binOfMaximum_(binofmax),
      tcpFormat_(tcpFormat),
      barrelOnly_(barrelOnly),
      debug_(debug),
      famos_(famos)
 
{
  maxNrSamples_ = 10;
  this->init(setup);
}
 
//----------------------------------------------------------------------
void EcalEBTrigPrimTestAlgo::init(const edm::EventSetup &setup) {
  if (!barrelOnly_) {
    //edm::ESHandle<CaloGeometry> theGeometry;
    //    edm::ESHandle<CaloSubdetectorGeometry> theEndcapGeometry_handle;
    setup.get<CaloGeometryRecord>().get(theGeometry);
    setup.get<IdealGeometryRecord>().get(eTTmap_);
  }
 
  // initialise data structures
  initStructures(towerMapEB_);
  hitTowers_.resize(maxNrTowers_);
 
  linearizer_.resize(nbMaxXtals_);
  for (int i = 0; i < nbMaxXtals_; i++)
    linearizer_[i] = new EcalFenixLinearizer(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 EcalFenixAmplitudeFilter();
  filt_out_.resize(maxNrSamples_);
  peak_out_.resize(maxNrSamples_);
  // these two are dummy
  fgvb_out_.resize(maxNrSamples_);
  fgvb_out_temp_.resize(maxNrSamples_);
  //
  peak_finder_ = new EcalFenixPeakFinder();
  fenixFormatterEB_ = new EcalFenixStripFormatEB();
  format_out_.resize(maxNrSamples_);
  //
  fenixTcpFormat_ = new EcalFenixTcpFormat(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_;
}
 
void EcalEBTrigPrimTestAlgo::run(const edm::EventSetup &setup,
                                 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;
}