SiPhase2Digitizer/test/PixelTestBeamValidation.cc

0001 //
0002 // Access Digi collection and creates histograms accumulating
0003 // data from the module to different pixel cells, 1x1 cell,
0004 // 2x2 cell, etc..
0005 //
0006 // Author:  J.Duarte-Campderros (IFCA)
0007 // Created: 2019-10-02
0008 //
0009
0010 // CMSSW Framework
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/Framework/interface/ESWatcher.h"
0013
0014 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0015
0016 // CMSSW Data formats
0017 #include "DataFormats/Math/interface/CMSUnits.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0020
0021 // system
0022 #include <algorithm>
0023
0024 // XXX - Be careful the relative position
0025 #include "PixelTestBeamValidation.h"
0026
0027 // Some needed units (um are more suited to the pixel size of the sensors)
0028 using cms_units::operators::operator""_deg;
0029 // Analogously to CMSUnits (no um defined, using inverse)
0030 constexpr double operator""_inv_um(long double length) { return length * 1e4; }
0031 // Energy (keV) -- to be used with the PSimHit::energyLoss with energy in GeV
0032 constexpr double operator""_inv_keV(long double energy_in_GeV) { return energy_in_GeV * 1e6; }
0033
0034 using Phase2TrackerGeomDetUnit = PixelGeomDetUnit;
0035
0036 PixelTestBeamValidation::PixelTestBeamValidation(const edm::ParameterSet& iConfig)
0037     : config_(iConfig),
0038       geomType_(iConfig.getParameter<std::string>("GeometryType")),
0039       //phiValues(iConfig.getParameter<std::vector<double> >("PhiAngles")),
0040       thresholdInElectrons_(iConfig.getParameter<double>("ThresholdInElectrons")),
0041       electronsPerADC_(iConfig.getParameter<double>("ElectronsPerADC")),
0042       tracksEntryAngleX_(
0043           iConfig.getUntrackedParameter<std::vector<double>>("TracksEntryAngleX", std::vector<double>())),
0044       tracksEntryAngleY_(
0045           iConfig.getUntrackedParameter<std::vector<double>>("TracksEntryAngleY", std::vector<double>())),
0046       digiToken_(consumes<edm::DetSetVector<PixelDigi>>(iConfig.getParameter<edm::InputTag>("PixelDigiSource"))),
0047       digiSimLinkToken_(
0048           consumes<edm::DetSetVector<PixelDigiSimLink>>(iConfig.getParameter<edm::InputTag>("PixelDigiSimSource"))),
0049       simTrackToken_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("SimTrackSource"))),
0050       topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0051       geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0052       topoBToken_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
0053       geomBToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()) {
0054   LogDebug("PixelTestBeamValidation") << ">>> Construct PixelTestBeamValidation ";
0055
0056   // The value to be used for ToT == 0 in electrons
0057   electronsAtToT0_ = 0.5 * (thresholdInElectrons_ + electronsPerADC_);
0058
0059   const std::vector<edm::InputTag> psimhit_v(config_.getParameter<std::vector<edm::InputTag>>("PSimHitSource"));
0060
0061   for (const auto& itag : psimhit_v) {
0062     simHitTokens_.push_back(consumes<edm::PSimHitContainer>(itag));
0063   }
0064
0065   // Parse the entry angles parameter. Remember the angles are defined in
0066   // radians and defined as 0 when perpendicular to the detector plane
0067   // ---------------------------------------------------------------------
0068   // Helper map to build up the active entry angles
0069   std::map<std::string, std::vector<double>*> prov_ref_m;
0070   // Get the range of entry angles for the tracks on the detector surfaces, if any
0071   if (tracksEntryAngleX_.size() != 0) {
0072     prov_ref_m["X"] = &tracksEntryAngleX_;
0073   }
0074   if (tracksEntryAngleY_.size() != 0) {
0075     prov_ref_m["Y"] = &tracksEntryAngleY_;
0076   }
0077
0078   // translation from string to int
0079   std::map<std::string, unsigned int> conversor({{"X", 0}, {"Y", 1}});
0080
0081   // For each range vector do some consistency checks
0082   for (const auto& label_v : prov_ref_m) {
0083     // If provided 2
0084     if (label_v.second->size() == 2) {
0085       // Just order the ranges, lower index the minimum
0086       std::sort(label_v.second->begin(), label_v.second->end());
0087     } else if (label_v.second->size() == 1) {
0088       // Create the range with  a +- 1 deg
0089       label_v.second->push_back((*label_v.second)[0] + 1.0_deg);
0090       (*label_v.second)[0] -= 1.0_deg;
0091     } else {
0092       // Not valid,
0093       throw cms::Exception("Configuration") << "Setup TrackEntryAngle parameters"
0094                                             << "Invalid number of elements in 'TracksEntryAngle" << label_v.first
0095                                             << ".size()' = " << label_v.second->size() << ". Valid sizes are 1 or 2.";
0096     }
0097     // convert the angles into its tangent (the check is going to be done
0098     // against the tangent dx/dz and dy/dz)
0099     (*label_v.second)[0] = std::tan((*label_v.second)[0]);
0100     (*label_v.second)[1] = std::tan((*label_v.second)[1]);
0101     active_entry_angles_[conversor[label_v.first]] =
0102         std::pair<double, double>({(*label_v.second)[0], (*label_v.second)[1]});
0103   }
0104
0105   if (prov_ref_m.size() != 0) {
0106     // The algorithm is defined in the implementation of _check_input_angles_
0107     use_this_track_ = std::bind(&PixelTestBeamValidation::_check_input_angles_, this, std::placeholders::_1);
0108     edm::LogInfo("Configuration") << "Considering hits from tracks entering the detectors between\n "
0109                                   << "Considering hits from tracks entering the detectors between\n "
0110                                   << "\tX-plane: (" << std::atan(active_entry_angles_[0].first) << ","
0111                                   << std::atan(active_entry_angles_[0].second) << ") rad. "
0112                                   << "\tY-plane: (" << std::atan(active_entry_angles_[1].first) << ","
0113                                   << std::atan(active_entry_angles_[1].second) << ") rad. ";
0114   } else {
0115     // There is no requiremnt, so always process it
0116     use_this_track_ = [](const PSimHit*) -> bool { return true; };
0117   }
0118 }
0119
0120 //
0121 // destructor
0122 //
0123 PixelTestBeamValidation::~PixelTestBeamValidation() {
0124   LogDebug("PixelTestBeamValidation") << ">>> Destroy PixelTestBeamValidation ";
0125 }
0126
0127 // -- DQM Begin Run
0128 //
0129 void PixelTestBeamValidation::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0130   edm::LogInfo("PixelTestBeamValidation") << "Initialize PixelTestBeamValidation ";
0131 }
0132 //
0133 // -- Analyze
0134 //
0135 void PixelTestBeamValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0136   // First clear the memoizers
0137   m_tId_det_simhits_.clear();
0138   m_illuminated_pixels_.clear();
0139
0140   // Get digis and the simhit links and the simhits
0141   edm::Handle<edm::DetSetVector<PixelDigiSimLink>> digiSimLinkHandle;
0142   iEvent.getByToken(digiSimLinkToken_, digiSimLinkHandle);
0143   const edm::DetSetVector<PixelDigiSimLink>* simdigis = digiSimLinkHandle.product();
0144
0145   edm::Handle<edm::DetSetVector<PixelDigi>> digiHandle;
0146   iEvent.getByToken(digiToken_, digiHandle);
0147   const edm::DetSetVector<PixelDigi>* digis = digiHandle.product();
0148
0149   // Vector of simHits
0150   // XXX : NOt like that, just an example
0151   std::vector<const edm::PSimHitContainer*> simhits;
0152   simhits.reserve(simHitTokens_.size());
0153   for (const auto& sh_token : simHitTokens_) {
0154     edm::Handle<edm::PSimHitContainer> simHitHandle;
0155     iEvent.getByToken(sh_token, simHitHandle);
0156     if (!simHitHandle.isValid()) {
0157       continue;
0158     }
0159     simhits.push_back(simHitHandle.product());
0160   }
0161
0162   // Get SimTrack
0163   /*edm::Handle<edm::SimTrackContainer> simTrackHandle;
0164     iEvent.getByToken(simTrackToken_, simTrackHandle);
0165     const edm::SimTrackContainer * simtracks = simTrackHandle.product();*/
0166
0167   // Geometry description
0168   edm::ESWatcher<TrackerDigiGeometryRecord> tkDigiGeomWatcher;
0169   if (!tkDigiGeomWatcher.check(iSetup)) {
0170     // XXX -- Should raise a Warning??
0171     return;
0172   }
0173
0174   // Tracker geometry and topology
0175   const TrackerTopology* topo = &iSetup.getData(topoToken_);
0176   const TrackerGeometry* tkGeom = &iSetup.getData(geomToken_);
0177
0178   // Let's loop over the detectors
0179   for (auto const& dunit : tkGeom->detUnits()) {
0180     if (!isPixelSystem_(dunit)) {
0181       continue;
0182     }
0183
0184     // --------------------------------------------
0185     // Get more info about the detector unit
0186     // -- layer, topology (nrows,ncolumns,pitch...)
0187     const Phase2TrackerGeomDetUnit* tkDetUnit = dynamic_cast<const Phase2TrackerGeomDetUnit*>(dunit);
0188     //const PixelTopology * topo = tkDetUnit->specificTopology();
0189
0190     const auto detId = dunit->geographicalId();
0191     const int layer = topo->layer(detId);
0192     // Get the relevant histo key
0193     const auto& me_unit = meUnit_(tkDetUnit->type().isBarrel(), layer, topo->side(detId));
0194
0195     // Get the id of the detector unit
0196     const unsigned int detId_raw = detId.rawId();
0197
0198     // Loop over the simhits to obtain the list of PSimHits created
0199     // in this detector unit
0200     std::vector<const PSimHit*> it_simhits;
0201
0202     for (const auto* sh_c : simhits) {
0203       for (const auto& sh : *sh_c) {
0204         if (sh.detUnitId() == detId_raw) {
0205           it_simhits.push_back(&sh);  //insert and/or reserve (?)
0206         }
0207       }
0208     }
0209
0210     if (it_simhits.size() == 0) {
0211       continue;
0212     }
0213
0214     // Find RAW digis (digis) created in this det unit
0215     const auto& it_digis = digis->find(detId);
0216
0217     //std::cout << "DETECTOR: " << tkDetUnit->type().name() << " ME UNIT: " << me_unit << std::endl;
0218
0219     // Loop over the list of PSimHits (i.e. the deposits created
0220     // by the primary+secundaries) to check if they are associated with
0221     // some digi, that is, if the simdigi link exists and to obtain the list
0222     // of channels illuminated
0223     for (const auto* psh : it_simhits) {
0224       // Check user conditions to accept the hits
0225       if (!use_this_track_(psh)) {
0226         continue;
0227       }
0228       // Fill some sim histograms
0229       const GlobalPoint tk_ep_gbl(dunit->surface().toGlobal(psh->entryPoint()));
0230       vME_track_XYMap_->Fill(tk_ep_gbl.x(), tk_ep_gbl.y());
0231       vME_track_RZMap_->Fill(tk_ep_gbl.z(), std::hypot(tk_ep_gbl.x(), tk_ep_gbl.y()));
0232       vME_track_dxdzAngle_[me_unit]->Fill(std::atan2(psh->momentumAtEntry().x(), psh->momentumAtEntry().z()));
0233       vME_track_dydzAngle_[me_unit]->Fill(std::atan2(psh->momentumAtEntry().y(), psh->momentumAtEntry().z()));
0234
0235       // Obtain the detected position of the sim particle:
0236       // the middle point between the entry and the exit
0237       const auto psh_pos = tkDetUnit->specificTopology().measurementPosition(psh->localPosition());
0238
0239       // MC Cluster finding: get the channels illuminated during digitization by this PSimHit
0240       // based on MC truth info (simdigi links)
0241       const auto psh_channels = get_illuminated_channels_(*psh, detId, simdigis);
0242
0243       // Get the total charge for this cluster size
0244       // and obtain the center of the cluster using a charge-weighted mean
0245       int cluster_tot = 0;
0246       double cluster_tot_elec = 0.0;
0247       int cluster_size = 0;
0248       std::pair<std::set<int>, std::set<int>> cluster_size_xy;
0249       std::pair<double, double> cluster_position({0.0, 0.0});
0250       std::set<int> used_channel;
0251       for (const auto& ch : psh_channels) {
0252         // Not re-using the digi XXX
0253         if (used_channel.find(ch) != used_channel.end()) {
0254           continue;
0255         }
0256         const PixelDigi& current_digi = get_digi_from_channel_(ch, it_digis);
0257         used_channel.insert(ch);
0258         // Fill the digi histograms
0259         vME_digi_charge1D_[me_unit]->Fill(current_digi.adc());
0260         // Fill maps: get the position in the sensor local frame to convert into global
0261         const LocalPoint digi_local_pos(
0262             tkDetUnit->specificTopology().localPosition(MeasurementPoint(current_digi.row(), current_digi.column())));
0263         const GlobalPoint digi_global_pos(dunit->surface().toGlobal(digi_local_pos));
0264         vME_digi_XYMap_->Fill(digi_global_pos.x(), digi_global_pos.y());
0265         vME_digi_RZMap_->Fill(digi_global_pos.z(), std::hypot(digi_global_pos.x(), digi_global_pos.y()));
0266         // Create the MC-cluster
0267         cluster_tot += current_digi.adc();
0268         // Assign the middle value between the threshold and the ElectronPerADC parameter
0269         // to the first bin in order to allow ToT = 0 (valid value)
0270         const double pixel_charge_elec =
0271             (bool(current_digi.adc()) ? (current_digi.adc() * electronsPerADC_) : electronsAtToT0_);
0272         cluster_tot_elec += pixel_charge_elec;
0273         // Use the center of the pixel
0274         cluster_position.first += pixel_charge_elec * (current_digi.row() + 0.5);
0275         cluster_position.second += pixel_charge_elec * (current_digi.column() + 0.5);
0276         // Size
0277         cluster_size_xy.first.insert(current_digi.row());
0278         cluster_size_xy.second.insert(current_digi.column());
0279         ++cluster_size;
0280       }
0281
0282       // Be careful here, there is 1 entry per each simhit
0283       vME_clsize1D_[me_unit]->Fill(cluster_size);
0284       vME_clsize1Dx_[me_unit]->Fill(cluster_size_xy.first.size());
0285       vME_clsize1Dy_[me_unit]->Fill(cluster_size_xy.second.size());
0286
0287       // mean weighted
0288       cluster_position.first /= cluster_tot_elec;
0289       cluster_position.second /= cluster_tot_elec;
0290
0291       // -- XXX Be careful, secondaries with already used the digis
0292       //        are going the be lost (then lost on efficiency)
0293       // Efficiency --> It was found a cluster of digis?
0294       const bool is_cluster_present = (cluster_size > 0);
0295
0296       // Get topology info of the module sensor
0297       //-const int n_rows = tkDetUnit->specificTopology().nrows();
0298       //-const int n_cols = tkDetUnit->specificTopology().ncolumns();
0299       const auto pitch = tkDetUnit->specificTopology().pitch();
0300       // Residuals, convert them to longitud units (so far, in units of row, col)
0301       const double dx_um = (psh_pos.x() - cluster_position.first) * pitch.first * 1.0_inv_um;
0302       const double dy_um = (psh_pos.y() - cluster_position.second) * pitch.second * 1.0_inv_um;
0303       if (is_cluster_present) {
0304         vME_charge1D_[me_unit]->Fill(cluster_tot);
0305         vME_charge_elec1D_[me_unit]->Fill(cluster_tot_elec);
0306         vME_dx1D_[me_unit]->Fill(dx_um);
0307         vME_dy1D_[me_unit]->Fill(dy_um);
0308         vME_dxy2D_[me_unit]->Fill(dx_um, dy_um);
0309         // The track energy loss corresponding to that cluster
0310         vME_sim_cluster_charge_[me_unit]->Fill(psh->energyLoss() * 1.0_inv_keV, cluster_tot_elec);
0311       }
0312       // Histograms per cell
0313       for (unsigned int i = 0; i < vME_position_cell_[me_unit].size(); ++i) {
0314         // Convert the PSimHit center position to the IxI-cell
0315         const std::pair<double, double> icell_psh = pixel_cell_transformation_(psh_pos, i, pitch);
0316         // Efficiency: (PSimHit matched to a digi-cluster)/PSimHit
0317         vME_eff_cell_[me_unit][i]->Fill(
0318             icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um, is_cluster_present);
0319         vME_pshpos_cell_[me_unit][i]->Fill(icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um);
0320         // Digi clusters related histoos
0321         if (is_cluster_present) {
0322           // Convert to the i-cell
0323           //const std::pair<double,double> icell_digi_cluster   = pixel_cell_transformation_(cluster_position,i,pitch);
0324           // Position
0325           vME_position_cell_[me_unit][i]->Fill(icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um);
0326           // Residuals
0327           vME_dx_cell_[me_unit][i]->Fill(icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um, dx_um);
0328           vME_dy_cell_[me_unit][i]->Fill(icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um, dy_um);
0329           // Charge
0330           vME_charge_cell_[me_unit][i]->Fill(icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um, cluster_tot);
0331           vME_charge_elec_cell_[me_unit][i]->Fill(
0332               icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um, cluster_tot_elec);
0333           // Cluster size
0334           vME_clsize_cell_[me_unit][i]->Fill(icell_psh.first * 1.0_inv_um, icell_psh.second * 1.0_inv_um, cluster_size);
0335         }
0336       }
0337     }
0338   }
0339 }
0340
0341 //
0342 // -- Book Histograms
0343 //
0344 void PixelTestBeamValidation::bookHistograms(DQMStore::IBooker& ibooker,
0345                                              edm::Run const& iRun,
0346                                              edm::EventSetup const& iSetup) {
0347   // Get Geometry to associate a folder to each Pixel subdetector
0348   const TrackerGeometry* tkGeom = &iSetup.getData(geomBToken_);
0349
0350   // Tracker Topology (layers, modules, side, etc..)
0351   const TrackerTopology* topo = &iSetup.getData(topoBToken_);
0352
0353   const std::string top_folder = config_.getParameter<std::string>("TopFolderName");
0354
0355   // Histograms independent of the subdetector units
0356   ibooker.cd();
0357   ibooker.setCurrentFolder(top_folder);
0358   vME_track_XYMap_ =
0359       setupH2D_(ibooker, "TrackXY", "Track entering position in the tracker system;x [cm];y [cm];N_{tracksId}");
0360   vME_track_RZMap_ =
0361       setupH2D_(ibooker, "TrackRZ", "Track entering position in the tracker system;z [cm];r [cm];N_{tracksId}");
0362   vME_digi_XYMap_ = setupH2D_(ibooker, "DigiXY", "Digi position;x [cm];y [cm];N_{digi}");
0363   vME_digi_RZMap_ = setupH2D_(ibooker, "DigiRZ", "Digi position;z [cm];r [cm];N_{digi}");
0364
0365   // Get all pixel subdetector, create histogram by layers
0366   // -- More granularity can be accomplished by modules (1 central + 5 modules: in z,
0367   //    each one composed by 12 pixel modules (1 long pixel + 2 ROCs) in Phi = 108 pixel-modules
0368   for (auto const& dunit : tkGeom->detUnits()) {
0369     if (!isPixelSystem_(dunit)) {
0370       continue;
0371     }
0372     //if( ! dtype.isInnerTracker() || ! dtype.isTrackerPixel() )
0373     //{
0374     //    continue;
0375     //}
0376     const auto& dtype = dunit->type();
0377
0378     const auto detId = dunit->geographicalId();
0379     const int layer = topo->layer(detId);
0380     // Create the sub-detector histogram if needed
0381     // Barrel_Layer/Endcap_Layer_Side
0382     const int me_unit = meUnit_(dtype.isBarrel(), layer, topo->side(detId));
0383     if (vME_position_cell_.find(me_unit) == vME_position_cell_.end()) {
0384       std::string folder_name(top_folder + "/");
0385       if (dtype.isBarrel()) {
0386         folder_name += "Barrel";
0387       } else if (dtype.isEndcap()) {
0388         folder_name += "Endcap";
0389       } else {
0390         cms::Exception("Geometry") << "Tracker subdetector '" << dtype.subDetector()
0391                                    << "' malformed: does not belong to the Endcap neither the Barrel";
0392       }
0393       folder_name += "/Layer" + std::to_string(layer);
0394
0395       if (topo->side(detId)) {
0396         folder_name += "/Side" + std::to_string(topo->side(detId));
0397       }
0398       // Go the proper folder
0399       ibooker.cd();
0400       ibooker.setCurrentFolder(folder_name);
0401
0402       const TrackerGeomDet* geomDetUnit = tkGeom->idToDetUnit(detId);
0403       const Phase2TrackerGeomDetUnit* tkDetUnit = dynamic_cast<const Phase2TrackerGeomDetUnit*>(geomDetUnit);
0404       const int nrows = tkDetUnit->specificTopology().nrows();
0405       const int ncols = tkDetUnit->specificTopology().ncolumns();
0406       const auto pitch = tkDetUnit->specificTopology().pitch();
0407       //const double x_size = nrows*pitch.first;
0408       //const double y_size = ncols*pitch.second;
0409
0410       // And create the histos
0411       // Per detector unit histos
0412       vME_clsize1D_[me_unit] =
0413           setupH1D_(ibooker, "ClusterSize1D", "MC-truth DIGI cluster size;Cluster size;N_{clusters}");
0414       vME_clsize1Dx_[me_unit] =
0415           setupH1D_(ibooker, "ClusterSize1Dx", "MC-truth DIGI cluster size in X;Cluster size;N_{clusters}");
0416       vME_clsize1Dy_[me_unit] =
0417           setupH1D_(ibooker, "ClusterSize1Dy", "MC-truth DIGI cluster size in Y;Cluster size;N_{clusters}");
0418       vME_charge1D_[me_unit] =
0419           setupH1D_(ibooker, "Charge1D", "MC-truth DIGI cluster charge;Cluster charge [ToT];N_{clusters}");
0420       vME_charge_elec1D_[me_unit] =
0421           setupH1D_(ibooker, "ChargeElec1D", "MC-truth DIGI cluster charge;Cluster charge [Electrons];N_{clusters}");
0422       vME_track_dxdzAngle_[me_unit] = setupH1D_(
0423           ibooker,
0424           "TrackAngleDxdz",
0425           "Angle between the track-momentum and detector surface (X-plane);#pi/2-#theta_{x} [rad];N_{tracks}");
0426       vME_track_dydzAngle_[me_unit] = setupH1D_(
0427           ibooker,
0428           "TrackAngleDydz",
0429           "Angle between the track-momentum and detector surface (Y-plane);#pi/2-#theta_{y} [rad];N_{tracks}");
0430       vME_dx1D_[me_unit] = setupH1D_(
0431           ibooker, "Dx1D", "MC-truth DIGI cluster residuals X;x_{PSimHit}-x^{cluster}_{digi} [#mum];N_{digi clusters}");
0432       vME_dy1D_[me_unit] = setupH1D_(
0433           ibooker, "Dy1D", "MC-truth DIGI cluster residual Y;y_{PSimHit}-y^{cluster}_{digi} [#mum];N_{digi clusters}");
0434       vME_dxy2D_[me_unit] = setupH2D_(ibooker,
0435                                       "Dxy2D",
0436                                       "MC-truth DIGI cluster residuals;x_{PSimHit}-x^{cluster}_{digi} "
0437                                       "[#mum];y_{PSimHit}-y^{cluster}_{digi} [#mum];N_{digi clusters}");
0438       vME_digi_charge1D_[me_unit] = setupH1D_(ibooker, "DigiCharge1D", "Digi charge;digi charge [ToT];N_{digi}");
0439       vME_sim_cluster_charge_[me_unit] =
0440           setupH2D_(ibooker,
0441                     "SimClusterCharge",
0442                     "PSimHit Energy deposit vs. Cluster Charge;deposited E_{sim} [keV];cluster_{charge} [Electrons];");
0443
0444       // The histos per cell
0445       // Prepare the ranges: 0- whole sensor, 1- cell 1x1, 2-cell 2x2,
0446       const std::vector<std::pair<double, double>> xranges = {
0447           std::make_pair<double, double>(0, (nrows - 1) * pitch.first * 1.0_inv_um),
0448           std::make_pair<double, double>(0, pitch.first * 1.0_inv_um),
0449           std::make_pair<double, double>(0, 2.0 * pitch.first * 1.0_inv_um)};
0450       const std::vector<std::pair<double, double>> yranges = {
0451           std::make_pair<double, double>(0, (ncols - 1) * pitch.second * 1.0_inv_um),
0452           std::make_pair<double, double>(0, pitch.second * 1.0_inv_um),
0453           std::make_pair<double, double>(0, 2.0 * pitch.second * 1.0_inv_um)};
0454       for (unsigned int i = 0; i < xranges.size(); ++i) {
0455         const std::string cell("Cell " + std::to_string(i) + "x" + std::to_string(i) + ": ");
0456         vME_pshpos_cell_[me_unit].push_back(
0457             setupH2D_(ibooker,
0458                       "Position_" + std::to_string(i),
0459                       cell + "PSimHit middle point position;x [#mum];y [#mum];N_{clusters}",
0460                       xranges[i],
0461                       yranges[i]));
0462         vME_position_cell_[me_unit].push_back(
0463             setupH2D_(ibooker,
0464                       "MatchedPosition_" + std::to_string(i),
0465                       cell + "PSimHit matched to a DIGI cluster position ;x [#mum];y [#mum];N_{clusters}",
0466                       xranges[i],
0467                       yranges[i]));
0468         vME_eff_cell_[me_unit].push_back(setupProf2D_(ibooker,
0469                                                       "Efficiency_" + std::to_string(i),
0470                                                       cell + "MC-truth efficiency;x [#mum];y [#mum];<#varepsilon>",
0471                                                       xranges[i],
0472                                                       yranges[i]));
0473         vME_clsize_cell_[me_unit].push_back(
0474             setupProf2D_(ibooker,
0475                          "ClusterSize_" + std::to_string(i),
0476                          cell + "MC-truth cluster size;x [#mum];y [#mum];<Cluster size>",
0477                          xranges[i],
0478                          yranges[i]));
0479         vME_charge_cell_[me_unit].push_back(setupProf2D_(ibooker,
0480                                                          "Charge_" + std::to_string(i),
0481                                                          cell + "MC-truth charge;x [#mum];y [#mum];<ToT>",
0482                                                          xranges[i],
0483                                                          yranges[i]));
0484         vME_charge_elec_cell_[me_unit].push_back(setupProf2D_(ibooker,
0485                                                               "Charge_elec_" + std::to_string(i),
0486                                                               cell + "MC-truth charge;x [#mum];y [#mum];<Electrons>",
0487                                                               xranges[i],
0488                                                               yranges[i]));
0489         vME_dx_cell_[me_unit].push_back(setupProf2D_(ibooker,
0490                                                      "Dx_" + std::to_string(i),
0491                                                      cell + "MC-truth residuals;x [#mum];y [#mum];<#Deltax [#mum]>",
0492                                                      xranges[i],
0493                                                      yranges[i]));
0494         vME_dy_cell_[me_unit].push_back(setupProf2D_(ibooker,
0495                                                      "Dy_" + std::to_string(i),
0496                                                      cell + "MC-truth residuals;x [#mum];y [#mum];<#Deltay [#mum]>",
0497                                                      xranges[i],
0498                                                      yranges[i]));
0499       }
0500
0501       edm::LogInfo("PixelTestBeamValidation") << "Booking Histograms in: " << folder_name << " ME UNIT:" << me_unit;
0502       edm::LogInfo("PixelTestBeamValidation") << "Booking Histograms in: " << folder_name;
0503     }
0504   }
0505 }
0506
0507 bool PixelTestBeamValidation::isPixelSystem_(const GeomDetUnit* dunit) const {
0508   const auto& dtype = dunit->type();
0509   if (!dtype.isInnerTracker() || !dtype.isTrackerPixel()) {
0510     return false;
0511   }
0512   return true;
0513 }
0514
0515 // Helper functions to setup
0516 int PixelTestBeamValidation::meUnit_(bool isBarrel, int layer, int side) const {
0517   // [isBarrel][layer][side]
0518   // [X][XXX][XX]
0519   return (static_cast<int>(isBarrel) << 6) | (layer << 3) | side;
0520 }
0521
0522 PixelTestBeamValidation::MonitorElement* PixelTestBeamValidation::setupH1D_(DQMStore::IBooker& ibooker,
0523                                                                             const std::string& histoname,
0524                                                                             const std::string& title) {
0525   // Config need to have exactly the same histo name
0526   edm::ParameterSet params = config_.getParameter<edm::ParameterSet>(histoname);
0527   return ibooker.book1D(histoname,
0528                         title,
0529                         params.getParameter<int32_t>("Nxbins"),
0530                         params.getParameter<double>("xmin"),
0531                         params.getParameter<double>("xmax"));
0532 }
0533
0534 PixelTestBeamValidation::MonitorElement* PixelTestBeamValidation::setupH2D_(DQMStore::IBooker& ibooker,
0535                                                                             const std::string& histoname,
0536                                                                             const std::string& title) {
0537   // Config need to have exactly the same histo name
0538   edm::ParameterSet params = config_.getParameter<edm::ParameterSet>(histoname);
0539   return ibooker.book2D(histoname,
0540                         title,
0541                         params.getParameter<int32_t>("Nxbins"),
0542                         params.getParameter<double>("xmin"),
0543                         params.getParameter<double>("xmax"),
0544                         params.getParameter<int32_t>("Nybins"),
0545                         params.getParameter<double>("ymin"),
0546                         params.getParameter<double>("ymax"));
0547 }
0548
0549 PixelTestBeamValidation::MonitorElement* PixelTestBeamValidation::setupH2D_(DQMStore::IBooker& ibooker,
0550                                                                             const std::string& histoname,
0551                                                                             const std::string& title,
0552                                                                             const std::pair<double, double>& xranges,
0553                                                                             const std::pair<double, double>& yranges) {
0554   // Config need to have exactly the same histo name
0555   edm::ParameterSet params = config_.getParameter<edm::ParameterSet>(histoname);
0556   return ibooker.book2D(histoname,
0557                         title,
0558                         params.getParameter<int32_t>("Nxbins"),
0559                         xranges.first,
0560                         xranges.second,
0561                         params.getParameter<int32_t>("Nybins"),
0562                         yranges.first,
0563                         yranges.second);
0564 }
0565
0566 PixelTestBeamValidation::MonitorElement* PixelTestBeamValidation::setupProf2D_(
0567     DQMStore::IBooker& ibooker,
0568     const std::string& histoname,
0569     const std::string& title,
0570     const std::pair<double, double>& xranges,
0571     const std::pair<double, double>& yranges) {
0572   // Config need to have exactly the same histo name
0573   edm::ParameterSet params = config_.getParameter<edm::ParameterSet>(histoname);
0574   return ibooker.bookProfile2D(histoname,
0575                                title,
0576                                params.getParameter<int32_t>("Nxbins"),
0577                                xranges.first,
0578                                xranges.second,
0579                                params.getParameter<int32_t>("Nybins"),
0580                                yranges.first,
0581                                yranges.second,
0582                                params.getParameter<double>("zmin"),
0583                                params.getParameter<double>("zmax"));
0584 }
0585
0586 const PixelDigi& PixelTestBeamValidation::get_digi_from_channel_(
0587     int ch, const edm::DetSetVector<PixelDigi>::const_iterator& itdigis) {
0588   for (const auto& dh : *itdigis) {
0589     if (dh.channel() == ch) {
0590       return dh;
0591     }
0592   }
0593   // MAke sense not to find a digi?
0594   throw cms::Exception("DIGI Pixel Validation") << "Not found a PixelDig for the given channel: " << ch;
0595 }
0596
0597 const SimTrack* PixelTestBeamValidation::get_simtrack_from_id_(unsigned int idx, const edm::SimTrackContainer* stc) {
0598   for (const auto& st : *stc) {
0599     if (st.trackId() == idx) {
0600       return &st;
0601     }
0602   }
0603   // Any simtrack correspond to this trackid index
0604   //edm::LogWarning("PixelTestBeamValidation::get_simtrack_from_id_")
0605   edm::LogInfo("PixelTestBeamValidation::get_simtrack_from_id_") << "Not found any SimTrack with trackId: " << idx;
0606   return nullptr;
0607 }
0608
0609 const std::vector<const PSimHit*> PixelTestBeamValidation::get_simhits_from_trackid_(
0610     unsigned int tid, unsigned int detid_raw, const std::vector<const edm::PSimHitContainer*>& psimhits) {
0611   // It was already found?
0612   if (m_tId_det_simhits_.find(tid) != m_tId_det_simhits_.end()) {
0613     // Note that if there were no simhits in detid_raw, now it
0614     // is creating an empty std::vector<const edm::PSimHitContainer*>
0615     return m_tId_det_simhits_[tid][detid_raw];
0616   }
0617
0618   // Otherwise,
0619   // Create the new map for the track
0620   m_tId_det_simhits_[tid] = std::map<unsigned int, std::vector<const PSimHit*>>();
0621   // and search for it, all the PsimHit found in all detectors
0622   // are going to be seeked once, therefore memoizing already
0623   for (const auto* sh_c : psimhits) {
0624     for (const auto& sh : *sh_c) {
0625       if (sh.trackId() == tid) {
0626         m_tId_det_simhits_[tid][sh.detUnitId()].push_back(&sh);
0627       }
0628     }
0629   }
0630   // And returning what was requested
0631   return m_tId_det_simhits_[tid][detid_raw];
0632 }
0633
0634 const std::pair<double, double> PixelTestBeamValidation::pixel_cell_transformation_(
0635     const MeasurementPoint& pos, unsigned int icell, const std::pair<double, double>& pitch) {
0636   return pixel_cell_transformation_(std::pair<double, double>({pos.x(), pos.y()}), icell, pitch);
0637 }
0638
0639 const std::pair<double, double> PixelTestBeamValidation::pixel_cell_transformation_(
0640     const std::pair<double, double>& pos, unsigned int icell, const std::pair<double, double>& pitch) {
0641   // Get the position modulo icell
0642   // Case the whole detector (icell==0)
0643   double xmod = pos.first;
0644   double ymod = pos.second;
0645
0646   // Actual pixel cells
0647   if (icell != 0) {
0648     xmod = std::fmod(pos.first, icell);
0649     ymod = std::fmod(pos.second, icell);
0650   }
0651
0652   const double xcell = xmod * pitch.first;
0653   const double ycell = ymod * pitch.second;
0654
0655   return std::pair<double, double>({xcell, ycell});
0656 }
0657
0658 //bool PixelTestBeamValidation::channel_iluminated_by_(const MeasurementPoint & localpos,int channel, double tolerance) const
0659 //{
0660 //    const auto pos_channel(PixelDigi::channelToPixel(channel));
0661 //    if( std::fabs(localpos.x()-pos_channel.first) <= tolerance
0662 //            && std::fabs(localpos.y()-pos_channel.second) <= tolerance )
0663 //    {
0664 //        return true;
0665 //    }
0666 //    return false;
0667 //}
0668
0669 bool PixelTestBeamValidation::channel_iluminated_by_(const PSimHit& ps,
0670                                                      int channel,
0671                                                      const PixelGeomDetUnit* tkDetUnit) {
0672   // Get the list of pixels illuminated by the PSimHit
0673   const auto pixel_list = get_illuminated_pixels_(ps, tkDetUnit);
0674   // Get the digi position
0675   const auto pos_channel(PixelDigi::channelToPixel(channel));
0676
0677   for (const auto& px_py : pixel_list) {
0678     if (px_py.first == pos_channel.first && px_py.second == pos_channel.second) {
0679       return true;
0680     }
0681   }
0682   return false;
0683 }
0684
0685 std::set<int> PixelTestBeamValidation::get_illuminated_channels_(const PSimHit& ps,
0686                                                                  const DetId& detid,
0687                                                                  const edm::DetSetVector<PixelDigiSimLink>* simdigis) {
0688   // Find simulated digi links (simdigis) created in this det unit
0689   const auto& it_simdigilink = simdigis->find(detid);
0690
0691   if (it_simdigilink == simdigis->end()) {
0692     return std::set<int>();
0693   }
0694
0695   std::set<int> channels;
0696   for (const auto& hdsim : *it_simdigilink) {
0697     if (ps.trackId() == hdsim.SimTrackId()) {
0698       channels.insert(hdsim.channel());
0699     }
0700   }
0701   return channels;
0702 }
0703
0704 std::set<std::pair<int, int>> PixelTestBeamValidation::get_illuminated_pixels_(const PSimHit& ps,
0705                                                                                const PixelGeomDetUnit* tkDetUnit) {
0706   auto ps_key = reinterpret_cast<std::uintptr_t>(&ps);
0707
0708   //  -- Check if was already memoized
0709   if (m_illuminated_pixels_.find(ps_key) != m_illuminated_pixels_.end()) {
0710     return m_illuminated_pixels_[ps_key];
0711   }
0712
0713   // Get the entry point - exit point position
0714   const double min_x = std::min(ps.entryPoint().x(), ps.exitPoint().x());
0715   const double min_y = std::min(ps.entryPoint().y(), ps.exitPoint().y());
0716   const double max_x = std::max(ps.entryPoint().x(), ps.exitPoint().x());
0717   const double max_y = std::max(ps.entryPoint().y(), ps.exitPoint().y());
0718   // Get the position in readout units for each point
0719   const auto min_pos_pre = tkDetUnit->specificTopology().measurementPosition(LocalPoint(min_x, min_y));
0720   const auto max_pos_pre = tkDetUnit->specificTopology().measurementPosition(LocalPoint(max_x, max_y));
0721   // Adding a unit at each side in order to include charge migration
0722   const MeasurementPoint min_pos(std::max(0.0, min_pos_pre.x() - 1.0), std::max(0.0, min_pos_pre.y() - 1.0));
0723   const MeasurementPoint max_pos(std::min(max_pos_pre.x() + 1.0, tkDetUnit->specificTopology().nrows() - 1.0),
0724                                  std::min(max_pos_pre.y() + 1.0, tkDetUnit->specificTopology().ncolumns() - 1.0));
0725   // Count how many cells has passed. Use the most conservative rounding:
0726   // round for maximums and int (floor) for minimums
0727   //const int ncells_x = std::round(max_pos.x())-std::floor(min_pos.x());
0728   //const int ncells_y = std::round(max_pos.y())-std::floor(min_pos.y());
0729
0730   //std::cout << "ENTRANDO --- PSimHit (" << ps_key << ") entry point: " << ps.entryPoint() << " exitPoint: " << ps.exitPoint()
0731   //  << " Entry (pixel units): " << min_pos
0732   //  << " Exit  (pixel units): " << max_pos << ")"
0733   //  << " ILLUMINATED PIXELS: ";
0734
0735   std::set<std::pair<int, int>> illuminated_pixels;
0736   for (unsigned int px = std::floor(min_pos.x()); px <= std::round(max_pos.x()); ++px) {
0737     for (unsigned int py = std::floor(min_pos.y()); py <= std::round(max_pos.y()); ++py) {
0738       //std::cout << " [" << px << "," << py << "]";
0739       illuminated_pixels.insert(std::pair<int, int>(px, py));
0740     }
0741   }
0742   //std::cout << " TOTAL PIXELS: " << illuminated_pixels.size() << std::endl;
0743   // Memoize, and return what expected.
0744   m_illuminated_pixels_[ps_key] = illuminated_pixels;
0745   return m_illuminated_pixels_[ps_key];
0746 }
0747
0748 bool PixelTestBeamValidation::_check_input_angles_(const PSimHit* psimhit) {
0749   // Create a vector to check against the range map where
0750   // X axis is in the key 0, Y axis in the key 1
0751   const std::vector<double> entry_tan({psimhit->momentumAtEntry().x() / psimhit->momentumAtEntry().z(),
0752                                        psimhit->momentumAtEntry().y() / psimhit->momentumAtEntry().z()});
0753
0754   for (const auto& axis_ranges : active_entry_angles_) {
0755     if (axis_ranges.second.first > entry_tan[axis_ranges.first] ||
0756         axis_ranges.second.second < entry_tan[axis_ranges.first]) {
0757       return false;
0758     }
0759   }
0760   return true;
0761 }
0762
0763 //define this as a plug-in
0764 DEFINE_FWK_MODULE(PixelTestBeamValidation);