CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFRecHitProducerHCAL.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducerHCAL.h"
2 
3 #include <memory>
4 
5 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
6 
7 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
8 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
9 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
10 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
11 
12 #include "DataFormats/DetId/interface/DetId.h"
13 #include "DataFormats/EcalDetId/interface/EEDetId.h"
14 #include "DataFormats/EcalDetId/interface/EBDetId.h"
15 
16 #include "FWCore/MessageLogger/interface/MessageLogger.h"
17 #include "FWCore/Framework/interface/ESHandle.h"
18 #include "FWCore/Framework/interface/EventSetup.h"
19 
20 
21 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
22 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
23 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
24 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
25 #include "Geometry/Records/interface/CaloGeometryRecord.h"
26 
27 #include "Geometry/CaloTopology/interface/HcalTopology.h"
28 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
29 
30 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
31 
32 #include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
33 #include "RecoCaloTools/Navigation/interface/CaloTowerNavigator.h"
34 
35 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
36 
37 
38 using namespace std;
39 using namespace edm;
40 
41 PFRecHitProducerHCAL::PFRecHitProducerHCAL(const edm::ParameterSet& iConfig)
42  : PFRecHitProducer( iConfig )
43 {
44 
45 
46 
47  // access to the collections of rechits
48 
49 
50  inputTagHcalRecHitsHBHE_ =
51  iConfig.getParameter<InputTag>("hcalRecHitsHBHE");
52 
53  inputTagHcalRecHitsHF_ =
54  iConfig.getParameter<InputTag>("hcalRecHitsHF");
55 
56 
57  inputTagCaloTowers_ =
58  iConfig.getParameter<InputTag>("caloTowers");
59 
60  thresh_HF_ =
61  iConfig.getParameter<double>("thresh_HF");
62  navigation_HF_ =
63  iConfig.getParameter<bool>("navigation_HF");
64  weight_HFem_ =
65  iConfig.getParameter<double>("weight_HFem");
66  weight_HFhad_ =
67  iConfig.getParameter<double>("weight_HFhad");
68 
69  HCAL_Calib_ =
70  iConfig.getParameter<bool>("HCAL_Calib");
71  HF_Calib_ =
72  iConfig.getParameter<bool>("HF_Calib");
73  max_Calib_ =
74  iConfig.getParameter<double>("Max_Calib");
75 
76  shortFibre_Cut = iConfig.getParameter<double>("ShortFibre_Cut");
77  longFibre_Fraction = iConfig.getParameter<double>("LongFibre_Fraction");
78 
79  longFibre_Cut = iConfig.getParameter<double>("LongFibre_Cut");
80  shortFibre_Fraction = iConfig.getParameter<double>("ShortFibre_Fraction");
81 
82  applyLongShortDPG_ = iConfig.getParameter<bool>("ApplyLongShortDPG");
83 
84  longShortFibre_Cut = iConfig.getParameter<double>("LongShortFibre_Cut");
85  minShortTiming_Cut = iConfig.getParameter<double>("MinShortTiming_Cut");
86  maxShortTiming_Cut = iConfig.getParameter<double>("MaxShortTiming_Cut");
87  minLongTiming_Cut = iConfig.getParameter<double>("MinLongTiming_Cut");
88  maxLongTiming_Cut = iConfig.getParameter<double>("MaxLongTiming_Cut");
89 
90  applyTimeDPG_ = iConfig.getParameter<bool>("ApplyTimeDPG");
91  applyPulseDPG_ = iConfig.getParameter<bool>("ApplyPulseDPG");
92 
93  ECAL_Compensate_ = iConfig.getParameter<bool>("ECAL_Compensate");
94  ECAL_Threshold_ = iConfig.getParameter<double>("ECAL_Threshold");
95  ECAL_Compensation_ = iConfig.getParameter<double>("ECAL_Compensation");
96  ECAL_Dead_Code_ = iConfig.getParameter<unsigned int>("ECAL_Dead_Code");
97 
98  EM_Depth_ = iConfig.getParameter<double>("EM_Depth");
99  HAD_Depth_ = iConfig.getParameter<double>("HAD_Depth");
100 
101  //--ab
102  produces<reco::PFRecHitCollection>("HFHAD").setBranchAlias("HFHADRecHits");
103  produces<reco::PFRecHitCollection>("HFEM").setBranchAlias("HFEMRecHits");
104  //--ab
105 }
106 
107 
108 
109 PFRecHitProducerHCAL::~PFRecHitProducerHCAL() {}
110 
111 
112 
113 void PFRecHitProducerHCAL::createRecHits(vector<reco::PFRecHit>& rechits,
114  vector<reco::PFRecHit>& rechitsCleaned,
115  edm::Event& iEvent,
116  const edm::EventSetup& iSetup ) {
117 
118 
119  // this map is necessary to find the rechit neighbours efficiently
120  //C but I should think about using Florian's hashed index to do this.
121  //C in which case the map might not be necessary anymore
122  //C however the hashed index does not seem to be implemented for HCAL
123  //
124  // the key of this map is detId.
125  // the value is the index in the rechits vector
126  map<unsigned, unsigned > idSortedRecHits;
127  map<unsigned, unsigned > idSortedRecHitsHFEM;
128  map<unsigned, unsigned > idSortedRecHitsHFHAD;
129  typedef map<unsigned, unsigned >::iterator IDH;
130 
131 
132  edm::ESHandle<CaloGeometry> geoHandle;
133  iSetup.get<CaloGeometryRecord>().get(geoHandle);
134 
135  // get the hcalBarrel geometry
136  const CaloSubdetectorGeometry *hcalBarrelGeometry =
137  geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
138 
139  // get the endcap geometry
140  const CaloSubdetectorGeometry *hcalEndcapGeometry =
141  geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalEndcap);
142 
143  //--ab
144  auto_ptr< vector<reco::PFRecHit> > HFHADRecHits( new vector<reco::PFRecHit> );
145  auto_ptr< vector<reco::PFRecHit> > HFEMRecHits( new vector<reco::PFRecHit> );
146  //--ab
147 
148  // 2 possibilities to make HCAL clustering :
149  // - from the HCAL rechits
150  // - from the CaloTowers.
151  // ultimately, clustering will be done taking CaloTowers as an
152  // input. This possibility is currently under investigation, and
153  // was thus made optional.
154 
155  // in the first step, we will fill the map of PFRecHits hcalrechits
156  // either from CaloTowers or from HCAL rechits.
157 
158  // in the second step, we will perform clustering on this map.
159 
160  if( !(inputTagCaloTowers_ == InputTag()) ) {
161 
162  edm::Handle<CaloTowerCollection> caloTowers;
163  CaloTowerTopology caloTowerTopology;
164  const CaloSubdetectorGeometry *caloTowerGeometry = 0;
165  // = geometry_->getSubdetectorGeometry(id)
166 
167  // get calotowers
168  bool found = iEvent.getByLabel(inputTagCaloTowers_,
169  caloTowers);
170 
171  if(!found) {
172  ostringstream err;
173  err<<"could not find rechits "<<inputTagCaloTowers_;
174  LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
175 
176  throw cms::Exception( "MissingProduct", err.str());
177  }
178  else {
179  assert( caloTowers.isValid() );
180 
181  // get HF rechits
182  edm::Handle<HFRecHitCollection> hfHandle;
183  found = iEvent.getByLabel(inputTagHcalRecHitsHF_,
184  hfHandle);
185 
186  if(!found) {
187  ostringstream err;
188  err<<"could not find HF rechits "<<inputTagHcalRecHitsHF_;
189  LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
190 
191  throw cms::Exception( "MissingProduct", err.str());
192  }
193  else {
194  assert( hfHandle.isValid() );
195  }
196 
197  // get HBHE rechits
198  edm::Handle<HBHERecHitCollection> hbheHandle;
199  found = iEvent.getByLabel(inputTagHcalRecHitsHBHE_,
200  hbheHandle);
201 
202  if(!found) {
203  ostringstream err;
204  err<<"could not find HBHE rechits "<<inputTagHcalRecHitsHBHE_;
205  LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
206 
207  throw cms::Exception( "MissingProduct", err.str());
208  }
209  else {
210  assert( hbheHandle.isValid() );
211  }
212 
213  // create rechits
214  typedef CaloTowerCollection::const_iterator ICT;
215 
216  for(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
217 
218  const CaloTower& ct = (*ict);
219 
220  //C
221  if(!caloTowerGeometry)
222  caloTowerGeometry = geoHandle->getSubdetectorGeometry(ct.id());
223 
224 
225  // get the hadronic energy.
226 
227  // Mike: Just ask for the Hadronic part only now!
228  // Patrick : ARGH ! While this is ok for the HCAL, this is
229  // just wrong for the HF (in which em/had are artificially
230  // separated.
231  double energy = ct.hadEnergy();
232  //Auguste: Photons in HF have no hadEnergy in fastsim: -> all RecHit collections are empty with photons.
233  double energyEM = ct.emEnergy(); // For HF !
234  //so test the total energy to deal with the photons in HF:
235  if( (energy+energyEM) < 1e-9 ) continue;
236 
237  assert( ct.constituentsSize() );
238  //Mike: The DetId will be taken by the first Hadronic constituent
239  // of the tower. That is only what we need
240 
241 
242  //get the constituents of the tower
243  const std::vector<DetId>& hits = ct.constituents();
244  const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id());
245 
246  /*
247  for(unsigned int i=0;i< hits.size();++i) {
248  if(hits[i].det()==DetId::Hcal) {
249  HcalDetId did = hits[i];
250  if ( did.subdet()==HcalEndcap || did.subdet()==HcalForward ) {
251  //double en = hits[i].energy();
252  int ieta = did.ieta();
253  const CaloCellGeometry *thisCell = hcalEndcapGeometry->getGeometry(did);
254  const GlobalPoint& position = thisCell->getPosition();
255  if ( abs(ieta) > 27 && abs(ieta) < 33 && energy > 10. ) {
256  std::cout << "HE/HF hit " << i << " at eta = " << ieta
257  << " with CT energy = " << energy
258  << " at eta, z (hit) = " << position.eta() << " " << position.z()
259  << " at eta, z (cte) = " << ct.emPosition().eta() << " " << ct.emPosition().z()
260  << " at eta, z (cth) = " << ct.hadPosition().eta() << " " << ct.hadPosition().z()
261  << " at eta, z (cto) = " << ct.eta() << " " << ct.vz()
262  << std::endl;
263  }
264  }
265  }
266  }
267  */
268 
269  //Reserve the DetId we are looking for:
270 
271  HcalDetId detid;
272  // EcalDetId edetid;
273  bool foundHCALConstituent = false;
274 
275 
276  //Loop on the calotower constituents and search for HCAL
277  double dead = 0.;
278  double alive = 0.;
279  for(unsigned int i=0;i< hits.size();++i) {
280  if(hits[i].det()==DetId::Hcal) {
281  foundHCALConstituent = true;
282  detid = hits[i];
283  // An HCAL tower was found: Look for dead ECAL channels in the same CaloTower.
284  if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
285  for(unsigned int j=0;j<allConstituents.size();++j) {
286  if ( allConstituents[j].det()==DetId::Ecal ) {
287  alive += 1.;
288  EcalChannelStatus::const_iterator chIt = theEcalChStatus->find(allConstituents[j]);
289  unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
290  if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
291  }
292  }
293  }
294  // Protection: tower 29 in HF is merged with tower 30.
295  // Just take the position of tower 30 in that case.
296  if ( detid.subdet() == HcalForward && abs(detid.ieta()) == 29 ) continue;
297  break;
298  }
299  }
300 
301  // In case of dead ECAL channel, rescale the HCAL energy...
302  double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
303 
304  reco::PFRecHit* pfrh = 0;
305  reco::PFRecHit* pfrhCleaned = 0;
306  //---ab: need 2 rechits for the HF:
307  reco::PFRecHit* pfrhHFEM = 0;
308  reco::PFRecHit* pfrhHFHAD = 0;
309  reco::PFRecHit* pfrhHFEMCleaned = 0;
310  reco::PFRecHit* pfrhHFHADCleaned = 0;
311  reco::PFRecHit* pfrhHFEMCleaned29 = 0;
312  reco::PFRecHit* pfrhHFHADCleaned29 = 0;
313 
314  if(foundHCALConstituent)
315  {
316  // std::cout << ", new Energy = " << energy << std::endl;
317  switch( detid.subdet() ) {
318  case HcalBarrel:
319  {
320  if(energy < thresh_Barrel_ ) continue;
321 
322  /*
323  // Check the timing
324  if ( energy > 5. ) {
325  for(unsigned int i=0;i< hits.size();++i) {
326  if( hits[i].det() != DetId::Hcal ) continue;
327  HcalDetId theDetId = hits[i];
328  typedef HBHERecHitCollection::const_iterator iHBHE;
329  iHBHE theHit = hbheHandle->find(theDetId);
330  if ( theHit != hbheHandle->end() )
331  std::cout << "HCAL hit : "
332  << theDetId.ieta() << " " << theDetId.iphi() << " "
333  << theHit->energy() << " " << theHit->time() << std::endl;
334  }
335  }
336  */
337 
338 
339  if ( HCAL_Calib_ ) energy *= std::min(max_Calib_,myPFCorr->getValues(detid)->getValue());
340  //if ( rescaleFactor > 1. )
341  // std::cout << "Barrel HCAL energy rescaled from = " << energy << " to " << energy*rescaleFactor << std::endl;
342  if ( rescaleFactor > 1. ) {
343  pfrhCleaned = createHcalRecHit( detid,
344  energy,
345  PFLayer::HCAL_BARREL1,
346  hcalBarrelGeometry,
347  ct.id().rawId() );
348  pfrhCleaned->setRescale(rescaleFactor);
349  energy *= rescaleFactor;
350  }
351  pfrh = createHcalRecHit( detid,
352  energy,
353  PFLayer::HCAL_BARREL1,
354  hcalBarrelGeometry,
355  ct.id().rawId() );
356  pfrh->setRescale(rescaleFactor);
357  }
358  break;
359  case HcalEndcap:
360  {
361  if(energy < thresh_Endcap_ ) continue;
362 
363  /*
364  // Check the timing
365  if ( energy > 5. ) {
366  for(unsigned int i=0;i< hits.size();++i) {
367  if( hits[i].det() != DetId::Hcal ) continue;
368  HcalDetId theDetId = hits[i];
369  typedef HBHERecHitCollection::const_iterator iHBHE;
370  iHBHE theHit = hbheHandle->find(theDetId);
371  if ( theHit != hbheHandle->end() )
372  std::cout << "HCAL hit : "
373  << theDetId.ieta() << " " << theDetId.iphi() << " "
374  << theHit->energy() << " " << theHit->time() << std::endl;
375  }
376  }
377  */
378 
379  if ( HCAL_Calib_ ) energy *= std::min(max_Calib_,myPFCorr->getValues(detid)->getValue());
380  //if ( rescaleFactor > 1. )
381  // std::cout << "End-cap HCAL energy rescaled from = " << energy << " to " << energy*rescaleFactor << std::endl;
382  if ( rescaleFactor > 1. ) {
383  pfrhCleaned = createHcalRecHit( detid,
384  energy,
385  PFLayer::HCAL_ENDCAP,
386  hcalEndcapGeometry,
387  ct.id().rawId() );
388  pfrhCleaned->setRescale(rescaleFactor);
389  energy *= rescaleFactor;
390  }
391  pfrh = createHcalRecHit( detid,
392  energy,
393  PFLayer::HCAL_ENDCAP,
394  hcalEndcapGeometry,
395  ct.id().rawId() );
396  pfrh->setRescale(rescaleFactor);
397  }
398  break;
399  case HcalOuter:
400  {
401  }
402  break;
403  case HcalForward:
404  {
405  //---ab: 2 rechits for HF:
406  //double energyemHF = weight_HFem_*ct.emEnergy();
407  //double energyhadHF = weight_HFhad_*ct.hadEnergy();
408  double energyemHF = weight_HFem_ * energyEM;
409  double energyhadHF = weight_HFhad_ * energy;
410  // Some energy in the tower !
411  if((energyemHF+energyhadHF) < thresh_HF_ ) continue;
412 
413  // Some cleaning in the HF
414  double longFibre = energyemHF + energyhadHF/2.;
415  double shortFibre = energyhadHF/2.;
416  int ieta = detid.ieta();
417  int iphi = detid.iphi();
418  HcalDetId theLongDetId (HcalForward, ieta, iphi, 1);
419  HcalDetId theShortDetId (HcalForward, ieta, iphi, 2);
420  typedef HFRecHitCollection::const_iterator iHF;
421  iHF theLongHit = hfHandle->find(theLongDetId);
422  iHF theShortHit = hfHandle->find(theShortDetId);
423  //
424  double theLongHitEnergy = 0.;
425  double theShortHitEnergy = 0.;
426  bool flagShortDPG = false;
427  bool flagLongDPG = false;
428  bool flagShortTimeDPG = false;
429  bool flagLongTimeDPG = false;
430  bool flagShortPulseDPG = false;
431  bool flagLongPulseDPG = false;
432  //
433  if ( theLongHit != hfHandle->end() ) {
434  theLongHitEnergy = theLongHit->energy();
435  flagLongDPG = applyLongShortDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
436  flagLongTimeDPG = applyTimeDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
437  flagLongPulseDPG = applyPulseDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
438  }
439  //
440  if ( theShortHit != hfHandle->end() ) {
441  theShortHitEnergy = theShortHit->energy();
442  flagShortDPG = applyLongShortDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
443  flagShortTimeDPG = applyTimeDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
444  flagShortPulseDPG = applyPulseDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
445  }
446 
447  // Then check the timing in short and long fibres in all other towers.
448  if ( theShortHitEnergy > longShortFibre_Cut &&
449  ( theShortHit->time() < minShortTiming_Cut ||
450  theShortHit->time() > maxShortTiming_Cut ||
451  flagShortTimeDPG || flagShortPulseDPG ) ) {
452  // rescaleFactor = 0. ;
453  pfrhHFHADCleaned = createHcalRecHit( detid,
454  theShortHitEnergy,
455  PFLayer::HF_HAD,
456  hcalEndcapGeometry,
457  ct.id().rawId() );
458  pfrhHFHADCleaned->setRescale(theShortHit->time());
459  /*
460  std::cout << "ieta/iphi = " << ieta << " " << iphi
461  << ", Energy em/had/long/short = "
462  << energyemHF << " " << energyhadHF << " "
463  << theLongHitEnergy << " " << theShortHitEnergy << " "
464  << ". Time = " << theShortHit->time() << " "
465  << ". The short and long flags : "
466  << flagShortDPG << " " << flagLongDPG << " "
467  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
468  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
469  << ". Short fibres were cleaned." << std::endl;
470  */
471  shortFibre -= theShortHitEnergy;
472  theShortHitEnergy = 0.;
473  }
474 
475 
476  if ( theLongHitEnergy > longShortFibre_Cut &&
477  ( theLongHit->time() < minLongTiming_Cut ||
478  theLongHit->time() > maxLongTiming_Cut ||
479  flagLongTimeDPG || flagLongPulseDPG ) ) {
480  //rescaleFactor = 0. ;
481  pfrhHFEMCleaned = createHcalRecHit( detid,
482  theLongHitEnergy,
483  PFLayer::HF_EM,
484  hcalEndcapGeometry,
485  ct.id().rawId() );
486  pfrhHFEMCleaned->setRescale(theLongHit->time());
487  /*
488  std::cout << "ieta/iphi = " << ieta << " " << iphi
489  << ", Energy em/had/long/short = "
490  << energyemHF << " " << energyhadHF << " "
491  << theLongHitEnergy << " " << theShortHitEnergy << " "
492  << ". Time = " << theLongHit->time() << " "
493  << ". The short and long flags : "
494  << flagShortDPG << " " << flagLongDPG << " "
495  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
496  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
497  << ". Long fibres were cleaned." << std::endl;
498  */
499  longFibre -= theLongHitEnergy;
500  theLongHitEnergy = 0.;
501  }
502 
503  // Some energy must be in the long fibres is there is some energy in the short fibres !
504  if ( theShortHitEnergy > shortFibre_Cut &&
505  ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
506  flagShortDPG ) ) {
507  // Check if the long-fibre hit was not cleaned already (because hot)
508  // In this case don't apply the cleaning
509  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId);
510  unsigned theStatusValue = theStatus->getValue();
511  // The channel is killed
512  if ( !theStatusValue ) {
513  // rescaleFactor = 0. ;
514  pfrhHFHADCleaned = createHcalRecHit( detid,
515  theShortHitEnergy,
516  PFLayer::HF_HAD,
517  hcalEndcapGeometry,
518  ct.id().rawId() );
519  pfrhHFHADCleaned->setRescale(theShortHit->time());
520  /*
521  std::cout << "ieta/iphi = " << ieta << " " << iphi
522  << ", Energy em/had/long/short = "
523  << energyemHF << " " << energyhadHF << " "
524  << theLongHitEnergy << " " << theShortHitEnergy << " "
525  << ". Time = " << theShortHit->time() << " "
526  << ". The status value is " << theStatusValue
527  << ". The short and long flags : "
528  << flagShortDPG << " " << flagLongDPG << " "
529  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
530  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
531  << ". Short fibres were cleaned." << std::endl;
532  */
533  shortFibre -= theShortHitEnergy;
534  theShortHitEnergy = 0.;
535  }
536  }
537 
538  if ( theLongHitEnergy > longFibre_Cut &&
539  ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
540  flagLongDPG ) ) {
541  // Check if the long-fibre hit was not cleaned already (because hot)
542  // In this case don't apply the cleaning
543  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId);
544  unsigned theStatusValue = theStatus->getValue();
545  // The channel is killed
546  if ( !theStatusValue ) {
547  //rescaleFactor = 0. ;
548  pfrhHFEMCleaned = createHcalRecHit( detid,
549  theLongHitEnergy,
550  PFLayer::HF_EM,
551  hcalEndcapGeometry,
552  ct.id().rawId() );
553  pfrhHFEMCleaned->setRescale(theLongHit->time());
554  /*
555  std::cout << "ieta/iphi = " << ieta << " " << iphi
556  << ", Energy em/had/long/short = "
557  << energyemHF << " " << energyhadHF << " "
558  << theLongHitEnergy << " " << theShortHitEnergy << " "
559  << ". The status value is " << theStatusValue
560  << ". Time = " << theLongHit->time() << " "
561  << ". The short and long flags : "
562  << flagShortDPG << " " << flagLongDPG << " "
563  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
564  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
565  << ". Long fibres were cleaned." << std::endl;
566  */
567  longFibre -= theLongHitEnergy;
568  theLongHitEnergy = 0.;
569  }
570  }
571 
572  // Special treatment for tower 29
573  // A tower with energy only at ieta = +/- 29 is not physical -> Clean
574  if ( abs(ieta) == 29 ) {
575  // rescaleFactor = 0. ;
576  // Clean long fibres
577  if ( theLongHitEnergy > shortFibre_Cut/2. ) {
578  pfrhHFEMCleaned29 = createHcalRecHit( detid,
579  theLongHitEnergy,
580  PFLayer::HF_EM,
581  hcalEndcapGeometry,
582  ct.id().rawId() );
583  pfrhHFEMCleaned29->setRescale(theLongHit->time());
584  /*
585  std::cout << "ieta/iphi = " << ieta << " " << iphi
586  << ", Energy em/had/long/short = "
587  << energyemHF << " " << energyhadHF << " "
588  << theLongHitEnergy << " " << theShortHitEnergy << " "
589  << ". Long fibres were cleaned." << std::endl;
590  */
591  longFibre -= theLongHitEnergy;
592  theLongHitEnergy = 0.;
593  }
594  // Clean short fibres
595  if ( theShortHitEnergy > shortFibre_Cut/2. ) {
596  pfrhHFHADCleaned29 = createHcalRecHit( detid,
597  theShortHitEnergy,
598  PFLayer::HF_HAD,
599  hcalEndcapGeometry,
600  ct.id().rawId() );
601  pfrhHFHADCleaned29->setRescale(theShortHit->time());
602  /*
603  std::cout << "ieta/iphi = " << ieta << " " << iphi
604  << ", Energy em/had/long/short = "
605  << energyemHF << " " << energyhadHF << " "
606  << theLongHitEnergy << " " << theShortHitEnergy << " "
607  << ". Short fibres were cleaned." << std::endl;
608  */
609  shortFibre -= theShortHitEnergy;
610  theShortHitEnergy = 0.;
611  }
612  }
613  // Check the timing of the long and short fibre rechits
614 
615  // First, check the timing of long and short fibre in eta = 29 if tower 30.
616  else if ( abs(ieta) == 30 ) {
617  int ieta29 = ieta > 0 ? 29 : -29;
618  HcalDetId theLongDetId29 (HcalForward, ieta29, iphi, 1);
619  HcalDetId theShortDetId29 (HcalForward, ieta29, iphi, 2);
620  iHF theLongHit29 = hfHandle->find(theLongDetId29);
621  iHF theShortHit29 = hfHandle->find(theShortDetId29);
622  //
623  double theLongHitEnergy29 = 0.;
624  double theShortHitEnergy29 = 0.;
625  bool flagShortDPG29 = false;
626  bool flagLongDPG29 = false;
627  bool flagShortTimeDPG29 = false;
628  bool flagLongTimeDPG29 = false;
629  bool flagShortPulseDPG29 = false;
630  bool flagLongPulseDPG29 = false;
631  //
632  if ( theLongHit29 != hfHandle->end() ) {
633  theLongHitEnergy29 = theLongHit29->energy() ;
634  flagLongDPG29 = applyLongShortDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
635  flagLongTimeDPG29 = applyTimeDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
636  flagLongPulseDPG29 = applyPulseDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
637  }
638  //
639  if ( theShortHit29 != hfHandle->end() ) {
640  theShortHitEnergy29 = theShortHit29->energy();
641  flagShortDPG29 = applyLongShortDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
642  flagShortTimeDPG29 = applyTimeDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
643  flagShortPulseDPG29 = applyPulseDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
644  }
645 
646  if ( theLongHitEnergy29 > longShortFibre_Cut &&
647  ( theLongHit29->time() < minLongTiming_Cut ||
648  theLongHit29->time() > maxLongTiming_Cut ||
649  flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
650  //rescaleFactor = 0. ;
651  pfrhHFEMCleaned29 = createHcalRecHit( detid,
652  theLongHitEnergy29,
653  PFLayer::HF_EM,
654  hcalEndcapGeometry,
655  ct.id().rawId() );
656  pfrhHFEMCleaned29->setRescale(theLongHit29->time());
657  /*
658  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
659  << ", Energy em/had/long/short = "
660  << energyemHF << " " << energyhadHF << " "
661  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
662  << ". Time = " << theLongHit29->time() << " "
663  << ". The short and long flags : "
664  << flagShortDPG29 << " " << flagLongDPG29 << " "
665  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
666  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
667  << ". Long fibres were cleaned." << std::endl;
668  */
669  longFibre -= theLongHitEnergy29;
670  theLongHitEnergy29 = 0;
671  }
672 
673  if ( theShortHitEnergy29 > longShortFibre_Cut &&
674  ( theShortHit29->time() < minShortTiming_Cut ||
675  theShortHit29->time() > maxShortTiming_Cut ||
676  flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
677  //rescaleFactor = 0. ;
678  pfrhHFHADCleaned29 = createHcalRecHit( detid,
679  theShortHitEnergy29,
680  PFLayer::HF_HAD,
681  hcalEndcapGeometry,
682  ct.id().rawId() );
683  pfrhHFHADCleaned29->setRescale(theShortHit29->time());
684  /*
685  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
686  << ", Energy em/had/long/short = "
687  << energyemHF << " " << energyhadHF << " "
688  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
689  << ". Time = " << theShortHit29->time() << " "
690  << ". The short and long flags : "
691  << flagShortDPG29 << " " << flagLongDPG29 << " "
692  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
693  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
694  << ". Short fibres were cleaned." << std::endl;
695  */
696  shortFibre -= theShortHitEnergy29;
697  theShortHitEnergy29 = 0.;
698  }
699 
700  // Some energy must be in the long fibres is there is some energy in the short fibres !
701  if ( theShortHitEnergy29 > shortFibre_Cut &&
702  ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
703  flagShortDPG29 ) ) {
704  // Check if the long-fibre hit was not cleaned already (because hot)
705  // In this case don't apply the cleaning
706  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId29);
707  unsigned theStatusValue = theStatus->getValue();
708  // The channel is killed
709  if ( !theStatusValue ) {
710  //rescaleFactor = 0. ;
711  pfrhHFHADCleaned29 = createHcalRecHit( detid,
712  theShortHitEnergy29,
713  PFLayer::HF_HAD,
714  hcalEndcapGeometry,
715  ct.id().rawId() );
716  pfrhHFHADCleaned29->setRescale(theShortHit29->time());
717  /*
718  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
719  << ", Energy em/had/long/short = "
720  << energyemHF << " " << energyhadHF << " "
721  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
722  << ". Time = " << theShortHit29->time() << " "
723  << ". The status value is " << theStatusValue
724  << ". The short and long flags : "
725  << flagShortDPG29 << " " << flagLongDPG29 << " "
726  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
727  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
728  << ". Short fibres were cleaned." << std::endl;
729  */
730  shortFibre -= theShortHitEnergy29;
731  theShortHitEnergy29 = 0.;
732  }
733  }
734 
735  // Some energy must be in the short fibres is there is some energy in the long fibres !
736  if ( theLongHitEnergy29 > longFibre_Cut &&
737  ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
738  flagLongDPG29 ) ) {
739  // Check if the long-fibre hit was not cleaned already (because hot)
740  // In this case don't apply the cleaning
741  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
742  unsigned theStatusValue = theStatus->getValue();
743  // The channel is killed
744  if ( !theStatusValue ) {
745  //rescaleFactor = 0. ;
746  pfrhHFEMCleaned29 = createHcalRecHit( detid,
747  theLongHitEnergy29,
748  PFLayer::HF_EM,
749  hcalEndcapGeometry,
750  ct.id().rawId() );
751  pfrhHFEMCleaned29->setRescale(theLongHit29->time());
752  /*
753  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
754  << ", Energy em/had/long/short = "
755  << energyemHF << " " << energyhadHF << " "
756  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
757  << ". The status value is " << theStatusValue
758  << ". Time = " << theLongHit29->time() << " "
759  << ". The short and long flags : "
760  << flagShortDPG29 << " " << flagLongDPG29 << " "
761  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
762  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
763  << ". Long fibres were cleaned." << std::endl;
764  */
765  longFibre -= theLongHitEnergy29;
766  theLongHitEnergy29 = 0.;
767  }
768  }
769 
770  // Check that the energy in tower 29 is smaller than in tower 30
771  // First in long fibres
772  if ( theLongHitEnergy29 > std::max(theLongHitEnergy,shortFibre_Cut/2) ) {
773  pfrhHFEMCleaned29 = createHcalRecHit( detid,
774  theLongHitEnergy29,
775  PFLayer::HF_EM,
776  hcalEndcapGeometry,
777  ct.id().rawId() );
778  pfrhHFEMCleaned29->setRescale(theLongHit29->time());
779  /*
780  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
781  << ", Energy L29/S29/L30/S30 = "
782  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
783  << theLongHitEnergy << " " << theShortHitEnergy << " "
784  << ". Long fibres were cleaned." << std::endl;
785  */
786  longFibre -= theLongHitEnergy29;
787  theLongHitEnergy29 = 0.;
788  }
789  // Second in short fibres
790  if ( theShortHitEnergy29 > std::max(theShortHitEnergy,shortFibre_Cut/2.) ) {
791  pfrhHFHADCleaned29 = createHcalRecHit( detid,
792  theShortHitEnergy29,
793  PFLayer::HF_HAD,
794  hcalEndcapGeometry,
795  ct.id().rawId() );
796  pfrhHFHADCleaned29->setRescale(theShortHit29->time());
797  /*
798  std::cout << "ieta/iphi = " << ieta << " " << iphi
799  << ", Energy L29/S29/L30/S30 = "
800  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
801  << theLongHitEnergy << " " << theShortHitEnergy << " "
802  << ". Short fibres were cleaned." << std::endl;
803  */
804  shortFibre -= theShortHitEnergy29;
805  theShortHitEnergy29 = 0.;
806  }
807  }
808 
809 
810  // Determine EM and HAD after cleaning of short and long fibres
811  energyhadHF = 2.*shortFibre;
812  energyemHF = longFibre - shortFibre;
813 
814  // The EM energy might be negative, as it amounts to Long - Short
815  // In that case, put the EM "energy" in the HAD energy
816  // Just to avoid systematic positive bias due to "Short" high fluctuations
817  if ( energyemHF < thresh_HF_ ) {
818  energyhadHF += energyemHF;
819  energyemHF = 0.;
820  }
821 
822  // Apply HCAL DPG calibration factors, if requested
823  if ( HF_Calib_ ) {
824  energyhadHF *= myPFCorr->getValues(detid)->getValue();
825  energyemHF *= myPFCorr->getValues(detid)->getValue();
826  }
827 
828 
829  // Create an EM and a HAD rechit if above threshold.
830  if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) {
831  pfrhHFEM = createHcalRecHit( detid,
832  energyemHF,
833  PFLayer::HF_EM,
834  hcalEndcapGeometry,
835  ct.id().rawId() );
836  pfrhHFHAD = createHcalRecHit( detid,
837  energyhadHF,
838  PFLayer::HF_HAD,
839  hcalEndcapGeometry,
840  ct.id().rawId() );
841  pfrhHFEM->setEnergyUp(energyhadHF);
842  pfrhHFHAD->setEnergyUp(energyemHF);
843  }
844 
845  }
846  break;
847  default:
848  LogError("PFRecHitProducerHCAL")
849  <<"CaloTower constituent: unknown layer : "
850  <<detid.subdet()<<endl;
851  }
852 
853  if(pfrh) {
854  rechits.push_back( *pfrh );
855  delete pfrh;
856  idSortedRecHits.insert( make_pair(ct.id().rawId(),
857  rechits.size()-1 ) );
858  }
859  if(pfrhCleaned) {
860  rechitsCleaned.push_back( *pfrhCleaned );
861  delete pfrhCleaned;
862  }
863  //---ab: 2 rechits for HF:
864  if(pfrhHFEM) {
865  HFEMRecHits->push_back( *pfrhHFEM );
866  delete pfrhHFEM;
867  idSortedRecHitsHFEM.insert( make_pair(ct.id().rawId(),
868  HFEMRecHits->size()-1 ) );
869  }
870  if(pfrhHFHAD) {
871  HFHADRecHits->push_back( *pfrhHFHAD );
872  delete pfrhHFHAD;
873  idSortedRecHitsHFHAD.insert( make_pair(ct.id().rawId(),
874  HFHADRecHits->size()-1 ) );
875  }
876  //---ab
877  if(pfrhHFEMCleaned) {
878  rechitsCleaned.push_back( *pfrhHFEMCleaned );
879  delete pfrhHFEMCleaned;
880  }
881  if(pfrhHFHADCleaned) {
882  rechitsCleaned.push_back( *pfrhHFHADCleaned );
883  delete pfrhHFHADCleaned;
884  }
885  if(pfrhHFEMCleaned29) {
886  rechitsCleaned.push_back( *pfrhHFEMCleaned29 );
887  delete pfrhHFEMCleaned29;
888  }
889  if(pfrhHFHADCleaned29) {
890  rechitsCleaned.push_back( *pfrhHFHADCleaned29 );
891  delete pfrhHFHADCleaned29;
892  }
893  }
894  }
895  // do navigation
896  for(unsigned i=0; i<rechits.size(); i++ ) {
897  findRecHitNeighboursCT( rechits[i],
898  idSortedRecHits,
899  caloTowerTopology);
900  }
901  for(unsigned i=0; i<HFEMRecHits->size(); i++ ) {
902  findRecHitNeighboursCT( (*HFEMRecHits)[i],
903  idSortedRecHitsHFEM,
904  caloTowerTopology);
905  }
906  for(unsigned i=0; i<HFHADRecHits->size(); i++ ) {
907  findRecHitNeighboursCT( (*HFHADRecHits)[i],
908  idSortedRecHitsHFHAD,
909  caloTowerTopology);
910  }
911  iEvent.put( HFHADRecHits,"HFHAD" );
912  iEvent.put( HFEMRecHits,"HFEM" );
913  }
914  }
915  else if( !(inputTagHcalRecHitsHBHE_ == InputTag()) ) {
916  // clustering is not done on CaloTowers but on HCAL rechits.
917 
918 
919  // get the hcal topology
920  HcalTopology hcalTopology;
921 
922  // HCAL rechits
923  // vector<edm::Handle<HBHERecHitCollection> > hcalHandles;
924  edm::Handle<HBHERecHitCollection> hcalHandle;
925 
926 
927  bool found = iEvent.getByLabel(inputTagHcalRecHitsHBHE_,
928  hcalHandle );
929 
930  if(!found) {
931  ostringstream err;
932  err<<"could not find rechits "<<inputTagHcalRecHitsHBHE_;
933  LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
934 
935  throw cms::Exception( "MissingProduct", err.str());
936  }
937  else {
938  assert( hcalHandle.isValid() );
939 
940  const edm::Handle<HBHERecHitCollection>& handle = hcalHandle;
941  for(unsigned irechit=0; irechit<handle->size(); irechit++) {
942  const HBHERecHit& hit = (*handle)[irechit];
943 
944  double energy = hit.energy();
945 
946  reco::PFRecHit* pfrh = 0;
947 
948 
949  const HcalDetId& detid = hit.detid();
950  switch( detid.subdet() ) {
951  case HcalBarrel:
952  {
953  if(energy < thresh_Barrel_ ) continue;
954  pfrh = createHcalRecHit( detid,
955  energy,
956  PFLayer::HCAL_BARREL1,
957  hcalBarrelGeometry );
958  }
959  break;
960  case HcalEndcap:
961  {
962  if(energy < thresh_Endcap_ ) continue;
963  pfrh = createHcalRecHit( detid,
964  energy,
965  PFLayer::HCAL_ENDCAP,
966  hcalEndcapGeometry );
967  }
968  break;
969  case HcalForward:
970  {
971  if(energy < thresh_HF_ ) continue;
972  pfrh = createHcalRecHit( detid,
973  energy,
974  PFLayer::HF_HAD,
975  hcalEndcapGeometry );
976  }
977  break;
978  default:
979  LogError("PFRecHitProducerHCAL")
980  <<"HCAL rechit: unknown layer : "<<detid.subdet()<<endl;
981  continue;
982  }
983 
984  if(pfrh) {
985  rechits.push_back( *pfrh );
986  delete pfrh;
987  idSortedRecHits.insert( make_pair(detid.rawId(),
988  rechits.size()-1 ) );
989  }
990  }
991 
992 
993  // do navigation:
994  for(unsigned i=0; i<rechits.size(); i++ ) {
995 
996  findRecHitNeighbours( rechits[i], idSortedRecHits,
997  hcalTopology,
998  *hcalBarrelGeometry,
999  hcalTopology,
1000  *hcalEndcapGeometry);
1001  } // loop for navigation
1002  } // endif hcal rechits were found
1003  } // endif clustering on rechits in hcal
1004 }
1005 
1006 
1007 
1008 
1009 
1010 
1011 reco::PFRecHit*
1012 PFRecHitProducerHCAL::createHcalRecHit( const DetId& detid,
1013  double energy,
1014  PFLayer::Layer layer,
1015  const CaloSubdetectorGeometry* geom,
1016  unsigned newDetId ) {
1017 
1018  const CaloCellGeometry *thisCell = geom->getGeometry(detid);
1019  if(!thisCell) {
1020  edm::LogError("PFRecHitProducerHCAL")
1021  <<"warning detid "<<detid.rawId()<<" not found in layer "
1022  <<layer<<endl;
1023  return 0;
1024  }
1025 
1026  const GlobalPoint& position = thisCell->getPosition();
1027 
1028  double depth_correction = 0.;
1029  switch ( layer ) {
1030  case PFLayer::HF_EM:
1031  depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
1032  break;
1033  case PFLayer::HF_HAD:
1034  depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
1035  break;
1036  default:
1037  break;
1038  }
1039 
1040  unsigned id = detid;
1041  if(newDetId) id = newDetId;
1042  reco::PFRecHit *rh =
1043  new reco::PFRecHit( id, layer, energy,
1044  position.x(), position.y(), position.z()+depth_correction,
1045  0,0,0 );
1046 
1047 
1048 
1049 
1050  // set the corners
1051  const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
1052 
1053  assert( corners.size() == 8 );
1054 
1055  rh->setNECorner( corners[0].x(), corners[0].y(), corners[0].z()+depth_correction );
1056  rh->setSECorner( corners[1].x(), corners[1].y(), corners[1].z()+depth_correction );
1057  rh->setSWCorner( corners[2].x(), corners[2].y(), corners[2].z()+depth_correction );
1058  rh->setNWCorner( corners[3].x(), corners[3].y(), corners[3].z()+depth_correction );
1059 
1060  return rh;
1061 }
1062 
1063 
1064 
1065 
1066 void
1067 PFRecHitProducerHCAL::findRecHitNeighbours
1068 ( reco::PFRecHit& rh,
1069  const map<unsigned,unsigned >& sortedHits,
1070  const CaloSubdetectorTopology& barrelTopology,
1071  const CaloSubdetectorGeometry& barrelGeometry,
1072  const CaloSubdetectorTopology& endcapTopology,
1073  const CaloSubdetectorGeometry& endcapGeometry ) {
1074 
1075  //cout<<"------PFRecHitProducerHcaL:findRecHitNeighbours navigation value "<<navigation_HF_<<endl;
1076  if(navigation_HF_ == false){
1077  if( rh.layer() == PFLayer::HF_HAD )
1078  return;
1079  if( rh.layer() == PFLayer::HF_EM )
1080  return;
1081  }
1082  DetId detid( rh.detId() );
1083 
1084  const CaloSubdetectorTopology* topology = 0;
1085  const CaloSubdetectorGeometry* geometry = 0;
1086  const CaloSubdetectorGeometry* othergeometry = 0;
1087 
1088  switch( rh.layer() ) {
1089  case PFLayer::ECAL_ENDCAP:
1090  topology = &endcapTopology;
1091  geometry = &endcapGeometry;
1092  break;
1093  case PFLayer::ECAL_BARREL:
1094  topology = &barrelTopology;
1095  geometry = &barrelGeometry;
1096  break;
1097  case PFLayer::HCAL_ENDCAP:
1098  topology = &endcapTopology;
1099  geometry = &endcapGeometry;
1100  othergeometry = &barrelGeometry;
1101  break;
1102  case PFLayer::HCAL_BARREL1:
1103  topology = &barrelTopology;
1104  geometry = &barrelGeometry;
1105  othergeometry = &endcapGeometry;
1106  break;
1107  case PFLayer::PS1:
1108  case PFLayer::PS2:
1109  topology = &barrelTopology;
1110  geometry = &barrelGeometry;
1111  othergeometry = &endcapGeometry;
1112  break;
1113  default:
1114  assert(0);
1115  }
1116 
1117  assert( topology && geometry );
1118 
1119  CaloNavigator<DetId> navigator(detid, topology);
1120 
1121  DetId north = navigator.north();
1122 
1123  DetId northeast(0);
1124  if( north != DetId(0) ) {
1125  northeast = navigator.east();
1126  }
1127  navigator.home();
1128 
1129 
1130  DetId south = navigator.south();
1131 
1132 
1133 
1134  DetId southwest(0);
1135  if( south != DetId(0) ) {
1136  southwest = navigator.west();
1137  }
1138  navigator.home();
1139 
1140 
1141  DetId east = navigator.east();
1142  DetId southeast;
1143  if( east != DetId(0) ) {
1144  southeast = navigator.south();
1145  }
1146  navigator.home();
1147  DetId west = navigator.west();
1148  DetId northwest;
1149  if( west != DetId(0) ) {
1150  northwest = navigator.north();
1151  }
1152  navigator.home();
1153 
1154  IDH i = sortedHits.find( north.rawId() );
1155  if(i != sortedHits.end() )
1156  rh.add4Neighbour( i->second );
1157 
1158  i = sortedHits.find( northeast.rawId() );
1159  if(i != sortedHits.end() )
1160  rh.add8Neighbour( i->second );
1161 
1162  i = sortedHits.find( south.rawId() );
1163  if(i != sortedHits.end() )
1164  rh.add4Neighbour( i->second );
1165 
1166  i = sortedHits.find( southwest.rawId() );
1167  if(i != sortedHits.end() )
1168  rh.add8Neighbour( i->second );
1169 
1170  i = sortedHits.find( east.rawId() );
1171  if(i != sortedHits.end() )
1172  rh.add4Neighbour( i->second );
1173 
1174  i = sortedHits.find( southeast.rawId() );
1175  if(i != sortedHits.end() )
1176  rh.add8Neighbour( i->second );
1177 
1178  i = sortedHits.find( west.rawId() );
1179  if(i != sortedHits.end() )
1180  rh.add4Neighbour( i->second );
1181 
1182  i = sortedHits.find( northwest.rawId() );
1183  if(i != sortedHits.end() )
1184  rh.add8Neighbour( i->second );
1185 
1186 
1187 }
1188 
1189 
1190 void
1191 PFRecHitProducerHCAL::findRecHitNeighboursCT
1192 ( reco::PFRecHit& rh,
1193  const map<unsigned, unsigned >& sortedHits,
1194  const CaloSubdetectorTopology& topology ) {
1195  //cout<<"------PFRecHitProducerHcaL:findRecHitNeighboursCT navigation value "<<navigation_HF_<<endl;
1196  // cout<<"----------- rechit print out"<<endl;
1197  // if(( rh.layer() == PFLayer::HF_HAD )||(rh.layer() == PFLayer::HF_EM)) {
1198 
1199  // cout<<rh<<endl;
1200  // }
1201  if(navigation_HF_ == false){
1202  if( rh.layer() == PFLayer::HF_HAD )
1203  return;
1204  if( rh.layer() == PFLayer::HF_EM )
1205  return;
1206  }
1207  CaloTowerDetId ctDetId( rh.detId() );
1208 
1209 
1210  vector<DetId> northids = topology.north(ctDetId);
1211  vector<DetId> westids = topology.west(ctDetId);
1212  vector<DetId> southids = topology.south(ctDetId);
1213  vector<DetId> eastids = topology.east(ctDetId);
1214 
1215 
1216  CaloTowerDetId badId;
1217 
1218  // all the following detids will be CaloTowerDetId
1219  CaloTowerDetId north;
1220  CaloTowerDetId northwest;
1221  CaloTowerDetId northwest2;
1222  CaloTowerDetId west;
1223  CaloTowerDetId west2;
1224  CaloTowerDetId southwest;
1225  CaloTowerDetId southwest2;
1226  CaloTowerDetId south;
1227  CaloTowerDetId southeast;
1228  CaloTowerDetId southeast2;
1229  CaloTowerDetId east;
1230  CaloTowerDetId east2;
1231  CaloTowerDetId northeast;
1232  CaloTowerDetId northeast2;
1233 
1234  // for north and south, there is no ambiguity : 1 or 0 neighbours
1235  stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours ");
1236 
1237  switch( northids.size() ) {
1238  case 0:
1239  break;
1240  case 1:
1241  north = northids[0];
1242  break;
1243  default:
1244  err<<"north: "<<northids.size();
1245  throw( err.str() );
1246  }
1247 
1248  switch( southids.size() ) {
1249  case 0:
1250  break;
1251  case 1:
1252  south = southids[0];
1253  break;
1254  default:
1255  err<<"south: "<<southids.size();
1256  throw( err.str() );
1257  }
1258 
1259  // for east and west, one must take care
1260  // of the pitch change in HCAL endcap.
1261 
1262  switch( eastids.size() ) {
1263  case 0:
1264  break;
1265  case 1:
1266  east = eastids[0];
1267  northeast = getNorth(east, topology);
1268  southeast = getSouth(east, topology);
1269  break;
1270  case 2:
1271  // in this case, 0 is more on the north than 1
1272  east = eastids[0];
1273  east2 = eastids[1];
1274  northeast = getNorth(east, topology );
1275  southeast = getSouth(east2, topology);
1276  northeast2 = getNorth(northeast, topology );
1277  southeast2 = getSouth(southeast, topology);
1278  break;
1279  default:
1280  err<<"eastids: "<<eastids.size();
1281  throw( err.str() );
1282  }
1283 
1284 
1285  switch( westids.size() ) {
1286  case 0:
1287  break;
1288  case 1:
1289  west = westids[0];
1290  northwest = getNorth(west, topology);
1291  southwest = getSouth(west, topology);
1292  break;
1293  case 2:
1294  // in this case, 0 is more on the north than 1
1295  west = westids[0];
1296  west2 = westids[1];
1297  northwest = getNorth(west, topology );
1298  southwest = getSouth(west2, topology );
1299  northwest2 = getNorth(northwest, topology );
1300  southwest2 = getSouth(southwest, topology );
1301  break;
1302  default:
1303  err<<"westids: "<< westids.size();
1304  throw( err.str() );
1305  }
1306 
1307 
1308 
1309 
1310  // find and set neighbours
1311 
1312  IDH i = sortedHits.find( north.rawId() );
1313  if(i != sortedHits.end() )
1314  rh.add4Neighbour( i->second );
1315 
1316  i = sortedHits.find( northeast.rawId() );
1317  if(i != sortedHits.end() )
1318  rh.add8Neighbour( i->second );
1319 
1320  i = sortedHits.find( northeast2.rawId() );
1321  if(i != sortedHits.end() )
1322  rh.add8Neighbour( i->second );
1323 
1324  i = sortedHits.find( south.rawId() );
1325  if(i != sortedHits.end() )
1326  rh.add4Neighbour( i->second );
1327 
1328  i = sortedHits.find( southwest.rawId() );
1329  if(i != sortedHits.end() )
1330  rh.add8Neighbour( i->second );
1331 
1332  i = sortedHits.find( southwest2.rawId() );
1333  if(i != sortedHits.end() )
1334  rh.add8Neighbour( i->second );
1335 
1336  i = sortedHits.find( east.rawId() );
1337  if(i != sortedHits.end() )
1338  rh.add4Neighbour( i->second );
1339 
1340  i = sortedHits.find( east2.rawId() );
1341  if(i != sortedHits.end() )
1342  rh.add4Neighbour( i->second );
1343 
1344  i = sortedHits.find( southeast.rawId() );
1345  if(i != sortedHits.end() )
1346  rh.add8Neighbour( i->second );
1347 
1348  i = sortedHits.find( southeast2.rawId() );
1349  if(i != sortedHits.end() )
1350  rh.add8Neighbour( i->second );
1351 
1352  i = sortedHits.find( west.rawId() );
1353  if(i != sortedHits.end() )
1354  rh.add4Neighbour( i->second );
1355 
1356  i = sortedHits.find( west2.rawId() );
1357  if(i != sortedHits.end() )
1358  rh.add4Neighbour( i->second );
1359 
1360  i = sortedHits.find( northwest.rawId() );
1361  if(i != sortedHits.end() )
1362  rh.add8Neighbour( i->second );
1363 
1364  i = sortedHits.find( northwest2.rawId() );
1365  if(i != sortedHits.end() )
1366  rh.add8Neighbour( i->second );
1367 
1368  // cout<<"----------- rechit print out"<<endl;
1369  // if(( rh.layer() == PFLayer::HF_HAD )||(rh.layer() == PFLayer::HF_EM)) {
1370 
1371  // cout<<rh<<endl;
1372  // }
1373 }
1374 
1375 
1376 
1377 DetId
1378 PFRecHitProducerHCAL::getSouth(const DetId& id,
1379  const CaloSubdetectorTopology& topology) {
1380 
1381  DetId south;
1382  vector<DetId> sids = topology.south(id);
1383  if(sids.size() == 1)
1384  south = sids[0];
1385 
1386  return south;
1387 }
1388 
1389 
1390 
1391 DetId
1392 PFRecHitProducerHCAL::getNorth(const DetId& id,
1393  const CaloSubdetectorTopology& topology) {
1394 
1395  DetId north;
1396  vector<DetId> nids = topology.north(id);
1397  if(nids.size() == 1)
1398  north = nids[0];
1399 
1400  return north;
1401 }
1402 
1403 
reco::PFRecHit::setSECorner
void setSECorner(double posx, double posy, double posz)
Definition: PFRecHit.cc:226
ExpressReco_HICollisions_FallBack.e
tuple e
Definition: ExpressReco_HICollisions_FallBack.py:1068
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
reco::PFRecHit::add4Neighbour
void add4Neighbour(unsigned index)
Definition: PFRecHit.cc:176
i
int i
Definition: DBlmapReader.cc:9
PFRecHitProducer::myPFCorr
const HcalPFCorrs * myPFCorr
Definition: PFRecHitProducer.h:69
reco::PFRecHit::add8Neighbour
void add8Neighbour(unsigned index)
Definition: PFRecHit.cc:181
CaloTower::constituentsSize
size_t constituentsSize() const
Definition: CaloTower.h:74
reco::PFRecHit::setNECorner
void setNECorner(double posx, double posy, double posz)
Definition: PFRecHit.cc:231
HcalDetId::subdet
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:32
CaloNavigator::home
void home() const
move the navigator back to the starting point
Definition: CaloNavigator.h:143
CaloRecHit::detid
const DetId & detid() const
Definition: CaloRecHit.h:21
reco::PFRecHit::detId
unsigned detId() const
rechit detId
Definition: PFRecHit.h:99
CaloNavigator::west
virtual T west() const
move the navigator west
Definition: CaloNavigator.h:81
edm::SortedCollection< CaloTower >::const_iterator
std::vector< CaloTower >::const_iterator const_iterator
Definition: SortedCollection.h:84
PV3DBase::y
T y() const
Definition: PV3DBase.h:57
CaloSubdetectorTopology::west
virtual std::vector< DetId > west(const DetId &id) const =0
abs
#define abs(x)
Definition: mlp_lapack.h:159
min
#define min(a, b)
Definition: mlp_lapack.h:161
PFRecHitProducerHCAL::createHcalRecHit
reco::PFRecHit * createHcalRecHit(const DetId &detid, double energy, PFLayer::Layer layer, const CaloSubdetectorGeometry *geom, unsigned newDetId=0)
Definition: PFRecHitProducerHCAL.cc:1012
PFRecHitProducer::theEcalChStatus
const EcalChannelStatus * theEcalChStatus
Definition: PFRecHitProducer.h:71
CaloNavigator::north
virtual T north() const
move the navigator north
Definition: CaloNavigator.h:51
HcalPFCorr::getValue
float getValue() const
Definition: HcalPFCorr.h:18
position
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
CaloSubdetectorGeometry::getGeometry
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
Definition: CaloSubdetectorGeometry.cc:76
CaloTowerConstituentsMap::constituentsOf
std::vector< DetId > constituentsOf(const CaloTowerDetId &id) const
Get the constituent detids for this tower id ( not yet implemented )
Definition: CaloTowerConstituentsMap.cc:57
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
reco::PFRecHit::layer
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:102
reco::PFRecHit::setSWCorner
void setSWCorner(double posx, double posy, double posz)
Definition: PFRecHit.cc:221
PFRecHitProducerHCAL::inputTagCaloTowers_
edm::InputTag inputTagCaloTowers_
Definition: PFRecHitProducerHCAL.h:86
CaloNavigator::east
virtual T east() const
move the navigator east
Definition: CaloNavigator.h:71
iEvent
int iEvent
Definition: GenABIO.cc:243
reco::PFRecHit
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
PFRecHitProducer
Base producer for particle flow rechits (PFRecHit)
Definition: PFRecHitProducer.h:33
z
Definition: DDAxes.h:10
ExpressReco_HICollisions_FallBack.caloTowers
tuple caloTowers
Definition: ExpressReco_HICollisions_FallBack.py:7904
CaloTower::emEnergy
double emEnergy() const
Definition: CaloTower.h:79
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
ExpressReco_HICollisions_FallBack.y
tuple y
Definition: ExpressReco_HICollisions_FallBack.py:4645
edm::Event::put
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
CaloRecHit::energy
float energy() const
Definition: CaloRecHit.h:19
reco::PFRecHit::setEnergyUp
void setEnergyUp(double eUp)
Definition: PFRecHit.h:73
PFRecHitProducerHCAL::inputTagHcalRecHitsHF_
edm::InputTag inputTagHcalRecHitsHF_
Definition: PFRecHitProducerHCAL.h:85
PV3DBase::z
T z() const
Definition: PV3DBase.h:58
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
patZpeak.handle
tuple handle
Definition: patZpeak.py:22
j
int j
Definition: DBlmapReader.cc:9
CaloTower::constituents
const std::vector< DetId > & constituents() const
Definition: CaloTower.h:73
PFRecHitProducerHCAL::getSouth
DetId getSouth(const DetId &id, const CaloSubdetectorTopology &topology)
Definition: PFRecHitProducerHCAL.cc:1378
reco::PFRecHit::setRescale
void setRescale(double factor)
Definition: PFRecHit.h:74
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:76
reco::PFRecHit::setNWCorner
void setNWCorner(double posx, double posy, double posz)
search for pointers to neighbours, using neighbours&#39; DetId.
Definition: PFRecHit.cc:216
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
PFRecHitProducerHCAL::inputTagHcalRecHitsHBHE_
edm::InputTag inputTagHcalRecHitsHBHE_
Definition: PFRecHitProducerHCAL.h:84
PFRecHitProducerHCAL::findRecHitNeighbours
void findRecHitNeighbours(reco::PFRecHit &rh, const std::map< unsigned, unsigned > &sortedHits, const CaloSubdetectorTopology &barrelTopo, const CaloSubdetectorGeometry &barrelGeom, const CaloSubdetectorTopology &endcapTopo, const CaloSubdetectorGeometry &endcapGeom)
Definition: PFRecHitProducerHCAL.cc:1068
PFRecHitProducerHCAL::findRecHitNeighboursCT
void findRecHitNeighboursCT(reco::PFRecHit &rh, const std::map< unsigned, unsigned > &sortedHits, const CaloSubdetectorTopology &topology)
Definition: PFRecHitProducerHCAL.cc:1192
CaloTower::hadEnergy
double hadEnergy() const
Definition: CaloTower.h:80
CaloSubdetectorTopology::north
virtual std::vector< DetId > north(const DetId &id) const =0
PFLayer::Layer
Layer
layer definition
Definition: PFLayer.h:31
HcalDetId::iphi
int iphi() const
get the cell iphi
Definition: HcalDetId.h:40
DetId
Definition: DetId.h:20
CaloTower::id
CaloTowerDetId id() const
Definition: CaloTower.h:72
PFRecHitProducer::theHcalChStatus
const HcalChannelQuality * theHcalChStatus
Definition: PFRecHitProducer.h:70
PFRecHitProducerHCAL::createRecHits
void createRecHits(std::vector< reco::PFRecHit > &rechits, std::vector< reco::PFRecHit > &rechitsCleaned, edm::Event &, const edm::EventSetup &)
Definition: PFRecHitProducerHCAL.cc:113
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:15
PFRecHitProducer::theTowerConstituentsMap
const CaloTowerConstituentsMap * theTowerConstituentsMap
Definition: PFRecHitProducer.h:72
geometry
ESHandle< TrackerGeometry > geometry
Definition: TkLasBeamFitter.cc:178
ExpressReco_HICollisions_FallBack.x
tuple x
Definition: ExpressReco_HICollisions_FallBack.py:4646
PFRecHitProducer::IDH
std::map< unsigned, unsigned >::const_iterator IDH
Definition: PFRecHitProducer.h:48
CaloSubdetectorTopology::south
virtual std::vector< DetId > south(const DetId &id) const =0
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:58
PFRecHitProducer::thresh_Barrel_
double thresh_Barrel_
rechits with E &lt; threshold will not give rise to a PFRecHit
Definition: PFRecHitProducer.h:64
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:88
PFRecHitProducerHCAL::thresh_HF_
double thresh_HF_
threshold for HF
Definition: PFRecHitProducerHCAL.h:89
PFRecHitProducerHCAL::getNorth
DetId getNorth(const DetId &id, const CaloSubdetectorTopology &topology)
Definition: PFRecHitProducerHCAL.cc:1392
HcalChannelStatus::getValue
uint32_t getValue() const
Definition: HcalChannelStatus.h:60
CaloNavigator::south
virtual T south() const
move the navigator south
Definition: CaloNavigator.h:61
HcalCondObjectContainer::getValues
const Item * getValues(DetId fId) const
Definition: HcalCondObjectContainer.h:95
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Definition: CaloCellGeometry.h:46
PV3DBase::x
T x() const
Definition: PV3DBase.h:56
CaloSubdetectorTopology::east
virtual std::vector< DetId > east(const DetId &id) const =0
PFRecHitProducerHCAL::PFRecHitProducerHCAL
PFRecHitProducerHCAL(const edm::ParameterSet &)
Definition: PFRecHitProducerHCAL.cc:41
CaloCellGeometry::getCorners
virtual const CornersVec & getCorners() const =0
Definition: CaloCellGeometry.cc:7
Generated for CMSSW Reference Manual by doxygen 1.8.5