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_ = 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(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
166  const CaloTower& ct = (*ict);
167 
168 
169  // get the hadronic energy.
170 
171  // Mike: Just ask for the Hadronic part only now!
172  // Patrick : ARGH ! While this is ok for the HCAL, this is
173  // just wrong for the HF (in which em/had are artificially
174  // separated.
175  double energy = ct.hadEnergy();
176  //Auguste: Photons in HF have no hadEnergy in fastsim: -> all RecHit collections are empty with photons.
177  double energyEM = ct.emEnergy(); // For HF !
178  //so test the total energy to deal with the photons in HF:
179  if( (energy+energyEM) < 1e-9 ) continue;
180 
181  //get the constituents of the tower
182  const std::vector<DetId>& hits = ct.constituents();
183  const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id());
184  HcalDetId detid;
185  bool foundHCALConstituent = false;
186  //Loop on the calotower constituents and search for HCAL
187  double dead = 0.;
188  double alive = 0.;
189  for(unsigned int i=0;i< hits.size();++i) {
190  if(hits[i].det()==DetId::Hcal) {
191  foundHCALConstituent = true;
192  detid = hits[i];
193  if (theHcalTopology->withSpecialRBXHBHE() &&
194  detid.subdet() == HcalEndcap) {
195  detid = theHcalTopology->idFront(detid);
196  }
197  // An HCAL tower was found: Look for dead ECAL channels in the same CaloTower.
198  if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
199  for(unsigned int j=0;j<allConstituents.size();++j) {
200  if ( allConstituents[j].det()==DetId::Ecal ) {
201  alive += 1.;
202  EcalChannelStatus::const_iterator chIt = theEcalChStatus->find(allConstituents[j]);
203  unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
204  if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
205  }
206  }
207  }
208  // Protection: tower 29 in HF is merged with tower 30.
209  // Just take the position of tower 30 in that case.
210  if ( detid.subdet() == HcalForward && abs(detid.ieta()) == 29 ) continue;
211  break;
212  }
213  }
214 
215  // In case of dead ECAL channel, rescale the HCAL energy...
216  double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
217 
218  reco::PFRecHit* pfrh = 0;
219  reco::PFRecHit* pfrhCleaned = 0;
220  //---ab: need 2 rechits for the HF:
221  reco::PFRecHit* pfrhHFEM = 0;
222  reco::PFRecHit* pfrhHFHAD = 0;
223  reco::PFRecHit* pfrhHFEMCleaned = 0;
224  reco::PFRecHit* pfrhHFHADCleaned = 0;
225  reco::PFRecHit* pfrhHFEMCleaned29 = 0;
226  reco::PFRecHit* pfrhHFHADCleaned29 = 0;
227 
228  if(foundHCALConstituent)
229  {
230  // std::cout << ", new Energy = " << energy << std::endl;
231  switch( detid.subdet() ) {
232  case HcalBarrel:
233  {
234  if(energy < thresh_Barrel_ ) continue;
235  if ( rescaleFactor > 1. ) {
236  pfrhCleaned = createHcalRecHit( detid,
237  energy,
239  hcalBarrelGeometry,
240  ct.id() );
241  pfrhCleaned->setTime(rescaleFactor);
242  energy *= rescaleFactor;
243  }
244  pfrh = createHcalRecHit( detid,
245  energy,
247  hcalBarrelGeometry,
248  ct.id() );
249  pfrh->setTime(rescaleFactor);
250  }
251  break;
252  case HcalEndcap:
253  {
254  if(energy < thresh_Endcap_ ) continue;
255  // Apply tower 29 calibration
256  if ( HCAL_Calib_ && abs(detid.ieta()) == 29 ) energy *= HCAL_Calib_29;
257  if ( rescaleFactor > 1. ) {
258  pfrhCleaned = createHcalRecHit( detid,
259  energy,
261  hcalEndcapGeometry,
262  ct.id() );
263  pfrhCleaned->setTime(rescaleFactor);
264  energy *= rescaleFactor;
265  }
266  pfrh = createHcalRecHit( detid,
267  energy,
269  hcalEndcapGeometry,
270  ct.id() );
271  pfrh->setTime(rescaleFactor);
272  }
273  break;
274  case HcalOuter:
275  {
276  }
277  break;
278  case HcalForward:
279  {
280  //---ab: 2 rechits for HF:
281  //double energyemHF = weight_HFem_*ct.emEnergy();
282  //double energyhadHF = weight_HFhad_*ct.hadEnergy();
283  double energyemHF = weight_HFem_ * energyEM;
284  double energyhadHF = weight_HFhad_ * energy;
285  // Some energy in the tower !
286  if((energyemHF+energyhadHF) < thresh_HF_ ) continue;
287 
288  // Some cleaning in the HF
289  double longFibre = energyemHF + energyhadHF/2.;
290  double shortFibre = energyhadHF/2.;
291  int ieta = detid.ieta();
292  int iphi = detid.iphi();
293  HcalDetId theLongDetId (HcalForward, ieta, iphi, 1);
294  HcalDetId theShortDetId (HcalForward, ieta, iphi, 2);
296  iHF theLongHit = hfHandle->find(theLongDetId);
297  iHF theShortHit = hfHandle->find(theShortDetId);
298  //
299  double theLongHitEnergy = 0.;
300  double theShortHitEnergy = 0.;
301  bool flagShortDPG = false;
302  bool flagLongDPG = false;
303  bool flagShortTimeDPG = false;
304  bool flagLongTimeDPG = false;
305  bool flagShortPulseDPG = false;
306  bool flagLongPulseDPG = false;
307  //
308  if ( theLongHit != hfHandle->end() ) {
309  int theLongFlag = theLongHit->flags();
310  theLongHitEnergy = theLongHit->energy();
311  flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
312  flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
313  flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
314 
315  //flagLongDPG = applyLongShortDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
316  //flagLongTimeDPG = applyTimeDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
317  //flagLongPulseDPG = applyPulseDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
318  }
319  //
320  if ( theShortHit != hfHandle->end() ) {
321  int theShortFlag = theShortHit->flags();
322  theShortHitEnergy = theShortHit->energy();
323  flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
324  flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
325  flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
326  //flagShortDPG = applyLongShortDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
327  //flagShortTimeDPG = applyTimeDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
328  //flagShortPulseDPG = applyPulseDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
329  }
330 
331  // Then check the timing in short and long fibres in all other towers.
332  if ( theShortHitEnergy > longShortFibre_Cut &&
333  ( theShortHit->time() < minShortTiming_Cut ||
334  theShortHit->time() > maxShortTiming_Cut ||
335  flagShortTimeDPG || flagShortPulseDPG ) ) {
336  // rescaleFactor = 0. ;
337  pfrhHFHADCleaned = createHcalRecHit( detid,
338  theShortHitEnergy,
340  hcalEndcapGeometry,
341  ct.id() );
342  pfrhHFHADCleaned->setTime(theShortHit->time());
343  /*
344  std::cout << "ieta/iphi = " << ieta << " " << iphi
345  << ", Energy em/had/long/short = "
346  << energyemHF << " " << energyhadHF << " "
347  << theLongHitEnergy << " " << theShortHitEnergy << " "
348  << ". Time = " << theShortHit->time() << " "
349  << ". The short and long flags : "
350  << flagShortDPG << " " << flagLongDPG << " "
351  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
352  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
353  << ". Short fibres were cleaned." << std::endl;
354  */
355  shortFibre -= theShortHitEnergy;
356  theShortHitEnergy = 0.;
357  }
358 
359 
360  if ( theLongHitEnergy > longShortFibre_Cut &&
361  ( theLongHit->time() < minLongTiming_Cut ||
362  theLongHit->time() > maxLongTiming_Cut ||
363  flagLongTimeDPG || flagLongPulseDPG ) ) {
364  //rescaleFactor = 0. ;
365  pfrhHFEMCleaned = createHcalRecHit( detid,
366  theLongHitEnergy,
368  hcalEndcapGeometry,
369  ct.id());
370  pfrhHFEMCleaned->setTime(theLongHit->time());
371  /*
372  std::cout << "ieta/iphi = " << ieta << " " << iphi
373  << ", Energy em/had/long/short = "
374  << energyemHF << " " << energyhadHF << " "
375  << theLongHitEnergy << " " << theShortHitEnergy << " "
376  << ". Time = " << theLongHit->time() << " "
377  << ". The short and long flags : "
378  << flagShortDPG << " " << flagLongDPG << " "
379  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
380  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
381  << ". Long fibres were cleaned." << std::endl;
382  */
383  longFibre -= theLongHitEnergy;
384  theLongHitEnergy = 0.;
385  }
386 
387  // Some energy must be in the long fibres is there is some energy in the short fibres !
388  if ( theShortHitEnergy > shortFibre_Cut &&
389  ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
390  flagShortDPG ) ) {
391  // Check if the long-fibre hit was not cleaned already (because hot)
392  // In this case don't apply the cleaning
393  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId);
394  unsigned theStatusValue = theStatus->getValue();
395  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
396  // The channel is killed
398  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
399  // rescaleFactor = 0. ;
400  pfrhHFHADCleaned = createHcalRecHit( detid,
401  theShortHitEnergy,
403  hcalEndcapGeometry,
404  ct.id() );
405  pfrhHFHADCleaned->setTime(theShortHit->time());
406  /*
407  std::cout << "ieta/iphi = " << ieta << " " << iphi
408  << ", Energy em/had/long/short = "
409  << energyemHF << " " << energyhadHF << " "
410  << theLongHitEnergy << " " << theShortHitEnergy << " "
411  << ". Time = " << theShortHit->time() << " "
412  << ". The status value is " << theStatusValue
413  << ". The short and long flags : "
414  << flagShortDPG << " " << flagLongDPG << " "
415  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
416  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
417  << ". Short fibres were cleaned." << std::endl;
418  */
419  shortFibre -= theShortHitEnergy;
420  theShortHitEnergy = 0.;
421  }
422  }
423 
424  if ( theLongHitEnergy > longFibre_Cut &&
425  ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
426  flagLongDPG ) ) {
427  // Check if the long-fibre hit was not cleaned already (because hot)
428  // In this case don't apply the cleaning
429  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId);
430  unsigned theStatusValue = theStatus->getValue();
431 
432  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
433  // The channel is killed
435  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
436 
437  //rescaleFactor = 0. ;
438  pfrhHFEMCleaned = createHcalRecHit( detid,
439  theLongHitEnergy,
441  hcalEndcapGeometry,
442  ct.id() );
443  pfrhHFEMCleaned->setTime(theLongHit->time());
444  /*
445  std::cout << "ieta/iphi = " << ieta << " " << iphi
446  << ", Energy em/had/long/short = "
447  << energyemHF << " " << energyhadHF << " "
448  << theLongHitEnergy << " " << theShortHitEnergy << " "
449  << ". The status value is " << theStatusValue
450  << ". Time = " << theLongHit->time() << " "
451  << ". The short and long flags : "
452  << flagShortDPG << " " << flagLongDPG << " "
453  << flagShortTimeDPG << " " << flagLongTimeDPG << " "
454  << flagShortPulseDPG << " " << flagLongPulseDPG << " "
455  << ". Long fibres were cleaned." << std::endl;
456  */
457  longFibre -= theLongHitEnergy;
458  theLongHitEnergy = 0.;
459  }
460  }
461 
462  // Special treatment for tower 29
463  // A tower with energy only at ieta = +/- 29 is not physical -> Clean
464  if ( abs(ieta) == 29 ) {
465  // rescaleFactor = 0. ;
466  // Clean long fibres
467  if ( theLongHitEnergy > shortFibre_Cut/2. ) {
468  pfrhHFEMCleaned29 = createHcalRecHit( detid,
469  theLongHitEnergy,
471  hcalEndcapGeometry,
472  ct.id() );
473  pfrhHFEMCleaned29->setTime(theLongHit->time());
474  /*
475  std::cout << "ieta/iphi = " << ieta << " " << iphi
476  << ", Energy em/had/long/short = "
477  << energyemHF << " " << energyhadHF << " "
478  << theLongHitEnergy << " " << theShortHitEnergy << " "
479  << ". Long fibres were cleaned." << std::endl;
480  */
481  longFibre -= theLongHitEnergy;
482  theLongHitEnergy = 0.;
483  }
484  // Clean short fibres
485  if ( theShortHitEnergy > shortFibre_Cut/2. ) {
486  pfrhHFHADCleaned29 = createHcalRecHit( detid,
487  theShortHitEnergy,
489  hcalEndcapGeometry,
490  ct.id());
491  pfrhHFHADCleaned29->setTime(theShortHit->time());
492  /*
493  std::cout << "ieta/iphi = " << ieta << " " << iphi
494  << ", Energy em/had/long/short = "
495  << energyemHF << " " << energyhadHF << " "
496  << theLongHitEnergy << " " << theShortHitEnergy << " "
497  << ". Short fibres were cleaned." << std::endl;
498  */
499  shortFibre -= theShortHitEnergy;
500  theShortHitEnergy = 0.;
501  }
502  }
503  // Check the timing of the long and short fibre rechits
504 
505  // First, check the timing of long and short fibre in eta = 29 if tower 30.
506  else if ( abs(ieta) == 30 ) {
507  int ieta29 = ieta > 0 ? 29 : -29;
508  HcalDetId theLongDetId29 (HcalForward, ieta29, iphi, 1);
509  HcalDetId theShortDetId29 (HcalForward, ieta29, iphi, 2);
510  iHF theLongHit29 = hfHandle->find(theLongDetId29);
511  iHF theShortHit29 = hfHandle->find(theShortDetId29);
512  //
513  double theLongHitEnergy29 = 0.;
514  double theShortHitEnergy29 = 0.;
515  bool flagShortDPG29 = false;
516  bool flagLongDPG29 = false;
517  bool flagShortTimeDPG29 = false;
518  bool flagLongTimeDPG29 = false;
519  bool flagShortPulseDPG29 = false;
520  bool flagLongPulseDPG29 = false;
521  //
522  if ( theLongHit29 != hfHandle->end() ) {
523  int theLongFlag29 = theLongHit29->flags();
524  theLongHitEnergy29 = theLongHit29->energy() ;
525  flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
526  flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
527  flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
528 
529  //flagLongDPG29 = applyLongShortDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
530  //flagLongTimeDPG29 = applyTimeDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
531  //flagLongPulseDPG29 = applyPulseDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
532  }
533  //
534  if ( theShortHit29 != hfHandle->end() ) {
535  int theShortFlag29 = theShortHit29->flags();
536  theShortHitEnergy29 = theShortHit29->energy();
537  flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
538  flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
539  flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
540 
541  //flagShortDPG29 = applyLongShortDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
542  //flagShortTimeDPG29 = applyTimeDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
543  //flagShortPulseDPG29 = applyPulseDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
544  }
545 
546  if ( theLongHitEnergy29 > longShortFibre_Cut &&
547  ( theLongHit29->time() < minLongTiming_Cut ||
548  theLongHit29->time() > maxLongTiming_Cut ||
549  flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
550  //rescaleFactor = 0. ;
551  pfrhHFEMCleaned29 = createHcalRecHit( detid,
552  theLongHitEnergy29,
554  hcalEndcapGeometry,
555  ct.id() );
556  pfrhHFEMCleaned29->setTime(theLongHit29->time());
557  /*
558  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
559  << ", Energy em/had/long/short = "
560  << energyemHF << " " << energyhadHF << " "
561  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
562  << ". Time = " << theLongHit29->time() << " "
563  << ". The short and long flags : "
564  << flagShortDPG29 << " " << flagLongDPG29 << " "
565  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
566  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
567  << ". Long fibres were cleaned." << std::endl;
568  */
569  longFibre -= theLongHitEnergy29;
570  theLongHitEnergy29 = 0;
571  }
572 
573  if ( theShortHitEnergy29 > longShortFibre_Cut &&
574  ( theShortHit29->time() < minShortTiming_Cut ||
575  theShortHit29->time() > maxShortTiming_Cut ||
576  flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
577  //rescaleFactor = 0. ;
578  pfrhHFHADCleaned29 = createHcalRecHit( detid,
579  theShortHitEnergy29,
581  hcalEndcapGeometry,
582  ct.id() );
583  pfrhHFHADCleaned29->setTime(theShortHit29->time());
584  /*
585  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
586  << ", Energy em/had/long/short = "
587  << energyemHF << " " << energyhadHF << " "
588  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
589  << ". Time = " << theShortHit29->time() << " "
590  << ". The short and long flags : "
591  << flagShortDPG29 << " " << flagLongDPG29 << " "
592  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
593  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
594  << ". Short fibres were cleaned." << std::endl;
595  */
596  shortFibre -= theShortHitEnergy29;
597  theShortHitEnergy29 = 0;
598  }
599 
600  // Some energy must be in the long fibres is there is some energy in the short fibres !
601  if ( theShortHitEnergy29 > shortFibre_Cut &&
602  ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
603  flagShortDPG29 ) ) {
604  // Check if the long-fibre hit was not cleaned already (because hot)
605  // In this case don't apply the cleaning
606  const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId29);
607  unsigned theStatusValue = theStatus->getValue();
608 
609  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
610  // The channel is killed
612  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
613  //rescaleFactor = 0. ;
614  pfrhHFHADCleaned29 = createHcalRecHit( detid,
615  theShortHitEnergy29,
617  hcalEndcapGeometry,
618  ct.id() );
619  pfrhHFHADCleaned29->setTime(theShortHit29->time());
620  /*
621  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
622  << ", Energy em/had/long/short = "
623  << energyemHF << " " << energyhadHF << " "
624  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
625  << ". Time = " << theShortHit29->time() << " "
626  << ". The status value is " << theStatusValue
627  << ". The short and long flags : "
628  << flagShortDPG29 << " " << flagLongDPG29 << " "
629  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
630  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
631  << ". Short fibres were cleaned." << std::endl;
632  */
633  shortFibre -= theShortHitEnergy29;
634  theShortHitEnergy29 = 0.;
635  }
636  }
637 
638  // Some energy must be in the short fibres is there is some energy in the long fibres !
639  if ( theLongHitEnergy29 > longFibre_Cut &&
640  ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
641  flagLongDPG29 ) ) {
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(theShortDetId29);
645  unsigned theStatusValue = theStatus->getValue();
646  int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
647  // The channel is killed
649  if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
650 
651  //rescaleFactor = 0. ;
652  pfrhHFEMCleaned29 = createHcalRecHit( detid,
653  theLongHitEnergy29,
655  hcalEndcapGeometry,
656  ct.id() );
657  pfrhHFEMCleaned29->setTime(theLongHit29->time());
658  /*
659  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
660  << ", Energy em/had/long/short = "
661  << energyemHF << " " << energyhadHF << " "
662  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
663  << ". The status value is " << theStatusValue
664  << ". Time = " << theLongHit29->time() << " "
665  << ". The short and long flags : "
666  << flagShortDPG29 << " " << flagLongDPG29 << " "
667  << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "
668  << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "
669  << ". Long fibres were cleaned." << std::endl;
670  */
671  longFibre -= theLongHitEnergy29;
672  theLongHitEnergy29 = 0.;
673  }
674  }
675 
676 
677  // Check that the energy in tower 29 is smaller than in tower 30
678  // First in long fibres
679  if ( theLongHitEnergy29 > std::max(theLongHitEnergy,shortFibre_Cut/2) ) {
680  pfrhHFEMCleaned29 = createHcalRecHit( detid,
681  theLongHitEnergy29,
683  hcalEndcapGeometry,
684  ct.id() );
685  pfrhHFEMCleaned29->setTime(theLongHit29->time());
686  /*
687  std::cout << "ieta/iphi = " << ieta29 << " " << iphi
688  << ", Energy L29/S29/L30/S30 = "
689  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
690  << theLongHitEnergy << " " << theShortHitEnergy << " "
691  << ". Long fibres were cleaned." << std::endl;
692  */
693  longFibre -= theLongHitEnergy29;
694  theLongHitEnergy29 = 0.;
695  }
696  // Second in short fibres
697  if ( theShortHitEnergy29 > std::max(theShortHitEnergy,shortFibre_Cut/2.) ) {
698  pfrhHFHADCleaned29 = createHcalRecHit( detid,
699  theShortHitEnergy29,
701  hcalEndcapGeometry,
702  ct.id() );
703  pfrhHFHADCleaned29->setTime(theShortHit29->time());
704  /*
705  std::cout << "ieta/iphi = " << ieta << " " << iphi
706  << ", Energy L29/S29/L30/S30 = "
707  << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
708  << theLongHitEnergy << " " << theShortHitEnergy << " "
709  << ". Short fibres were cleaned." << std::endl;
710  */
711  shortFibre -= theShortHitEnergy29;
712  theShortHitEnergy29 = 0.;
713  }
714  }
715 
716 
717  // Determine EM and HAD after cleaning of short and long fibres
718  energyhadHF = 2.*shortFibre;
719  energyemHF = longFibre - shortFibre;
720 
721  // The EM energy might be negative, as it amounts to Long - Short
722  // In that case, put the EM "energy" in the HAD energy
723  // Just to avoid systematic positive bias due to "Short" high fluctuations
724  if ( energyemHF < thresh_HF_ ) {
725  energyhadHF += energyemHF;
726  energyemHF = 0.;
727  }
728 
729  // Apply HCAL calibration factors flor towers close to 29, if requested
730  if ( HF_Calib_ && abs(detid.ieta()) <= 32 ) {
731  energyhadHF *= HF_Calib_29;
732  energyemHF *= HF_Calib_29;
733  }
734 
735  // Create an EM and a HAD rechit if above threshold.
736  if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) {
737  pfrhHFEM = createHcalRecHit( detid,
738  energyemHF,
740  hcalEndcapGeometry,
741  ct.id() );
742  pfrhHFHAD = createHcalRecHit( detid,
743  energyhadHF,
745  hcalEndcapGeometry,
746  ct.id() );
747 
748  }
749 
750  }
751  break;
752  default:
753  LogError("PFCTRecHitProducerHCAL")
754  <<"CaloTower constituent: unknown layer : "
755  <<detid.subdet()<<endl;
756  }
757 
758 
759  if(pfrh) {
760  rechits->push_back( *pfrh );
761  delete pfrh;
762  }
763  if(pfrhCleaned) {
764  rechitsCleaned->push_back( *pfrhCleaned );
765  delete pfrhCleaned;
766  }
767  if(pfrhHFEM) {
768  HFEMRecHits->push_back( *pfrhHFEM );
769  delete pfrhHFEM;
770  }
771  if(pfrhHFHAD) {
772  HFHADRecHits->push_back( *pfrhHFHAD );
773  delete pfrhHFHAD;
774  }
775  if(pfrhHFEMCleaned) {
776  rechitsCleaned->push_back( *pfrhHFEMCleaned );
777  delete pfrhHFEMCleaned;
778  }
779  if(pfrhHFHADCleaned) {
780  rechitsCleaned->push_back( *pfrhHFHADCleaned );
781  delete pfrhHFHADCleaned;
782  }
783  if(pfrhHFEMCleaned29) {
784  rechitsCleaned->push_back( *pfrhHFEMCleaned29 );
785  delete pfrhHFEMCleaned29;
786  }
787  if(pfrhHFHADCleaned29) {
788  rechitsCleaned->push_back( *pfrhHFHADCleaned29 );
789  delete pfrhHFHADCleaned29;
790  }
791  }
792  }
793 
794  //ok now do navigation
795  //create a refprod here
796 
799 
800 
801  for( unsigned int i=0;i<rechits->size();++i) {
802  navigator_->associateNeighbours(rechits->at(i),rechits,refProd);
803  }
804 
805  if (navigation_HF_) {
806 
809 
810 
811  for( unsigned int i=0;i<HFEMRecHits->size();++i) {
812  navigator_->associateNeighbours(HFEMRecHits->at(i),HFEMRecHits,refProdEM);
813  }
814 
816  iEvent.getRefBeforePut<reco::PFRecHitCollection>("HFHAD");
817 
818 
819  for( unsigned int i=0;i<HFHADRecHits->size();++i) {
820  navigator_->associateNeighbours(HFHADRecHits->at(i),HFHADRecHits,refProdHAD);
821  }
822  }
823 
824  iEvent.put(std::move(rechits),"");
825  iEvent.put(std::move(rechitsCleaned),"Cleaned");
826  iEvent.put(std::move(HFEMRecHits),"HFEM");
827  iEvent.put(std::move(HFHADRecHits),"HFHAD");
828 
829 }
830 
832 
833 // ------------ method called once each job just before starting event loop ------------
834 void
836  const EventSetup& es) {
837 
838  // Get cleaned channels in the HCAL and HF
839  // HCAL channel status map ****************************************
840  edm::ESHandle<HcalChannelQuality> hcalChStatus;
841  es.get<HcalChannelQualityRcd>().get( "withTopo", hcalChStatus );
842  theHcalChStatus = hcalChStatus.product();
843 
844  // Retrieve the good/bad ECAL channels from the DB
846  es.get<EcalChannelStatusRcd>().get(ecalChStatus);
847  theEcalChStatus = ecalChStatus.product();
848 
850  es.get<CaloGeometryRecord>().get(cttopo);
851  theTowerConstituentsMap = cttopo.product();
852 }
853 
854 
857  double energy,
858  PFLayer::Layer layer,
860  const CaloTowerDetId& newDetId ) {
861 
862  const CaloCellGeometry *thisCell = geom->getGeometry(detid);
863  if(!thisCell) {
864  edm::LogError("PFRecHitProducerHCAL")
865  <<"warning detid "<<detid.rawId()<<" not found in layer "
866  <<layer<<endl;
867  return 0;
868  }
869 
870  switch ( layer ) {
871  case PFLayer::HF_EM:
872  case PFLayer::HF_HAD:
873  {
874  auto zp = dynamic_cast<IdealZPrism const*>(thisCell);
875  assert(zp);
876  thisCell = zp->forPF();
877  }
878  break;
879  default:
880  break;
881  }
882 
883  reco::PFRecHit *rh =
884  new reco::PFRecHit(thisCell, newDetId.rawId(), layer, energy);
885 
886  return rh;
887 }
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:45
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:122
virtual void beginLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &es) override
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:49
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
PFCTRecHitProducer(const edm::ParameterSet &)
std::vector< CaloTower >::const_iterator const_iterator
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)
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
int iEvent
Definition: GenABIO.cc:230
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
double emEnergy() const
Definition: CaloTower.h:110
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const std::vector< DetId > & constituents() const
Definition: CaloTower.h:104
void setTime(double time)
Definition: PFRecHit.h:79
RefProd< PROD > getRefBeforePut()
Definition: Event.h:134
double hadEnergy() const
Definition: CaloTower.h:111
const_iterator end() const
Layer
layer definition
Definition: PFLayer.h:31
int iphi() const
get the cell iphi
Definition: HcalDetId.cc:103
Definition: DetId.h:18
CaloTowerDetId id() const
Definition: CaloTower.h:103
const T & get() const
Definition: EventSetup.h:56
std::vector< Item >::const_iterator const_iterator
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
iterator find(key_type k)
HLT enums.
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:510
T get(const Candidate &c)
Definition: component.h:55