test
CMS 3D CMS Logo

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