HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc

0001 #include <memory>
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/stream/EDProducer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "Math/Vector3D.h"
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/Framework/interface/ESWatcher.h"
0013 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0014
0015 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
0016
0017 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/VertexReco/interface/Vertex.h"
0021 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0022
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0025 #include "CondFormats/DataRecord/interface/HeavyIonRPRcd.h"
0026 #include "CondFormats/DataRecord/interface/HeavyIonRcd.h"
0027 #include "CondFormats/HIObjects/interface/CentralityTable.h"
0028 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0029 #include "CondFormats/HIObjects/interface/RPFlatParams.h"
0030
0031 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h"
0032 #include <ctime>
0033 #include <cstdlib>
0034
0035 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h"
0036 #include "RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h"
0037
0038 using namespace std;
0039 using namespace hi;
0040
0041 #include <vector>
0042 using std::vector;
0043
0044 //
0045 // class declaration
0046 //
0047
0048 class HiEvtPlaneFlatProducer : public edm::stream::EDProducer<> {
0049 public:
0050   explicit HiEvtPlaneFlatProducer(const edm::ParameterSet&);
0051   ~HiEvtPlaneFlatProducer() override;
0052
0053 private:
0054   void produce(edm::Event&, const edm::EventSetup&) override;
0055
0056   // ----------member data ---------------------------
0057
0058   std::string centralityVariable_;
0059   std::string centralityMC_;
0060
0061   edm::InputTag centralityBinTag_;
0062   edm::EDGetTokenT<int> centralityBinToken_;
0063
0064   edm::InputTag centralityTag_;
0065   edm::EDGetTokenT<reco::Centrality> centralityToken_;
0066
0067   edm::InputTag vertexTag_;
0068   edm::EDGetTokenT<std::vector<reco::Vertex>> vertexToken_;
0069
0070   edm::InputTag inputPlanesTag_;
0071   edm::EDGetTokenT<reco::EvtPlaneCollection> inputPlanesToken_;
0072
0073   edm::InputTag trackTag_;
0074   edm::EDGetTokenT<reco::TrackCollection> trackToken_;
0075   edm::Handle<reco::TrackCollection> trackCollection_;
0076
0077   edm::ESWatcher<HeavyIonRcd> hiWatcher;
0078   edm::ESWatcher<HeavyIonRPRcd> hirpWatcher;
0079   edm::ESGetToken<CentralityTable, HeavyIonRcd> centralityESToken_;
0080   edm::ESGetToken<RPFlatParams, HeavyIonRPRcd> flatparmsToken_;
0081
0082   const int FlatOrder_;
0083   int NumFlatBins_;
0084   int flatnvtxbins_;
0085   double flatminvtx_;
0086   double flatdelvtx_;
0087   double caloCentRef_;
0088   double caloCentRefWidth_;
0089   int CentBinCompression_;
0090   HiEvtPlaneFlatten* flat[NumEPNames];
0091   bool useOffsetPsi_;
0092   double nCentBins_;
0093 };
0094 //
0095 // constructors and destructor
0096 //
0097 HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig)
0098     : centralityVariable_(iConfig.getParameter<std::string>("centralityVariable")),
0099       centralityBinTag_(iConfig.getParameter<edm::InputTag>("centralityBinTag")),
0100       centralityTag_(iConfig.getParameter<edm::InputTag>("centralityTag")),
0101       vertexTag_(iConfig.getParameter<edm::InputTag>("vertexTag")),
0102       inputPlanesTag_(iConfig.getParameter<edm::InputTag>("inputPlanesTag")),
0103       trackTag_(iConfig.getParameter<edm::InputTag>("trackTag")),
0104       FlatOrder_(iConfig.getParameter<int>("FlatOrder")),
0105       NumFlatBins_(iConfig.getParameter<int>("NumFlatBins")),
0106       flatnvtxbins_(iConfig.getParameter<int>("flatnvtxbins")),
0107       flatminvtx_(iConfig.getParameter<double>("flatminvtx")),
0108       flatdelvtx_(iConfig.getParameter<double>("flatdelvtx")),
0109       caloCentRef_(iConfig.getParameter<double>("caloCentRef")),
0110       caloCentRefWidth_(iConfig.getParameter<double>("caloCentRefWidth")),
0111       CentBinCompression_(iConfig.getParameter<int>("CentBinCompression")),
0112       useOffsetPsi_(iConfig.getParameter<bool>("useOffsetPsi")) {
0113   nCentBins_ = 200.;
0114
0115   if (iConfig.exists("nonDefaultGlauberModel")) {
0116     centralityMC_ = iConfig.getParameter<std::string>("nonDefaultGlauberModel");
0117   }
0118   centralityESToken_ = esConsumes(edm::ESInputTag("", centralityVariable_ + centralityMC_));
0119   flatparmsToken_ = esConsumes();
0120
0121   centralityBinToken_ = consumes<int>(centralityBinTag_);
0122
0123   centralityToken_ = consumes<reco::Centrality>(centralityTag_);
0124
0125   vertexToken_ = consumes<std::vector<reco::Vertex>>(vertexTag_);
0126
0127   trackToken_ = consumes<reco::TrackCollection>(trackTag_);
0128
0129   inputPlanesToken_ = consumes<reco::EvtPlaneCollection>(inputPlanesTag_);
0130
0131   //register your products
0132   produces<reco::EvtPlaneCollection>();
0133   //now do what ever other initialization is needed
0134   for (int i = 0; i < NumEPNames; i++) {
0135     flat[i] = new HiEvtPlaneFlatten();
0136     flat[i]->init(FlatOrder_, NumFlatBins_, flatnvtxbins_, flatminvtx_, flatdelvtx_, EPNames[i], EPOrder[i]);
0137   }
0138 }
0139
0140 HiEvtPlaneFlatProducer::~HiEvtPlaneFlatProducer() {
0141   // do anything here that needs to be done at desctruction time
0142   // (e.g. close files, deallocate resources etc.)
0143   for (int i = 0; i < NumEPNames; i++) {
0144     delete flat[i];
0145   }
0146 }
0147
0148 //
0149 // member functions
0150 //
0151
0152 // ------------ method called to produce the data  ------------
0153 void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0154   using namespace edm;
0155   using namespace std;
0156   using namespace reco;
0157
0158   if (hiWatcher.check(iSetup)) {
0159     //
0160     //Get Size of Centrality Table
0161     //
0162     auto const& centDB = iSetup.getData(centralityESToken_);
0163     nCentBins_ = centDB.m_table.size();
0164     for (int i = 0; i < NumEPNames; i++) {
0165       if (caloCentRef_ > 0) {
0166         int minbin = (caloCentRef_ - caloCentRefWidth_ / 2.) * nCentBins_ / 100.;
0167         int maxbin = (caloCentRef_ + caloCentRefWidth_ / 2.) * nCentBins_ / 100.;
0168         minbin /= CentBinCompression_;
0169         maxbin /= CentBinCompression_;
0170         if (minbin > 0 && maxbin >= minbin) {
0171           if (EPDet[i] == HF || EPDet[i] == Castor)
0172             flat[i]->setCaloCentRefBins(minbin, maxbin);
0173         }
0174       }
0175     }
0176   }
0177   //
0178   //Get flattening parameter file.
0179   //
0180   if (hirpWatcher.check(iSetup)) {
0181     LoadEPDB db(iSetup.getData(flatparmsToken_), flat);
0182   }  //rp record change
0183
0184   //
0185   //Get Centrality
0186   //
0187   int bin = 0;
0188   int cbin = 0;
0189   cbin = iEvent.get(centralityBinToken_);
0190   bin = cbin / CentBinCompression_;
0191   //
0192   //Get Vertex
0193   //
0194
0195   //best vertex
0196   double bestvz = -999.9;
0197   const reco::Vertex& vtx = iEvent.get(vertexToken_)[0];
0198   bestvz = vtx.z();
0199
0200   //
0201   //Get Event Planes
0202   //
0203
0204   auto const& evtPlanes = iEvent.get(inputPlanesToken_);
0205
0206   auto evtplaneOutput = std::make_unique<EvtPlaneCollection>();
0207   EvtPlane* ep[NumEPNames];
0208   for (int i = 0; i < NumEPNames; i++) {
0209     ep[i] = nullptr;
0210   }
0211   int indx = 0;
0212   for (auto&& rp : (evtPlanes)) {
0213     double s = rp.sumSin(0);
0214     double c = rp.sumCos(0);
0215     uint m = rp.mult();
0216     double soff = s;
0217     double coff = c;
0218     double psiOffset = -10;
0219     double psiFlat = -10;
0220     if (rp.angle(0) > -5) {
0221       if (useOffsetPsi_) {
0222         soff = flat[indx]->soffset(s, bestvz, bin);
0223         coff = flat[indx]->coffset(c, bestvz, bin);
0224         psiOffset = flat[indx]->offsetPsi(soff, coff);
0225       }
0226       psiFlat = flat[indx]->getFlatPsi(psiOffset, bestvz, bin);
0227     }
0228     ep[indx] = new EvtPlane(indx, 2, psiFlat, soff, coff, rp.sumw(), rp.sumw2(), rp.sumPtOrEt(), rp.sumPtOrEt2(), m);
0229     ep[indx]->addLevel(0, rp.angle(0), s, c);
0230     ep[indx]->addLevel(3, 0., rp.sumSin(3), rp.sumCos(3));
0231     if (useOffsetPsi_)
0232       ep[indx]->addLevel(1, psiOffset, soff, coff);
0233     ++indx;
0234   }
0235
0236   for (int i = 0; i < NumEPNames; i++) {
0237     if (ep[i] != nullptr)
0238       evtplaneOutput->push_back(*ep[i]);
0239   }
0240   iEvent.put(std::move(evtplaneOutput));
0241   for (int i = 0; i < indx; i++)
0242     delete ep[i];
0243 }
0244
0245 //define this as a plug-in
0246 DEFINE_FWK_MODULE(HiEvtPlaneFlatProducer);