EcalSelectiveReadoutValidation.h

Go to the documentation of this file.

#ifndef EcalSelectiveReadoutValidation_H
#define EcalSelectiveReadoutValidation_H

/*
 * \file EcalSelectiveReadoutValidation.h
 *
 * $Date: 2012/01/31 22:28:20 $
 * $Revision: 1.14 $
 *
 */

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "SimDataFormats/CaloHit/interface/PCaloHit.h"


#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
#include "DataFormats/EcalDetId/interface/EcalScDetId.h"
#include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"

#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
//#include "DataFormats/EcalDetId/interface/EEDetId.h"
//#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/FEDRawData/interface/FEDRawDataCollection.h"

#include "Validation/EcalDigis/src/CollHandle.h"

#include <string>
#include <set>
#include <utility>
#include <fstream>
#include <inttypes.h>
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"

class DQMStore;
class MonitorElement;
class EBDetId;
class EEDetId;
class EcalTPParameters;
class EcalElectronicsMapping;

class EcalSelectiveReadoutValidation: public edm::EDAnalyzer{

  //typedef EcalUncalibratedRecHitCollection RecHitCollection;
  //typedef EcalUncalibratedRecHit RecHit;
  typedef EcalRecHitCollection RecHitCollection;
  typedef EcalRecHit RecHit;
  
public:

  EcalSelectiveReadoutValidation(const edm::ParameterSet& ps);

  ~EcalSelectiveReadoutValidation();

protected:

  void analyze(edm::Event const & e, edm::EventSetup const & c);

  void beginRun(const edm::Run& r, const edm::EventSetup& c); 
  
  void endRun(const edm::Run& r, const edm::EventSetup& c); 

private:
  enum subdet_t {EB, EE};

  template<class T, class U>
  void anaDigi(const T& frame, const U& srFlagColl);

  void anaDigiInit();

  void analyzeDataVolume(const edm::Event& e, const edm::EventSetup& es);

  void analyzeEB(const edm::Event& event, const edm::EventSetup& es);

  void analyzeEE(const edm::Event& event, const edm::EventSetup& es);

  void analyzeTP(const edm::Event& event, const edm::EventSetup& es); 
  
  //  /** Selective Readout decisions Validation
  //    * @param event EDM event
  //    * @param es event setup
  //    */ 
  //   void SRFlagValidation(const edm::Event& event, const edm::EventSetup& es);

  void SRFlagValidation(const edm::Event& event, const edm::EventSetup& es);

  double frame2Energy(const EcalDataFrame& frame) const;

  template<class T>
  double frame2EnergyForTp(const T& frame, int offset = 0) const;
  
  //   double getEcalEventSize(double nReadXtals) const{
  //     return getDccOverhead(EB)*nEbDccs+getDccOverhead(EE)*nEeDccs
  //       + nReadXtals*getBytesPerCrystal()
  //       + (nEeRus+nEbRus)*8;
  //   }

  double getEbEventSize(double nReadXtals) const;

  double getEeEventSize(double nReadXtals) const;

  double getDccOverhead(subdet_t subdet) const{
    //  return (subdet==EB?34:25)*8;
    return (subdet==EB?34:52)*8;
  }

  double getBytesPerCrystal() const{
    return 3*8;
  }

  double getDccEventSize(int iDcc0, double nReadXtals) const{
    subdet_t subdet;
     if(iDcc0<9 || iDcc0>=45){
       subdet = EE;
     } else{
       subdet = EB;
     }
     //     return getDccOverhead(subdet)+nReadXtals*getBytesPerCrystal()
     //       + getRuCount(iDcc0)*8;
    return getDccOverhead(subdet)
      + getDccSrDependentPayload(iDcc0, getRuCount(iDcc0), nReadXtals);
  }

  double getDccSrDependentPayload(int iDcc0, double nReadRus,
                                  double nReadXtals) const{
    return nReadXtals*getBytesPerCrystal() + nReadRus*8;
  }

  int getRuCount(int iDcc0) const;

  void readAllCollections(const edm::Event& e);

  void setTtEtSums(const edm::EventSetup& es,
                   const EBDigiCollection& ebDigis,
                   const EEDigiCollection& eeDigis);

  //   /** Retrieves the logical number of the DCC reading a given crystal channel.
  //    * @param xtarId crystal channel identifier
  //    * @return the DCC logical number starting from 1.
  //    */
  //   unsigned dccNum(const DetId& xtalId) const;

  std::pair<int,int> dccCh(const DetId& xtalId) const;

  
  int iEta2cIndex(int iEta) const{
    return (iEta<0)?iEta+85:iEta+84;
  }
  
  int iPhi2cIndex(int iPhi) const{
    //    return iPhi-1;
    int iPhi0 = iPhi - 11;
    if(iPhi0<0) iPhi0 += 360;
    return iPhi0;
  }

  int iXY2cIndex(int iX) const{
    return iX-1;
  }

  int cIndex2iXY(int iX0) const{
    return iX0+1;
  }

  int cIndex2iEta(int i) const{
    return (i<85)?i-85:i-84;
  }


  int cIndex2iPhi(int i) const {
    return (i+11) % 360;
  }

  int iTtEta2cIndex(int iEta) const{
    return (iEta<0)?iEta+28:iEta+27;
  }
  
  int iTtPhi2cIndex(int iPhi) const{
    return iPhi-1;
    //int iPhi0 = iPhi - 3;
    //if(iPhi0<0) iPhi0 += 72;
    //return iPhi0;
  }

  int cIndex2iTtEta(int i) const{
    return (i<27)?i-28:i-27;
  }

  int cIndex2iTtPhi(int i) const{
    return i + 1;
  }
  
  EcalScDetId superCrystalOf(const EEDetId& xtalId) const;


  EcalTrigTowerDetId readOutUnitOf(const EBDetId& xtalId) const;
  
  EcalScDetId readOutUnitOf(const EEDetId& xtalId) const;

  static int dccZsFIR(const EcalDataFrame& frame,
              const std::vector<int>& firWeights,
              int firstFIRSample,
              bool* saturated = 0);
  

  static std::vector<int> getFIRWeights(const std::vector<double>&
                    normalizedWeights);
  

  MonitorElement* bookFloat(const std::string& name);
  
  MonitorElement* book1D(const std::string& name,
             const std::string& title, int nbins,
             double xmin, double xmax);
 
  MonitorElement* book2D(const std::string& name,
             const std::string& title,
             int nxbins, double xmin, double xmax,
             int nybins, double ymin, double ymax);

  MonitorElement* bookProfile(const std::string& name,
                  const std::string& title,
                  int nbins, double xmin, double xmax);
  
  MonitorElement* bookProfile2D(const std::string& name,
                const std::string& title,
                int nbinx, double xmin, double xmax,
                int nbiny, double ymin, double ymax,
                const char* option = "");


  void fill(MonitorElement* me, float x){
    if(me) me->Fill(x);
  }
  void fill(MonitorElement* me,float x, float yw){
    if(me) me->Fill(x, yw);
  }
  void fill(MonitorElement* me,float x, float y, float zw){
    if(me) me->Fill(x, y, zw);
  }
  void fill(MonitorElement* me,float x, float y, float z, float w){
    if(me) me->Fill(x, y, z, w);
  }

  void initAsciiFile();

  void updateL1aRate(const edm::Event& event);

  double getL1aRate() const;
  
private:
  struct energiesEb_t{
    double simE;     
    double noZsRecE; 
    double recE;     
    //    EBDigiCollection::const_iterator itNoZsFrame; //
    int simHit;      
    double phi;      
    double eta;      
    bool gain12;     //all MGPA samples at gain 12?
  };

  struct energiesEe_t{
    double simE;     
    double noZsRecE; 
    double recE;     
    //    EEDigiCollection::const_iterator itNoZsFrame;
    int simHit;      
    double phi;      
    double eta;      
    bool gain12;     //all MGPA samples at gain 12?
  };

  static const int kByte_ = 1024;
  
  static const unsigned nDccs_ = 54;

  // = maximum number of RUs read by a DCC
  static const unsigned nDccChs_ = 68;
  
  //Lower bound of DCC ID range
  static const int minDccId_ = 1 ;

  //Upper bound of DCC ID range
  static const int maxDccId_  = minDccId_ + nDccs_ -1;
  
  static const int nEbDccs = 36;

  static const int nEeDccs = 18;

  static const int nEbRus = 36*68;

  static const int nEeRus = 2*(34+32+33+33+32+34+33+34+33);

  static const int nDccRus_[nDccs_];

  static const int nEndcaps = 2;

  static const int nEbEta = 170;

  static const int nEbPhi = 360;

  static const int nEeX = 100;

  static const int nEeY = 100;

  static const int ebTtEdge = 5;
  
  static const int scEdge = 5;

  static const int nOneEeTtEta = 11;
  
  static const int nEbTtEta = 34;
  
  static const int nTtEta = 2*nOneEeTtEta + nEbTtEta;

  static const int nTtPhi = 72;

  static const int nMaxXtalPerRu = 25;
  static const double rad2deg;
  
  bool verbose_;

  DQMStore* dbe_;

  std::string outputFile_;

  bool collNotFoundWarn_;
  
  std::string srpAlgoErrorLogFileName_;

  std::string srApplicationErrorLogFileName_;
  
  std::ofstream srpAlgoErrorLog_; 

  std::ofstream srApplicationErrorLog_; 

  std::ofstream zsErrorLog_;
  

  CollHandle<EBDigiCollection>           ebDigis_;
  CollHandle<EEDigiCollection>           eeDigis_;
  CollHandle<EBDigiCollection>           ebNoZsDigis_;
  CollHandle<EEDigiCollection>           eeNoZsDigis_;
  CollHandle<EBSrFlagCollection>         ebSrFlags_;
  CollHandle<EESrFlagCollection>         eeSrFlags_;
  CollHandle<EBSrFlagCollection>         ebComputedSrFlags_;
  CollHandle<EESrFlagCollection>         eeComputedSrFlags_;
  CollHandle<std::vector<PCaloHit> >     ebSimHits_;
  CollHandle<std::vector<PCaloHit> >     eeSimHits_;
  CollHandle<EcalTrigPrimDigiCollection> tps_;
  CollHandle<RecHitCollection>       ebRecHits_;
  CollHandle<RecHitCollection>       eeRecHits_;
  CollHandle<FEDRawDataCollection>       fedRaw_;


  int64_t tmax;
  int64_t tmin;
  int64_t l1aOfTmin;
  int64_t l1aOfTmax;
  bool l1aRateErr;
  

  MonitorElement* meDccVol_;
  MonitorElement* meDccLiVol_;
  MonitorElement* meDccHiVol_;
  MonitorElement* meDccVolFromData_;
  MonitorElement* meVol_;
  MonitorElement* meVolB_;
  MonitorElement* meVolE_;
  MonitorElement* meVolBLI_;
  MonitorElement* meVolELI_;
  MonitorElement* meVolLI_;
  MonitorElement* meVolBHI_;
  MonitorElement* meVolEHI_;
  MonitorElement* meVolHI_;
  MonitorElement* meChOcc_;
  
  MonitorElement* meTp_;
  MonitorElement* meTtf_;
  MonitorElement* meTtfVsTp_;
  MonitorElement* meTtfVsEtSum_;
  MonitorElement* meTpVsEtSum_;

  MonitorElement* meEbRecE_;
  MonitorElement* meEbEMean_;
  MonitorElement* meEbNoise_;
  MonitorElement* meEbSimE_;
  MonitorElement* meEbRecEHitXtal_;
  MonitorElement* meEbRecVsSimE_;
  MonitorElement* meEbNoZsRecVsSimE_;
  
  MonitorElement* meEeRecE_;
  MonitorElement* meEeEMean_;
  MonitorElement* meEeNoise_;
  MonitorElement* meEeSimE_;
  MonitorElement* meEeRecEHitXtal_;
  MonitorElement* meEeRecVsSimE_;
  MonitorElement* meEeNoZsRecVsSimE_;

  MonitorElement* meFullRoRu_;
  MonitorElement* meZs1Ru_;
  MonitorElement* meForcedRu_;

  MonitorElement* meLiTtf_;
  MonitorElement* meMiTtf_;
  MonitorElement* meHiTtf_;
  MonitorElement* meForcedTtf_;

  MonitorElement* meTpMap_;

  MonitorElement* meFullRoCnt_;
  MonitorElement* meEbFullRoCnt_;
  MonitorElement* meEeFullRoCnt_;
  
  MonitorElement* meEbLiZsFir_;
  MonitorElement* meEbHiZsFir_;
  MonitorElement* meEbIncompleteRUZsFir_;
  
  MonitorElement* meEeLiZsFir_;
  MonitorElement* meEeHiZsFir_;
  MonitorElement* meSRFlagsFromData_; 
  MonitorElement* meSRFlagsComputed_; 
  MonitorElement* meSRFlagsConsistency_;
    
  MonitorElement* meIncompleteFRO_;
  MonitorElement* meDroppedFRO_;
  MonitorElement* meCompleteZS_;
  
  MonitorElement* meIncompleteFROMap_;
  MonitorElement* meDroppedFROMap_;
  MonitorElement* meCompleteZSMap_;

  MonitorElement* meIncompleteFRORateMap_;
  MonitorElement* meDroppedFRORateMap_;
  MonitorElement* meCompleteZSRateMap_;

  MonitorElement* meIncompleteFROCnt_;
  MonitorElement* meDroppedFROCnt_;
  MonitorElement* meCompleteZSCnt_;
  MonitorElement* meEbZsErrCnt_;
  MonitorElement* meEeZsErrCnt_;
  MonitorElement* meZsErrCnt_;
  MonitorElement* meEbZsErrType1Cnt_;
  MonitorElement* meEeZsErrType1Cnt_;
  MonitorElement* meZsErrType1Cnt_;


  MonitorElement* meEbFixedPayload_;
  MonitorElement* meEeFixedPayload_;
  MonitorElement* meFixedPayload_;
  
  MonitorElement* meL1aRate_;
  
  int nDroppedFRO_;

  int nIncompleteFRO_;

  int nCompleteZS_;

  int nEbFROCnt_;

  int nEeFROCnt_;

  int nEbZsErrors_;

  int nEeZsErrors_;

  int nEbZsErrorsType1_;

  int nEeZsErrorsType1_;

  
  const EcalTrigTowerConstituentsMap * triggerTowerMap_;

  const EcalTPParameters* tpParam_;

  const EcalElectronicsMapping* elecMap_;
  
  bool localReco_;
  
  std::vector<double> weights_;

  std::vector<int> firWeights_;

  int ebZsThr_;

  int eeZsThr_;
  
  bool tpInGeV_;

  int firstFIRSample_;

  bool useEventRate_;

  std::vector<bool> logErrForDccs_;
  
  int nEb_;

  int nEe_;

  int nEeLI_;

  int nEeHI_;

  int nEbLI_;

  int nEbHI_;

  int nPerDcc_[nDccs_];

  int nLiPerDcc_[nDccs_];

  int nHiPerDcc_[nDccs_];

  
  int nRuPerDcc_[nDccs_];

  int nLiRuPerDcc_[nDccs_];

  int nHiRuPerDcc_[nDccs_];

  

  bool ebRuActive_[nEbEta/ebTtEdge][nEbPhi/ebTtEdge];
  bool eeRuActive_[nEndcaps][nEeX/scEdge][nEeY/scEdge];
  
  bool isRuComplete_[nDccs_][nDccChs_];

  int nPerRu_[nDccs_][nDccChs_];  
  
  int ievt_;

  double ttEtSums[nTtEta][nTtPhi];

  energiesEb_t ebEnergies[nEbEta][nEbPhi];

  energiesEe_t eeEnergies[nEndcaps][nEeX][nEeY];

  bool SkipInnerSC_;

  std::set<std::string> histList_;

  bool allHists_;

  std::string histDir_;
  
  std::map<std::string, std::string> availableHistList_;


  bool withEeSimHit_;

  bool withEbSimHit_;
  
  bool registerHist(const std::string& name, const std::string& title);

  void printAvailableHists();

  void normalizeHists(double eventCount);

  void configFirWeights(std::vector<double> weightsForZsFIR);

  bool logSrpAlgoErrors_;

  bool logSrApplicationErrors_;

  template<class T>
  void compareSrfColl(const edm::Event& event, T& srfFromData, T& computedSrf);

  template<class T>
  void checkSrApplication(const edm::Event& event, T& srfs);

  
  int ruGraphX(const EcalScDetId& id) const{
    return id.ix() + (id.zside()>0?20:-40);
  }
  
  int ruGraphY(const EcalScDetId& id) const{
    return id.iy();
  }
  
  int ruGraphX(const EcalTrigTowerDetId& id) const{
    return id.ieta();
  }
  
  int ruGraphY(const EcalTrigTowerDetId& id) const{
    return id.iphi();
  }
  
  int xtalGraphX(const EEDetId& id) const{
    return id.ix() + (id.zside()>0?100:-200);
  }
  
  int xtalGraphY(const EEDetId& id) const{
    return id.iy();
  }
  
  int xtalGraphX(const EBDetId& id) const{
    return id.ieta();
  }
  
  int xtalGraphY(const EBDetId& id) const{
    return id.iphi();
  }


  int dccId(const EcalScDetId& detId) const;
  int dccId(const EcalTrigTowerDetId& detId) const;

  void selectFedsForLog();

  int getCrystalCount(int iDcc, int iDccCh);

  
private:
  class Sorter{
    EcalSelectiveReadoutValidation* validation;
  public:
    Sorter(EcalSelectiveReadoutValidation* v): validation(v){};
    bool operator()(std::pair<int,int>a, std::pair<int,int>b){
      return (validation->ebEnergies[a.first][a.second].simE
          > validation->ebEnergies[b.first][b.second].simE);
    }
  };

  void myAna();
};

#endif //EcalSelectiveReadoutValidation_H not defined