PixelVertexFinding/plugins/FastPrimaryVertexWithWeightsProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FastPrimaryVertexWithWeightsProducer
0004 // Class:      FastPrimaryVertexWithWeightsProducer
0005 //
0006 /**\class FastPrimaryVertexWithWeightsProducer FastPrimaryVertexWithWeightsProducer.cc RecoBTag/FastPrimaryVertexWithWeightsProducer/src/FastPrimaryVertexWithWeightsProducer.cc
0007
0008  Description: The FastPrimaryVertex is an algorithm used to find the primary vertex of an event @HLT. It takes the pixel clusters compabible with a jet and project it to the beamSpot with the eta-angle of the jet. The peak on the z-projected clusters distribution is our FastPrimaryVertex.
0009
0010
0011  Implementation:
0012      [Notes on implementation]
0013 */
0014 //
0015 // Original Author:  Silvio DONATO
0016 //         Created:  Wed Dec 18 10:05:40 CET 2013
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022 #include <vector>
0023 #include <cmath>
0024 // user include files
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0033
0034 #include "DataFormats/JetReco/interface/CaloJet.h"
0035 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0036 #include "DataFormats/Math/interface/deltaPhi.h"
0037 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0038 #include "DataFormats/VertexReco/interface/Vertex.h"
0039 #include "DataFormats/TrackReco/interface/Track.h"
0040 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0041 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0042 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0043 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0044
0045 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0046 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0047 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0048 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0049
0050 #include "RecoVertex/PixelVertexFinding/interface/FindPeakFastPV.h"
0051
0052 #include "FWCore/Utilities/interface/InputTag.h"
0053
0054 using namespace std;
0055
0056 class FastPrimaryVertexWithWeightsProducer : public edm::stream::EDProducer<> {
0057 public:
0058   explicit FastPrimaryVertexWithWeightsProducer(const edm::ParameterSet&);
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062   void produce(edm::Event&, const edm::EventSetup&) override;
0063
0064   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> const m_geomToken;
0065   edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> const m_pixelCPEToken;
0066
0067   const double m_maxZ;             // Use only pixel clusters with |z| < maxZ
0068   const edm::InputTag m_clusters;  // PixelClusters InputTag
0069   edm::EDGetTokenT<SiPixelClusterCollectionNew> clustersToken;
0070   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken;
0071   edm::EDGetTokenT<edm::View<reco::Jet> > jetsToken;
0072
0073   // PARAMETERS USED IN THE BARREL PIXEL CLUSTERS PROJECTION
0074   const int m_njets;                  // Use only the first njets
0075   const double m_maxJetEta;           // Use only jets with |eta| < maxJetEta
0076   const double m_minJetPt;            // Use only jets with Pt > minJetPt
0077   const bool m_barrel;                // Use clusters from pixel endcap
0078   const double m_maxSizeX;            // Use only pixel clusters with sizeX <= maxSizeX
0079   const double m_maxDeltaPhi;         // Use only pixel clusters with DeltaPhi(Jet,Cluster) < maxDeltaPhi
0080   const double m_weight_charge_down;  // Use only pixel clusters with ClusterCharge > weight_charge_down
0081   const double m_weight_charge_up;    // Use only pixel clusters with ClusterCharge < weight_charge_up
0082       //It is the ratio between pixel cell height and width along z coordinate about 285µm/150µm=1.9
0083   const double m_PixelCellHeightOverWidth;
0084   // Use only pixel clusters with sizeY > PixelCellHeightOverWidth * |jetZOverRho| + minSizeY_q
0085   const double m_minSizeY_q;
0086   // Use only pixel clusters with sizeY < PixelCellHeightOverWidth * |jetZOverRho| + maxSizeY_q
0087   const double m_maxSizeY_q;
0088
0089   // PARAMETERS USED TO WEIGHT THE BARREL PIXEL CLUSTERS
0090   // The cluster weight is defined as weight = weight_dPhi * weight_sizeY  * weight_rho * weight_sizeX1 * weight_charge
0091
0092   const double m_weight_dPhi;         // used in weight_dPhi = exp(-|DeltaPhi(JetCluster)|/m_weight_dPhi)
0093   const double m_weight_SizeX1;       // used in weight_SizeX1 = (ClusterSizeX==2)*1+(ClusterSizeX==1)*m_weight_SizeX1;
0094   const double m_weight_rho_up;       // used in weight_rho = (m_weight_rho_up - ClusterRho)/m_weight_rho_up
0095   const double m_weight_charge_peak;  // Give the maximum weight_charge for a cluster with Charge = m_weight_charge_peak
0096       // Give the maximum weight_sizeY for a cluster with sizeY = PixelCellHeightOverWidth * |jetZOverRho| + peakSizeY_q
0097   const double m_peakSizeY_q;
0098
0099   // PARAMETERS USED IN THE ENDCAP PIXEL CLUSTERS PROJECTION
0100   const bool m_endCap;            // Use clusters from pixel endcap
0101   const double m_minJetEta_EC;    // Use only jets with |eta| > minJetEta_EC
0102   const double m_maxJetEta_EC;    // Use only jets with |eta| < maxJetEta_EC
0103   const double m_maxDeltaPhi_EC;  // Use only pixel clusters with DeltaPhi(Jet,Cluster) < maxDeltaPhi_EC
0104
0105   // PARAMETERS USED TO WEIGHT THE ENDCAP PIXEL CLUSTERS
0106   const double m_EC_weight;       // In EndCap the weight is defined as weight = m_EC_weight*(weight_dPhi)
0107   const double m_weight_dPhi_EC;  // Used in weight_dPhi = exp(-|DeltaPhi|/m_weight_dPhi_EC )
0108
0109   // PARAMETERS USED TO FIND THE FASTPV AS PEAK IN THE Z-PROJECTIONS DISTRIBUTION
0110   // First Iteration: look for a cluster with a width = m_zClusterWidth_step1
0111   const double m_zClusterWidth_step1;  // cluster width in step1
0112
0113   // Second Iteration: use only z-projections with weight > weightCut_step2 and look for a cluster with a width = m_zClusterWidth_step2, within of weightCut_step2 of the previous result
0114   const double m_zClusterWidth_step2;       // cluster width in step2
0115   const double m_zClusterSearchArea_step2;  // cluster width in step2
0116   const double m_weightCut_step2;           // minimum z-projections weight required in step2
0117
0118   // Third Iteration: use only z-projections with weight > weightCut_step3 and look for a cluster with a width = m_zClusterWidth_step3, within of weightCut_step3 of the previous result
0119   const double m_zClusterWidth_step3;       // cluster width in step3
0120   const double m_zClusterSearchArea_step3;  // cluster width in step3
0121   const double m_weightCut_step3;           // minimum z-projections weight required in step3
0122
0123   // use the jetPt weighting
0124   const bool m_ptWeighting;           // use weight=weight*pt_weigth ?;
0125   const double m_ptWeighting_slope;   // pt_weigth= pt*ptWeighting_slope +m_ptWeighting_offset;
0126   const double m_ptWeighting_offset;  //
0127 };
0128
0129 FastPrimaryVertexWithWeightsProducer::FastPrimaryVertexWithWeightsProducer(const edm::ParameterSet& iConfig)
0130     : m_geomToken(esConsumes()),
0131       m_pixelCPEToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("pixelCPE")))),
0132       m_maxZ(iConfig.getParameter<double>("maxZ")),
0133       m_njets(iConfig.getParameter<int>("njets")),
0134       m_maxJetEta(iConfig.getParameter<double>("maxJetEta")),
0135       m_minJetPt(iConfig.getParameter<double>("minJetPt")),
0136
0137       m_barrel(iConfig.getParameter<bool>("barrel")),
0138       m_maxSizeX(iConfig.getParameter<double>("maxSizeX")),
0139       m_maxDeltaPhi(iConfig.getParameter<double>("maxDeltaPhi")),
0140       m_weight_charge_down(iConfig.getParameter<double>("weight_charge_down")),
0141       m_weight_charge_up(iConfig.getParameter<double>("weight_charge_up")),
0142       m_PixelCellHeightOverWidth(iConfig.getParameter<double>("PixelCellHeightOverWidth")),
0143       m_minSizeY_q(iConfig.getParameter<double>("minSizeY_q")),
0144       m_maxSizeY_q(iConfig.getParameter<double>("maxSizeY_q")),
0145
0146       m_weight_dPhi(iConfig.getParameter<double>("weight_dPhi")),
0147       m_weight_SizeX1(iConfig.getParameter<double>("weight_SizeX1")),
0148       m_weight_rho_up(iConfig.getParameter<double>("weight_rho_up")),
0149       m_weight_charge_peak(iConfig.getParameter<double>("weight_charge_peak")),
0150       m_peakSizeY_q(iConfig.getParameter<double>("peakSizeY_q")),
0151
0152       m_endCap(iConfig.getParameter<bool>("endCap")),
0153       m_minJetEta_EC(iConfig.getParameter<double>("minJetEta_EC")),
0154       m_maxJetEta_EC(iConfig.getParameter<double>("maxJetEta_EC")),
0155       m_maxDeltaPhi_EC(iConfig.getParameter<double>("maxDeltaPhi_EC")),
0156       m_EC_weight(iConfig.getParameter<double>("EC_weight")),
0157       m_weight_dPhi_EC(iConfig.getParameter<double>("weight_dPhi_EC")),
0158
0159       m_zClusterWidth_step1(iConfig.getParameter<double>("zClusterWidth_step1")),
0160
0161       m_zClusterWidth_step2(iConfig.getParameter<double>("zClusterWidth_step2")),
0162       m_zClusterSearchArea_step2(iConfig.getParameter<double>("zClusterSearchArea_step2")),
0163       m_weightCut_step2(iConfig.getParameter<double>("weightCut_step2")),
0164
0165       m_zClusterWidth_step3(iConfig.getParameter<double>("zClusterWidth_step3")),
0166       m_zClusterSearchArea_step3(iConfig.getParameter<double>("zClusterSearchArea_step3")),
0167       m_weightCut_step3(iConfig.getParameter<double>("weightCut_step3")),
0168
0169       m_ptWeighting(iConfig.getParameter<bool>("ptWeighting")),
0170       m_ptWeighting_slope(iConfig.getParameter<double>("ptWeighting_slope")),
0171       m_ptWeighting_offset(iConfig.getParameter<double>("ptWeighting_offset")) {
0172   clustersToken = consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("clusters"));
0173   beamSpotToken = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"));
0174   jetsToken = consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("jets"));
0175
0176   produces<reco::VertexCollection>();
0177   produces<float>();
0178 }
0179
0180 void FastPrimaryVertexWithWeightsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0181   edm::ParameterSetDescription desc;
0182   desc.add<edm::InputTag>("clusters", edm::InputTag("hltSiPixelClusters"));
0183   desc.add<edm::InputTag>("beamSpot", edm::InputTag("hltOnlineBeamSpot"));
0184   desc.add<edm::InputTag>("jets", edm::InputTag("hltCaloJetL1FastJetCorrected"));
0185   desc.add<std::string>("pixelCPE", "hltESPPixelCPEGeneric");
0186   desc.add<double>("maxZ", 19.0);
0187   desc.add<int>("njets", 999);
0188   desc.add<double>("maxJetEta", 2.6);
0189   desc.add<double>("minJetPt", 40.);
0190   desc.add<bool>("barrel", true);
0191   desc.add<double>("maxSizeX", 2.1);
0192   desc.add<double>("maxDeltaPhi", 0.21);
0193   desc.add<double>("PixelCellHeightOverWidth", 1.8);
0194   desc.add<double>("weight_charge_down", 11. * 1000.);
0195   desc.add<double>("weight_charge_up", 190. * 1000.);
0196   desc.add<double>("maxSizeY_q", 2.0);
0197   desc.add<double>("minSizeY_q", -0.6);
0198   desc.add<double>("weight_dPhi", 0.13888888);
0199   desc.add<double>("weight_SizeX1", 0.88);
0200   desc.add<double>("weight_rho_up", 22.);
0201   desc.add<double>("weight_charge_peak", 22. * 1000.);
0202   desc.add<double>("peakSizeY_q", 1.0);
0203   desc.add<bool>("endCap", true);
0204   desc.add<double>("minJetEta_EC", 1.3);
0205   desc.add<double>("maxJetEta_EC", 2.6);
0206   desc.add<double>("maxDeltaPhi_EC", 0.14);
0207   desc.add<double>("EC_weight", 0.008);
0208   desc.add<double>("weight_dPhi_EC", 0.064516129);
0209   desc.add<double>("zClusterWidth_step1", 2.0);
0210   desc.add<double>("zClusterWidth_step2", 0.65);
0211   desc.add<double>("zClusterSearchArea_step2", 3.0);
0212   desc.add<double>("weightCut_step2", 0.05);
0213   desc.add<double>("zClusterWidth_step3", 0.3);
0214   desc.add<double>("zClusterSearchArea_step3", 0.55);
0215   desc.add<double>("weightCut_step3", 0.1);
0216   desc.add<bool>("ptWeighting", false);  // <---- newMethod
0217   desc.add<double>("ptWeighting_slope", 1 / 20.);
0218   desc.add<double>("ptWeighting_offset", -1);
0219   descriptions.add("fastPrimaryVertexWithWeightsProducer", desc);
0220 }
0221
0222 void FastPrimaryVertexWithWeightsProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0223   using namespace edm;
0224   using namespace reco;
0225   using namespace std;
0226
0227   const float barrel_lenght = 30;  //half lenght of the pixel barrel 30 cm
0228
0229   //get pixel cluster
0230   Handle<SiPixelClusterCollectionNew> cH;
0231   iEvent.getByToken(clustersToken, cH);
0232   const SiPixelClusterCollectionNew& pixelClusters = *cH.product();
0233
0234   //get jets
0235   Handle<edm::View<reco::Jet> > jH;
0236   iEvent.getByToken(jetsToken, jH);
0237   const edm::View<reco::Jet>& jets = *jH.product();
0238
0239   vector<const reco::Jet*> selectedJets;
0240   int countjet = 0;
0241   for (edm::View<reco::Jet>::const_iterator it = jets.begin(); it != jets.end() && countjet < m_njets; it++) {
0242     if (  //select jets used in barrel or endcap pixel cluster projections
0243         ((it->pt() >= m_minJetPt) && std::abs(it->eta()) <= m_maxJetEta) ||  //barrel
0244         ((it->pt() >= m_minJetPt) && std::abs(it->eta()) <= m_maxJetEta_EC &&
0245          std::abs(it->eta()) >= m_minJetEta_EC)  //endcap
0246     ) {
0247       selectedJets.push_back(&(*it));
0248       countjet++;
0249     }
0250   }
0251
0252   //get PixelClusterParameterEstimator
0253   const PixelClusterParameterEstimator* pp = &iSetup.getData(m_pixelCPEToken);
0254
0255   //get beamSpot
0256   edm::Handle<BeamSpot> beamSpot;
0257   iEvent.getByToken(beamSpotToken, beamSpot);
0258
0259   //get TrackerGeometry
0260   const TrackerGeometry* trackerGeometry = &iSetup.getData(m_geomToken);
0261
0262   // PART I: get z-projections with z-weights
0263   std::vector<float> zProjections;
0264   std::vector<float> zWeights;
0265   int jet_count = 0;
0266   for (vector<const reco::Jet*>::iterator jit = selectedJets.begin(); jit != selectedJets.end();
0267        jit++) {  //loop on selected jets
0268     float px = (*jit)->px();
0269     float py = (*jit)->py();
0270     float pz = (*jit)->pz();
0271     float pt = (*jit)->pt();
0272     float eta = (*jit)->eta();
0273     float jetZOverRho = (*jit)->momentum().Z() / (*jit)->momentum().Rho();
0274     float pt_weight = pt * m_ptWeighting_slope + m_ptWeighting_offset;
0275     for (SiPixelClusterCollectionNew::const_iterator it = pixelClusters.begin(); it != pixelClusters.end();
0276          it++)  //Loop on pixel modules with clusters
0277     {           //loop on pixel modules
0278       DetId id = it->detId();
0279       const edmNew::DetSet<SiPixelCluster>& detset = (*it);
0280       Point3DBase<float, GlobalTag> modulepos = trackerGeometry->idToDet(id)->position();
0281       float zmodule = modulepos.z() -
0282                       ((modulepos.x() - beamSpot->x0()) * px + (modulepos.y() - beamSpot->y0()) * py) / pt * pz / pt;
0283       if ((std::abs(deltaPhi((*jit)->momentum().Phi(), modulepos.phi())) < m_maxDeltaPhi * 2) &&
0284           (std::abs(zmodule) < (m_maxZ + barrel_lenght))) {  //if it is a compatible module
0285         for (size_t j = 0; j < detset.size(); j++) {         //loop on pixel clusters on this module
0286           const SiPixelCluster& aCluster = detset[j];
0287           if (   //it is a cluster to project
0288               (  // barrel
0289                   m_barrel && std::abs(modulepos.z()) < barrel_lenght && pt >= m_minJetPt && jet_count < m_njets &&
0290                   std::abs(eta) <= m_maxJetEta && aCluster.sizeX() <= m_maxSizeX &&
0291                   aCluster.sizeY() >= m_PixelCellHeightOverWidth * std::abs(jetZOverRho) + m_minSizeY_q &&
0292                   aCluster.sizeY() <= m_PixelCellHeightOverWidth * std::abs(jetZOverRho) + m_maxSizeY_q) ||
0293               (  // EC
0294                   m_endCap && std::abs(modulepos.z()) > barrel_lenght && pt > m_minJetPt && jet_count < m_njets &&
0295                   std::abs(eta) <= m_maxJetEta_EC && std::abs(eta) >= m_minJetEta_EC &&
0296                   aCluster.sizeX() <= m_maxSizeX)) {
0297             Point3DBase<float, GlobalTag> v = trackerGeometry->idToDet(id)->surface().toGlobal(
0298                 pp->localParametersV(aCluster, (*trackerGeometry->idToDetUnit(id)))[0].first);
0299             GlobalPoint v_bs(v.x() - beamSpot->x0(), v.y() - beamSpot->y0(), v.z());
0300             if (  //it pass DeltaPhi(Jet,Cluster) requirements
0301                 (m_barrel && std::abs(modulepos.z()) < barrel_lenght &&
0302                  std::abs(deltaPhi((*jit)->momentum().Phi(), v_bs.phi())) <= m_maxDeltaPhi) ||  //barrel
0303                 (m_endCap && std::abs(modulepos.z()) > barrel_lenght &&
0304                  std::abs(deltaPhi((*jit)->momentum().Phi(), v_bs.phi())) <= m_maxDeltaPhi_EC)  //EC
0305             ) {
0306               //calculate z-projection
0307               float z = v.z() - ((v.x() - beamSpot->x0()) * px + (v.y() - beamSpot->y0()) * py) / pt * pz / pt;
0308               if (std::abs(z) < m_maxZ) {
0309                 zProjections.push_back(z);  //add z-projection in zProjections
0310                 float weight = 0;
0311                 //calculate zWeight
0312                 if (std::abs(modulepos.z()) < barrel_lenght) {  //barrel
0313                                                                 //calculate weight_sizeY
0314                   float sizeY = aCluster.sizeY();
0315                   float sizeY_up = m_PixelCellHeightOverWidth * std::abs(jetZOverRho) + m_maxSizeY_q;
0316                   float sizeY_peak = m_PixelCellHeightOverWidth * std::abs(jetZOverRho) + m_peakSizeY_q;
0317                   float sizeY_down = m_PixelCellHeightOverWidth * std::abs(jetZOverRho) + m_minSizeY_q;
0318                   float weight_sizeY_up = (sizeY_up - sizeY) / (sizeY_up - sizeY_peak);
0319                   float weight_sizeY_down = (sizeY - sizeY_down) / (sizeY_peak - sizeY_down);
0320                   weight_sizeY_down = weight_sizeY_down * (weight_sizeY_down > 0) * (weight_sizeY_down < 1);
0321                   weight_sizeY_up = weight_sizeY_up * (weight_sizeY_up > 0) * (weight_sizeY_up < 1);
0322                   float weight_sizeY = weight_sizeY_up + weight_sizeY_down;
0323
0324                   //calculate weight_rho
0325                   float rho = sqrt(v_bs.x() * v_bs.x() + v_bs.y() * v_bs.y());
0326                   float weight_rho = ((m_weight_rho_up - rho) / m_weight_rho_up);
0327
0328                   //calculate weight_dPhi
0329                   float weight_dPhi = exp(-std::abs(deltaPhi((*jit)->momentum().Phi(), v_bs.phi())) / m_weight_dPhi);
0330
0331                   //calculate weight_sizeX1
0332                   float weight_sizeX1 = (aCluster.sizeX() == 2) + (aCluster.sizeX() == 1) * m_weight_SizeX1;
0333
0334                   //calculate weight_charge
0335                   float charge = aCluster.charge();
0336                   float weightCluster_up = (m_weight_charge_up - charge) / (m_weight_charge_up - m_weight_charge_peak);
0337                   float weightCluster_down =
0338                       (charge - m_weight_charge_down) / (m_weight_charge_peak - m_weight_charge_down);
0339                   weightCluster_down = weightCluster_down * (weightCluster_down > 0) * (weightCluster_down < 1);
0340                   weightCluster_up = weightCluster_up * (weightCluster_up > 0) * (weightCluster_up < 1);
0341                   float weight_charge = weightCluster_up + weightCluster_down;
0342
0343                   //calculate the final weight
0344                   weight = weight_dPhi * weight_sizeY * weight_rho * weight_sizeX1 * weight_charge;
0345                 } else if (std::abs(modulepos.z()) > barrel_lenght)  // EC
0346                 {                                                    // EC
0347                                                                      //calculate weight_dPhi
0348                   float weight_dPhi = exp(-std::abs(deltaPhi((*jit)->momentum().Phi(), v_bs.phi())) / m_weight_dPhi_EC);
0349                   //calculate the final weight
0350                   weight = m_EC_weight * (weight_dPhi);
0351                 }
0352                 if (m_ptWeighting)
0353                   weight = weight * pt_weight;
0354                 zWeights.push_back(weight);  //add the weight to zWeights
0355               }
0356             }  //if it pass DeltaPhi(Jet,Cluster) requirements
0357           }  ///if it is a cluster to project
0358         }  //loop on pixel clusters on this module
0359       }  //if it is a compatible module
0360     }  //loop on pixel modules
0361     jet_count++;
0362   }  //loop on selected jets
0363
0364   //order zProjections and zWeights by z
0365   std::multimap<float, float> zWithW;
0366   size_t i = 0;
0367   for (i = 0; i < zProjections.size(); i++)
0368     zWithW.insert(std::pair<float, float>(zProjections[i], zWeights[i]));
0369   i = 0;
0370   for (std::multimap<float, float>::iterator it = zWithW.begin(); it != zWithW.end(); it++, i++) {
0371     zProjections[i] = it->first;
0372     zWeights[i] = it->second;
0373   }  //order zProjections and zWeights by z
0374
0375   //calculate zWeightsSquared
0376   std::vector<float> zWeightsSquared;
0377   for (std::vector<float>::iterator it = zWeights.begin(); it != zWeights.end(); it++) {
0378     zWeightsSquared.push_back((*it) * (*it));
0379   }
0380
0381   //do multi-step peak searching
0382   float res_step1 = FindPeakFastPV(zProjections, zWeights, 0.0, m_zClusterWidth_step1, 999.0, -1.0);
0383   float res_step2 = FindPeakFastPV(
0384       zProjections, zWeights, res_step1, m_zClusterWidth_step2, m_zClusterSearchArea_step2, m_weightCut_step2);
0385   float res_step3 = FindPeakFastPV(zProjections,
0386                                    zWeightsSquared,
0387                                    res_step2,
0388                                    m_zClusterWidth_step3,
0389                                    m_zClusterSearchArea_step3,
0390                                    m_weightCut_step3 * m_weightCut_step3);
0391
0392   float centerWMax = res_step3;
0393   //End of PART II
0394
0395   //Make the output
0396   float res = 0;
0397   if (zProjections.size() > 2) {
0398     res = centerWMax;
0399     Vertex::Error e;
0400     e(0, 0) = 0.0015 * 0.0015;
0401     e(1, 1) = 0.0015 * 0.0015;
0402     e(2, 2) = 1.5 * 1.5;
0403     Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
0404     Vertex thePV(p, e, 1, 1, 0);
0405     auto pOut = std::make_unique<reco::VertexCollection>();
0406     pOut->push_back(thePV);
0407     iEvent.put(std::move(pOut));
0408   } else {
0409     Vertex::Error e;
0410     e(0, 0) = 0.0015 * 0.0015;
0411     e(1, 1) = 0.0015 * 0.0015;
0412     e(2, 2) = 1.5 * 1.5;
0413     Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
0414     Vertex thePV(p, e, 0, 0, 0);
0415     auto pOut = std::make_unique<reco::VertexCollection>();
0416     pOut->push_back(thePV);
0417     iEvent.put(std::move(pOut));
0418   }
0419
0420   //Finally, calculate the zClusterQuality as Sum(weights near the fastPV)/sqrt(Sum(total weights)) [a kind of zCluster significance]
0421
0422   const float half_width_peak = 1;
0423   float nWeightedTot = 0;
0424   float nWeightedTotPeak = 0;
0425   for (std::vector<float>::iterator it = zProjections.begin(); it != zProjections.end(); it++) {
0426     nWeightedTot += zWeights[it - zProjections.begin()];
0427     if ((res - half_width_peak) <= (*it) && (*it) <= (res + half_width_peak)) {
0428       nWeightedTotPeak += zWeights[it - zProjections.begin()];
0429     }
0430   }
0431
0432   auto zClusterQuality = std::make_unique<float>();
0433   *zClusterQuality = -1;
0434   if (nWeightedTot != 0) {
0435     // where 30 is the beam spot lenght
0436     *zClusterQuality = nWeightedTotPeak / sqrt(nWeightedTot / (2 * half_width_peak));
0437     iEvent.put(std::move(zClusterQuality));
0438   } else
0439     iEvent.put(std::move(zClusterQuality));
0440 }
0441
0442 //define this as a plug-in
0443 DEFINE_FWK_MODULE(FastPrimaryVertexWithWeightsProducer);