Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 
0002 #include "RecoMET/METAlgorithms/interface/CSCHaloAlgo.h"
0003 #include "FWCore/Common/interface/TriggerNames.h"
0004 /*
0005   [class]:  CSCHaloAlgo
0006   [authors]: R. Remington, The University of Florida
0007   [description]: See CSCHaloAlgo.h
0008   [date]: October 15, 2009
0009 */
0010 using namespace reco;
0011 using namespace std;
0012 using namespace edm;
0013 #include "TMath.h"
0014 namespace {
0015   constexpr float c_cm_per_ns = 29.9792458;
0016 };
0017 CSCHaloAlgo::CSCHaloAlgo(edm::ConsumesCollector iC) : geoToken_(iC.esConsumes()), geo_(nullptr), hgeo_(nullptr) {
0018   deta_threshold = 0.;
0019   min_inner_radius = 0.;
0020   max_inner_radius = 9999.;
0021   min_outer_radius = 0.;
0022   max_outer_radius = 9999.;
0023   dphi_threshold = 999.;
0024   norm_chi2_threshold = 999.;
0025   recHit_t0 = 0.;
0026   recHit_twindow = 25.;
0027   expected_BX = 3;
0028   max_dt_muon_segment = -10.0;
0029   max_free_inverse_beta = 0.0;
0030 
0031   min_outer_theta = 0.;
0032   max_outer_theta = TMath::Pi();
0033 
0034   matching_dphi_threshold = 0.18;  //radians
0035   matching_deta_threshold = 0.4;
0036   matching_dwire_threshold = 5.;
0037 
0038   et_thresh_rh_hbhe = 25;  //GeV
0039   et_thresh_rh_ee = 10;
0040   et_thresh_rh_eb = 10;
0041 
0042   dphi_thresh_segvsrh_hbhe = 0.05;  //radians
0043   dphi_thresh_segvsrh_eb = 0.05;
0044   dphi_thresh_segvsrh_ee = 0.05;
0045 
0046   dr_lowthresh_segvsrh_hbhe = -20;  //cm
0047   dr_lowthresh_segvsrh_eb = -20;
0048   dr_lowthresh_segvsrh_ee = -20;
0049 
0050   dr_highthresh_segvsrh_hbhe = 20;  //cm
0051   dr_highthresh_segvsrh_eb = 20;
0052   dr_highthresh_segvsrh_ee = 20;
0053 
0054   dt_lowthresh_segvsrh_hbhe = 5;  //ns
0055   dt_lowthresh_segvsrh_eb = 5;
0056   dt_lowthresh_segvsrh_ee = 5;
0057 
0058   dt_highthresh_segvsrh_hbhe = 30;  //ns
0059   dt_highthresh_segvsrh_eb = 30;
0060   dt_highthresh_segvsrh_ee = 30;
0061 }
0062 
0063 reco::CSCHaloData CSCHaloAlgo::Calculate(const CSCGeometry& TheCSCGeometry,
0064                                          edm::Handle<reco::MuonCollection>& TheCosmicMuons,
0065                                          const edm::Handle<reco::MuonTimeExtraMap> TheCSCTimeMap,
0066                                          edm::Handle<reco::MuonCollection>& TheMuons,
0067                                          edm::Handle<CSCSegmentCollection>& TheCSCSegments,
0068                                          edm::Handle<CSCRecHit2DCollection>& TheCSCRecHits,
0069                                          edm::Handle<L1MuGMTReadoutCollection>& TheL1GMTReadout,
0070                                          edm::Handle<HBHERecHitCollection>& hbhehits,
0071                                          edm::Handle<EcalRecHitCollection>& ecalebhits,
0072                                          edm::Handle<EcalRecHitCollection>& ecaleehits,
0073                                          edm::Handle<edm::TriggerResults>& TheHLTResults,
0074                                          const edm::TriggerNames* triggerNames,
0075                                          const edm::Handle<CSCALCTDigiCollection>& TheALCTs,
0076                                          MuonSegmentMatcher* TheMatcher,
0077                                          const edm::Event& TheEvent,
0078                                          const edm::EventSetup& TheSetup) {
0079   reco::CSCHaloData TheCSCHaloData;
0080   int imucount = 0;
0081 
0082   bool calomatched = false;
0083   bool ECALBmatched = false;
0084   bool ECALEmatched = false;
0085   bool HCALmatched = false;
0086 
0087   geo_ = &TheSetup.getData(geoToken_);
0088   hgeo_ = dynamic_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, 1));
0089 
0090   //}
0091   bool trkmuunvetoisdefault = false;  //Pb with low pt tracker muons that veto good csc segments/halo triggers.
0092   //Test to "unveto" low pt trk muons.
0093   //For now, we just recalculate everything without the veto and add an extra set of variables to the class CSCHaloData.
0094   //If this is satisfactory, these variables can become the default ones by setting trkmuunvetoisdefault to true.
0095   if (TheCosmicMuons.isValid()) {
0096     short int n_tracks_small_beta = 0;
0097     short int n_tracks_small_dT = 0;
0098     short int n_tracks_small_dT_and_beta = 0;
0099     for (reco::MuonCollection::const_iterator iMuon = TheCosmicMuons->begin(); iMuon != TheCosmicMuons->end();
0100          iMuon++, imucount++) {
0101       reco::TrackRef Track = iMuon->outerTrack();
0102       if (!Track)
0103         continue;
0104 
0105       bool StoreTrack = false;
0106       // Calculate global phi coordinate for central most rechit in the track
0107       float innermost_global_z = 1500.;
0108       float outermost_global_z = 0.;
0109       GlobalPoint InnerMostGlobalPosition(0., 0., 0.);  // smallest abs(z)
0110       GlobalPoint OuterMostGlobalPosition(0., 0., 0.);  // largest abs(z)
0111       int nCSCHits = 0;
0112       for (unsigned int j = 0; j < Track->extra()->recHitsSize(); j++) {
0113         auto hit = Track->extra()->recHitRef(j);
0114         if (!hit->isValid())
0115           continue;
0116         DetId TheDetUnitId(hit->geographicalId());
0117         if (TheDetUnitId.det() != DetId::Muon)
0118           continue;
0119         if (TheDetUnitId.subdetId() != MuonSubdetId::CSC)
0120           continue;
0121 
0122         const GeomDetUnit* TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
0123         LocalPoint TheLocalPosition = hit->localPosition();
0124         const BoundPlane& TheSurface = TheUnit->surface();
0125         const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);
0126 
0127         float z = TheGlobalPosition.z();
0128         if (abs(z) < innermost_global_z) {
0129           innermost_global_z = abs(z);
0130           InnerMostGlobalPosition = GlobalPoint(TheGlobalPosition);
0131         }
0132         if (abs(z) > outermost_global_z) {
0133           outermost_global_z = abs(z);
0134           OuterMostGlobalPosition = GlobalPoint(TheGlobalPosition);
0135         }
0136         nCSCHits++;
0137       }
0138 
0139       std::vector<const CSCSegment*> MatchedSegments = TheMatcher->matchCSC(*Track, TheEvent);
0140       // Find the inner and outer segments separately in case they don't agree completely with recHits
0141       // Plan for the possibility segments in both endcaps
0142       float InnerSegmentTime[2] = {0, 0};
0143       float OuterSegmentTime[2] = {0, 0};
0144       float innermost_seg_z[2] = {1500, 1500};
0145       float outermost_seg_z[2] = {0, 0};
0146       for (std::vector<const CSCSegment*>::const_iterator segment = MatchedSegments.begin();
0147            segment != MatchedSegments.end();
0148            ++segment) {
0149         CSCDetId TheCSCDetId((*segment)->cscDetId());
0150         const CSCChamber* TheCSCChamber = TheCSCGeometry.chamber(TheCSCDetId);
0151         LocalPoint TheLocalPosition = (*segment)->localPosition();
0152         const GlobalPoint TheGlobalPosition = TheCSCChamber->toGlobal(TheLocalPosition);
0153         float z = TheGlobalPosition.z();
0154         int TheEndcap = TheCSCDetId.endcap();
0155         if (abs(z) < innermost_seg_z[TheEndcap - 1]) {
0156           innermost_seg_z[TheEndcap - 1] = abs(z);
0157           InnerSegmentTime[TheEndcap - 1] = (*segment)->time();
0158         }
0159         if (abs(z) > outermost_seg_z[TheEndcap - 1]) {
0160           outermost_seg_z[TheEndcap - 1] = abs(z);
0161           OuterSegmentTime[TheEndcap - 1] = (*segment)->time();
0162         }
0163       }
0164 
0165       if (nCSCHits < 3)
0166         continue;  // This needs to be optimized, but is the minimum
0167 
0168       float dT_Segment = 0;  // default safe value, looks like collision muon
0169 
0170       if (innermost_seg_z[0] < outermost_seg_z[0])  // two segments in ME+
0171         dT_Segment = OuterSegmentTime[0] - InnerSegmentTime[0];
0172       if (innermost_seg_z[1] < outermost_seg_z[1])  // two segments in ME-
0173       {
0174         // replace the measurement if there weren't segments in ME+ or
0175         // if the track in ME- has timing more consistent with an incoming particle
0176         if (dT_Segment == 0.0 || OuterSegmentTime[1] - InnerSegmentTime[1] < dT_Segment)
0177           dT_Segment = OuterSegmentTime[1] - InnerSegmentTime[1];
0178       }
0179 
0180       if (OuterMostGlobalPosition.x() == 0. || OuterMostGlobalPosition.y() == 0. || OuterMostGlobalPosition.z() == 0.)
0181         continue;
0182       if (InnerMostGlobalPosition.x() == 0. || InnerMostGlobalPosition.y() == 0. || InnerMostGlobalPosition.z() == 0.)
0183         continue;
0184 
0185       //Its a CSC Track,store it if it passes halo selection
0186       StoreTrack = true;
0187 
0188       float deta = abs(OuterMostGlobalPosition.eta() - InnerMostGlobalPosition.eta());
0189       float dphi = abs(deltaPhi(OuterMostGlobalPosition.barePhi(), InnerMostGlobalPosition.barePhi()));
0190       float theta = Track->outerMomentum().theta();
0191       float innermost_x = InnerMostGlobalPosition.x();
0192       float innermost_y = InnerMostGlobalPosition.y();
0193       float outermost_x = OuterMostGlobalPosition.x();
0194       float outermost_y = OuterMostGlobalPosition.y();
0195       float innermost_r = TMath::Sqrt(innermost_x * innermost_x + innermost_y * innermost_y);
0196       float outermost_r = TMath::Sqrt(outermost_x * outermost_x + outermost_y * outermost_y);
0197 
0198       if (deta < deta_threshold)
0199         StoreTrack = false;
0200       if (theta > min_outer_theta && theta < max_outer_theta)
0201         StoreTrack = false;
0202       if (dphi > dphi_threshold)
0203         StoreTrack = false;
0204       if (innermost_r < min_inner_radius)
0205         StoreTrack = false;
0206       if (innermost_r > max_inner_radius)
0207         StoreTrack = false;
0208       if (outermost_r < min_outer_radius)
0209         StoreTrack = false;
0210       if (outermost_r > max_outer_radius)
0211         StoreTrack = false;
0212       if (Track->normalizedChi2() > norm_chi2_threshold)
0213         StoreTrack = false;
0214 
0215       if (StoreTrack) {
0216         TheCSCHaloData.GetCSCTrackImpactPositions().push_back(InnerMostGlobalPosition);
0217         TheCSCHaloData.GetTracks().push_back(Track);
0218       }
0219 
0220       // Analyze the MuonTimeExtra information
0221       if (TheCSCTimeMap.isValid()) {
0222         reco::MuonRef muonR(TheCosmicMuons, imucount);
0223         const reco::MuonTimeExtraMap& timeMapCSC = *TheCSCTimeMap;
0224         reco::MuonTimeExtra timecsc = timeMapCSC[muonR];
0225         float freeInverseBeta = timecsc.freeInverseBeta();
0226 
0227         if (dT_Segment < max_dt_muon_segment)
0228           n_tracks_small_dT++;
0229         if (freeInverseBeta < max_free_inverse_beta)
0230           n_tracks_small_beta++;
0231         if ((dT_Segment < max_dt_muon_segment) && (freeInverseBeta < max_free_inverse_beta))
0232           n_tracks_small_dT_and_beta++;
0233       } else {
0234         static std::atomic<bool> MuonTimeFail{false};
0235         bool expected = false;
0236         if (MuonTimeFail.compare_exchange_strong(expected, true, std::memory_order_acq_rel)) {
0237           edm::LogWarning("InvalidInputTag")
0238               << "The MuonTimeExtraMap does not appear to be in the event. Some beam halo "
0239               << " identification variables will be empty";
0240         }
0241       }
0242     }
0243     TheCSCHaloData.SetNIncomingTracks(n_tracks_small_dT, n_tracks_small_beta, n_tracks_small_dT_and_beta);
0244   } else  // collection is invalid
0245   {
0246     static std::atomic<bool> CosmicFail{false};
0247     bool expected = false;
0248     if (CosmicFail.compare_exchange_strong(expected, true, std::memory_order_acq_rel)) {
0249       edm::LogWarning("InvalidInputTag")
0250           << " The Cosmic Muon collection does not appear to be in the event. These beam halo "
0251           << " identification variables will be empty";
0252     }
0253   }
0254 
0255   // Loop over the CSCRecHit2D collection to look for out-of-time recHits
0256   // Out-of-time is defined as tpeak outside [t_0 + TOF - t_window, t_0 + TOF + t_window]
0257   // where t_0 and t_window are configurable parameters
0258   short int n_recHitsP = 0;
0259   short int n_recHitsM = 0;
0260   if (TheCSCRecHits.isValid()) {
0261     CSCRecHit2DCollection::const_iterator dRHIter;
0262     for (dRHIter = TheCSCRecHits->begin(); dRHIter != TheCSCRecHits->end(); dRHIter++) {
0263       if (!((*dRHIter).isValid()))
0264         continue;  // only interested in valid hits
0265       CSCDetId idrec = (CSCDetId)(*dRHIter).cscDetId();
0266       float RHTime = (*dRHIter).tpeak();
0267       LocalPoint rhitlocal = (*dRHIter).localPosition();
0268       const CSCChamber* chamber = TheCSCGeometry.chamber(idrec);
0269       GlobalPoint globalPosition = chamber->toGlobal(rhitlocal);
0270       float globZ = globalPosition.z();
0271       if (RHTime < (recHit_t0 - recHit_twindow)) {
0272         if (globZ > 0)
0273           n_recHitsP++;
0274         else
0275           n_recHitsM++;
0276       }
0277 
0278       /*
0279 
0280        float globX = globalPosition.x();
0281        float globY = globalPosition.y();
0282        float globZ = globalPosition.z();
0283        float TOF = (sqrt(globX*globX+ globY*globY + globZ*globZ))/29.9792458 ; //cm -> ns
0284        if ( (RHTime < (recHit_t0 + TOF - recHit_twindow)) || (RHTime > (recHit_t0 + TOF + recHit_twindow)) )
0285          {
0286            if( globZ > 0 ) 
0287          n_recHitsP++;
0288            else
0289          n_recHitsM++;
0290          }
0291        */
0292     }
0293   } else {
0294     static std::atomic<bool> RecHitFail{false};
0295     bool expected = false;
0296     if (RecHitFail.compare_exchange_strong(expected, true, std::memory_order_acq_rel)) {
0297       edm::LogWarning("InvalidInputTag")
0298           << "The requested CSCRecHit2DCollection does not appear to be in the event. The CSC RecHit "
0299           << " variables used for halo identification will not be calculated or stored";
0300     }
0301   }
0302   TheCSCHaloData.SetNOutOfTimeHits(n_recHitsP + n_recHitsM);
0303   // MLR
0304   // Loop through CSCSegments and count the number of "flat" segments with the same (r,phi),
0305   // saving the value in TheCSCHaloData.
0306   short int maxNSegments = 0;
0307   bool plus_endcap = false;
0308   bool minus_endcap = false;
0309   bool both_endcaps = false;
0310   bool both_endcaps_loose = false;
0311   //   bool both_endcaps_dtcut = false;
0312 
0313   short int maxNSegments_alt = 0;
0314   bool both_endcaps_alt = false;
0315   bool both_endcaps_loose_alt = false;
0316   bool both_endcaps_loose_dtcut_alt = false;
0317 
0318   //float r = 0., phi = 0.;
0319   if (TheCSCSegments.isValid()) {
0320     for (CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end();
0321          iSegment++) {
0322       CSCDetId iCscDetID = iSegment->cscDetId();
0323       bool Segment1IsGood = true;
0324       bool Segment1IsGood_alt = true;
0325 
0326       //avoid segments from collision muons
0327       if (TheMuons.isValid()) {
0328         for (reco::MuonCollection::const_iterator mu = TheMuons->begin();
0329              mu != TheMuons->end() && (Segment1IsGood || !trkmuunvetoisdefault);
0330              mu++) {
0331           bool lowpttrackmu = false;
0332           if (!mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon())
0333             continue;
0334           if (!mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon() && trkmuunvetoisdefault)
0335             continue;
0336           if (!mu->isGlobalMuon() && mu->isTrackerMuon() && mu->pt() < 3)
0337             lowpttrackmu = true;
0338           const std::vector<MuonChamberMatch> chambers = mu->matches();
0339           for (std::vector<MuonChamberMatch>::const_iterator kChamber = chambers.begin(); kChamber != chambers.end();
0340                kChamber++) {
0341             if (kChamber->detector() != MuonSubdetId::CSC)
0342               continue;
0343             for (std::vector<reco::MuonSegmentMatch>::const_iterator kSegment = kChamber->segmentMatches.begin();
0344                  kSegment != kChamber->segmentMatches.end();
0345                  kSegment++) {
0346               edm::Ref<CSCSegmentCollection> cscSegRef = kSegment->cscSegmentRef;
0347               CSCDetId kCscDetID = cscSegRef->cscDetId();
0348 
0349               if (kCscDetID == iCscDetID) {
0350                 Segment1IsGood = false;
0351                 if (!lowpttrackmu)
0352                   Segment1IsGood_alt = false;
0353               }
0354             }
0355           }
0356         }
0357       }
0358       if (!Segment1IsGood && !Segment1IsGood_alt)
0359         continue;
0360 
0361       // Get local direction vector; if direction runs parallel to beamline,
0362       // count this segment as beam halo candidate.
0363       LocalPoint iLocalPosition = iSegment->localPosition();
0364       LocalVector iLocalDirection = iSegment->localDirection();
0365 
0366       GlobalPoint iGlobalPosition = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalPosition);
0367       GlobalVector iGlobalDirection = TheCSCGeometry.chamber(iCscDetID)->toGlobal(iLocalDirection);
0368 
0369       float iTheta = iGlobalDirection.theta();
0370       if (iTheta > max_segment_theta && iTheta < TMath::Pi() - max_segment_theta)
0371         continue;
0372 
0373       float iPhi = iGlobalPosition.phi();
0374       float iR = TMath::Sqrt(iGlobalPosition.x() * iGlobalPosition.x() + iGlobalPosition.y() * iGlobalPosition.y());
0375       float iZ = iGlobalPosition.z();
0376       float iT = iSegment->time();
0377 
0378       //       if(abs(iZ)<650&& TheEvent.id().run()< 251737) iT-= 25;
0379       //Calo matching:
0380 
0381       bool hbhematched = HCALSegmentMatching(hbhehits,
0382                                              et_thresh_rh_hbhe,
0383                                              dphi_thresh_segvsrh_hbhe,
0384                                              dr_lowthresh_segvsrh_hbhe,
0385                                              dr_highthresh_segvsrh_hbhe,
0386                                              dt_lowthresh_segvsrh_hbhe,
0387                                              dt_highthresh_segvsrh_hbhe,
0388                                              iZ,
0389                                              iR,
0390                                              iT,
0391                                              iPhi);
0392       bool ebmatched = ECALSegmentMatching(ecalebhits,
0393                                            et_thresh_rh_eb,
0394                                            dphi_thresh_segvsrh_eb,
0395                                            dr_lowthresh_segvsrh_eb,
0396                                            dr_highthresh_segvsrh_eb,
0397                                            dt_lowthresh_segvsrh_eb,
0398                                            dt_highthresh_segvsrh_eb,
0399                                            iZ,
0400                                            iR,
0401                                            iT,
0402                                            iPhi);
0403       bool eematched = ECALSegmentMatching(ecaleehits,
0404                                            et_thresh_rh_ee,
0405                                            dphi_thresh_segvsrh_ee,
0406                                            dr_lowthresh_segvsrh_ee,
0407                                            dr_highthresh_segvsrh_ee,
0408                                            dt_lowthresh_segvsrh_ee,
0409                                            dt_highthresh_segvsrh_ee,
0410                                            iZ,
0411                                            iR,
0412                                            iT,
0413                                            iPhi);
0414       calomatched |= (hbhematched || ebmatched || eematched);
0415       HCALmatched |= hbhematched;
0416       ECALBmatched |= ebmatched;
0417       ECALEmatched |= eematched;
0418 
0419       short int nSegs = 0;
0420       short int nSegs_alt = 0;
0421       // Changed to loop over all Segments (so N^2) to catch as many segments as possible.
0422       for (CSCSegmentCollection::const_iterator jSegment = TheCSCSegments->begin(); jSegment != TheCSCSegments->end();
0423            jSegment++) {
0424         if (jSegment == iSegment)
0425           continue;
0426         bool Segment2IsGood = true;
0427         bool Segment2IsGood_alt = true;
0428         LocalPoint jLocalPosition = jSegment->localPosition();
0429         LocalVector jLocalDirection = jSegment->localDirection();
0430         CSCDetId jCscDetID = jSegment->cscDetId();
0431         GlobalPoint jGlobalPosition = TheCSCGeometry.chamber(jCscDetID)->toGlobal(jLocalPosition);
0432         GlobalVector jGlobalDirection = TheCSCGeometry.chamber(jCscDetID)->toGlobal(jLocalDirection);
0433         float jTheta = jGlobalDirection.theta();
0434         float jPhi = jGlobalPosition.phi();
0435         float jR = TMath::Sqrt(jGlobalPosition.x() * jGlobalPosition.x() + jGlobalPosition.y() * jGlobalPosition.y());
0436         float jZ = jGlobalPosition.z();
0437         float jT = jSegment->time();
0438         //   if(abs(jZ)<650&& TheEvent.id().run() < 251737)jT-= 25;
0439         if (std::abs(deltaPhi(jPhi, iPhi)) <= max_segment_phi_diff &&
0440             (std::abs(jR - iR) < 0.05 * std::abs(jZ - iZ) || jZ * iZ > 0) &&
0441             ((jR - iR) > -0.02 * std::abs(jZ - iZ) || iT > jT || jZ * iZ > 0) &&
0442             ((iR - jR) > -0.02 * std::abs(jZ - iZ) || iT < jT || jZ * iZ > 0) &&
0443             (std::abs(jR - iR) <= max_segment_r_diff || jZ * iZ < 0) &&
0444             (jTheta < max_segment_theta || jTheta > TMath::Pi() - max_segment_theta)) {
0445           //// Check if Segment matches to a colision muon
0446           if (TheMuons.isValid()) {
0447             for (reco::MuonCollection::const_iterator mu = TheMuons->begin();
0448                  mu != TheMuons->end() && (Segment2IsGood || !trkmuunvetoisdefault);
0449                  mu++) {
0450               bool lowpttrackmu = false;
0451               if (!mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon())
0452                 continue;
0453               if (!mu->isGlobalMuon() && mu->isTrackerMuon() && mu->pt() < 3)
0454                 lowpttrackmu = true;
0455               const std::vector<MuonChamberMatch> chambers = mu->matches();
0456               for (std::vector<MuonChamberMatch>::const_iterator kChamber = chambers.begin();
0457                    kChamber != chambers.end();
0458                    kChamber++) {
0459                 if (kChamber->detector() != MuonSubdetId::CSC)
0460                   continue;
0461                 for (std::vector<reco::MuonSegmentMatch>::const_iterator kSegment = kChamber->segmentMatches.begin();
0462                      kSegment != kChamber->segmentMatches.end();
0463                      kSegment++) {
0464                   edm::Ref<CSCSegmentCollection> cscSegRef = kSegment->cscSegmentRef;
0465                   CSCDetId kCscDetID = cscSegRef->cscDetId();
0466 
0467                   if (kCscDetID == jCscDetID) {
0468                     Segment2IsGood = false;
0469                     if (!lowpttrackmu)
0470                       Segment2IsGood_alt = false;
0471                   }
0472                 }
0473               }
0474             }
0475           }
0476           if (Segment1IsGood && Segment2IsGood) {
0477             nSegs++;
0478             minus_endcap = iGlobalPosition.z() < 0 || jGlobalPosition.z() < 0;
0479             plus_endcap = iGlobalPosition.z() > 0 || jGlobalPosition.z() > 0;
0480             //       if( abs(jT-iT)/sqrt( (jR-iR)*(jR-iR)+(jZ-iZ)*(jZ-iZ) )<0.05 && abs(jT-iT)/sqrt( (jR-iR)*(jR-iR)+(jZ-iZ)*(jZ-iZ) )>0.02 && minus_endcap&&plus_endcap ) both_endcaps_dtcut =true;
0481           }
0482           if (Segment1IsGood_alt && Segment2IsGood_alt) {
0483             nSegs_alt++;
0484             minus_endcap = iGlobalPosition.z() < 0 || jGlobalPosition.z() < 0;
0485             plus_endcap = iGlobalPosition.z() > 0 || jGlobalPosition.z() > 0;
0486             double iTBX = iT + sqrt(iGlobalPosition.mag2()) / c_cm_per_ns;
0487             double jTBX = jT + sqrt(jGlobalPosition.mag2()) / c_cm_per_ns;
0488             double truedt = (iTBX > jTBX) ? iTBX - jTBX - std::abs(iZ - jZ) / c_cm_per_ns
0489                                           : jTBX - iTBX - std::abs(iZ - jZ) / c_cm_per_ns;
0490             if (std::abs(truedt) < 15 && (iT > -15 || jT > -15) && minus_endcap && plus_endcap)
0491               both_endcaps_loose_dtcut_alt = true;
0492           }
0493         }
0494       }
0495       // Correct the fact that the way nSegs counts will always be short by 1
0496       if (nSegs > 0)
0497         nSegs++;
0498 
0499       // The opposite endcaps segments do not need to belong to the longest chain.
0500       if (nSegs > 0)
0501         both_endcaps_loose = both_endcaps_loose ? both_endcaps_loose : minus_endcap && plus_endcap;
0502       if (nSegs_alt > 0)
0503         nSegs_alt++;
0504       if (nSegs_alt > 0)
0505         both_endcaps_loose_alt = both_endcaps_loose_alt ? both_endcaps_loose_alt : minus_endcap && plus_endcap;
0506 
0507       //       if (nSegs > 0) both_endcaps_dt20ns = both_endcaps_dt20ns ? both_endcaps_dt20ns : minus_endcap && plus_endcap &&dt20ns;
0508       if (nSegs > maxNSegments) {
0509         // Use value of r, phi to collect halo CSCSegments for examining timing (not coded yet...)
0510         //r = iR;
0511         //phi = iPhi;
0512         maxNSegments = nSegs;
0513         both_endcaps = both_endcaps ? both_endcaps : minus_endcap && plus_endcap;
0514       }
0515 
0516       if (nSegs_alt > maxNSegments_alt) {
0517         maxNSegments_alt = nSegs_alt;
0518         both_endcaps_alt = both_endcaps_alt ? both_endcaps_alt : minus_endcap && plus_endcap;
0519       }
0520     }
0521   }
0522 
0523   //Deprecated methods, kept for backward compatibility
0524   TheCSCHaloData.SetHLTBit(false);
0525   TheCSCHaloData.SetNumberOfHaloTriggers(0, 0);
0526   TheCSCHaloData.SetNumberOfHaloTriggers_TrkMuUnVeto(0, 0);
0527   TheCSCHaloData.SetNOutOfTimeTriggers(0, 0);
0528 
0529   //Current methods used
0530   TheCSCHaloData.SetNFlatHaloSegments(maxNSegments);
0531   TheCSCHaloData.SetSegmentsBothEndcaps(both_endcaps);
0532   TheCSCHaloData.SetNFlatHaloSegments_TrkMuUnVeto(maxNSegments_alt);
0533   TheCSCHaloData.SetSegmentsBothEndcaps_Loose_TrkMuUnVeto(both_endcaps_loose_alt);
0534   TheCSCHaloData.SetSegmentsBothEndcaps_Loose_dTcut_TrkMuUnVeto(both_endcaps_loose_dtcut_alt);
0535   TheCSCHaloData.SetSegmentIsCaloMatched(calomatched);
0536   TheCSCHaloData.SetSegmentIsHCaloMatched(HCALmatched);
0537   TheCSCHaloData.SetSegmentIsEBCaloMatched(ECALBmatched);
0538   TheCSCHaloData.SetSegmentIsEECaloMatched(ECALEmatched);
0539 
0540   return TheCSCHaloData;
0541 }
0542 
0543 math::XYZPoint CSCHaloAlgo::getPosition(const DetId& id, reco::Vertex::Point vtx) {
0544   const GlobalPoint pos = ((id.det() == DetId::Hcal) ? hgeo_->getPosition(id) : GlobalPoint(geo_->getPosition(id)));
0545   math::XYZPoint posV(pos.x() - vtx.x(), pos.y() - vtx.y(), pos.z() - vtx.z());
0546   return posV;
0547 }
0548 
0549 bool CSCHaloAlgo::HCALSegmentMatching(edm::Handle<HBHERecHitCollection>& rechitcoll,
0550                                       float et_thresh_rh,
0551                                       float dphi_thresh_segvsrh,
0552                                       float dr_lowthresh_segvsrh,
0553                                       float dr_highthresh_segvsrh,
0554                                       float dt_lowthresh_segvsrh,
0555                                       float dt_highthresh_segvsrh,
0556                                       float iZ,
0557                                       float iR,
0558                                       float iT,
0559                                       float iPhi) {
0560   reco::Vertex::Point vtx(0, 0, 0);
0561   for (size_t ihit = 0; ihit < rechitcoll->size(); ++ihit) {
0562     const HBHERecHit& rechit = (*rechitcoll)[ihit];
0563     math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
0564     double rhet = rechit.energy() / cosh(rhpos.eta());
0565     double dphi_rhseg = abs(deltaPhi(rhpos.phi(), iPhi));
0566     double dr_rhseg = sqrt(rhpos.x() * rhpos.x() + rhpos.y() * rhpos.y()) - iR;
0567     double dtcorr_rhseg = rechit.time() - abs(rhpos.z() - iZ) / 30 - iT;
0568     if ((rechit.time() < -3) && (rhpos.z() * iZ < 0 || abs(rhpos.z()) < 200) && rhet > et_thresh_rh &&
0569         dphi_rhseg < dphi_thresh_segvsrh && dr_rhseg < dr_highthresh_segvsrh &&
0570         dr_rhseg > dr_lowthresh_segvsrh &&  //careful: asymmetric cut might not be the most appropriate thing
0571         dtcorr_rhseg > dt_lowthresh_segvsrh && dtcorr_rhseg < dt_highthresh_segvsrh)
0572       return true;
0573   }
0574   return false;
0575 }
0576 
0577 bool CSCHaloAlgo::ECALSegmentMatching(edm::Handle<EcalRecHitCollection>& rechitcoll,
0578                                       float et_thresh_rh,
0579                                       float dphi_thresh_segvsrh,
0580                                       float dr_lowthresh_segvsrh,
0581                                       float dr_highthresh_segvsrh,
0582                                       float dt_lowthresh_segvsrh,
0583                                       float dt_highthresh_segvsrh,
0584                                       float iZ,
0585                                       float iR,
0586                                       float iT,
0587                                       float iPhi) {
0588   reco::Vertex::Point vtx(0, 0, 0);
0589   for (size_t ihit = 0; ihit < rechitcoll->size(); ++ihit) {
0590     const EcalRecHit& rechit = (*rechitcoll)[ihit];
0591     math::XYZPoint rhpos = getPosition(rechit.id(), vtx);
0592     double rhet = rechit.energy() / cosh(rhpos.eta());
0593     double dphi_rhseg = abs(deltaPhi(rhpos.phi(), iPhi));
0594     double dr_rhseg = sqrt(rhpos.x() * rhpos.x() + rhpos.y() * rhpos.y()) - iR;
0595     double dtcorr_rhseg = rechit.time() - abs(rhpos.z() - iZ) / 30 - iT;
0596     if ((rechit.time() < -1) && (rhpos.z() * iZ < 0 || abs(rhpos.z()) < 200) && rhet > et_thresh_rh &&
0597         dphi_rhseg < dphi_thresh_segvsrh && dr_rhseg < dr_highthresh_segvsrh &&
0598         dr_rhseg > dr_lowthresh_segvsrh &&  //careful: asymmetric cut might not be the most appropriate thing
0599         dtcorr_rhseg > dt_lowthresh_segvsrh && dtcorr_rhseg < dr_highthresh_segvsrh)
0600       return true;
0601   }
0602   return false;
0603 }
[ 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!