CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoLocalCalo/HcalRecAlgos/src/HcalHFStatusBitFromRecHits.cc

Go to the documentation of this file.
00001 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHFStatusBitFromRecHits.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryData.h" // for eta bounds
00004 
00005 #include <algorithm> // for "max"
00006 #include <cmath>
00007 #include <iostream>
00008 
00009 HcalHFStatusBitFromRecHits::HcalHFStatusBitFromRecHits()
00010 {
00011   // use simple values in default constructor
00012   long_HFlongshortratio_   = 0.99;
00013   short_HFlongshortratio_  = 0.99;
00014   long_thresholdET_    = 0.5; // default energy requirement
00015   short_thresholdET_   = 0.5;
00016   long_thresholdEnergy_ = 100;
00017   short_thresholdEnergy_ = 100;
00018 }
00019 
00020 HcalHFStatusBitFromRecHits::HcalHFStatusBitFromRecHits(double shortR, double shortET, double shortE,
00021                                                        double longR, double longET, double longE)
00022 {
00023   long_HFlongshortratio_   = longR; 
00024   short_HFlongshortratio_  = shortR;
00025   long_thresholdET_        = longET;
00026   short_thresholdET_       = shortET;
00027   long_thresholdEnergy_    = longE;
00028   short_thresholdEnergy_   = shortE;
00029 }
00030 
00031 HcalHFStatusBitFromRecHits::~HcalHFStatusBitFromRecHits(){}
00032 
00033 void HcalHFStatusBitFromRecHits::hfSetFlagFromRecHits(HFRecHitCollection& rec, 
00034                                                       HcalChannelQuality* myqual, 
00035                                                       const HcalSeverityLevelComputer* mySeverity)
00036 {
00037   // Compares energies from long & short fibers
00038   int status;
00039   float en, en2;
00040   int ieta, iphi, depth;
00041   float ratio; // ratio of (L-S)/(L+S) energy magnitudes
00042   double coshEta;
00043 
00044   // Is there a faster way to do this than a double loop?
00045   for (HFRecHitCollection::iterator iHF=rec.begin(); iHF!=rec.end();++iHF)
00046     {
00047       // skip cells that have already been tagged -- shouldn't happen in current algorithm
00048       //if (iHF->flagField(HcalCaloFlagLabels::HFLongShort, HcalCaloFlagLabels::HFLongShort+1)) continue;
00049 
00050       ieta =iHF->id().ieta();  // int between 29-41
00051       // eta = average value between cell eta bounds
00052       coshEta=fabs(cosh(0.5*(theHFEtaBounds[abs(ieta)-29]+theHFEtaBounds[abs(ieta)-28])));
00053 
00054       status=0; // status bit for rechit
00055 
00056       en2=-999; // dummy starting value for partner energy
00057       depth=iHF->id().depth();
00058       en=iHF->energy(); // energy of current rechit
00059 
00060       if (depth==1) // check long fiber
00061         {
00062           if (en<1.2) continue;  // never flag long rechits < 1.2 GeV
00063           if (long_thresholdEnergy_>0. && en<long_thresholdEnergy_) continue;
00064           if (long_thresholdET_>0. && en<long_thresholdET_*coshEta) continue;
00065         }
00066 
00067       else if (depth==2) // check short fiber
00068         {
00069           if (en<1.8) continue;  // never flag short rechits < 1.8 GeV
00070           if (short_thresholdEnergy_>0. && en<short_thresholdEnergy_) continue;
00071           if (short_thresholdET_>0. && en<short_thresholdET_*coshEta) continue;
00072         }
00073       
00074       iphi =iHF->id().iphi();
00075 
00076       // Check for cells whose partners have been excluded from the rechit collections
00077       // Such cells will not get flagged (since we can't make an L vs. S comparison)
00078 
00079       HcalDetId partner(HcalForward, ieta, iphi, 3-depth); // if depth=1, 3-depth =2, and vice versa
00080       DetId detpartner=DetId(partner);
00081       const HcalChannelStatus* partnerstatus=myqual->getValues(detpartner.rawId());
00082       if (mySeverity->dropChannel(partnerstatus->getValue() ) ) continue;  // partner was dropped; don't set flag
00083 
00084       // inner loop will find 'partner' channel (same ieta, iphi, different depth)
00085       for (HFRecHitCollection::iterator iHF2=rec.begin(); iHF2!=rec.end();++iHF2)
00086         {
00087           if (iHF2->id().ieta()!=ieta) continue; // require ieta match
00088           if (iHF2->id().iphi()!=iphi) continue; // require iphi match
00089           if (iHF2->id().depth()==depth) continue;  // require short/long combo
00090 
00091           en2=iHF2->energy(); 
00092 
00093           /* 
00094              We used to use absolute values of energies for ratios, but I don't think we want to do this any more.
00095              For example, for a threshold of 0.995, if en=50 and en2=0, the flag would be set.
00096              But if en=50 and en2<-0.125, the threshold would not be set if using the absolute values.
00097              I don't think we want to have a range of en2 below which the flag is not set.
00098              This does mean that we need to be careful not to set the en energy threshold too low,
00099              so as to not falsely flag fluctuations (en=2, en2=-0.01, for example), but we should never be setting our
00100              thresholds that low.
00101           */
00102           
00103           ratio = (en - en2)/(en + en2);
00104           
00105           if (depth==1 && ratio>long_HFlongshortratio_)
00106             status=1;
00107           else if (depth==2 && ratio>short_HFlongshortratio_)
00108             status=1;
00109           break; // once partner channel found, break out of loop
00110         } // inner loop
00111 
00112       // Consider the case where only one depth present
00113       if (en2==-999) // outer rechit has already passed energy, ET thresholds above; no partner cell means this looks like isolated energy in one channel -- flag it!
00114         status=1;
00115 
00116       // flag rechit as potential PMT window hit
00117       if (status==1) 
00118         iHF->setFlagField(status,HcalCaloFlagLabels::HFLongShort, 1);
00119     } // outer loop
00120   return;
00121 } // void HcalHFStatusBitFromRecHits::hfSetFlagFromRecHits(...)