Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-23 03:25:22

0001 // -*- C++ -*-
0002 //
0003 // Package:    EvtPlaneProducer
0004 // Class:      EvtPlaneProducer
0005 //
0006 /**\class EvtPlaneProducer EvtPlaneProducer.cc RecoHI/EvtPlaneProducer/src/EvtPlaneProducer.cc
0007 
0008 Description: <one line class summary>
0009 
0010 Implementation:
0011 <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Sergey Petrushanko
0015 //         Created:  Fri Jul 11 10:05:00 2008
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <iostream>
0022 #include <ctime>
0023 #include <cmath>
0024 #include <cstdlib>
0025 
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/stream/EDProducer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/EventSetup.h"
0031 
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/Framework/interface/ESWatcher.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 
0036 #include "DataFormats/Candidate/interface/Candidate.h"
0037 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
0038 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0039 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0040 
0041 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0042 #include "DataFormats/Common/interface/Ref.h"
0043 
0044 #include "DataFormats/Common/interface/Handle.h"
0045 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0046 #include "DataFormats/CastorReco/interface/CastorTower.h"
0047 
0048 #include "FWCore/ServiceRegistry/interface/Service.h"
0049 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/VertexReco/interface/Vertex.h"
0052 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0053 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h"
0054 #include "CondFormats/HIObjects/interface/RPFlatParams.h"
0055 #include "CondFormats/DataRecord/interface/HeavyIonRPRcd.h"
0056 #include "CondFormats/DataRecord/interface/HeavyIonRcd.h"
0057 #include "CondFormats/HIObjects/interface/CentralityTable.h"
0058 
0059 #include "DataFormats/Math/interface/Point3D.h"
0060 #include "DataFormats/Common/interface/RefProd.h"
0061 #include "DataFormats/Common/interface/Ref.h"
0062 #include "DataFormats/Common/interface/RefVector.h"
0063 
0064 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0065 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0066 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h"
0067 #include "RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h"
0068 #include "RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h"
0069 
0070 using namespace std;
0071 using namespace hi;
0072 
0073 //
0074 // class decleration
0075 //
0076 
0077 namespace hi {
0078   class GenPlane {
0079   public:
0080     GenPlane(string name, double etaminval1, double etamaxval1, double etaminval2, double etamaxval2, int orderval) {
0081       epname = name;
0082       etamin1 = etaminval1;
0083       etamax1 = etamaxval1;
0084       etamin2 = etaminval2;
0085       etamax2 = etamaxval2;
0086       sumsin = 0;
0087       sumcos = 0;
0088       sumsinNoWgt = 0;
0089       sumcosNoWgt = 0;
0090 
0091       mult = 0;
0092       order = (double)orderval;
0093     }
0094     ~GenPlane() { ; }
0095     void addParticle(double w, double PtOrEt, double s, double c, double eta) {
0096       if ((eta >= etamin1 && eta < etamax1) || (etamin2 != etamax2 && eta >= etamin2 && eta < etamax2)) {
0097         sumsin += w * s;
0098         sumcos += w * c;
0099         sumsinNoWgt += s;
0100         sumcosNoWgt += c;
0101 
0102         sumw += fabs(w);
0103         sumw2 += w * w;
0104         sumPtOrEt += PtOrEt;
0105         sumPtOrEt2 += PtOrEt * PtOrEt;
0106         ++mult;
0107       }
0108     }
0109 
0110     double getAngle(double &ang,
0111                     double &sv,
0112                     double &cv,
0113                     double &svNoWgt,
0114                     double &cvNoWgt,
0115                     double &w,
0116                     double &w2,
0117                     double &PtOrEt,
0118                     double &PtOrEt2,
0119                     uint &epmult) {
0120       ang = -10;
0121       sv = sumsin;
0122       cv = sumcos;
0123       svNoWgt = sumsinNoWgt;
0124       cvNoWgt = sumcosNoWgt;
0125       w = sumw;
0126       w2 = sumw2;
0127       PtOrEt = sumPtOrEt;
0128       PtOrEt2 = sumPtOrEt2;
0129       epmult = mult;
0130       double q = sv * sv + cv * cv;
0131       if (q > 0)
0132         ang = atan2(sv, cv) / order;
0133       return ang;
0134     }
0135     void reset() {
0136       sumsin = 0;
0137       sumcos = 0;
0138       sumsinNoWgt = 0;
0139       sumcosNoWgt = 0;
0140       sumw = 0;
0141       sumw2 = 0;
0142       mult = 0;
0143       sumPtOrEt = 0;
0144       sumPtOrEt2 = 0;
0145     }
0146 
0147   private:
0148     string epname;
0149     double etamin1;
0150     double etamax1;
0151 
0152     double etamin2;
0153     double etamax2;
0154     double sumsin;
0155     double sumcos;
0156     double sumsinNoWgt;
0157     double sumcosNoWgt;
0158     uint mult;
0159     double sumw;
0160     double sumw2;
0161     double sumPtOrEt;
0162     double sumPtOrEt2;
0163     double order;
0164   };
0165 }  // namespace hi
0166 
0167 class EvtPlaneProducer : public edm::stream::EDProducer<> {
0168 public:
0169   explicit EvtPlaneProducer(const edm::ParameterSet &);
0170   ~EvtPlaneProducer() override;
0171 
0172 private:
0173   GenPlane *rp[NumEPNames];
0174 
0175   void produce(edm::Event &, const edm::EventSetup &) override;
0176 
0177   // ----------member data ---------------------------
0178   EPCuts cuts_;
0179 
0180   std::string centralityVariable_;
0181   std::string centralityMC_;
0182 
0183   edm::InputTag centralityBinTag_;
0184   edm::EDGetTokenT<int> centralityBinToken_;
0185 
0186   edm::InputTag vertexTag_;
0187   edm::EDGetTokenT<std::vector<reco::Vertex>> vertexToken_;
0188   edm::Handle<std::vector<reco::Vertex>> vertex_;
0189 
0190   edm::InputTag caloTag_;
0191   edm::EDGetTokenT<CaloTowerCollection> caloToken_;
0192   edm::Handle<CaloTowerCollection> caloCollection_;
0193   edm::EDGetTokenT<reco::PFCandidateCollection> caloTokenPF_;
0194 
0195   edm::InputTag castorTag_;
0196   edm::EDGetTokenT<std::vector<reco::CastorTower>> castorToken_;
0197   edm::Handle<std::vector<reco::CastorTower>> castorCollection_;
0198 
0199   edm::InputTag trackTag_;
0200   edm::EDGetTokenT<reco::TrackCollection> trackToken_;
0201   edm::InputTag losttrackTag_;
0202   edm::Handle<reco::TrackCollection> trackCollection_;
0203   bool bStrack_packedPFCandidates_;
0204   bool bScalo_particleFlow_;
0205   edm::EDGetTokenT<pat::PackedCandidateCollection> packedToken_;
0206   edm::EDGetTokenT<pat::PackedCandidateCollection> lostToken_;
0207 
0208   edm::InputTag chi2MapTag_;
0209   edm::EDGetTokenT<edm::ValueMap<float>> chi2MapToken_;
0210   edm::InputTag chi2MapLostTag_;
0211   edm::EDGetTokenT<edm::ValueMap<float>> chi2MapLostToken_;
0212 
0213   edm::ESGetToken<CentralityTable, HeavyIonRcd> centralityToken_;
0214   edm::ESGetToken<RPFlatParams, HeavyIonRPRcd> flatparmsToken_;
0215 
0216   bool loadDB_;
0217   double minet_;
0218   double maxet_;
0219   double minpt_;
0220   double maxpt_;
0221   int flatnvtxbins_;
0222   double flatminvtx_;
0223   double flatdelvtx_;
0224   double dzdzerror_;
0225   double d0d0error_;
0226   double pterror_;
0227   double chi2perlayer_;
0228   double dzerr_;
0229   double dzdzerror_pix_;
0230   double chi2_;
0231   int nhitsValid_;
0232   int FlatOrder_;
0233   int NumFlatBins_;
0234   double nCentBins_;
0235   double caloCentRef_;
0236   double caloCentRefWidth_;
0237   int CentBinCompression_;
0238   int cutEra_;
0239   HiEvtPlaneFlatten *flat[NumEPNames];
0240   TrackStructure track_;
0241 
0242   edm::ESWatcher<HeavyIonRcd> hiWatcher_;
0243   edm::ESWatcher<HeavyIonRPRcd> hirpWatcher_;
0244 
0245   void fillHF(const TrackStructure &track, double vz, int bin) {
0246     double minet = minet_;
0247     double maxet = maxet_;
0248     for (int i = 0; i < NumEPNames; i++) {
0249       if (EPDet[i] != HF)
0250         continue;
0251       if (minet_ < 0)
0252         minet = minTransverse[i];
0253       if (maxet_ < 0)
0254         maxet = maxTransverse[i];
0255       if (track.et < minet)
0256         continue;
0257       if (track.et > maxet)
0258         continue;
0259       if (not passEta(track.eta, i))
0260         continue;
0261       double w = track.et;
0262       if (loadDB_)
0263         w = track.et * flat[i]->etScale(vz, bin);
0264       if (EPOrder[i] == 1) {
0265         if (MomConsWeight[i][0] == 'y' && loadDB_) {
0266           w = flat[i]->getW(track.et, vz, bin);
0267         }
0268       }
0269       rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta);
0270     }
0271   };
0272 
0273   void fillCastor(const TrackStructure &track, double vz, int bin) {
0274     double minet = minet_;
0275     double maxet = maxet_;
0276     for (int i = 0; i < NumEPNames; i++) {
0277       if (EPDet[i] == Castor) {
0278         if (minet_ < 0)
0279           minet = minTransverse[i];
0280         if (maxet_ < 0)
0281           maxet = maxTransverse[i];
0282         if (track.et < minet)
0283           continue;
0284         if (track.et > maxet)
0285           continue;
0286         if (not passEta(track.eta, i))
0287           continue;
0288         double w = track.et;
0289         if (EPOrder[i] == 1) {
0290           if (MomConsWeight[i][0] == 'y' && loadDB_) {
0291             w = flat[i]->getW(track.et, vz, bin);
0292           }
0293         }
0294         rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta);
0295       }
0296     }
0297   }
0298 
0299   bool passEta(float eta, int i) {
0300     if (EPEtaMin2[i] == EPEtaMax2[i]) {
0301       if (eta < EPEtaMin1[i])
0302         return false;
0303       if (eta > EPEtaMax1[i])
0304         return false;
0305     } else {
0306       if (eta < EPEtaMin1[i])
0307         return false;
0308       if (eta > EPEtaMax2[i])
0309         return false;
0310       if (eta > EPEtaMax1[i] && eta < EPEtaMin2[i])
0311         return false;
0312     }
0313     return true;
0314   }
0315 
0316   void fillTracker(const TrackStructure &track, double vz, int bin) {
0317     double minpt = minpt_;
0318     double maxpt = maxpt_;
0319     for (int i = 0; i < NumEPNames; i++) {
0320       if (EPDet[i] == Tracker) {
0321         if (minpt_ < 0)
0322           minpt = minTransverse[i];
0323         if (maxpt_ < 0)
0324           maxpt = maxTransverse[i];
0325         if (track.pt < minpt)
0326           continue;
0327         if (track.pt > maxpt)
0328           continue;
0329         if (not passEta(track.eta, i))
0330           continue;
0331         double w = track.pt;
0332         if (w > 2.5)
0333           w = 2.0;  //v2 starts decreasing above ~2.5 GeV/c
0334         if (EPOrder[i] == 1) {
0335           if (MomConsWeight[i][0] == 'y' && loadDB_) {
0336             w = flat[i]->getW(track.pt, vz, bin);
0337           }
0338         }
0339         rp[i]->addParticle(w, track.pt, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta);
0340       }
0341     }
0342   };
0343 };
0344 
0345 EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig)
0346     : centralityVariable_(iConfig.getParameter<std::string>("centralityVariable")),
0347       centralityBinTag_(iConfig.getParameter<edm::InputTag>("centralityBinTag")),
0348       vertexTag_(iConfig.getParameter<edm::InputTag>("vertexTag")),
0349       caloTag_(iConfig.getParameter<edm::InputTag>("caloTag")),
0350       castorTag_(iConfig.getParameter<edm::InputTag>("castorTag")),
0351       trackTag_(iConfig.getParameter<edm::InputTag>("trackTag")),
0352       losttrackTag_(iConfig.getParameter<edm::InputTag>("lostTag")),
0353       chi2MapTag_(iConfig.getParameter<edm::InputTag>("chi2MapTag")),
0354       chi2MapLostTag_(iConfig.getParameter<edm::InputTag>("chi2MapLostTag")),
0355       loadDB_(iConfig.getParameter<bool>("loadDB")),
0356       minet_(iConfig.getParameter<double>("minet")),
0357       maxet_(iConfig.getParameter<double>("maxet")),
0358       minpt_(iConfig.getParameter<double>("minpt")),
0359       maxpt_(iConfig.getParameter<double>("maxpt")),
0360       flatnvtxbins_(iConfig.getParameter<int>("flatnvtxbins")),
0361       flatminvtx_(iConfig.getParameter<double>("flatminvtx")),
0362       flatdelvtx_(iConfig.getParameter<double>("flatdelvtx")),
0363       dzdzerror_(iConfig.getParameter<double>("dzdzerror")),
0364       d0d0error_(iConfig.getParameter<double>("d0d0error")),
0365       pterror_(iConfig.getParameter<double>("pterror")),
0366       chi2perlayer_(iConfig.getParameter<double>("chi2perlayer")),
0367       dzdzerror_pix_(iConfig.getParameter<double>("dzdzerror_pix")),
0368       chi2_(iConfig.getParameter<double>("chi2")),
0369       nhitsValid_(iConfig.getParameter<int>("nhitsValid")),
0370       FlatOrder_(iConfig.getParameter<int>("FlatOrder")),
0371       NumFlatBins_(iConfig.getParameter<int>("NumFlatBins")),
0372       caloCentRef_(iConfig.getParameter<double>("caloCentRef")),
0373       caloCentRefWidth_(iConfig.getParameter<double>("caloCentRefWidth")),
0374       CentBinCompression_(iConfig.getParameter<int>("CentBinCompression")),
0375       cutEra_(iConfig.getParameter<int>("cutEra"))
0376 
0377 {
0378   if (cutEra_ > 3)
0379     throw edm::Exception(edm::errors::Configuration) << "wrong range in cutEra parameter";
0380   cuts_ = EPCuts(
0381       static_cast<EP_ERA>(cutEra_), pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_);
0382   nCentBins_ = 200.;
0383 
0384   if (iConfig.exists("nonDefaultGlauberModel")) {
0385     centralityMC_ = iConfig.getParameter<std::string>("nonDefaultGlauberModel");
0386   }
0387   centralityToken_ = esConsumes(edm::ESInputTag("", centralityVariable_ + centralityMC_));
0388   if (loadDB_) {
0389     flatparmsToken_ = esConsumes();
0390   }
0391 
0392   centralityBinToken_ = consumes<int>(centralityBinTag_);
0393 
0394   vertexToken_ = consumes<std::vector<reco::Vertex>>(vertexTag_);
0395 
0396   bStrack_packedPFCandidates_ = (trackTag_.label().find("packedPFCandidates") != std::string::npos);
0397   bScalo_particleFlow_ = (caloTag_.label().find("particleFlow") != std::string::npos);
0398   if (bStrack_packedPFCandidates_) {
0399     packedToken_ = consumes<pat::PackedCandidateCollection>(trackTag_);
0400     lostToken_ = consumes<pat::PackedCandidateCollection>(losttrackTag_);
0401     chi2MapToken_ = consumes<edm::ValueMap<float>>(chi2MapTag_);
0402     chi2MapLostToken_ = consumes<edm::ValueMap<float>>(chi2MapLostTag_);
0403 
0404   } else {
0405     if (bScalo_particleFlow_) {
0406       caloTokenPF_ = consumes<reco::PFCandidateCollection>(caloTag_);
0407     } else {
0408       caloToken_ = consumes<CaloTowerCollection>(caloTag_);
0409     }
0410     castorToken_ = consumes<std::vector<reco::CastorTower>>(castorTag_);
0411     trackToken_ = consumes<reco::TrackCollection>(trackTag_);
0412   }
0413 
0414   produces<reco::EvtPlaneCollection>();
0415   for (int i = 0; i < NumEPNames; i++) {
0416     rp[i] = new GenPlane(EPNames[i], EPEtaMin1[i], EPEtaMax1[i], EPEtaMin2[i], EPEtaMax2[i], EPOrder[i]);
0417   }
0418   for (int i = 0; i < NumEPNames; i++) {
0419     flat[i] = new HiEvtPlaneFlatten();
0420     flat[i]->init(FlatOrder_, NumFlatBins_, flatnvtxbins_, flatminvtx_, flatdelvtx_, EPNames[i], EPOrder[i]);
0421   }
0422 }
0423 
0424 EvtPlaneProducer::~EvtPlaneProducer() {
0425   // do anything here that needs to be done at desctruction time
0426   // (e.g. close files, deallocate resources etc.)
0427   for (int i = 0; i < NumEPNames; i++) {
0428     delete flat[i];
0429   }
0430 }
0431 
0432 //
0433 // member functions
0434 //
0435 
0436 // ------------ method called to produce the data  ------------
0437 void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0438   using namespace edm;
0439   using namespace std;
0440   using namespace reco;
0441   if (hiWatcher_.check(iSetup)) {
0442     //
0443     //Get Size of Centrality Table
0444     //
0445     auto const &centDB = iSetup.getData(centralityToken_);
0446     nCentBins_ = centDB.m_table.size();
0447     for (int i = 0; i < NumEPNames; i++) {
0448       if (caloCentRef_ > 0) {
0449         int minbin = (caloCentRef_ - caloCentRefWidth_ / 2.) * nCentBins_ / 100.;
0450         int maxbin = (caloCentRef_ + caloCentRefWidth_ / 2.) * nCentBins_ / 100.;
0451         minbin /= CentBinCompression_;
0452         maxbin /= CentBinCompression_;
0453         if (minbin > 0 && maxbin >= minbin) {
0454           if (EPDet[i] == HF || EPDet[i] == Castor)
0455             flat[i]->setCaloCentRefBins(minbin, maxbin);
0456         }
0457       }
0458     }
0459   }
0460   //
0461   //Get flattening parameter file.
0462   //
0463   if (loadDB_ && hirpWatcher_.check(iSetup)) {
0464     LoadEPDB db(iSetup.getData(flatparmsToken_), flat);
0465     if (!db.IsSuccess()) {
0466       loadDB_ = kFALSE;
0467     }
0468   }
0469   //
0470   //Get Centrality
0471   //
0472   int bin = 0;
0473   int cbin = 0;
0474   if (loadDB_) {
0475     cbin = iEvent.get(centralityBinToken_);
0476     bin = cbin / CentBinCompression_;
0477   }
0478   //
0479   //Get Vertex
0480   //
0481   //best vertex
0482   const reco::Vertex &vtx = iEvent.get(vertexToken_)[0];
0483   double bestvz = vtx.z();
0484   double bestvx = vtx.x();
0485   double bestvy = vtx.y();
0486   double bestvzError = vtx.zError();
0487   math::XYZPoint bestvtx(bestvx, bestvy, bestvz);
0488   math::Error<3>::type vtx_cov = vtx.covariance();
0489 
0490   for (int i = 0; i < NumEPNames; i++)
0491     rp[i]->reset();
0492   edm::Handle<edm::ValueMap<float>> chi2Map;
0493   edm::Handle<pat::PackedCandidateCollection> cands;
0494   edm::Handle<reco::PFCandidateCollection> calocands;
0495   if (bStrack_packedPFCandidates_) {
0496     for (int idx = 1; idx < 3; idx++) {
0497       if (idx == 1) {
0498         iEvent.getByToken(packedToken_, cands);
0499         iEvent.getByToken(chi2MapToken_, chi2Map);
0500       }
0501       if (idx == 2) {
0502         iEvent.getByToken(lostToken_, cands);
0503         iEvent.getByToken(chi2MapLostToken_, chi2Map);
0504       }
0505       for (unsigned int i = 0, n = cands->size(); i < n; ++i) {
0506         track_ = {};
0507         track_.centbin = cbin;
0508         const pat::PackedCandidate &pf = (*cands)[i];
0509         track_.et = pf.et();
0510         track_.eta = pf.eta();
0511         track_.phi = pf.phi();
0512         track_.pdgid = pf.pdgId();
0513         if ((idx == 1) and cuts_.isGoodHF(track_)) {
0514           fillHF(track_, bestvz, bin);
0515         }
0516         if (!pf.hasTrackDetails())
0517           continue;
0518         const reco::Track &trk = pf.pseudoTrack();
0519         track_.highPurity = pf.trackHighPurity();
0520         track_.charge = trk.charge();
0521         if (!track_.highPurity || track_.charge == 0)
0522           continue;
0523         track_.collection = idx;
0524         track_.eta = trk.eta();
0525         track_.phi = trk.phi();
0526         track_.pt = trk.pt();
0527         track_.ptError = trk.ptError();
0528         track_.numberOfValidHits = trk.numberOfValidHits();
0529         track_.algos = trk.algo();
0530         track_.dz = std::abs(trk.dz(bestvtx));
0531         track_.dxy = std::abs(trk.dxy(bestvtx));
0532         track_.dzError = std::hypot(trk.dzError(), bestvzError);
0533         track_.dxyError = trk.dxyError(bestvtx, vtx_cov);
0534         track_.dzSig = track_.dz / track_.dzError;
0535         track_.dxySig = track_.dxy / track_.dxyError;
0536         const reco::HitPattern &hit_pattern = trk.hitPattern();
0537         track_.normalizedChi2 = (*chi2Map)[pat::PackedCandidateRef(cands, i)];
0538         track_.chi2layer = (*chi2Map)[pat::PackedCandidateRef(cands, i)] / hit_pattern.trackerLayersWithMeasurement();
0539         if (cuts_.isGoodTrack(track_)) {
0540           fillTracker(track_, bestvz, bin);
0541         }
0542       }
0543     }
0544   } else {
0545     //calorimetry part
0546     if (bScalo_particleFlow_) {
0547       iEvent.getByToken(caloTokenPF_, calocands);
0548       for (unsigned int i = 0, n = calocands->size(); i < n; ++i) {
0549         track_ = {};
0550         track_.centbin = cbin;
0551         const reco::PFCandidate &pf = (*calocands)[i];
0552         track_.et = pf.et();
0553         track_.eta = pf.eta();
0554         track_.phi = pf.phi();
0555         track_.pdgid = pf.pdgId();
0556         if (cuts_.isGoodHF(track_)) {
0557           fillHF(track_, bestvz, bin);
0558         }
0559       }
0560     } else {
0561       iEvent.getByToken(caloToken_, caloCollection_);
0562       for (const auto &tower : *caloCollection_) {
0563         track_.eta = tower.eta();
0564         track_.phi = tower.phi();
0565         track_.et = tower.emEt() + tower.hadEt();
0566         track_.pdgid = 1;
0567         if (cuts_.isGoodHF(track_))
0568           fillHF(track_, bestvz, bin);
0569       }
0570     }
0571 
0572     //Castor part
0573     iEvent.getByToken(castorToken_, castorCollection_);
0574     for (const auto &tower : *castorCollection_) {
0575       track_.eta = tower.eta();
0576       track_.phi = tower.phi();
0577       track_.et = tower.et();
0578       track_.pdgid = 1;
0579       if (cuts_.isGoodCastor(track_))
0580         fillCastor(track_, bestvz, bin);
0581     }
0582     //Tracking part
0583     iEvent.getByToken(trackToken_, trackCollection_);
0584     for (const auto &trk : *trackCollection_) {
0585       track_.highPurity = trk.quality(reco::TrackBase::highPurity);
0586       track_.charge = trk.charge();
0587       if (!track_.highPurity || track_.charge == 0)
0588         continue;
0589       track_.centbin = cbin;
0590       track_.collection = 0;
0591       track_.eta = trk.eta();
0592       track_.phi = trk.phi();
0593       track_.pt = trk.pt();
0594       track_.ptError = trk.ptError();
0595       track_.numberOfValidHits = trk.numberOfValidHits();
0596       track_.algos = trk.algo();
0597       track_.dz = std::abs(trk.dz(bestvtx));
0598       track_.dxy = std::abs(trk.dxy(bestvtx));
0599       track_.dzError = std::hypot(trk.dzError(), bestvzError);
0600       track_.dxyError = trk.dxyError(bestvtx, vtx_cov);
0601       track_.dzSig = track_.dz / track_.dzError;
0602       track_.dxySig = track_.dxy / track_.dxyError;
0603       track_.normalizedChi2 = trk.normalizedChi2();
0604       track_.chi2layer = track_.normalizedChi2 / trk.hitPattern().trackerLayersWithMeasurement();
0605       if (cuts_.isGoodTrack(track_))
0606         fillTracker(track_, bestvz, bin);
0607     }
0608   }
0609 
0610   auto evtplaneOutput = std::make_unique<EvtPlaneCollection>();
0611 
0612   double ang = -10;
0613   double sv = 0;
0614   double cv = 0;
0615   double svNoWgt = 0;
0616   double cvNoWgt = 0;
0617 
0618   double wv = 0;
0619   double wv2 = 0;
0620   double pe = 0;
0621   double pe2 = 0;
0622   uint epmult = 0;
0623 
0624   for (int i = 0; i < NumEPNames; i++) {
0625     rp[i]->getAngle(ang, sv, cv, svNoWgt, cvNoWgt, wv, wv2, pe, pe2, epmult);
0626     evtplaneOutput->push_back(EvtPlane(i, 0, ang, sv, cv, wv, wv2, pe, pe2, epmult));
0627     evtplaneOutput->back().addLevel(3, 0., svNoWgt, cvNoWgt);
0628   }
0629 
0630   iEvent.put(std::move(evtplaneOutput));
0631 }
0632 
0633 //define this as a plug-in
0634 DEFINE_FWK_MODULE(EvtPlaneProducer);