CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/HLTrigger/Muon/src/HLTMuonDimuonL2Filter.cc

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/TrajectorySeed/interface/TrajectorySeed.h"
00021 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
00022 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
00023 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
00024 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
00025 
00026 #include "HLTrigger/Muon/interface/HLTMuonDimuonL2Filter.h"
00027 #include "HLTrigger/Muon/interface/HLTMuonL2ToL1Map.h"
00028 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00029 
00030 using namespace edm;
00031 using namespace std;
00032 using namespace reco;
00033 using namespace trigger;
00034 using namespace l1extra;
00035 //
00036 // constructors and destructor
00037 //
00038 HLTMuonDimuonL2Filter::HLTMuonDimuonL2Filter(const edm::ParameterSet& iConfig):
00039    beamspotTag_   (iConfig.getParameter< edm::InputTag > ("BeamSpotTag")),
00040    candTag_     (iConfig.getParameter< edm::InputTag > ("CandTag")),
00041    previousCandTag_   (iConfig.getParameter<InputTag > ("PreviousCandTag")),
00042    seedMapTag_( iConfig.getParameter<InputTag >("SeedMapTag") ),
00043    fast_Accept_ (iConfig.getParameter<bool> ("FastAccept")),
00044    max_Eta_     (iConfig.getParameter<double> ("MaxEta")),
00045    min_Nhits_   (iConfig.getParameter<int> ("MinNhits")),
00046    max_Dr_      (iConfig.getParameter<double> ("MaxDr")),
00047    max_Dz_      (iConfig.getParameter<double> ("MaxDz")),
00048    chargeOpt_   (iConfig.getParameter<int> ("ChargeOpt")),
00049    min_PtPair_  (iConfig.getParameter<double> ("MinPtPair")),
00050    min_PtMax_   (iConfig.getParameter<double> ("MinPtMax")),
00051    min_PtMin_   (iConfig.getParameter<double> ("MinPtMin")),
00052    min_InvMass_ (iConfig.getParameter<double> ("MinInvMass")),
00053    max_InvMass_ (iConfig.getParameter<double> ("MaxInvMass")),
00054    min_Acop_    (iConfig.getParameter<double> ("MinAcop")),
00055    max_Acop_    (iConfig.getParameter<double> ("MaxAcop")),
00056    min_PtBalance_ (iConfig.getParameter<double> ("MinPtBalance")),
00057    max_PtBalance_ (iConfig.getParameter<double> ("MaxPtBalance")),
00058    nsigma_Pt_   (iConfig.getParameter<double> ("NSigmaPt")),
00059    saveTag_  (iConfig.getUntrackedParameter<bool> ("SaveTag",false)) 
00060 {
00061 
00062    LogDebug("HLTMuonDimuonL2Filter")
00063       << " CandTag/MinN/MaxEta/MinNhits/MaxDr/MaxDz/MinPt1/MinPt2/MinInvMass/MaxInvMass/MinAcop/MaxAcop/MinPtBalance/MaxPtBalance/NSigmaPt : " 
00064       << candTag_.encode()
00065       << " " << fast_Accept_
00066       << " " << max_Eta_
00067       << " " << min_Nhits_
00068       << " " << max_Dr_
00069       << " " << max_Dz_
00070       << " " << chargeOpt_ << " " << min_PtPair_
00071       << " " << min_PtMax_ << " " << min_PtMin_
00072       << " " << min_InvMass_ << " " << max_InvMass_
00073       << " " << min_Acop_ << " " << max_Acop_
00074       << " " << min_PtBalance_ << " " << max_PtBalance_
00075       << " " << nsigma_Pt_;
00076 
00077    //register your products
00078    produces<trigger::TriggerFilterObjectWithRefs>();
00079 }
00080 
00081 HLTMuonDimuonL2Filter::~HLTMuonDimuonL2Filter()
00082 {
00083 }
00084 
00085 
00086 //
00087 // member functions
00088 //
00089 
00090 // ------------ method called to produce the data  ------------
00091 bool
00092 HLTMuonDimuonL2Filter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
00093 {
00094 
00095    double const MuMass = 0.106;
00096    double const MuMass2 = MuMass*MuMass;
00097    // All HLT filters must create and fill an HLT filter object,
00098    // recording any reconstructed physics objects satisfying (or not)
00099    // this HLT filter, and place it in the Event.
00100 
00101    // The filter object
00102    auto_ptr<TriggerFilterObjectWithRefs>
00103      filterproduct (new TriggerFilterObjectWithRefs(path(),module()));
00104 
00105    // Ref to Candidate object to be recorded in filter object
00106    RecoChargedCandidateRef ref1;
00107    RecoChargedCandidateRef ref2;
00108 
00109    // get hold of trks
00110    Handle<RecoChargedCandidateCollection> mucands;
00111    iEvent.getByLabel (candTag_,mucands);
00112    if(saveTag_)filterproduct->addCollectionTag(candTag_);
00113 
00114    BeamSpot beamSpot;
00115    Handle<BeamSpot> recoBeamSpotHandle;
00116    iEvent.getByLabel(beamspotTag_,recoBeamSpotHandle);
00117    beamSpot = *recoBeamSpotHandle;
00118 
00119    // get the L2 to L1 map object for this event
00120    HLTMuonL2ToL1Map mapL2ToL1(previousCandTag_, seedMapTag_, iEvent);
00121 
00122    // look at all mucands,  check cuts and add to filter object
00123    int n = 0;
00124    double e1,e2;
00125    Particle::LorentzVector p,p1,p2;
00126 
00127    RecoChargedCandidateCollection::const_iterator cand1;
00128    RecoChargedCandidateCollection::const_iterator cand2;
00129    for (cand1=mucands->begin(); cand1!=mucands->end(); cand1++) {
00130       TrackRef tk1 = cand1->get<TrackRef>();
00131       // eta cut
00132       LogDebug("HLTMuonDimuonL2Filter") << " 1st muon in loop: q*pt= "
00133             << tk1->charge()*tk1->pt() << ", eta= " << tk1->eta() << ", hits= " << tk1->numberOfValidHits();
00134       // find the L1 Particle corresponding to the L2 Track
00135       if (!mapL2ToL1.isTriggeredByL1(tk1)) continue;
00136  
00137       if (fabs(tk1->eta())>max_Eta_) continue;
00138 
00139       // cut on number of hits
00140       if (tk1->numberOfValidHits()<min_Nhits_) continue;
00141 
00142       //dr cut
00143       //if (fabs(tk1->d0())>max_Dr_) continue;
00144       if (fabs(tk1->dxy(beamSpot.position()))>max_Dr_) continue;
00145 
00146       //dz cut
00147       if (fabs(tk1->dz())>max_Dz_) continue;
00148 
00149       // Pt threshold cut
00150       double pt1 = tk1->pt();
00151       double err1 = tk1->error(0);
00152       double abspar1 = fabs(tk1->parameter(0));
00153       double ptLx1 = pt1;
00154       // convert 50% efficiency threshold to 90% efficiency threshold
00155       if (abspar1>0) ptLx1 += nsigma_Pt_*err1/abspar1*pt1;
00156       LogDebug("HLTMuonDimuonL2Filter") << " ... 1st muon in loop, pt1= "
00157             << pt1 << ", ptLx1= " << ptLx1;
00158 
00159       cand2 = cand1; cand2++;
00160       for (; cand2!=mucands->end(); cand2++) {
00161             TrackRef tk2 = cand2->get<TrackRef>();
00162 
00163             // eta cut
00164             LogDebug("HLTMuonDimuonL2Filter") << " 2nd muon in loop: q*pt= " << tk2->charge()*tk2->pt() << ", eta= " << tk2->eta() << ", hits= " << tk2->numberOfValidHits() << ", d0= " << tk2->d0();
00165             if (!mapL2ToL1.isTriggeredByL1(tk2)) continue;
00166 
00167             if (fabs(tk2->eta())>max_Eta_) continue;
00168 
00169             // cut on number of hits
00170             if (tk2->numberOfValidHits()<min_Nhits_) continue;
00171 
00172             //dr cut
00173             //if (fabs(tk2->d0())>max_Dr_) continue;
00174             if (fabs(tk2->dxy(beamSpot.position()))>max_Dr_) continue;
00175 
00176             //dz cut
00177             if (fabs(tk2->dz())>max_Dz_) continue;
00178 
00179             // Pt threshold cut
00180             double pt2 = tk2->pt();
00181             double err2 = tk2->error(0);
00182             double abspar2 = fabs(tk2->parameter(0));
00183             double ptLx2 = pt2;
00184             // convert 50% efficiency threshold to 90% efficiency threshold
00185             if (abspar2>0) ptLx2 += nsigma_Pt_*err2/abspar2*pt2;
00186             LogDebug("HLTMuonDimuonL2Filter") << " ... 2nd muon in loop, pt2= "
00187                   << pt2 << ", ptLx2= " << ptLx2;
00188 
00189             if (ptLx1>ptLx2) {
00190                   if (ptLx1<min_PtMax_) continue;
00191                   if (ptLx2<min_PtMin_) continue;
00192             } else {
00193                   if (ptLx2<min_PtMax_) continue;
00194                   if (ptLx1<min_PtMin_) continue;
00195             }
00196 
00197             if (chargeOpt_<0) {
00198                   if (tk1->charge()*tk2->charge()>0) continue;
00199             } else if (chargeOpt_>0) {
00200                   if (tk1->charge()*tk2->charge()<0) continue;
00201             }
00202 
00203             // Acoplanarity
00204             double acop = fabs(tk1->phi()-tk2->phi());
00205             if (acop>M_PI) acop = 2*M_PI - acop;
00206             acop = M_PI - acop;
00207             LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 acop= " << acop;
00208             if (acop<min_Acop_) continue;
00209             if (acop>max_Acop_) continue;
00210 
00211             // Pt balance
00212             double ptbalance = fabs(tk1->pt()-tk2->pt());
00213             if (ptbalance<min_PtBalance_) continue;
00214             if (ptbalance>max_PtBalance_) continue;
00215 
00216             // Combined dimuon system
00217             e1 = sqrt(tk1->momentum().Mag2()+MuMass2);
00218             e2 = sqrt(tk2->momentum().Mag2()+MuMass2);
00219             p1 = Particle::LorentzVector(tk1->px(),tk1->py(),tk1->pz(),e1);
00220             p2 = Particle::LorentzVector(tk2->px(),tk2->py(),tk2->pz(),e2);
00221             p = p1+p2;
00222 
00223             double pt12 = p.pt();
00224             LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 pt12= " << pt12;
00225             if (pt12<min_PtPair_) continue;
00226 
00227             double invmass = abs(p.mass());
00228          // if (invmass>0) invmass = sqrt(invmass); else invmass = 0;
00229             LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 invmass= " << invmass;
00230             if (invmass<min_InvMass_) continue;
00231             if (invmass>max_InvMass_) continue;
00232 
00233             // Add this pair
00234             n++;
00235             LogDebug("HLTMuonDimuonL2Filter") << " Track1 passing filter: pt= " << tk1->pt() << ", eta: " << tk1->eta();
00236             LogDebug("HLTMuonDimuonL2Filter") << " Track2 passing filter: pt= " << tk2->pt() << ", eta: " << tk2->eta();
00237             LogDebug("HLTMuonDimuonL2Filter") << " Invmass= " << invmass;
00238 
00239             bool i1done = false;
00240             bool i2done = false;
00241             vector<RecoChargedCandidateRef> vref;
00242             filterproduct->getObjects(TriggerMuon,vref);
00243             for (unsigned int i=0; i<vref.size(); i++) {
00244               RecoChargedCandidateRef candref =  RecoChargedCandidateRef(vref[i]);
00245               TrackRef tktmp = candref->get<TrackRef>();
00246               if (tktmp==tk1) {
00247                         i1done = true;
00248               } else if (tktmp==tk2) {
00249                 i2done = true;
00250               }
00251               if (i1done && i2done) break;
00252             }
00253             
00254             if (!i1done) { 
00255               ref1=RecoChargedCandidateRef( Ref<RecoChargedCandidateCollection> (mucands,distance(mucands->begin(), cand1)));
00256               filterproduct->addObject(TriggerMuon,ref1);
00257             }
00258             if (!i2done) { 
00259               ref2=RecoChargedCandidateRef( Ref<RecoChargedCandidateCollection> (mucands,distance(mucands->begin(),cand2 )));
00260             filterproduct->addObject(TriggerMuon,ref2);
00261             }
00262 
00263             if (fast_Accept_) break;
00264       }
00265 
00266    }
00267 
00268    // filter decision
00269    const bool accept (n >= 1);
00270 
00271    // put filter object into the Event
00272    iEvent.put(filterproduct);
00273 
00274    LogDebug("HLTMuonDimuonL2Filter") << " >>>>> Result of HLTMuonDimuonL2Filter is "<< accept << ", number of muon pairs passing thresholds= " << n; 
00275 
00276    return accept;
00277 }
00278