SiStripChannelChargeFilter/src/ClusterMTCCFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     SiStripChannelChargeFilter
0004 // Class  :     ClusterMTCCFilter
0005 //
0006 //
0007 // Original Author:  dkcira
0008
0009 #include "EventFilter/SiStripChannelChargeFilter/interface/ClusterMTCCFilter.h"
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/Common/interface/DetSetVector.h"
0012 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0013 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 using namespace std;
0017
0018 namespace cms {
0019
0020   ClusterMTCCFilter::ClusterMTCCFilter(const edm::ParameterSet& ps) {
0021     //
0022     ModulesToBeExcluded.clear();
0023     ModulesToBeExcluded = ps.getParameter<std::vector<unsigned>>("ModulesToBeExcluded");
0024     edm::LogInfo("ClusterMTCCFilter") << "Clusters from " << ModulesToBeExcluded.size()
0025                                       << " modules will be ignored in the filter:";
0026     for (std::vector<uint32_t>::const_iterator imod = ModulesToBeExcluded.begin(); imod != ModulesToBeExcluded.end();
0027          imod++) {
0028       edm::LogInfo("ClusterMTCCFilter") << *imod;
0029     }
0030     //
0031     ChargeThresholdTIB = ps.getParameter<int>("ChargeThresholdTIB");
0032     ChargeThresholdTOB = ps.getParameter<int>("ChargeThresholdTOB");
0033     ChargeThresholdTEC = ps.getParameter<int>("ChargeThresholdTEC");
0034     MinClustersDiffComponents = ps.getParameter<int>("MinClustersDiffComponents");
0035     clusterProducer = ps.getParameter<string>("ClusterProducer");
0036     tTopoToken_ = esConsumes();
0037     //
0038     produces<int>();
0039     produces<unsigned int>();
0040     produces<map<unsigned int, vector<SiStripCluster>>>();
0041   }
0042
0043   bool ClusterMTCCFilter::filter(edm::Event& e, edm::EventSetup const& c) {
0044     const auto& tTopo = c.getData(tTopoToken_);
0045
0046     //get SiStripCluster
0047     edm::Handle<edm::DetSetVector<SiStripCluster>> h;
0048     e.getByLabel(clusterProducer, h);
0049
0050     //
0051     unsigned int sum_of_cluster_charges = 0;
0052     clusters_in_subcomponents.clear();
0053     // first find all clusters that are over the threshold
0054     for (edm::DetSetVector<SiStripCluster>::const_iterator it = h->begin(); it != h->end(); it++) {
0055       for (vector<SiStripCluster>::const_iterator vit = (it->data).begin(); vit != (it->data).end(); vit++) {
0056         // calculate sum of amplitudes
0057         unsigned int amplclus = 0;
0058         for (auto ia = vit->amplitudes().begin(); ia != vit->amplitudes().end(); ia++) {
0059           if ((*ia) > 0)
0060             amplclus += (*ia);  // why should this be negative?
0061         }
0062         sum_of_cluster_charges += amplclus;
0063         DetId thedetId = DetId(it->detId());
0064         unsigned int generalized_layer = 0;
0065         bool exclude_this_detid = false;
0066         for (std::vector<uint32_t>::const_iterator imod = ModulesToBeExcluded.begin();
0067              imod != ModulesToBeExcluded.end();
0068              imod++) {
0069           if (*imod == thedetId.rawId())
0070             exclude_this_detid = true;  // found in exclusion list
0071         }
0072         // apply different thresholds for TIB/TOB/TEC
0073         if (!exclude_this_detid) {  // only consider if not in exclusion list
0074           if ((thedetId.subdetId() == StripSubdetector::TIB && amplclus > ChargeThresholdTIB) ||
0075               (thedetId.subdetId() == StripSubdetector::TOB && amplclus > ChargeThresholdTOB) ||
0076               (thedetId.subdetId() == StripSubdetector::TEC && amplclus > ChargeThresholdTEC)) {
0077             // calculate generalized_layer:  31 = TIB1, 32 = TIB2, 33 = TIB3, 50 = TOB, 60 = TEC
0078             if (thedetId.subdetId() == StripSubdetector::TIB) {
0079               generalized_layer =
0080                   10 * thedetId.subdetId() + tTopo.tibLayer(thedetId.rawId()) + tTopo.tibStereo(thedetId.rawId());
0081               if (tTopo.tibLayer(thedetId.rawId()) == 2) {
0082                 generalized_layer++;
0083                 if (tTopo.tibGlued(thedetId.rawId()))
0084                   edm::LogError("ClusterMTCCFilter") << "WRONGGGG" << endl;
0085               }
0086             } else {
0087               generalized_layer = 10 * thedetId.subdetId();
0088               if (thedetId.subdetId() == StripSubdetector::TOB) {
0089                 generalized_layer += tTopo.tobLayer(thedetId.rawId());
0090               }
0091             }
0092             // fill clusters_in_subcomponents
0093             map<unsigned int, vector<SiStripCluster>>::iterator layer_it =
0094                 clusters_in_subcomponents.find(generalized_layer);
0095             if (layer_it == clusters_in_subcomponents
0096                                 .end()) {  // if layer not found yet, create DATA vector and generate map KEY + DATA
0097               vector<SiStripCluster> local_vector;
0098               local_vector.push_back(*vit);
0099               clusters_in_subcomponents.insert(std::make_pair(generalized_layer, local_vector));
0100             } else {  // push into already existing vector
0101               (layer_it->second).push_back(*vit);
0102             }
0103           }
0104         }
0105       }
0106     }
0107
0108     bool decision = false;  // default value, only accept if set true in this loop
0109     unsigned int nr_of_subcomps_with_clusters = 0;
0110     // dk: 2006.08.24 - change filter decision as proposed by V. Ciulli. || TIB1 TIB2 counted as 1, TEC excluded
0111     //  if( clusters_in_subcomponents[31].size()>0 ) nr_of_subcomps_with_clusters++; // TIB1
0112     //  if( clusters_in_subcomponents[32].size()>0 ) nr_of_subcomps_with_clusters++; // TIB2
0113     //  if( clusters_in_subcomponents[60].size()>0 ) nr_of_subcomps_with_clusters++; // TEC
0114     if (!clusters_in_subcomponents[31].empty() || !clusters_in_subcomponents[32].empty())
0115       nr_of_subcomps_with_clusters++;  // TIB1 || TIB2
0116     if (!clusters_in_subcomponents[33].empty())
0117       nr_of_subcomps_with_clusters++;  // TIB3
0118     if (!clusters_in_subcomponents[51].empty())
0119       nr_of_subcomps_with_clusters++;  // TOB1
0120     if (!clusters_in_subcomponents[52].empty())
0121       nr_of_subcomps_with_clusters++;  // TOB2
0122     if (nr_of_subcomps_with_clusters >=
0123         MinClustersDiffComponents  // more than 'MinClustersDiffComponents' components have at least 1 cluster
0124     ) {
0125       decision = true;  // accept event
0126     }
0127
0128     e.put(std::make_unique<int>(decision));
0129
0130     e.put(std::make_unique<unsigned int>(sum_of_cluster_charges));
0131
0132     e.put(std::make_unique<std::map<unsigned int, std::vector<SiStripCluster>>>(clusters_in_subcomponents));
0133
0134     return decision;
0135   }
0136 }  // namespace cms