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