EcalEBTrigPrimAlgos/src/EcalEBTrigPrimTestAlgo.cc

0001 /** \class EcalEBTrigPrimTestAlgo
0002  *
0003  * EcalEBTrigPrimTestAlgo
0004  * starting point for Phase II: build TPs out of Phase I digis to start building the
0005  * infrastructures
0006  *
0007  *
0008  ************************************************************/
0009 #include <string>
0010 #include <algorithm>
0011 #include <numeric>
0012 #include <functional>
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 //#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0017 //#include "Geometry/Records/interface/CaloGeometryRecord.h"
0018 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
0019 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
0020
0021 #include "SimCalorimetry/EcalEBTrigPrimAlgos/interface/EcalEBTrigPrimTestAlgo.h"
0022
0023 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0024 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0025 #include "DataFormats/EcalDigi/interface/EEDataFrame.h"
0026 #include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
0027 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0028 #include "DataFormats/EcalDetId/interface/EcalTriggerElectronicsId.h"
0029
0030 #include "CondFormats/EcalObjects/interface/EcalTPGPedestals.h"
0031 #include "CondFormats/DataRecord/interface/EcalTPGPedestalsRcd.h"
0032
0033 #include <TTree.h>
0034 #include <TMath.h>
0035
0036 //----------------------------------------------------------------------
0037
0038 const unsigned int EcalEBTrigPrimTestAlgo::nrSamples_ = 5;
0039 const unsigned int EcalEBTrigPrimTestAlgo::maxNrTowers_ = 2448;
0040 const unsigned int EcalEBTrigPrimTestAlgo::maxNrSamplesOut_ = 10;
0041
0042 // not BarrelOnly
0043 EcalEBTrigPrimTestAlgo::EcalEBTrigPrimTestAlgo(const EcalTrigTowerConstituentsMap *eTTmap,
0044                                                const CaloGeometry *theGeometry,
0045                                                int nSam,
0046                                                int binofmax,
0047                                                bool tcpFormat,
0048                                                bool debug,
0049                                                bool famos)
0050     : eTTmap_(eTTmap),
0051       theGeometry_(theGeometry),
0052       nSamples_(nSam),
0053       binOfMaximum_(binofmax),
0054       tcpFormat_(tcpFormat),
0055       barrelOnly_(false),
0056       debug_(debug),
0057       famos_(famos) {
0058   maxNrSamples_ = 10;
0059   this->init();
0060 }
0061
0062 //barrel only
0063 EcalEBTrigPrimTestAlgo::EcalEBTrigPrimTestAlgo(int nSam, int binofmax, bool tcpFormat, bool debug, bool famos)
0064     : nSamples_(nSam),
0065       binOfMaximum_(binofmax),
0066       tcpFormat_(tcpFormat),
0067       barrelOnly_(true),
0068       debug_(debug),
0069       famos_(famos)
0070
0071 {
0072   maxNrSamples_ = 10;
0073   this->init();
0074 }
0075
0076 //----------------------------------------------------------------------
0077 void EcalEBTrigPrimTestAlgo::init() {
0078   // initialise data structures
0079   theMapping_ = new EcalElectronicsMapping();
0080   initStructures(towerMapEB_);
0081   hitTowers_.resize(maxNrTowers_);
0082
0083   linearizer_.resize(nbMaxXtals_);
0084   for (int i = 0; i < nbMaxXtals_; i++)
0085     linearizer_[i] = new EcalEBFenixLinearizer(famos_);
0086
0087   //
0088   std::vector<int> v;
0089   v.resize(maxNrSamples_);
0090   lin_out_.resize(nbMaxXtals_);
0091   for (int i = 0; i < 5; i++)
0092     lin_out_[i] = v;
0093   //
0094   amplitude_filter_ = new EcalEBFenixAmplitudeFilter();
0095   filt_out_.resize(maxNrSamples_);
0096   peak_out_.resize(maxNrSamples_);
0097   // these two are dummy
0098   fgvb_out_.resize(maxNrSamples_);
0099   fgvb_out_temp_.resize(maxNrSamples_);
0100   //
0101   peak_finder_ = new EcalEBFenixPeakFinder();
0102   fenixFormatterEB_ = new EcalEBFenixStripFormatEB();
0103   format_out_.resize(maxNrSamples_);
0104   //
0105   fenixTcpFormat_ = new EcalEBFenixTcpFormat(tcpFormat_, debug_, famos_, binOfMaximum_);
0106   tcpformat_out_.resize(maxNrSamples_);
0107 }
0108 //----------------------------------------------------------------------
0109
0110 EcalEBTrigPrimTestAlgo::~EcalEBTrigPrimTestAlgo() {
0111   for (int i = 0; i < nbMaxXtals_; i++)
0112     delete linearizer_[i];
0113   delete amplitude_filter_;
0114   delete peak_finder_;
0115   delete fenixFormatterEB_;
0116   delete fenixTcpFormat_;
0117   delete theMapping_;
0118 }
0119
0120 void EcalEBTrigPrimTestAlgo::run(EBDigiCollection const *digi,
0121                                  EcalEBTrigPrimDigiCollection &result,
0122                                  EcalEBTrigPrimDigiCollection &resultTcp) {
0123   //typedef typename Coll::Digi Digi;
0124   if (debug_) {
0125     std::cout << "  EcalEBTrigPrimTestAlgo: Testing that the algorythm with digis is well plugged " << std::endl;
0126     std::cout << "  EcalEBTrigPrimTestAlgo: digi size " << digi->size() << std::endl;
0127   }
0128
0129   uint16_t etInADC;
0130   EcalEBTriggerPrimitiveDigi tp;
0131   int firstSample = binOfMaximum_ - 1 - nrSamples_ / 2;
0132   int lastSample = binOfMaximum_ - 1 + nrSamples_ / 2;
0133
0134   if (debug_) {
0135     std::cout << "  binOfMaximum_ " << binOfMaximum_ << " nrSamples_" << nrSamples_ << std::endl;
0136     std::cout << " first sample " << firstSample << " last " << lastSample << std::endl;
0137   }
0138
0139   clean(towerMapEB_);
0140   fillMap(digi, towerMapEB_);
0141
0142   for (int itow = 0; itow < nrTowers_; ++itow) {
0143     int index = hitTowers_[itow].first;
0144     const EcalTrigTowerDetId &thisTower = hitTowers_[itow].second;
0145     if (debug_)
0146       std::cout << " Data for TOWER num " << itow << " index " << index << " TowerId " << thisTower << " size "
0147                 << towerMapEB_[itow].size() << std::endl;
0148     // loop over all strips assigned to this trigger tower
0149     int nxstals = 0;
0150     for (unsigned int iStrip = 0; iStrip < towerMapEB_[itow].size(); ++iStrip) {
0151       if (debug_)
0152         std::cout << " Data for STRIP num " << iStrip << std::endl;
0153       std::vector<EBDataFrame> &dataFrames =
0154           (towerMapEB_[index])[iStrip].second;  //vector of dataframes for this strip, size; nr of crystals/strip
0155
0156       nxstals = (towerMapEB_[index])[iStrip].first;
0157       if (nxstals <= 0)
0158         continue;
0159       if (debug_)
0160         std::cout << " Number of xTals " << nxstals << std::endl;
0161
0162       const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(dataFrames[0].id());
0163       uint32_t stripid = elId.rawId() & 0xfffffff8;
0164
0165       // loop over the xstals in a strip
0166       for (int iXstal = 0; iXstal < nxstals; iXstal++) {
0167         const EBDetId &myid = dataFrames[iXstal].id();
0168         tp = EcalEBTriggerPrimitiveDigi(myid);
0169         tp.setSize(nrSamples_);
0170
0171         if (debug_) {
0172           std::cout << std::endl;
0173           std::cout << "iXstal= " << iXstal << " id " << dataFrames[iXstal].id() << " EBDataFrame is: " << std::endl;
0174           for (int i = 0; i < dataFrames[iXstal].size(); i++) {
0175             std::cout << " " << std::dec << dataFrames[iXstal][i].adc();
0176           }
0177           std::cout << std::endl;
0178         }
0179         //   Call the linearizer
0180         this->getLinearizer(iXstal)->setParameters(
0181             dataFrames[iXstal].id().rawId(), ecaltpPed_, ecaltpLin_, ecaltpgBadX_);
0182         this->getLinearizer(iXstal)->process(dataFrames[iXstal], lin_out_[iXstal]);
0183
0184         for (unsigned int i = 0; i < lin_out_[iXstal].size(); i++) {
0185           if ((lin_out_[iXstal])[i] > 0X3FFFF)
0186             (lin_out_[iXstal])[i] = 0X3FFFF;
0187         }
0188
0189         if (debug_) {
0190           std::cout << "output of linearizer for channel " << iXstal << std::endl;
0191           std::cout << " lin_out[iXstal].size()= " << std::dec << lin_out_[iXstal].size() << std::endl;
0192           for (unsigned int i = 0; i < lin_out_[iXstal].size(); i++) {
0193             std::cout << " " << std::dec << (lin_out_[iXstal])[i];
0194           }
0195           std::cout << std::endl;
0196         }
0197
0198         // Call the amplitude filter
0199         this->getFilter()->setParameters(stripid, ecaltpgWeightMap_, ecaltpgWeightGroup_);
0200         this->getFilter()->process(lin_out_[iXstal], filt_out_, fgvb_out_temp_, fgvb_out_);
0201
0202         if (debug_) {
0203           std::cout << "output of filter is a vector of size: " << std::dec << filt_out_.size() << std::endl;
0204           for (unsigned int ix = 0; ix < filt_out_.size(); ix++) {
0205             std::cout << std::dec << filt_out_[ix] << " ";
0206           }
0207           std::cout << std::endl;
0208         }
0209
0210         // call peakfinder
0211         this->getPeakFinder()->process(filt_out_, peak_out_);
0212
0213         if (debug_) {
0214           std::cout << "output of peakfinder is a vector of size: " << std::dec << peak_out_.size() << std::endl;
0215           for (unsigned int ix = 0; ix < peak_out_.size(); ix++) {
0216             std::cout << std::dec << peak_out_[ix] << " ";
0217           }
0218           std::cout << std::endl;
0219         }
0220
0221         // call formatter
0222         this->getFormatterEB()->setParameters(stripid, ecaltpgSlidW_);
0223         this->getFormatterEB()->process(fgvb_out_, peak_out_, filt_out_, format_out_);
0224
0225         if (debug_) {
0226           std::cout << "output of formatter is a vector of size: " << format_out_.size() << std::endl;
0227           for (unsigned int i = 0; i < format_out_.size(); i++) {
0228             std::cout << " " << std::dec << format_out_[i] << " ";
0229           }
0230           std::cout << std::endl;
0231         }
0232
0233         // call final tcp formatter
0234         this->getFormatter()->setParameters(
0235             thisTower.rawId(), ecaltpgLutGroup_, ecaltpgLut_, ecaltpgBadTT_, ecaltpgSpike_);
0236         this->getFormatter()->process(format_out_, tcpformat_out_);
0237         // loop over the time samples and fill the TP
0238         int nSam = 0;
0239         for (int iSample = firstSample; iSample <= lastSample; ++iSample) {
0240           etInADC = tcpformat_out_[iSample];
0241           if (debug_)
0242             std::cout << " format_out " << tcpformat_out_[iSample] << " etInADC " << etInADC << std::endl;
0243
0244           bool isASpike = false;  // no spikes for now
0245           int timing = 0;         //  set to 0  value for now
0246           tp.setSample(nSam, EcalEBTriggerPrimitiveSample(etInADC, isASpike, timing));
0247
0248           nSam++;
0249           if (debug_)
0250             std::cout << "in TestAlgo"
0251                       << " tp size " << tp.size() << std::endl;
0252         }
0253
0254         if (!tcpFormat_)
0255           result.push_back(tp);
0256         else
0257           resultTcp.push_back(tp);
0258
0259       }  // Loop over the xStals
0260
0261     }  //loop over strips in one tower
0262   }
0263
0264   /*
0265   for (unsigned int i=0;i<digi->size();i++) {
0266     EBDataFrame myFrame((*digi)[i]);
0267     const EBDetId & myid1 = myFrame.id();
0268     tp=  EcalTriggerPrimitiveDigi(  myid1);
0269     tp.setSize( myFrame.size());
0270     int nSam=0;
0271
0272     if (debug_) {
0273       std::cout << " data frame size " << myFrame.size() << " Id " <<  myFrame.id()  << std::endl;
0274       std::cout << " Sample data ADC: " << std::endl;
0275       for (int iSample=0; iSample<myFrame.size(); iSample++) {
0276     std::cout << " " << std::dec<< myFrame.sample(iSample).adc() ;
0277       }
0278       std::cout<<std::endl;
0279     }
0280
0281
0282     this->getLinearizer(i)->setParameters( myFrame.id().rawId(),ecaltpPed_,ecaltpLin_,ecaltpgBadX_) ;
0283     //this->getLinearizer(i)->process( myFrame,lin_out_[i]);
0284
0285     if (debug_) {
0286       std::cout<< "cryst: "<< i <<"  value : "<<std::dec<<std::endl;
0287       std::cout<<" lin_out[i].size()= "<<std::dec<<lin_out_[i].size()<<std::endl;
0288       for (unsigned int j =0; j<lin_out_[i].size();j++){
0289     std::cout <<" "<<std::dec<<(lin_out_[i])[j];
0290       }
0291       std::cout<<std::endl;
0292     }
0293
0294
0295     for (int iSample=0; iSample<myFrame.size(); iSample++) {
0296       etInADC= myFrame.sample(iSample).adc();
0297       EcalEBTriggerPrimitiveSample mysam(etInADC);
0298       tp.setSample(nSam, mysam );
0299       nSam++;
0300       if (debug_) std::cout << "in TestAlgo" <<" tp size "<<tp.size() << std::endl;
0301     }
0302
0303     if (!tcpFormat_)
0304       result.push_back(tp);
0305     else
0306       resultTcp.push_back(tp);
0307
0308
0309     if (debug_) std::cout << " result size " << result.size() << std::endl;
0310
0311
0312
0313   }
0314   */
0315 }
0316
0317 //----------------------------------------------------------------------
0318
0319 int EcalEBTrigPrimTestAlgo::findStripNr(const EBDetId &id) {
0320   int stripnr;
0321   int n = ((id.ic() - 1) % 100) / 20;  //20 corresponds to 4 * ecal_barrel_crystals_per_strip FIXME!!
0322   if (id.ieta() < 0)
0323     stripnr = n + 1;
0324   else
0325     stripnr = nbMaxStrips_ - n;
0326   return stripnr;
0327 }