Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:27

0001 #include "RecoMET/METAlgorithms/interface/GlobalHaloAlgo.h"
0002 namespace {
0003   constexpr float c_cm_per_ns = 29.9792458;
0004 };
0005 /*
0006   [class]:  GlobalHaloAlgo
0007   [authors]: R. Remington, The University of Florida
0008   [description]: See GlobalHaloAlgo.h
0009   [date]: October 15, 2009
0010 */
0011 using namespace std;
0012 using namespace edm;
0013 using namespace reco;
0014 
0015 enum detectorregion { EB, EE, HB, HE };
0016 int Phi_To_HcaliPhi(float phi) {
0017   phi = phi < 0 ? phi + 2. * TMath::Pi() : phi;
0018   float phi_degrees = phi * (360.) / (2. * TMath::Pi());
0019   int iPhi = (int)((phi_degrees / 5.) + 1.);
0020 
0021   return iPhi < 73 ? iPhi : 73;
0022 }
0023 
0024 int Phi_To_EcaliPhi(float phi) {
0025   phi = phi < 0 ? phi + 2. * TMath::Pi() : phi;
0026   float phi_degrees = phi * (360.) / (2. * TMath::Pi());
0027   int iPhi = (int)(phi_degrees + 1.);
0028 
0029   return iPhi < 361 ? iPhi : 360;
0030 }
0031 
0032 GlobalHaloAlgo::GlobalHaloAlgo() {
0033   // Defaults are "loose"
0034   Ecal_R_Min = 110.;  // Tight: 200.
0035   Ecal_R_Max = 330.;  // Tight: 250.
0036   Hcal_R_Min = 110.;  // Tight: 220.
0037   Hcal_R_Max = 490.;  // Tight: 350.
0038 }
0039 
0040 reco::GlobalHaloData GlobalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry,
0041                                                const CSCGeometry& TheCSCGeometry,
0042                                                const reco::CaloMET& TheCaloMET,
0043                                                edm::Handle<edm::View<Candidate> >& TheCaloTowers,
0044                                                edm::Handle<CSCSegmentCollection>& TheCSCSegments,
0045                                                edm::Handle<CSCRecHit2DCollection>& TheCSCRecHits,
0046                                                edm::Handle<reco::MuonCollection>& TheMuons,
0047                                                const CSCHaloData& TheCSCHaloData,
0048                                                const EcalHaloData& TheEcalHaloData,
0049                                                const HcalHaloData& TheHcalHaloData,
0050                                                bool ishlt) {
0051   GlobalHaloData TheGlobalHaloData;
0052   float METOverSumEt = TheCaloMET.sumEt() ? TheCaloMET.pt() / TheCaloMET.sumEt() : 0;
0053   TheGlobalHaloData.SetMETOverSumEt(METOverSumEt);
0054 
0055   int EcalOverlapping_CSCRecHits[361] = {};
0056   int EcalOverlapping_CSCSegments[361] = {};
0057   int HcalOverlapping_CSCRecHits[73] = {};
0058   int HcalOverlapping_CSCSegments[73] = {};
0059 
0060   if (TheCSCSegments.isValid()) {
0061     for (CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end();
0062          iSegment++) {
0063       bool EcalOverlap[361];
0064       bool HcalOverlap[73];
0065       for (int i = 0; i < 361; i++) {
0066         EcalOverlap[i] = false;
0067         if (i < 73)
0068           HcalOverlap[i] = false;
0069       }
0070 
0071       std::vector<CSCRecHit2D> Hits = iSegment->specificRecHits();
0072       for (std::vector<CSCRecHit2D>::iterator iHit = Hits.begin(); iHit != Hits.end(); iHit++) {
0073         DetId TheDetUnitId(iHit->geographicalId());
0074         if (TheDetUnitId.det() != DetId::Muon)
0075           continue;
0076         if (TheDetUnitId.subdetId() != MuonSubdetId::CSC)
0077           continue;
0078 
0079         const GeomDetUnit* TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
0080         LocalPoint TheLocalPosition = iHit->localPosition();
0081         const BoundPlane& TheSurface = TheUnit->surface();
0082         const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);
0083 
0084         int Hcal_iphi = Phi_To_HcaliPhi(TheGlobalPosition.phi());
0085         int Ecal_iphi = Phi_To_EcaliPhi(TheGlobalPosition.phi());
0086         float x = TheGlobalPosition.x();
0087         float y = TheGlobalPosition.y();
0088 
0089         float r = TMath::Sqrt(x * x + y * y);
0090 
0091         if (r < Ecal_R_Max && r > Ecal_R_Min)
0092           EcalOverlap[Ecal_iphi] = true;
0093         if (r < Hcal_R_Max && r > Hcal_R_Max)
0094           HcalOverlap[Hcal_iphi] = true;
0095       }
0096       for (int i = 0; i < 361; i++) {
0097         if (EcalOverlap[i])
0098           EcalOverlapping_CSCSegments[i]++;
0099         if (i < 73 && HcalOverlap[i])
0100           HcalOverlapping_CSCSegments[i]++;
0101       }
0102     }
0103   }
0104   if (TheCSCRecHits.isValid()) {
0105     for (CSCRecHit2DCollection::const_iterator iCSCRecHit = TheCSCRecHits->begin(); iCSCRecHit != TheCSCRecHits->end();
0106          iCSCRecHit++) {
0107       DetId TheDetUnitId(iCSCRecHit->geographicalId());
0108       if (TheDetUnitId.det() != DetId::Muon)
0109         continue;
0110       if (TheDetUnitId.subdetId() != MuonSubdetId::CSC)
0111         continue;
0112 
0113       const GeomDetUnit* TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
0114       LocalPoint TheLocalPosition = iCSCRecHit->localPosition();
0115       const BoundPlane& TheSurface = TheUnit->surface();
0116       const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);
0117 
0118       int Hcaliphi = Phi_To_HcaliPhi(TheGlobalPosition.phi());
0119       int Ecaliphi = Phi_To_EcaliPhi(TheGlobalPosition.phi());
0120       float x = TheGlobalPosition.x();
0121       float y = TheGlobalPosition.y();
0122 
0123       float r = TMath::Sqrt(x * x + y * y);
0124 
0125       if (r < Ecal_R_Max && r > Ecal_R_Min)
0126         EcalOverlapping_CSCRecHits[Ecaliphi]++;
0127       if (r < Hcal_R_Max && r > Hcal_R_Max)
0128         HcalOverlapping_CSCRecHits[Hcaliphi]++;
0129     }
0130   }
0131 
0132   // In development....
0133   // Get Ecal Wedges
0134   std::vector<PhiWedge> EcalWedges = TheEcalHaloData.GetPhiWedges();
0135 
0136   // Get Hcal Wedges
0137   std::vector<PhiWedge> HcalWedges = TheHcalHaloData.GetPhiWedges();
0138 
0139   //Get Ref to CSC Tracks
0140   //edm::RefVector<reco::TrackCollection> TheCSCTracks = TheCSCHaloData.GetTracks();
0141   //for(unsigned int i = 0 ; i < TheCSCTracks.size() ; i++ )
0142   //edm::Ref<reco::TrackCollection> iTrack( TheCSCTracks, i );
0143 
0144   // Get global positions of central most rechit of CSC Halo tracks
0145   std::vector<GlobalPoint> TheGlobalPositions = TheCSCHaloData.GetCSCTrackImpactPositions();
0146 
0147   // Container to store Ecal/Hcal iPhi values matched to impact point of CSC tracks
0148   std::vector<int> vEcaliPhi, vHcaliPhi;
0149 
0150   // Keep track of number of calo pointing CSC halo tracks that do not match to Phi wedges
0151   int N_Unmatched_Tracks = 0;
0152 
0153   for (std::vector<GlobalPoint>::iterator Pos = TheGlobalPositions.begin(); Pos != TheGlobalPositions.end(); Pos++) {
0154     // Calculate global phi coordinate for central most rechit in the track
0155     float global_phi = Pos->phi();
0156     float global_r = TMath::Sqrt(Pos->x() * Pos->x() + Pos->y() * Pos->y());
0157 
0158     // Convert global phi to iPhi
0159     int global_EcaliPhi = Phi_To_EcaliPhi(global_phi);
0160     int global_HcaliPhi = Phi_To_HcaliPhi(global_phi);
0161     bool MATCHED = false;
0162 
0163     //Loop over Ecal Phi Wedges
0164     for (std::vector<PhiWedge>::iterator iWedge = EcalWedges.begin(); iWedge != EcalWedges.end(); iWedge++) {
0165       if ((TMath::Abs(global_EcaliPhi - iWedge->iPhi()) <= 5) && (global_r > Ecal_R_Min && global_r < Ecal_R_Max)) {
0166         bool StoreWedge = true;
0167         for (unsigned int i = 0; i < vEcaliPhi.size(); i++)
0168           if (vEcaliPhi[i] == iWedge->iPhi())
0169             StoreWedge = false;
0170 
0171         if (StoreWedge) {
0172           PhiWedge NewWedge(*iWedge);
0173           NewWedge.SetOverlappingCSCSegments(EcalOverlapping_CSCSegments[iWedge->iPhi()]);
0174           NewWedge.SetOverlappingCSCRecHits(EcalOverlapping_CSCRecHits[iWedge->iPhi()]);
0175           vEcaliPhi.push_back(iWedge->iPhi());
0176           TheGlobalHaloData.GetMatchedEcalPhiWedges().push_back(NewWedge);
0177         }
0178         MATCHED = true;
0179       }
0180     }
0181     //Loop over Hcal Phi Wedges
0182     for (std::vector<PhiWedge>::iterator iWedge = HcalWedges.begin(); iWedge != HcalWedges.end(); iWedge++) {
0183       if ((TMath::Abs(global_HcaliPhi - iWedge->iPhi()) <= 2) && (global_r > Hcal_R_Min && global_r < Hcal_R_Max)) {
0184         bool StoreWedge = true;
0185         for (unsigned int i = 0; i < vHcaliPhi.size(); i++)
0186           if (vHcaliPhi[i] == iWedge->iPhi())
0187             StoreWedge = false;
0188 
0189         if (StoreWedge) {
0190           vHcaliPhi.push_back(iWedge->iPhi());
0191           PhiWedge NewWedge(*iWedge);
0192           NewWedge.SetOverlappingCSCSegments(HcalOverlapping_CSCSegments[iWedge->iPhi()]);
0193           NewWedge.SetOverlappingCSCRecHits(HcalOverlapping_CSCRecHits[iWedge->iPhi()]);
0194           PhiWedge wedge(*iWedge);
0195           TheGlobalHaloData.GetMatchedHcalPhiWedges().push_back(NewWedge);
0196         }
0197         MATCHED = true;
0198       }
0199     }
0200     if (!MATCHED)
0201       N_Unmatched_Tracks++;
0202   }
0203 
0204   // Corrections to MEx, MEy
0205   float dMEx = 0.;
0206   float dMEy = 0.;
0207   // Loop over calotowers and correct the MET for the towers that lie in the trajectory of the CSC Halo Tracks
0208   for (edm::View<Candidate>::const_iterator iCandidate = TheCaloTowers->begin(); iCandidate != TheCaloTowers->end();
0209        iCandidate++) {
0210     const Candidate* c = &(*iCandidate);
0211     if (c) {
0212       const CaloTower* iTower = dynamic_cast<const CaloTower*>(c);
0213       if (iTower->et() < TowerEtThreshold)
0214         continue;
0215       if (abs(iTower->ieta()) > 24)
0216         continue;  // not in barrel/endcap
0217       int iphi = iTower->iphi();
0218       for (unsigned int x = 0; x < vEcaliPhi.size(); x++) {
0219         if (iphi == vEcaliPhi[x]) {
0220           dMEx += (TMath::Cos(iTower->phi()) * iTower->emEt());
0221           dMEy += (TMath::Sin(iTower->phi()) * iTower->emEt());
0222         }
0223       }
0224       for (unsigned int x = 0; x < vHcaliPhi.size(); x++) {
0225         if (iphi == vHcaliPhi[x]) {
0226           dMEx += (TMath::Cos(iTower->phi()) * iTower->hadEt());
0227           dMEy += (TMath::Sin(iTower->phi()) * iTower->hadEt());
0228         }
0229       }
0230     }
0231   }
0232 
0233   TheGlobalHaloData.SetMETCorrections(dMEx, dMEy);
0234 
0235   std::vector<HaloClusterCandidateECAL> hccandEB = TheEcalHaloData.getHaloClusterCandidatesEB();
0236   std::vector<HaloClusterCandidateECAL> hccandEE = TheEcalHaloData.getHaloClusterCandidatesEE();
0237   std::vector<HaloClusterCandidateHCAL> hccandHB = TheHcalHaloData.getHaloClusterCandidatesHB();
0238   std::vector<HaloClusterCandidateHCAL> hccandHE = TheHcalHaloData.getHaloClusterCandidatesHE();
0239 
0240   //CSC-calo matching
0241   bool ECALBmatched(false), ECALEmatched(false), HCALBmatched(false), HCALEmatched(false);
0242 
0243   if (TheCSCSegments.isValid()) {
0244     for (CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end();
0245          iSegment++) {
0246       CSCDetId iCscDetID = iSegment->cscDetId();
0247       bool Segment1IsGood = true;
0248 
0249       //avoid segments from collision muons
0250       if (TheMuons.isValid()) {
0251         for (reco::MuonCollection::const_iterator mu = TheMuons->begin(); mu != TheMuons->end() && (Segment1IsGood);
0252              mu++) {
0253           if (!mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon())
0254             continue;
0255           if (!mu->isGlobalMuon() && mu->isTrackerMuon() && mu->pt() < 3)
0256             continue;
0257           const std::vector<MuonChamberMatch> chambers = mu->matches();
0258           for (std::vector<MuonChamberMatch>::const_iterator kChamber = chambers.begin(); kChamber != chambers.end();
0259                kChamber++) {
0260             if (kChamber->detector() != MuonSubdetId::CSC)
0261               continue;
0262             for (std::vector<reco::MuonSegmentMatch>::const_iterator kSegment = kChamber->segmentMatches.begin();
0263                  kSegment != kChamber->segmentMatches.end();
0264                  kSegment++) {
0265               edm::Ref<CSCSegmentCollection> cscSegRef = kSegment->cscSegmentRef;
0266               CSCDetId kCscDetID = cscSegRef->cscDetId();
0267 
0268               if (kCscDetID == iCscDetID) {
0269                 Segment1IsGood = false;
0270               }
0271             }
0272           }
0273         }
0274       }
0275       if (!Segment1IsGood)
0276         continue;
0277 
0278       // Get local direction vector; if direction runs parallel to beamline,
0279       // count this segment as beam halo candidate.
0280       LocalPoint iLocalPosition = iSegment->localPosition();
0281       LocalVector iLocalDirection = iSegment->localDirection();
0282 
0283       GlobalPoint iGlobalPosition = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalPosition);
0284       GlobalVector iGlobalDirection = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalDirection);
0285 
0286       float iTheta = iGlobalDirection.theta();
0287       if (iTheta > max_segment_theta && iTheta < TMath::Pi() - max_segment_theta)
0288         continue;
0289 
0290       float iPhi = iGlobalPosition.phi();
0291       float iR = sqrt(iGlobalPosition.perp2());
0292       float iZ = iGlobalPosition.z();
0293       float iT = iSegment->time();
0294 
0295       //CSC-calo matching:
0296       //Here, one checks if any halo cluster can be matched to a CSC segment.
0297       //The matching uses both geometric (dphi, dR) and timing information (dt).
0298       //The cut values depend on the subdetector considered (e.g. in HB, Rcalo-Rsegment is allowed to be very negative)
0299 
0300       bool ebmatched = SegmentMatchingEB(TheGlobalHaloData, hccandEB, iZ, iR, iT, iPhi, ishlt);
0301       bool eematched = SegmentMatchingEE(TheGlobalHaloData, hccandEE, iZ, iR, iT, iPhi, ishlt);
0302       bool hbmatched = SegmentMatchingHB(TheGlobalHaloData, hccandHB, iZ, iR, iT, iPhi, ishlt);
0303       bool hematched = SegmentMatchingHE(TheGlobalHaloData, hccandHE, iZ, iR, iT, iPhi, ishlt);
0304 
0305       ECALBmatched |= ebmatched;
0306       ECALEmatched |= eematched;
0307       HCALBmatched |= hbmatched;
0308       HCALEmatched |= hematched;
0309     }
0310   }
0311 
0312   TheGlobalHaloData.SetSegmentIsEBCaloMatched(ECALBmatched);
0313   TheGlobalHaloData.SetSegmentIsEECaloMatched(ECALEmatched);
0314   TheGlobalHaloData.SetSegmentIsHBCaloMatched(HCALBmatched);
0315   TheGlobalHaloData.SetSegmentIsHECaloMatched(HCALEmatched);
0316 
0317   //Now checking patterns from EcalHaloData and HcalHaloData:
0318   //Simply check whether any cluster has a halo pattern
0319   //In that case store the rhits in GlobalHaloData
0320 
0321   bool HaloPatternFoundInEB = false;
0322   for (auto& hcand : hccandEB) {
0323     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0324       HaloPatternFoundInEB = true;
0325       edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0326       AddtoBeamHaloEBEERechits(bhrhcandidates, TheGlobalHaloData, true);
0327     }
0328   }
0329 
0330   bool HaloPatternFoundInEE = false;
0331   for (auto& hcand : hccandEE) {
0332     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0333       HaloPatternFoundInEE = true;
0334       edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0335       AddtoBeamHaloEBEERechits(bhrhcandidates, TheGlobalHaloData, false);
0336     }
0337   }
0338 
0339   bool HaloPatternFoundInHB = false;
0340   for (auto& hcand : hccandHB) {
0341     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0342       HaloPatternFoundInHB = true;
0343       edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0344       AddtoBeamHaloHBHERechits(bhrhcandidates, TheGlobalHaloData);
0345     }
0346   }
0347 
0348   bool HaloPatternFoundInHE = false;
0349   for (auto& hcand : hccandHE) {
0350     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0351       HaloPatternFoundInHE = true;
0352       edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0353       AddtoBeamHaloHBHERechits(bhrhcandidates, TheGlobalHaloData);
0354     }
0355   }
0356   TheGlobalHaloData.SetHaloPatternFoundEB(HaloPatternFoundInEB);
0357   TheGlobalHaloData.SetHaloPatternFoundEE(HaloPatternFoundInEE);
0358   TheGlobalHaloData.SetHaloPatternFoundHB(HaloPatternFoundInHB);
0359   TheGlobalHaloData.SetHaloPatternFoundHE(HaloPatternFoundInHE);
0360 
0361   return TheGlobalHaloData;
0362 }
0363 
0364 bool GlobalHaloAlgo::SegmentMatchingEB(reco::GlobalHaloData& thehalodata,
0365                                        const std::vector<HaloClusterCandidateECAL>& haloclustercands,
0366                                        float iZ,
0367                                        float iR,
0368                                        float iT,
0369                                        float iPhi,
0370                                        bool ishlt) {
0371   bool rhmatchingfound = false;
0372 
0373   for (auto& hcand : haloclustercands) {
0374     if (!ApplyMatchingCuts(EB,
0375                            ishlt,
0376                            hcand.getSeedEt(),
0377                            iZ,
0378                            hcand.getSeedZ(),
0379                            iR,
0380                            hcand.getSeedR(),
0381                            iT,
0382                            hcand.getSeedTime(),
0383                            iPhi,
0384                            hcand.getSeedPhi()))
0385       continue;
0386 
0387     rhmatchingfound = true;
0388 
0389     edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0390 
0391     AddtoBeamHaloEBEERechits(bhrhcandidates, thehalodata, true);
0392   }
0393 
0394   return rhmatchingfound;
0395 }
0396 
0397 bool GlobalHaloAlgo::SegmentMatchingEE(reco::GlobalHaloData& thehalodata,
0398                                        const std::vector<HaloClusterCandidateECAL>& haloclustercands,
0399                                        float iZ,
0400                                        float iR,
0401                                        float iT,
0402                                        float iPhi,
0403                                        bool ishlt) {
0404   bool rhmatchingfound = false;
0405 
0406   for (auto& hcand : haloclustercands) {
0407     if (!ApplyMatchingCuts(EE,
0408                            ishlt,
0409                            hcand.getSeedEt(),
0410                            iZ,
0411                            hcand.getSeedZ(),
0412                            iR,
0413                            hcand.getSeedR(),
0414                            iT,
0415                            hcand.getSeedTime(),
0416                            iPhi,
0417                            hcand.getSeedPhi()))
0418       continue;
0419 
0420     rhmatchingfound = true;
0421 
0422     edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0423 
0424     AddtoBeamHaloEBEERechits(bhrhcandidates, thehalodata, false);
0425   }
0426 
0427   return rhmatchingfound;
0428 }
0429 
0430 bool GlobalHaloAlgo::SegmentMatchingHB(reco::GlobalHaloData& thehalodata,
0431                                        const std::vector<HaloClusterCandidateHCAL>& haloclustercands,
0432                                        float iZ,
0433                                        float iR,
0434                                        float iT,
0435                                        float iPhi,
0436                                        bool ishlt) {
0437   bool rhmatchingfound = false;
0438 
0439   for (auto& hcand : haloclustercands) {
0440     if (!ApplyMatchingCuts(HB,
0441                            ishlt,
0442                            hcand.getSeedEt(),
0443                            iZ,
0444                            hcand.getSeedZ(),
0445                            iR,
0446                            hcand.getSeedR(),
0447                            iT,
0448                            hcand.getSeedTime(),
0449                            iPhi,
0450                            hcand.getSeedPhi()))
0451       continue;
0452 
0453     rhmatchingfound = true;
0454 
0455     edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0456 
0457     AddtoBeamHaloHBHERechits(bhrhcandidates, thehalodata);
0458   }
0459 
0460   return rhmatchingfound;
0461 }
0462 
0463 bool GlobalHaloAlgo::SegmentMatchingHE(reco::GlobalHaloData& thehalodata,
0464                                        const std::vector<HaloClusterCandidateHCAL>& haloclustercands,
0465                                        float iZ,
0466                                        float iR,
0467                                        float iT,
0468                                        float iPhi,
0469                                        bool ishlt) {
0470   bool rhmatchingfound = false;
0471 
0472   for (auto& hcand : haloclustercands) {
0473     if (!ApplyMatchingCuts(HE,
0474                            ishlt,
0475                            hcand.getSeedEt(),
0476                            iZ,
0477                            hcand.getSeedZ(),
0478                            iR,
0479                            hcand.getSeedR(),
0480                            iT,
0481                            hcand.getSeedTime(),
0482                            iPhi,
0483                            hcand.getSeedPhi()))
0484       continue;
0485 
0486     rhmatchingfound = true;
0487 
0488     edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0489 
0490     AddtoBeamHaloHBHERechits(bhrhcandidates, thehalodata);
0491   }
0492 
0493   return rhmatchingfound;
0494 }
0495 
0496 bool GlobalHaloAlgo::ApplyMatchingCuts(int subdet,
0497                                        bool ishlt,
0498                                        double rhet,
0499                                        double segZ,
0500                                        double rhZ,
0501                                        double segR,
0502                                        double rhR,
0503                                        double segT,
0504                                        double rhT,
0505                                        double segPhi,
0506                                        double rhPhi) {
0507   //Std::Absolute time wrt BX
0508   double tBXrh = rhT + sqrt(rhR * rhR + rhZ * rhZ) / c_cm_per_ns;
0509   double tBXseg = segT + sqrt(segR * segR + segZ * segZ) / c_cm_per_ns;
0510   //Time at z=0, under beam halo hypothesis
0511   double tcorseg = tBXseg - std::abs(segZ) / c_cm_per_ns;       //Outgoing beam halo
0512   double tcorsegincbh = tBXseg + std::abs(segZ) / c_cm_per_ns;  //Ingoing beam halo
0513   double truedt[4] = {1000, 1000, 1000, 1000};
0514   //There are four types of segments associated to beam halo, test each hypothesis:
0515   //IT beam halo, ingoing track
0516   double twindow_seg = 15;
0517   if (std::abs(tcorsegincbh) < twindow_seg)
0518     truedt[0] = tBXrh - tBXseg - std::abs(rhZ - segZ) / c_cm_per_ns;
0519   //IT beam halo, outgoing track
0520   if (std::abs(tcorseg) < twindow_seg)
0521     truedt[1] = tBXseg - tBXrh - std::abs(rhZ - segZ) / c_cm_per_ns;
0522   //OT beam halo (from next BX), ingoing track
0523   if (tcorsegincbh > 25 - twindow_seg && std::abs(tcorsegincbh) < 25 + twindow_seg)
0524     truedt[2] = tBXrh - tBXseg - std::abs(rhZ - segZ) / c_cm_per_ns;
0525   //OT beam halo (from next BX), outgoing track
0526   if (tcorseg > 25 - twindow_seg && tcorseg < 25 + twindow_seg)
0527     truedt[3] = tBXseg - tBXrh - std::abs(rhZ - segZ) / c_cm_per_ns;
0528 
0529   if (subdet == EB) {
0530     if (rhet < et_thresh_rh_eb)
0531       return false;
0532     if (rhet < 20 && ishlt)
0533       return false;
0534     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_eb)
0535       return false;
0536     if (rhR - segR < dr_lowthresh_segvsrh_eb)
0537       return false;
0538     if (rhR - segR > dr_highthresh_segvsrh_eb)
0539       return false;
0540     if (std::abs(truedt[0]) > dt_segvsrh_eb && std::abs(truedt[1]) > dt_segvsrh_eb &&
0541         std::abs(truedt[2]) > dt_segvsrh_eb && std::abs(truedt[3]) > dt_segvsrh_eb)
0542       return false;
0543     return true;
0544   }
0545 
0546   if (subdet == EE) {
0547     if (rhet < et_thresh_rh_ee)
0548       return false;
0549     if (rhet < 20 && ishlt)
0550       return false;
0551     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_ee)
0552       return false;
0553     if (rhR - segR < dr_lowthresh_segvsrh_ee)
0554       return false;
0555     if (rhR - segR > dr_highthresh_segvsrh_ee)
0556       return false;
0557     if (std::abs(truedt[0]) > dt_segvsrh_ee && std::abs(truedt[1]) > dt_segvsrh_ee &&
0558         std::abs(truedt[2]) > dt_segvsrh_ee && std::abs(truedt[3]) > dt_segvsrh_ee)
0559       return false;
0560     return true;
0561   }
0562 
0563   if (subdet == HB) {
0564     if (rhet < et_thresh_rh_hb)
0565       return false;
0566     if (rhet < 20 && ishlt)
0567       return false;
0568     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_hb)
0569       return false;
0570     if (rhR - segR < dr_lowthresh_segvsrh_hb)
0571       return false;
0572     if (rhR - segR > dr_highthresh_segvsrh_hb)
0573       return false;
0574     if (std::abs(truedt[0]) > dt_segvsrh_hb && std::abs(truedt[1]) > dt_segvsrh_hb &&
0575         std::abs(truedt[2]) > dt_segvsrh_hb && std::abs(truedt[3]) > dt_segvsrh_hb)
0576       return false;
0577     return true;
0578   }
0579 
0580   if (subdet == HE) {
0581     if (rhet < et_thresh_rh_he)
0582       return false;
0583     if (rhet < 20 && ishlt)
0584       return false;
0585     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_he)
0586       return false;
0587     if (rhR - segR < dr_lowthresh_segvsrh_he)
0588       return false;
0589     if (rhR - segR > dr_highthresh_segvsrh_he)
0590       return false;
0591     if (std::abs(truedt[0]) > dt_segvsrh_he && std::abs(truedt[1]) > dt_segvsrh_he &&
0592         std::abs(truedt[2]) > dt_segvsrh_he && std::abs(truedt[3]) > dt_segvsrh_he)
0593       return false;
0594     return true;
0595   }
0596 
0597   return false;
0598 }
0599 
0600 void GlobalHaloAlgo::AddtoBeamHaloEBEERechits(edm::RefVector<EcalRecHitCollection>& bhtaggedrechits,
0601                                               reco::GlobalHaloData& thehalodata,
0602                                               bool isbarrel) {
0603   for (size_t ihit = 0; ihit < bhtaggedrechits.size(); ++ihit) {
0604     bool alreadyincl = false;
0605     edm::Ref<EcalRecHitCollection> rhRef(bhtaggedrechits[ihit]);
0606     edm::RefVector<EcalRecHitCollection> refrhcoll;
0607     if (isbarrel)
0608       refrhcoll = thehalodata.GetEBRechits();
0609     else
0610       refrhcoll = thehalodata.GetEERechits();
0611     for (size_t jhit = 0; jhit < refrhcoll.size(); jhit++) {
0612       edm::Ref<EcalRecHitCollection> rhitRef(refrhcoll[jhit]);
0613       if (rhitRef->detid() == rhRef->detid())
0614         alreadyincl = true;
0615       if (rhitRef->detid() == rhRef->detid())
0616         break;
0617     }
0618     if (!alreadyincl && isbarrel)
0619       thehalodata.GetEBRechits().push_back(rhRef);
0620     if (!alreadyincl && !isbarrel)
0621       thehalodata.GetEERechits().push_back(rhRef);
0622   }
0623 }
0624 
0625 void GlobalHaloAlgo::AddtoBeamHaloHBHERechits(edm::RefVector<HBHERecHitCollection>& bhtaggedrechits,
0626                                               reco::GlobalHaloData& thehalodata) {
0627   for (size_t ihit = 0; ihit < bhtaggedrechits.size(); ++ihit) {
0628     bool alreadyincl = false;
0629     edm::Ref<HBHERecHitCollection> rhRef(bhtaggedrechits[ihit]);
0630     edm::RefVector<HBHERecHitCollection> refrhcoll;
0631     refrhcoll = thehalodata.GetHBHERechits();
0632     for (size_t jhit = 0; jhit < refrhcoll.size(); jhit++) {
0633       edm::Ref<HBHERecHitCollection> rhitRef(refrhcoll[jhit]);
0634       if (rhitRef->detid() == rhRef->detid())
0635         alreadyincl = true;
0636       if (rhitRef->detid() == rhRef->detid())
0637         break;
0638     }
0639     if (!alreadyincl)
0640       thehalodata.GetHBHERechits().push_back(rhRef);
0641   }
0642 }