HiEvtPlaneFlatProducer.cc

Go to the documentation of this file.

#include <memory>
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "Math/Vector3D.h"

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"

#include "DataFormats/HeavyIonEvent/interface/Centrality.h"

#include "DataFormats/HeavyIonEvent/interface/EvtPlane.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 "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Track/interface/CoreSimTrack.h"
#include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
#include "SimDataFormats/Vertex/interface/SimVertex.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
#include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "CondFormats/DataRecord/interface/HeavyIonRPRcd.h"
#include "CondFormats/DataRecord/interface/HeavyIonRcd.h"
#include "CondFormats/HIObjects/interface/CentralityTable.h"
#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
#include "CondFormats/HIObjects/interface/RPFlatParams.h"

#include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h"
#include "TList.h"
#include "TString.h"
#include <ctime>
#include <cstdlib>

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

using namespace std;
using namespace hi;

#include <vector>
using std::vector;

//
// class declaration
//

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

private:
  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 centralityTag_;
  edm::EDGetTokenT<reco::Centrality> centralityToken;

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

  edm::InputTag inputPlanesTag_;
  edm::EDGetTokenT<reco::EvtPlaneCollection> inputPlanesToken;

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

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

  const int FlatOrder_;
  int NumFlatBins_;
  double caloCentRef_;
  double caloCentRefWidth_;
  int CentBinCompression_;
  int Noffmin_;
  int Noffmax_;
  HiEvtPlaneFlatten* flat[NumEPNames];
  bool useOffsetPsi_;
  double nCentBins_;
};
//
// constants, enums and typedefs
//

typedef std::vector<TrackingParticle> TrackingParticleCollection;
typedef TrackingParticleRefVector::iterator tp_iterator;

//
// static data member definitions
//

//
// constructors and destructor
//
HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig)
    : centralityVariable_(iConfig.getParameter<std::string>("centralityVariable")),
      centralityBinTag_(iConfig.getParameter<edm::InputTag>("centralityBinTag")),
      centralityTag_(iConfig.getParameter<edm::InputTag>("centralityTag")),
      vertexTag_(iConfig.getParameter<edm::InputTag>("vertexTag")),
      inputPlanesTag_(iConfig.getParameter<edm::InputTag>("inputPlanesTag")),
      trackTag_(iConfig.getParameter<edm::InputTag>("trackTag")),
      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")),
      Noffmin_(iConfig.getParameter<int>("Noffmin")),
      Noffmax_(iConfig.getParameter<int>("Noffmax")),
      useOffsetPsi_(iConfig.getParameter<bool>("useOffsetPsi")) {
  //  UseEtHF = kFALSE;
  nCentBins_ = 200.;

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

  centralityBinToken = consumes<int>(centralityBinTag_);

  centralityToken = consumes<reco::Centrality>(centralityTag_);

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

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

  inputPlanesToken = consumes<reco::EvtPlaneCollection>(inputPlanesTag_);

  //register your products
  produces<reco::EvtPlaneCollection>();
  //now do what ever other initialization is needed
  for (int i = 0; i < NumEPNames; i++) {
    flat[i] = new HiEvtPlaneFlatten();
    flat[i]->init(FlatOrder_, NumFlatBins_, EPNames[i], EPOrder[i]);
  }
}

HiEvtPlaneFlatProducer::~HiEvtPlaneFlatProducer() {
  // 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 HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  using namespace edm;
  using namespace std;
  using namespace reco;

  //
  //Get Flattening Parameters
  //
  if (hiWatcher.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++) {
      flat[i]->setCaloCentRefBins(-1, -1);
      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);
        }
      }
    }
  }

  if (hirpWatcher.check(iSetup)) {
    edm::ESHandle<RPFlatParams> flatparmsDB_;
    iSetup.get<HeavyIonRPRcd>().get(flatparmsDB_);
    LoadEPDB db(flatparmsDB_, flat);
    if (!db.IsSuccess())
      return;
  }
  //
  //Get Centrality
  //

  int bin = 0;
  edm::Handle<int> cbin_;
  iEvent.getByToken(centralityBinToken, cbin_);
  int cbin = *cbin_;
  bin = cbin / CentBinCompression_;

  if (Noffmin_ >= 0) {
    edm::Handle<reco::Centrality> centrality_;
    iEvent.getByToken(centralityToken, centrality_);
    int Noff = centrality_->Ntracks();
    if ((Noff < Noffmin_) or (Noff >= Noffmax_)) {
      return;
    }
  }
  //
  //Get Vertex
  //
  int vs_sell;  // vertex collection size
  float vzr_sell;
  edm::Handle<std::vector<reco::Vertex>> vertex_;
  iEvent.getByToken(vertexToken, vertex_);
  const reco::VertexCollection* vertices3 = vertex_.product();
  vs_sell = vertices3->size();
  if (vs_sell > 0) {
    vzr_sell = vertices3->begin()->z();
  } else
    vzr_sell = -999.9;

  //
  //Get Event Planes
  //

  edm::Handle<reco::EvtPlaneCollection> evtPlanes_;
  iEvent.getByToken(inputPlanesToken, evtPlanes_);

  if (!evtPlanes_.isValid()) {
    //    cout << "Error! Can't get hiEvtPlane product!" << endl;
    return;
  }

  auto evtplaneOutput = std::make_unique<EvtPlaneCollection>();
  EvtPlane* ep[NumEPNames];
  for (int i = 0; i < NumEPNames; i++) {
    ep[i] = nullptr;
  }
  int indx = 0;
  for (EvtPlaneCollection::const_iterator rp = evtPlanes_->begin(); rp != evtPlanes_->end(); rp++) {
    double psiOffset = rp->angle(0);
    double s = rp->sumSin(0);
    double c = rp->sumCos(0);
    uint m = rp->mult();

    double soff = s;
    double coff = c;
    if (useOffsetPsi_) {
      soff = flat[indx]->getSoffset(s, vzr_sell, bin);
      coff = flat[indx]->getCoffset(c, vzr_sell, bin);
      psiOffset = flat[indx]->getOffsetPsi(soff, coff);
    }
    double psiFlat = flat[indx]->getFlatPsi(psiOffset, vzr_sell, bin);
    ep[indx] =
        new EvtPlane(indx, 2, psiFlat, soff, coff, rp->sumw(), rp->sumw2(), rp->sumPtOrEt(), rp->sumPtOrEt2(), m);
    ep[indx]->addLevel(0, rp->angle(0), s, c);
    ep[indx]->addLevel(3, 0., rp->sumSin(3), rp->sumCos(3));
    if (useOffsetPsi_)
      ep[indx]->addLevel(1, psiOffset, soff, coff);
    ++indx;
  }

  for (int i = 0; i < NumEPNames; i++) {
    if (ep[i] != nullptr)
      evtplaneOutput->push_back(*ep[i]);
  }
  iEvent.put(std::move(evtplaneOutput));
  for (int i = 0; i < indx; i++)
    delete ep[i];
}

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