PreshowerPhiClusterAlgo.cc

Go to the documentation of this file.

#include <vector>
#include <map>
#include <iostream>

#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "RecoEcal/EgammaClusterAlgos/interface/PreshowerPhiClusterAlgo.h"
#include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include "TH1.h"

reco::PreshowerCluster PreshowerPhiClusterAlgo::makeOneCluster(ESDetId strip,
                                                               HitsID* used_strips,
                                                               RecHitsMap* the_rechitsMap_p,
                                                               const CaloSubdetectorGeometry* geometry_p,
                                                               double deltaEta,
                                                               double minDeltaPhi,
                                                               double maxDeltaPhi) {
  rechits_map = the_rechitsMap_p;
  used_s = used_strips;

  int plane = strip.plane();

  // create null-cluster
  std::vector<std::pair<DetId, float> > dummy;
  Point posi(0, 0, 0);
  reco::PreshowerCluster nullcluster = reco::PreshowerCluster(0., posi, dummy, plane);
  if (strip == ESDetId(0))
    return nullcluster;  // works in case of no intersected strip found (e.g. in the Barrel)

  auto refCell = geometry_p->getGeometry(strip);
  const GlobalPoint& refpos = refCell->getPosition();
  double refEta = refpos.eta();
  double refPhi = refpos.phi();

  // Collection of cluster strips
  EcalRecHitCollection clusterRecHits;

  double x_pos = 0;
  double y_pos = 0;
  double z_pos = 0;

  RecHitsMap::iterator strip_it;
  for (strip_it = rechits_map->begin(); strip_it != rechits_map->end(); ++strip_it) {
    if (!goodStrip(strip_it))
      continue;

    ESDetId mystrip = (strip_it->first == DetId(0)) ? ESDetId(0) : ESDetId(strip_it->first);
    if (mystrip.plane() != plane)
      continue;

    auto thisCell = geometry_p->getGeometry(strip_it->first);
    const GlobalPoint& position = thisCell->getPosition();

    if (fabs(position.eta() - refEta) < deltaEta) {
      //std::cout<<"all strips : "<<mystrip.plane()<<" "<<position.phi()<<" "<<reco::deltaPhi(position.phi(), refPhi)<<std::endl;

      if (reco::deltaPhi(position.phi(), refPhi) > 0 && reco::deltaPhi(position.phi(), refPhi) > maxDeltaPhi)
        continue;
      if (reco::deltaPhi(position.phi(), refPhi) < 0 && reco::deltaPhi(position.phi(), refPhi) < minDeltaPhi)
        continue;

      //std::cout<<mystrip.zside()<<" "<<mystrip.plane()<<" "<<mystrip.six()<<" "<<mystrip.siy()<<" "<<mystrip.strip()<<" "<<position.eta()<<" "<<position.phi()<<" "<<strip_it->second.energy()<<" "<<strip_it->second.recoFlag()<<" "<<refEta<<" "<<refPhi<<" "<<reco::deltaPhi(position.phi(), refPhi)<<std::endl;

      clusterRecHits.push_back(strip_it->second);
      used_s->insert(strip_it->first);

      x_pos += strip_it->second.energy() * position.x();
      y_pos += strip_it->second.energy() * position.y();
      z_pos += strip_it->second.energy() * position.z();
    }
  }

  EcalRecHitCollection::iterator it;
  double Eclust = 0;

  std::vector<std::pair<DetId, float> > usedHits;
  for (it = clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
    Eclust += it->energy();
    usedHits.push_back(std::pair<DetId, float>(it->id(), 1.));
  }

  if (Eclust > 0.) {
    x_pos /= Eclust;
    y_pos /= Eclust;
    z_pos /= Eclust;
  }
  Point pos(x_pos, y_pos, z_pos);

  reco::PreshowerCluster cluster = reco::PreshowerCluster(Eclust, pos, usedHits, plane);

  return cluster;
}

// returns true if the candidate strip fulfills the requirements to be added to the cluster:
bool PreshowerPhiClusterAlgo::goodStrip(RecHitsMap::iterator candidate_it) {
  if (used_s->find(candidate_it->first) != used_s->end())
    LogTrace("PreShowerPhiClusterAlgo") << " This strip is in use";
  if (candidate_it == rechits_map->end())
    LogTrace("PreShowerPhiClusterAlgo") << " No such a strip in rechits_map";
  if (candidate_it->second.energy() <= esStripEnergyCut_)
    LogTrace("PreShowerPhiClusterAlgo") << "Strip energy" << candidate_it->second.energy() << "is below threshold";

  // crystal should not be included...
  if ((used_s->find(candidate_it->first) != used_s->end()) ||  // ...if it already belongs to a cluster
      (candidate_it == rechits_map->end()) ||                  // ...if it corresponds to a hit
      (candidate_it->second.energy() <= esStripEnergyCut_))    // ...if it has a negative or zero energy
  {
    return false;
  }
  return true;
}