EvtPlaneProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    EvtPlaneProducer
// Class:      EvtPlaneProducer
//
//
// Original Author:  Sergey Petrushanko
//         Created:  Fri Jul 11 10:05:00 2008
//
//

// system include files
#include <memory>
#include <iostream>
#include <ctime>
#include <cmath>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"

#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/Common/interface/Ref.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/CastorReco/interface/CastorTower.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include <cstdlib>
#include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h"
#include "CondFormats/HIObjects/interface/RPFlatParams.h"
#include "CondFormats/DataRecord/interface/HeavyIonRPRcd.h"
#include "CondFormats/DataRecord/interface/HeavyIonRcd.h"
#include "CondFormats/HIObjects/interface/CentralityTable.h"

#include "DataFormats/Math/interface/Point3D.h"
#include "DataFormats/Common/interface/RefProd.h"
#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/Common/interface/RefVector.h"

#include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h"
#include "RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h"

using namespace std;
using namespace hi;

//
// class decleration
//

namespace hi {
  class GenPlane {
  public:
    GenPlane(string name, double etaminval1, double etamaxval1, double etaminval2, double etamaxval2, int orderval) {
      epname = name;
      etamin1 = etaminval1;
      etamax1 = etamaxval1;
      etamin2 = etaminval2;
      etamax2 = etamaxval2;
      sumsin = 0;
      sumcos = 0;
      sumsinNoWgt = 0;
      sumcosNoWgt = 0;

      mult = 0;
      order = (double)orderval;
    }
    ~GenPlane() { ; }
    void addParticle(double w, double PtOrEt, double s, double c, double eta) {
      if ((eta >= etamin1 && eta < etamax1) || (etamin2 != etamax2 && eta >= etamin2 && eta < etamax2)) {
        sumsin += w * s;
        sumcos += w * c;
        sumsinNoWgt += s;
        sumcosNoWgt += c;

        sumw += fabs(w);
        sumw2 += w * w;
        sumPtOrEt += PtOrEt;
        sumPtOrEt2 += PtOrEt * PtOrEt;
        ++mult;
      }
    }

    double getAngle(double &ang,
                    double &sv,
                    double &cv,
                    double &svNoWgt,
                    double &cvNoWgt,
                    double &w,
                    double &w2,
                    double &PtOrEt,
                    double &PtOrEt2,
                    uint &epmult) {
      ang = -10;
      sv = 0;
      cv = 0;
      sv = sumsin;
      cv = sumcos;
      svNoWgt = sumsinNoWgt;
      cvNoWgt = sumcosNoWgt;
      w = sumw;
      w2 = sumw2;
      PtOrEt = sumPtOrEt;
      PtOrEt2 = sumPtOrEt2;
      epmult = mult;
      double q = sv * sv + cv * cv;
      if (q > 0)
        ang = atan2(sv, cv) / order;
      return ang;
    }
    void reset() {
      sumsin = 0;
      sumcos = 0;
      sumsinNoWgt = 0;
      sumcosNoWgt = 0;
      sumw = 0;
      sumw2 = 0;
      mult = 0;
      sumPtOrEt = 0;
      sumPtOrEt2 = 0;
    }

  private:
    string epname;
    double etamin1;
    double etamax1;

    double etamin2;
    double etamax2;
    double sumsin;
    double sumcos;
    double sumsinNoWgt;
    double sumcosNoWgt;
    uint mult;
    double sumw;
    double sumw2;
    double sumPtOrEt;
    double sumPtOrEt2;
    double order;
  };
}  // namespace hi

class EvtPlaneProducer : public edm::stream::EDProducer<> {
public:
  explicit EvtPlaneProducer(const edm::ParameterSet &);
  ~EvtPlaneProducer() override;

private:
  GenPlane *rp[NumEPNames];

  void produce(edm::Event &, const edm::EventSetup &) override;

  // ----------member data ---------------------------

  std::string centralityVariable_;
  std::string centralityLabel_;
  std::string centralityMC_;

  edm::InputTag centralityBinTag_;
  edm::EDGetTokenT<int> centralityBinToken;

  edm::InputTag vertexTag_;
  edm::EDGetTokenT<std::vector<reco::Vertex>> vertexToken;
  edm::Handle<std::vector<reco::Vertex>> vertex_;

  edm::InputTag caloTag_;
  edm::EDGetTokenT<CaloTowerCollection> caloToken;
  edm::Handle<CaloTowerCollection> caloCollection_;

  edm::InputTag castorTag_;
  edm::EDGetTokenT<std::vector<reco::CastorTower>> castorToken;
  edm::Handle<std::vector<reco::CastorTower>> castorCollection_;

  edm::InputTag trackTag_;
  edm::EDGetTokenT<reco::TrackCollection> trackToken;
  edm::Handle<reco::TrackCollection> trackCollection_;

  edm::ESWatcher<HeavyIonRcd> hiWatcher;
  edm::ESWatcher<HeavyIonRPRcd> hirpWatcher;

  bool loadDB_;
  double minet_;
  double maxet_;
  double minpt_;
  double maxpt_;
  double minvtx_;
  double maxvtx_;
  double dzerr_;
  double chi2_;
  int FlatOrder_;
  int NumFlatBins_;
  double nCentBins_;
  double caloCentRef_;
  double caloCentRefWidth_;
  int CentBinCompression_;
  HiEvtPlaneFlatten *flat[NumEPNames];
};

EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig)
    : centralityVariable_(iConfig.getParameter<std::string>("centralityVariable")),
      centralityBinTag_(iConfig.getParameter<edm::InputTag>("centralityBinTag")),
      vertexTag_(iConfig.getParameter<edm::InputTag>("vertexTag")),
      caloTag_(iConfig.getParameter<edm::InputTag>("caloTag")),
      castorTag_(iConfig.getParameter<edm::InputTag>("castorTag")),
      trackTag_(iConfig.getParameter<edm::InputTag>("trackTag")),
      loadDB_(iConfig.getParameter<bool>("loadDB")),
      minet_(iConfig.getParameter<double>("minet")),
      maxet_(iConfig.getParameter<double>("maxet")),
      minpt_(iConfig.getParameter<double>("minpt")),
      maxpt_(iConfig.getParameter<double>("maxpt")),
      minvtx_(iConfig.getParameter<double>("minvtx")),
      maxvtx_(iConfig.getParameter<double>("maxvtx")),
      dzerr_(iConfig.getParameter<double>("dzerr")),
      chi2_(iConfig.getParameter<double>("chi2")),
      FlatOrder_(iConfig.getParameter<int>("FlatOrder")),
      NumFlatBins_(iConfig.getParameter<int>("NumFlatBins")),
      caloCentRef_(iConfig.getParameter<double>("caloCentRef")),
      caloCentRefWidth_(iConfig.getParameter<double>("caloCentRefWidth")),
      CentBinCompression_(iConfig.getParameter<int>("CentBinCompression")) {
  nCentBins_ = 200.;

  if (iConfig.exists("nonDefaultGlauberModel")) {
    centralityMC_ = iConfig.getParameter<std::string>("nonDefaultGlauberModel");
  }
  centralityLabel_ = centralityVariable_ + centralityMC_;

  centralityBinToken = consumes<int>(centralityBinTag_);

  vertexToken = consumes<std::vector<reco::Vertex>>(vertexTag_);

  caloToken = consumes<CaloTowerCollection>(caloTag_);

  castorToken = consumes<std::vector<reco::CastorTower>>(castorTag_);

  trackToken = consumes<reco::TrackCollection>(trackTag_);

  produces<reco::EvtPlaneCollection>();
  for (int i = 0; i < NumEPNames; i++) {
    rp[i] = new GenPlane(EPNames[i], EPEtaMin1[i], EPEtaMax1[i], EPEtaMin2[i], EPEtaMax2[i], EPOrder[i]);
  }
  for (int i = 0; i < NumEPNames; i++) {
    flat[i] = new HiEvtPlaneFlatten();
    flat[i]->init(FlatOrder_, NumFlatBins_, EPNames[i], EPOrder[i]);
  }
}

EvtPlaneProducer::~EvtPlaneProducer() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  for (int i = 0; i < NumEPNames; i++) {
    delete flat[i];
  }
}

//
// member functions
//

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

  if (loadDB_ && (hiWatcher.check(iSetup) || hirpWatcher.check(iSetup))) {
    //
    //Get Size of Centrality Table
    //
    edm::ESHandle<CentralityTable> centDB_;
    iSetup.get<HeavyIonRcd>().get(centralityLabel_, centDB_);
    nCentBins_ = centDB_->m_table.size();
    for (int i = 0; i < NumEPNames; i++) {
      if (caloCentRef_ > 0) {
        int minbin = (caloCentRef_ - caloCentRefWidth_ / 2.) * nCentBins_ / 100.;
        int maxbin = (caloCentRef_ + caloCentRefWidth_ / 2.) * nCentBins_ / 100.;
        minbin /= CentBinCompression_;
        maxbin /= CentBinCompression_;
        if (minbin > 0 && maxbin >= minbin) {
          if (EPDet[i] == HF || EPDet[i] == Castor)
            flat[i]->setCaloCentRefBins(minbin, maxbin);
        }
      }
    }

    //
    //Get flattening parameter file.
    //
    if (loadDB_) {
      edm::ESHandle<RPFlatParams> flatparmsDB_;
      iSetup.get<HeavyIonRPRcd>().get(flatparmsDB_);
      LoadEPDB db(flatparmsDB_, flat);
      if (!db.IsSuccess()) {
        loadDB_ = kFALSE;
      }
    }

  }  //rp record change

  //
  //Get Centrality
  //
  int bin = 0;
  if (loadDB_) {
    edm::Handle<int> cbin_;
    iEvent.getByToken(centralityBinToken, cbin_);
    int cbin = *cbin_;
    bin = cbin / CentBinCompression_;
  }
  //
  //Get Vertex
  //
  int vs_sell = 0.;
  float vzr_sell;
  iEvent.getByToken(vertexToken, vertex_);
  const reco::VertexCollection *vertices3 = nullptr;
  if (vertex_.isValid()) {
    vertices3 = vertex_.product();
    vs_sell = vertices3->size();
  }
  if (vs_sell > 0) {
    vzr_sell = vertices3->begin()->z();
  } else
    vzr_sell = -999.9;
  //
  for (int i = 0; i < NumEPNames; i++)
    rp[i]->reset();
  if (vzr_sell < minvtx_ or vzr_sell > maxvtx_)
    return;

  //calorimetry part

  double tower_eta, tower_phi;
  double tower_energyet, tower_energyet_e, tower_energyet_h;

  iEvent.getByToken(caloToken, caloCollection_);

  if (caloCollection_.isValid()) {
    for (CaloTowerCollection::const_iterator j = caloCollection_->begin(); j != caloCollection_->end(); j++) {
      tower_eta = j->eta();
      tower_phi = j->phi();
      tower_energyet_e = j->emEt();
      tower_energyet_h = j->hadEt();
      tower_energyet = tower_energyet_e + tower_energyet_h;
      double minet = minet_;
      double maxet = maxet_;
      for (int i = 0; i < NumEPNames; i++) {
        if (minet_ < 0)
          minet = minTransverse[i];
        if (maxet_ < 0)
          maxet = maxTransverse[i];
        if (tower_energyet < minet)
          continue;
        if (tower_energyet > maxet)
          continue;
        if (EPDet[i] == HF) {
          double w = tower_energyet;
          if (loadDB_)
            w = tower_energyet * flat[i]->getEtScale(vzr_sell, bin);
          if (EPOrder[i] == 1) {
            if (MomConsWeight[i][0] == 'y' && loadDB_) {
              w = flat[i]->getW(tower_energyet, vzr_sell, bin);
            }
            if (tower_eta < 0)
              w = -w;
          }
          rp[i]->addParticle(w, tower_energyet, sin(EPOrder[i] * tower_phi), cos(EPOrder[i] * tower_phi), tower_eta);
        }
      }
    }
  }

  //Castor part

  iEvent.getByToken(castorToken, castorCollection_);

  if (castorCollection_.isValid()) {
    for (std::vector<reco::CastorTower>::const_iterator j = castorCollection_->begin(); j != castorCollection_->end();
         j++) {
      tower_eta = j->eta();
      tower_phi = j->phi();
      tower_energyet = j->et();
      double minet = minet_;
      double maxet = maxet_;
      for (int i = 0; i < NumEPNames; i++) {
        if (EPDet[i] == Castor) {
          if (minet_ < 0)
            minet = minTransverse[i];
          if (maxet_ < 0)
            maxet = maxTransverse[i];
          if (tower_energyet < minet)
            continue;
          if (tower_energyet > maxet)
            continue;
          double w = tower_energyet;
          if (EPOrder[i] == 1) {
            if (MomConsWeight[i][0] == 'y' && loadDB_) {
              w = flat[i]->getW(tower_energyet, vzr_sell, bin);
            }
            if (tower_eta < 0)
              w = -w;
          }
          rp[i]->addParticle(w, tower_energyet, sin(EPOrder[i] * tower_phi), cos(EPOrder[i] * tower_phi), tower_eta);
        }
      }
    }
  }

  //Tracking part

  double track_eta;
  double track_phi;
  double track_pt;

  double vzErr2 = 0.0, vxyErr = 0.0;
  math::XYZPoint vtxPoint(0.0, 0.0, 0.0);
  if (vertex_.isValid() && !vertex_->empty()) {
    vtxPoint = vertex_->begin()->position();
    vzErr2 = (vertex_->begin()->zError()) * (vertex_->begin()->zError());
    vxyErr = vertex_->begin()->xError() * vertex_->begin()->yError();
  }

  iEvent.getByToken(trackToken, trackCollection_);
  if (trackCollection_.isValid()) {
    for (reco::TrackCollection::const_iterator j = trackCollection_->begin(); j != trackCollection_->end(); j++) {
      bool accepted = true;
      bool isPixel = false;
      // determine if the track is a pixel track
      if (j->numberOfValidHits() < 7)
        isPixel = true;

      // determine the vertex significance
      double d0 = 0.0, dz = 0.0, d0sigma = 0.0, dzsigma = 0.0;
      d0 = -1. * j->dxy(vtxPoint);
      dz = j->dz(vtxPoint);
      d0sigma = sqrt(j->d0Error() * j->d0Error() + vxyErr);
      dzsigma = sqrt(j->dzError() * j->dzError() + vzErr2);

      // cuts for pixel tracks
      if (isPixel) {
        // dz significance cut
        if (fabs(dz / dzsigma) > dzerr_)
          accepted = false;
        // chi2/ndof cut
        if (j->normalizedChi2() > chi2_)
          accepted = false;
      }
      // cuts for full tracks
      if (!isPixel) {
        // dz and d0 significance cuts
        if (fabs(dz / dzsigma) > 3)
          accepted = false;
        if (fabs(d0 / d0sigma) > 3)
          accepted = false;
        // pt resolution cut
        if (j->ptError() / j->pt() > 0.1)
          accepted = false;
        // number of valid hits cut
        if (j->numberOfValidHits() < 12)
          accepted = false;
      }
      if (accepted) {
        track_eta = j->eta();
        track_phi = j->phi();
        track_pt = j->pt();
        double minpt = minpt_;
        double maxpt = maxpt_;
        for (int i = 0; i < NumEPNames; i++) {
          if (minpt_ < 0)
            minpt = minTransverse[i];
          if (maxpt_ < 0)
            maxpt = maxTransverse[i];
          if (track_pt < minpt)
            continue;
          if (track_pt > maxpt)
            continue;
          if (EPDet[i] == Tracker) {
            double w = track_pt;
            if (w > 2.5)
              w = 2.0;  //v2 starts decreasing above ~2.5 GeV/c
            if (EPOrder[i] == 1) {
              if (MomConsWeight[i][0] == 'y' && loadDB_) {
                w = flat[i]->getW(track_pt, vzr_sell, bin);
              }
              if (track_eta < 0)
                w = -w;
            }
            rp[i]->addParticle(w, track_pt, sin(EPOrder[i] * track_phi), cos(EPOrder[i] * track_phi), track_eta);
          }
        }
      }
    }  //end for
  }

  auto evtplaneOutput = std::make_unique<EvtPlaneCollection>();

  double ang = -10;
  double sv = 0;
  double cv = 0;
  double svNoWgt = 0;
  double cvNoWgt = 0;

  double wv = 0;
  double wv2 = 0;
  double pe = 0;
  double pe2 = 0;
  uint epmult = 0;

  for (int i = 0; i < NumEPNames; i++) {
    rp[i]->getAngle(ang, sv, cv, svNoWgt, cvNoWgt, wv, wv2, pe, pe2, epmult);
    evtplaneOutput->push_back(EvtPlane(i, 0, ang, sv, cv, wv, wv2, pe, pe2, epmult));
    evtplaneOutput->back().addLevel(3, 0., svNoWgt, cvNoWgt);
  }

  iEvent.put(std::move(evtplaneOutput));
}

//define this as a plug-in
DEFINE_FWK_MODULE(EvtPlaneProducer);