Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:41 

0001 #include "FWCore/Framework/interface/ConsumesCollector.h"
0002 #include "RecoMET/METAlgorithms/interface/HcalHaloAlgo.h"
0003 #include <map>
0004 
0005 /*
0006   [class]:  HcalHaloAlgo
0007   [authors]: R. Remington, The University of Florida
0008   [description]: See HcalHaloAlgo.h
0009   [date]: October 15, 2009
0010 */
0011 namespace {
0012   constexpr float c_cm_per_ns = 29.9792458;
0013   constexpr float zseparation_HBHE = 380.;
0014 };  // namespace
0015 
0016 using namespace reco;
0017 
0018 #include <iomanip>
0019 bool CompareTime(const HBHERecHit* x, const HBHERecHit* y) { return x->time() < y->time(); }
0020 bool CompareTowers(const CaloTower* x, const CaloTower* y) {
0021   return x->iphi() * 1000 + x->ieta() < y->iphi() * 1000 + y->ieta();
0022 }
0023 
0024 HcalHaloAlgo::HcalHaloAlgo(edm::ConsumesCollector iC) : geoToken_(iC.esConsumes()), geo_(nullptr), hgeo_(nullptr) {
0025   HBRecHitEnergyThreshold = 0.;
0026   HERecHitEnergyThreshold = 0.;
0027   SumEnergyThreshold = 0.;
0028   NHitsThreshold = 0;
0029 }
0030 
0031 HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
0032                                      edm::Handle<HBHERecHitCollection>& TheHBHERecHits,
0033                                      edm::Handle<EBRecHitCollection>& TheEBRecHits,
0034                                      edm::Handle<EERecHitCollection>& TheEERecHits,
0035                                      const edm::EventSetup& TheSetup) {
0036   edm::Handle<CaloTowerCollection> TheCaloTowers;
0037   return Calculate(TheCaloGeometry, TheHBHERecHits, TheCaloTowers, TheEBRecHits, TheEERecHits, TheSetup);
0038 }
0039 
0040 HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
0041                                      edm::Handle<HBHERecHitCollection>& TheHBHERecHits,
0042                                      edm::Handle<CaloTowerCollection>& TheCaloTowers,
0043                                      edm::Handle<EBRecHitCollection>& TheEBRecHits,
0044                                      edm::Handle<EERecHitCollection>& TheEERecHits,
0045                                      const edm::EventSetup& TheSetup) {
0046   HcalHaloData TheHcalHaloData;
0047   // ieta overlap geometrically w/ HB
0048   const int iEtaOverlap = 22;
0049   const int nPhiMax = 73;
0050   // Store Energy sum of rechits as a function of iPhi (iPhi goes from 1 to 72)
0051   float SumE[nPhiMax];
0052   // Store Number of rechits as a function of iPhi
0053   int NumHits[nPhiMax];
0054   // Store minimum time of rechit as a function of iPhi
0055   float MinTimeHits[nPhiMax];
0056   // Store maximum time of rechit as a function of iPhi
0057   float MaxTimeHits[nPhiMax];
0058   for (unsigned int i = 0; i < nPhiMax; i++) {
0059     SumE[i] = 0;
0060     NumHits[i] = 0;
0061     MinTimeHits[i] = 0.;
0062     MaxTimeHits[i] = 0.;
0063   }
0064 
0065   for (const auto& hit : (*TheHBHERecHits)) {
0066     HcalDetId id = HcalDetId(hit.id());
0067     switch (id.subdet()) {
0068       case HcalBarrel:
0069         if (hit.energy() < HBRecHitEnergyThreshold)
0070           continue;
0071         break;
0072       case HcalEndcap:
0073         if (hit.energy() < HERecHitEnergyThreshold)
0074           continue;
0075         break;
0076       default:
0077         continue;
0078     }
0079 
0080     int iEta = id.ieta();
0081     int iPhi = id.iphi();
0082     if (iPhi < nPhiMax && std::abs(iEta) <= iEtaOverlap) {
0083       SumE[iPhi] += hit.energy();
0084       NumHits[iPhi]++;
0085 
0086       float time = hit.time();
0087       MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
0088       MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
0089     }
0090   }
0091 
0092   for (int iPhi = 1; iPhi < nPhiMax; iPhi++) {
0093     if (SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold) {
0094       // Build PhiWedge and store to HcalHaloData if energy or #hits pass thresholds
0095       PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
0096 
0097       // Loop over rechits again to calculate direction based on timing info
0098       std::vector<const HBHERecHit*> Hits;
0099       for (const auto& hit : (*TheHBHERecHits)) {
0100         HcalDetId id = HcalDetId(hit.id());
0101         if (id.iphi() != iPhi)
0102           continue;
0103         if (std::abs(id.ieta()) > iEtaOverlap)
0104           continue;  // has to overlap geometrically w/ HB
0105         switch (id.subdet()) {
0106           case HcalBarrel:
0107             if (hit.energy() < HBRecHitEnergyThreshold)
0108               continue;
0109             break;
0110           case HcalEndcap:
0111             if (hit.energy() < HERecHitEnergyThreshold)
0112               continue;
0113             break;
0114           default:
0115             continue;
0116         }
0117         Hits.push_back(&(hit));
0118       }
0119 
0120       std::sort(Hits.begin(), Hits.end(), CompareTime);
0121       float MinusToPlus = 0.;
0122       float PlusToMinus = 0.;
0123       for (unsigned int i = 0; i < Hits.size(); i++) {
0124         HcalDetId id_i = HcalDetId(Hits[i]->id());
0125         int ieta_i = id_i.ieta();
0126         for (unsigned int j = (i + 1); j < Hits.size(); j++) {
0127           HcalDetId id_j = HcalDetId(Hits[j]->id());
0128           int ieta_j = id_j.ieta();
0129           if (ieta_i > ieta_j)
0130             PlusToMinus += std::abs(ieta_i - ieta_j);
0131           else
0132             MinusToPlus += std::abs(ieta_i - ieta_j);
0133         }
0134       }
0135       float PlusZOriginConfidence = (PlusToMinus + MinusToPlus) ? PlusToMinus / (PlusToMinus + MinusToPlus) : -1.;
0136       wedge.SetPlusZOriginConfidence(PlusZOriginConfidence);
0137       TheHcalHaloData.GetPhiWedges().push_back(wedge);
0138     }
0139   }
0140 
0141   // Don't use HF.
0142   int maxAbsIEta = 29;
0143 
0144   std::map<int, float> iPhiHadEtMap;
0145   std::vector<const CaloTower*> sortedCaloTowers;
0146   for (const auto& tower : (*TheCaloTowers)) {
0147     if (std::abs(tower.ieta()) > maxAbsIEta)
0148       continue;
0149 
0150     int iPhi = tower.iphi();
0151     if (!iPhiHadEtMap.count(iPhi))
0152       iPhiHadEtMap[iPhi] = 0.0;
0153     iPhiHadEtMap[iPhi] += tower.hadEt();
0154 
0155     if (tower.numProblematicHcalCells() > 0)
0156       sortedCaloTowers.push_back(&(tower));
0157   }
0158 
0159   // Sort towers such that lowest iphi and ieta are first, highest last, and towers
0160   // with same iphi value are consecutive. Then we can do everything else in one loop.
0161   std::sort(sortedCaloTowers.begin(), sortedCaloTowers.end(), CompareTowers);
0162 
0163   HaloTowerStrip strip;
0164 
0165   int prevIEta = -99, prevIPhi = -99;
0166   float prevHadEt = 0.;
0167   float prevEmEt = 0.;
0168   std::pair<uint8_t, CaloTowerDetId> prevPair, towerPair;
0169   bool wasContiguous = true;
0170 
0171   // Loop through and store a vector of pairs (problematicCells, DetId) for each contiguous strip we find
0172   for (unsigned int i = 0; i < sortedCaloTowers.size(); i++) {
0173     const CaloTower* tower = sortedCaloTowers[i];
0174 
0175     towerPair = std::make_pair((uint8_t)tower->numProblematicHcalCells(), tower->id());
0176 
0177     bool newIPhi = tower->iphi() != prevIPhi;
0178     bool isContiguous = tower->ieta() == 1 ? tower->ieta() - 2 == prevIEta : tower->ieta() - 1 == prevIEta;
0179 
0180     isContiguous = isContiguous || (tower->ieta() == -maxAbsIEta);
0181     if (newIPhi)
0182       isContiguous = false;
0183 
0184     if (!wasContiguous && isContiguous) {
0185       strip.cellTowerIds.push_back(prevPair);
0186       strip.cellTowerIds.push_back(towerPair);
0187       strip.hadEt += prevHadEt + tower->hadEt();
0188       strip.emEt += prevEmEt + tower->emEt();
0189     }
0190 
0191     if (wasContiguous && isContiguous) {
0192       strip.cellTowerIds.push_back(towerPair);
0193       strip.hadEt += tower->hadEt();
0194       strip.emEt += tower->emEt();
0195     }
0196 
0197     if ((wasContiguous && !isContiguous) || i == sortedCaloTowers.size() - 1) {  //ended the strip, so flush it
0198 
0199       if (strip.cellTowerIds.size() > 3) {
0200         int iPhi = strip.cellTowerIds.at(0).second.iphi();
0201         int iPhiLower = (iPhi == 1) ? 72 : iPhi - 1;
0202         int iPhiUpper = (iPhi == 72) ? 1 : iPhi + 1;
0203 
0204         float energyRatio = 0.0;
0205         if (iPhiHadEtMap.count(iPhiLower))
0206           energyRatio += iPhiHadEtMap[iPhiLower];
0207         if (iPhiHadEtMap.count(iPhiUpper))
0208           energyRatio += iPhiHadEtMap[iPhiUpper];
0209         iPhiHadEtMap[iPhi] = std::max(iPhiHadEtMap[iPhi], 0.001F);
0210 
0211         energyRatio /= iPhiHadEtMap[iPhi];
0212         strip.energyRatio = energyRatio;
0213 
0214         TheHcalHaloData.getProblematicStrips().push_back(strip);
0215       }
0216       strip = HaloTowerStrip();
0217     }
0218 
0219     wasContiguous = isContiguous;
0220     prevPair = towerPair;
0221     prevEmEt = tower->emEt();
0222     prevIPhi = tower->iphi();
0223     prevIEta = tower->ieta();
0224     prevHadEt = tower->hadEt();
0225   }
0226 
0227   geo_ = &TheSetup.getData(geoToken_);
0228   hgeo_ = dynamic_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, 1));
0229 
0230   //Halo cluster building:
0231   //Various clusters are built, depending on the subdetector.
0232   //In barrel, one looks for deposits narrow in phi.
0233   //In endcaps, one looks for localized deposits (dr condition in EE where r =sqrt(dphi*dphi+deta*deta)
0234   //E/H condition is also applied.
0235   //The halo cluster building step targets a large efficiency (ideally >99%) for beam halo deposits.
0236   //These clusters are used as input for the halo pattern finding methods in HcalHaloAlgo and for the CSC-calo matching methods in GlobalHaloAlgo.
0237 
0238   //Et threshold hardcoded for now. Might one to get it from config
0239 
0240   std::vector<HaloClusterCandidateHCAL> haloclustercands_HB;
0241   haloclustercands_HB = GetHaloClusterCandidateHB(TheEBRecHits, TheHBHERecHits, 5);
0242 
0243   std::vector<HaloClusterCandidateHCAL> haloclustercands_HE;
0244   haloclustercands_HE = GetHaloClusterCandidateHE(TheEERecHits, TheHBHERecHits, 10);
0245 
0246   TheHcalHaloData.setHaloClusterCandidatesHB(haloclustercands_HB);
0247   TheHcalHaloData.setHaloClusterCandidatesHE(haloclustercands_HE);
0248 
0249   return TheHcalHaloData;
0250 }
0251 
0252 std::vector<HaloClusterCandidateHCAL> HcalHaloAlgo::GetHaloClusterCandidateHB(
0253     edm::Handle<EcalRecHitCollection>& ecalrechitcoll,
0254     edm::Handle<HBHERecHitCollection>& hbherechitcoll,
0255     float et_thresh_seedrh) {
0256   std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHB;
0257 
0258   reco::Vertex::Point vtx(0, 0, 0);
0259 
0260   for (size_t ihit = 0; ihit < hbherechitcoll->size(); ++ihit) {
0261     HaloClusterCandidateHCAL clustercand;
0262 
0263     const HBHERecHit& rechit = (*hbherechitcoll)[ihit];
0264     math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
0265     //Et condition
0266     double rhet = rechit.energy() * sqrt(rhpos.perp2() / rhpos.mag2());
0267     if (rhet < et_thresh_seedrh)
0268       continue;
0269     if (std::abs(rhpos.z()) > zseparation_HBHE)
0270       continue;
0271     double eta = rhpos.eta();
0272     double phi = rhpos.phi();
0273 
0274     bool isiso = true;
0275     double etcluster(0);
0276     int nbtowerssameeta(0);
0277     double timediscriminatorITBH(0), timediscriminatorOTBH(0);
0278     double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
0279 
0280     //Building the cluster
0281     edm::RefVector<HBHERecHitCollection> bhrhcandidates;
0282     for (size_t jhit = 0; jhit < hbherechitcoll->size(); ++jhit) {
0283       const HBHERecHit& rechitj = (*hbherechitcoll)[jhit];
0284       HBHERecHitRef rhRef(hbherechitcoll, jhit);
0285       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0286       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0287       if (rhetj < 2)
0288         continue;
0289       if (std::abs(rhposj.z()) > zseparation_HBHE)
0290         continue;
0291       double etaj = rhposj.eta();
0292       double phij = rhposj.phi();
0293       double deta = eta - etaj;
0294       double dphi = deltaPhi(phi, phij);
0295       if (std::abs(deta) > 0.4)
0296         continue;  //This means +/-4 towers in eta
0297       if (std::abs(dphi) > 0.2)
0298         continue;  //This means +/-2 towers in phi
0299       if (std::abs(dphi) > 0.1 && std::abs(deta) < 0.2) {
0300         isiso = false;
0301         break;
0302       }  //The strip should be isolated
0303       if (std::abs(dphi) > 0.1)
0304         continue;
0305       if (std::abs(dphi) < 0.05)
0306         nbtowerssameeta++;
0307       if (dphi > 0.05)
0308         etstrip_phiseedplus1 += rhetj;
0309       if (dphi < -0.05)
0310         etstrip_phiseedminus1 += rhetj;
0311 
0312       etcluster += rhetj;
0313       //Timing discriminator
0314       //We assign a weight to the rechit defined as:
0315       //Log10(Et)*f(T,R,Z)
0316       //where f(T,R,Z) is the separation curve between halo-like and IP-like times.
0317       //The time difference between a deposit from a outgoing IT halo and a deposit coming from a particle emitted at the IP is given by:
0318       //dt= ( - sqrt(R^2+z^2) + |z| )/c
0319       // For OT beam halo, the time difference is:
0320       //dt= ( 25 + sqrt(R^2+z^2) + |z| )/c
0321       //only consider the central part of HB as things get hard at large z.
0322       //The best fitted value for R leads to 240 cm (IT) and 330 cm (OT)
0323       double rhtj = rechitj.time();
0324       timediscriminatorITBH +=
0325           std::log10(rhetj) *
0326           (rhtj + 0.5 * (sqrt(240. * 240. + rhposj.z() * rhposj.z()) - std::abs(rhposj.z())) / c_cm_per_ns);
0327       if (std::abs(rhposj.z()) < 300)
0328         timediscriminatorOTBH +=
0329             std::log10(rhetj) *
0330             (rhtj - 0.5 * (25 - (sqrt(330. * 330. + rhposj.z() * rhposj.z()) + std::abs(rhposj.z())) / c_cm_per_ns));
0331       bhrhcandidates.push_back(rhRef);
0332     }
0333     //Isolation conditions
0334     if (!isiso)
0335       continue;
0336     if (etstrip_phiseedplus1 / etcluster > 0.2 && etstrip_phiseedminus1 / etcluster > 0.2)
0337       continue;
0338 
0339     //Calculate E/H
0340     double eoh(0);
0341     for (size_t jhit = 0; jhit < ecalrechitcoll->size(); ++jhit) {
0342       const EcalRecHit& rechitj = (*ecalrechitcoll)[jhit];
0343       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0344       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0345       if (rhetj < 2)
0346         continue;
0347       double etaj = rhposj.eta();
0348       double phij = rhposj.phi();
0349       if (std::abs(eta - etaj) > 0.2)
0350         continue;
0351       if (std::abs(deltaPhi(phi, phij)) > 0.2)
0352         continue;
0353       eoh += rhetj / etcluster;
0354     }
0355     //E/H condition
0356     if (eoh > 0.1)
0357       continue;
0358 
0359     clustercand.setClusterEt(etcluster);
0360     clustercand.setSeedEt(rhet);
0361     clustercand.setSeedEta(eta);
0362     clustercand.setSeedPhi(phi);
0363     clustercand.setSeedZ(rhpos.Z());
0364     clustercand.setSeedR(sqrt(rhpos.perp2()));
0365     clustercand.setSeedTime(rechit.time());
0366     clustercand.setEoverH(eoh);
0367     clustercand.setNbTowersInEta(nbtowerssameeta);
0368     clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
0369     clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
0370     clustercand.setTimeDiscriminatorITBH(timediscriminatorITBH);
0371     clustercand.setTimeDiscriminatorOTBH(timediscriminatorOTBH);
0372     clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
0373 
0374     bool isbeamhalofrompattern = HBClusterShapeandTimeStudy(clustercand, false);
0375     clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
0376     bool isbeamhalofrompattern_hlt = HBClusterShapeandTimeStudy(clustercand, true);
0377     clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
0378 
0379     TheHaloClusterCandsHB.push_back(clustercand);
0380   }
0381 
0382   return TheHaloClusterCandsHB;
0383 }
0384 
0385 std::vector<HaloClusterCandidateHCAL> HcalHaloAlgo::GetHaloClusterCandidateHE(
0386     edm::Handle<EcalRecHitCollection>& ecalrechitcoll,
0387     edm::Handle<HBHERecHitCollection>& hbherechitcoll,
0388     float et_thresh_seedrh) {
0389   std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHE;
0390 
0391   reco::Vertex::Point vtx(0, 0, 0);
0392 
0393   for (size_t ihit = 0; ihit < hbherechitcoll->size(); ++ihit) {
0394     HaloClusterCandidateHCAL clustercand;
0395 
0396     const HBHERecHit& rechit = (*hbherechitcoll)[ihit];
0397     math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
0398     //Et condition
0399     double rhet = rechit.energy() * sqrt(rhpos.perp2() / rhpos.mag2());
0400     if (rhet < et_thresh_seedrh)
0401       continue;
0402     if (std::abs(rhpos.z()) < zseparation_HBHE)
0403       continue;
0404     double eta = rhpos.eta();
0405     double phi = rhpos.phi();
0406     double rhr = sqrt(rhpos.perp2());
0407     bool isiso = true;
0408     double etcluster(0), hdepth1(0);
0409     int clustersize(0);
0410     double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
0411 
0412     //Building the cluster
0413     edm::RefVector<HBHERecHitCollection> bhrhcandidates;
0414     for (size_t jhit = 0; jhit < hbherechitcoll->size(); ++jhit) {
0415       const HBHERecHit& rechitj = (*hbherechitcoll)[jhit];
0416       HBHERecHitRef rhRef(hbherechitcoll, jhit);
0417       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0418       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0419       if (rhetj < 2)
0420         continue;
0421       if (std::abs(rhposj.z()) < zseparation_HBHE)
0422         continue;
0423       if (rhpos.z() * rhposj.z() < 0)
0424         continue;
0425       double phij = rhposj.phi();
0426       double dphi = deltaPhi(phi, phij);
0427       if (std::abs(dphi) > 0.4)
0428         continue;
0429       double rhrj = sqrt(rhposj.perp2());
0430       if (std::abs(rhr - rhrj) > 50)
0431         continue;
0432       if (std::abs(dphi) > 0.2 || std::abs(rhr - rhrj) > 20) {
0433         isiso = false;
0434         break;
0435       }  //The deposit should be isolated
0436       if (dphi > 0.05)
0437         etstrip_phiseedplus1 += rhetj;
0438       if (dphi < -0.05)
0439         etstrip_phiseedminus1 += rhetj;
0440       clustersize++;
0441       etcluster += rhetj;
0442       if (std::abs(rhposj.z()) < 405)
0443         hdepth1 += rhetj;
0444       //No timing condition for now in HE
0445       bhrhcandidates.push_back(rhRef);
0446     }
0447     //Isolation conditions
0448     if (!isiso)
0449       continue;
0450     if (etstrip_phiseedplus1 / etcluster > 0.1 && etstrip_phiseedminus1 / etcluster > 0.1)
0451       continue;
0452 
0453     //Calculate E/H
0454     double eoh(0);
0455     for (size_t jhit = 0; jhit < ecalrechitcoll->size(); ++jhit) {
0456       const EcalRecHit& rechitj = (*ecalrechitcoll)[jhit];
0457       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0458       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0459       if (rhetj < 2)
0460         continue;
0461       if (rhpos.z() * rhposj.z() < 0)
0462         continue;
0463       double etaj = rhposj.eta();
0464       double phij = rhposj.phi();
0465       double dr = sqrt((eta - etaj) * (eta - etaj) + deltaPhi(phi, phij) * deltaPhi(phi, phij));
0466       if (dr > 0.3)
0467         continue;
0468 
0469       eoh += rhetj / etcluster;
0470     }
0471     //E/H condition
0472     if (eoh > 0.1)
0473       continue;
0474 
0475     clustercand.setClusterEt(etcluster);
0476     clustercand.setSeedEt(rhet);
0477     clustercand.setSeedEta(eta);
0478     clustercand.setSeedPhi(phi);
0479     clustercand.setSeedZ(rhpos.Z());
0480     clustercand.setSeedR(sqrt(rhpos.perp2()));
0481     clustercand.setSeedTime(rechit.time());
0482     clustercand.setEoverH(eoh);
0483     clustercand.setH1overH123(hdepth1 / etcluster);
0484     clustercand.setClusterSize(clustersize);
0485     clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
0486     clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
0487     clustercand.setTimeDiscriminator(0);
0488     clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
0489 
0490     bool isbeamhalofrompattern = HEClusterShapeandTimeStudy(clustercand, false);
0491     clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
0492     bool isbeamhalofrompattern_hlt = HEClusterShapeandTimeStudy(clustercand, true);
0493     clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
0494 
0495     TheHaloClusterCandsHE.push_back(clustercand);
0496   }
0497 
0498   return TheHaloClusterCandsHE;
0499 }
0500 
0501 bool HcalHaloAlgo::HBClusterShapeandTimeStudy(HaloClusterCandidateHCAL hcand, bool ishlt) {
0502   //Conditions on the central strip size in eta.
0503   //For low size, extra conditions on seed et, isolation and cluster timing
0504   //Here we target both IT and OT beam halo. Two separate discriminators were built for the two cases.
0505 
0506   if (hcand.getSeedEt() < 10)
0507     return false;
0508 
0509   if (hcand.getNbTowersInEta() < 3)
0510     return false;
0511   //Isolation criteria for very short eta strips
0512   if (hcand.getNbTowersInEta() == 3 && (hcand.getEtStripPhiSeedPlus1() > 0.1 || hcand.getEtStripPhiSeedMinus1() > 0.1))
0513     return false;
0514   if (hcand.getNbTowersInEta() <= 5 && (hcand.getEtStripPhiSeedPlus1() > 0.1 && hcand.getEtStripPhiSeedMinus1() > 0.1))
0515     return false;
0516 
0517   //Timing conditions for short eta strips
0518   if (hcand.getNbTowersInEta() == 3 && hcand.getTimeDiscriminatorITBH() >= 0.)
0519     return false;
0520   if (hcand.getNbTowersInEta() <= 6 && hcand.getTimeDiscriminatorITBH() >= 5. && hcand.getTimeDiscriminatorOTBH() < 0.)
0521     return false;
0522 
0523   //For HLT, only use conditions without timing
0524   if (ishlt && hcand.getNbTowersInEta() < 7)
0525     return false;
0526 
0527   hcand.setIsHaloFromPattern(true);
0528 
0529   return true;
0530 }
0531 
0532 bool HcalHaloAlgo::HEClusterShapeandTimeStudy(HaloClusterCandidateHCAL hcand, bool ishlt) {
0533   //Conditions on H1/H123 to spot halo interacting only in one HCAL layer.
0534   //For R> about 170cm, HE has only one layer and this condition cannot be applied
0535   //Note that for R>170 cm, the halo is in CSC acceptance and will most likely be spotted by the CSC-calo matching method
0536   //A method to identify halos interacting in both H1 and H2/H3 at low R is still missing.
0537 
0538   if (hcand.getSeedEt() < 20)
0539     return false;
0540   if (hcand.getSeedR() > 170)
0541     return false;
0542 
0543   if (hcand.getH1overH123() > 0.02 && hcand.getH1overH123() < 0.98)
0544     return false;
0545 
0546   //This method is one of the ones with the highest fake rate: in JetHT dataset, it happens in around 0.1% of the cases that a low pt jet (pt= 20) leaves all of its energy in only one HCAL layer.
0547   //At HLT, one only cares about large deposits from BH that would lead to a MET/SinglePhoton trigger to be fired.
0548   //Rising the seed Et threshold at HLT has therefore little impact on the HLT performances but ensures that possible controversial events are still recorded.
0549   if (ishlt && hcand.getSeedEt() < 50)
0550     return false;
0551 
0552   hcand.setIsHaloFromPattern(true);
0553 
0554   return true;
0555 }
0556 
0557 math::XYZPoint HcalHaloAlgo::getPosition(const DetId& id, reco::Vertex::Point vtx) {
0558   const GlobalPoint pos = ((id.det() == DetId::Hcal) ? hgeo_->getPosition(id) : GlobalPoint(geo_->getPosition(id)));
0559   math::XYZPoint posV(pos.x() - vtx.x(), pos.y() - vtx.y(), pos.z() - vtx.z());
0560   return posV;
0561 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!