Reference Manual

Go to the documentation of this file.
00001 
00010 #include "DataFormats/Common/interface/Handle.h"
00011 
00012 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
00013 #include "DataFormats/HLTReco/interface/TriggerRefsCollections.h"
00014 
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 
00017 #include "DataFormats/TrackReco/interface/Track.h"
00018 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
00019 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
00020 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00021 #include "HLTrigger/Muon/interface/HLTMuonDimuonL3Filter.h"
00022 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00023 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
00024 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeedCollection.h"
00025 
00026 using namespace edm;
00027 using namespace std;
00028 using namespace reco;
00029 using namespace trigger;
00030 
00031 //
00032 // constructors and destructor
00033 //
00034 HLTMuonDimuonL3Filter::HLTMuonDimuonL3Filter(const edm::ParameterSet& iConfig) :   beamspotTag_   (iConfig.getParameter< edm::InputTag > ("BeamSpotTag")),
00035    candTag_     (iConfig.getParameter< edm::InputTag > ("CandTag")),
00036    previousCandTag_   (iConfig.getParameter<InputTag > ("PreviousCandTag")),
00037     fast_Accept_ (iConfig.getParameter<bool> ("FastAccept")),
00038    max_Eta_     (iConfig.getParameter<double> ("MaxEta")),
00039    min_Nhits_   (iConfig.getParameter<int> ("MinNhits")),
00040    max_Dr_      (iConfig.getParameter<double> ("MaxDr")),
00041    max_Dz_      (iConfig.getParameter<double> ("MaxDz")),
00042    chargeOpt_   (iConfig.getParameter<int> ("ChargeOpt")),
00043    min_PtPair_  (iConfig.getParameter<double> ("MinPtPair")),
00044    min_PtMax_   (iConfig.getParameter<double> ("MinPtMax")),
00045    min_PtMin_   (iConfig.getParameter<double> ("MinPtMin")),
00046    min_InvMass_ (iConfig.getParameter<double> ("MinInvMass")),
00047    max_InvMass_ (iConfig.getParameter<double> ("MaxInvMass")),
00048    min_Acop_    (iConfig.getParameter<double> ("MinAcop")),
00049    max_Acop_    (iConfig.getParameter<double> ("MaxAcop")),
00050    min_PtBalance_ (iConfig.getParameter<double> ("MinPtBalance")),
00051    max_PtBalance_ (iConfig.getParameter<double> ("MaxPtBalance")),
00052    nsigma_Pt_   (iConfig.getParameter<double> ("NSigmaPt")), 
00053    max_DzMuMu_  (iConfig.getParameter<double>("MaxDzMuMu")),
00054    max_YPair_   (iConfig.getParameter<double>("MaxRapidityPair")),
00055    saveTag_  (iConfig.getUntrackedParameter<bool> ("SaveTag",false)) 
00056 {
00057 
00058    LogDebug("HLTMuonDimuonL3Filter")
00059       << " CandTag/MinN/MaxEta/MinNhits/MaxDr/MaxDz/MinPt1/MinPt2/MinInvMass/MaxInvMass/MinAcop/MaxAcop/MinPtBalance/MaxPtBalance/NSigmaPt/MaxDzMuMu/MaxRapidityPair : " 
00060       << candTag_.encode()
00061       << " " << fast_Accept_
00062       << " " << max_Eta_
00063       << " " << min_Nhits_
00064       << " " << max_Dr_
00065       << " " << max_Dz_
00066       << " " << chargeOpt_ << " " << min_PtPair_
00067       << " " << min_PtMax_ << " " << min_PtMin_
00068       << " " << min_InvMass_ << " " << max_InvMass_
00069       << " " << min_Acop_ << " " << max_Acop_
00070       << " " << min_PtBalance_ << " " << max_PtBalance_
00071       << " " << nsigma_Pt_
00072       << " " << max_DzMuMu_
00073       << " " << max_YPair_;
00074 
00075    //register your products
00076    produces<trigger::TriggerFilterObjectWithRefs>();
00077 }
00078 
00079 HLTMuonDimuonL3Filter::~HLTMuonDimuonL3Filter()
00080 {
00081 }
00082 
00083 //
00084 // member functions
00085 //
00086 
00087 // ------------ method called to produce the data  ------------
00088 bool
00089 HLTMuonDimuonL3Filter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
00090 {
00091 
00092    double const MuMass = 0.106;
00093    double const MuMass2 = MuMass*MuMass;
00094    // All HLT filters must create and fill an HLT filter object,
00095    // recording any reconstructed physics objects satisfying (or not)
00096    // this HLT filter, and place it in the Event.
00097 
00098    // The filter object
00099    auto_ptr<TriggerFilterObjectWithRefs>
00100      filterproduct (new TriggerFilterObjectWithRefs(path(),module()));
00101 
00102    // get hold of trks
00103    Handle<RecoChargedCandidateCollection> mucands;
00104    if(saveTag_)filterproduct->addCollectionTag(candTag_);
00105    iEvent.getByLabel (candTag_,mucands);
00106    // sort them by L2Track
00107    std::map<reco::TrackRef, std::vector<RecoChargedCandidateRef> > L2toL3s;
00108    unsigned int maxI = mucands->size();
00109    for (unsigned int i=0;i!=maxI;i++){
00110      TrackRef tk = (*mucands)[i].track();
00111      edm::Ref<L3MuonTrajectorySeedCollection> l3seedRef = tk->seedRef().castTo<edm::Ref<L3MuonTrajectorySeedCollection> >();
00112      TrackRef staTrack = l3seedRef->l2Track();
00113      L2toL3s[staTrack].push_back(RecoChargedCandidateRef(mucands,i));
00114    }
00115 
00116    Handle<TriggerFilterObjectWithRefs> previousLevelCands;
00117    iEvent.getByLabel (previousCandTag_,previousLevelCands);
00118    BeamSpot beamSpot;
00119    Handle<BeamSpot> recoBeamSpotHandle;
00120    iEvent.getByLabel(beamspotTag_,recoBeamSpotHandle);
00121    beamSpot = *recoBeamSpotHandle;
00122   
00123    // needed to compare to L2
00124    vector<RecoChargedCandidateRef> vl2cands;
00125    previousLevelCands->getObjects(TriggerMuon,vl2cands);
00126 
00127    // look at all mucands,  check cuts and add to filter object
00128    int n = 0;
00129    double e1,e2;
00130    Particle::LorentzVector p,p1,p2;
00131 
00132    std::map<reco::TrackRef, std::vector<RecoChargedCandidateRef> > ::iterator L2toL3s_it1 = L2toL3s.begin();
00133    std::map<reco::TrackRef, std::vector<RecoChargedCandidateRef> > ::iterator L2toL3s_end = L2toL3s.end();
00134    bool atLeastOnePair=false;
00135    for (; L2toL3s_it1!=L2toL3s_end; ++L2toL3s_it1){
00136 
00137      if (!triggeredByLevel2(L2toL3s_it1->first,vl2cands)) continue;
00138      
00139      //loop over the L3Tk reconstructed for this L2.
00140      unsigned int iTk1=0;
00141      unsigned int maxItk1=L2toL3s_it1->second.size();
00142      for (; iTk1!=maxItk1; iTk1++){
00143        bool thisL3Index1isDone=false;
00144        RecoChargedCandidateRef & cand1=L2toL3s_it1->second[iTk1];
00145        TrackRef tk1 = cand1->get<TrackRef>();
00146        // eta cut
00147        LogDebug("HLTMuonDimuonL3Filter") << " 1st muon in loop: q*pt= "
00148                                          << tk1->charge()*tk1->pt() << ", eta= " << tk1->eta() << ", hits= " << tk1->numberOfValidHits();
00149 
00150        if (fabs(tk1->eta())>max_Eta_) continue;
00151        
00152        // cut on number of hits
00153        if (tk1->numberOfValidHits()<min_Nhits_) continue;
00154        
00155        //dr cut
00156        //      if (fabs(tk1->d0())>max_Dr_) continue;
00157        if (fabs(tk1->dxy(beamSpot.position()))>max_Dr_) continue;
00158        
00159        //dz cut
00160        if (fabs(tk1->dz())>max_Dz_) continue;
00161        
00162        // Pt threshold cut
00163        double pt1 = tk1->pt();
00164        double err1 = tk1->error(0);
00165        double abspar1 = fabs(tk1->parameter(0));
00166        double ptLx1 = pt1;
00167        // convert 50% efficiency threshold to 90% efficiency threshold
00168        if (abspar1>0) ptLx1 += nsigma_Pt_*err1/abspar1*pt1;
00169        LogDebug("HLTMuonDimuonL3Filter") << " ... 1st muon in loop, pt1= "
00170                                          << pt1 << ", ptLx1= " << ptLx1;
00171        std::map<reco::TrackRef, std::vector<RecoChargedCandidateRef> > ::iterator L2toL3s_it2 = L2toL3s_it1;
00172        L2toL3s_it2++;
00173        for (; L2toL3s_it2!=L2toL3s_end; ++L2toL3s_it2){
00174          if (!triggeredByLevel2(L2toL3s_it2->first,vl2cands)) continue;
00175          
00176             //loop over the L3Tk reconstructed for this L2.
00177             unsigned int iTk2=0;
00178             unsigned int maxItk2=L2toL3s_it2->second.size();
00179             for (; iTk2!=maxItk2; iTk2++){
00180               RecoChargedCandidateRef & cand2=L2toL3s_it2->second[iTk2];
00181               TrackRef tk2 = cand2->get<TrackRef>();
00182               
00183               // eta cut
00184               LogDebug("HLTMuonDimuonL3Filter") << " 2nd muon in loop: q*pt= " << tk2->charge()*tk2->pt() << ", eta= " << tk2->eta() << ", hits= " << tk2->numberOfValidHits() << ", d0= " << tk2->d0();
00185               if (fabs(tk2->eta())>max_Eta_) continue;
00186               
00187               // cut on number of hits
00188               if (tk2->numberOfValidHits()<min_Nhits_) continue;
00189               
00190               //dr cut
00191               // if (fabs(tk2->d0())>max_Dr_) continue;
00192               if (fabs(tk2->dxy(beamSpot.position()))>max_Dr_) continue;
00193               
00194               //dz cut
00195               if (fabs(tk2->dz())>max_Dz_) continue;
00196               
00197               // Pt threshold cut
00198               double pt2 = tk2->pt();
00199               double err2 = tk2->error(0);
00200               double abspar2 = fabs(tk2->parameter(0));
00201               double ptLx2 = pt2;
00202               // convert 50% efficiency threshold to 90% efficiency threshold
00203               if (abspar2>0) ptLx2 += nsigma_Pt_*err2/abspar2*pt2;
00204               LogDebug("HLTMuonDimuonL3Filter") << " ... 2nd muon in loop, pt2= "
00205                                                 << pt2 << ", ptLx2= " << ptLx2;
00206               
00207               if (ptLx1>ptLx2) {
00208                 if (ptLx1<min_PtMax_) continue;
00209                 if (ptLx2<min_PtMin_) continue;
00210               } else {
00211                 if (ptLx2<min_PtMax_) continue;
00212                 if (ptLx1<min_PtMin_) continue;
00213               }
00214               
00215               if (chargeOpt_<0) {
00216                 if (tk1->charge()*tk2->charge()>0) continue;
00217               } else if (chargeOpt_>0) {
00218                 if (tk1->charge()*tk2->charge()<0) continue;
00219               }
00220               
00221               // Acoplanarity
00222               double acop = fabs(tk1->phi()-tk2->phi());
00223               if (acop>M_PI) acop = 2*M_PI - acop;
00224               acop = M_PI - acop;
00225               LogDebug("HLTMuonDimuonL3Filter") << " ... 1-2 acop= " << acop;
00226               if (acop<min_Acop_) continue;
00227               if (acop>max_Acop_) continue;
00228 
00229               // Pt balance
00230               double ptbalance = fabs(tk1->pt()-tk2->pt());
00231               if (ptbalance<min_PtBalance_) continue;
00232               if (ptbalance>max_PtBalance_) continue;
00233 
00234               // Combined dimuon system
00235               e1 = sqrt(tk1->momentum().Mag2()+MuMass2);
00236               e2 = sqrt(tk2->momentum().Mag2()+MuMass2);
00237               p1 = Particle::LorentzVector(tk1->px(),tk1->py(),tk1->pz(),e1);
00238               p2 = Particle::LorentzVector(tk2->px(),tk2->py(),tk2->pz(),e2);
00239               p = p1+p2;
00240               
00241               double pt12 = p.pt();
00242               LogDebug("HLTMuonDimuonL3Filter") << " ... 1-2 pt12= " << pt12;
00243               if (pt12<min_PtPair_) continue;
00244               
00245               double invmass = abs(p.mass());
00246               // if (invmass>0) invmass = sqrt(invmass); else invmass = 0;
00247               LogDebug("HLTMuonDimuonL3Filter") << " ... 1-2 invmass= " << invmass;
00248               if (invmass<min_InvMass_) continue;
00249               if (invmass>max_InvMass_) continue;
00250 
00251               // Delta Z between the two muons
00252               double DeltaZMuMu = fabs(tk2->dz(beamSpot.position())-tk1->dz(beamSpot.position()));
00253               if ( DeltaZMuMu > max_DzMuMu_) continue;
00254               
00255               // Max dimuon |rapidity|
00256               double rapidity = fabs(p.Rapidity());
00257               if ( rapidity > max_YPair_) continue;
00258               
00259               // Add this pair
00260               n++;
00261               LogDebug("HLTMuonDimuonL3Filter") << " Track1 passing filter: pt= " << tk1->pt() << ", eta: " << tk1->eta();
00262               LogDebug("HLTMuonDimuonL3Filter") << " Track2 passing filter: pt= " << tk2->pt() << ", eta: " << tk2->eta();
00263               LogDebug("HLTMuonDimuonL3Filter") << " Invmass= " << invmass;
00264 
00265               bool i1done = false;
00266               bool i2done = false;
00267               vector<RecoChargedCandidateRef> vref;
00268               filterproduct->getObjects(TriggerMuon,vref);
00269               for (unsigned int i=0; i<vref.size(); i++) {
00270                 RecoChargedCandidateRef candref =  RecoChargedCandidateRef(vref[i]);
00271                 TrackRef tktmp = candref->get<TrackRef>();
00272                 if (tktmp==tk1) {
00273                   i1done = true;
00274                 } else if (tktmp==tk2) {
00275                   i2done = true;
00276                 }
00277                 if (i1done && i2done) break;
00278               }
00279             
00280               if (!i1done) { 
00281                 filterproduct->addObject(TriggerMuon,cand1);
00282               }
00283               if (!i2done) { 
00284                 filterproduct->addObject(TriggerMuon,cand2);
00285               }
00286               
00287               //break anyway since a L3 track pair has been found matching the criteria
00288               thisL3Index1isDone=true;
00289               atLeastOnePair=true;
00290               break;
00291             }//loop on the track of the second L2
00292             //break the loop if fast accept.
00293             if (atLeastOnePair && fast_Accept_) break;
00294        }//loop on the second L2
00295 
00296        
00297        //break the loop if fast accept.
00298        if (atLeastOnePair && fast_Accept_) break;
00299        if (thisL3Index1isDone) break;
00300      }//loop on tracks for first L2
00301      //break the loop if fast accept.
00302      if (atLeastOnePair && fast_Accept_) break;
00303    }//loop on the first L2
00304    
00305    // filter decision
00306    const bool accept (n >= 1);
00307 
00308    // put filter object into the Event
00309    iEvent.put(filterproduct);
00310 
00311    LogDebug("HLTMuonDimuonL3Filter") << " >>>>> Result of HLTMuonDimuonL3Filter is "<< accept << ", number of muon pairs passing thresholds= " << n; 
00312 
00313    return accept;
00314 }
00315 
00316 
00317 bool
00318 HLTMuonDimuonL3Filter::triggeredByLevel2(const TrackRef& staTrack,vector<RecoChargedCandidateRef>& vcands)
00319 {
00320   bool ok=false;
00321   for (unsigned int i=0; i<vcands.size(); i++) {
00322     if ( vcands[i]->get<TrackRef>() == staTrack ) {
00323       ok=true;
00324       LogDebug("HLTMuonL3PreFilter") << "The L2 track triggered";
00325       break;
00326     }
00327   }
00328   return ok;
00329 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/HLTrigger/Muon/src/HLTMuonDimuonL3Filter.cc
