Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:00:56 

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   for (std::vector<GlobalPoint>::iterator Pos = TheGlobalPositions.begin(); Pos != TheGlobalPositions.end(); Pos++) {
0151     // Calculate global phi coordinate for central most rechit in the track
0152     float global_phi = Pos->phi();
0153     float global_r = TMath::Sqrt(Pos->x() * Pos->x() + Pos->y() * Pos->y());
0154 
0155     // Convert global phi to iPhi
0156     int global_EcaliPhi = Phi_To_EcaliPhi(global_phi);
0157     int global_HcaliPhi = Phi_To_HcaliPhi(global_phi);
0158 
0159     //Loop over Ecal Phi Wedges
0160     for (std::vector<PhiWedge>::iterator iWedge = EcalWedges.begin(); iWedge != EcalWedges.end(); iWedge++) {
0161       if ((TMath::Abs(global_EcaliPhi - iWedge->iPhi()) <= 5) && (global_r > Ecal_R_Min && global_r < Ecal_R_Max)) {
0162         bool StoreWedge = true;
0163         for (unsigned int i = 0; i < vEcaliPhi.size(); i++)
0164           if (vEcaliPhi[i] == iWedge->iPhi())
0165             StoreWedge = false;
0166 
0167         if (StoreWedge) {
0168           PhiWedge NewWedge(*iWedge);
0169           NewWedge.SetOverlappingCSCSegments(EcalOverlapping_CSCSegments[iWedge->iPhi()]);
0170           NewWedge.SetOverlappingCSCRecHits(EcalOverlapping_CSCRecHits[iWedge->iPhi()]);
0171           vEcaliPhi.push_back(iWedge->iPhi());
0172           TheGlobalHaloData.GetMatchedEcalPhiWedges().push_back(NewWedge);
0173         }
0174       }
0175     }
0176     //Loop over Hcal Phi Wedges
0177     for (std::vector<PhiWedge>::iterator iWedge = HcalWedges.begin(); iWedge != HcalWedges.end(); iWedge++) {
0178       if ((TMath::Abs(global_HcaliPhi - iWedge->iPhi()) <= 2) && (global_r > Hcal_R_Min && global_r < Hcal_R_Max)) {
0179         bool StoreWedge = true;
0180         for (unsigned int i = 0; i < vHcaliPhi.size(); i++)
0181           if (vHcaliPhi[i] == iWedge->iPhi())
0182             StoreWedge = false;
0183 
0184         if (StoreWedge) {
0185           vHcaliPhi.push_back(iWedge->iPhi());
0186           PhiWedge NewWedge(*iWedge);
0187           NewWedge.SetOverlappingCSCSegments(HcalOverlapping_CSCSegments[iWedge->iPhi()]);
0188           NewWedge.SetOverlappingCSCRecHits(HcalOverlapping_CSCRecHits[iWedge->iPhi()]);
0189           PhiWedge wedge(*iWedge);
0190           TheGlobalHaloData.GetMatchedHcalPhiWedges().push_back(NewWedge);
0191         }
0192       }
0193     }
0194   }
0195 
0196   // Corrections to MEx, MEy
0197   float dMEx = 0.;
0198   float dMEy = 0.;
0199   // Loop over calotowers and correct the MET for the towers that lie in the trajectory of the CSC Halo Tracks
0200   for (edm::View<Candidate>::const_iterator iCandidate = TheCaloTowers->begin(); iCandidate != TheCaloTowers->end();
0201        iCandidate++) {
0202     const Candidate* c = &(*iCandidate);
0203     if (c) {
0204       const CaloTower* iTower = dynamic_cast<const CaloTower*>(c);
0205       if (iTower->et() < TowerEtThreshold)
0206         continue;
0207       if (abs(iTower->ieta()) > 24)
0208         continue;  // not in barrel/endcap
0209       int iphi = iTower->iphi();
0210       for (unsigned int x = 0; x < vEcaliPhi.size(); x++) {
0211         if (iphi == vEcaliPhi[x]) {
0212           dMEx += (TMath::Cos(iTower->phi()) * iTower->emEt());
0213           dMEy += (TMath::Sin(iTower->phi()) * iTower->emEt());
0214         }
0215       }
0216       for (unsigned int x = 0; x < vHcaliPhi.size(); x++) {
0217         if (iphi == vHcaliPhi[x]) {
0218           dMEx += (TMath::Cos(iTower->phi()) * iTower->hadEt());
0219           dMEy += (TMath::Sin(iTower->phi()) * iTower->hadEt());
0220         }
0221       }
0222     }
0223   }
0224 
0225   TheGlobalHaloData.SetMETCorrections(dMEx, dMEy);
0226 
0227   std::vector<HaloClusterCandidateECAL> hccandEB = TheEcalHaloData.getHaloClusterCandidatesEB();
0228   std::vector<HaloClusterCandidateECAL> hccandEE = TheEcalHaloData.getHaloClusterCandidatesEE();
0229   std::vector<HaloClusterCandidateHCAL> hccandHB = TheHcalHaloData.getHaloClusterCandidatesHB();
0230   std::vector<HaloClusterCandidateHCAL> hccandHE = TheHcalHaloData.getHaloClusterCandidatesHE();
0231 
0232   //CSC-calo matching
0233   bool ECALBmatched(false), ECALEmatched(false), HCALBmatched(false), HCALEmatched(false);
0234 
0235   if (TheCSCSegments.isValid()) {
0236     for (CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end();
0237          iSegment++) {
0238       CSCDetId iCscDetID = iSegment->cscDetId();
0239       bool Segment1IsGood = true;
0240 
0241       //avoid segments from collision muons
0242       if (TheMuons.isValid()) {
0243         for (reco::MuonCollection::const_iterator mu = TheMuons->begin(); mu != TheMuons->end() && (Segment1IsGood);
0244              mu++) {
0245           if (!mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon())
0246             continue;
0247           if (!mu->isGlobalMuon() && mu->isTrackerMuon() && mu->pt() < 3)
0248             continue;
0249           const std::vector<MuonChamberMatch> chambers = mu->matches();
0250           for (std::vector<MuonChamberMatch>::const_iterator kChamber = chambers.begin(); kChamber != chambers.end();
0251                kChamber++) {
0252             if (kChamber->detector() != MuonSubdetId::CSC)
0253               continue;
0254             for (std::vector<reco::MuonSegmentMatch>::const_iterator kSegment = kChamber->segmentMatches.begin();
0255                  kSegment != kChamber->segmentMatches.end();
0256                  kSegment++) {
0257               edm::Ref<CSCSegmentCollection> cscSegRef = kSegment->cscSegmentRef;
0258               CSCDetId kCscDetID = cscSegRef->cscDetId();
0259 
0260               if (kCscDetID == iCscDetID) {
0261                 Segment1IsGood = false;
0262               }
0263             }
0264           }
0265         }
0266       }
0267       if (!Segment1IsGood)
0268         continue;
0269 
0270       // Get local direction vector; if direction runs parallel to beamline,
0271       // count this segment as beam halo candidate.
0272       LocalPoint iLocalPosition = iSegment->localPosition();
0273       LocalVector iLocalDirection = iSegment->localDirection();
0274 
0275       GlobalPoint iGlobalPosition = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalPosition);
0276       GlobalVector iGlobalDirection = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalDirection);
0277 
0278       float iTheta = iGlobalDirection.theta();
0279       if (iTheta > max_segment_theta && iTheta < TMath::Pi() - max_segment_theta)
0280         continue;
0281 
0282       float iPhi = iGlobalPosition.phi();
0283       float iR = sqrt(iGlobalPosition.perp2());
0284       float iZ = iGlobalPosition.z();
0285       float iT = iSegment->time();
0286 
0287       //CSC-calo matching:
0288       //Here, one checks if any halo cluster can be matched to a CSC segment.
0289       //The matching uses both geometric (dphi, dR) and timing information (dt).
0290       //The cut values depend on the subdetector considered (e.g. in HB, Rcalo-Rsegment is allowed to be very negative)
0291 
0292       bool ebmatched = SegmentMatchingEB(TheGlobalHaloData, hccandEB, iZ, iR, iT, iPhi, ishlt);
0293       bool eematched = SegmentMatchingEE(TheGlobalHaloData, hccandEE, iZ, iR, iT, iPhi, ishlt);
0294       bool hbmatched = SegmentMatchingHB(TheGlobalHaloData, hccandHB, iZ, iR, iT, iPhi, ishlt);
0295       bool hematched = SegmentMatchingHE(TheGlobalHaloData, hccandHE, iZ, iR, iT, iPhi, ishlt);
0296 
0297       ECALBmatched |= ebmatched;
0298       ECALEmatched |= eematched;
0299       HCALBmatched |= hbmatched;
0300       HCALEmatched |= hematched;
0301     }
0302   }
0303 
0304   TheGlobalHaloData.SetSegmentIsEBCaloMatched(ECALBmatched);
0305   TheGlobalHaloData.SetSegmentIsEECaloMatched(ECALEmatched);
0306   TheGlobalHaloData.SetSegmentIsHBCaloMatched(HCALBmatched);
0307   TheGlobalHaloData.SetSegmentIsHECaloMatched(HCALEmatched);
0308 
0309   //Now checking patterns from EcalHaloData and HcalHaloData:
0310   //Simply check whether any cluster has a halo pattern
0311   //In that case store the rhits in GlobalHaloData
0312 
0313   bool HaloPatternFoundInEB = false;
0314   for (auto& hcand : hccandEB) {
0315     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0316       HaloPatternFoundInEB = true;
0317       edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0318       AddtoBeamHaloEBEERechits(bhrhcandidates, TheGlobalHaloData, true);
0319     }
0320   }
0321 
0322   bool HaloPatternFoundInEE = false;
0323   for (auto& hcand : hccandEE) {
0324     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0325       HaloPatternFoundInEE = true;
0326       edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0327       AddtoBeamHaloEBEERechits(bhrhcandidates, TheGlobalHaloData, false);
0328     }
0329   }
0330 
0331   bool HaloPatternFoundInHB = false;
0332   for (auto& hcand : hccandHB) {
0333     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0334       HaloPatternFoundInHB = true;
0335       edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0336       AddtoBeamHaloHBHERechits(bhrhcandidates, TheGlobalHaloData);
0337     }
0338   }
0339 
0340   bool HaloPatternFoundInHE = false;
0341   for (auto& hcand : hccandHE) {
0342     if ((hcand.getIsHaloFromPattern() && !ishlt) || (hcand.getIsHaloFromPattern_HLT() && ishlt)) {
0343       HaloPatternFoundInHE = true;
0344       edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0345       AddtoBeamHaloHBHERechits(bhrhcandidates, TheGlobalHaloData);
0346     }
0347   }
0348   TheGlobalHaloData.SetHaloPatternFoundEB(HaloPatternFoundInEB);
0349   TheGlobalHaloData.SetHaloPatternFoundEE(HaloPatternFoundInEE);
0350   TheGlobalHaloData.SetHaloPatternFoundHB(HaloPatternFoundInHB);
0351   TheGlobalHaloData.SetHaloPatternFoundHE(HaloPatternFoundInHE);
0352 
0353   return TheGlobalHaloData;
0354 }
0355 
0356 bool GlobalHaloAlgo::SegmentMatchingEB(reco::GlobalHaloData& thehalodata,
0357                                        const std::vector<HaloClusterCandidateECAL>& haloclustercands,
0358                                        float iZ,
0359                                        float iR,
0360                                        float iT,
0361                                        float iPhi,
0362                                        bool ishlt) {
0363   bool rhmatchingfound = false;
0364 
0365   for (auto& hcand : haloclustercands) {
0366     if (!ApplyMatchingCuts(EB,
0367                            ishlt,
0368                            hcand.getSeedEt(),
0369                            iZ,
0370                            hcand.getSeedZ(),
0371                            iR,
0372                            hcand.getSeedR(),
0373                            iT,
0374                            hcand.getSeedTime(),
0375                            iPhi,
0376                            hcand.getSeedPhi()))
0377       continue;
0378 
0379     rhmatchingfound = true;
0380 
0381     edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0382 
0383     AddtoBeamHaloEBEERechits(bhrhcandidates, thehalodata, true);
0384   }
0385 
0386   return rhmatchingfound;
0387 }
0388 
0389 bool GlobalHaloAlgo::SegmentMatchingEE(reco::GlobalHaloData& thehalodata,
0390                                        const std::vector<HaloClusterCandidateECAL>& haloclustercands,
0391                                        float iZ,
0392                                        float iR,
0393                                        float iT,
0394                                        float iPhi,
0395                                        bool ishlt) {
0396   bool rhmatchingfound = false;
0397 
0398   for (auto& hcand : haloclustercands) {
0399     if (!ApplyMatchingCuts(EE,
0400                            ishlt,
0401                            hcand.getSeedEt(),
0402                            iZ,
0403                            hcand.getSeedZ(),
0404                            iR,
0405                            hcand.getSeedR(),
0406                            iT,
0407                            hcand.getSeedTime(),
0408                            iPhi,
0409                            hcand.getSeedPhi()))
0410       continue;
0411 
0412     rhmatchingfound = true;
0413 
0414     edm::RefVector<EcalRecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0415 
0416     AddtoBeamHaloEBEERechits(bhrhcandidates, thehalodata, false);
0417   }
0418 
0419   return rhmatchingfound;
0420 }
0421 
0422 bool GlobalHaloAlgo::SegmentMatchingHB(reco::GlobalHaloData& thehalodata,
0423                                        const std::vector<HaloClusterCandidateHCAL>& haloclustercands,
0424                                        float iZ,
0425                                        float iR,
0426                                        float iT,
0427                                        float iPhi,
0428                                        bool ishlt) {
0429   bool rhmatchingfound = false;
0430 
0431   for (auto& hcand : haloclustercands) {
0432     if (!ApplyMatchingCuts(HB,
0433                            ishlt,
0434                            hcand.getSeedEt(),
0435                            iZ,
0436                            hcand.getSeedZ(),
0437                            iR,
0438                            hcand.getSeedR(),
0439                            iT,
0440                            hcand.getSeedTime(),
0441                            iPhi,
0442                            hcand.getSeedPhi()))
0443       continue;
0444 
0445     rhmatchingfound = true;
0446 
0447     edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0448 
0449     AddtoBeamHaloHBHERechits(bhrhcandidates, thehalodata);
0450   }
0451 
0452   return rhmatchingfound;
0453 }
0454 
0455 bool GlobalHaloAlgo::SegmentMatchingHE(reco::GlobalHaloData& thehalodata,
0456                                        const std::vector<HaloClusterCandidateHCAL>& haloclustercands,
0457                                        float iZ,
0458                                        float iR,
0459                                        float iT,
0460                                        float iPhi,
0461                                        bool ishlt) {
0462   bool rhmatchingfound = false;
0463 
0464   for (auto& hcand : haloclustercands) {
0465     if (!ApplyMatchingCuts(HE,
0466                            ishlt,
0467                            hcand.getSeedEt(),
0468                            iZ,
0469                            hcand.getSeedZ(),
0470                            iR,
0471                            hcand.getSeedR(),
0472                            iT,
0473                            hcand.getSeedTime(),
0474                            iPhi,
0475                            hcand.getSeedPhi()))
0476       continue;
0477 
0478     rhmatchingfound = true;
0479 
0480     edm::RefVector<HBHERecHitCollection> bhrhcandidates = hcand.getBeamHaloRecHitsCandidates();
0481 
0482     AddtoBeamHaloHBHERechits(bhrhcandidates, thehalodata);
0483   }
0484 
0485   return rhmatchingfound;
0486 }
0487 
0488 bool GlobalHaloAlgo::ApplyMatchingCuts(int subdet,
0489                                        bool ishlt,
0490                                        double rhet,
0491                                        double segZ,
0492                                        double rhZ,
0493                                        double segR,
0494                                        double rhR,
0495                                        double segT,
0496                                        double rhT,
0497                                        double segPhi,
0498                                        double rhPhi) {
0499   //Std::Absolute time wrt BX
0500   double tBXrh = rhT + sqrt(rhR * rhR + rhZ * rhZ) / c_cm_per_ns;
0501   double tBXseg = segT + sqrt(segR * segR + segZ * segZ) / c_cm_per_ns;
0502   //Time at z=0, under beam halo hypothesis
0503   double tcorseg = tBXseg - std::abs(segZ) / c_cm_per_ns;       //Outgoing beam halo
0504   double tcorsegincbh = tBXseg + std::abs(segZ) / c_cm_per_ns;  //Ingoing beam halo
0505   double truedt[4] = {1000, 1000, 1000, 1000};
0506   //There are four types of segments associated to beam halo, test each hypothesis:
0507   //IT beam halo, ingoing track
0508   double twindow_seg = 15;
0509   if (std::abs(tcorsegincbh) < twindow_seg)
0510     truedt[0] = tBXrh - tBXseg - std::abs(rhZ - segZ) / c_cm_per_ns;
0511   //IT beam halo, outgoing track
0512   if (std::abs(tcorseg) < twindow_seg)
0513     truedt[1] = tBXseg - tBXrh - std::abs(rhZ - segZ) / c_cm_per_ns;
0514   //OT beam halo (from next BX), ingoing track
0515   if (tcorsegincbh > 25 - twindow_seg && std::abs(tcorsegincbh) < 25 + twindow_seg)
0516     truedt[2] = tBXrh - tBXseg - std::abs(rhZ - segZ) / c_cm_per_ns;
0517   //OT beam halo (from next BX), outgoing track
0518   if (tcorseg > 25 - twindow_seg && tcorseg < 25 + twindow_seg)
0519     truedt[3] = tBXseg - tBXrh - std::abs(rhZ - segZ) / c_cm_per_ns;
0520 
0521   if (subdet == EB) {
0522     if (rhet < et_thresh_rh_eb)
0523       return false;
0524     if (rhet < 20 && ishlt)
0525       return false;
0526     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_eb)
0527       return false;
0528     if (rhR - segR < dr_lowthresh_segvsrh_eb)
0529       return false;
0530     if (rhR - segR > dr_highthresh_segvsrh_eb)
0531       return false;
0532     if (std::abs(truedt[0]) > dt_segvsrh_eb && std::abs(truedt[1]) > dt_segvsrh_eb &&
0533         std::abs(truedt[2]) > dt_segvsrh_eb && std::abs(truedt[3]) > dt_segvsrh_eb)
0534       return false;
0535     return true;
0536   }
0537 
0538   if (subdet == EE) {
0539     if (rhet < et_thresh_rh_ee)
0540       return false;
0541     if (rhet < 20 && ishlt)
0542       return false;
0543     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_ee)
0544       return false;
0545     if (rhR - segR < dr_lowthresh_segvsrh_ee)
0546       return false;
0547     if (rhR - segR > dr_highthresh_segvsrh_ee)
0548       return false;
0549     if (std::abs(truedt[0]) > dt_segvsrh_ee && std::abs(truedt[1]) > dt_segvsrh_ee &&
0550         std::abs(truedt[2]) > dt_segvsrh_ee && std::abs(truedt[3]) > dt_segvsrh_ee)
0551       return false;
0552     return true;
0553   }
0554 
0555   if (subdet == HB) {
0556     if (rhet < et_thresh_rh_hb)
0557       return false;
0558     if (rhet < 20 && ishlt)
0559       return false;
0560     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_hb)
0561       return false;
0562     if (rhR - segR < dr_lowthresh_segvsrh_hb)
0563       return false;
0564     if (rhR - segR > dr_highthresh_segvsrh_hb)
0565       return false;
0566     if (std::abs(truedt[0]) > dt_segvsrh_hb && std::abs(truedt[1]) > dt_segvsrh_hb &&
0567         std::abs(truedt[2]) > dt_segvsrh_hb && std::abs(truedt[3]) > dt_segvsrh_hb)
0568       return false;
0569     return true;
0570   }
0571 
0572   if (subdet == HE) {
0573     if (rhet < et_thresh_rh_he)
0574       return false;
0575     if (rhet < 20 && ishlt)
0576       return false;
0577     if (std::abs(deltaPhi(rhPhi, segPhi)) > dphi_thresh_segvsrh_he)
0578       return false;
0579     if (rhR - segR < dr_lowthresh_segvsrh_he)
0580       return false;
0581     if (rhR - segR > dr_highthresh_segvsrh_he)
0582       return false;
0583     if (std::abs(truedt[0]) > dt_segvsrh_he && std::abs(truedt[1]) > dt_segvsrh_he &&
0584         std::abs(truedt[2]) > dt_segvsrh_he && std::abs(truedt[3]) > dt_segvsrh_he)
0585       return false;
0586     return true;
0587   }
0588 
0589   return false;
0590 }
0591 
0592 void GlobalHaloAlgo::AddtoBeamHaloEBEERechits(edm::RefVector<EcalRecHitCollection>& bhtaggedrechits,
0593                                               reco::GlobalHaloData& thehalodata,
0594                                               bool isbarrel) {
0595   for (size_t ihit = 0; ihit < bhtaggedrechits.size(); ++ihit) {
0596     bool alreadyincl = false;
0597     edm::Ref<EcalRecHitCollection> rhRef(bhtaggedrechits[ihit]);
0598     edm::RefVector<EcalRecHitCollection> refrhcoll;
0599     if (isbarrel)
0600       refrhcoll = thehalodata.GetEBRechits();
0601     else
0602       refrhcoll = thehalodata.GetEERechits();
0603     for (size_t jhit = 0; jhit < refrhcoll.size(); jhit++) {
0604       edm::Ref<EcalRecHitCollection> rhitRef(refrhcoll[jhit]);
0605       if (rhitRef->detid() == rhRef->detid())
0606         alreadyincl = true;
0607       if (rhitRef->detid() == rhRef->detid())
0608         break;
0609     }
0610     if (!alreadyincl && isbarrel)
0611       thehalodata.GetEBRechits().push_back(rhRef);
0612     if (!alreadyincl && !isbarrel)
0613       thehalodata.GetEERechits().push_back(rhRef);
0614   }
0615 }
0616 
0617 void GlobalHaloAlgo::AddtoBeamHaloHBHERechits(edm::RefVector<HBHERecHitCollection>& bhtaggedrechits,
0618                                               reco::GlobalHaloData& thehalodata) {
0619   for (size_t ihit = 0; ihit < bhtaggedrechits.size(); ++ihit) {
0620     bool alreadyincl = false;
0621     edm::Ref<HBHERecHitCollection> rhRef(bhtaggedrechits[ihit]);
0622     edm::RefVector<HBHERecHitCollection> refrhcoll;
0623     refrhcoll = thehalodata.GetHBHERechits();
0624     for (size_t jhit = 0; jhit < refrhcoll.size(); jhit++) {
0625       edm::Ref<HBHERecHitCollection> rhitRef(refrhcoll[jhit]);
0626       if (rhitRef->detid() == rhRef->detid())
0627         alreadyincl = true;
0628       if (rhitRef->detid() == rhRef->detid())
0629         break;
0630     }
0631     if (!alreadyincl)
0632       thehalodata.GetHBHERechits().push_back(rhRef);
0633   }
0634 }
[ 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!