CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducerHCAL.cc

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