Back to home page

Project CMSSW displayed by LXR
 
 

    

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