CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalDigiSelector.cc
Go to the documentation of this file.
1 #include <vector>
17 #include <TMath.h>
18 #include <iostream>
19 #include <set>
20 
21 
22 using namespace reco;
24 {
25  selectedEcalEBDigiCollection_ = ps.getParameter<std::string>("selectedEcalEBDigiCollection");
26  selectedEcalEEDigiCollection_ = ps.getParameter<std::string>("selectedEcalEEDigiCollection");
27 
28  barrelSuperClusterProducer_ =
29  consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("barrelSuperClusterProducer"));
30  endcapSuperClusterProducer_ =
31  consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("endcapSuperClusterProducer"));
32 
33 
34  EcalEBRecHitToken_ =
35  consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("EcalEBRecHitTag"));
36  EcalEERecHitToken_ =
37  consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("EcalEERecHitTag"));
38 
39  EcalEBDigiToken_ =
40  consumes<EBDigiCollection>(ps.getParameter<edm::InputTag>("EcalEBDigiTag"));
41  EcalEEDigiToken_ =
42  consumes<EEDigiCollection>(ps.getParameter<edm::InputTag>("EcalEEDigiTag"));
43 
44  cluster_pt_thresh_ = ps.getParameter<double>("cluster_pt_thresh");
45  single_cluster_thresh_ = ps.getParameter<double>("single_cluster_thresh");
46 
47  nclus_sel_ = ps.getParameter<int>("nclus_sel");
48  produces<EBDigiCollection>(selectedEcalEBDigiCollection_);
49  produces<EEDigiCollection>(selectedEcalEEDigiCollection_);
50 
51 }
52 
54 {
55 
56  //Get BarrelSuperClusters to start.
57  edm::Handle<reco::SuperClusterCollection> pBarrelSuperClusters;
58 
59  evt.getByToken(barrelSuperClusterProducer_, pBarrelSuperClusters);
60 
61  const reco::SuperClusterCollection & BarrelSuperClusters = *pBarrelSuperClusters;
62  //Got BarrelSuperClusters
63 
64  //Get BarrelSuperClusters to start.
65  edm::Handle<reco::SuperClusterCollection> pEndcapSuperClusters;
66 
67  evt.getByToken(endcapSuperClusterProducer_, pEndcapSuperClusters);
68 
69  const reco::SuperClusterCollection & EndcapSuperClusters = *pEndcapSuperClusters;
70  //Got EndcapSuperClusters
71 
72  reco::SuperClusterCollection saveBarrelSuperClusters;
73  reco::SuperClusterCollection saveEndcapSuperClusters;
74  bool meet_single_thresh = false;
75  //Loop over barrel superclusters, and apply threshold
76  for (int loop=0;loop<int(BarrelSuperClusters.size());loop++){
77  SuperCluster clus1 = BarrelSuperClusters[loop];
78  float eta1 = clus1.eta();
79  float energy1 = clus1.energy();
80  float theta1 = 2*atan(exp(-1.*eta1));
81  float cluspt1 = energy1 * sin(theta1);
82  if (cluspt1 > cluster_pt_thresh_){
83  saveBarrelSuperClusters.push_back(clus1);
84  if (cluspt1 > single_cluster_thresh_)
85  meet_single_thresh = true;
86 
87  }
88  }
89 
90  //Loop over endcap superclusters, and apply threshold
91  for (int loop=0;loop<int(EndcapSuperClusters.size());loop++){
92  SuperCluster clus1 = EndcapSuperClusters[loop];
93  float eta1 = clus1.eta();
94  float energy1 = clus1.energy();
95  float theta1 = 2*atan(exp(-1.*eta1));
96  float cluspt1 = energy1 * sin(theta1);
97  if (cluspt1 > cluster_pt_thresh_){
98  saveEndcapSuperClusters.push_back(clus1);
99  if (cluspt1 > single_cluster_thresh_)
100  meet_single_thresh = true;
101  }
102  }
103 
104  std::auto_ptr<EBDigiCollection> SEBDigiCol(new EBDigiCollection);
105  std::auto_ptr<EEDigiCollection> SEEDigiCol(new EEDigiCollection);
106  int TotClus = saveBarrelSuperClusters.size() + saveEndcapSuperClusters.size();
107 
108  if (TotClus >= nclus_sel_ || meet_single_thresh){
109 
110  if (saveBarrelSuperClusters.size() > 0){
111 
112 
113  edm::ESHandle<CaloTopology> pTopology;
114  es.get<CaloTopologyRecord>().get(pTopology);
115  const CaloTopology *topology = pTopology.product();
116 
117  //get barrel digi collection
119  const EBDigiCollection* digis=0;
120  evt.getByToken(EcalEBDigiToken_,pdigis);
121  digis = pdigis.product(); // get a ptr to the product
122 
124  const EcalRecHitCollection* rechits=0;
125  evt.getByToken(EcalEBRecHitToken_,prechits);
126  rechits = prechits.product(); // get a ptr to the product
127 
128 
129  if ( digis ) {
130  std::vector<DetId> saveTheseDetIds;
131  //pick out the detids for the 3x3 in each of the selected superclusters
132  for (int loop = 0;loop < int(saveBarrelSuperClusters.size());loop++){
133  SuperCluster clus1 = saveBarrelSuperClusters[loop];
134  CaloClusterPtr bcref = clus1.seed();
135  const BasicCluster *bc = bcref.get();
136  //Get the maximum detid
137  std::pair<DetId, float> EDetty =
138  EcalClusterTools::getMaximum(*bc,rechits);
139  //get the 3x3 array centered on maximum detid.
140  std::vector<DetId> detvec =
141  EcalClusterTools::matrixDetId(topology,EDetty.first, -1, 1, -1, 1);
142  //Loop over the 3x3
143  for (int ik = 0;ik<int(detvec.size());++ik)
144  saveTheseDetIds.push_back(detvec[ik]);
145 
146  }
147  for (int detloop=0; detloop < int(saveTheseDetIds.size());++detloop){
148  EBDetId detL = EBDetId(saveTheseDetIds[detloop]);
149 
150  for (EBDigiCollection::const_iterator blah = digis->begin();
151  blah!=digis->end();blah++){
152 
153  if (detL == blah->id()){
154  EBDataFrame myDigi = (*blah);
155  SEBDigiCol->push_back(detL);
156 
157  EBDataFrame df( SEBDigiCol->back());
158  for (int iq =0;iq<myDigi.size();++iq){
159  df.setSample(iq, myDigi.sample(iq).raw());
160  }
161  //ebcounter++;
162  }
163  }
164  //if (ebcounter >= int(saveTheseDetIds.size())) break;
165  }//loop over dets
166 
167  }
168 
169  }//If barrel superclusters need saving.
170 
171 
172  if (saveEndcapSuperClusters.size() > 0){
173 
174  edm::ESHandle<CaloTopology> pTopology;
175  es.get<CaloTopologyRecord>().get(pTopology);
176  const CaloTopology *topology = pTopology.product();
177 
178  //Get endcap rec hit collection
179  //get endcap digi collection
181  const EEDigiCollection* digis=0;
182  evt.getByToken(EcalEEDigiToken_,pdigis);
183  digis = pdigis.product(); // get a ptr to the product
184 
186  const EcalRecHitCollection* rechits=0;
187  evt.getByToken(EcalEERecHitToken_,prechits);
188  rechits = prechits.product(); // get a ptr to the product
189 
190  if ( digis ) {
191  //std::vector<DetId> saveTheseDetIds;
192  std::set<DetId> saveTheseDetIds;
193  //pick out the digis for the 3x3 in each of the selected superclusters
194  for (int loop = 0;loop < int(saveEndcapSuperClusters.size());loop++){
195  SuperCluster clus1 = saveEndcapSuperClusters[loop];
196  CaloClusterPtr bcref = clus1.seed();
197  const BasicCluster *bc = bcref.get();
198  //Get the maximum detid
199  std::pair<DetId, float> EDetty = EcalClusterTools::getMaximum(*bc,rechits);
200  //get the 3x3 array centered on maximum detid.
201  std::vector<DetId> detvec =
202  EcalClusterTools::matrixDetId(topology,EDetty.first, -1, 1, -1, 1);
203  //Loop over the 3x3
204  for (int ik = 0;ik<int(detvec.size());++ik)
205  saveTheseDetIds.insert(detvec[ik]);
206  }
207  int eecounter=0;
208  for (EEDigiCollection::const_iterator blah = digis->begin();
209  blah!=digis->end();blah++){
210  std::set<DetId>::const_iterator finder = saveTheseDetIds.find(blah->id());
211  if (finder!=saveTheseDetIds.end()){
212  EEDetId detL = EEDetId(*finder);
213 
214  if (detL == blah->id()){
215  EEDataFrame myDigi = (*blah);
216  SEEDigiCol->push_back(detL);
217  EEDataFrame df( SEEDigiCol->back());
218  for (int iq =0;iq<myDigi.size();++iq){
219  df.setSample(iq, myDigi.sample(iq).raw());
220  }
221  eecounter++;
222 
223  }
224  }
225  if (eecounter >= int(saveTheseDetIds.size())) break;
226  }//loop over digis
227  }
228  }//If endcap superclusters need saving.
229 
230  }//If we're actually saving stuff
231 
232  //Okay, either my collections have been filled with the requisite Digis, or they haven't.
233 
234 
235  //Empty collection, or full, still put in event.
236  SEBDigiCol->sort();
237  SEEDigiCol->sort();
238  evt.put(SEBDigiCol, selectedEcalEBDigiCollection_);
239  evt.put(SEEDigiCol, selectedEcalEEDigiCollection_);
240 
241 }
T getParameter(std::string const &) const
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
CaloTopology const * topology(0)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:160
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
EcalMGPASample sample(int i) const
Definition: EcalDataFrame.h:29
int loop
CMSSW
const_iterator begin() const
uint16_t raw() const
get the raw word
int size() const
Definition: EcalDataFrame.h:26
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:163
EcalDigiSelector(const edm::ParameterSet &ps)
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:121
double energy() const
cluster energy
Definition: CaloCluster.h:121
void produce(edm::Event &, const edm::EventSetup &)
T const * product() const
Definition: Handle.h:81
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
const_iterator end() const
void setSample(int i, EcalMGPASample sam)
Definition: EcalDataFrame.h:43
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:66