Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217
#include "Calibration/EcalAlCaRecoProducers/plugins/AlCaECALRecHitReducer.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"

//#define ALLrecHits
//#define DEBUG

//#define QUICK -> if commented loop over the recHits of the SC and add them to the list of recHits to be saved
//                 comment it if you want a faster module but be sure the window is large enough

/** This module reduces the recHitCollections and the caloCalusterCollections in input 
 * keeping only those associated to the given electrons/photons 
 */

/// \todo make sure that the new caloClusterCollection has no duplicates

AlCaECALRecHitReducer::AlCaECALRecHitReducer(const edm::ParameterSet& iConfig) {
  ebRecHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ebRecHitsLabel"));
  eeRecHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("eeRecHitsLabel"));
  //  esRecHitsLabel_ = iConfig.getParameter< edm::InputTag > ("esRecHitsLabel");

  std::vector<edm::InputTag> srcLabels = iConfig.getParameter<std::vector<edm::InputTag> >("srcLabels");
  for (auto inputTag = srcLabels.begin(); inputTag != srcLabels.end(); ++inputTag) {
    eleViewTokens_.push_back(consumes<edm::View<reco::RecoCandidate> >(*inputTag));
  }

  //eleViewToken_ = consumes<edm::View <reco::RecoCandidate> > (iConfig.getParameter< edm::InputTag > ("electronLabel"));
  photonToken_ = consumes<reco::PhotonCollection>(iConfig.getParameter<edm::InputTag>("photonLabel"));
  EESuperClusterToken_ =
      consumes<reco::SuperClusterCollection>(iConfig.getParameter<edm::InputTag>("EESuperClusterCollection"));

  caloTopologyToken_ = esConsumes<CaloTopology, CaloTopologyRecord>();

  minEta_highEtaSC_ = iConfig.getParameter<double>("minEta_highEtaSC");

  alcaBarrelHitsCollection_ = iConfig.getParameter<std::string>("alcaBarrelHitCollection");
  alcaEndcapHitsCollection_ = iConfig.getParameter<std::string>("alcaEndcapHitCollection");
  alcaCaloClusterCollection_ = iConfig.getParameter<std::string>("alcaCaloClusterCollection");

  //  alcaPreshowerHitsCollection_ = iConfig.getParameter<std::string>("alcaPreshowerHitCollection");

  etaSize_ = iConfig.getParameter<int>("etaSize");
  phiSize_ = iConfig.getParameter<int>("phiSize");
  // FIXME: minimum size of etaSize_ and phiSize_
  if (phiSize_ % 2 == 0 || etaSize_ % 2 == 0)
    edm::LogError("AlCaECALRecHitReducerError") << "Size of eta/phi should be odd numbers";

  //  esNstrips_  = iConfig.getParameter<int> ("esNstrips");
  //  esNcolumns_ = iConfig.getParameter<int> ("esNcolumns");

  //register your products
  produces<EBRecHitCollection>(alcaBarrelHitsCollection_);
  produces<EERecHitCollection>(alcaEndcapHitsCollection_);
  produces<reco::CaloClusterCollection>(alcaCaloClusterCollection_);
  //  produces< ESRecHitCollection > (alcaPreshowerHitsCollection_) ;
}

AlCaECALRecHitReducer::~AlCaECALRecHitReducer() {}

// ------------ method called to produce the data  ------------
void AlCaECALRecHitReducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  using namespace edm;
  //using namespace std;

  const CaloTopology* caloTopology = &iSetup.getData(caloTopologyToken_);

  // Get Photons
  Handle<reco::PhotonCollection> phoHandle;
  iEvent.getByToken(photonToken_, phoHandle);

  // get RecHits
  Handle<EBRecHitCollection> barrelRecHitsHandle;
  iEvent.getByToken(ebRecHitsToken_, barrelRecHitsHandle);
  const EBRecHitCollection* barrelHitsCollection = barrelRecHitsHandle.product();

  // get RecHits
  Handle<EERecHitCollection> endcapRecHitsHandle;
  iEvent.getByToken(eeRecHitsToken_, endcapRecHitsHandle);
  const EERecHitCollection* endcapHitsCollection = endcapRecHitsHandle.product();

  //   // get ES RecHits
  //   Handle<ESRecHitCollection> preshowerRecHitsHandle;
  //   iEvent.getByToken(esRecHitsToken_,preshowerRecHitsHandle);

  //   const ESRecHitCollection * preshowerHitsCollection = 0 ;
  //   if (preshowerIsFull)
  //     preshowerHitsCollection = preshowerRecHitsHandle.product () ;

  //   // make a vector to store the used ES rechits:
  //   set<ESDetId> used_strips;
  //   used_strips.clear();

  // for Z->ele+SC
  Handle<reco::SuperClusterCollection> EESCHandle;
  iEvent.getByToken(EESuperClusterToken_, EESCHandle);

  //Create empty output collections
  auto miniEBRecHitCollection = std::make_unique<EBRecHitCollection>();
  auto miniEERecHitCollection = std::make_unique<EERecHitCollection>();
  //  auto miniESRecHitCollection = std::make_unique<ESRecHitCollection>();

  std::set<DetId> reducedRecHit_EBmap;
  std::set<DetId> reducedRecHit_EEmap;

  //  std::set< edm::Ref<reco::CaloCluster> > reducedCaloClusters_map;

  auto reducedCaloClusterCollection = std::make_unique<reco::CaloClusterCollection>();

  //Photons:
#ifdef shervin
  for (reco::PhotonCollection::const_iterator phoIt = phoHandle->begin(); phoIt != phoHandle->end(); phoIt++) {
    const reco::SuperCluster& sc = *(phoIt->superCluster());

    if (phoIt->isEB()) {
      AddMiniRecHitCollection(sc, reducedRecHit_EBmap, caloTopology);
    } else {  // endcap
      AddMiniRecHitCollection(sc, reducedRecHit_EEmap, caloTopology);
    }  // end of endcap

    /// \todo check if this works when you ask sc->seed(), I suspect that the references have to be updated
    reco::CaloCluster_iterator it = sc.clustersBegin();
    reco::CaloCluster_iterator itend = sc.clustersEnd();
    for (; it != itend; ++it) {
      reco::CaloCluster caloClus(**it);
      reducedCaloClusterCollection->push_back(caloClus);
    }
  }
#endif

  Handle<edm::View<reco::RecoCandidate> > eleViewHandle;
  for (auto iToken = eleViewTokens_.begin(); iToken != eleViewTokens_.end(); iToken++) {
    iEvent.getByToken(*iToken, eleViewHandle);

    //Electrons:
    for (auto eleIt = eleViewHandle->begin(); eleIt != eleViewHandle->end(); eleIt++) {
      const reco::SuperCluster& sc = *(eleIt->superCluster());

      if (fabs(sc.eta()) < 1.479) {
        AddMiniRecHitCollection(sc, reducedRecHit_EBmap, caloTopology);
      } else {  // endcap
        AddMiniRecHitCollection(sc, reducedRecHit_EEmap, caloTopology);
      }  // end of endcap

      reco::CaloCluster_iterator it = sc.clustersBegin();
      reco::CaloCluster_iterator itend = sc.clustersEnd();
      for (; it != itend; ++it) {
        reco::CaloCluster caloClus(**it);
        reducedCaloClusterCollection->push_back(caloClus);
      }
    }
  }

  //saving recHits for highEta SCs for highEta studies
  for (reco::SuperClusterCollection::const_iterator SC_iter = EESCHandle->begin(); SC_iter != EESCHandle->end();
       SC_iter++) {
    if (fabs(SC_iter->eta()) < minEta_highEtaSC_)
      continue;
    AddMiniRecHitCollection(*SC_iter, reducedRecHit_EEmap, caloTopology);

    const reco::SuperCluster& sc = *(SC_iter);
    reco::CaloCluster_iterator it = sc.clustersBegin();
    reco::CaloCluster_iterator itend = sc.clustersEnd();
    for (; it != itend; ++it) {
      reco::CaloCluster caloClus(**it);
      reducedCaloClusterCollection->push_back(caloClus);
    }
  }

  //------------------------------ fill the alcareco reduced recHit collection
  for (std::set<DetId>::const_iterator itr = reducedRecHit_EBmap.begin(); itr != reducedRecHit_EBmap.end(); itr++) {
    if (barrelHitsCollection->find(*itr) != barrelHitsCollection->end())
      miniEBRecHitCollection->push_back(*(barrelHitsCollection->find(*itr)));
  }

  for (std::set<DetId>::const_iterator itr = reducedRecHit_EEmap.begin(); itr != reducedRecHit_EEmap.end(); itr++) {
    if (endcapHitsCollection->find(*itr) != endcapHitsCollection->end())
      miniEERecHitCollection->push_back(*(endcapHitsCollection->find(*itr)));
  }

  //--------------------------------------- Put selected information in the event
  iEvent.put(std::move(miniEBRecHitCollection), alcaBarrelHitsCollection_);
  iEvent.put(std::move(miniEERecHitCollection), alcaEndcapHitsCollection_);
  //  iEvent.put(std::move(miniESRecHitCollection),alcaPreshowerHitsCollection_ );
  iEvent.put(std::move(reducedCaloClusterCollection), alcaCaloClusterCollection_);
}

void AlCaECALRecHitReducer::AddMiniRecHitCollection(const reco::SuperCluster& sc,
                                                    std::set<DetId>& reducedRecHitMap,
                                                    const CaloTopology* caloTopology) const {
  DetId seed = (sc.seed()->seed());
  int phiSize = phiSize_, etaSize = etaSize_;
  if (seed.subdetId() != EcalBarrel) {  // if not EB, take a square window
    etaSize = std::max(phiSize_, etaSize_);
    phiSize = etaSize;
  }

  std::vector<DetId> recHit_window = caloTopology->getWindow(seed, phiSize, etaSize);
  for (unsigned int i = 0; i < recHit_window.size(); i++) {
    reducedRecHitMap.insert(recHit_window[i]);
  }

  const std::vector<std::pair<DetId, float> >& scHits = sc.hitsAndFractions();
  for (std::vector<std::pair<DetId, float> >::const_iterator scHit_itr = scHits.begin(); scHit_itr != scHits.end();
       scHit_itr++) {
    // the map fills just one time (avoiding double insert of recHits)
    reducedRecHitMap.insert(scHit_itr->first);
  }

  return;
}

DEFINE_FWK_MODULE(AlCaECALRecHitReducer);