text: #include "DataFormats/EcalDetId/interface/EBDetId.h"
  #include "FWCore/MessageLogger/interface/MessageLogger.h"
  #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
  #include "EcalBasicClusterLocalContCorrection.h"
  
  #include "TVector2.h"
  
  EcalBasicClusterLocalContCorrection::EcalBasicClusterLocalContCorrection(edm::ConsumesCollector &&cc)
      : paramsToken_{cc.esConsumes()}, caloGeometryToken_{cc.esConsumes()} {}
  
  void EcalBasicClusterLocalContCorrection::init(const edm::EventSetup &es) {
    params_ = &es.getData(paramsToken_);
    caloGeometry_ = &es.getData(caloGeometryToken_);
  }
  
  void EcalBasicClusterLocalContCorrection::checkInit() const {
    if (!params_) {
      // non-initialized function parameters: throw exception
      throw cms::Exception("EcalBasicClusterLocalContCorrection::checkInit()")
          << "Trying to access an uninitialized crack correction function.\n"
             "Please call `init( edm::EventSetup &)' before any use of the function.\n";
    }
  }
  
  using namespace std;
  using namespace edm;
  
  float EcalBasicClusterLocalContCorrection::operator()(const reco::BasicCluster &basicCluster,
                                                        const EcalRecHitCollection &recHit) const {
    checkInit();
  
    // number of parameters needed by this parametrization
    constexpr size_t nparams = 24;
  
    //correction factor to be returned, and to be calculated in this present function:
    double correction_factor = 1.;
    double fetacor = 1.;  //eta dependent part of the correction factor
    double fphicor = 1.;  //phi dependent part of the correction factor
  
    //--------------if barrel calculate local position wrt xtal center -------------------
    const CaloSubdetectorGeometry *geom =
        caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);  //EcalBarrel = 1
  
    const math::XYZPoint &position_ = basicCluster.position();
    double Theta = -position_.theta() + 0.5 * M_PI;
    double Eta = position_.eta();
    double Phi = TVector2::Phi_mpi_pi(position_.phi());
  
    //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
    // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
    constexpr float X0 = 0.89;
    constexpr float T0 = 7.4;
    double depth = X0 * (T0 + log(basicCluster.energy()));
  
    //search which crystal is closest to the cluster position and call it crystalseed:
    //std::vector<DetId> crystals_vector = *scRef.getHitsByDetId();   //deprecated
    std::vector<std::pair<DetId, float> > crystals_vector = basicCluster.hitsAndFractions();
    float dphimin = 999.;
    float detamin = 999.;
    int ietaclosest = 0;
    int iphiclosest = 0;
  
    for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
      EBDetId crystal(crystals_vector[icry].first);
      auto cell = geom->getGeometry(crystal);  // problema qui
      GlobalPoint center_pos = cell->getPosition(depth);
      double EtaCentr = center_pos.eta();
      double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
      if (std::abs(EtaCentr - Eta) < detamin) {
        detamin = std::abs(EtaCentr - Eta);
        ietaclosest = crystal.ieta();
      }
      if (std::abs(TVector2::Phi_mpi_pi(PhiCentr - Phi)) < dphimin) {
        dphimin = std::abs(TVector2::Phi_mpi_pi(PhiCentr - Phi));
        iphiclosest = crystal.iphi();
      }
    }
  
    EBDetId crystalseed(ietaclosest, iphiclosest);
  
    // Get center cell position from shower depth
    auto cell = geom->getGeometry(crystalseed);
    GlobalPoint center_pos = cell->getPosition(depth);
  
    //PHI
    double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
    double PhiWidth = (M_PI / 180.);
    double PhiCry = (TVector2::Phi_mpi_pi(Phi - PhiCentr)) / PhiWidth;
    if (PhiCry > 0.5)
      PhiCry = 0.5;
    if (PhiCry < -0.5)
      PhiCry = -0.5;
    //flip to take into account ECAL barrel symmetries:
    if (ietaclosest < 0)
      PhiCry *= -1.;
  
    //ETA
    double ThetaCentr = -center_pos.theta() + 0.5 * M_PI;
    double ThetaWidth = (M_PI / 180.) * std::cos(ThetaCentr);
    double EtaCry = (Theta - ThetaCentr) / ThetaWidth;
    if (EtaCry > 0.5)
      EtaCry = 0.5;
    if (EtaCry < -0.5)
      EtaCry = -0.5;
    //flip to take into account ECAL barrel symmetries:
    if (ietaclosest < 0)
      EtaCry *= -1.;
  
    //-------------- end calculate local position -------------
  
    size_t payloadsize = params_->params().size();
  
    if (payloadsize < nparams)
      edm::LogError("Invalid Payload") << "Parametrization requires " << nparams << " parameters but only " << payloadsize
                                       << " are found in DB. Perhaps incompatible Global Tag" << std::endl;
  
    if (payloadsize > nparams)
      edm::LogWarning("Size mismatch ") << "Parametrization requires " << nparams << " parameters but " << payloadsize
                                        << " are found in DB. Perhaps incompatible Global Tag" << std::endl;
  
    std::pair<double, double> localPosition(EtaCry, PhiCry);
  
    //--- local cluster coordinates
    float localEta = localPosition.first;
    float localPhi = localPosition.second;
  
    //--- ecal module
    int imod = getEcalModule(basicCluster.seed());
  
    //-- corrections parameters
    float pe[3], pp[3];
    pe[0] = (params_->params())[0 + imod * 3];
    pe[1] = (params_->params())[1 + imod * 3];
    pe[2] = (params_->params())[2 + imod * 3];
    pp[0] = (params_->params())[12 + imod * 3];
    pp[1] = (params_->params())[13 + imod * 3];
    pp[2] = (params_->params())[14 + imod * 3];
  
    //--- correction vs local eta
    fetacor = pe[0] + pe[1] * localEta + pe[2] * localEta * localEta;
  
    //--- correction vs local phi
    fphicor = pp[0] + pp[1] * localPhi + pp[2] * localPhi * localPhi;
  
    //if the seed crystal is neighbourgh of a supermodule border, don't apply the phi dependent  containment corrections, but use the larger crack corrections instead.
    int iphimod20 = std::abs(iphiclosest % 20);
    if (iphimod20 <= 1)
      fphicor = 1.;
  
    correction_factor = (1. / fetacor) * (1. / fphicor);
  
    //return the correction factor. Use it to multiply the cluster energy.
    return correction_factor;
  }
  
  //------------------------------------------------------------------------------------------------------
  int EcalBasicClusterLocalContCorrection::getEcalModule(DetId id) const {
    int mod = 0;
    int ieta = (EBDetId(id)).ieta();
  
    if (fabs(ieta) <= 25)
      mod = 0;
    if (fabs(ieta) > 25 && fabs(ieta) <= 45)
      mod = 1;
    if (fabs(ieta) > 45 && fabs(ieta) <= 65)
      mod = 2;
    if (fabs(ieta) > 65 && fabs(ieta) <= 85)
      mod = 3;
  
    return (mod);
  }
  