CMS 3D CMS Logo

PFCTRecHitProducer.cc
Go to the documentation of this file.
2 
3 #include <memory>
4 
6 
14 
18 
22 
30 
32 
33 
36 
40 
41 
42 using namespace std;
43 using namespace edm;
44 
45 
47 {
48  thresh_Barrel_ =
49  iConfig.getParameter<double>("thresh_Barrel");
50  thresh_Endcap_ =
51  iConfig.getParameter<double>("thresh_Endcap");
52 
53 
54  thresh_HF_ =
55  iConfig.getParameter<double>("thresh_HF");
56  navigation_HF_ =
57  iConfig.getParameter<bool>("navigation_HF");
58  weight_HFem_ =
59  iConfig.getParameter<double>("weight_HFem");
60  weight_HFhad_ =
61  iConfig.getParameter<double>("weight_HFhad");
62 
63  HCAL_Calib_ =
64  iConfig.getParameter<bool>("HCAL_Calib");
65  HF_Calib_ =
66  iConfig.getParameter<bool>("HF_Calib");
67  HCAL_Calib_29 =
68  iConfig.getParameter<double>("HCAL_Calib_29");
69  HF_Calib_29 =
70  iConfig.getParameter<double>("HF_Calib_29");
71 
72  shortFibre_Cut = iConfig.getParameter<double>("ShortFibre_Cut");
73  longFibre_Fraction = iConfig.getParameter<double>("LongFibre_Fraction");
74 
75  longFibre_Cut = iConfig.getParameter<double>("LongFibre_Cut");
76  shortFibre_Fraction = iConfig.getParameter<double>("ShortFibre_Fraction");
77 
78  applyLongShortDPG_ = iConfig.getParameter<bool>("ApplyLongShortDPG");
79 
80  longShortFibre_Cut = iConfig.getParameter<double>("LongShortFibre_Cut");
81  minShortTiming_Cut = iConfig.getParameter<double>("MinShortTiming_Cut");
82  maxShortTiming_Cut = iConfig.getParameter<double>("MaxShortTiming_Cut");
83  minLongTiming_Cut = iConfig.getParameter<double>("MinLongTiming_Cut");
84  maxLongTiming_Cut = iConfig.getParameter<double>("MaxLongTiming_Cut");
85 
86  applyTimeDPG_ = iConfig.getParameter<bool>("ApplyTimeDPG");
87  applyPulseDPG_ = iConfig.getParameter<bool>("ApplyPulseDPG");
88  HcalMaxAllowedHFLongShortSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFLongShortSev");
89  HcalMaxAllowedHFDigiTimeSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFDigiTimeSev");
90  HcalMaxAllowedHFInTimeWindowSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFInTimeWindowSev");
91  HcalMaxAllowedChannelStatusSev_ = iConfig.getParameter<int>("HcalMaxAllowedChannelStatusSev");
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  //Get integer values of individual HCAL HF flags
102  hcalHFLongShortFlagValue_=1<<HcalCaloFlagLabels::HFLongShort;
103  hcalHFDigiTimeFlagValue_=1<<HcalCaloFlagLabels::HFDigiTime;
104  hcalHFInTimeWindowFlagValue_=1<<HcalCaloFlagLabels::HFInTimeWindow;
105 
106  hcalToken_ = consumes<HBHERecHitCollection>(iConfig.getParameter<InputTag>("hcalRecHitsHBHE"));
107  hfToken_ = consumes<HFRecHitCollection>(iConfig.getParameter<InputTag>("hcalRecHitsHF"));
108  towersToken_ = consumes<CaloTowerCollection>(iConfig.getParameter<InputTag>("caloTowers"));
109 
110 
111  edm::ParameterSet navSet = iConfig.getParameter<edm::ParameterSet>("navigator");
112  navigator_ = std::unique_ptr<PFRecHitNavigatorBase>{PFRecHitNavigationFactory::get()->create(navSet.getParameter<std::string>("name"),navSet)};
113 
114  //--ab
115  produces<reco::PFRecHitCollection>();
116  produces<reco::PFRecHitCollection>("Cleaned");
117  produces<reco::PFRecHitCollection>("HFHAD").setBranchAlias("HFHADRecHits");
118  produces<reco::PFRecHitCollection>("HFEM").setBranchAlias("HFEMRecHits");
119  //--ab
120 
121 }
122 
123 
125  const edm::EventSetup& iSetup) {
126 
127  navigator_->beginEvent(iSetup);
128 
129  edm::ESHandle<CaloGeometry> geoHandle;
130  iSetup.get<CaloGeometryRecord>().get(geoHandle);
131 
132  // get the hcalBarrel geometry
133  const CaloSubdetectorGeometry *hcalBarrelGeometry =
135 
136  // get the endcap geometry
137  const CaloSubdetectorGeometry *hcalEndcapGeometry =
139 
140  // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined
141  edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
142  iSetup.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
143  const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
144 
145  // Get Hcal Topology
146  edm::ESHandle<HcalTopology> hcalTopology;
147  iSetup.get<HcalRecNumberingRecord>().get( hcalTopology );
148  theHcalTopology = hcalTopology.product();
149 
150  auto rechits = std::make_unique<std::vector<reco::PFRecHit>>();
151  auto rechitsCleaned = std::make_unique<std::vector<reco::PFRecHit>>();
152  auto HFHADRecHits = std::make_unique<std::vector<reco::PFRecHit>>();
153  auto HFEMRecHits = std::make_unique<std::vector<reco::PFRecHit>>();
154 
156  iEvent.getByToken(towersToken_,caloTowers);
158  iEvent.getByToken(hfToken_,hfHandle);
159 
161  iEvent.getByToken(hcalToken_,hbheHandle);
162  // create rechits
164 
165  for(auto const & ct : *caloTowers) {
166  // get the hadronic energy.
167 
168  // Mike: Just ask for the Hadronic part only now!
169  // Patrick : ARGH ! While this is ok for the HCAL, this is
170  // just wrong for the HF (in which em/had are artificially
171  // separated.
172  double energy = ct.hadEnergy();
173  //Auguste: Photons in HF have no hadEnergy in fastsim: -> all RecHit collections are empty with photons.
174  double energyEM = ct.emEnergy(); // For HF !
175  //so test the total energy to deal with the photons in HF:
176  if( (energy+energyEM) < 1e-9 ) continue;
177 
178  //get the constituents of the tower
179  const std::vector<DetId>& hits = ct.constituents();
180  const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id());
181  HcalDetId detid;
182  bool foundHCALConstituent = false;
183  //Loop on the calotower constituents and search for HCAL
184  double dead = 0.;
185  double alive = 0.;
186  for(auto hit : hits) {
187  if(hit.det()==DetId::Hcal) {
188  foundHCALConstituent = true;
189  detid = hit;
190  if (theHcalTopology->getMergePositionFlag() &&
191  detid.subdet() == HcalEndcap) {
192  detid = theHcalTopology->idFront(detid);
193  }
194  // An HCAL tower was found: Look for dead ECAL channels in the same CaloTower.
195  if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
196  for(auto allConstituent : allConstituents) {
197  if ( allConstituent.det()==DetId::Ecal ) {
198  alive += 1.;
199  auto chIt = theEcalChStatus->find(allConstituent);
200  unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
201  if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
202  }
203  }
204  }
205  // Protection: tower 29 in HF is merged with tower 30.
206  // Just take the position of tower 30 in that case.
207  if ( detid.subdet() == HcalForward && abs(detid.ieta()) == 29 ) continue;
208  break;
209  }
210  }
211 
212  // In case of dead ECAL channel, rescale the HCAL energy...
213  double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
214 
215  reco::PFRecHit* pfrh = nullptr;
216  reco::PFRecHit* pfrhCleaned = nullptr;
217  //---ab: need 2 rechits for the HF:
218  reco::PFRecHit* pfrhHFEM = nullptr;
219  reco::PFRecHit* pfrhHFHAD = nullptr;
220  reco::PFRecHit* pfrhHFEMCleaned = nullptr;
221  reco::PFRecHit* pfrhHFHADCleaned = nullptr;
222  reco::PFRecHit* pfrhHFEMCleaned29 = nullptr;
223  reco::PFRecHit* pfrhHFHADCleaned29 = nullptr;
224 
225  if(foundHCALConstituent)
226  {
227  // std::cout << ", new Energy = " << energy << std::endl;
228  switch( detid.subdet() ) {
229  case HcalBarrel:
230  {
231  if(energy < thresh_Barrel_ ) continue;
232  if ( rescaleFactor > 1. ) {
233  pfrhCleaned = createHcalRecHit( detid,
234  energy,
236  hcalBarrelGeometry,
237  ct.id() );
238  pfrhCleaned->setTime(rescaleFactor);
239  energy *= rescaleFactor;
240  }
241  pfrh = createHcalRecHit( detid,
242  energy,
244  hcalBarrelGeometry,
245  ct.id() );
246  pfrh->setTime(rescaleFactor);
247  }
248  break;
249  case HcalEndcap:
250  {
251  if(energy < thresh_Endcap_ ) continue;
252  // Apply tower 29 calibration
253  if ( HCAL_Calib_ && abs(detid.ieta()) == 29 ) energy *= HCAL_Calib_29;
254  if ( rescaleFactor > 1. ) {
255  pfrhCleaned = createHcalRecHit( detid,
256  energy,
258  hcalEndcapGeometry,
259  ct.id() );
260  pfrhCleaned->setTime(rescaleFactor);
261  energy *= rescaleFactor;
262  }
263  pfrh = createHcalRecHit( detid,
264  energy,
266  hcalEndcapGeometry,
267  ct.id() );
268  pfrh->setTime(rescaleFactor);
269  }
270  break;
271  case HcalOuter:
272  {
273  }
274  break;
275  case HcalForward:
276  {
277  //---ab: 2 rechits for HF:
278  //double energyemHF = weight_HFem_*ct.emEnergy();
279  //double energyhadHF = weight_HFhad_*ct.hadEnergy();
280  double energyemHF = weight_HFem_ * energyEM;
281  double energyhadHF = weight_HFhad_ * energy;
282  // Some energy in the tower !
283  if((energyemHF+energyhadHF) < thresh_HF_ ) continue;
284 
285  // Some cleaning in the HF
286  double longFibre = energyemHF + energyhadHF/2.;
287  double shortFibre = energyhadHF/2.;
288  int ieta = detid.ieta();
289  int iphi = detid.iphi();
290  HcalDetId theLongDetId (HcalForward, ieta, iphi, 1);
291  HcalDetId theShortDetId (HcalForward, ieta, iphi, 2);
293  auto theLongHit = hfHandle->find(theLongDetId);
294  auto theShortHit = hfHandle->find(theShortDetId);
295  //
296  double theLongHitEnergy = 0.;
297  double theShortHitEnergy = 0.;
298  bool flagShortDPG = false;
299  bool flagLongDPG = false;
300  bool flagShortTimeDPG = false;
301  bool flagLongTimeDPG = false;
302  bool flagShortPulseDPG = false;
303  bool flagLongPulseDPG = false;
304  //
305  if ( theLongHit != hfHandle->end() ) {
306  int theLongFlag = theLongHit->flags();
307  theLongHitEnergy = theLongHit->energy();
308  flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
309  flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
310  flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
311 
312  //flagLongDPG = applyLongShortDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
313  //flagLongTimeDPG = applyTimeDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
314  //flagLongPulseDPG = applyPulseDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
315  }
316  //
317  if ( theShortHit != hfHandle->end() ) {
318  int theShortFlag = theShortHit->flags();
319  theShortHitEnergy = theShortHit->energy();
320  flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
321  flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
322  flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
323  //flagShortDPG = applyLongShortDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
324  //flagShortTimeDPG = applyTimeDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
325  //flagShortPulseDPG = applyPulseDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
326  }
327 
328  // Then check the timing in short and long fibres in all other towers.
329  if ( theShortHitEnergy > longShortFibre_Cut &&
330  ( theShortHit->time() < minShortTiming_Cut ||
331  theShortHit->time() > maxShortTiming_Cut ||
332  flagShortTimeDPG || flagShortPulseDPG ) ) {
333  // rescaleFactor = 0. ;
334  pfrhHFHADCleaned = createHcalRecHit( detid,
335  theShortHitEnergy,
337  hcalEndcapGeometry,
338  ct.id() );
339  pfrhHFHADCleaned->setTime(theShortHit->time());
340  /*
341  std::cout << "ieta/iphi = " << ieta << " " << iphi
342  << ", Energy em/had/long/short = "
343  << energyemHF << " " << energyhadHF << " "
344  << theLongHitEnergy << " " << theShortHitEnergy << " "
345  << ". Time = " << theShortHit->time() << " "
346  << ". The short and long flags : "
347  << flagShortDPG << " " << flagLongDPG << " "
348  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
349  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
350  << ". Short fibres were cleaned." << std::endl;
351  */
352  shortFibre -= theShortHitEnergy;
353  theShortHitEnergy = 0.;
354  }
355 
356 
357  if ( theLongHitEnergy > longShortFibre_Cut &&
358  ( theLongHit->time() < minLongTiming_Cut ||
359  theLongHit->time() > maxLongTiming_Cut ||
360  flagLongTimeDPG || flagLongPulseDPG ) ) {
361  //rescaleFactor = 0. ;
362  pfrhHFEMCleaned = createHcalRecHit( detid,
363  theLongHitEnergy,
365  hcalEndcapGeometry,
366  ct.id());
367  pfrhHFEMCleaned->setTime(theLongHit->time());
368  /*
369  std::cout << "ieta/iphi = " << ieta << " " << iphi
370  << ", Energy em/had/long/short = "
371  << energyemHF << " " << energyhadHF << " "
372  << theLongHitEnergy << " " << theShortHitEnergy << " "
373  << ". Time = " << theLongHit->time() << " "
374  << ". The short and long flags : "
375  << flagShortDPG << " " << flagLongDPG << " "
376  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
377  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
378  << ". Long fibres were cleaned." << std::endl;
379  */
380  longFibre -= theLongHitEnergy;
381  theLongHitEnergy = 0.;
382  }
383 
384  // Some energy must be in the long fibres is there is some energy in the short fibres !
385  if ( theShortHitEnergy > shortFibre_Cut &&
386  ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
387  flagShortDPG ) ) {
388  // Check if the long-fibre hit was not cleaned already (because hot)
389  // In this case don't apply the cleaning
390  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId);
391  unsigned theStatusValue = theStatus->getValue();
392  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
393  // The channel is killed
395  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
396  // rescaleFactor = 0. ;
397  pfrhHFHADCleaned = createHcalRecHit( detid,
398  theShortHitEnergy,
400  hcalEndcapGeometry,
401  ct.id() );
402  pfrhHFHADCleaned->setTime(theShortHit->time());
403  /*
404  std::cout << "ieta/iphi = " << ieta << " " << iphi
405  << ", Energy em/had/long/short = "
406  << energyemHF << " " << energyhadHF << " "
407  << theLongHitEnergy << " " << theShortHitEnergy << " "
408  << ". Time = " << theShortHit->time() << " "
409  << ". The status value is " << theStatusValue
410  << ". The short and long flags : "
411  << flagShortDPG << " " << flagLongDPG << " "
412  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
413  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
414  << ". Short fibres were cleaned." << std::endl;
415  */
416  shortFibre -= theShortHitEnergy;
417  theShortHitEnergy = 0.;
418  }
419  }
420 
421  if ( theLongHitEnergy > longFibre_Cut &&
422  ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
423  flagLongDPG ) ) {
424  // Check if the long-fibre hit was not cleaned already (because hot)
425  // In this case don't apply the cleaning
426  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId);
427  unsigned theStatusValue = theStatus->getValue();
428 
429  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
430  // The channel is killed
432  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
433 
434  //rescaleFactor = 0. ;
435  pfrhHFEMCleaned = createHcalRecHit( detid,
436  theLongHitEnergy,
438  hcalEndcapGeometry,
439  ct.id() );
440  pfrhHFEMCleaned->setTime(theLongHit->time());
441  /*
442  std::cout << "ieta/iphi = " << ieta << " " << iphi
443  << ", Energy em/had/long/short = "
444  << energyemHF << " " << energyhadHF << " "
445  << theLongHitEnergy << " " << theShortHitEnergy << " "
446  << ". The status value is " << theStatusValue
447  << ". Time = " << theLongHit->time() << " "
448  << ". The short and long flags : "
449  << flagShortDPG << " " << flagLongDPG << " "
450  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
451  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
452  << ". Long fibres were cleaned." << std::endl;
453  */
454  longFibre -= theLongHitEnergy;
455  theLongHitEnergy = 0.;
456  }
457  }
458 
459  // Special treatment for tower 29
460  // A tower with energy only at ieta = +/- 29 is not physical -> Clean
461  if ( abs(ieta) == 29 ) {
462  // rescaleFactor = 0. ;
463  // Clean long fibres
464  if ( theLongHitEnergy > shortFibre_Cut/2. ) {
465  pfrhHFEMCleaned29 = createHcalRecHit( detid,
466  theLongHitEnergy,
468  hcalEndcapGeometry,
469  ct.id() );
470  pfrhHFEMCleaned29->setTime(theLongHit->time());
471  /*
472  std::cout << "ieta/iphi = " << ieta << " " << iphi
473  << ", Energy em/had/long/short = "
474  << energyemHF << " " << energyhadHF << " "
475  << theLongHitEnergy << " " << theShortHitEnergy << " "
476  << ". Long fibres were cleaned." << std::endl;
477  */
478  longFibre -= theLongHitEnergy;
479  theLongHitEnergy = 0.;
480  }
481  // Clean short fibres
482  if ( theShortHitEnergy > shortFibre_Cut/2. ) {
483  pfrhHFHADCleaned29 = createHcalRecHit( detid,
484  theShortHitEnergy,
486  hcalEndcapGeometry,
487  ct.id());
488  pfrhHFHADCleaned29->setTime(theShortHit->time());
489  /*
490  std::cout << "ieta/iphi = " << ieta << " " << iphi
491  << ", Energy em/had/long/short = "
492  << energyemHF << " " << energyhadHF << " "
493  << theLongHitEnergy << " " << theShortHitEnergy << " "
494  << ". Short fibres were cleaned." << std::endl;
495  */
496  shortFibre -= theShortHitEnergy;
497  theShortHitEnergy = 0.;
498  }
499  }
500  // Check the timing of the long and short fibre rechits
501 
502  // First, check the timing of long and short fibre in eta = 29 if tower 30.
503  else if ( abs(ieta) == 30 ) {
504  int ieta29 = ieta > 0 ? 29 : -29;
505  HcalDetId theLongDetId29 (HcalForward, ieta29, iphi, 1);
506  HcalDetId theShortDetId29 (HcalForward, ieta29, iphi, 2);
507  auto theLongHit29 = hfHandle->find(theLongDetId29);
508  auto theShortHit29 = hfHandle->find(theShortDetId29);
509  //
510  double theLongHitEnergy29 = 0.;
511  double theShortHitEnergy29 = 0.;
512  bool flagShortDPG29 = false;
513  bool flagLongDPG29 = false;
514  bool flagShortTimeDPG29 = false;
515  bool flagLongTimeDPG29 = false;
516  bool flagShortPulseDPG29 = false;
517  bool flagLongPulseDPG29 = false;
518  //
519  if ( theLongHit29 != hfHandle->end() ) {
520  int theLongFlag29 = theLongHit29->flags();
521  theLongHitEnergy29 = theLongHit29->energy() ;
522  flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
523  flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
524  flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
525 
526  //flagLongDPG29 = applyLongShortDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
527  //flagLongTimeDPG29 = applyTimeDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
528  //flagLongPulseDPG29 = applyPulseDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
529  }
530  //
531  if ( theShortHit29 != hfHandle->end() ) {
532  int theShortFlag29 = theShortHit29->flags();
533  theShortHitEnergy29 = theShortHit29->energy();
534  flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
535  flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
536  flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
537 
538  //flagShortDPG29 = applyLongShortDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
539  //flagShortTimeDPG29 = applyTimeDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
540  //flagShortPulseDPG29 = applyPulseDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
541  }
542 
543  if ( theLongHitEnergy29 > longShortFibre_Cut &&
544  ( theLongHit29->time() < minLongTiming_Cut ||
545  theLongHit29->time() > maxLongTiming_Cut ||
546  flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
547  //rescaleFactor = 0. ;
548  pfrhHFEMCleaned29 = createHcalRecHit( detid,
549  theLongHitEnergy29,
551  hcalEndcapGeometry,
552  ct.id() );
553  pfrhHFEMCleaned29->setTime(theLongHit29->time());
554  /*
555  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
556  << ", Energy em/had/long/short = "
557  << energyemHF << " " << energyhadHF << " "
558  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
559  << ". Time = " << theLongHit29->time() << " "
560  << ". The short and long flags : "
561  << flagShortDPG29 << " " << flagLongDPG29 << " "
562  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
563  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
564  << ". Long fibres were cleaned." << std::endl;
565  */
566  longFibre -= theLongHitEnergy29;
567  theLongHitEnergy29 = 0;
568  }
569 
570  if ( theShortHitEnergy29 > longShortFibre_Cut &&
571  ( theShortHit29->time() < minShortTiming_Cut ||
572  theShortHit29->time() > maxShortTiming_Cut ||
573  flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
574  //rescaleFactor = 0. ;
575  pfrhHFHADCleaned29 = createHcalRecHit( detid,
576  theShortHitEnergy29,
578  hcalEndcapGeometry,
579  ct.id() );
580  pfrhHFHADCleaned29->setTime(theShortHit29->time());
581  /*
582  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
583  << ", Energy em/had/long/short = "
584  << energyemHF << " " << energyhadHF << " "
585  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
586  << ". Time = " << theShortHit29->time() << " "
587  << ". The short and long flags : "
588  << flagShortDPG29 << " " << flagLongDPG29 << " "
589  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
590  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
591  << ". Short fibres were cleaned." << std::endl;
592  */
593  shortFibre -= theShortHitEnergy29;
594  theShortHitEnergy29 = 0;
595  }
596 
597  // Some energy must be in the long fibres is there is some energy in the short fibres !
598  if ( theShortHitEnergy29 > shortFibre_Cut &&
599  ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
600  flagShortDPG29 ) ) {
601  // Check if the long-fibre hit was not cleaned already (because hot)
602  // In this case don't apply the cleaning
603  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId29);
604  unsigned theStatusValue = theStatus->getValue();
605 
606  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
607  // The channel is killed
609  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
610  //rescaleFactor = 0. ;
611  pfrhHFHADCleaned29 = createHcalRecHit( detid,
612  theShortHitEnergy29,
614  hcalEndcapGeometry,
615  ct.id() );
616  pfrhHFHADCleaned29->setTime(theShortHit29->time());
617  /*
618  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
619  << ", Energy em/had/long/short = "
620  << energyemHF << " " << energyhadHF << " "
621  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
622  << ". Time = " << theShortHit29->time() << " "
623  << ". The status value is " << theStatusValue
624  << ". The short and long flags : "
625  << flagShortDPG29 << " " << flagLongDPG29 << " "
626  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
627  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
628  << ". Short fibres were cleaned." << std::endl;
629  */
630  shortFibre -= theShortHitEnergy29;
631  theShortHitEnergy29 = 0.;
632  }
633  }
634 
635  // Some energy must be in the short fibres is there is some energy in the long fibres !
636  if ( theLongHitEnergy29 > longFibre_Cut &&
637  ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
638  flagLongDPG29 ) ) {
639  // Check if the long-fibre hit was not cleaned already (because hot)
640  // In this case don't apply the cleaning
641  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
642  unsigned theStatusValue = theStatus->getValue();
643  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
644  // The channel is killed
646  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
647 
648  //rescaleFactor = 0. ;
649  pfrhHFEMCleaned29 = createHcalRecHit( detid,
650  theLongHitEnergy29,
652  hcalEndcapGeometry,
653  ct.id() );
654  pfrhHFEMCleaned29->setTime(theLongHit29->time());
655  /*
656  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
657  << ", Energy em/had/long/short = "
658  << energyemHF << " " << energyhadHF << " "
659  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
660  << ". The status value is " << theStatusValue
661  << ". Time = " << theLongHit29->time() << " "
662  << ". The short and long flags : "
663  << flagShortDPG29 << " " << flagLongDPG29 << " "
664  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
665  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
666  << ". Long fibres were cleaned." << std::endl;
667  */
668  longFibre -= theLongHitEnergy29;
669  theLongHitEnergy29 = 0.;
670  }
671  }
672 
673 
674  // Check that the energy in tower 29 is smaller than in tower 30
675  // First in long fibres
676  if ( theLongHitEnergy29 > std::max(theLongHitEnergy,shortFibre_Cut/2) ) {
677  pfrhHFEMCleaned29 = createHcalRecHit( detid,
678  theLongHitEnergy29,
680  hcalEndcapGeometry,
681  ct.id() );
682  pfrhHFEMCleaned29->setTime(theLongHit29->time());
683  /*
684  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
685  << ", Energy L29/S29/L30/S30 = "
686  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
687  << theLongHitEnergy << " " << theShortHitEnergy << " "
688  << ". Long fibres were cleaned." << std::endl;
689  */
690  longFibre -= theLongHitEnergy29;
691  theLongHitEnergy29 = 0.;
692  }
693  // Second in short fibres
694  if ( theShortHitEnergy29 > std::max(theShortHitEnergy,shortFibre_Cut/2.) ) {
695  pfrhHFHADCleaned29 = createHcalRecHit( detid,
696  theShortHitEnergy29,
698  hcalEndcapGeometry,
699  ct.id() );
700  pfrhHFHADCleaned29->setTime(theShortHit29->time());
701  /*
702  std::cout << "ieta/iphi = " << ieta << " " << iphi
703  << ", Energy L29/S29/L30/S30 = "
704  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
705  << theLongHitEnergy << " " << theShortHitEnergy << " "
706  << ". Short fibres were cleaned." << std::endl;
707  */
708  shortFibre -= theShortHitEnergy29;
709  theShortHitEnergy29 = 0.;
710  }
711  }
712 
713 
714  // Determine EM and HAD after cleaning of short and long fibres
715  energyhadHF = 2.*shortFibre;
716  energyemHF = longFibre - shortFibre;
717 
718  // The EM energy might be negative, as it amounts to Long - Short
719  // In that case, put the EM "energy" in the HAD energy
720  // Just to avoid systematic positive bias due to "Short" high fluctuations
721  if ( energyemHF < thresh_HF_ ) {
722  energyhadHF += energyemHF;
723  energyemHF = 0.;
724  }
725 
726  // Apply HCAL calibration factors flor towers close to 29, if requested
727  if ( HF_Calib_ && abs(detid.ieta()) <= 32 ) {
728  energyhadHF *= HF_Calib_29;
729  energyemHF *= HF_Calib_29;
730  }
731 
732  // Create an EM and a HAD rechit if above threshold.
733  if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) {
734  pfrhHFEM = createHcalRecHit( detid,
735  energyemHF,
737  hcalEndcapGeometry,
738  ct.id() );
739  pfrhHFHAD = createHcalRecHit( detid,
740  energyhadHF,
742  hcalEndcapGeometry,
743  ct.id() );
744 
745  }
746 
747  }
748  break;
749  default:
750  LogError("PFCTRecHitProducerHCAL")
751  <<"CaloTower constituent: unknown layer : "
752  <<detid.subdet()<<endl;
753  }
754 
755 
756  if(pfrh) {
757  rechits->push_back( *pfrh );
758  delete pfrh;
759  }
760  if(pfrhCleaned) {
761  rechitsCleaned->push_back( *pfrhCleaned );
762  delete pfrhCleaned;
763  }
764  if(pfrhHFEM) {
765  HFEMRecHits->push_back( *pfrhHFEM );
766  delete pfrhHFEM;
767  }
768  if(pfrhHFHAD) {
769  HFHADRecHits->push_back( *pfrhHFHAD );
770  delete pfrhHFHAD;
771  }
772  if(pfrhHFEMCleaned) {
773  rechitsCleaned->push_back( *pfrhHFEMCleaned );
774  delete pfrhHFEMCleaned;
775  }
776  if(pfrhHFHADCleaned) {
777  rechitsCleaned->push_back( *pfrhHFHADCleaned );
778  delete pfrhHFHADCleaned;
779  }
780  if(pfrhHFEMCleaned29) {
781  rechitsCleaned->push_back( *pfrhHFEMCleaned29 );
782  delete pfrhHFEMCleaned29;
783  }
784  if(pfrhHFHADCleaned29) {
785  rechitsCleaned->push_back( *pfrhHFHADCleaned29 );
786  delete pfrhHFHADCleaned29;
787  }
788  }
789  }
790 
791  //ok now do navigation
792  //create a refprod here
793 
796 
797 
798  for( unsigned int i=0;i<rechits->size();++i) {
799  navigator_->associateNeighbours(rechits->at(i),rechits,refProd);
800  }
801 
802  if (navigation_HF_) {
803 
806 
807 
808  for( unsigned int i=0;i<HFEMRecHits->size();++i) {
809  navigator_->associateNeighbours(HFEMRecHits->at(i),HFEMRecHits,refProdEM);
810  }
811 
813  iEvent.getRefBeforePut<reco::PFRecHitCollection>("HFHAD");
814 
815 
816  for( unsigned int i=0;i<HFHADRecHits->size();++i) {
817  navigator_->associateNeighbours(HFHADRecHits->at(i),HFHADRecHits,refProdHAD);
818  }
819  }
820 
821  iEvent.put(std::move(rechits),"");
822  iEvent.put(std::move(rechitsCleaned),"Cleaned");
823  iEvent.put(std::move(HFEMRecHits),"HFEM");
824  iEvent.put(std::move(HFHADRecHits),"HFHAD");
825 
826 }
827 
829 
830 // ------------ method called once each job just before starting event loop ------------
831 void
833  const EventSetup& es) {
834 
835  // Get cleaned channels in the HCAL and HF
836  // HCAL channel status map ****************************************
837  edm::ESHandle<HcalChannelQuality> hcalChStatus;
838  es.get<HcalChannelQualityRcd>().get( "withTopo", hcalChStatus );
839  theHcalChStatus = hcalChStatus.product();
840 
841  // Retrieve the good/bad ECAL channels from the DB
843  es.get<EcalChannelStatusRcd>().get(ecalChStatus);
844  theEcalChStatus = ecalChStatus.product();
845 
847  es.get<CaloGeometryRecord>().get(cttopo);
848  theTowerConstituentsMap = cttopo.product();
849 }
850 
851 
854  double energy,
855  PFLayer::Layer layer,
857  const CaloTowerDetId& newDetId ) {
858 
859  std::shared_ptr<const CaloCellGeometry> thisCell = geom->getGeometry(detid);
860  if(!thisCell) {
861  edm::LogError("PFRecHitProducerHCAL")
862  <<"warning detid "<<detid.rawId()<<" not found in layer "
863  <<layer<<endl;
864  return nullptr;
865  }
866 
867  switch ( layer ) {
868  case PFLayer::HF_EM:
869  case PFLayer::HF_HAD:
870  {
871  auto zp = dynamic_cast<IdealZPrism const*>(thisCell.get());
872  assert(zp);
873  thisCell = zp->forPF();
874  }
875  break;
876  default:
877  break;
878  }
879 
880  auto rh =
881  new reco::PFRecHit(thisCell, newDetId.rawId(), layer, energy);
882 
883  return rh;
884 }
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:49
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
void beginLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &es) override
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:146
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
PFCTRecHitProducer(const edm::ParameterSet &)
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
std::vector< CaloTower >::const_iterator const_iterator
~PFCTRecHitProducer() override
std::vector< PFRecHit > PFRecHitCollection
collection of PFRecHit objects
Definition: PFRecHitFwd.h:9
reco::PFRecHit * createHcalRecHit(const DetId &detid, double energy, PFLayer::Layer layer, const CaloSubdetectorGeometry *geom, const CaloTowerDetId &newDetId)
int iEvent
Definition: GenABIO.cc:224
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:32
int ieta() const
get the cell ieta
Definition: HcalDetId.h:159
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void setTime(double time)
Definition: PFRecHit.h:79
RefProd< PROD > getRefBeforePut()
Definition: Event.h:150
const_iterator end() const
Layer
layer definition
Definition: PFLayer.h:31
int iphi() const
get the cell iphi
Definition: HcalDetId.h:161
Definition: DetId.h:18
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
iterator find(key_type k)
HLT enums.
T get() const
Definition: EventSetup.h:71
uint32_t getValue() const
T const * product() const
Definition: ESHandle.h:86
void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override
def move(src, dest)
Definition: eostools.py:511
T get(const Candidate &c)
Definition: component.h:55