TrackingMCTruth/plugins/SimDoubletsAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Validation/TrackingMCTruth
0004 // Class:      SimDoubletsAnalyzer
0005 //
0006
0007 // user include files
0008 #include "Validation/TrackingMCTruth/plugins/SimDoubletsAnalyzer.h"
0009
0010 #include "DataFormats/Histograms/interface/MonitorElementCollection.h"
0011 #include "DataFormats/Math/interface/deltaPhi.h"
0012 #include "DataFormats/Math/interface/approx_atan2.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0014 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0015 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018
0019 #include <cstddef>
0020
0021 namespace simdoublets {
0022   // helper function that takes the layerPairId and returns two strings with the
0023   // inner and outer layer id
0024   std::pair<std::string, std::string> getInnerOuterLayerNames(int const layerPairId) {
0025     // make a string from the Id (int)
0026     std::string index = std::to_string(layerPairId);
0027     // determine inner and outer layer name
0028     std::string innerLayerName;
0029     std::string outerLayerName;
0030     if (index.size() < 3) {
0031       innerLayerName = "0";
0032       outerLayerName = index;
0033     } else if (index.size() == 3) {
0034       innerLayerName = index.substr(0, 1);
0035       outerLayerName = index.substr(1, 3);
0036     } else {
0037       innerLayerName = index.substr(0, 2);
0038       outerLayerName = index.substr(2, 4);
0039     }
0040     if (outerLayerName[0] == '0') {
0041       outerLayerName = outerLayerName.substr(1, 2);
0042     }
0043
0044     return {innerLayerName, outerLayerName};
0045   }
0046
0047   // make bins logarithmic
0048   void BinLogX(TH1* h) {
0049     TAxis* axis = h->GetXaxis();
0050     int bins = axis->GetNbins();
0051
0052     float from = axis->GetXmin();
0053     float to = axis->GetXmax();
0054     float width = (to - from) / bins;
0055     std::vector<float> new_bins(bins + 1, 0);
0056
0057     for (int i = 0; i <= bins; i++) {
0058       new_bins[i] = TMath::Power(10, from + i * width);
0059     }
0060     axis->Set(bins, new_bins.data());
0061   }
0062
0063   // function to produce histogram with log scale on x (taken from MultiTrackValidator)
0064   template <typename... Args>
0065   dqm::reco::MonitorElement* make1DLogX(dqm::reco::DQMStore::IBooker& ibook, Args&&... args) {
0066     auto h = std::make_unique<TH1F>(std::forward<Args>(args)...);
0067     BinLogX(h.get());
0068     const auto& name = h->GetName();
0069     return ibook.book1D(name, h.release());
0070   }
0071
0072   // function to produce profile with log scale on x (taken from MultiTrackValidator)
0073   template <typename... Args>
0074   dqm::reco::MonitorElement* makeProfileLogX(dqm::reco::DQMStore::IBooker& ibook, Args&&... args) {
0075     auto h = std::make_unique<TProfile>(std::forward<Args>(args)...);
0076     BinLogX(h.get());
0077     const auto& name = h->GetName();
0078     return ibook.bookProfile(name, h.release());
0079   }
0080
0081   // function that checks if two vector share a common element
0082   template <typename T>
0083   bool haveCommonElement(std::vector<T> const& v1, std::vector<T> const& v2) {
0084     return std::find_first_of(v1.begin(), v1.end(), v2.begin(), v2.end()) != v1.end();
0085   }
0086
0087   // fillDescriptionsCommon: description that is identical for Phase 1 and 2
0088   template <typename TrackerTraits>
0089   void fillDescriptionsCommon(edm::ParameterSetDescription& desc) {
0090     desc.add<std::string>("folder", "Tracking/TrackingMCTruth/SimDoublets");
0091
0092     // cut parameters with scalar values
0093     desc.add<int>("cellMaxDYSize12", TrackerTraits::maxDYsize12)
0094         ->setComment("Maximum difference in cluster size for B1/B2");
0095     desc.add<int>("cellMaxDYSize", TrackerTraits::maxDYsize)->setComment("Maximum difference in cluster size");
0096     desc.add<int>("cellMaxDYPred", TrackerTraits::maxDYPred)
0097         ->setComment("Maximum difference between actual and expected cluster size of inner RecHit");
0098   }
0099 }  // namespace simdoublets
0100
0101 // -------------------------------------------------------------------------------------------------------------
0102 // constructors and destructor
0103 // -------------------------------------------------------------------------------------------------------------
0104
0105 template <typename TrackerTraits>
0106 SimDoubletsAnalyzer<TrackerTraits>::SimDoubletsAnalyzer(const edm::ParameterSet& iConfig)
0107     : topology_getToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0108       simDoublets_getToken_(consumes(iConfig.getParameter<edm::InputTag>("simDoubletsSrc"))),
0109       cellMinz_(iConfig.getParameter<std::vector<double>>("cellMinz")),
0110       cellMaxz_(iConfig.getParameter<std::vector<double>>("cellMaxz")),
0111       cellPhiCuts_(iConfig.getParameter<std::vector<int>>("cellPhiCuts")),
0112       cellMaxr_(iConfig.getParameter<std::vector<double>>("cellMaxr")),
0113       cellMinYSizeB1_(iConfig.getParameter<int>("cellMinYSizeB1")),
0114       cellMinYSizeB2_(iConfig.getParameter<int>("cellMinYSizeB2")),
0115       cellMaxDYSize12_(iConfig.getParameter<int>("cellMaxDYSize12")),
0116       cellMaxDYSize_(iConfig.getParameter<int>("cellMaxDYSize")),
0117       cellMaxDYPred_(iConfig.getParameter<int>("cellMaxDYPred")),
0118       cellZ0Cut_(iConfig.getParameter<double>("cellZ0Cut")),
0119       cellPtCut_(iConfig.getParameter<double>("cellPtCut")),
0120       folder_(iConfig.getParameter<std::string>("folder")) {
0121   // get layer pairs from configuration
0122   std::vector<int> layerPairs{iConfig.getParameter<std::vector<int>>("layerPairs")};
0123
0124   // number of configured layer pairs
0125   size_t numLayerPairs = layerPairs.size() / 2;
0126
0127   // fill the map of layer pairs
0128   for (size_t i{0}; i < numLayerPairs; i++) {
0129     int layerPairId = 100 * layerPairs[2 * i] + layerPairs[2 * i + 1];
0130     layerPairId2Index_.insert({layerPairId, i});
0131   }
0132
0133   // resize all histogram vectors, so that we can fill them according to the
0134   // layerPairIndex saved in the map that we just created
0135   hVector_dr_.resize(numLayerPairs);
0136   hVector_dphi_.resize(numLayerPairs);
0137   hVector_idphi_.resize(numLayerPairs);
0138   hVector_innerZ_.resize(numLayerPairs);
0139   hVector_Ysize_.resize(numLayerPairs);
0140   hVector_DYsize_.resize(numLayerPairs);
0141   hVector_DYPred_.resize(numLayerPairs);
0142   hVector_pass_dr_.resize(numLayerPairs);
0143   hVector_pass_idphi_.resize(numLayerPairs);
0144   hVector_pass_innerZ_.resize(numLayerPairs);
0145 }
0146
0147 template <typename TrackerTraits>
0148 SimDoubletsAnalyzer<TrackerTraits>::~SimDoubletsAnalyzer() {}
0149
0150 // -------------------------------------------------------------------------------------------------------------
0151 // member functions
0152 // -------------------------------------------------------------------------------------------------------------
0153
0154 template <typename TrackerTraits>
0155 void SimDoubletsAnalyzer<TrackerTraits>::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0156
0157 template <typename TrackerTraits>
0158 void SimDoubletsAnalyzer<TrackerTraits>::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0159   using namespace edm;
0160
0161   // get tracker topology
0162   trackerTopology_ = &iSetup.getData(topology_getToken_);
0163
0164   // get simDoublets
0165   SimDoubletsCollection const& simDoubletsCollection = iEvent.get(simDoublets_getToken_);
0166
0167   // create vectors for inner and outer RecHits of SimDoublets passing all cuts
0168   std::vector<SiPixelRecHitRef> innerRecHitsPassing;
0169   std::vector<SiPixelRecHitRef> outerRecHitsPassing;
0170
0171   // initialize a bunch of variables that we will use in the coming for loops
0172   double true_pT, true_eta, inner_r, inner_z, inner_phi, outer_r, outer_z, outer_phi, dz, dr, dphi, z0, curvature, pT;
0173   int numSimDoublets, pass_numSimDoublets, inner_iphi, outer_iphi, idphi, layerPairId, layerPairIdIndex,
0174       innerClusterSizeY;
0175   bool doubletGetsCut, subjectToYsizeB1, subjectToYsizeB2, subjectToDYsize, subjectToDYsize12, subjectToDYPred;
0176
0177   // loop over SimDoublets (= loop over TrackingParticles)
0178   for (auto const& simDoublets : simDoubletsCollection) {
0179     // get true pT of the TrackingParticle
0180     true_pT = simDoublets.trackingParticle()->pt();
0181     true_eta = simDoublets.trackingParticle()->eta();
0182
0183     // create the true RecHit doublets of the TrackingParticle
0184     auto doublets = simDoublets.getSimDoublets(trackerTopology_);
0185
0186     // number of SimDoublets of the Tracking Particle
0187     numSimDoublets = doublets.size();
0188     // number of SimDoublets of the Tracking Particle passing all cuts
0189     pass_numSimDoublets = 0;
0190
0191     // fill histograms for number of SimDoublets
0192     h_numSimDoubletsPerTrackingParticle_->Fill(numSimDoublets);
0193     h_numLayersPerTrackingParticle_->Fill(simDoublets.numLayers());
0194
0195     // fill histograms for number of TrackingParticles
0196     h_numTPVsPt_->Fill(true_pT);
0197     h_numTPVsEta_->Fill(true_eta);
0198
0199     // clear passing inner and outer RecHits
0200     innerRecHitsPassing.clear();
0201     outerRecHitsPassing.clear();
0202
0203     // loop over those doublets
0204     for (auto const& doublet : doublets) {
0205       // RecHit properties
0206       inner_r = doublet.innerGlobalPos().perp();
0207       inner_z = doublet.innerGlobalPos().z();
0208       inner_phi = doublet.innerGlobalPos().barePhi();  // returns float, whereas .phi() returns phi object
0209       inner_iphi = unsafe_atan2s<7>(doublet.innerGlobalPos().y(), doublet.innerGlobalPos().x());
0210       outer_r = doublet.outerGlobalPos().perp();
0211       outer_z = doublet.outerGlobalPos().z();
0212       outer_phi = doublet.outerGlobalPos().barePhi();
0213       outer_iphi = unsafe_atan2s<7>(doublet.outerGlobalPos().y(), doublet.outerGlobalPos().x());
0214
0215       dz = outer_z - inner_z;
0216       dr = outer_r - inner_r;
0217       dphi = reco::deltaPhi(inner_phi, outer_phi);
0218       idphi = std::min(std::abs(int16_t(outer_iphi - inner_iphi)), std::abs(int16_t(inner_iphi - outer_iphi)));
0219
0220       // ----------------------------------------------------------
0221       // general plots (general folder)
0222       // ----------------------------------------------------------
0223
0224       // outer layer vs inner layer of SimDoublets
0225       h_layerPairs_->Fill(doublet.innerLayerId(), doublet.outerLayerId());
0226
0227       // number of skipped layers by SimDoublets
0228       h_numSkippedLayers_->Fill(doublet.numSkippedLayers());
0229
0230       // ----------------------------------------------------------
0231       // layer pair independent cuts (global folder)
0232       // ----------------------------------------------------------
0233
0234       // longitudinal impact parameter with respect to the beamspot
0235       z0 = std::abs(inner_r * outer_z - inner_z * outer_r) / dr;
0236       h_z0_->Fill(z0);
0237
0238       // radius of the circle defined by the two RecHits and the beamspot
0239       curvature = 1.f / 2.f * std::sqrt((dr / dphi) * (dr / dphi) + (inner_r * outer_r));
0240       h_curvatureR_->Fill(curvature);
0241
0242       // pT that this curvature radius corresponds to
0243       pT = curvature / 87.78f;
0244       h_pTFromR_->Fill(pT);
0245
0246       // ----------------------------------------------------------
0247       // layer pair dependent cuts (sub-folders for layer pairs)
0248       // ----------------------------------------------------------
0249
0250       // first, get layer pair Id and exclude layer pairs that are not considered
0251       layerPairId = doublet.layerPairId();
0252       if (layerPairId2Index_.find(layerPairId) == layerPairId2Index_.end()) {
0253         continue;
0254       }
0255
0256       // get the position of the layer pair in the vectors of histograms
0257       layerPairIdIndex = layerPairId2Index_.at(layerPairId);
0258
0259       // dr = (outer_r - inner_r) histogram
0260       hVector_dr_[layerPairIdIndex]->Fill(dr);
0261
0262       // dphi histogram
0263       hVector_dphi_[layerPairIdIndex]->Fill(dphi);
0264       hVector_idphi_[layerPairIdIndex]->Fill(idphi);
0265
0266       // z of the inner RecHit histogram
0267       hVector_innerZ_[layerPairIdIndex]->Fill(inner_z);
0268
0269       // ----------------------------------------------------------
0270       // cluster size cuts (global + sub-folders for layer pairs)
0271       // ----------------------------------------------------------
0272
0273       // cluster size in local y histogram
0274       innerClusterSizeY = doublet.innerRecHit()->cluster()->sizeY();
0275       hVector_Ysize_[layerPairIdIndex]->Fill(innerClusterSizeY);
0276
0277       // create bool that indicates if the doublet gets cut
0278       doubletGetsCut = false;
0279       // create bools that trace if doublet is subject to any clsuter size cut
0280       subjectToYsizeB1 = false;
0281       subjectToYsizeB2 = false;
0282       subjectToDYsize = false;
0283       subjectToDYsize12 = false;
0284       subjectToDYPred = false;
0285
0286       // apply all cuts that do not depend on the cluster size
0287       // z window cut
0288       if (inner_z < cellMinz_[layerPairIdIndex] || inner_z > cellMaxz_[layerPairIdIndex]) {
0289         doubletGetsCut = true;
0290       }
0291       // z0cutoff
0292       if (dr > cellMaxr_[layerPairIdIndex] || dr < 0 || z0 > cellZ0Cut_) {
0293         doubletGetsCut = true;
0294       }
0295       // ptcut
0296       if (pT < cellPtCut_) {
0297         doubletGetsCut = true;
0298       }
0299       // iphicut
0300       if (idphi > cellPhiCuts_[layerPairIdIndex]) {
0301         doubletGetsCut = true;
0302       }
0303
0304       // determine the moduleId
0305       const GeomDetUnit* geomDetUnit = doublet.innerRecHit()->det();
0306       const uint32_t moduleId = geomDetUnit->index();
0307
0308       // define bools needed to decide on cutting parameters
0309       const bool innerInB1 = (doublet.innerLayerId() == 0);
0310       const bool innerInB2 = (doublet.innerLayerId() == 1);
0311       const bool isOuterLadder = (0 == (moduleId / 8) % 2);  // check if this even makes sense in Phase-2
0312       const bool innerInBarrel = (doublet.innerLayerId() < 4);
0313       const bool outerInBarrel = (doublet.outerLayerId() < 4);
0314
0315       // histograms for clusterCut
0316       // cluster size in local y
0317       if (!outerInBarrel) {
0318         if (innerInB1 && isOuterLadder) {
0319           subjectToYsizeB1 = true;
0320           h_YsizeB1_->Fill(innerClusterSizeY);
0321           // apply the cut
0322           if (innerClusterSizeY < cellMinYSizeB1_) {
0323             doubletGetsCut = true;
0324           }
0325         }
0326         if (innerInB2) {
0327           subjectToYsizeB2 = true;
0328           h_YsizeB2_->Fill(innerClusterSizeY);
0329           // apply the cut
0330           if (innerClusterSizeY < cellMinYSizeB2_) {
0331             doubletGetsCut = true;
0332           }
0333         }
0334       }
0335
0336       // histograms for zSizeCut
0337       int DYsize{0}, DYPred{0};
0338       if (innerInBarrel) {
0339         if (outerInBarrel) {  // onlyBarrel
0340           DYsize = std::abs(innerClusterSizeY - doublet.outerRecHit()->cluster()->sizeY());
0341           if (innerInB1 && isOuterLadder) {
0342             subjectToDYsize12 = true;
0343             hVector_DYsize_[layerPairIdIndex]->Fill(DYsize);
0344             h_DYsize12_->Fill(DYsize);
0345             // apply the cut
0346             if (DYsize > cellMaxDYSize12_) {
0347               doubletGetsCut = true;
0348             }
0349           } else if (!innerInB1) {
0350             subjectToDYsize = true;
0351             hVector_DYsize_[layerPairIdIndex]->Fill(DYsize);
0352             h_DYsize_->Fill(DYsize);
0353             // apply the cut
0354             if (DYsize > cellMaxDYSize_) {
0355               doubletGetsCut = true;
0356             }
0357           }
0358         } else {  // not onlyBarrel
0359           subjectToDYPred = true;
0360           DYPred = std::abs(innerClusterSizeY - int(std::abs(dz / dr) * pixelTopology::Phase2::dzdrFact + 0.5f));
0361           hVector_DYPred_[layerPairIdIndex]->Fill(DYPred);
0362           h_DYPred_->Fill(DYPred);
0363           // apply the cut
0364           if (DYPred > cellMaxDYPred_) {
0365             doubletGetsCut = true;
0366           }
0367         }
0368       }
0369
0370       // ----------------------------------------------------------
0371       // all kinds of plots for doublets passing all cuts
0372       // ----------------------------------------------------------
0373
0374       // fill the number histograms
0375       // histogram of all valid doublets
0376       h_numVsPt_->Fill(true_pT);
0377       h_numVsEta_->Fill(true_eta);
0378
0379       // if the doublet passes all cuts
0380       if (!doubletGetsCut) {
0381         // increment number of SimDoublets passing all cuts
0382         pass_numSimDoublets++;
0383
0384         // fill histogram of doublets that pass all cuts
0385         h_pass_layerPairs_->Fill(doublet.innerLayerId(), doublet.outerLayerId());
0386         h_pass_numVsPt_->Fill(true_pT);
0387         h_pass_numVsEta_->Fill(true_eta);
0388
0389         // also put the inner/outer RecHit in the respective vector
0390         innerRecHitsPassing.push_back(doublet.innerRecHit());
0391         outerRecHitsPassing.push_back(doublet.outerRecHit());
0392
0393         // fill pass_ histograms
0394         h_pass_z0_->Fill(z0);
0395         h_pass_pTFromR_->Fill(pT);
0396         hVector_pass_dr_[layerPairIdIndex]->Fill(dr);
0397         hVector_pass_idphi_[layerPairIdIndex]->Fill(idphi);
0398         hVector_pass_innerZ_[layerPairIdIndex]->Fill(inner_z);
0399         if (subjectToDYPred) {
0400           h_pass_DYPred_->Fill(DYPred);
0401         }
0402         if (subjectToDYsize) {
0403           h_pass_DYsize_->Fill(DYsize);
0404         }
0405         if (subjectToDYsize12) {
0406           h_pass_DYsize12_->Fill(DYsize);
0407         }
0408         if (subjectToYsizeB1) {
0409           h_pass_YsizeB1_->Fill(innerClusterSizeY);
0410         }
0411         if (subjectToYsizeB2) {
0412           h_pass_YsizeB2_->Fill(innerClusterSizeY);
0413         }
0414       }
0415     }  // end loop over those doublets
0416
0417     // Now check if the TrackingParticle is reconstructable by at least two conencted SimDoublets surviving the cuts
0418     if (simdoublets::haveCommonElement<SiPixelRecHitRef>(innerRecHitsPassing, outerRecHitsPassing)) {
0419       h_pass_numTPVsPt_->Fill(true_pT);
0420       h_pass_numTPVsEta_->Fill(true_eta);
0421     }
0422
0423     // Fill the efficiency profile per Tracking Particle only if the TP has at least one SimDoublet
0424     if (numSimDoublets > 0) {
0425       h_effSimDoubletsPerTPVsEta_->Fill(true_eta, pass_numSimDoublets / numSimDoublets);
0426       h_effSimDoubletsPerTPVsPt_->Fill(true_pT, pass_numSimDoublets / numSimDoublets);
0427     }
0428   }  // end loop over SimDoublets (= loop over TrackingParticles)
0429 }
0430
0431 // booking the histograms
0432 template <typename TrackerTraits>
0433 void SimDoubletsAnalyzer<TrackerTraits>::bookHistograms(DQMStore::IBooker& ibook,
0434                                                         edm::Run const& run,
0435                                                         edm::EventSetup const& iSetup) {
0436   // set some common parameters
0437   int pTNBins = 50;
0438   double pTmin = log10(0.01);
0439   double pTmax = log10(1000);
0440   int etaNBins = 80;
0441   double etamin = -4.;
0442   double etamax = 4.;
0443
0444   // ----------------------------------------------------------
0445   // booking general histograms (general folder)
0446   // ----------------------------------------------------------
0447
0448   ibook.setCurrentFolder(folder_ + "/general");
0449
0450   // overview histograms and profiles
0451   h_effSimDoubletsPerTPVsPt_ =
0452       simdoublets::makeProfileLogX(ibook,
0453                                    "efficiencyPerTP_vs_pT",
0454                                    "SimDoublets efficiency per TP vs p_{T}; TP transverse momentum p_{T} [GeV]; "
0455                                    "Average fraction of SimDoublets per TP passing all cuts",
0456                                    pTNBins,
0457                                    pTmin,
0458                                    pTmax,
0459                                    0,
0460                                    1,
0461                                    " ");
0462   h_effSimDoubletsPerTPVsEta_ = ibook.bookProfile("efficiencyPerTP_vs_eta",
0463                                                   "SimDoublets efficiency per TP vs #eta; TP transverse momentum #eta; "
0464                                                   "Average fraction of SimDoublets per TP passing all cuts",
0465                                                   etaNBins,
0466                                                   etamin,
0467                                                   etamax,
0468                                                   0,
0469                                                   1,
0470                                                   " ");
0471   h_layerPairs_ = ibook.book2D("layerPairs",
0472                                "Layer pairs in SimDoublets; Inner layer ID; Outer layer ID",
0473                                TrackerTraits::numberOfLayers,
0474                                -0.5,
0475                                -0.5 + TrackerTraits::numberOfLayers,
0476                                TrackerTraits::numberOfLayers,
0477                                -0.5,
0478                                -0.5 + TrackerTraits::numberOfLayers);
0479   h_pass_layerPairs_ = ibook.book2D("pass_layerPairs",
0480                                     "Layer pairs in SimDoublets passing all cuts; Inner layer ID; Outer layer ID",
0481                                     TrackerTraits::numberOfLayers,
0482                                     -0.5,
0483                                     -0.5 + TrackerTraits::numberOfLayers,
0484                                     TrackerTraits::numberOfLayers,
0485                                     -0.5,
0486                                     -0.5 + TrackerTraits::numberOfLayers);
0487   h_numSkippedLayers_ = ibook.book1D(
0488       "numSkippedLayers", "Number of skipped layers; Number of skipped layers; Number of SimDoublets", 16, -1.5, 14.5);
0489   h_numSimDoubletsPerTrackingParticle_ =
0490       ibook.book1D("numSimDoubletsPerTrackingParticle",
0491                    "Number of SimDoublets per Tracking Particle; Number of SimDoublets; Number of Tracking Particles",
0492                    31,
0493                    -0.5,
0494                    30.5);
0495   h_numLayersPerTrackingParticle_ =
0496       ibook.book1D("numLayersPerTrackingParticle",
0497                    "Number of layers hit by Tracking Particle; Number of layers; Number of Tracking Particles",
0498                    29,
0499                    -0.5,
0500                    28.5);
0501   h_numTPVsPt_ = simdoublets::make1DLogX(
0502       ibook,
0503       "numTPVsPt",
0504       "Total number of TrackingParticles; True transverse momentum p_{T} [GeV]; Total number of TrackingParticles",
0505       pTNBins,
0506       pTmin,
0507       pTmax);
0508   h_pass_numTPVsPt_ = simdoublets::make1DLogX(ibook,
0509                                               "pass_numTPVsPt",
0510                                               "Reconstructable TrackingParticles (two or more connected SimDoublets "
0511                                               "pass cuts); True transverse momentum p_{T} [GeV]; "
0512                                               "Number of reconstructable TrackingParticles",
0513                                               pTNBins,
0514                                               pTmin,
0515                                               pTmax);
0516   h_numTPVsEta_ =
0517       ibook.book1D("numTPVsEta",
0518                    "Total number of TrackingParticles; True pseudorapidity #eta; Total number of TrackingParticles",
0519                    etaNBins,
0520                    etamin,
0521                    etamax);
0522   h_pass_numTPVsEta_ = ibook.book1D("pass_numTPVsEta",
0523                                     "Reconstructable TrackingParticles (two or more connected SimDoublets "
0524                                     "pass cuts); True pseudorapidity #eta; Number of reconstructable TrackingParticles",
0525                                     etaNBins,
0526                                     etamin,
0527                                     etamax);
0528   h_numVsPt_ = simdoublets::make1DLogX(
0529       ibook,
0530       "numVsPt",
0531       "Total number of SimDoublets; True transverse momentum p_{T} [GeV]; Total number of SimDoublets",
0532       pTNBins,
0533       pTmin,
0534       pTmax);
0535   h_pass_numVsPt_ = simdoublets::make1DLogX(ibook,
0536                                             "pass_numVsPt",
0537                                             "Number of passing SimDoublets; True transverse momentum p_{T} [GeV]; "
0538                                             "Number of SimDoublets passing all cuts",
0539                                             pTNBins,
0540                                             pTmin,
0541                                             pTmax);
0542   h_numVsEta_ = ibook.book1D("numVsEta",
0543                              "Total number of SimDoublets; True pseudorapidity #eta; Total number of SimDoublets",
0544                              etaNBins,
0545                              etamin,
0546                              etamax);
0547   h_pass_numVsEta_ =
0548       ibook.book1D("pass_numVsEta",
0549                    "Number of SimDoublets; True pseudorapidity #eta; Number of SimDoublets passing all cuts",
0550                    etaNBins,
0551                    etamin,
0552                    etamax);
0553
0554   // -------------------------------------------------------------
0555   // booking layer pair independent cut histograms (global folder)
0556   // -------------------------------------------------------------
0557
0558   ibook.setCurrentFolder(folder_ + "/cutParameters/global");
0559
0560   // histogram for z0cutoff  (z0Cut)
0561   h_z0_ = ibook.book1D("z0", "z_{0}; Longitudinal impact parameter z_{0} [cm]; Number of SimDoublets", 51, -1, 50);
0562   h_pass_z0_ = ibook.book1D(
0563       "pass_z0",
0564       "z_{0} of SimDoublets passing all cuts; Longitudinal impact parameter z_{0} [cm]; Number of SimDoublets",
0565       51,
0566       -1,
0567       50);
0568
0569   // histograms for ptcut  (ptCut)
0570   h_curvatureR_ = ibook.book1D(
0571       "curvatureR", "Curvature from SimDoublet+beamspot; Curvature radius [cm] ; Number of SimDoublets", 100, 0, 1000);
0572   h_pTFromR_ = simdoublets::make1DLogX(
0573       ibook,
0574       "pTFromR",
0575       "Transverse momentum from curvature; Transverse momentum p_{T} [GeV]; Number of SimDoublets",
0576       pTNBins,
0577       pTmin,
0578       pTmax);
0579   h_pass_pTFromR_ = simdoublets::make1DLogX(ibook,
0580                                             "pass_pTFromR",
0581                                             "Transverse momentum from curvature of SimDoublets passing all cuts; "
0582                                             "Transverse momentum p_{T} [GeV]; Number of SimDoublets",
0583                                             pTNBins,
0584                                             pTmin,
0585                                             pTmax);
0586
0587   // histograms for clusterCut  (minYsizeB1 and minYsizeB2)
0588   h_YsizeB1_ = ibook.book1D(
0589       "YsizeB1",
0590       "Cluster size along z (inner from B1); Size along z of inner cluster [num of pixels]; Number of SimDoublets",
0591       51,
0592       -1,
0593       50);
0594   h_YsizeB2_ = ibook.book1D(
0595       "YsizeB2",
0596       "Cluster size along z (inner not from B1); Size along z of inner cluster [num of pixels]; Number of SimDoublets",
0597       51,
0598       -1,
0599       50);
0600   h_pass_YsizeB1_ = ibook.book1D("pass_YsizeB1",
0601                                  "Cluster size along z of SimDoublets passing all cuts (inner from B1); Size along z "
0602                                  "of inner cluster [num of pixels]; Number of SimDoublets",
0603                                  51,
0604                                  -1,
0605                                  50);
0606   h_pass_YsizeB2_ = ibook.book1D("pass_YsizeB2",
0607                                  "Cluster size along z of SimDoublets passing all cuts (inner not from B1); Size along "
0608                                  "z of inner cluster [num of pixels]; Number of SimDoublets",
0609                                  51,
0610                                  -1,
0611                                  50);
0612
0613   // histograms for zSizeCut  (maxDYsize12, maxDYsize and maxDYPred)
0614   h_DYsize12_ =
0615       ibook.book1D("DYsize12",
0616                    "Difference in cluster size along z (inner from B1); Absolute difference in cluster size along z of "
0617                    "the two RecHits [num of pixels]; Number of SimDoublets",
0618                    31,
0619                    -1,
0620                    30);
0621   h_DYsize_ = ibook.book1D("DYsize",
0622                            "Difference in cluster size along z; Absolute difference in cluster size along z of the two "
0623                            "RecHits [num of pixels]; Number of SimDoublets",
0624                            31,
0625                            -1,
0626                            30);
0627   h_DYPred_ = ibook.book1D("DYPred",
0628                            "Difference between actual and predicted cluster size along z of inner cluster; Absolute "
0629                            "difference [num of pixels]; Number of SimDoublets",
0630                            201,
0631                            -1,
0632                            200);
0633   h_pass_DYsize12_ = ibook.book1D("pass_DYsize12",
0634                                   "Difference in cluster size along z of SimDoublets passing all cuts (inner from B1); "
0635                                   "Absolute difference in cluster size along z of "
0636                                   "the two RecHits [num of pixels]; Number of SimDoublets",
0637                                   31,
0638                                   -1,
0639                                   30);
0640   h_pass_DYsize_ =
0641       ibook.book1D("pass_DYsize",
0642                    "Difference in cluster size along z of SimDoublets passing all cuts; Absolute difference in "
0643                    "cluster size along z of the two RecHits [num of pixels]; Number of SimDoublets",
0644                    31,
0645                    -1,
0646                    30);
0647   h_pass_DYPred_ =
0648       ibook.book1D("pass_DYPred",
0649                    "Difference between actual and predicted cluster size along z of inner cluster of SimDoublets "
0650                    "passing all cuts; Absolute difference [num of pixels]; Number of SimDoublets",
0651                    201,
0652                    -1,
0653                    200);
0654
0655   // -----------------------------------------------------------------------
0656   // booking layer pair dependent histograms (sub-folders for layer pairs)
0657   // -----------------------------------------------------------------------
0658
0659   // loop through valid layer pairs and add for each one booked hist per vector
0660   for (auto id = layerPairId2Index_.begin(); id != layerPairId2Index_.end(); ++id) {
0661     // get the position of the layer pair in the histogram vectors
0662     int layerPairIdIndex = id->second;
0663
0664     // get layer names from the layer pair Id
0665     auto layerNames = simdoublets::getInnerOuterLayerNames(id->first);
0666     std::string innerLayerName = layerNames.first;
0667     std::string outerLayerName = layerNames.second;
0668
0669     // name the sub-folder for the layer pair "lp_${innerLayerId}_${outerLayerId}"
0670     std::string subFolderName = "/cutParameters/lp_" + innerLayerName + "_" + outerLayerName;
0671
0672     // layer mentioning in histogram titles
0673     std::string layerTitle = "(layers (" + innerLayerName + "," + outerLayerName + "))";
0674
0675     // set folder to the sub-folder for the layer pair
0676     ibook.setCurrentFolder(folder_ + subFolderName);
0677
0678     // histogram for z0cutoff  (maxr)
0679     hVector_dr_.at(layerPairIdIndex) = ibook.book1D(
0680         "dr",
0681         "dr of RecHit pair " + layerTitle + "; dr between outer and inner RecHit [cm]; Number of SimDoublets",
0682         31,
0683         -1,
0684         30);
0685     hVector_pass_dr_.at(layerPairIdIndex) = ibook.book1D(
0686         "pass_dr",
0687         "dr of RecHit pair " + layerTitle +
0688             " for SimDoublets passing all cuts; dr between outer and inner RecHit [cm]; Number of SimDoublets",
0689         31,
0690         -1,
0691         30);
0692
0693     // histograms for iphicut  (phiCuts)
0694     hVector_dphi_.at(layerPairIdIndex) = ibook.book1D(
0695         "dphi",
0696         "dphi of RecHit pair " + layerTitle + "; d#phi between outer and inner RecHit [rad]; Number of SimDoublets",
0697         50,
0698         -M_PI,
0699         M_PI);
0700     hVector_idphi_.at(layerPairIdIndex) =
0701         ibook.book1D("idphi",
0702                      "idphi of RecHit pair " + layerTitle +
0703                          "; Absolute int d#phi between outer and inner RecHit; Number of SimDoublets",
0704                      50,
0705                      0,
0706                      1000);
0707     hVector_pass_idphi_.at(layerPairIdIndex) = ibook.book1D("pass_idphi",
0708                                                             "idphi of RecHit pair " + layerTitle +
0709                                                                 " for SimDoublets passing all cuts; Absolute int d#phi "
0710                                                                 "between outer and inner RecHit; Number of SimDoublets",
0711                                                             50,
0712                                                             0,
0713                                                             1000);
0714
0715     // histogram for z window  (minz and maxz)
0716     hVector_innerZ_.at(layerPairIdIndex) =
0717         ibook.book1D("innerZ",
0718                      "z of the inner RecHit " + layerTitle + "; z of inner RecHit [cm]; Number of SimDoublets",
0719                      100,
0720                      -300,
0721                      300);
0722     hVector_pass_innerZ_.at(layerPairIdIndex) =
0723         ibook.book1D("pass_innerZ",
0724                      "z of the inner RecHit " + layerTitle +
0725                          " for SimDoublets passing all cuts; z of inner RecHit [cm]; Number of SimDoublets",
0726                      100,
0727                      -300,
0728                      300);
0729
0730     // histograms for cluster size and size differences
0731     hVector_DYsize_.at(layerPairIdIndex) =
0732         ibook.book1D("DYsize",
0733                      "Difference in cluster size along z between outer and inner RecHit " + layerTitle +
0734                          "; Absolute difference in cluster size along z of the two "
0735                          "RecHits [num of pixels]; Number of SimDoublets",
0736                      51,
0737                      -1,
0738                      50);
0739     hVector_DYPred_.at(layerPairIdIndex) =
0740         ibook.book1D("DYPred",
0741                      "Difference between actual and predicted cluster size along z of inner cluster " + layerTitle +
0742                          "; Absolute difference [num of pixels]; Number of SimDoublets",
0743                      51,
0744                      -1,
0745                      50);
0746     hVector_Ysize_.at(layerPairIdIndex) = ibook.book1D(
0747         "Ysize",
0748         "Cluster size along z " + layerTitle + "; Size along z of inner cluster [num of pixels]; Number of SimDoublets",
0749         51,
0750         -1,
0751         50);
0752   }
0753 }
0754
0755 // -------------------------------------------------------------------------------------------------------------
0756 // fillDescriptions
0757 // -------------------------------------------------------------------------------------------------------------
0758
0759 // dummy default fillDescription
0760 template <typename TrackerTraits>
0761 void SimDoubletsAnalyzer<TrackerTraits>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0762   edm::ParameterSetDescription desc;
0763   descriptions.addWithDefaultLabel(desc);
0764 }
0765
0766 // fillDescription for Phase 1
0767 template <>
0768 void SimDoubletsAnalyzer<pixelTopology::Phase1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0769   edm::ParameterSetDescription desc;
0770   simdoublets::fillDescriptionsCommon<pixelTopology::Phase1>(desc);
0771
0772   // input source for SimDoublets
0773   desc.add<edm::InputTag>("simDoubletsSrc", edm::InputTag("simDoubletsProducerPhase1"));
0774
0775   // layer pairs in reconstruction
0776   desc.add<std::vector<int>>(
0777           "layerPairs",
0778           std::vector<int>(std::begin(phase1PixelTopology::layerPairs), std::end(phase1PixelTopology::layerPairs)))
0779       ->setComment(
0780           "Array of length 2*NumberOfPairs where the elements at 2i and 2i+1 are the inner and outer layers of layer "
0781           "pair i");
0782
0783   // cutting parameters
0784   desc.add<int>("cellMinYSizeB1", 36)->setComment("Minimum cluster size for inner RecHit from B1");
0785   desc.add<int>("cellMinYSizeB2", 28)->setComment("Minimum cluster size for inner RecHit not from B1");
0786   desc.add<double>("cellZ0Cut", 12.0)->setComment("Maximum longitudinal impact parameter");
0787   desc.add<double>("cellPtCut", 0.5)->setComment("Minimum tranverse momentum");
0788   desc.add<std::vector<double>>(
0789           "cellMinz", std::vector<double>(std::begin(phase1PixelTopology::minz), std::end(phase1PixelTopology::minz)))
0790       ->setComment("Minimum z of inner RecHit for each layer pair");
0791   desc.add<std::vector<double>>(
0792           "cellMaxz", std::vector<double>(std::begin(phase1PixelTopology::maxz), std::end(phase1PixelTopology::maxz)))
0793       ->setComment("Maximum z of inner RecHit for each layer pair");
0794   desc.add<std::vector<int>>(
0795           "cellPhiCuts",
0796           std::vector<int>(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts)))
0797       ->setComment("Cuts in delta phi for cells");
0798   desc.add<std::vector<double>>(
0799           "cellMaxr", std::vector<double>(std::begin(phase1PixelTopology::maxr), std::end(phase1PixelTopology::maxr)))
0800       ->setComment("Cut for dr of cells");
0801
0802   descriptions.addWithDefaultLabel(desc);
0803 }
0804
0805 // fillDescription for Phase 2
0806 template <>
0807 void SimDoubletsAnalyzer<pixelTopology::Phase2>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0808   edm::ParameterSetDescription desc;
0809   simdoublets::fillDescriptionsCommon<pixelTopology::Phase2>(desc);
0810
0811   // input source for SimDoublets
0812   desc.add<edm::InputTag>("simDoubletsSrc", edm::InputTag("simDoubletsProducerPhase2"));
0813
0814   // layer pairs in reconstruction
0815   desc.add<std::vector<int>>(
0816           "layerPairs",
0817           std::vector<int>(std::begin(phase2PixelTopology::layerPairs), std::end(phase2PixelTopology::layerPairs)))
0818       ->setComment(
0819           "Array of length 2*NumberOfPairs where the elements at 2i and 2i+1 are the inner and outer layers of layer "
0820           "pair i");
0821
0822   // cutting parameters
0823   desc.add<int>("cellMinYSizeB1", 25)->setComment("Minimum cluster size for inner RecHit from B1");
0824   desc.add<int>("cellMinYSizeB2", 15)->setComment("Minimum cluster size for inner RecHit not from B1");
0825   desc.add<double>("cellZ0Cut", 7.5)->setComment("Maximum longitudinal impact parameter");
0826   desc.add<double>("cellPtCut", 0.85)->setComment("Minimum tranverse momentum");
0827   desc.add<std::vector<double>>(
0828           "cellMinz", std::vector<double>(std::begin(phase2PixelTopology::minz), std::end(phase2PixelTopology::minz)))
0829       ->setComment("Minimum z of inner RecHit for each layer pair");
0830   desc.add<std::vector<double>>(
0831           "cellMaxz", std::vector<double>(std::begin(phase2PixelTopology::maxz), std::end(phase2PixelTopology::maxz)))
0832       ->setComment("Maximum z of inner RecHit for each layer pair");
0833   desc.add<std::vector<int>>(
0834           "cellPhiCuts",
0835           std::vector<int>(std::begin(phase2PixelTopology::phicuts), std::end(phase2PixelTopology::phicuts)))
0836       ->setComment("Cuts in delta phi for cells");
0837   desc.add<std::vector<double>>(
0838           "cellMaxr", std::vector<double>(std::begin(phase2PixelTopology::maxr), std::end(phase2PixelTopology::maxr)))
0839       ->setComment("Cut for dr of cells");
0840
0841   descriptions.addWithDefaultLabel(desc);
0842 }
0843
0844 // define two plugins for Phase 1 and 2
0845 using SimDoubletsAnalyzerPhase1 = SimDoubletsAnalyzer<pixelTopology::Phase1>;
0846 using SimDoubletsAnalyzerPhase2 = SimDoubletsAnalyzer<pixelTopology::Phase2>;
0847
0848 // define this as a plug-in
0849 DEFINE_FWK_MODULE(SimDoubletsAnalyzerPhase1);
0850 DEFINE_FWK_MODULE(SimDoubletsAnalyzerPhase2);