Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:22

0001 #include <string>
0002 #include <algorithm>
0003 #include <numeric>
0004 #include <functional>
0005 
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
0008 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
0009 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0010 
0011 #include "SimCalorimetry/EcalEBTrigPrimAlgos/interface/EcalEBPhase2TrigPrimAlgo.h"
0012 
0013 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0014 #include "DataFormats/EcalDigi/interface/EBDataFrame_Ph2.h"
0015 #include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0017 #include "DataFormats/EcalDetId/interface/EcalTriggerElectronicsId.h"
0018 
0019 #include "CondFormats/EcalObjects/interface/EcalTPGPedestals.h"
0020 #include "CondFormats/DataRecord/interface/EcalTPGPedestalsRcd.h"
0021 
0022 #include "DataFormats/EcalDigi/interface/EcalConstants.h"
0023 
0024 #include <TTree.h>
0025 #include <TMath.h>
0026 
0027 //----------------------------------------------------------------------
0028 
0029 const unsigned int EcalEBPhase2TrigPrimAlgo::nrSamples_ =
0030     ecalPh2::sampleSize;  // set to 16 samples, might change (less than 16) in the future
0031 const unsigned int EcalEBPhase2TrigPrimAlgo::maxNrTowers_ = 2448;  // number of towers in EB
0032 
0033 EcalEBPhase2TrigPrimAlgo::EcalEBPhase2TrigPrimAlgo(const EcalTrigTowerConstituentsMap *eTTmap,
0034                                                    const CaloGeometry *theGeometry,
0035                                                    int binofmax,
0036                                                    bool debug)
0037     : eTTmap_(eTTmap),
0038       theGeometry_(theGeometry),
0039       binOfMaximum_(binofmax),
0040       debug_(debug)
0041 
0042 {
0043   maxNrSamples_ = ecalPh2::sampleSize;
0044   this->init();
0045 }
0046 
0047 void EcalEBPhase2TrigPrimAlgo::init() {
0048   theMapping_ = new EcalElectronicsMapping();
0049   // initialise data structures
0050   initStructures(towerMapEB_);
0051   hitTowers_.resize(maxNrTowers_);
0052 
0053   linearizer_ = new EcalEBPhase2Linearizer(debug_);
0054   lin_out_.resize(maxNrSamples_);
0055 
0056   amplitude_reconstructor_ = new EcalEBPhase2AmplitudeReconstructor(debug_);
0057   filt_out_.resize(maxNrSamples_);
0058 
0059   tpFormatter_ = new EcalEBPhase2TPFormatter(debug_);
0060   outEt_.resize(maxNrSamples_);
0061   outTime_.resize(maxNrSamples_);
0062 
0063   //
0064 
0065   time_reconstructor_ = new EcalEBPhase2TimeReconstructor(debug_);
0066   time_out_.resize(maxNrSamples_);
0067   spike_tagger_ = new EcalEBPhase2SpikeTagger(debug_);
0068 }
0069 //----------------------------------------------------------------------
0070 
0071 EcalEBPhase2TrigPrimAlgo::~EcalEBPhase2TrigPrimAlgo() {
0072   delete linearizer_;
0073   delete amplitude_reconstructor_;
0074   delete time_reconstructor_;
0075   delete spike_tagger_;
0076   delete tpFormatter_;
0077   delete theMapping_;
0078 }
0079 
0080 void EcalEBPhase2TrigPrimAlgo::run(EBDigiCollectionPh2 const *digi, EcalEBPhase2TrigPrimDigiCollection &result) {
0081   if (debug_)
0082     LogDebug("") << "  EcalEBPhase2TrigPrimAlgo: digi size " << digi->size() << std::endl;
0083 
0084   EcalEBPhase2TriggerPrimitiveDigi tp;
0085   int firstSample = binOfMaximum_ - 1 - nrSamples_ / 2;
0086   int lastSample = binOfMaximum_ - 1 + nrSamples_ / 2;
0087 
0088   if (debug_) {
0089     LogDebug("") << "  binOfMaximum_ " << binOfMaximum_ << " nrSamples_" << nrSamples_ << std::endl;
0090     LogDebug("") << " first sample " << firstSample << " last " << lastSample << std::endl;
0091   }
0092 
0093   clean(towerMapEB_);
0094   fillMap(digi, towerMapEB_);
0095 
0096   int iChannel = 0;
0097   int nXinBCP = 0;
0098   for (int itow = 0; itow < nrTowers_; ++itow) {
0099     int index = hitTowers_[itow].first;
0100     const EcalTrigTowerDetId &thisTower = hitTowers_[itow].second;
0101     if (debug_)
0102       LogDebug("") << " Data for TOWER num " << itow << " index " << index << " TowerId " << thisTower << " zside "
0103                    << thisTower.zside() << " ieta " << thisTower.ieta() << " iphi " << thisTower.iphi() << " size "
0104                    << towerMapEB_[itow].size() << std::endl;
0105 
0106     // loop over all strips assigned to this trigger tower
0107     int nxstals = 0;
0108     for (unsigned int iStrip = 0; iStrip < towerMapEB_[itow].size(); ++iStrip) {
0109       if (debug_)
0110         LogDebug("") << " Data for STRIP num " << iStrip << std::endl;
0111       std::vector<EBDataFrame_Ph2> &dataFrames =
0112           (towerMapEB_[index])[iStrip].second;  //vector of dataframes for this strip, size; nr of crystals/strip
0113 
0114       nxstals = (towerMapEB_[index])[iStrip].first;
0115       if (nxstals <= 0)
0116         continue;
0117       if (debug_)
0118         LogDebug("") << " Number of xTals " << nxstals << std::endl;
0119 
0120       //const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(dataFrames[0].id());
0121 
0122       // loop over the xstals in a strip
0123 
0124       for (int iXstal = 0; iXstal < nxstals; iXstal++) {
0125         const EBDetId &myid = dataFrames[iXstal].id();
0126 
0127         nXinBCP++;
0128         if (debug_) {
0129           LogDebug("") << " Data for TOWER num " << itow << " index " << index << " TowerId " << thisTower << " size "
0130                        << towerMapEB_[itow].size() << std::endl;
0131           LogDebug("") << "nXinBCP " << nXinBCP << " myid rawId " << myid.rawId() << " xTal iEta " << myid.ieta()
0132                        << " iPhi " << myid.iphi() << std::endl;
0133         }
0134 
0135         tp = EcalEBPhase2TriggerPrimitiveDigi(myid);
0136         tp.setSize(nrSamples_);
0137 
0138         iChannel++;
0139         if (debug_) {
0140           LogDebug("") << " " << std::endl;
0141           LogDebug("") << " ******  iChannel " << iChannel << std::endl;
0142           for (int i = 0; i < dataFrames[iXstal].size(); i++) {
0143             LogDebug("") << " " << dataFrames[iXstal][i].adc();
0144           }
0145           LogDebug("") << " " << std::endl;
0146         }
0147 
0148         if (debug_) {
0149           LogDebug("") << std::endl;
0150           EBDetId id = dataFrames[iXstal].id();
0151           LogDebug("") << "iXstal= " << iXstal << std::endl;
0152           LogDebug("") << "iXstal= " << iXstal << " id " << id << " EcalDataFrame_Ph2 is: " << std::endl;
0153           for (int i = 0; i < dataFrames[iXstal].size(); i++) {
0154             LogDebug("") << " " << std::dec << dataFrames[iXstal][i].adc();
0155           }
0156           LogDebug("") << std::endl;
0157         }
0158 
0159         //   Call the linearizer
0160         this->getLinearizer()->setParameters(dataFrames[iXstal].id(), ecaltpPed_, ecaltpLin_, ecaltpgBadX_);
0161         this->getLinearizer()->process(dataFrames[iXstal], lin_out_);
0162 
0163         for (unsigned int i = 0; i < lin_out_.size(); i++) {
0164           if (lin_out_[i] > 0X3FFF)
0165             lin_out_[i] = 0X3FFF;
0166         }
0167 
0168         if (debug_) {
0169           LogDebug("") << "EcalEBPhase2TrigPrimAlgo output of linearize for channel " << iXstal << std::endl;
0170           for (unsigned int i = 0; i < lin_out_.size(); i++) {
0171             LogDebug("") << " " << std::dec << lin_out_[i];
0172           }
0173           LogDebug("") << std::endl;
0174         }
0175 
0176         // call spike finder right after the linearizer
0177         this->getSpikeTagger()->setParameters(dataFrames[iXstal].id(), ecaltpPed_, ecaltpLin_, ecaltpgBadX_);
0178         bool isASpike = this->getSpikeTagger()->process(lin_out_);
0179 
0180         //if (!isASpike) {
0181 
0182         // Call the amplitude reconstructor
0183         this->getAmplitudeFinder()->setParameters(myid.rawId(), ecaltpgAmplWeightMap_, ecaltpgWeightGroup_);
0184         this->getAmplitudeFinder()->process(lin_out_, filt_out_);
0185 
0186         if (debug_) {
0187           LogDebug("") << "EcalEBPhase2TrigPrimAlgo output of amp finder is a vector of size: " << std::dec
0188                        << time_out_.size() << std::endl;
0189           for (unsigned int ix = 0; ix < filt_out_.size(); ix++) {
0190             LogDebug("") << std::dec << filt_out_[ix] << " ";
0191           }
0192           LogDebug("") << std::endl;
0193         }
0194 
0195         if (debug_) {
0196           LogDebug("") << " Ampl "
0197                        << " ";
0198           for (unsigned int ix = 0; ix < filt_out_.size(); ix++) {
0199             LogDebug("") << std::dec << filt_out_[ix] << " ";
0200           }
0201           LogDebug("") << std::endl;
0202         }
0203 
0204         // call time finder
0205         this->getTimeFinder()->setParameters(myid.rawId(), ecaltpgTimeWeightMap_, ecaltpgWeightGroup_);
0206         this->getTimeFinder()->process(lin_out_, filt_out_, time_out_);
0207 
0208         if (debug_) {
0209           LogDebug("") << " Time "
0210                        << " ";
0211           for (unsigned int ix = 0; ix < time_out_.size(); ix++) {
0212             LogDebug("") << std::dec << time_out_[ix] << " ";
0213           }
0214           LogDebug("") << std::endl;
0215         }
0216 
0217         if (debug_) {
0218           LogDebug("") << "EcalEBPhase2TrigPrimAlgo output of timefinder is a vector of size: " << std::dec
0219                        << time_out_.size() << std::endl;
0220           for (unsigned int ix = 0; ix < time_out_.size(); ix++) {
0221             LogDebug("") << std::dec << time_out_[ix] << " ";
0222           }
0223           LogDebug("") << std::endl;
0224         }
0225 
0226         this->getTPFormatter()->process(filt_out_, time_out_, outEt_, outTime_);
0227 
0228         if (debug_) {
0229           LogDebug("") << " compressed Et "
0230                        << " ";
0231           for (unsigned int iSample = 0; iSample < outEt_.size(); ++iSample) {
0232             LogDebug("") << outEt_[iSample] << " ";
0233           }
0234           LogDebug("") << std::endl;
0235 
0236           LogDebug("") << " compressed time "
0237                        << " ";
0238           for (unsigned int iSample = 0; iSample < outEt_.size(); ++iSample) {
0239             LogDebug("") << outTime_[iSample] << " ";
0240           }
0241           LogDebug("") << std::endl;
0242         }
0243 
0244         if (debug_) {
0245           LogDebug("") << " EcalEBPhase2TrigPrimAlgo  after getting the formatter " << std::endl;
0246           for (unsigned int iSample = 0; iSample < outEt_.size(); ++iSample) {
0247             LogDebug("") << " outEt " << outEt_[iSample] << " outTime " << outTime_[iSample] << " ";
0248           }
0249           LogDebug("") << std::endl;
0250         }
0251 
0252         // } not a spike
0253 
0254         // create the final TP samples
0255         int etInADC = 0;
0256         ;
0257         int64_t time = -999;
0258         int nSam = 0;
0259         for (int iSample = 0; iSample < 16; ++iSample) {
0260           etInADC = outEt_[iSample];
0261           time = outTime_[iSample];
0262           if (debug_) {
0263             LogDebug("") << "TrigPrimAlgo   outEt " << outEt_[iSample] << " outTime " << outTime_[iSample] << std::endl;
0264             LogDebug("") << "TrigPrimAlgo etInADCt " << outEt_[iSample] << " outTime " << time << std::endl;
0265           }
0266 
0267           tp.setSample(nSam, EcalEBPhase2TriggerPrimitiveSample(etInADC, isASpike, time));
0268           nSam++;
0269         }
0270 
0271         result.push_back(tp);
0272 
0273       }  // Loop over the xStals
0274 
0275     }  //loop over strips in one tower
0276 
0277     if (debug_) {
0278       if (nXinBCP > 0)
0279         LogDebug("") << " Accepted xTals " << nXinBCP << std::endl;
0280     }
0281   }
0282 }
0283 
0284 //----------------------------------------------------------------------
0285 
0286 int EcalEBPhase2TrigPrimAlgo::findStripNr(const EBDetId &id) {
0287   int stripnr;
0288   int n = ((id.ic() - 1) % 100) / 20;  //20 corresponds to 4 * ecal_barrel_crystals_per_strip FIXME!!
0289   if (id.ieta() < 0)
0290     stripnr = n + 1;
0291   else
0292     stripnr = nbMaxStrips_ - n;
0293   return stripnr;
0294 }