METAlgorithms/src/EcalHaloAlgo.cc

0001 #include "FWCore/Framework/interface/ConsumesCollector.h"
0002 #include "RecoMET/METAlgorithms/interface/EcalHaloAlgo.h"
0003 #include "DataFormats/Common/interface/ValueMap.h"
0004
0005 /*
0006   [class]:  EcalHaloAlgo
0007   [authors]: R. Remington, The University of Florida
0008   [description]: See EcalHaloAlgo.h
0009   [date]: October 15, 2009
0010 */
0011 namespace {
0012   constexpr float c_cm_per_ns = 29.9792458;
0013 };
0014 using namespace std;
0015 using namespace reco;
0016 using namespace edm;
0017
0018 bool CompareTime(const EcalRecHit* x, const EcalRecHit* y) { return x->time() < y->time(); }
0019
0020 EcalHaloAlgo::EcalHaloAlgo(edm::ConsumesCollector iC) : geoToken_(iC.esConsumes()) {
0021   RoundnessCut = 0;
0022   AngleCut = 0;
0023   EBRecHitEnergyThreshold = 0.;
0024   EERecHitEnergyThreshold = 0.;
0025   ESRecHitEnergyThreshold = 0.;
0026   SumEnergyThreshold = 0.;
0027   NHitsThreshold = 0;
0028
0029   geo = nullptr;
0030 }
0031
0032 EcalHaloData EcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
0033                                      edm::Handle<reco::PhotonCollection>& ThePhotons,
0034                                      edm::Handle<reco::SuperClusterCollection>& TheSuperClusters,
0035                                      edm::Handle<EBRecHitCollection>& TheEBRecHits,
0036                                      edm::Handle<EERecHitCollection>& TheEERecHits,
0037                                      edm::Handle<ESRecHitCollection>& TheESRecHits,
0038                                      edm::Handle<HBHERecHitCollection>& TheHBHERecHits,
0039                                      const edm::EventSetup& TheSetup) {
0040   EcalHaloData TheEcalHaloData;
0041
0042   // Store energy sum of rechits as a function of iPhi (iphi goes from 1 to 72)
0043   float SumE[361];
0044   // Store number of rechits as a function of iPhi
0045   int NumHits[361];
0046   // Store minimum time of rechit as a function of iPhi
0047   float MinTimeHits[361];
0048   // Store maximum time of rechit as a function of iPhi
0049   float MaxTimeHits[361];
0050
0051   // initialize
0052   for (int i = 0; i < 361; i++) {
0053     SumE[i] = 0.;
0054     NumHits[i] = 0;
0055     MinTimeHits[i] = 9999.;
0056     MaxTimeHits[i] = -9999.;
0057   }
0058
0059   // Loop over EB RecHits
0060   for (EBRecHitCollection::const_iterator hit = TheEBRecHits->begin(); hit != TheEBRecHits->end(); hit++) {
0061     // Arbitrary threshold to kill noise (needs to be optimized with data)
0062     if (hit->energy() < EBRecHitEnergyThreshold)
0063       continue;
0064
0065     // Get Det Id of the rechit
0066     DetId id = DetId(hit->id());
0067
0068     // Get EB geometry
0069     const CaloSubdetectorGeometry* TheSubGeometry = TheCaloGeometry.getSubdetectorGeometry(DetId::Ecal, 1);
0070     EBDetId EcalID(id.rawId());
0071     auto cell = (TheSubGeometry) ? (TheSubGeometry->getGeometry(id)) : decltype(TheSubGeometry->getGeometry(id))();
0072
0073     if (cell) {
0074       // GlobalPoint globalpos = cell->getPosition();
0075       //      float r = TMath::Sqrt ( globalpos.y()*globalpos.y() + globalpos.x()*globalpos.x());
0076       int iPhi = EcalID.iphi();
0077
0078       if (iPhi < 361)  // just to be safe
0079       {
0080         //iPhi = (iPhi-1)/5 +1;  // convert ecal iphi to phiwedge iphi  (e.g. there are 5 crystal per phi wedge, as in calotowers )
0081         SumE[iPhi] += hit->energy();
0082         NumHits[iPhi]++;
0083
0084         float time = hit->time();
0085         MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
0086         MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
0087       }
0088     }
0089   }
0090
0091   //for( int iPhi = 1 ; iPhi < 73; iPhi++ )
0092   for (int iPhi = 1; iPhi < 361; iPhi++) {
0093     if (SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold) {
0094       // Build PhiWedge and store to EcalHaloData 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
0099       // Loop over EB RecHits
0100       std::vector<const EcalRecHit*> Hits;
0101       for (EBRecHitCollection::const_iterator hit = TheEBRecHits->begin(); hit != TheEBRecHits->end(); hit++) {
0102         if (hit->energy() < EBRecHitEnergyThreshold)
0103           continue;
0104
0105         // Get Det Id of the rechit
0106         DetId id = DetId(hit->id());
0107         EBDetId EcalID(id.rawId());
0108         int Hit_iPhi = EcalID.iphi();
0109         //Hit_iPhi = (Hit_iPhi-1)/5 +1; // convert ecal iphi to phiwedge iphi
0110         if (Hit_iPhi != iPhi)
0111           continue;
0112         Hits.push_back(&(*hit));
0113       }
0114       std::sort(Hits.begin(), Hits.end(), CompareTime);
0115       float MinusToPlus = 0.;
0116       float PlusToMinus = 0.;
0117       for (unsigned int i = 0; i < Hits.size(); i++) {
0118         DetId id_i = DetId(Hits[i]->id());
0119         EBDetId EcalID_i(id_i.rawId());
0120         int ieta_i = EcalID_i.ieta();
0121         for (unsigned int j = (i + 1); j < Hits.size(); j++) {
0122           DetId id_j = DetId(Hits[j]->id());
0123           EBDetId EcalID_j(id_j.rawId());
0124           int ieta_j = EcalID_j.ieta();
0125           if (ieta_i > ieta_j)
0126             PlusToMinus += TMath::Abs(ieta_j - ieta_i);
0127           else
0128             MinusToPlus += TMath::Abs(ieta_j - ieta_i);
0129         }
0130       }
0131
0132       float PlusZOriginConfidence = (PlusToMinus + MinusToPlus) ? PlusToMinus / (PlusToMinus + MinusToPlus) : -1.;
0133       wedge.SetPlusZOriginConfidence(PlusZOriginConfidence);
0134       TheEcalHaloData.GetPhiWedges().push_back(wedge);
0135     }
0136   }
0137
0138   std::vector<float> vShowerShapes_Roundness;
0139   std::vector<float> vShowerShapes_Angle;
0140   if (TheSuperClusters.isValid()) {
0141     for (reco::SuperClusterCollection::const_iterator cluster = TheSuperClusters->begin();
0142          cluster != TheSuperClusters->end();
0143          cluster++) {
0144       if (abs(cluster->eta()) <= 1.48) {
0145         vector<float> shapes = EcalClusterTools::roundnessBarrelSuperClusters(*cluster, (*TheEBRecHits.product()));
0146         float roundness = shapes[0];
0147         float angle = shapes[1];
0148
0149         // Check if supercluster belongs to photon and passes the cuts on roundness and angle, if so store the reference to it
0150         if ((roundness >= 0 && roundness < GetRoundnessCut()) && angle >= 0 && angle < GetAngleCut()) {
0151           edm::Ref<SuperClusterCollection> TheClusterRef(TheSuperClusters, cluster - TheSuperClusters->begin());
0152           bool BelongsToPhoton = false;
0153           if (ThePhotons.isValid()) {
0154             for (reco::PhotonCollection::const_iterator iPhoton = ThePhotons->begin(); iPhoton != ThePhotons->end();
0155                  iPhoton++) {
0156               if (iPhoton->isEB())
0157                 if (TheClusterRef == iPhoton->superCluster()) {
0158                   BelongsToPhoton = true;
0159                   break;
0160                 }
0161             }
0162           }
0163           //Only store refs to suspicious EB SuperClusters which belong to Photons
0164           //Showershape variables are more discriminating for these cases
0165           if (BelongsToPhoton) {
0166             TheEcalHaloData.GetSuperClusters().push_back(TheClusterRef);
0167           }
0168         }
0169         vShowerShapes_Roundness.push_back(shapes[0]);
0170         vShowerShapes_Angle.push_back(shapes[1]);
0171       } else {
0172         vShowerShapes_Roundness.push_back(-1.);
0173         vShowerShapes_Angle.push_back(-1.);
0174       }
0175     }
0176
0177     edm::ValueMap<float>::Filler TheRoundnessFiller(TheEcalHaloData.GetShowerShapesRoundness());
0178     TheRoundnessFiller.insert(TheSuperClusters, vShowerShapes_Roundness.begin(), vShowerShapes_Roundness.end());
0179     TheRoundnessFiller.fill();
0180
0181     edm::ValueMap<float>::Filler TheAngleFiller(TheEcalHaloData.GetShowerShapesAngle());
0182     TheAngleFiller.insert(TheSuperClusters, vShowerShapes_Angle.begin(), vShowerShapes_Angle.end());
0183     TheAngleFiller.fill();
0184   }
0185
0186   geo = &TheSetup.getData(geoToken_);
0187
0188   //Halo cluster building:
0189   //Various clusters are built, depending on the subdetector.
0190   //In barrel, one looks for deposits narrow in phi.
0191   //In endcaps, one looks for localized deposits (dr condition in EE where r =sqrt(dphi*dphi+deta*deta)
0192   //H/E condition is also applied in EB.
0193   //The halo cluster building step targets a large efficiency (ideally >99%) for beam halo deposits.
0194   //These clusters are used as input for the halo pattern finding methods in EcalHaloAlgo and for the CSC-calo matching methods in GlobalHaloAlgo.
0195
0196   //Et threshold hardcoded for now. Might one to get it from config
0197   std::vector<HaloClusterCandidateECAL> haloclustercands_EB;
0198   haloclustercands_EB = GetHaloClusterCandidateEB(TheEBRecHits, TheHBHERecHits, 5);
0199
0200   std::vector<HaloClusterCandidateECAL> haloclustercands_EE;
0201   haloclustercands_EE = GetHaloClusterCandidateEE(TheEERecHits, TheHBHERecHits, 10);
0202
0203   TheEcalHaloData.setHaloClusterCandidatesEB(haloclustercands_EB);
0204   TheEcalHaloData.setHaloClusterCandidatesEE(haloclustercands_EE);
0205
0206   return TheEcalHaloData;
0207 }
0208
0209 std::vector<HaloClusterCandidateECAL> EcalHaloAlgo::GetHaloClusterCandidateEB(
0210     edm::Handle<EcalRecHitCollection>& ecalrechitcoll,
0211     edm::Handle<HBHERecHitCollection>& hbherechitcoll,
0212     float et_thresh_seedrh) {
0213   std::vector<HaloClusterCandidateECAL> TheHaloClusterCandsEB;
0214   reco::Vertex::Point vtx(0, 0, 0);
0215
0216   for (size_t ihit = 0; ihit < ecalrechitcoll->size(); ++ihit) {
0217     HaloClusterCandidateECAL clustercand;
0218
0219     const EcalRecHit& rechit = (*ecalrechitcoll)[ihit];
0220     math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
0221     //Et condition
0222
0223     double rhet = rechit.energy() * sqrt(rhpos.perp2() / rhpos.mag2());
0224     if (rhet < et_thresh_seedrh)
0225       continue;
0226     double eta = rhpos.eta();
0227     double phi = rhpos.phi();
0228
0229     bool isiso = true;
0230     double etcluster(0);
0231     int nbcrystalsameeta(0);
0232     double timediscriminator(0);
0233     double etstrip_iphiseedplus1(0), etstrip_iphiseedminus1(0);
0234
0235     //Building the cluster
0236     edm::RefVector<EcalRecHitCollection> bhrhcandidates;
0237     for (size_t jhit = 0; jhit < ecalrechitcoll->size(); ++jhit) {
0238       const EcalRecHit& rechitj = (*ecalrechitcoll)[jhit];
0239       EcalRecHitRef rhRef(ecalrechitcoll, jhit);
0240       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0241
0242       double etaj = rhposj.eta();
0243       double phij = rhposj.phi();
0244
0245       double deta = eta - etaj;
0246       double dphi = deltaPhi(phi, phij);
0247       if (std::abs(deta) > 0.2)
0248         continue;  //This means +/-11 crystals in eta
0249       if (std::abs(dphi) > 0.08)
0250         continue;  //This means +/-4 crystals in phi
0251
0252       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0253       //Rechits with et between 1 and 2 GeV are saved in the rh list but not used in the calculation of the halocluster variables
0254       if (rhetj < 1)
0255         continue;
0256       bhrhcandidates.push_back(rhRef);
0257       if (rhetj < 2)
0258         continue;
0259
0260       if (std::abs(dphi) > 0.03) {
0261         isiso = false;
0262         break;
0263       }  //The strip should be isolated
0264       if (std::abs(dphi) < 0.01)
0265         nbcrystalsameeta++;
0266       if (dphi > 0.01)
0267         etstrip_iphiseedplus1 += rhetj;
0268       if (dphi < -0.01)
0269         etstrip_iphiseedminus1 += rhetj;
0270       etcluster += rhetj;
0271       //Timing discriminator
0272       //We assign a weight to the rechit defined as:
0273       //Log10(Et)*f(T,R,Z)
0274       //where f(T,R,Z) is the separation curve between halo-like and IP-like times.
0275       //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:
0276       //dt= ( - sqrt(R^2+z^2) + |z| )/c
0277       //Here we take R to be 130 cm.
0278       //For EB, the function was parametrized as a function of ieta instead of Z.
0279       double rhtj = rechitj.time();
0280       EBDetId detj = rechitj.id();
0281       int rhietaj = detj.ieta();
0282       timediscriminator += std::log10(rhetj) *
0283                            (rhtj + 0.5 * (sqrt(16900 + 9 * rhietaj * rhietaj) - 3 * std::abs(rhietaj)) / c_cm_per_ns);
0284     }
0285     //Isolation condition
0286     if (!isiso)
0287       continue;
0288
0289     //Calculate H/E
0290     double hoe(0);
0291     for (size_t jhit = 0; jhit < hbherechitcoll->size(); ++jhit) {
0292       const HBHERecHit& rechitj = (*hbherechitcoll)[jhit];
0293       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0294       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0295       if (rhetj < 2)
0296         continue;
0297       double etaj = rhposj.eta();
0298       double phij = rhposj.phi();
0299       double deta = eta - etaj;
0300       double dphi = deltaPhi(phi, phij);
0301       if (std::abs(deta) > 0.2)
0302         continue;
0303       if (std::abs(dphi) > 0.2)
0304         continue;
0305       hoe += rhetj / etcluster;
0306     }
0307     //H/E condition
0308     if (hoe > 0.1)
0309       continue;
0310
0311     clustercand.setClusterEt(etcluster);
0312     clustercand.setSeedEt(rhet);
0313     clustercand.setSeedEta(eta);
0314     clustercand.setSeedPhi(phi);
0315     clustercand.setSeedZ(rhpos.Z());
0316     clustercand.setSeedR(sqrt(rhpos.perp2()));
0317     clustercand.setSeedTime(rechit.time());
0318     clustercand.setHoverE(hoe);
0319     clustercand.setNbofCrystalsInEta(nbcrystalsameeta);
0320     clustercand.setEtStripIPhiSeedPlus1(etstrip_iphiseedplus1);
0321     clustercand.setEtStripIPhiSeedMinus1(etstrip_iphiseedminus1);
0322     clustercand.setTimeDiscriminator(timediscriminator);
0323     clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
0324
0325     bool isbeamhalofrompattern = EBClusterShapeandTimeStudy(clustercand, false);
0326     clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
0327
0328     bool isbeamhalofrompattern_hlt = EBClusterShapeandTimeStudy(clustercand, true);
0329     clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
0330
0331     TheHaloClusterCandsEB.push_back(clustercand);
0332   }
0333
0334   return TheHaloClusterCandsEB;
0335 }
0336
0337 std::vector<HaloClusterCandidateECAL> EcalHaloAlgo::GetHaloClusterCandidateEE(
0338     edm::Handle<EcalRecHitCollection>& ecalrechitcoll,
0339     edm::Handle<HBHERecHitCollection>& hbherechitcoll,
0340     float et_thresh_seedrh) {
0341   std::vector<HaloClusterCandidateECAL> TheHaloClusterCandsEE;
0342
0343   reco::Vertex::Point vtx(0, 0, 0);
0344
0345   for (size_t ihit = 0; ihit < ecalrechitcoll->size(); ++ihit) {
0346     HaloClusterCandidateECAL clustercand;
0347
0348     const EcalRecHit& rechit = (*ecalrechitcoll)[ihit];
0349     math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
0350     //Et condition
0351     double rhet = rechit.energy() * sqrt(rhpos.perp2() / rhpos.mag2());
0352     if (rhet < et_thresh_seedrh)
0353       continue;
0354     double eta = rhpos.eta();
0355     double phi = rhpos.phi();
0356     double rhr = sqrt(rhpos.perp2());
0357
0358     bool isiso = true;
0359     double etcluster(0);
0360     double timediscriminator(0);
0361     int clustersize(0);
0362     int nbcrystalssmallt(0);
0363     int nbcrystalshight(0);
0364     //Building the cluster
0365     edm::RefVector<EcalRecHitCollection> bhrhcandidates;
0366     for (size_t jhit = 0; jhit < ecalrechitcoll->size(); ++jhit) {
0367       const EcalRecHit& rechitj = (*ecalrechitcoll)[jhit];
0368       EcalRecHitRef rhRef(ecalrechitcoll, jhit);
0369       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0370
0371       //Ask the hits to be in the same endcap
0372       if (rhposj.z() * rhpos.z() < 0)
0373         continue;
0374
0375       double etaj = rhposj.eta();
0376       double phij = rhposj.phi();
0377       double dr = sqrt((eta - etaj) * (eta - etaj) + deltaPhi(phi, phij) * deltaPhi(phi, phij));
0378
0379       //Outer cone
0380       if (dr > 0.3)
0381         continue;
0382
0383       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0384       //Rechits with et between 1 and 2 GeV are saved in the rh list but not used in the calculation of the halocluster variables
0385       if (rhetj < 1)
0386         continue;
0387       bhrhcandidates.push_back(rhRef);
0388       if (rhetj < 2)
0389         continue;
0390
0391       //Isolation between outer and inner cone
0392       if (dr > 0.05) {
0393         isiso = false;
0394         break;
0395       }  //The deposit should be isolated
0396
0397       etcluster += rhetj;
0398
0399       //Timing infos:
0400       //Here we target both IT and OT beam halo
0401       double rhtj = rechitj.time();
0402
0403       //Discriminating variables for OT beam halo:
0404       if (rhtj > 1)
0405         nbcrystalshight++;
0406       if (rhtj < 0)
0407         nbcrystalssmallt++;
0408       //Timing test (likelihood ratio), only for seeds with large R (100 cm) and for crystals with et>5,
0409       //This targets IT beam halo (t around - 1ns)
0410       if (rhtj > 5) {
0411         double corrt_j = rhtj + sqrt(rhposj.x() * rhposj.x() + rhposj.y() * rhposj.y() + 320. * 320.) / c_cm_per_ns -
0412                          320. / c_cm_per_ns;
0413         //BH is modeled by a Gaussian peaking at 0.
0414         //Collisions is modeled by a Gaussian peaking at 0.3
0415         //The width is similar and taken to be 0.4
0416         timediscriminator += 0.5 * (pow((corrt_j - 0.3) / 0.4, 2) - pow((corrt_j - 0.) / 0.4, 2));
0417         clustersize++;
0418       }
0419     }
0420     //Isolation condition
0421     if (!isiso)
0422       continue;
0423
0424     //Calculate H2/E
0425     //Only second hcal layer is considered as it can happen that a shower initiated in EE reaches HCAL first layer
0426     double h2oe(0);
0427     for (size_t jhit = 0; jhit < hbherechitcoll->size(); ++jhit) {
0428       const HBHERecHit& rechitj = (*hbherechitcoll)[jhit];
0429       math::XYZPoint rhposj = getPosition(rechitj.id(), vtx);
0430
0431       //Ask the hits to be in the same endcap
0432       if (rhposj.z() * rhpos.z() < 0)
0433         continue;
0434       //Selects only second HCAL layer
0435       if (std::abs(rhposj.z()) < 425)
0436         continue;
0437
0438       double rhetj = rechitj.energy() * sqrt(rhposj.perp2() / rhposj.mag2());
0439       if (rhetj < 2)
0440         continue;
0441
0442       double phij = rhposj.phi();
0443       if (std::abs(deltaPhi(phi, phij)) > 0.4)
0444         continue;
0445
0446       double rhrj = sqrt(rhposj.perp2());
0447       if (std::abs(rhr - rhrj) > 50)
0448         continue;
0449
0450       h2oe += rhetj / etcluster;
0451     }
0452     //H/E condition
0453     if (h2oe > 0.1)
0454       continue;
0455
0456     clustercand.setClusterEt(etcluster);
0457     clustercand.setSeedEt(rhet);
0458     clustercand.setSeedEta(eta);
0459     clustercand.setSeedPhi(phi);
0460     clustercand.setSeedZ(rhpos.Z());
0461     clustercand.setSeedR(sqrt(rhpos.perp2()));
0462     clustercand.setSeedTime(rechit.time());
0463     clustercand.setH2overE(h2oe);
0464     clustercand.setNbEarlyCrystals(nbcrystalssmallt);
0465     clustercand.setNbLateCrystals(nbcrystalshight);
0466     clustercand.setClusterSize(clustersize);
0467     clustercand.setTimeDiscriminator(timediscriminator);
0468     clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
0469
0470     bool isbeamhalofrompattern =
0471         EEClusterShapeandTimeStudy_ITBH(clustercand, false) || EEClusterShapeandTimeStudy_OTBH(clustercand, false);
0472     clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
0473
0474     bool isbeamhalofrompattern_hlt =
0475         EEClusterShapeandTimeStudy_ITBH(clustercand, true) || EEClusterShapeandTimeStudy_OTBH(clustercand, true);
0476     clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
0477
0478     TheHaloClusterCandsEE.push_back(clustercand);
0479   }
0480
0481   return TheHaloClusterCandsEE;
0482 }
0483
0484 bool EcalHaloAlgo::EBClusterShapeandTimeStudy(HaloClusterCandidateECAL hcand, bool ishlt) {
0485   //Conditions on the central strip size in eta.
0486   //For low size, extra conditions on seed et, isolation and cluster timing
0487   //The time condition only targets IT beam halo.
0488   //EB rechits from OT beam halos are typically too late (around 5 ns or more) and seem therefore already cleaned by the reconstruction.
0489
0490   if (hcand.getSeedEt() < 5)
0491     return false;
0492   if (hcand.getNbofCrystalsInEta() < 4)
0493     return false;
0494   if (hcand.getNbofCrystalsInEta() == 4 && hcand.getSeedEt() < 10)
0495     return false;
0496   if (hcand.getNbofCrystalsInEta() == 4 && hcand.getEtStripIPhiSeedPlus1() > 0.1 &&
0497       hcand.getEtStripIPhiSeedMinus1() > 0.1)
0498     return false;
0499   if (hcand.getNbofCrystalsInEta() <= 5 && hcand.getTimeDiscriminator() >= 0.)
0500     return false;
0501
0502   //For HLT, only use conditions without timing and tighten seed et condition
0503   if (ishlt && hcand.getNbofCrystalsInEta() <= 5)
0504     return false;
0505   if (ishlt && hcand.getSeedEt() < 10)
0506     return false;
0507
0508   hcand.setIsHaloFromPattern(true);
0509
0510   return true;
0511 }
0512
0513 bool EcalHaloAlgo::EEClusterShapeandTimeStudy_OTBH(HaloClusterCandidateECAL hcand, bool ishlt) {
0514   //Separate conditions targeting IT and OT beam halos
0515   //For OT beam halos, just require enough crystals with large T
0516   if (hcand.getSeedEt() < 20)
0517     return false;
0518   if (hcand.getSeedTime() < 0.5)
0519     return false;
0520   if (hcand.getNbLateCrystals() - hcand.getNbEarlyCrystals() < 2)
0521     return false;
0522
0523   //The use of time information does not allow this method to work at HLT
0524   if (ishlt)
0525     return false;
0526
0527   hcand.setIsHaloFromPattern(true);
0528
0529   return true;
0530 }
0531
0532 bool EcalHaloAlgo::EEClusterShapeandTimeStudy_ITBH(HaloClusterCandidateECAL hcand, bool ishlt) {
0533   //Separate conditions targeting IT and OT beam halos
0534   //For IT beam halos, fakes from collisions are higher => require the cluster size to be small.
0535   //Only halos with R>100 cm are considered here.
0536   //For lower values, the time difference with particles from collisions is too small
0537   //IT outgoing beam halos that interact in EE at low R is probably the most difficult category to deal with:
0538   //Their signature is very close to the one of photon from collisions (similar cluster shape and timing)
0539   if (hcand.getSeedEt() < 20)
0540     return false;
0541   if (hcand.getSeedR() < 100)
0542     return false;
0543   if (hcand.getTimeDiscriminator() < 1)
0544     return false;
0545   if (hcand.getClusterSize() < 2)
0546     return false;
0547   if (hcand.getClusterSize() > 4)
0548     return false;
0549
0550   //The use of time information does not allow this method to work at HLT
0551   if (ishlt)
0552     return false;
0553
0554   hcand.setIsHaloFromPattern(true);
0555
0556   return true;
0557 }
0558
0559 math::XYZPoint EcalHaloAlgo::getPosition(const DetId& id, reco::Vertex::Point vtx) {
0560   const GlobalPoint& pos = geo->getPosition(id);
0561   math::XYZPoint posV(pos.x() - vtx.x(), pos.y() - vtx.y(), pos.z() - vtx.z());
0562   return posV;
0563 }