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
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
00078 produces<trigger::TriggerFilterObjectWithRefs>();
00079 }
00080
00081 HLTMuonDimuonL2Filter::~HLTMuonDimuonL2Filter()
00082 {
00083 }
00084
00085
00086
00087
00088
00089
00090
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
00098
00099
00100
00101
00102 auto_ptr<TriggerFilterObjectWithRefs>
00103 filterproduct (new TriggerFilterObjectWithRefs(path(),module()));
00104
00105
00106 RecoChargedCandidateRef ref1;
00107 RecoChargedCandidateRef ref2;
00108
00109
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
00120 HLTMuonL2ToL1Map mapL2ToL1(previousCandTag_, seedMapTag_, iEvent);
00121
00122
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
00132 LogDebug("HLTMuonDimuonL2Filter") << " 1st muon in loop: q*pt= "
00133 << tk1->charge()*tk1->pt() << ", eta= " << tk1->eta() << ", hits= " << tk1->numberOfValidHits();
00134
00135 if (!mapL2ToL1.isTriggeredByL1(tk1)) continue;
00136
00137 if (fabs(tk1->eta())>max_Eta_) continue;
00138
00139
00140 if (tk1->numberOfValidHits()<min_Nhits_) continue;
00141
00142
00143
00144 if (fabs(tk1->dxy(beamSpot.position()))>max_Dr_) continue;
00145
00146
00147 if (fabs(tk1->dz())>max_Dz_) continue;
00148
00149
00150 double pt1 = tk1->pt();
00151 double err1 = tk1->error(0);
00152 double abspar1 = fabs(tk1->parameter(0));
00153 double ptLx1 = pt1;
00154
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
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
00170 if (tk2->numberOfValidHits()<min_Nhits_) continue;
00171
00172
00173
00174 if (fabs(tk2->dxy(beamSpot.position()))>max_Dr_) continue;
00175
00176
00177 if (fabs(tk2->dz())>max_Dz_) continue;
00178
00179
00180 double pt2 = tk2->pt();
00181 double err2 = tk2->error(0);
00182 double abspar2 = fabs(tk2->parameter(0));
00183 double ptLx2 = pt2;
00184
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
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
00212 double ptbalance = fabs(tk1->pt()-tk2->pt());
00213 if (ptbalance<min_PtBalance_) continue;
00214 if (ptbalance>max_PtBalance_) continue;
00215
00216
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
00229 LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 invmass= " << invmass;
00230 if (invmass<min_InvMass_) continue;
00231 if (invmass>max_InvMass_) continue;
00232
00233
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
00269 const bool accept (n >= 1);
00270
00271
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