MtdAssociatorProducers/plugins/MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include "MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl.h"
0004 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0005
0006 using namespace reco;
0007 using namespace std;
0008
0009 /* Constructor */
0010
0011 MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl::MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl(
0012     edm::EDProductGetter const& productGetter, double energyCut, double timeCut, mtd::MTDGeomUtil& geomTools)
0013     : productGetter_(&productGetter), energyCut_(energyCut), timeCut_(timeCut), geomTools_(geomTools) {}
0014
0015 //
0016 //---member functions
0017 //
0018
0019 reco::RecoToSimCollectionMtd MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl::associateRecoToSim(
0020     const edm::Handle<FTLClusterCollection>& btlRecoClusH,
0021     const edm::Handle<FTLClusterCollection>& etlRecoClusH,
0022     const edm::Handle<MtdSimLayerClusterCollection>& simClusH) const {
0023   RecoToSimCollectionMtd outputCollection;
0024
0025   // -- get the collections
0026   std::array<edm::Handle<FTLClusterCollection>, 2> inputRecoClusH{{btlRecoClusH, etlRecoClusH}};
0027
0028   const auto& simClusters = *simClusH.product();
0029
0030   // -- create temporary map  DetId, SimClusterRef
0031   std::map<uint32_t, std::vector<MtdSimLayerClusterRef>> simClusIdsMap;
0032   for (auto simClusIt = simClusters.begin(); simClusIt != simClusters.end(); simClusIt++) {
0033     const auto& simClus = *simClusIt;
0034
0035     edm::Ref<MtdSimLayerClusterCollection> simClusterRef =
0036         edm::Ref<MtdSimLayerClusterCollection>(simClusH, simClusIt - simClusters.begin());
0037
0038     std::vector<std::pair<uint32_t, std::pair<uint8_t, uint8_t>>> detIdsAndRows = simClus.detIds_and_rows();
0039     std::vector<uint32_t> detIds(detIdsAndRows.size());
0040     std::transform(detIdsAndRows.begin(),
0041                    detIdsAndRows.end(),
0042                    detIds.begin(),
0043                    [](const std::pair<uint32_t, std::pair<uint8_t, uint8_t>>& pair) { return pair.first; });
0044
0045     // -- get the sensor module id from the first hit id of the MtdSimCluster
0046     DetId id(detIds[0]);
0047     uint32_t modId = geomTools_.sensorModuleId(id);
0048     simClusIdsMap[modId].push_back(simClusterRef);
0049
0050     LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl") << "Sim cluster  detId = " << modId << std::endl;
0051   }
0052
0053   for (auto const& recoClusH : inputRecoClusH) {
0054     // -- loop over detSetVec
0055     for (const auto& detSet : *recoClusH) {
0056       // -- loop over reco clusters
0057       for (const auto& recoClus : detSet) {
0058         MTDDetId clusId = recoClus.id();
0059
0060         LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl") << "Reco cluster : " << clusId;
0061
0062         std::vector<uint64_t> recoClusHitIds;
0063
0064         // -- loop over hits in the reco cluster and find their unique ids
0065         for (int ihit = 0; ihit < recoClus.size(); ++ihit) {
0066           int hit_row = recoClus.minHitRow() + recoClus.hitOffset()[ihit * 2];
0067           int hit_col = recoClus.minHitCol() + recoClus.hitOffset()[ihit * 2 + 1];
0068
0069           // -- Get an unique id from sensor module detId , row, column
0070           uint64_t uniqueId = static_cast<uint64_t>(clusId.rawId()) << 32;
0071           uniqueId |= hit_row << 16;
0072           uniqueId |= hit_col;
0073           recoClusHitIds.push_back(uniqueId);
0074
0075           LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0076               << " ======= cluster raw Id : " << clusId.rawId() << "  row : " << hit_row << "   column: " << hit_col
0077               << " recHit uniqueId : " << uniqueId;
0078         }
0079
0080         // -- loop over sim clusters and check if this reco clus shares some hits
0081         std::vector<MtdSimLayerClusterRef> simClusterRefs;
0082
0083         for (const auto& simClusterRef : simClusIdsMap[clusId.rawId()]) {
0084           const auto& simClus = *simClusterRef;
0085
0086           // get the hit ids of hits in the SimCluster
0087           std::vector<std::pair<uint32_t, std::pair<uint8_t, uint8_t>>> detIdsAndRows = simClus.detIds_and_rows();
0088           std::vector<uint64_t> simClusHitIds(detIdsAndRows.size());
0089           for (unsigned int i = 0; i < detIdsAndRows.size(); i++) {
0090             DetId id(detIdsAndRows[i].first);
0091             uint32_t modId = geomTools_.sensorModuleId(id);
0092             // build the unique id from sensor module detId , row, column
0093             uint64_t uniqueId = static_cast<uint64_t>(modId) << 32;
0094             uniqueId |= detIdsAndRows[i].second.first << 16;
0095             uniqueId |= detIdsAndRows[i].second.second;
0096             simClusHitIds.push_back(uniqueId);
0097           }
0098
0099           // -- Get shared hits
0100           std::vector<uint64_t> sharedHitIds;
0101           std::set_intersection(recoClusHitIds.begin(),
0102                                 recoClusHitIds.end(),
0103                                 simClusHitIds.begin(),
0104                                 simClusHitIds.end(),
0105                                 std::back_inserter(sharedHitIds));
0106
0107           if (sharedHitIds.empty())
0108             continue;
0109
0110           float dE = recoClus.energy() * 0.001 / simClus.simLCEnergy();  // reco cluster energy is in MeV!
0111           float dtSig = std::abs((recoClus.time() - simClus.simLCTime()) / recoClus.timeError());
0112
0113           // -- If the sim and reco clusters have common hits, fill the std:vector of sim clusters refs
0114           if (!sharedHitIds.empty() &&
0115               ((clusId.mtdSubDetector() == MTDDetId::BTL && dE < energyCut_ && dtSig < timeCut_) ||
0116                (clusId.mtdSubDetector() == MTDDetId::ETL && dtSig < timeCut_))) {
0117             simClusterRefs.push_back(simClusterRef);
0118
0119             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0120                 << "RecoToSim --> Found " << sharedHitIds.size() << " shared hits";
0121             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0122                 << "E_recoClus = " << recoClus.energy() << "   E_simClus = " << simClus.simLCEnergy()
0123                 << "   E_recoClus/E_simClus = " << dE;
0124             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0125                 << "(t_recoClus-t_simClus)/sigma_t = " << dtSig;
0126           }
0127         }  // -- end loop over sim clus refs
0128
0129         // -- Now fill the output collection
0130         edm::Ref<edmNew::DetSetVector<FTLCluster>, FTLCluster> recoClusterRef = edmNew::makeRefTo(recoClusH, &recoClus);
0131         outputCollection.emplace_back(recoClusterRef, simClusterRefs);
0132
0133       }  // -- end loop over reco clus
0134     }  // -- end loop over detsetclus
0135   }
0136
0137   outputCollection.post_insert();
0138   return outputCollection;
0139 }
0140
0141 reco::SimToRecoCollectionMtd MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl::associateSimToReco(
0142     const edm::Handle<FTLClusterCollection>& btlRecoClusH,
0143     const edm::Handle<FTLClusterCollection>& etlRecoClusH,
0144     const edm::Handle<MtdSimLayerClusterCollection>& simClusH) const {
0145   SimToRecoCollectionMtd outputCollection;
0146
0147   // -- get the collections
0148   const auto& simClusters = *simClusH.product();
0149
0150   std::array<edm::Handle<FTLClusterCollection>, 2> inputH{{btlRecoClusH, etlRecoClusH}};
0151
0152   // -- loop over MtdSimLayerClusters
0153   for (auto simClusIt = simClusters.begin(); simClusIt != simClusters.end(); simClusIt++) {
0154     const auto& simClus = *simClusIt;
0155
0156     // get the hit ids of hits in the SimCluster
0157     std::vector<std::pair<uint32_t, std::pair<uint8_t, uint8_t>>> detIdsAndRows = simClus.detIds_and_rows();
0158     std::vector<uint64_t> simClusHitIds(detIdsAndRows.size());
0159     uint32_t modId = 0;
0160     for (unsigned int i = 0; i < detIdsAndRows.size(); i++) {
0161       DetId id(detIdsAndRows[i].first);
0162       modId = geomTools_.sensorModuleId(id);
0163       // build the unique id from sensor module detId , row, column
0164       uint64_t uniqueId = static_cast<uint64_t>(modId) << 32;
0165       uniqueId |= detIdsAndRows[i].second.first << 16;
0166       uniqueId |= detIdsAndRows[i].second.second;
0167       simClusHitIds.push_back(uniqueId);
0168     }
0169
0170     DetId simClusId(modId);
0171
0172     std::vector<FTLClusterRef> recoClusterRefs;
0173
0174     // -- loop over reco clusters
0175     for (auto const& recoClusH : inputH) {
0176       // -- loop over detSetVec
0177       for (const auto& detSet : *recoClusH) {
0178         if (detSet.id() != simClusId.rawId())
0179           continue;
0180
0181         // -- loop over reco clusters
0182         for (const auto& recoClus : detSet) {
0183           MTDDetId clusId = recoClus.id();
0184
0185           std::vector<uint64_t> recoClusHitIds;
0186
0187           // -- loop over hits in the reco cluster and find their ids
0188           for (int ihit = 0; ihit < recoClus.size(); ++ihit) {
0189             int hit_row = recoClus.minHitRow() + recoClus.hitOffset()[ihit * 2];
0190             int hit_col = recoClus.minHitCol() + recoClus.hitOffset()[ihit * 2 + 1];
0191
0192             // -- Get an unique id from sensor module detId , row, column
0193             uint64_t uniqueId = static_cast<uint64_t>(clusId.rawId()) << 32;
0194             uniqueId |= hit_row << 16;
0195             uniqueId |= hit_col;
0196             recoClusHitIds.push_back(uniqueId);
0197
0198             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0199                 << " ======= cluster raw Id : " << clusId.rawId() << "  row : " << hit_row << "   column: " << hit_col
0200                 << " recHit uniqueId : " << uniqueId;
0201           }
0202
0203           // -- Get shared hits
0204           std::vector<uint64_t> sharedHitIds;
0205           std::set_intersection(simClusHitIds.begin(),
0206                                 simClusHitIds.end(),
0207                                 recoClusHitIds.begin(),
0208                                 recoClusHitIds.end(),
0209                                 std::back_inserter(sharedHitIds));
0210           if (sharedHitIds.empty())
0211             continue;
0212
0213           float dE = recoClus.energy() * 0.001 / simClus.simLCEnergy();  // reco cluster energy is in MeV
0214           float dtSig = std::abs((recoClus.time() - simClus.simLCTime()) / recoClus.timeError());
0215
0216           // -- If the sim and reco clusters have common hits, fill the std:vector of reco clusters refs
0217           if (!sharedHitIds.empty() &&
0218               ((clusId.mtdSubDetector() == MTDDetId::BTL && dE < energyCut_ && dtSig < timeCut_) ||
0219                (clusId.mtdSubDetector() == MTDDetId::ETL && dtSig < timeCut_))) {
0220             // Create a persistent edm::Ref to the cluster
0221             edm::Ref<edmNew::DetSetVector<FTLCluster>, FTLCluster> recoClusterRef =
0222                 edmNew::makeRefTo(recoClusH, &recoClus);
0223             recoClusterRefs.push_back(recoClusterRef);
0224
0225             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0226                 << "SimToReco --> Found " << sharedHitIds.size() << " shared hits";
0227             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0228                 << "E_recoClus = " << recoClus.energy() << "   E_simClus = " << simClus.simLCEnergy()
0229                 << "   E_recoClus/E_simClus = " << dE;
0230             LogDebug("MtdRecoClusterToSimLayerClusterAssociatorByHitsImpl")
0231                 << "(t_recoClus-t_simClus)/sigma_t = " << dtSig;
0232           }
0233
0234         }  // end loop ove reco clus
0235       }  // end loop over detsets
0236     }
0237
0238     // -- Now fill the output collection
0239     edm::Ref<MtdSimLayerClusterCollection> simClusterRef =
0240         edm::Ref<MtdSimLayerClusterCollection>(simClusH, simClusIt - simClusters.begin());
0241     outputCollection.emplace_back(simClusterRef, recoClusterRefs);
0242
0243   }  // -- end loop over sim clusters
0244
0245   outputCollection.post_insert();
0246   return outputCollection;
0247 }