Reference Manual

Go to the documentation of this file.
00001 #include "RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducerHCAL.h"
00002 
00003 #include <memory>
00004 
00005 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
00006 
00007 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
00008 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
00009 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00010 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00011 
00012 #include "DataFormats/DetId/interface/DetId.h"
00013 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00014 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00015 
00016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 
00020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00021 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00022 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00023 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
00024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00025 
00026 #include "Geometry/CaloTopology/interface/HcalTopology.h"
00027 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
00028 
00029 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00030 
00031 #include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
00032 #include "RecoCaloTools/Navigation/interface/CaloTowerNavigator.h"
00033 
00034 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
00035 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
00036 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
00037 
00038 using namespace std;
00039 using namespace edm;
00040 
00041 PFRecHitProducerHCAL::PFRecHitProducerHCAL(const edm::ParameterSet& iConfig)
00042   : PFRecHitProducer( iConfig )
00043 {
00044   // access to the collections of rechits 
00045   
00046   inputTagHcalRecHitsHBHE_ =
00047     iConfig.getParameter<InputTag>("hcalRecHitsHBHE");
00048     
00049   inputTagHcalRecHitsHF_ =
00050     iConfig.getParameter<InputTag>("hcalRecHitsHF");
00051     
00052  
00053   inputTagCaloTowers_ = 
00054     iConfig.getParameter<InputTag>("caloTowers");
00055    
00056   thresh_HF_ = 
00057     iConfig.getParameter<double>("thresh_HF");
00058   navigation_HF_ = 
00059     iConfig.getParameter<bool>("navigation_HF");
00060   weight_HFem_ =
00061     iConfig.getParameter<double>("weight_HFem");
00062   weight_HFhad_ =
00063     iConfig.getParameter<double>("weight_HFhad");
00064 
00065   HCAL_Calib_ =
00066     iConfig.getParameter<bool>("HCAL_Calib");
00067   HF_Calib_ =
00068     iConfig.getParameter<bool>("HF_Calib");
00069   HCAL_Calib_29 = 
00070     iConfig.getParameter<double>("HCAL_Calib_29");
00071   HF_Calib_29 = 
00072     iConfig.getParameter<double>("HF_Calib_29");
00073 
00074   shortFibre_Cut = iConfig.getParameter<double>("ShortFibre_Cut");
00075   longFibre_Fraction = iConfig.getParameter<double>("LongFibre_Fraction");
00076 
00077   longFibre_Cut = iConfig.getParameter<double>("LongFibre_Cut");
00078   shortFibre_Fraction = iConfig.getParameter<double>("ShortFibre_Fraction");
00079 
00080   applyLongShortDPG_ = iConfig.getParameter<bool>("ApplyLongShortDPG");
00081 
00082   longShortFibre_Cut = iConfig.getParameter<double>("LongShortFibre_Cut");
00083   minShortTiming_Cut = iConfig.getParameter<double>("MinShortTiming_Cut");
00084   maxShortTiming_Cut = iConfig.getParameter<double>("MaxShortTiming_Cut");
00085   minLongTiming_Cut = iConfig.getParameter<double>("MinLongTiming_Cut");
00086   maxLongTiming_Cut = iConfig.getParameter<double>("MaxLongTiming_Cut");
00087 
00088   applyTimeDPG_ = iConfig.getParameter<bool>("ApplyTimeDPG");
00089   applyPulseDPG_ = iConfig.getParameter<bool>("ApplyPulseDPG");
00090   HcalMaxAllowedHFLongShortSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFLongShortSev");
00091   HcalMaxAllowedHFDigiTimeSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFDigiTimeSev");
00092   HcalMaxAllowedHFInTimeWindowSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFInTimeWindowSev");
00093   HcalMaxAllowedChannelStatusSev_ = iConfig.getParameter<int>("HcalMaxAllowedChannelStatusSev");
00094 
00095   ECAL_Compensate_ = iConfig.getParameter<bool>("ECAL_Compensate");
00096   ECAL_Threshold_ = iConfig.getParameter<double>("ECAL_Threshold");
00097   ECAL_Compensation_ = iConfig.getParameter<double>("ECAL_Compensation");
00098   ECAL_Dead_Code_ = iConfig.getParameter<unsigned int>("ECAL_Dead_Code");
00099 
00100   EM_Depth_ = iConfig.getParameter<double>("EM_Depth");
00101   HAD_Depth_ = iConfig.getParameter<double>("HAD_Depth");
00102 
00103   //Get integer values of individual HCAL HF flags
00104   hcalHFLongShortFlagValue_=1<<HcalCaloFlagLabels::HFLongShort;
00105   hcalHFDigiTimeFlagValue_=1<<HcalCaloFlagLabels::HFDigiTime;
00106   hcalHFInTimeWindowFlagValue_=1<<HcalCaloFlagLabels::HFInTimeWindow;
00107 
00108 
00109   //--ab
00110   produces<reco::PFRecHitCollection>("HFHAD").setBranchAlias("HFHADRecHits");
00111   produces<reco::PFRecHitCollection>("HFEM").setBranchAlias("HFEMRecHits");
00112   //--ab
00113 }
00114 
00115 
00116 
00117 PFRecHitProducerHCAL::~PFRecHitProducerHCAL() {}
00118 
00119 
00120 
00121 void PFRecHitProducerHCAL::createRecHits(vector<reco::PFRecHit>& rechits,
00122                                          vector<reco::PFRecHit>& rechitsCleaned,
00123                                          edm::Event& iEvent, 
00124                                          const edm::EventSetup& iSetup ) {
00125 
00126   
00127   // this map is necessary to find the rechit neighbours efficiently
00128   //C but I should think about using Florian's hashed index to do this.
00129   //C in which case the map might not be necessary anymore
00130   //C however the hashed index does not seem to be implemented for HCAL
00131   // 
00132   // the key of this map is detId. 
00133   // the value is the index in the rechits vector
00134   map<unsigned,  unsigned > idSortedRecHits;
00135   map<unsigned,  unsigned > idSortedRecHitsHFEM;
00136   map<unsigned,  unsigned > idSortedRecHitsHFHAD;
00137   typedef map<unsigned, unsigned >::iterator IDH;  
00138 
00139 
00140   edm::ESHandle<CaloGeometry> geoHandle;
00141   iSetup.get<CaloGeometryRecord>().get(geoHandle);
00142   
00143   // get the hcalBarrel geometry
00144   const CaloSubdetectorGeometry *hcalBarrelGeometry = 
00145     geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
00146 
00147   // get the endcap geometry
00148   const CaloSubdetectorGeometry *hcalEndcapGeometry = 
00149     geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalEndcap);
00150 
00151   // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined
00152   edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
00153   iSetup.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
00154   const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
00155 
00156   //--ab
00157   auto_ptr< vector<reco::PFRecHit> > HFHADRecHits( new vector<reco::PFRecHit> ); 
00158   auto_ptr< vector<reco::PFRecHit> > HFEMRecHits( new vector<reco::PFRecHit> ); 
00159   //--ab
00160 
00161   // 2 possibilities to make HCAL clustering :
00162   // - from the HCAL rechits
00163   // - from the CaloTowers. 
00164   // ultimately, clustering will be done taking CaloTowers as an 
00165   // input. This possibility is currently under investigation, and 
00166   // was thus made optional.
00167 
00168   // in the first step, we will fill the map of PFRecHits hcalrechits
00169   // either from CaloTowers or from HCAL rechits. 
00170 
00171   // in the second step, we will perform clustering on this map.
00172 
00173   if( !(inputTagCaloTowers_ == InputTag()) ) {
00174       
00175     edm::Handle<CaloTowerCollection> caloTowers; 
00176     CaloTowerTopology caloTowerTopology; 
00177     const CaloSubdetectorGeometry *caloTowerGeometry = 0; 
00178     // = geometry_->getSubdetectorGeometry(id)
00179 
00180     // get calotowers
00181     bool found = iEvent.getByLabel(inputTagCaloTowers_,
00182                                    caloTowers);
00183 
00184     if(!found) {
00185       ostringstream err;
00186       err<<"could not find rechits "<<inputTagCaloTowers_;
00187       LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
00188     
00189       throw cms::Exception( "MissingProduct", err.str());
00190     }
00191     else {
00192       assert( caloTowers.isValid() );
00193 
00194       // get HF rechits
00195       edm::Handle<HFRecHitCollection>  hfHandle;  
00196       found = iEvent.getByLabel(inputTagHcalRecHitsHF_,
00197                                 hfHandle);
00198       
00199       if(!found) {
00200         ostringstream err;
00201         err<<"could not find HF rechits "<<inputTagHcalRecHitsHF_;
00202         LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
00203         
00204         throw cms::Exception( "MissingProduct", err.str());
00205       }
00206       else {
00207         assert( hfHandle.isValid() );
00208       }
00209       
00210       // get HBHE rechits
00211       edm::Handle<HBHERecHitCollection>  hbheHandle;  
00212       found = iEvent.getByLabel(inputTagHcalRecHitsHBHE_,
00213                                 hbheHandle);
00214       
00215       if(!found) {
00216         ostringstream err;
00217         err<<"could not find HBHE rechits "<<inputTagHcalRecHitsHBHE_;
00218         LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
00219         
00220         throw cms::Exception( "MissingProduct", err.str());
00221       }
00222       else {
00223         assert( hbheHandle.isValid() );
00224       }
00225       
00226       // create rechits
00227       typedef CaloTowerCollection::const_iterator ICT;
00228     
00229       for(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
00230           
00231         const CaloTower& ct = (*ict);
00232           
00233         //C     
00234         if(!caloTowerGeometry) 
00235           caloTowerGeometry = geoHandle->getSubdetectorGeometry(ct.id());
00236 
00237           
00238         // get the hadronic energy.
00239         
00240         // Mike: Just ask for the Hadronic part only now!
00241         // Patrick : ARGH ! While this is ok for the HCAL, this is 
00242         // just wrong for the HF (in which em/had are artificially 
00243         // separated. 
00244         double energy = ct.hadEnergy();
00245         //Auguste: Photons in HF have no hadEnergy in fastsim: -> all RecHit collections are empty with photons.
00246         double energyEM = ct.emEnergy(); // For HF !
00247         //so test the total energy to deal with the photons in  HF:
00248         if( (energy+energyEM) < 1e-9 ) continue;
00249           
00250         assert( ct.constituentsSize() );          
00251         //Mike: The DetId will be taken by the first Hadronic constituent
00252         //         of the tower. That is only what we need
00253 
00254         
00255         //get the constituents of the tower
00256         const std::vector<DetId>& hits = ct.constituents();
00257         const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id());
00258 
00259         /*
00260         for(unsigned int i=0;i< hits.size();++i) {
00261           if(hits[i].det()==DetId::Hcal) {
00262             HcalDetId did = hits[i];
00263             if ( did.subdet()==HcalEndcap || did.subdet()==HcalForward ) { 
00264               //double en = hits[i].energy();
00265               int ieta = did.ieta();
00266               const CaloCellGeometry *thisCell = hcalEndcapGeometry->getGeometry(did);
00267               const GlobalPoint& position = thisCell->getPosition();
00268               if ( abs(ieta) > 27 && abs(ieta) < 33 && energy > 10. ) { 
00269                 std::cout << "HE/HF hit " << i << " at eta = " << ieta 
00270                           << " with CT energy = " << energy 
00271                           << " at eta, z (hit) = " << position.eta() << " " << position.z()
00272                           << " at eta, z (cte) = " << ct.emPosition().eta() << " " << ct.emPosition().z()
00273                           << " at eta, z (cth) = " << ct.hadPosition().eta() << " " << ct.hadPosition().z()
00274                           << " at eta, z (cto) = " << ct.eta() << " " << ct.vz() 
00275                           << std::endl;
00276               }
00277             }
00278           }
00279         }
00280         */
00281         
00282         //Reserve the DetId we are looking for:
00283 
00284         HcalDetId detid;
00285         // EcalDetId edetid;
00286         bool foundHCALConstituent = false;
00287 
00288  
00289         //Loop on the calotower constituents and search for HCAL
00290         double dead = 0.;
00291         double alive = 0.;
00292         for(unsigned int i=0;i< hits.size();++i) {
00293           if(hits[i].det()==DetId::Hcal) { 
00294             foundHCALConstituent = true;
00295             detid = hits[i];
00296             // An HCAL tower was found: Look for dead ECAL channels in the same CaloTower.
00297             if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
00298               for(unsigned int j=0;j<allConstituents.size();++j) { 
00299                 if ( allConstituents[j].det()==DetId::Ecal ) { 
00300                   alive += 1.;
00301                   EcalChannelStatus::const_iterator chIt = theEcalChStatus->find(allConstituents[j]);
00302                   unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
00303                   if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
00304                 }
00305               }
00306             } 
00307             // Protection: tower 29 in HF is merged with tower 30. 
00308             // Just take the position of tower 30 in that case. 
00309             if ( detid.subdet() == HcalForward && abs(detid.ieta()) == 29 ) continue; 
00310             break;
00311           }
00312         }
00313 
00314         // In case of dead ECAL channel, rescale the HCAL energy...
00315         double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
00316           
00317         reco::PFRecHit* pfrh = 0;
00318         reco::PFRecHit* pfrhCleaned = 0;
00319         //---ab: need 2 rechits for the HF:
00320         reco::PFRecHit* pfrhHFEM = 0;
00321         reco::PFRecHit* pfrhHFHAD = 0;
00322         reco::PFRecHit* pfrhHFEMCleaned = 0;
00323         reco::PFRecHit* pfrhHFHADCleaned = 0;
00324         reco::PFRecHit* pfrhHFEMCleaned29 = 0;
00325         reco::PFRecHit* pfrhHFHADCleaned29 = 0;
00326 
00327         if(foundHCALConstituent)
00328           {
00329             // std::cout << ", new Energy = " << energy << std::endl;
00330             switch( detid.subdet() ) {
00331             case HcalBarrel: 
00332               {
00333                 if(energy < thresh_Barrel_ ) continue;
00334 
00335                 /*
00336                 // Check the timing
00337                 if ( energy > 5. ) { 
00338                   for(unsigned int i=0;i< hits.size();++i) {
00339                     if( hits[i].det() != DetId::Hcal ) continue; 
00340                     HcalDetId theDetId = hits[i]; 
00341                     typedef HBHERecHitCollection::const_iterator iHBHE;
00342                     iHBHE theHit = hbheHandle->find(theDetId); 
00343                     if ( theHit != hbheHandle->end() ) 
00344                       std::cout << "HCAL hit : " 
00345                                 << theDetId.ieta() << " " << theDetId.iphi() << " " 
00346                                 << theHit->energy() << " " << theHit->time() << std::endl;
00347                   }
00348                 }
00349                 */
00350 
00351 
00352                 // if ( HCAL_Calib_ ) energy   *= std::min(max_Calib_,myPFCorr->getValues(detid)->getValue());
00353                 //if ( rescaleFactor > 1. ) 
00354                 // std::cout << "Barrel HCAL energy rescaled from = " << energy << " to " << energy*rescaleFactor << std::endl;
00355                 if ( rescaleFactor > 1. ) { 
00356                   pfrhCleaned = createHcalRecHit( detid, 
00357                                                   energy, 
00358                                                   PFLayer::HCAL_BARREL1, 
00359                                                   hcalBarrelGeometry,
00360                                                   ct.id().rawId() );
00361                   pfrhCleaned->setRescale(rescaleFactor);
00362                   energy *= rescaleFactor;
00363                 }
00364                 pfrh = createHcalRecHit( detid, 
00365                                          energy, 
00366                                  PFLayer::HCAL_BARREL1, 
00367                                          hcalBarrelGeometry,
00368                                          ct.id().rawId() );
00369                 pfrh->setRescale(rescaleFactor);
00370               }
00371               break;
00372             case HcalEndcap:
00373               {
00374                 if(energy < thresh_Endcap_ ) continue;
00375 
00376                 /*
00377                 // Check the timing
00378                 if ( energy > 5. ) { 
00379                   for(unsigned int i=0;i< hits.size();++i) {
00380                     if( hits[i].det() != DetId::Hcal ) continue; 
00381                     HcalDetId theDetId = hits[i]; 
00382                     typedef HBHERecHitCollection::const_iterator iHBHE;
00383                     iHBHE theHit = hbheHandle->find(theDetId); 
00384                     if ( theHit != hbheHandle->end() ) 
00385                       std::cout << "HCAL hit : " 
00386                                 << theDetId.ieta() << " " << theDetId.iphi() << " " 
00387                                 << theHit->energy() << " " << theHit->time() << std::endl;
00388                   }
00389                 }
00390                 */
00391 
00392                 // Apply tower 29 calibration
00393                 if ( HCAL_Calib_ && abs(detid.ieta()) == 29 ) energy *= HCAL_Calib_29;
00394                 //if ( rescaleFactor > 1. ) 
00395                 // std::cout << "End-cap HCAL energy rescaled from = " << energy << " to " << energy*rescaleFactor << std::endl;
00396                 if ( rescaleFactor > 1. ) { 
00397                   pfrhCleaned = createHcalRecHit( detid, 
00398                                                   energy, 
00399                                                   PFLayer::HCAL_ENDCAP, 
00400                                                   hcalEndcapGeometry,
00401                                                   ct.id().rawId() );
00402                   pfrhCleaned->setRescale(rescaleFactor);
00403                   energy *= rescaleFactor;
00404                 }
00405                 pfrh = createHcalRecHit( detid, 
00406                                          energy, 
00407                                          PFLayer::HCAL_ENDCAP, 
00408                                          hcalEndcapGeometry,
00409                                          ct.id().rawId() );
00410                 pfrh->setRescale(rescaleFactor);
00411               }
00412               break;
00413             case HcalOuter:
00414               {
00415               }
00416               break;
00417             case HcalForward:
00418               {
00419                 //---ab: 2 rechits for HF:
00420                 //double energyemHF = weight_HFem_*ct.emEnergy();
00421                 //double energyhadHF = weight_HFhad_*ct.hadEnergy();
00422                 double energyemHF = weight_HFem_ * energyEM;
00423                 double energyhadHF = weight_HFhad_ * energy;
00424                 // Some energy in the tower !
00425                 if((energyemHF+energyhadHF) < thresh_HF_ ) continue;
00426 
00427                 // Some cleaning in the HF 
00428                 double longFibre = energyemHF + energyhadHF/2.;
00429                 double shortFibre = energyhadHF/2.;
00430                 int ieta = detid.ieta();
00431                 int iphi = detid.iphi();
00432                 HcalDetId theLongDetId (HcalForward, ieta, iphi, 1);
00433                 HcalDetId theShortDetId (HcalForward, ieta, iphi, 2);
00434                 typedef HFRecHitCollection::const_iterator iHF;
00435                 iHF theLongHit = hfHandle->find(theLongDetId); 
00436                 iHF theShortHit = hfHandle->find(theShortDetId); 
00437                 // 
00438                 double theLongHitEnergy = 0.;
00439                 double theShortHitEnergy = 0.;
00440                 bool flagShortDPG =  false; 
00441                 bool flagLongDPG = false; 
00442                 bool flagShortTimeDPG = false; 
00443                 bool flagLongTimeDPG = false;
00444                 bool flagShortPulseDPG = false;
00445                 bool flagLongPulseDPG = false;
00446                 //
00447                 if ( theLongHit != hfHandle->end() ) { 
00448                   int theLongFlag = theLongHit->flags();
00449                   theLongHitEnergy = theLongHit->energy();
00450                   flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
00451                   flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
00452                   flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
00453 
00454                   //flagLongDPG =  applyLongShortDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
00455                   //flagLongTimeDPG = applyTimeDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
00456                   //flagLongPulseDPG = applyPulseDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
00457                 }
00458                 //
00459                 if ( theShortHit != hfHandle->end() ) { 
00460                   int theShortFlag = theShortHit->flags();
00461                   theShortHitEnergy = theShortHit->energy();
00462                   flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
00463                   flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
00464                   flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
00465                   //flagShortDPG =  applyLongShortDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
00466                   //flagShortTimeDPG = applyTimeDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
00467                   //flagShortPulseDPG = applyPulseDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
00468                 }
00469 
00470                 // Then check the timing in short and long fibres in all other towers.
00471                 if ( theShortHitEnergy > longShortFibre_Cut && 
00472                      ( theShortHit->time() < minShortTiming_Cut ||
00473                        theShortHit->time() > maxShortTiming_Cut || 
00474                        flagShortTimeDPG || flagShortPulseDPG ) ) { 
00475                   // rescaleFactor = 0. ;
00476                   pfrhHFHADCleaned = createHcalRecHit( detid, 
00477                                                        theShortHitEnergy, 
00478                                                        PFLayer::HF_HAD, 
00479                                                        hcalEndcapGeometry,
00480                                                        ct.id().rawId() );
00481                   pfrhHFHADCleaned->setRescale(theShortHit->time());
00482                   /*
00483                   std::cout << "ieta/iphi = " << ieta << " " << iphi 
00484                             << ", Energy em/had/long/short = " 
00485                             << energyemHF << " " << energyhadHF << " "
00486                             << theLongHitEnergy << " " << theShortHitEnergy << " " 
00487                             << ". Time = " << theShortHit->time() << " " 
00488                             << ". The short and long flags : " 
00489                             << flagShortDPG << " " << flagLongDPG << " "  
00490                             << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
00491                             << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
00492                             << ". Short fibres were cleaned." << std::endl;
00493                   */
00494                   shortFibre -= theShortHitEnergy;
00495                   theShortHitEnergy = 0.;
00496                 }
00497                 
00498                 
00499                 if ( theLongHitEnergy > longShortFibre_Cut && 
00500                      ( theLongHit->time() < minLongTiming_Cut ||
00501                        theLongHit->time() > maxLongTiming_Cut  || 
00502                        flagLongTimeDPG || flagLongPulseDPG ) ) { 
00503                   //rescaleFactor = 0. ;
00504                   pfrhHFEMCleaned = createHcalRecHit( detid, 
00505                                                       theLongHitEnergy, 
00506                                                       PFLayer::HF_EM, 
00507                                                       hcalEndcapGeometry,
00508                                                       ct.id().rawId() );
00509                   pfrhHFEMCleaned->setRescale(theLongHit->time());
00510                   /*
00511                   std::cout << "ieta/iphi = " << ieta << " " << iphi 
00512                             << ", Energy em/had/long/short = " 
00513                             << energyemHF << " " << energyhadHF << " "
00514                             << theLongHitEnergy << " " << theShortHitEnergy << " " 
00515                             << ". Time = " << theLongHit->time() << " " 
00516                             << ". The short and long flags : " 
00517                             << flagShortDPG << " " << flagLongDPG << " "  
00518                             << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
00519                             << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
00520                             << ". Long fibres were cleaned." << std::endl;
00521                   */
00522                   longFibre -= theLongHitEnergy;
00523                   theLongHitEnergy = 0.;
00524                 }
00525 
00526                 // Some energy must be in the long fibres is there is some energy in the short fibres ! 
00527                 if ( theShortHitEnergy > shortFibre_Cut && 
00528                      ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction || 
00529                        flagShortDPG ) ) {
00530                   // Check if the long-fibre hit was not cleaned already (because hot)
00531                   // In this case don't apply the cleaning
00532                   const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId);
00533                   unsigned theStatusValue = theStatus->getValue();
00534                   int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
00535                   // The channel is killed
00537                   if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
00538                     // rescaleFactor = 0. ;
00539                     pfrhHFHADCleaned = createHcalRecHit( detid, 
00540                                                          theShortHitEnergy, 
00541                                                          PFLayer::HF_HAD, 
00542                                                          hcalEndcapGeometry,
00543                                                          ct.id().rawId() );
00544                     pfrhHFHADCleaned->setRescale(theShortHit->time());
00545                     /*
00546                     std::cout << "ieta/iphi = " << ieta << " " << iphi 
00547                               << ", Energy em/had/long/short = " 
00548                               << energyemHF << " " << energyhadHF << " "
00549                               << theLongHitEnergy << " " << theShortHitEnergy << " " 
00550                               << ". Time = " << theShortHit->time() << " " 
00551                               << ". The status value is " << theStatusValue
00552                               << ". The short and long flags : " 
00553                               << flagShortDPG << " " << flagLongDPG << " "  
00554                               << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
00555                               << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
00556                               << ". Short fibres were cleaned." << std::endl;
00557                     */
00558                     shortFibre -= theShortHitEnergy;
00559                     theShortHitEnergy = 0.;
00560                   }
00561                 }
00562 
00563                 if ( theLongHitEnergy > longFibre_Cut && 
00564                      ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction || 
00565                        flagLongDPG ) ) {
00566                   // Check if the long-fibre hit was not cleaned already (because hot)
00567                   // In this case don't apply the cleaning
00568                   const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId);
00569                   unsigned theStatusValue = theStatus->getValue();
00570 
00571                   int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
00572                   // The channel is killed
00574                   if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
00575                     
00576                     //rescaleFactor = 0. ;
00577                     pfrhHFEMCleaned = createHcalRecHit( detid, 
00578                                                       theLongHitEnergy, 
00579                                                       PFLayer::HF_EM, 
00580                                                       hcalEndcapGeometry,
00581                                                       ct.id().rawId() );
00582                     pfrhHFEMCleaned->setRescale(theLongHit->time());
00583                     /*
00584                     std::cout << "ieta/iphi = " << ieta << " " << iphi 
00585                               << ", Energy em/had/long/short = " 
00586                               << energyemHF << " " << energyhadHF << " "
00587                               << theLongHitEnergy << " " << theShortHitEnergy << " " 
00588                               << ". The status value is " << theStatusValue
00589                               << ". Time = " << theLongHit->time() << " " 
00590                               << ". The short and long flags : " 
00591                               << flagShortDPG << " " << flagLongDPG << " "  
00592                               << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
00593                               << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
00594                               << ". Long fibres were cleaned." << std::endl;
00595                     */
00596                     longFibre -= theLongHitEnergy;
00597                     theLongHitEnergy = 0.;
00598                   }
00599                 }
00600 
00601                 // Special treatment for tower 29
00602                 // A tower with energy only at ieta = +/- 29 is not physical -> Clean
00603                 if ( abs(ieta) == 29 ) { 
00604                   // rescaleFactor = 0. ;
00605                   // Clean long fibres
00606                   if ( theLongHitEnergy > shortFibre_Cut/2. ) { 
00607                     pfrhHFEMCleaned29 = createHcalRecHit( detid, 
00608                                                           theLongHitEnergy, 
00609                                                           PFLayer::HF_EM, 
00610                                                           hcalEndcapGeometry,
00611                                                           ct.id().rawId() );
00612                     pfrhHFEMCleaned29->setRescale(theLongHit->time());
00613                     /*
00614                     std::cout << "ieta/iphi = " << ieta << " " << iphi 
00615                               << ", Energy em/had/long/short = " 
00616                               << energyemHF << " " << energyhadHF << " "
00617                               << theLongHitEnergy << " " << theShortHitEnergy << " " 
00618                               << ". Long fibres were cleaned." << std::endl;
00619                     */
00620                     longFibre -= theLongHitEnergy;
00621                     theLongHitEnergy = 0.;
00622                   }
00623                   // Clean short fibres
00624                   if ( theShortHitEnergy > shortFibre_Cut/2. ) { 
00625                     pfrhHFHADCleaned29 = createHcalRecHit( detid, 
00626                                                            theShortHitEnergy, 
00627                                                            PFLayer::HF_HAD, 
00628                                                            hcalEndcapGeometry,
00629                                                            ct.id().rawId() );
00630                     pfrhHFHADCleaned29->setRescale(theShortHit->time());
00631                     /*
00632                     std::cout << "ieta/iphi = " << ieta << " " << iphi 
00633                               << ", Energy em/had/long/short = " 
00634                               << energyemHF << " " << energyhadHF << " "
00635                               << theLongHitEnergy << " " << theShortHitEnergy << " " 
00636                               << ". Short fibres were cleaned." << std::endl;
00637                     */
00638                     shortFibre -= theShortHitEnergy;
00639                     theShortHitEnergy = 0.;
00640                   }
00641                 }
00642                 // Check the timing of the long and short fibre rechits
00643                 
00644                 // First, check the timing of long and short fibre in eta = 29 if tower 30.
00645                 else if ( abs(ieta) == 30 ) { 
00646                   int ieta29 = ieta > 0 ? 29 : -29;
00647                   HcalDetId theLongDetId29 (HcalForward, ieta29, iphi, 1);
00648                   HcalDetId theShortDetId29 (HcalForward, ieta29, iphi, 2);
00649                   iHF theLongHit29 = hfHandle->find(theLongDetId29); 
00650                   iHF theShortHit29 = hfHandle->find(theShortDetId29); 
00651                   // 
00652                   double theLongHitEnergy29 = 0.;
00653                   double theShortHitEnergy29 = 0.;
00654                   bool flagShortDPG29 =  false; 
00655                   bool flagLongDPG29 = false; 
00656                   bool flagShortTimeDPG29 = false; 
00657                   bool flagLongTimeDPG29 = false;
00658                   bool flagShortPulseDPG29 = false;
00659                   bool flagLongPulseDPG29 = false;
00660                   //
00661                   if ( theLongHit29 != hfHandle->end() ) { 
00662                     int theLongFlag29 = theLongHit29->flags();
00663                     theLongHitEnergy29 = theLongHit29->energy() ;
00664                     flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
00665                     flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
00666                     flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
00667 
00668                     //flagLongDPG29 = applyLongShortDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
00669                     //flagLongTimeDPG29 = applyTimeDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
00670                     //flagLongPulseDPG29 = applyPulseDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
00671                   }
00672                   //
00673                   if ( theShortHit29 != hfHandle->end() ) { 
00674                     int theShortFlag29 = theShortHit29->flags();
00675                     theShortHitEnergy29 = theShortHit29->energy();                
00676                     flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
00677                     flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
00678                     flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
00679 
00680                     //flagShortDPG29 = applyLongShortDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
00681                     //flagShortTimeDPG29 = applyTimeDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
00682                     //flagShortPulseDPG29 = applyPulseDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
00683                   }
00684 
00685                   if ( theLongHitEnergy29 > longShortFibre_Cut && 
00686                        ( theLongHit29->time() < minLongTiming_Cut ||
00687                          theLongHit29->time() > maxLongTiming_Cut ||
00688                          flagLongTimeDPG29 || flagLongPulseDPG29 ) ) { 
00689                     //rescaleFactor = 0. ;
00690                     pfrhHFEMCleaned29 = createHcalRecHit( detid, 
00691                                                           theLongHitEnergy29, 
00692                                                           PFLayer::HF_EM, 
00693                                                           hcalEndcapGeometry,
00694                                                           ct.id().rawId() );
00695                     pfrhHFEMCleaned29->setRescale(theLongHit29->time());
00696                     /*
00697                     std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
00698                               << ", Energy em/had/long/short = " 
00699                               << energyemHF << " " << energyhadHF << " "
00700                               << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
00701                               << ". Time = " << theLongHit29->time() << " " 
00702                               << ". The short and long flags : " 
00703                               << flagShortDPG29 << " " << flagLongDPG29 << " "  
00704                               << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
00705                               << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
00706                               << ". Long fibres were cleaned." << std::endl;
00707                     */
00708                     longFibre -= theLongHitEnergy29;
00709                     theLongHitEnergy29 = 0;
00710                   }
00711 
00712                   if ( theShortHitEnergy29 > longShortFibre_Cut && 
00713                        ( theShortHit29->time() < minShortTiming_Cut ||
00714                          theShortHit29->time() > maxShortTiming_Cut ||
00715                          flagShortTimeDPG29 || flagShortPulseDPG29 ) ) { 
00716                     //rescaleFactor = 0. ;
00717                     pfrhHFHADCleaned29 = createHcalRecHit( detid, 
00718                                                          theShortHitEnergy29, 
00719                                                          PFLayer::HF_HAD, 
00720                                                          hcalEndcapGeometry,
00721                                                          ct.id().rawId() );
00722                     pfrhHFHADCleaned29->setRescale(theShortHit29->time());
00723                     /*
00724                     std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
00725                               << ", Energy em/had/long/short = " 
00726                               << energyemHF << " " << energyhadHF << " "
00727                               << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
00728                               << ". Time = " << theShortHit29->time() << " " 
00729                               << ". The short and long flags : " 
00730                               << flagShortDPG29 << " " << flagLongDPG29 << " "  
00731                               << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
00732                               << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
00733                               << ". Short fibres were cleaned." << std::endl;
00734                     */
00735                     shortFibre -= theShortHitEnergy29;
00736                     theShortHitEnergy29 = 0.;
00737                   }
00738 
00739                   // Some energy must be in the long fibres is there is some energy in the short fibres ! 
00740                   if ( theShortHitEnergy29 > shortFibre_Cut && 
00741                        ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction || 
00742                          flagShortDPG29 ) ) {
00743                     // Check if the long-fibre hit was not cleaned already (because hot)
00744                     // In this case don't apply the cleaning
00745                     const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId29);
00746                     unsigned theStatusValue = theStatus->getValue();
00747 
00748                     int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
00749                     // The channel is killed
00751                     if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
00752                       //rescaleFactor = 0. ;
00753                       pfrhHFHADCleaned29 = createHcalRecHit( detid, 
00754                                                            theShortHitEnergy29, 
00755                                                            PFLayer::HF_HAD, 
00756                                                            hcalEndcapGeometry,
00757                                                            ct.id().rawId() );
00758                       pfrhHFHADCleaned29->setRescale(theShortHit29->time());
00759                       /*
00760                       std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
00761                                 << ", Energy em/had/long/short = " 
00762                                 << energyemHF << " " << energyhadHF << " "
00763                                 << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
00764                                 << ". Time = " << theShortHit29->time() << " " 
00765                                 << ". The status value is " << theStatusValue
00766                                 << ". The short and long flags : " 
00767                                 << flagShortDPG29 << " " << flagLongDPG29 << " "  
00768                                 << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
00769                                 << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
00770                                 << ". Short fibres were cleaned." << std::endl;
00771                       */
00772                       shortFibre -= theShortHitEnergy29;
00773                       theShortHitEnergy29 = 0.;
00774                     }       
00775                   }
00776                   
00777                   // Some energy must be in the short fibres is there is some energy in the long fibres ! 
00778                   if ( theLongHitEnergy29 > longFibre_Cut && 
00779                        ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction || 
00780                          flagLongDPG29 ) ) {
00781                     // Check if the long-fibre hit was not cleaned already (because hot)
00782                     // In this case don't apply the cleaning
00783                     const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
00784                     unsigned theStatusValue = theStatus->getValue();
00785                     int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
00786                     // The channel is killed
00788                     if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
00789 
00790                       //rescaleFactor = 0. ;
00791                       pfrhHFEMCleaned29 = createHcalRecHit( detid, 
00792                                                             theLongHitEnergy29, 
00793                                                             PFLayer::HF_EM, 
00794                                                             hcalEndcapGeometry,
00795                                                             ct.id().rawId() );
00796                       pfrhHFEMCleaned29->setRescale(theLongHit29->time());
00797                       /* 
00798                       std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
00799                                 << ", Energy em/had/long/short = " 
00800                                 << energyemHF << " " << energyhadHF << " "
00801                                 << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
00802                                 << ". The status value is " << theStatusValue
00803                                 << ". Time = " << theLongHit29->time() << " " 
00804                                 << ". The short and long flags : " 
00805                                 << flagShortDPG29 << " " << flagLongDPG29 << " "  
00806                                 << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
00807                                 << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
00808                                 << ". Long fibres were cleaned." << std::endl;
00809                       */
00810                       longFibre -= theLongHitEnergy29;
00811                       theLongHitEnergy29 = 0.;
00812                     }
00813                   }
00814 
00815                   // Check that the energy in tower 29 is smaller than in tower 30
00816                   // First in long fibres
00817                   if ( theLongHitEnergy29 > std::max(theLongHitEnergy,shortFibre_Cut/2) ) { 
00818                     pfrhHFEMCleaned29 = createHcalRecHit( detid, 
00819                                                           theLongHitEnergy29, 
00820                                                           PFLayer::HF_EM, 
00821                                                           hcalEndcapGeometry,
00822                                                           ct.id().rawId() );
00823                     pfrhHFEMCleaned29->setRescale(theLongHit29->time());
00824                     /*
00825                     std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
00826                               << ", Energy L29/S29/L30/S30 = " 
00827                               << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
00828                               << theLongHitEnergy << " " << theShortHitEnergy << " " 
00829                               << ". Long fibres were cleaned." << std::endl;
00830                     */
00831                     longFibre -= theLongHitEnergy29;
00832                     theLongHitEnergy29 = 0.;
00833                   }
00834                   // Second in short fibres
00835                   if ( theShortHitEnergy29 > std::max(theShortHitEnergy,shortFibre_Cut/2.) ) { 
00836                     pfrhHFHADCleaned29 = createHcalRecHit( detid, 
00837                                                            theShortHitEnergy29, 
00838                                                            PFLayer::HF_HAD, 
00839                                                            hcalEndcapGeometry,
00840                                                            ct.id().rawId() );
00841                     pfrhHFHADCleaned29->setRescale(theShortHit29->time());
00842                     /*
00843                     std::cout << "ieta/iphi = " << ieta << " " << iphi 
00844                               << ", Energy L29/S29/L30/S30 = " 
00845                               << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
00846                               << theLongHitEnergy << " " << theShortHitEnergy << " " 
00847                               << ". Short fibres were cleaned." << std::endl;
00848                     */
00849                     shortFibre -= theShortHitEnergy29;
00850                     theShortHitEnergy29 = 0.;
00851                   }
00852                 }
00853 
00854 
00855                 // Determine EM and HAD after cleaning of short and long fibres
00856                 energyhadHF = 2.*shortFibre;
00857                 energyemHF = longFibre - shortFibre;
00858 
00859                 // The EM energy might be negative, as it amounts to Long - Short
00860                 // In that case, put the EM "energy" in the HAD energy
00861                 // Just to avoid systematic positive bias due to "Short" high fluctuations
00862                 if ( energyemHF < thresh_HF_ ) { 
00863                   energyhadHF += energyemHF;
00864                   energyemHF = 0.;
00865                 }
00866 
00867                 // Apply HCAL calibration factors flor towers close to 29, if requested
00868                 if ( HF_Calib_ && abs(detid.ieta()) <= 32 ) { 
00869                   energyhadHF *= HF_Calib_29;
00870                   energyemHF *= HF_Calib_29;
00871                 }
00872                                 
00873                 // Create an EM and a HAD rechit if above threshold.
00874                 if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) { 
00875                   pfrhHFEM = createHcalRecHit( detid, 
00876                                                energyemHF, 
00877                                                PFLayer::HF_EM, 
00878                                                hcalEndcapGeometry,
00879                                                ct.id().rawId() );
00880                   pfrhHFHAD = createHcalRecHit( detid, 
00881                                                 energyhadHF, 
00882                                                 PFLayer::HF_HAD, 
00883                                                 hcalEndcapGeometry,
00884                                                 ct.id().rawId() );
00885                   pfrhHFEM->setEnergyUp(energyhadHF);
00886                   pfrhHFHAD->setEnergyUp(energyemHF);
00887                 }
00888                 
00889               }
00890               break;
00891             default:
00892               LogError("PFRecHitProducerHCAL")
00893                 <<"CaloTower constituent: unknown layer : "
00894                 <<detid.subdet()<<endl;
00895             } 
00896 
00897             if(pfrh) { 
00898               rechits.push_back( *pfrh );
00899               delete pfrh;
00900               idSortedRecHits.insert( make_pair(ct.id().rawId(), 
00901                                                 rechits.size()-1 ) ); 
00902             }
00903             if(pfrhCleaned) { 
00904               rechitsCleaned.push_back( *pfrhCleaned );
00905               delete pfrhCleaned;
00906             }
00907             //---ab: 2 rechits for HF:     
00908             if(pfrhHFEM) { 
00909               HFEMRecHits->push_back( *pfrhHFEM );
00910               delete pfrhHFEM;
00911               idSortedRecHitsHFEM.insert( make_pair(ct.id().rawId(), 
00912                                                 HFEMRecHits->size()-1 ) ); 
00913             }
00914             if(pfrhHFHAD) { 
00915               HFHADRecHits->push_back( *pfrhHFHAD );
00916               delete pfrhHFHAD;
00917               idSortedRecHitsHFHAD.insert( make_pair(ct.id().rawId(), 
00918                                                 HFHADRecHits->size()-1 ) ); 
00919             }
00920             //---ab        
00921             if(pfrhHFEMCleaned) { 
00922               rechitsCleaned.push_back( *pfrhHFEMCleaned );
00923               delete pfrhHFEMCleaned;
00924             }
00925             if(pfrhHFHADCleaned) { 
00926               rechitsCleaned.push_back( *pfrhHFHADCleaned );
00927               delete pfrhHFHADCleaned;
00928             }
00929             if(pfrhHFEMCleaned29) { 
00930               rechitsCleaned.push_back( *pfrhHFEMCleaned29 );
00931               delete pfrhHFEMCleaned29;
00932             }
00933             if(pfrhHFHADCleaned29) { 
00934               rechitsCleaned.push_back( *pfrhHFHADCleaned29 );
00935               delete pfrhHFHADCleaned29;
00936             }
00937           }
00938       }
00939       // do navigation 
00940       for(unsigned i=0; i<rechits.size(); i++ ) {
00941         findRecHitNeighboursCT( rechits[i], 
00942                                 idSortedRecHits, 
00943                                 caloTowerTopology);
00944       }
00945       for(unsigned i=0; i<HFEMRecHits->size(); i++ ) {
00946         findRecHitNeighboursCT( (*HFEMRecHits)[i], 
00947                                 idSortedRecHitsHFEM, 
00948                                 caloTowerTopology);
00949       }
00950       for(unsigned i=0; i<HFHADRecHits->size(); i++ ) {
00951         findRecHitNeighboursCT( (*HFHADRecHits)[i], 
00952                                 idSortedRecHitsHFHAD, 
00953                                 caloTowerTopology);
00954       }
00955       iEvent.put( HFHADRecHits,"HFHAD" );       
00956       iEvent.put( HFEMRecHits,"HFEM" ); 
00957     }   
00958   }
00959   else if( !(inputTagHcalRecHitsHBHE_ == InputTag()) ) { 
00960     // clustering is not done on CaloTowers but on HCAL rechits.
00961        
00962   
00963     // get the hcal topology
00964     edm::ESHandle<HcalTopology> hcalTopology;
00965     iSetup.get<IdealGeometryRecord>().get( hcalTopology );
00966     
00967     // HCAL rechits 
00968     //    vector<edm::Handle<HBHERecHitCollection> > hcalHandles;  
00969     edm::Handle<HBHERecHitCollection>  hcalHandle;  
00970 
00971     
00972     bool found = iEvent.getByLabel(inputTagHcalRecHitsHBHE_, 
00973                                    hcalHandle );
00974 
00975     if(!found) {
00976       ostringstream err;
00977       err<<"could not find rechits "<<inputTagHcalRecHitsHBHE_;
00978       LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
00979     
00980       throw cms::Exception( "MissingProduct", err.str());
00981     }
00982     else {
00983       assert( hcalHandle.isValid() );
00984       
00985       const edm::Handle<HBHERecHitCollection>& handle = hcalHandle;
00986       for(unsigned irechit=0; irechit<handle->size(); irechit++) {
00987         const HBHERecHit& hit = (*handle)[irechit];
00988         
00989         double energy = hit.energy();
00990         
00991         reco::PFRecHit* pfrh = 0;
00992         
00993 
00994         const HcalDetId& detid = hit.detid();
00995         switch( detid.subdet() ) {
00996         case HcalBarrel:
00997           {
00998             if(energy < thresh_Barrel_ ) continue;
00999             pfrh = createHcalRecHit( detid, 
01000                                      energy, 
01001                                      PFLayer::HCAL_BARREL1, 
01002                                      hcalBarrelGeometry );
01003           }
01004           break;
01005         case HcalEndcap:
01006           {
01007             if(energy < thresh_Endcap_ ) continue;
01008             pfrh = createHcalRecHit( detid, 
01009                                      energy, 
01010                                      PFLayer::HCAL_ENDCAP, 
01011                                      hcalEndcapGeometry );        
01012           }
01013           break;
01014         case HcalForward:
01015           {
01016             if(energy < thresh_HF_ ) continue;
01017             pfrh = createHcalRecHit( detid, 
01018                                      energy, 
01019                                      PFLayer::HF_HAD, 
01020                                      hcalEndcapGeometry );
01021           }
01022           break;
01023         default:
01024           LogError("PFRecHitProducerHCAL")
01025             <<"HCAL rechit: unknown layer : "<<detid.subdet()<<endl;
01026           continue;
01027         } 
01028 
01029         if(pfrh) { 
01030           rechits.push_back( *pfrh );
01031           delete pfrh;
01032           idSortedRecHits.insert( make_pair(detid.rawId(), 
01033                                             rechits.size()-1 ) ); 
01034         }
01035       }
01036       
01037       
01038       // do navigation:
01039       for(unsigned i=0; i<rechits.size(); i++ ) {
01040         
01041         findRecHitNeighbours( rechits[i], idSortedRecHits, 
01042                               *hcalTopology, 
01043                               *hcalBarrelGeometry, 
01044                               *hcalTopology,
01045                               *hcalEndcapGeometry);
01046       } // loop for navigation
01047     }  // endif hcal rechits were found
01048   } // endif clustering on rechits in hcal
01049 }
01050 
01051 
01052 
01053 
01054 
01055 
01056 reco::PFRecHit* 
01057 PFRecHitProducerHCAL::createHcalRecHit( const DetId& detid,
01058                                         double energy,
01059                                         PFLayer::Layer layer,
01060                                         const CaloSubdetectorGeometry* geom,
01061                                         unsigned newDetId ) {
01062   
01063   const CaloCellGeometry *thisCell = geom->getGeometry(detid);
01064   if(!thisCell) {
01065     edm::LogError("PFRecHitProducerHCAL")
01066       <<"warning detid "<<detid.rawId()<<" not found in layer "
01067       <<layer<<endl;
01068     return 0;
01069   }
01070   
01071   const GlobalPoint& position = thisCell->getPosition();
01072   
01073   double depth_correction = 0.;
01074   switch ( layer ) { 
01075   case PFLayer::HF_EM:
01076     depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
01077     break;
01078   case PFLayer::HF_HAD:
01079     depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
01080     break;
01081   default:
01082     break;
01083   }
01084 
01085   unsigned id = detid;
01086   if(newDetId) id = newDetId;
01087   reco::PFRecHit *rh = 
01088     new reco::PFRecHit( id,  layer, energy, 
01089                         position.x(), position.y(), position.z()+depth_correction, 
01090                         0,0,0 );
01091  
01092   
01093   
01094   
01095   // set the corners
01096   const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
01097 
01098   assert( corners.size() == 8 );
01099 
01100   rh->setNECorner( corners[0].x(), corners[0].y(),  corners[0].z()+depth_correction );
01101   rh->setSECorner( corners[1].x(), corners[1].y(),  corners[1].z()+depth_correction );
01102   rh->setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z()+depth_correction );
01103   rh->setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z()+depth_correction );
01104  
01105   return rh;
01106 }
01107 
01108 
01109 
01110 
01111 void 
01112 PFRecHitProducerHCAL::findRecHitNeighbours
01113 ( reco::PFRecHit& rh, 
01114   const map<unsigned,unsigned >& sortedHits, 
01115   const CaloSubdetectorTopology& barrelTopology, 
01116   const CaloSubdetectorGeometry& barrelGeometry, 
01117   const CaloSubdetectorTopology& endcapTopology, 
01118   const CaloSubdetectorGeometry& endcapGeometry ) {
01119   
01120   //cout<<"------PFRecHitProducerHcaL:findRecHitNeighbours navigation value "<<navigation_HF_<<endl;
01121  if(navigation_HF_ == false){
01122     if( rh.layer() == PFLayer::HF_HAD )
01123       return;
01124     if( rh.layer() == PFLayer::HF_EM )
01125       return;
01126   } 
01127   DetId detid( rh.detId() );
01128 
01129   const CaloSubdetectorTopology* topology = 0;
01130   const CaloSubdetectorGeometry* geometry = 0;
01131   
01132   switch( rh.layer() ) {
01133   case PFLayer::ECAL_ENDCAP: 
01134     topology = &endcapTopology;
01135     geometry = &endcapGeometry;
01136     break;
01137   case PFLayer::ECAL_BARREL: 
01138     topology = &barrelTopology;
01139     geometry = &barrelGeometry;
01140     break;
01141   case PFLayer::HCAL_ENDCAP:
01142     topology = &endcapTopology;
01143     geometry = &endcapGeometry;
01144     break;
01145   case PFLayer::HCAL_BARREL1:
01146     topology = &barrelTopology;
01147     geometry = &barrelGeometry;
01148     break;
01149   case PFLayer::PS1:
01150   case PFLayer::PS2:
01151     topology = &barrelTopology;
01152     geometry = &barrelGeometry;
01153     break;
01154   default:
01155     assert(0);
01156   }
01157   
01158   assert( topology && geometry );
01159 
01160   CaloNavigator<DetId> navigator(detid, topology);
01161 
01162   DetId north = navigator.north();  
01163   
01164   DetId northeast(0);
01165   if( north != DetId(0) ) {
01166     northeast = navigator.east();  
01167   }
01168   navigator.home();
01169 
01170 
01171   DetId south = navigator.south();
01172 
01173   
01174 
01175   DetId southwest(0); 
01176   if( south != DetId(0) ) {
01177     southwest = navigator.west();
01178   }
01179   navigator.home();
01180 
01181 
01182   DetId east = navigator.east();
01183   DetId southeast;
01184   if( east != DetId(0) ) {
01185     southeast = navigator.south(); 
01186   }
01187   navigator.home();
01188   DetId west = navigator.west();
01189   DetId northwest;
01190   if( west != DetId(0) ) {   
01191     northwest = navigator.north();  
01192   }
01193   navigator.home();
01194     
01195   IDH i = sortedHits.find( north.rawId() );
01196   if(i != sortedHits.end() ) 
01197     rh.add4Neighbour( i->second );
01198   
01199   i = sortedHits.find( northeast.rawId() );
01200   if(i != sortedHits.end() ) 
01201     rh.add8Neighbour( i->second );
01202   
01203   i = sortedHits.find( south.rawId() );
01204   if(i != sortedHits.end() ) 
01205     rh.add4Neighbour( i->second );
01206     
01207   i = sortedHits.find( southwest.rawId() );
01208   if(i != sortedHits.end() ) 
01209     rh.add8Neighbour( i->second );
01210     
01211   i = sortedHits.find( east.rawId() );
01212   if(i != sortedHits.end() ) 
01213     rh.add4Neighbour( i->second );
01214     
01215   i = sortedHits.find( southeast.rawId() );
01216   if(i != sortedHits.end() ) 
01217     rh.add8Neighbour( i->second );
01218     
01219   i = sortedHits.find( west.rawId() );
01220   if(i != sortedHits.end() ) 
01221      rh.add4Neighbour( i->second );
01222    
01223   i = sortedHits.find( northwest.rawId() );
01224   if(i != sortedHits.end() ) 
01225     rh.add8Neighbour( i->second );
01226     
01227 
01228 }
01229 
01230 
01231 void 
01232 PFRecHitProducerHCAL::findRecHitNeighboursCT
01233 ( reco::PFRecHit& rh, 
01234   const map<unsigned, unsigned >& sortedHits, 
01235   const CaloSubdetectorTopology& topology ) {
01236   //cout<<"------PFRecHitProducerHcaL:findRecHitNeighboursCT navigation value "<<navigation_HF_<<endl;
01237   //  cout<<"----------- rechit print out"<<endl;
01238   // if(( rh.layer() == PFLayer::HF_HAD )||(rh.layer() == PFLayer::HF_EM)) {  
01239     
01240   //    cout<<rh<<endl;
01241     //  }
01242   if(navigation_HF_ == false){
01243     if( rh.layer() == PFLayer::HF_HAD )
01244       return;
01245     if( rh.layer() == PFLayer::HF_EM )
01246       return;
01247   }
01248   CaloTowerDetId ctDetId( rh.detId() );
01249     
01250 
01251   vector<DetId> northids = topology.north(ctDetId);
01252   vector<DetId> westids = topology.west(ctDetId);
01253   vector<DetId> southids = topology.south(ctDetId);
01254   vector<DetId> eastids = topology.east(ctDetId);
01255 
01256 
01257   CaloTowerDetId badId;
01258 
01259   // all the following detids will be CaloTowerDetId
01260   CaloTowerDetId north;
01261   CaloTowerDetId northwest;
01262   CaloTowerDetId northwest2;
01263   CaloTowerDetId west;
01264   CaloTowerDetId west2;
01265   CaloTowerDetId southwest;
01266   CaloTowerDetId southwest2;
01267   CaloTowerDetId south;
01268   CaloTowerDetId southeast;
01269   CaloTowerDetId southeast2;
01270   CaloTowerDetId east;
01271   CaloTowerDetId east2;
01272   CaloTowerDetId northeast;
01273   CaloTowerDetId northeast2;
01274   
01275   // for north and south, there is no ambiguity : 1 or 0 neighbours
01276   
01277   switch( northids.size() ) {
01278   case 0: 
01279     break;
01280   case 1: 
01281     north = northids[0];
01282     break;
01283   default:
01284   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours north: "); 
01285     err<<northids.size();
01286     throw( err.str() ); 
01287   }
01288 
01289   switch( southids.size() ) {
01290   case 0: 
01291     break;
01292   case 1: 
01293     south = southids[0];
01294     break;
01295   default:
01296   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours south: "); 
01297     err<<southids.size();
01298     throw( err.str() ); 
01299   }
01300   
01301   // for east and west, one must take care 
01302   // of the pitch change in HCAL endcap.
01303 
01304   switch( eastids.size() ) {
01305   case 0: 
01306     break;
01307   case 1: 
01308     east = eastids[0];
01309     northeast = getNorth(east, topology);
01310     southeast = getSouth(east, topology);
01311     break;
01312   case 2:  
01313     // in this case, 0 is more on the north than 1
01314     east = eastids[0];
01315     east2 = eastids[1];
01316     northeast = getNorth(east, topology );
01317     southeast = getSouth(east2, topology);    
01318     northeast2 = getNorth(northeast, topology );
01319     southeast2 = getSouth(southeast, topology);    
01320     break;
01321   default:
01322   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours eastids: "); 
01323     err<<eastids.size();
01324     throw( err.str() ); 
01325   }
01326   
01327   
01328   switch( westids.size() ) {
01329   case 0: 
01330     break;
01331   case 1: 
01332     west = westids[0];
01333     northwest = getNorth(west, topology);
01334     southwest = getSouth(west, topology);
01335     break;
01336   case 2:  
01337     // in this case, 0 is more on the north than 1
01338     west = westids[0];
01339     west2 = westids[1];
01340     northwest = getNorth(west, topology );
01341     southwest = getSouth(west2, topology );    
01342     northwest2 = getNorth(northwest, topology );
01343     southwest2 = getSouth(southwest, topology );    
01344     break;
01345   default:
01346   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours westids: "); 
01347     err<< westids.size();
01348     throw( err.str() ); 
01349   }
01350 
01351 
01352 
01353 
01354   // find and set neighbours
01355     
01356   IDH i = sortedHits.find( north.rawId() );
01357   if(i != sortedHits.end() ) 
01358     rh.add4Neighbour( i->second );
01359   
01360   i = sortedHits.find( northeast.rawId() );
01361   if(i != sortedHits.end() ) 
01362     rh.add8Neighbour( i->second );
01363   
01364   i = sortedHits.find( northeast2.rawId() );
01365   if(i != sortedHits.end() ) 
01366     rh.add8Neighbour( i->second );
01367   
01368   i = sortedHits.find( south.rawId() );
01369   if(i != sortedHits.end() ) 
01370     rh.add4Neighbour( i->second );
01371     
01372   i = sortedHits.find( southwest.rawId() );
01373   if(i != sortedHits.end() ) 
01374     rh.add8Neighbour( i->second );
01375     
01376   i = sortedHits.find( southwest2.rawId() );
01377   if(i != sortedHits.end() ) 
01378     rh.add8Neighbour( i->second );
01379     
01380   i = sortedHits.find( east.rawId() );
01381   if(i != sortedHits.end() ) 
01382     rh.add4Neighbour( i->second );
01383     
01384   i = sortedHits.find( east2.rawId() );
01385   if(i != sortedHits.end() ) 
01386     rh.add4Neighbour( i->second );
01387     
01388   i = sortedHits.find( southeast.rawId() );
01389   if(i != sortedHits.end() ) 
01390     rh.add8Neighbour( i->second );
01391     
01392   i = sortedHits.find( southeast2.rawId() );
01393   if(i != sortedHits.end() ) 
01394     rh.add8Neighbour( i->second );
01395     
01396   i = sortedHits.find( west.rawId() );
01397   if(i != sortedHits.end() ) 
01398      rh.add4Neighbour( i->second );
01399    
01400   i = sortedHits.find( west2.rawId() );
01401   if(i != sortedHits.end() ) 
01402      rh.add4Neighbour( i->second );
01403    
01404   i = sortedHits.find( northwest.rawId() );
01405   if(i != sortedHits.end() ) 
01406     rh.add8Neighbour( i->second );
01407 
01408   i = sortedHits.find( northwest2.rawId() );
01409   if(i != sortedHits.end() ) 
01410     rh.add8Neighbour( i->second );
01411 
01412   //  cout<<"----------- rechit print out"<<endl;
01413   // if(( rh.layer() == PFLayer::HF_HAD )||(rh.layer() == PFLayer::HF_EM)) {  
01414     
01415   //   cout<<rh<<endl;
01416     //  }
01417 }
01418 
01419 
01420 
01421 DetId 
01422 PFRecHitProducerHCAL::getSouth(const DetId& id, 
01423                                const CaloSubdetectorTopology& topology) {
01424 
01425   DetId south;
01426   vector<DetId> sids = topology.south(id);
01427   if(sids.size() == 1)
01428     south = sids[0];
01429   
01430   return south;
01431 } 
01432 
01433 
01434 
01435 DetId 
01436 PFRecHitProducerHCAL::getNorth(const DetId& id, 
01437                                const CaloSubdetectorTopology& topology) {
01438 
01439   DetId north;
01440   vector<DetId> nids = topology.north(id);
01441   if(nids.size() == 1)
01442     north = nids[0];
01443   
01444   return north;
01445 } 
01446 
01447
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducerHCAL.cc
