CMS 3D CMS Logo

BPHMonitor.cc
Go to the documentation of this file.
2 
4 
6 
7 
8 // -----------------------------
9 // constructors and destructor
10 // -----------------------------
11 
13  folderName_ ( iConfig.getParameter<std::string>("FolderName") )
14  , muoToken_ ( mayConsume<reco::MuonCollection> (iConfig.getParameter<edm::InputTag>("muons") ) )
15  , bsToken_ ( mayConsume<reco::BeamSpot> (iConfig.getParameter<edm::InputTag>("beamSpot")))
16  , trToken_ ( mayConsume<reco::TrackCollection> (iConfig.getParameter<edm::InputTag>("tracks")))
17  , phToken_ ( mayConsume<reco::PhotonCollection> (iConfig.getParameter<edm::InputTag>("photons")))
18  , phi_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("phiPSet") ) )
19  , pt_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("ptPSet") ) )
20  , eta_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("etaPSet") ) )
21  , d0_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("d0PSet") ) )
22  , z0_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("z0PSet") ) )
23  , dR_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("dRPSet") ) )
24  , mass_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("massPSet") ) )
25  , dca_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("dcaPSet") ) )
26  , ds_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("dsPSet") ) )
27  , cos_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("cosPSet") ) )
28  , prob_binning_ ( getHistoPSet (iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet> ("probPSet") ) )
29  , num_genTriggerEventFlag_(new GenericTriggerEventFlag(iConfig.getParameter<edm::ParameterSet>("numGenericTriggerEventPSet"),consumesCollector(), *this))
30  , den_genTriggerEventFlag_(new GenericTriggerEventFlag(iConfig.getParameter<edm::ParameterSet>("denGenericTriggerEventPSet"),consumesCollector(), *this))
31  , prescaleWeightProvider_( new PrescaleWeightProvider( iConfig.getParameter<edm::ParameterSet>("PrescaleTriggerEventPSet"),consumesCollector(), *this))
32  , muoSelection_ ( iConfig.getParameter<std::string>("muoSelection") )
33  , muoSelection_ref ( iConfig.getParameter<std::string>("muoSelection_ref") )
34  , muoSelection_tag ( iConfig.getParameter<std::string>("muoSelection_tag") )
35  , muoSelection_probe ( iConfig.getParameter<std::string>("muoSelection_probe") )
36  , nmuons_ ( iConfig.getParameter<int>("nmuons" ) )
37  , tnp_ ( iConfig.getParameter<bool>("tnp" ) )
38  , L3_ ( iConfig.getParameter<int>("L3" ) )
39  , trOrMu_ ( iConfig.getParameter<int>("trOrMu" ) )
40  , Jpsi_ ( iConfig.getParameter<int>("Jpsi" ) )
41  , Upsilon_ ( iConfig.getParameter<int>("Upsilon" ) ) // if ==1 path with Upsilon constraint
42  , enum_ ( iConfig.getParameter<int>("enum" ) )
43  , seagull_ ( iConfig.getParameter<int>("seagull" ) )
44  , maxmass_ ( iConfig.getParameter<double>("maxmass" ) )
45  , minmass_ ( iConfig.getParameter<double>("minmass" ) )
46  , maxmassJpsi ( iConfig.getParameter<double>("maxmassJpsi" ) )
47  , minmassJpsi ( iConfig.getParameter<double>("minmassJpsi" ) )
48  , maxmassUpsilon ( iConfig.getParameter<double>("maxmassUpsilon" ) )
49  , minmassUpsilon ( iConfig.getParameter<double>("minmassUpsilon" ) )
50  , maxmassTkTk ( iConfig.getParameter<double>("maxmassTkTk" ) )
51  , minmassTkTk ( iConfig.getParameter<double>("minmassTkTk" ) )
52  , maxmassJpsiTk ( iConfig.getParameter<double>("maxmassJpsiTk" ) )
53  , minmassJpsiTk ( iConfig.getParameter<double>("minmassJpsiTk" ) )
54  , kaon_mass ( iConfig.getParameter<double>("kaon_mass" ) )
55  , mu_mass ( iConfig.getParameter<double>("mu_mass" ) )
56  , min_dR ( iConfig.getParameter<double>("min_dR" ) )
57  , minprob ( iConfig.getParameter<double>("minprob" ) )
58  , mincos ( iConfig.getParameter<double>("mincos" ) )
59  , minDS ( iConfig.getParameter<double>("minDS" ) )
60  , hltTrigResTag_ (mayConsume<edm::TriggerResults>( iConfig.getParameter<edm::ParameterSet>("denGenericTriggerEventPSet").getParameter<edm::InputTag>("hltInputTag")))
61  , hltInputTag_ (mayConsume<trigger::TriggerEvent>( iConfig.getParameter<edm::InputTag>("hltTriggerSummaryAOD")))
62  , hltpaths_num ( iConfig.getParameter<edm::ParameterSet>("numGenericTriggerEventPSet").getParameter<std::vector<std::string>>("hltPaths"))
63  , hltpaths_den ( iConfig.getParameter<edm::ParameterSet>("denGenericTriggerEventPSet").getParameter<std::vector<std::string>>("hltPaths"))
64  , trSelection_ ( iConfig.getParameter<std::string>("muoSelection") )
65  , trSelection_ref ( iConfig.getParameter<std::string>("trSelection_ref") )
66  , DMSelection_ref ( iConfig.getParameter<std::string>("DMSelection_ref") )
67 {
68 
69  muPhi_.numerator = nullptr;
70  muPhi_.denominator = nullptr;
71  muEta_.numerator = nullptr;
72  muEta_.denominator = nullptr;
73  muPt_.numerator = nullptr;
74  muPt_.denominator = nullptr;
75  mud0_.numerator = nullptr;
76  mud0_.denominator = nullptr;
77  muz0_.numerator = nullptr;
78  muz0_.denominator = nullptr;
79 
80  mu1Phi_.numerator = nullptr;
81  mu1Phi_.denominator = nullptr;
82  mu1Eta_.numerator = nullptr;
83  mu1Eta_.denominator = nullptr;
84  mu1Pt_.numerator = nullptr;
85  mu1Pt_.denominator = nullptr;
86 
87  mu2Phi_.numerator = nullptr;
88  mu2Phi_.denominator = nullptr;
89  mu2Eta_.numerator = nullptr;
90  mu2Eta_.denominator = nullptr;
91  mu2Pt_.numerator = nullptr;
92  mu2Pt_.denominator = nullptr;
93 
94  mu3Phi_.numerator = nullptr;
95  mu3Phi_.denominator = nullptr;
96  mu3Eta_.numerator = nullptr;
97  mu3Eta_.denominator = nullptr;
98  mu3Pt_.numerator = nullptr;
99  mu3Pt_.denominator = nullptr;
100 
101  phPhi_.numerator = nullptr;
102  phPhi_.denominator = nullptr;
103  phEta_.numerator = nullptr;
104  phEta_.denominator = nullptr;
105  phPt_.numerator = nullptr;
106  phPt_.denominator = nullptr;
107 
108 
109  DiMuPhi_.numerator = nullptr;
110  DiMuPhi_.denominator = nullptr;
111  DiMuEta_.numerator = nullptr;
112  DiMuEta_.denominator = nullptr;
113  DiMuPt_.numerator = nullptr;
114  DiMuPt_.denominator = nullptr;
115  DiMuPVcos_.numerator = nullptr;
116  DiMuPVcos_.denominator = nullptr;
117  DiMuProb_.numerator = nullptr;
118  DiMuProb_.denominator = nullptr;
119  DiMuDS_.numerator = nullptr;
120  DiMuDS_.denominator = nullptr;
121  DiMuDCA_.numerator = nullptr;
122  DiMuDCA_.denominator = nullptr;
123  DiMuMass_.numerator = nullptr;
124  DiMuMass_.denominator = nullptr;
125  DiMudR_.numerator = nullptr;
126  DiMudR_.denominator = nullptr;
127 
128 
129 }
130 
132 {
136 }
137 
139 {
140  // Due to the setup of the fillDescription only one of the
141  // two cases is possible at this point.
142  if (pset.existsAs<std::vector<double>>("edges")) {
143  return MEbinning{pset.getParameter<std::vector<double>>("edges")};
144  }
145 
146  return MEbinning {
147  pset.getParameter<int32_t>("nbins"),
148  pset.getParameter<double>("xmin"),
149  pset.getParameter<double>("xmax"),
150  };
151 }
152 
153 // MEbinning BPHMonitor::getHistoLSPSet(edm::ParameterSet pset)
154 // {
155 // return MEbinning{
156 // pset.getParameter<int32_t>("nbins"),
157 // 0.,
158 // double(pset.getParameter<int32_t>("nbins"))
159 // };
160 // }
161 
163 {
164  me.numerator->setAxisTitle(titleX,1);
165  me.numerator->setAxisTitle(titleY,2);
166  me.denominator->setAxisTitle(titleX,1);
167  me.denominator->setAxisTitle(titleY,2);
168 
169 }
170 
171 void BPHMonitor::bookME(DQMStore::IBooker &ibooker, METME& me, std::string& histname, std::string& histtitle, int& nbins, double& min, double& max)
172 {
173  me.numerator = ibooker.book1D(histname+"_numerator", histtitle+" (numerator)", nbins, min, max);
174  me.denominator = ibooker.book1D(histname+"_denominator", histtitle+" (denominator)", nbins, min, max);
175 }
176 void BPHMonitor::bookME(DQMStore::IBooker &ibooker, METME& me, std::string& histname, std::string& histtitle, std::vector<double> binning)
177 {
178  int nbins = binning.size()-1;
179  std::vector<float> fbinning(binning.begin(),binning.end());
180  float* arr = &fbinning[0];
181  me.numerator = ibooker.book1D(histname+"_numerator", histtitle+" (numerator)", nbins, arr);
182  me.denominator = ibooker.book1D(histname+"_denominator", histtitle+" (denominator)", nbins, arr);
183 }
184 void BPHMonitor::bookME(DQMStore::IBooker &ibooker, METME& me, std::string& histname, std::string& histtitle, int& nbinsX, double& xmin, double& xmax, double& ymin, double& ymax)
185 {
186  me.numerator = ibooker.bookProfile(histname+"_numerator", histtitle+" (numerator)", nbinsX, xmin, xmax, ymin, ymax);
187  me.denominator = ibooker.bookProfile(histname+"_denominator", histtitle+" (denominator)", nbinsX, xmin, xmax, ymin, ymax);
188 }
189 void BPHMonitor::bookME(DQMStore::IBooker &ibooker, METME& me, std::string& histname, std::string& histtitle, int& nbinsX, double& xmin, double& xmax, int& nbinsY, double& ymin, double& ymax)
190 {
191  me.numerator = ibooker.book2D(histname+"_numerator", histtitle+" (numerator)", nbinsX, xmin, xmax, nbinsY, ymin, ymax);
192  me.denominator = ibooker.book2D(histname+"_denominator", histtitle+" (denominator)", nbinsX, xmin, xmax, nbinsY, ymin, ymax);
193 }
194 void BPHMonitor::bookME(DQMStore::IBooker &ibooker, METME& me, std::string& histname, std::string& histtitle, std::vector<double> binningX, std::vector<double> binningY)
195 {
196  int nbinsX = binningX.size()-1;
197  std::vector<float> fbinningX(binningX.begin(),binningX.end());
198  float* arrX = &fbinningX[0];
199  int nbinsY = binningY.size()-1;
200  std::vector<float> fbinningY(binningY.begin(),binningY.end());
201  float* arrY = &fbinningY[0];
202 
203  me.numerator = ibooker.book2D(histname+"_numerator", histtitle+" (numerator)", nbinsX, arrX, nbinsY, arrY);
204  me.denominator = ibooker.book2D(histname+"_denominator", histtitle+" (denominator)", nbinsX, arrX, nbinsY, arrY);
205 }
206 
207 void BPHMonitor::bookME(DQMStore::IBooker &ibooker, METME &me, std::string &histname, std::string &histtitle, /*const*/ MEbinning& binning)
208 {
209  // If the vector in the binning is filled use the bins defined there
210  // otherwise use a linear binning between min and max
211  if (binning.edges.empty()) {
212  this->bookME(ibooker, me, histname, histtitle, binning.nbins, binning.xmin, binning.xmax);
213  } else {
214  this->bookME(ibooker, me, histname, histtitle, binning.edges);
215  }
216 }
217 
218 
220  edm::Run const & iRun,
221  edm::EventSetup const & iSetup)
222 {
223 
224  std::string histname, histtitle, istnp, trMuPh;
225  bool Ph_ = false; if (enum_ == 7) Ph_ = true;
226  if (tnp_) istnp = "Tag_and_Probe/"; else istnp = "";
227  std::string currentFolder = folderName_ + istnp;
228  ibooker.setCurrentFolder(currentFolder);
229  if (trOrMu_) trMuPh = "tr"; else if (Ph_) trMuPh = "ph"; else trMuPh = "mu";
230 
231  if (enum_ == 7 || enum_ == 1 || enum_ == 9 || enum_ == 10) {
232  histname = trMuPh+"Pt"; histtitle = trMuPh+"_P_{t}";
233  bookME(ibooker,muPt_,histname,histtitle, pt_binning_);
234  setMETitle(muPt_,trMuPh+"_Pt[GeV]","events / 1 GeV");
235 
236  histname = trMuPh+"Phi"; histtitle = trMuPh+"Phi";
237  bookME(ibooker,muPhi_,histname,histtitle, phi_binning_);
238  setMETitle(muPhi_,trMuPh+"_#phi","events / 0.1 rad");
239 
240  histname = trMuPh+"Eta"; histtitle = trMuPh+"_Eta";
241  bookME(ibooker,muEta_,histname,histtitle, eta_binning_);
242  setMETitle(muEta_,trMuPh+"_#eta","events / 0.2");
243  }
244  else {
245  histname = trMuPh+"1Pt"; histtitle = trMuPh+"1_P_{t}";
246  bookME(ibooker,mu1Pt_,histname,histtitle, pt_binning_);
247  setMETitle(mu1Pt_,trMuPh+"_Pt[GeV]","events / 1 GeV");
248 
249  histname = trMuPh+"1Phi"; histtitle = trMuPh+"1Phi";
250  bookME(ibooker,mu1Phi_,histname,histtitle, phi_binning_);
251  setMETitle(mu1Phi_,trMuPh+"_#phi","events / 0.1 rad");
252 
253  histname = trMuPh+"1Eta"; histtitle = trMuPh+"1_Eta";
254  bookME(ibooker,mu1Eta_,histname,histtitle, eta_binning_);
255  setMETitle(mu1Eta_,trMuPh+"_#eta","events / 0.2");
256 
257  histname = trMuPh+"2Pt"; histtitle = trMuPh+"2_P_{t}";
258  bookME(ibooker,mu2Pt_,histname,histtitle, pt_binning_);
259  setMETitle(mu2Pt_,trMuPh+"_Pt[GeV]","events / 1 GeV");
260 
261  histname = trMuPh+"2Phi"; histtitle = trMuPh+"2Phi";
262  bookME(ibooker,mu2Phi_,histname,histtitle, phi_binning_);
263  setMETitle(mu2Phi_,trMuPh+"_#phi","events / 0.1 rad");
264 
265  histname = trMuPh+"2Eta"; histtitle = trMuPh+"2_Eta";
266  bookME(ibooker,mu2Eta_,histname,histtitle, eta_binning_);
267  setMETitle(mu2Eta_,trMuPh+"_#eta","events / 0.2");
268 
269  if (enum_ == 6) {
270  histname = trMuPh+"3Eta"; histtitle = trMuPh+"3Eta";
271  bookME(ibooker,mu3Eta_,histname,histtitle, eta_binning_);
272  setMETitle(mu3Eta_,trMuPh+"3#eta","events / 0.2");
273 
274  histname = trMuPh+"3Pt"; histtitle = trMuPh+"3_P_{t}";
275  bookME(ibooker,mu3Pt_,histname,histtitle, pt_binning_);
276  setMETitle(mu3Pt_,trMuPh+"3_Pt[GeV]","events / 1 GeV");
277 
278  histname = trMuPh+"3Phi"; histtitle = trMuPh+"3Phi";
279  bookME(ibooker,mu3Phi_,histname,histtitle, phi_binning_);
280  setMETitle(mu3Phi_,trMuPh+"3_#phi","events / 0.1 rad");
281  }
282  else if (enum_ == 2 || enum_ == 4 || enum_ == 5 || enum_ == 8) {
283  histname = "DiMuEta"; histtitle = "DiMuEta";
284  bookME(ibooker,DiMuEta_,histname,histtitle, eta_binning_);
285  setMETitle(DiMuEta_,"DiMu#eta","events / 0.2");
286 
287  histname = "DiMuPt"; histtitle = "DiMu_P_{t}";
288  bookME(ibooker,DiMuPt_,histname,histtitle, pt_binning_);
289  setMETitle(DiMuPt_,"DiMu_Pt[GeV]","events / 1 GeV");
290 
291  histname = "DiMuPhi"; histtitle = "DiMuPhi";
292  bookME(ibooker,DiMuPhi_,histname,histtitle, phi_binning_);
293  setMETitle(DiMuPhi_,"DiMu_#phi","events / 0.1 rad");
294 
295  if (enum_ == 4 || enum_ == 5) {
296  histname = "DiMudR"; histtitle = "DiMudR";
297  bookME(ibooker,DiMudR_,histname,histtitle, dR_binning_);
298  setMETitle(DiMudR_,"DiMu_#dR","events /");
299  if (enum_ == 4) {
300  histname = "DiMuMass"; histtitle = "DiMuMass";
301  bookME(ibooker,DiMuMass_,histname,histtitle, mass_binning_);
302  setMETitle(DiMuMass_,"DiMu_#mass","events /");
303  }
304  } else if (enum_ == 8) {
305  histname = "DiMuProb"; histtitle = "DiMuProb";
306  bookME(ibooker,DiMuProb_,histname,histtitle, prob_binning_);
307  setMETitle(DiMuProb_,"DiMu_#prob","events /");
308 
309  histname = "DiMuPVcos"; histtitle = "DiMuPVcos";
310  bookME(ibooker,DiMuPVcos_,histname,histtitle, cos_binning_);
311  setMETitle(DiMuPVcos_,"DiMu_#cosPV","events /");
312 
313  histname = "DiMuDS"; histtitle = "DiMuDS";
314  bookME(ibooker,DiMuDS_,histname,histtitle, ds_binning_);
315  setMETitle(DiMuDS_,"DiMu_#ds","events /");
316 
317  histname = "DiMuDCA"; histtitle = "DiMuDCA";
318  bookME(ibooker,DiMuDCA_,histname,histtitle, dca_binning_);
319  setMETitle(DiMuDCA_,"DiMu_#dca","events /");
320  }
321  } // if (enum_ == 2 || enum_ == 4 || enum_ == 5 || enum_ == 8)
322  }
323 
324  if (trOrMu_) {
325  if (false) { // not filled anywhere
326  histname = trMuPh+"_d0"; histtitle = trMuPh+"_d0";
327  bookME(ibooker,mud0_,histname,histtitle, d0_binning_);
328  setMETitle(mud0_,trMuPh+"_d0","events / bin");
329 
330  histname = trMuPh+"_z0"; histtitle = trMuPh+"_z0";
331  bookME(ibooker,muz0_,histname,histtitle, z0_binning_);
332  setMETitle(muz0_,trMuPh+"_z0","events / bin");
333  }
334  }
335 
336  // Initialize the GenericTriggerEventFlag
339  prescaleWeightProvider_->initRun( iRun, iSetup );
340 }
341 
346 #include "TLorentzVector.h"
348 
350  iEvent.getByToken( bsToken_, beamSpot);
351 
353  iEvent.getByToken( muoToken_, muoHandle );
354 
356  iEvent.getByToken( trToken_, trHandle );
357 
359  iEvent.getByToken( phToken_, phHandle );
360 
361 
362  edm::Handle<edm::TriggerResults> handleTriggerTrigRes;
363 
364  edm::Handle<trigger::TriggerEvent> handleTriggerEvent;
365  edm::ESHandle<MagneticField> bFieldHandle;
366  // Filter out events if Trigger Filtering is requested
367  double PrescaleWeight = prescaleWeightProvider_->prescaleWeight( iEvent, iSetup );
368 
369  if (tnp_> 0) { // TnP method
370  if (den_genTriggerEventFlag_->on() && ! den_genTriggerEventFlag_->accept( iEvent, iSetup) ) return;
371  iEvent.getByToken( hltInputTag_, handleTriggerEvent);
372  if (handleTriggerEvent->sizeFilters()== 0) return;
373  const std::string & hltpath = hltpaths_num[0];
374  std::vector<reco::Muon> tagMuons;
375  for ( auto const & m : *muoHandle ) { // applying tag selection
376  if (false && !matchToTrigger(hltpath,m, handleTriggerEvent)) continue;
377  if ( muoSelection_ref( m ) ) tagMuons.push_back(m);
378  }
379  for (int i = 0; i<int(tagMuons.size());i++) {
380  for ( auto const & m : *muoHandle ) {
381  if (false && !matchToTrigger(hltpath,m, handleTriggerEvent)) continue;
382  if ((tagMuons[i].pt() == m.pt())) continue; // not the same
383  if ((tagMuons[i].p4()+m.p4()).M() >minmass_&& (tagMuons[i].p4()+m.p4()).M() <maxmass_) { // near to J/psi mass
384  muPhi_.denominator->Fill(m.phi());
385  muEta_.denominator->Fill(m.eta());
386  muPt_.denominator ->Fill(m.pt());
387  if (muoSelection_( m )) {
388  muPhi_.numerator->Fill(m.phi(),PrescaleWeight);
389  muEta_.numerator->Fill(m.eta(),PrescaleWeight);
390  muPt_.numerator ->Fill(m.pt(),PrescaleWeight);
391  }
392  }
393  }
394 
395  }
396 
397 
398  }
399  else { // reference method
400  if (den_genTriggerEventFlag_->on() && ! den_genTriggerEventFlag_->accept( iEvent, iSetup) ) return;
401  iEvent.getByToken( hltInputTag_, handleTriggerEvent);
402  if (handleTriggerEvent->sizeFilters()== 0) return;
403  const std::string & hltpath = hltpaths_den[0];
404  for (auto const & m : *muoHandle ) {
405  if (false && !matchToTrigger(hltpath,m, handleTriggerEvent)) continue;
406  if (!muoSelection_ref(m)) continue;
407  for (auto const & m1 : *muoHandle ) {
408  if (&m - &(*muoHandle)[0] >= &m1 - &(*muoHandle)[0]) continue;
409  //if (m1.pt() == m.pt()) continue; // probably not needed if using the above check
410  if (!muoSelection_ref(m1)) continue;
411  if (false && !matchToTrigger(hltpath,m1, handleTriggerEvent)) continue;
412  if (enum_ != 10) {
413  if (!DMSelection_ref(m1.p4() + m.p4())) continue;
414  if (m.charge()*m1.charge() > 0 ) continue;
415  }
416  iSetup.get<IdealMagneticFieldRecord>().get(bFieldHandle);
417  const reco::BeamSpot& vertexBeamSpot = *beamSpot;
418  std::vector<reco::TransientTrack> j_tks;
419  reco::TransientTrack mu1TT(m.track(), &(*bFieldHandle));
420  reco::TransientTrack mu2TT(m1.track(), &(*bFieldHandle));
421  j_tks.push_back(mu1TT);
422  j_tks.push_back(mu2TT);
423  KalmanVertexFitter jkvf;
424  TransientVertex jtv = jkvf.vertex(j_tks);
425  if (!jtv.isValid()) continue;
426  reco::Vertex jpsivertex = jtv;
427  float dimuonCL = 0;
428  if ( (jpsivertex.chi2() >= 0) && (jpsivertex.ndof() > 0) ) // I think these values are "unphysical"(no one will need to change them ever)so the can be fixed
429  dimuonCL = TMath::Prob(jpsivertex.chi2(), jpsivertex.ndof() );
430  math::XYZVector jpperp(m.px() + m1.px() ,
431  m.py() + m1.py() ,
432  0.);
433  GlobalPoint jVertex = jtv.position();
434  GlobalError jerr = jtv.positionError();
435  GlobalPoint displacementFromBeamspotJpsi( -1*((vertexBeamSpot.x0() - jVertex.x()) + (jVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
436  -1*((vertexBeamSpot.y0() - jVertex.y()) + (jVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
437  0);
438  reco::Vertex::Point vperpj(displacementFromBeamspotJpsi.x(), displacementFromBeamspotJpsi.y(), 0.);
439  float jpsi_cos = vperpj.Dot(jpperp) / (vperpj.R()*jpperp.R());
440  TrajectoryStateClosestToPoint mu1TS = mu1TT.impactPointTSCP();
441  TrajectoryStateClosestToPoint mu2TS = mu2TT.impactPointTSCP();
443  if (mu1TS.isValid() && mu2TS.isValid()) {
444  cApp.calculate(mu1TS.theState(), mu2TS.theState());
445  }
446  double DiMuMass = (m1.p4()+m.p4()).M();
447  switch(enum_) { // enum_ = 1...9, represents different sets of variables for different paths, we want to have different hists for different paths
448  case 1: tnp_=true; // already filled hists for tnp method
449  case 2:
450  if ((Jpsi_) && (!Upsilon_)) {
451  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
452  }
453  if ((!Jpsi_) && (Upsilon_)) {
454  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
455  }
456  if (dimuonCL < minprob) continue;
457  mu1Phi_.denominator->Fill(m.phi());
458  mu1Eta_.denominator->Fill(m.eta());
459  mu1Pt_.denominator ->Fill(m.pt());
460  mu2Phi_.denominator->Fill(m1.phi());
461  mu2Eta_.denominator->Fill(m1.eta());
462  mu2Pt_.denominator ->Fill(m1.pt());
463  DiMuPt_.denominator ->Fill((m1.p4()+m.p4()).Pt() );
464  DiMuEta_.denominator ->Fill((m1.p4()+m.p4()).Eta() );
465  DiMuPhi_.denominator ->Fill((m1.p4()+m.p4()).Phi());
466  break;
467  case 3:
468  if ((Jpsi_) && (!Upsilon_)) {
469  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
470  }
471  if ((!Jpsi_) && (Upsilon_)) {
472  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
473  }
474  if (dimuonCL < minprob) continue;
475  mu1Eta_.denominator->Fill(m.eta());
476  mu1Pt_.denominator ->Fill(m.pt());
477  mu2Eta_.denominator->Fill(m1.eta());
478  mu2Pt_.denominator ->Fill(m1.pt());
479  break;
480  case 4:
481  if (dimuonCL < minprob) continue;
482  DiMuMass_.denominator ->Fill(DiMuMass);
483  if ((Jpsi_) && (!Upsilon_)) {
484  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
485  }
486  if ((!Jpsi_) && (Upsilon_)) {
487  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
488  }
489  mu1Phi_.denominator->Fill(m.phi());
490  mu1Eta_.denominator->Fill(m.eta());
491  mu1Pt_.denominator ->Fill(m.pt());
492  mu2Phi_.denominator->Fill(m1.phi());
493  mu2Eta_.denominator->Fill(m1.eta());
494  mu2Pt_.denominator ->Fill(m1.pt());
495  DiMuPt_.denominator ->Fill((m1.p4()+m.p4()).Pt() );
496  DiMuEta_.denominator ->Fill((m1.p4()+m.p4()).Eta() );
497  DiMuPhi_.denominator ->Fill((m1.p4()+m.p4()).Phi());
499  break;
500  case 5:
501  if (dimuonCL < minprob) continue;
502  if ((Jpsi_) && (!Upsilon_)) {
503  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
504  }
505 
506  if ((!Jpsi_) && (Upsilon_)) {
507  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
508  }
509  mu1Phi_.denominator->Fill(m.phi());
510  mu1Eta_.denominator->Fill(m.eta());
511  mu1Pt_.denominator ->Fill(m.pt());
512  mu2Phi_.denominator->Fill(m1.phi());
513  mu2Eta_.denominator->Fill(m1.eta());
514  mu2Pt_.denominator ->Fill(m1.pt());
515  DiMuPt_.denominator ->Fill((m1.p4()+m.p4()).Pt() );
516  DiMuEta_.denominator ->Fill((m1.p4()+m.p4()).Eta() );
517  DiMuPhi_.denominator ->Fill((m1.p4()+m.p4()).Phi());
519  break;
520  case 6:
521  if (dimuonCL < minprob) continue;
522  if ((Jpsi_) && (!Upsilon_)) {
523  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
524  }
525  if ((!Jpsi_) && (Upsilon_)) {
526  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
527  }
528  for (auto const & m2 : *muoHandle) { // triple muon paths
529  if (false && !matchToTrigger(hltpath,m2, handleTriggerEvent)) continue;
530  if (m2.pt() == m.pt()) continue;
531  mu1Phi_.denominator->Fill(m.phi());
532  mu1Eta_.denominator->Fill(m.eta());
533  mu1Pt_.denominator ->Fill(m.pt());
534  mu2Phi_.denominator->Fill(m1.phi());
535  mu2Eta_.denominator->Fill(m1.eta());
536  mu2Pt_.denominator ->Fill(m1.pt());
537  mu3Phi_.denominator->Fill(m2.phi());
538  mu3Eta_.denominator->Fill(m2.eta());
539  mu3Pt_.denominator ->Fill(m2.pt());
540  }
541  break;
542 
543  case 7: // the hists for photon monitoring will be filled on 515 line
544  tnp_=false;
545  break;
546 
547  case 8: // vtx monitoring, filling probability, DS, DCA, cos of pointing angle to the PV, eta, pT of dimuon
548  if ((Jpsi_) && (!Upsilon_)) {
549  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
550  }
551 
552  if ((!Jpsi_) && (Upsilon_)) {
553  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
554  }
555 
556  DiMuProb_.denominator ->Fill( dimuonCL);
557  if (dimuonCL < minprob) continue;
558  DiMuDS_.denominator ->Fill( displacementFromBeamspotJpsi.perp() / sqrt(jerr.rerr(displacementFromBeamspotJpsi)));
559  DiMuPVcos_.denominator ->Fill(jpsi_cos );
560  DiMuPt_.denominator ->Fill((m1.p4()+m.p4()).Pt() );
561  DiMuEta_.denominator ->Fill((m1.p4()+m.p4()).Eta() );
562  DiMuDCA_.denominator ->Fill( cApp.distance());
563  break;
564  case 9:
565  if (dimuonCL < minprob) continue;
566  if (fabs(jpsi_cos) < mincos) continue;
567  if ((displacementFromBeamspotJpsi.perp() / sqrt(jerr.rerr(displacementFromBeamspotJpsi))) < minDS) continue;
568 
569  if (trHandle.isValid()) {
571 
572  for (auto const & t : *trHandle) {
573 
574  if (!trSelection_ref(t)) continue;
575  if (false && !matchToTrigger(hltpath,t, handleTriggerEvent)) continue;
576  const reco::Track& itrk1 = t ;
577 
578  if ((reco::deltaR(t,m1) <= min_dR)) continue; // checking overlapping
579  if ((reco::deltaR(t,m) <= min_dR)) continue;
580 
581  if (! itrk1.quality(reco::TrackBase::highPurity)) continue;
582 
584  double trackMass2 = kaon_mass * kaon_mass;
585  double MuMass2 = mu_mass * mu_mass; // 0.1056583745 *0.1056583745;
586  double e1 = sqrt(m.momentum().Mag2() + MuMass2 );
587  double e2 = sqrt(m1.momentum().Mag2() + MuMass2 );
588  double e3 = sqrt(itrk1.momentum().Mag2() + trackMass2 );
589 
590  p1 = reco::Particle::LorentzVector(m.px() , m.py() , m.pz() , e1 );
591  p2 = reco::Particle::LorentzVector(m1.px() , m1.py() , m1.pz() , e2 );
592  p3 = reco::Particle::LorentzVector(itrk1.px(), itrk1.py(), itrk1.pz(), e3 );
593  pB = p1 + p2 + p3;
594  if ( pB.mass() > maxmassJpsiTk || pB.mass() < minmassJpsiTk) continue;
595  reco::TransientTrack trTT(itrk1, &(*bFieldHandle));
596 
597  std::vector<reco::TransientTrack> t_tks;
598  t_tks.push_back(mu1TT);
599  t_tks.push_back(mu2TT);
600  t_tks.push_back(trTT);
601 
602  KalmanVertexFitter kvf;
603  TransientVertex tv = kvf.vertex(t_tks);
604  reco::Vertex vertex = tv;
605  if (!tv.isValid()) continue;
606  float JpsiTkCL = 0;
607  if ((vertex.chi2() >= 0.0) && (vertex.ndof() > 0) )
608  JpsiTkCL = TMath::Prob(vertex.chi2(), vertex.ndof() );
609  math::XYZVector pperp(m.px() + m1.px() + itrk1.px(),
610  m.py() + m1.py() + itrk1.py(),
611  0.);
612  GlobalPoint secondaryVertex = tv.position();
613  GlobalError err = tv.positionError();
614  GlobalPoint displacementFromBeamspot( -1*((vertexBeamSpot.x0() - secondaryVertex.x()) +
615  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
616  -1*((vertexBeamSpot.y0() - secondaryVertex.y()) +
617  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
618  0);
619  reco::Vertex::Point vperp(displacementFromBeamspot.x(),displacementFromBeamspot.y(),0.);
620  float jpsiKcos = vperp.Dot(pperp) / (vperp.R()*pperp.R());
621  if (JpsiTkCL < minprob) continue;
622  if (fabs(jpsiKcos) < mincos) continue;
623  if ((displacementFromBeamspot.perp() / sqrt(err.rerr(displacementFromBeamspot)))<minDS) continue;
624  muPhi_.denominator->Fill(t.phi());
625  muEta_.denominator->Fill(t.eta());
626  muPt_.denominator ->Fill(t.pt());
627 
629  }
630  }
631  break;
632  case 10:
633  if (trHandle.isValid()) {
634  for (auto const & t : *trHandle) {
635  if (!trSelection_ref(t)) continue;
636  if (false && !matchToTrigger(hltpath,t, handleTriggerEvent)) continue;
637  const reco::Track& itrk1 = t ;
638  if ((reco::deltaR(t,m1) <= min_dR)) continue; // checking overlapping
639  if ((reco::deltaR(t,m) <= min_dR)) continue;
640  if (! itrk1.quality(reco::TrackBase::highPurity)) continue;
642  double trackMass2 = kaon_mass * kaon_mass;
643  double MuMass2 = mu_mass * mu_mass; // 0.1056583745 *0.1056583745;
644  double e2 = sqrt(m1.momentum().Mag2() + MuMass2 );
645  double e3 = sqrt(itrk1.momentum().Mag2() + trackMass2 );
646  p2 = reco::Particle::LorentzVector(m1.px() , m1.py() , m1.pz() , e2 );
647  p3 = reco::Particle::LorentzVector(itrk1.px(), itrk1.py(), itrk1.pz(), e3 );
648  pB = p2 + p3;
649  if ( pB.mass() > maxmassJpsiTk || pB.mass()< minmassJpsiTk) continue;
650  reco::TransientTrack trTT(itrk1, &(*bFieldHandle));
651  std::vector<reco::TransientTrack> t_tks;
652  t_tks.push_back(mu2TT);
653  t_tks.push_back(trTT);
654  KalmanVertexFitter kvf;
655  TransientVertex tv = kvf.vertex(t_tks);
656  reco::Vertex vertex = tv;
657  if (!tv.isValid()) continue;
658  float JpsiTkCL = 0;
659  if ((vertex.chi2() >= 0.0) && (vertex.ndof() > 0) )
660  JpsiTkCL = TMath::Prob(vertex.chi2(), vertex.ndof() );
661  math::XYZVector pperp(m1.px() + itrk1.px(),
662  m1.py() + itrk1.py(),
663  0.);
664  GlobalPoint secondaryVertex = tv.position();
665  GlobalError err = tv.positionError();
666  GlobalPoint displacementFromBeamspot( -1*((vertexBeamSpot.x0() - secondaryVertex.x()) +
667  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
668  -1*((vertexBeamSpot.y0() - secondaryVertex.y()) +
669  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
670  0);
671  reco::Vertex::Point vperp(displacementFromBeamspot.x(),displacementFromBeamspot.y(),0.);
672  if (JpsiTkCL<minprob) continue;
673  muPhi_.denominator->Fill(m1.phi());
674  muEta_.denominator->Fill(m1.eta());
675  muPt_.denominator ->Fill(m1.pt());
676  }
677  }
678  break;
679  case 11:
680  case11_selection(dimuonCL, jpsi_cos, displacementFromBeamspotJpsi, jerr, trHandle, hltpath, handleTriggerEvent, m, m1, bFieldHandle, vertexBeamSpot, mu1Phi_.denominator, mu1Eta_.denominator, mu1Pt_.denominator, mu2Phi_.denominator, mu2Eta_.denominator, mu2Pt_.denominator);
681  break;
682  }
683  }
684  }
685 
686 
687  if (enum_ == 7) { // photons
688  const std::string & hltpath = hltpaths_den[0];
689  for (auto const & p : *phHandle) {
690  if (false && !matchToTrigger(hltpath,p, handleTriggerEvent)) continue;
691  phPhi_.denominator->Fill(p.phi());
692  phEta_.denominator->Fill(p.eta());
693  phPt_.denominator ->Fill(p.pt());
694  }
695 
696  }
697  //
699  // filling numerator hists
700  if (num_genTriggerEventFlag_->on() && ! num_genTriggerEventFlag_->accept( iEvent, iSetup) ) return;
701  iEvent.getByToken( hltInputTag_, handleTriggerEvent);
702  if (handleTriggerEvent->sizeFilters()== 0) return;
703  const std::string & hltpath1 = hltpaths_num[0];
704  for (auto const & m : *muoHandle ) {
705  if (false && !matchToTrigger(hltpath1,m, handleTriggerEvent)) continue;
706  if (!muoSelection_ref(m)) continue;
707  for (auto const & m1 : *muoHandle ) {
708  if (seagull_ && ((m.charge()* deltaPhi(m.phi(), m1.phi())) > 0.) ) continue;
709  if (m.charge()*m1.charge() > 0 ) continue;
710  if (m1.pt() == m.pt()) continue;
711  if (!muoSelection_ref(m1)) continue;
712  if (false && !matchToTrigger(hltpath1,m1, handleTriggerEvent)) continue;
713  if (!DMSelection_ref(m1.p4() + m.p4())) continue;
714  iSetup.get<IdealMagneticFieldRecord>().get(bFieldHandle);
715  const reco::BeamSpot& vertexBeamSpot = *beamSpot;
716  std::vector<reco::TransientTrack> j_tks;
717  reco::TransientTrack mu1TT(m.track(), &(*bFieldHandle));
718  reco::TransientTrack mu2TT(m1.track(), &(*bFieldHandle));
719  j_tks.push_back(mu1TT);
720  j_tks.push_back(mu2TT);
721  KalmanVertexFitter jkvf;
722  TransientVertex jtv = jkvf.vertex(j_tks);
723  if (!jtv.isValid()) continue;
724  reco::Vertex jpsivertex = jtv;
725  float dimuonCL = 0;
726  if ( (jpsivertex.chi2() >= 0) && (jpsivertex.ndof() > 0) ) // I think these values are "unphysical"(no one will need to change them ever)so the can be fixed
727  dimuonCL = TMath::Prob(jpsivertex.chi2(), jpsivertex.ndof() );
728  math::XYZVector jpperp(m.px() + m1.px() ,
729  m.py() + m1.py() ,
730  0.);
731  GlobalPoint jVertex = jtv.position();
732  GlobalError jerr = jtv.positionError();
733  GlobalPoint displacementFromBeamspotJpsi( -1*((vertexBeamSpot.x0() - jVertex.x()) + (jVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
734  -1*((vertexBeamSpot.y0() - jVertex.y()) + (jVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
735  0);
736  reco::Vertex::Point vperpj(displacementFromBeamspotJpsi.x(), displacementFromBeamspotJpsi.y(), 0.);
737  float jpsi_cos = vperpj.Dot(jpperp) / (vperpj.R()*jpperp.R());
738  TrajectoryStateClosestToPoint mu1TS = mu1TT.impactPointTSCP();
739  TrajectoryStateClosestToPoint mu2TS = mu2TT.impactPointTSCP();
741  if (mu1TS.isValid() && mu2TS.isValid()) {
742  cApp.calculate(mu1TS.theState(), mu2TS.theState());
743  }
744  double DiMuMass = (m1.p4()+m.p4()).M();
745  switch(enum_) { // enum_ = 1...9, represents different sets of variables for different paths, we want to have different hists for different paths
746  case 1: tnp_=true; // already filled hists for tnp method
747  case 2:
748  if ((Jpsi_) && (!Upsilon_)) {
749  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
750  }
751  if ((!Jpsi_) && (Upsilon_)) {
752  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
753  }
754  if (dimuonCL < minprob) continue;
755  mu1Phi_.numerator->Fill(m.phi(),PrescaleWeight);
756  mu1Eta_.numerator->Fill(m.eta(),PrescaleWeight);
757  mu1Pt_.numerator ->Fill(m.pt(),PrescaleWeight);
758  mu2Phi_.numerator->Fill(m1.phi(),PrescaleWeight);
759  mu2Eta_.numerator->Fill(m1.eta(),PrescaleWeight);
760  mu2Pt_.numerator ->Fill(m1.pt(),PrescaleWeight);
761  DiMuPt_.numerator ->Fill((m1.p4()+m.p4()).Pt() ,PrescaleWeight);
762  DiMuEta_.numerator ->Fill((m1.p4()+m.p4()).Eta() ,PrescaleWeight);
763  DiMuPhi_.numerator ->Fill((m1.p4()+m.p4()).Phi(),PrescaleWeight);
764  break;
765  case 3:
766  if ((Jpsi_) && (!Upsilon_)) {
767  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
768  }
769 
770  if ((!Jpsi_) && (Upsilon_)) {
771  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
772  }
773  if (dimuonCL < minprob) continue;
774  mu1Eta_.numerator->Fill(m.eta(),PrescaleWeight);
775  mu1Pt_.numerator ->Fill(m.pt(),PrescaleWeight);
776  mu2Eta_.numerator->Fill(m1.eta(),PrescaleWeight);
777  mu2Pt_.numerator ->Fill(m1.pt(),PrescaleWeight);
778  break;
779  case 4:
780  if (dimuonCL < minprob) continue;
781  DiMuMass_.numerator ->Fill(DiMuMass);
782  if ((Jpsi_) && (!Upsilon_)) {
783  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
784  }
785  if ((!Jpsi_) && (Upsilon_)) {
786  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
787  }
788  mu1Phi_.numerator->Fill(m.phi(),PrescaleWeight);
789  mu1Eta_.numerator->Fill(m.eta(),PrescaleWeight);
790  mu1Pt_.numerator ->Fill(m.pt(),PrescaleWeight);
791  mu2Phi_.numerator->Fill(m1.phi(),PrescaleWeight);
792  mu2Eta_.numerator->Fill(m1.eta(),PrescaleWeight);
793  mu2Pt_.numerator ->Fill(m1.pt(),PrescaleWeight);
794  DiMuPt_.numerator ->Fill((m1.p4()+m.p4()).Pt() ,PrescaleWeight);
795  DiMuEta_.numerator ->Fill((m1.p4()+m.p4()).Eta() ,PrescaleWeight);
796  DiMuPhi_.numerator ->Fill((m1.p4()+m.p4()).Phi(),PrescaleWeight);
797  DiMudR_.numerator ->Fill(reco::deltaR(m,m1),PrescaleWeight);
798  break;
799  case 5:
800  if (dimuonCL < minprob) continue;
801  if ((Jpsi_) && (!Upsilon_)) {
802  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
803  }
804  if ((!Jpsi_) && (Upsilon_)) {
805  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
806  }
807  mu1Phi_.numerator->Fill(m.phi(),PrescaleWeight);
808  mu1Eta_.numerator->Fill(m.eta(),PrescaleWeight);
809  mu1Pt_.numerator ->Fill(m.pt(),PrescaleWeight);
810  mu2Phi_.numerator->Fill(m1.phi(),PrescaleWeight);
811  mu2Eta_.numerator->Fill(m1.eta(),PrescaleWeight);
812  mu2Pt_.numerator ->Fill(m1.pt(),PrescaleWeight);
813  DiMuPt_.numerator ->Fill((m1.p4()+m.p4()).Pt() ,PrescaleWeight);
814  DiMuEta_.numerator ->Fill((m1.p4()+m.p4()).Eta() ,PrescaleWeight);
815  DiMuPhi_.numerator ->Fill((m1.p4()+m.p4()).Phi(),PrescaleWeight);
816  DiMudR_.numerator ->Fill(reco::deltaR(m,m1),PrescaleWeight);
817  break;
818  case 6:
819  if (dimuonCL < minprob) continue;
820  if ((Jpsi_) && (!Upsilon_)) {
821  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
822  }
823  if ((!Jpsi_) && (Upsilon_)) {
824  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
825  }
826  for (auto const & m2 : *muoHandle) { // triple muon paths
827  if (false && !matchToTrigger(hltpath1,m2, handleTriggerEvent)) continue;
828  if (m2.pt() == m.pt()) continue;
829  mu1Phi_.numerator->Fill(m.phi(),PrescaleWeight);
830  mu1Eta_.numerator->Fill(m.eta(),PrescaleWeight);
831  mu1Pt_.numerator ->Fill(m.pt(),PrescaleWeight);
832  mu2Phi_.numerator->Fill(m1.phi(),PrescaleWeight);
833  mu2Eta_.numerator->Fill(m1.eta(),PrescaleWeight);
834  mu2Pt_.numerator ->Fill(m1.pt(),PrescaleWeight);
835  mu3Phi_.numerator->Fill(m2.phi(),PrescaleWeight);
836  mu3Eta_.numerator->Fill(m2.eta(),PrescaleWeight);
837  mu3Pt_.numerator ->Fill(m2.pt(),PrescaleWeight);
838  }
839  break;
840  case 7: // the hists for photon monitoring will be filled on 515 line
841  tnp_=false;
842  break;
843  case 8: // vtx monitoring, filling probability, DS, DCA, cos of pointing angle to the PV, eta, pT of dimuon
844  if ((Jpsi_) && (!Upsilon_)) {
845  if (DiMuMass> maxmassJpsi || DiMuMass< minmassJpsi) continue;
846  }
847  if ((!Jpsi_) && (Upsilon_)) {
848  if (DiMuMass> maxmassUpsilon || DiMuMass< minmassUpsilon) continue;
849  }
850  DiMuProb_.numerator ->Fill( dimuonCL,PrescaleWeight);
851  if (dimuonCL < minprob) continue;
852  DiMuDS_.numerator ->Fill( displacementFromBeamspotJpsi.perp() / sqrt(jerr.rerr(displacementFromBeamspotJpsi)),PrescaleWeight);
853  DiMuPVcos_.numerator ->Fill(jpsi_cos ,PrescaleWeight);
854  DiMuPt_.numerator ->Fill((m1.p4()+m.p4()).Pt() ,PrescaleWeight);
855  DiMuEta_.numerator ->Fill((m1.p4()+m.p4()).Eta() ,PrescaleWeight);
856  DiMuDCA_.numerator ->Fill( cApp.distance(),PrescaleWeight);
857  break;
858  case 9:
859  if (dimuonCL < minprob) continue;
860  if (fabs(jpsi_cos) < mincos) continue;
861  if ((displacementFromBeamspotJpsi.perp() / sqrt(jerr.rerr(displacementFromBeamspotJpsi))) < minDS) continue;
862  if (trHandle.isValid()) {
863  for (auto const & t : *trHandle) {
864  if (!trSelection_ref(t)) continue;
865  if (false && !matchToTrigger(hltpath1,t, handleTriggerEvent)) continue;
866  const reco::Track& itrk1 = t ;
867  if ((reco::deltaR(t,m1) <= min_dR)) continue; // checking overlapping
868  if ((reco::deltaR(t,m) <= min_dR)) continue;
869  if (! itrk1.quality(reco::TrackBase::highPurity)) continue;
871  double trackMass2 = kaon_mass * kaon_mass;
872  double MuMass2 = mu_mass * mu_mass; // 0.1056583745 *0.1056583745;
873  double e1 = sqrt(m.momentum().Mag2() + MuMass2 );
874  double e2 = sqrt(m1.momentum().Mag2() + MuMass2 );
875  double e3 = sqrt(itrk1.momentum().Mag2() + trackMass2 );
876  p1 = reco::Particle::LorentzVector(m.px() , m.py() , m.pz() , e1 );
877  p2 = reco::Particle::LorentzVector(m1.px() , m1.py() , m1.pz() , e2 );
878  p3 = reco::Particle::LorentzVector(itrk1.px(), itrk1.py(), itrk1.pz(), e3 );
879  pB = p1 + p2 + p3;
880  if ( pB.mass() > maxmassJpsiTk || pB.mass()< minmassJpsiTk) continue;
881  reco::TransientTrack trTT(itrk1, &(*bFieldHandle));
882  std::vector<reco::TransientTrack> t_tks;
883  t_tks.push_back(mu1TT);
884  t_tks.push_back(mu2TT);
885  t_tks.push_back(trTT);
886  KalmanVertexFitter kvf;
887  TransientVertex tv = kvf.vertex(t_tks);
888  reco::Vertex vertex = tv;
889  if (!tv.isValid()) continue;
890  float JpsiTkCL = 0;
891  if ((vertex.chi2() >= 0.0) && (vertex.ndof() > 0) )
892  JpsiTkCL = TMath::Prob(vertex.chi2(), vertex.ndof() );
893  math::XYZVector pperp(m.px() + m1.px() + itrk1.px(),
894  m.py() + m1.py() + itrk1.py(),
895  0.);
896  GlobalPoint secondaryVertex = tv.position();
897  GlobalError err = tv.positionError();
898  GlobalPoint displacementFromBeamspot( -1*((vertexBeamSpot.x0() - secondaryVertex.x()) +
899  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
900  -1*((vertexBeamSpot.y0() - secondaryVertex.y()) +
901  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
902  0);
903  reco::Vertex::Point vperp(displacementFromBeamspot.x(),displacementFromBeamspot.y(),0.);
904  float jpsiKcos = vperp.Dot(pperp) / (vperp.R()*pperp.R());
905  if (JpsiTkCL<minprob) continue;
906  if (fabs(jpsiKcos)<mincos) continue;
907  if ((displacementFromBeamspot.perp() / sqrt(err.rerr(displacementFromBeamspot)))<minDS) continue;
908  muPhi_.numerator->Fill(t.phi(),PrescaleWeight);
909  muEta_.numerator->Fill(t.eta(),PrescaleWeight);
910  muPt_.numerator ->Fill(t.pt(),PrescaleWeight);
911  }
912  }
913  break;
914 
915  case 10:
916  if (trHandle.isValid()) {
917  for (auto const & t : *trHandle) {
918  if (!trSelection_ref(t)) continue;
919  if (false && !matchToTrigger(hltpath1,t, handleTriggerEvent)) continue;
920  const reco::Track& itrk1 = t ;
921  if ((reco::deltaR(t,m1) <= min_dR)) continue; // checking overlapping
922  if ((reco::deltaR(t,m) <= min_dR)) continue;
923  if (! itrk1.quality(reco::TrackBase::highPurity)) continue;
925  double trackMass2 = kaon_mass * kaon_mass;
926  double MuMass2 = mu_mass * mu_mass; // 0.1056583745 *0.1056583745;
927  double e2 = sqrt(m1.momentum().Mag2() + MuMass2 );
928  double e3 = sqrt(itrk1.momentum().Mag2() + trackMass2 );
929  p2 = reco::Particle::LorentzVector(m1.px() , m1.py() , m1.pz() , e2 );
930  p3 = reco::Particle::LorentzVector(itrk1.px(), itrk1.py(), itrk1.pz(), e3 );
931  pB = p2 + p3;
932  if ( pB.mass() > maxmassJpsiTk || pB.mass()< minmassJpsiTk) continue;
933  reco::TransientTrack trTT(itrk1, &(*bFieldHandle));
934  std::vector<reco::TransientTrack> t_tks;
935  t_tks.push_back(mu2TT);
936  t_tks.push_back(trTT);
937  if (t_tks.size()!=2) continue;
938  KalmanVertexFitter kvf;
939  TransientVertex tv = kvf.vertex(t_tks);
940  reco::Vertex vertex = tv;
941  if (!tv.isValid()) continue;
942  float JpsiTkCL = 0;
943  if ((vertex.chi2() >= 0.0) && (vertex.ndof() > 0) )
944  JpsiTkCL = TMath::Prob(vertex.chi2(), vertex.ndof() );
945  math::XYZVector pperp(m1.px() + itrk1.px(),
946  m1.py() + itrk1.py(),
947  0.);
948  GlobalPoint secondaryVertex = tv.position();
949  GlobalError err = tv.positionError();
950  GlobalPoint displacementFromBeamspot( -1*((vertexBeamSpot.x0() - secondaryVertex.x()) +
951  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
952  -1*((vertexBeamSpot.y0() - secondaryVertex.y()) +
953  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
954  0);
955  reco::Vertex::Point vperp(displacementFromBeamspot.x(),displacementFromBeamspot.y(),0.);
956  if (JpsiTkCL<minprob) continue;
957  muPhi_.numerator->Fill(m1.phi(),PrescaleWeight);
958  muEta_.numerator->Fill(m1.eta(),PrescaleWeight);
959  muPt_.numerator ->Fill(m1.pt(),PrescaleWeight);
960  }
961  }
962  break;
963  case 11:
964  case11_selection(dimuonCL, jpsi_cos, displacementFromBeamspotJpsi, jerr, trHandle, hltpath, handleTriggerEvent, m, m1, bFieldHandle, vertexBeamSpot, mu1Phi_.numerator, mu1Eta_.numerator, mu1Pt_.numerator, mu2Phi_.numerator, mu2Eta_.numerator, mu2Pt_.numerator);
965  break;
966  }
967  }
968  }
969  if (enum_ == 7) { // photons
970  const std::string &hltpath = hltpaths_num[0];
971  for (auto const & p : *phHandle) {
972  if (false && !matchToTrigger(hltpath,p, handleTriggerEvent)) continue;
973  phPhi_.numerator->Fill(p.phi(),PrescaleWeight);
974  phEta_.numerator->Fill(p.eta(),PrescaleWeight);
975  phPt_.numerator ->Fill(p.pt(),PrescaleWeight);
976  }
977  }
978  }
979 }
980 
981 
982 
983 
985 {
986  pset.addNode((edm::ParameterDescription<int>("nbins", true) and
987  edm::ParameterDescription<double>("xmin", true) and
988  edm::ParameterDescription<double>("xmax", true)) xor
989  edm::ParameterDescription<std::vector<double>>("edges", true));
990 }
991 
993 {
994  pset.add<int> ( "nbins", 2500);
995 }
996 
998 {
1000  desc.add<std::string> ( "FolderName", "HLT/BPH/" );
1001  desc.add<edm::InputTag>( "tracks", edm::InputTag("generalTracks") );
1002  desc.add<edm::InputTag>( "photons", edm::InputTag("photons") );
1003  desc.add<edm::InputTag>( "offlinePVs", edm::InputTag("offlinePrimaryVertices") );
1004  desc.add<edm::InputTag>( "beamSpot",edm::InputTag("offlineBeamSpot") );
1005  desc.add<edm::InputTag>( "muons", edm::InputTag("muons") );
1006  desc.add<edm::InputTag>( "hltTriggerSummaryAOD", edm::InputTag("hltTriggerSummaryAOD","","HLT") );
1007  desc.add<std::string>("muoSelection", "abs(eta)<1.4 & isPFMuon & isGlobalMuon & innerTrack.hitPattern.trackerLayersWithMeasurement>5 & innerTrack.hitPattern.numberOfValidPixelHits> 0");
1008  desc.add<std::string>("muoSelection_ref", "isPFMuon & isGlobalMuon & innerTrack.hitPattern.trackerLayersWithMeasurement>5 & innerTrack.hitPattern.numberOfValidPixelHits> 0");
1009  desc.add<std::string>("muoSelection_tag", "isGlobalMuon && isPFMuon && isTrackerMuon && abs(eta) < 2.4 && innerTrack.hitPattern.numberOfValidPixelHits > 0 && innerTrack.hitPattern.trackerLayersWithMeasurement > 5 && globalTrack.hitPattern.numberOfValidMuonHits > 0 && globalTrack.normalizedChi2 < 10"); // tight selection for tag muon
1010  desc.add<std::string>("muoSelection_probe", "isPFMuon & isGlobalMuon & innerTrack.hitPattern.trackerLayersWithMeasurement>5 & innerTrack.hitPattern.numberOfValidPixelHits> 0");
1011  desc.add<std::string>("trSelection_ref", "");
1012  desc.add<std::string>("DMSelection_ref", "Pt>4 & abs(eta)");
1013 
1014  desc.add<int>("nmuons", 1);
1015  desc.add<bool>( "tnp", false );
1016  desc.add<int>( "L3", 0 );
1017  desc.add<int>( "trOrMu", 0 ); // if =0, track param monitoring
1018  desc.add<int>( "Jpsi", 0 );
1019  desc.add<int>( "Upsilon", 0 );
1020  desc.add<int>( "enum", 1 ); // 1...9, 9 sets of variables to be filled, depends on the hlt path
1021  desc.add<int>( "seagull", 1 ); // 1...9, 9 sets of variables to be filled, depends on the hlt path
1022  desc.add<double>( "maxmass", 3.596 );
1023  desc.add<double>( "minmass", 2.596 );
1024  desc.add<double>( "maxmassJpsi", 3.2 );
1025  desc.add<double>( "minmassJpsi", 3. );
1026  desc.add<double>( "maxmassUpsilon", 8.1 );
1027  desc.add<double>( "minmassUpsilon", 8. );
1028  desc.add<double>( "maxmassTkTk", 10);
1029  desc.add<double>( "minmassTkTk", 0);
1030  desc.add<double>( "maxmassJpsiTk", 5.46 );
1031  desc.add<double>( "minmassJpsiTk", 5.1 );
1032  desc.add<double>( "kaon_mass", 0.493677 );
1033  desc.add<double>( "mu_mass", 0.1056583745);
1034  desc.add<double>( "min_dR", 0.001);
1035  desc.add<double>( "minprob", 0.005 );
1036  desc.add<double>( "mincos", 0.95 );
1037  desc.add<double>( "minDS", 3. );
1038 
1039  edm::ParameterSetDescription genericTriggerEventPSet;
1040  genericTriggerEventPSet.add<bool>("andOr");
1041  genericTriggerEventPSet.add<edm::InputTag>("dcsInputTag", edm::InputTag("scalersRawToDigi") );
1042  genericTriggerEventPSet.add<edm::InputTag>("hltInputTag", edm::InputTag("TriggerResults::HLT") );
1043  genericTriggerEventPSet.add<std::vector<int> >("dcsPartitions",{});
1044  genericTriggerEventPSet.add<bool>("andOrDcs", false);
1045  genericTriggerEventPSet.add<bool>("errorReplyDcs", true);
1046  genericTriggerEventPSet.add<std::string>("dbLabel","");
1047  genericTriggerEventPSet.add<bool>("andOrHlt", true);
1048  genericTriggerEventPSet.add<bool>("andOrL1", true);
1049  genericTriggerEventPSet.add<std::vector<std::string> >("hltPaths",{});
1050  genericTriggerEventPSet.add<std::vector<std::string> >("l1Algorithms",{});
1051  genericTriggerEventPSet.add<std::string>("hltDBKey","");
1052  genericTriggerEventPSet.add<bool>("errorReplyHlt",false);
1053  genericTriggerEventPSet.add<bool>("errorReplyL1",true);
1054  genericTriggerEventPSet.add<bool>("l1BeforeMask",true);
1055  genericTriggerEventPSet.add<unsigned int>("verbosityLevel",0);
1056  desc.add<edm::ParameterSetDescription>("numGenericTriggerEventPSet", genericTriggerEventPSet);
1057  desc.add<edm::ParameterSetDescription>("denGenericTriggerEventPSet", genericTriggerEventPSet);
1058 
1059  edm::ParameterSetDescription PrescaleTriggerEventPSet;
1060  PrescaleTriggerEventPSet.add<unsigned int>("prescaleWeightVerbosityLevel",0);
1061  PrescaleTriggerEventPSet.add<edm::InputTag>("prescaleWeightTriggerResults",edm::InputTag("TriggerResults::HLT"));
1062  PrescaleTriggerEventPSet.add<edm::InputTag>("prescaleWeightL1GtTriggerMenuLite",edm::InputTag("l1GtTriggerMenuLite"));
1063  PrescaleTriggerEventPSet.add<std::vector<std::string>>("prescaleWeightHltPaths",{});
1064  desc.add<edm::ParameterSetDescription>("PrescaleTriggerEventPSet", PrescaleTriggerEventPSet);
1065 
1066  edm::ParameterSetDescription histoPSet;
1078  fillHistoPSetDescription(phiPSet);
1079  fillHistoPSetDescription(ptPSet);
1080  fillHistoPSetDescription(etaPSet);
1081  fillHistoPSetDescription(z0PSet);
1082  fillHistoPSetDescription(d0PSet);
1083  fillHistoPSetDescription(dRPSet);
1084  fillHistoPSetDescription(massPSet);
1085  fillHistoPSetDescription(dcaPSet);
1086  fillHistoPSetDescription(dsPSet);
1087  fillHistoPSetDescription(cosPSet);
1088  fillHistoPSetDescription(probPSet);
1089  histoPSet.add<edm::ParameterSetDescription>("d0PSet", d0PSet);
1090  histoPSet.add<edm::ParameterSetDescription>("etaPSet", etaPSet);
1091  histoPSet.add<edm::ParameterSetDescription>("phiPSet", phiPSet);
1092  histoPSet.add<edm::ParameterSetDescription>("ptPSet", ptPSet);
1093  histoPSet.add<edm::ParameterSetDescription>("z0PSet", z0PSet);
1094  histoPSet.add<edm::ParameterSetDescription>("dRPSet", dRPSet);
1095  histoPSet.add<edm::ParameterSetDescription>("massPSet", massPSet);
1096  histoPSet.add<edm::ParameterSetDescription>("dcaPSet", dcaPSet);
1097  histoPSet.add<edm::ParameterSetDescription>("dsPSet", dsPSet);
1098  histoPSet.add<edm::ParameterSetDescription>("cosPSet", cosPSet);
1099  histoPSet.add<edm::ParameterSetDescription>("probPSet", probPSet);
1100  desc.add<edm::ParameterSetDescription>("histoPSet",histoPSet);
1101 
1102  descriptions.add("bphMonitoring", desc);
1103 }
1104 template <typename T>
1105 bool BPHMonitor::matchToTrigger(const std::string &theTriggerName , T t, edm::Handle<trigger::TriggerEvent> handleTriggerEvent) {
1106  //bool BPHMonitor::matchToTrigger(std::string theTriggerName,T t, edm::Handle<trigger::TriggerEventWithRefs> handleTriggerEvent) {
1107 
1108  bool matchedToTrigger = false;
1109  if (handleTriggerEvent->sizeFilters() > 0) {
1110  const trigger::TriggerObjectCollection & toc(handleTriggerEvent->getObjects()); //Handle< trigger::TriggerEvent > handleTriggerEvent;
1111  for ( size_t ia = 0; ia < handleTriggerEvent->sizeFilters(); ++ ia) {
1112  std::string fullname = handleTriggerEvent->filterTag(ia).encode();
1113  std::string name;
1114  size_t p = fullname.find_first_of(':');
1115  if ( p != std::string::npos) {name = fullname.substr(0, p);}
1116  else {name = fullname;}
1117  const trigger::Keys & k = handleTriggerEvent->filterKeys(ia);
1118  for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
1119  reco::Particle theTriggerParticle = toc[*ki].particle();
1120  if (name.find(theTriggerName) != string::npos) {
1121  if ((reco::deltaR(t.eta(), t.phi(),theTriggerParticle.eta(),theTriggerParticle.phi()) <= 0.2)) {
1122  matchedToTrigger = true;
1123  }
1124  }
1125  }
1126  }
1127 
1128  return matchedToTrigger;
1129  }
1130  else {cout <<theTriggerName <<"\t\tNo HLT filters" <<endl; return false;}
1131 }
1132 
1133 void BPHMonitor::case11_selection(const float & dimuonCL, const float & jpsi_cos, const GlobalPoint & displacementFromBeamspotJpsi, const GlobalError & jerr, const edm::Handle<reco::TrackCollection> & trHandle, const std::string & hltpath, const edm::Handle<trigger::TriggerEvent> & handleTriggerEvent, const reco::Muon& m, const reco::Muon& m1, const edm::ESHandle<MagneticField> & bFieldHandle, const reco::BeamSpot & vertexBeamSpot, MonitorElement* phi1, MonitorElement* eta1, MonitorElement* pT1, MonitorElement* phi2, MonitorElement* eta2, MonitorElement* pT2) {
1134  //cout <<"\nInside case11_selection" <<endl;
1135  if (dimuonCL < minprob) return;
1136  if (fabs(jpsi_cos) < mincos) return;
1137  if ((displacementFromBeamspotJpsi.perp() / sqrt(jerr.rerr(displacementFromBeamspotJpsi))) < minDS) return;
1138  if (trHandle.isValid()) {
1140  for (auto const & t : *trHandle) {
1141  if (!trSelection_ref(t)) continue;
1142  if (false && !matchToTrigger(hltpath,t, handleTriggerEvent)) continue;
1143  if ((reco::deltaR(t,m) <= min_dR)) continue; // checking overlapping
1144  if ((reco::deltaR(t,m1) <= min_dR)) continue; // checking overlapping
1145  for (auto const & t1 : *trHandle) {
1146  if (&t - &(*trHandle)[0] >= &t1 - &(*trHandle)[0]) continue; // not enough, need the following DeltaR checks
1147  //if (t.pt() == t1.pt()) continue;
1148  if (!trSelection_ref(t1)) continue;
1149  if (false && !matchToTrigger(hltpath,t1, handleTriggerEvent)) continue;
1150  if ((reco::deltaR(t1,m) <= min_dR)) continue; // checking overlapping
1151  if ((reco::deltaR(t1,m1) <= min_dR)) continue; // checking overlapping
1152  if ((reco::deltaR(t,t1) <= min_dR)) continue; // checking overlapping
1153  const reco::Track& itrk1 = t ;
1154  const reco::Track& itrk2 = t1 ;
1155  if (! itrk1.quality(reco::TrackBase::highPurity)) continue;
1156  if (! itrk2.quality(reco::TrackBase::highPurity)) continue;
1157  reco::Particle::LorentzVector pB, pTkTk, p1, p2, p3, p4;
1158  double trackMass2 = kaon_mass * kaon_mass;
1159  double MuMass2 = mu_mass * mu_mass; // 0.1056583745 *0.1056583745;
1160  double e1 = sqrt(m.momentum().Mag2() + MuMass2 );
1161  double e2 = sqrt(m1.momentum().Mag2() + MuMass2 );
1162  double e3 = sqrt(itrk1.momentum().Mag2() + trackMass2 );
1163  double e4 = sqrt(itrk2.momentum().Mag2() + trackMass2 );
1164  p1 = reco::Particle::LorentzVector(m.px() , m.py() , m.pz() , e1 );
1165  p2 = reco::Particle::LorentzVector(m1.px() , m1.py() , m1.pz() , e2 );
1166  p3 = reco::Particle::LorentzVector(itrk1.px(), itrk1.py(), itrk1.pz(), e3 );
1167  p4 = reco::Particle::LorentzVector(itrk2.px(), itrk2.py(), itrk2.pz(), e4 );
1168  pTkTk = p3 + p4;
1169  if (pTkTk.mass() > maxmassTkTk || pTkTk.mass() < minmassTkTk) continue;
1170  pB = p1 + p2 + p3 + p4;
1171  if ( pB.mass() > maxmassJpsiTk || pB.mass()< minmassJpsiTk) continue;
1172  reco::TransientTrack mu1TT(m.track(), &(*bFieldHandle));
1173  reco::TransientTrack mu2TT(m1.track(), &(*bFieldHandle));
1174  reco::TransientTrack trTT(itrk1, &(*bFieldHandle));
1175  reco::TransientTrack tr1TT(itrk2, &(*bFieldHandle));
1176  std::vector<reco::TransientTrack> t_tks;
1177  t_tks.push_back(mu1TT);
1178  t_tks.push_back(mu2TT);
1179  t_tks.push_back(trTT);
1180  t_tks.push_back(tr1TT);
1181  KalmanVertexFitter kvf;
1182  TransientVertex tv = kvf.vertex(t_tks); // this will compare the tracks
1183  reco::Vertex vertex = tv;
1184  if (!tv.isValid()) continue;
1185  float JpsiTkCL = 0;
1186  if ((vertex.chi2() >= 0.0) && (vertex.ndof() > 0) )
1187  JpsiTkCL = TMath::Prob(vertex.chi2(), vertex.ndof() );
1188  math::XYZVector pperp(m.px() + m1.px() + itrk1.px() + itrk2.px(),
1189  m.py() + m1.py() + itrk1.py() + itrk2.py(),
1190  0.);
1191  GlobalPoint secondaryVertex = tv.position();
1192  GlobalError err = tv.positionError();
1193  GlobalPoint displacementFromBeamspot( -1*((vertexBeamSpot.x0() - secondaryVertex.x()) +
1194  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
1195  -1*((vertexBeamSpot.y0() - secondaryVertex.y()) +
1196  (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
1197  0);
1198  reco::Vertex::Point vperp(displacementFromBeamspot.x(),displacementFromBeamspot.y(),0.);
1199  float jpsiKcos = vperp.Dot(pperp) / (vperp.R()*pperp.R());
1200  if (JpsiTkCL < minprob) continue;
1201  if (fabs(jpsiKcos) < mincos) continue;
1202  if ((displacementFromBeamspot.perp() / sqrt(err.rerr(displacementFromBeamspot))) < minDS) continue;
1203 
1204  phi1->Fill(t.phi());
1205  eta1->Fill(t.eta());
1206  pT1->Fill(t.pt());
1207  phi2->Fill(t1.phi());
1208  eta2->Fill(t1.eta());
1209  pT2->Fill(t1.pt());
1210  } // for (auto const & t1 : *trHandle)
1211  } // for (auto const & t : *trHandle)
1212  } // if (trHandle.isValid())
1213 }
1214 
1215 // Define this as a plug-in
GlobalError positionError() const
MEbinning z0_binning_
Definition: BPHMonitor.h:114
int seagull_
Definition: BPHMonitor.h:176
MonitorElement * numerator
Definition: BPHMonitor.h:60
T getParameter(std::string const &) const
double z0() const
z coordinate
Definition: BeamSpot.h:68
METME mu2Phi_
Definition: BPHMonitor.h:133
int Upsilon_
Definition: BPHMonitor.h:174
METME DiMuDCA_
Definition: BPHMonitor.h:155
float distance() const override
METME phEta_
Definition: BPHMonitor.h:147
TrackRef track() const override
reference to a Track
Definition: Muon.h:49
METME phPt_
Definition: BPHMonitor.h:148
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:186
StringCutObjectSelector< reco::Candidate::LorentzVector, true > DMSelection_ref
Definition: BPHMonitor.h:200
T perp() const
Definition: PV3DBase.h:72
double maxmassJpsiTk
Definition: BPHMonitor.h:185
double kaon_mass
Definition: BPHMonitor.h:187
trigger::size_type sizeFilters() const
Definition: TriggerEvent.h:135
METME mu2Pt_
Definition: BPHMonitor.h:135
METME muPt_
Definition: BPHMonitor.h:124
int prescaleWeight(const edm::Event &event, const edm::EventSetup &setup)
METME mu1Eta_
Definition: BPHMonitor.h:129
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:160
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:519
METME phPhi_
Definition: BPHMonitor.h:146
const FreeTrajectoryState & theState() const
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
edm::EDGetTokenT< reco::MuonCollection > muoToken_
Definition: BPHMonitor.h:105
void bookME(DQMStore::IBooker &, METME &me, std::string &histname, std::string &histtitle, int &nbins, double &xmin, double &xmax)
Definition: BPHMonitor.cc:171
double px() const final
x coordinate of momentum vector
std::string folderName_
Definition: BPHMonitor.h:102
const Keys & filterKeys(trigger::size_type index) const
Definition: TriggerEvent.h:111
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
T y() const
Definition: PV3DBase.h:63
static void fillHistoLSPSetDescription(edm::ParameterSetDescription &pset)
Definition: BPHMonitor.cc:992
METME mu2Eta_
Definition: BPHMonitor.h:134
ParameterDescriptionNode * addNode(ParameterDescriptionNode const &node)
METME DiMuPhi_
Definition: BPHMonitor.h:149
METME mu3Pt_
Definition: BPHMonitor.h:140
Provides a code based selection for trigger and DCS information in order to have no failing filters i...
PrescaleWeightProvider * prescaleWeightProvider_
Definition: BPHMonitor.h:164
MEbinning eta_binning_
Definition: BPHMonitor.h:112
bool matchToTrigger(const std::string &theTriggerName, T t, edm::Handle< trigger::TriggerEvent > handleTriggerEvent)
Definition: BPHMonitor.cc:1105
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:627
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:675
edm::EDGetTokenT< trigger::TriggerEvent > hltInputTag_
Definition: BPHMonitor.h:195
double maxmassTkTk
Definition: BPHMonitor.h:183
std::string encode() const
Definition: InputTag.cc:166
edm::EDGetTokenT< reco::TrackCollection > trToken_
Definition: BPHMonitor.h:107
double minmass_
Definition: BPHMonitor.h:178
std::vector< Muon > MuonCollection
collection of Muon objects
Definition: MuonFwd.h:9
double minprob
Definition: BPHMonitor.h:191
void Fill(long long x)
METME mu1Phi_
Definition: BPHMonitor.h:128
METME DiMuEta_
Definition: BPHMonitor.h:150
double dydz() const
dydz slope
Definition: BeamSpot.h:84
MEbinning dca_binning_
Definition: BPHMonitor.h:117
Vector momentum() const final
spatial momentum vector
int iEvent
Definition: GenABIO.cc:230
void setMETitle(METME &me, std::string titleX, std::string titleY)
Definition: BPHMonitor.cc:162
MEbinning pt_binning_
Definition: BPHMonitor.h:111
double phi() const
momentum azimuthal angle
Definition: Particle.h:106
double mincos
Definition: BPHMonitor.h:192
double maxmassJpsi
Definition: BPHMonitor.h:179
const TriggerObjectCollection & getObjects() const
Definition: TriggerEvent.h:98
void case11_selection(const float &dimuonCL, const float &jpsi_cos, const GlobalPoint &displacementFromBeamspotJpsi, const GlobalError &jerr, const edm::Handle< reco::TrackCollection > &trHandle, const std::string &hltpath, const edm::Handle< trigger::TriggerEvent > &handleTriggerEvent, const reco::Muon &m, const reco::Muon &m1, const edm::ESHandle< MagneticField > &bFieldHandle, const reco::BeamSpot &vertexBeamSpot, MonitorElement *phi1, MonitorElement *eta1, MonitorElement *pT1, MonitorElement *phi2, MonitorElement *eta2, MonitorElement *pT2)
Definition: BPHMonitor.cc:1133
double pz() const final
z coordinate of momentum vector
bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb) override
T sqrt(T t)
Definition: SSEVec.h:18
double p4[4]
Definition: TauolaWrapper.h:92
MEbinning dR_binning_
Definition: BPHMonitor.h:115
bool accept(const edm::Event &event, const edm::EventSetup &setup)
To be called from analyze/filter() methods.
GlobalPoint position() const
T z() const
Definition: PV3DBase.h:64
double mu_mass
Definition: BPHMonitor.h:188
GenericTriggerEventFlag * den_genTriggerEventFlag_
Definition: BPHMonitor.h:163
BPHMonitor(const edm::ParameterSet &)
Definition: BPHMonitor.cc:12
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:118
double chi2() const
chi-squares
Definition: Vertex.h:98
MonitorElement * denominator
Definition: BPHMonitor.h:61
CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &tracks) const override
METME DiMuDS_
Definition: BPHMonitor.h:154
void initRun(const edm::Run &run, const edm::EventSetup &setup)
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:28
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
MEbinning phi_binning_
Definition: BPHMonitor.h:110
static void fillHistoPSetDescription(edm::ParameterSetDescription &pset)
Definition: BPHMonitor.cc:984
T min(T a, T b)
Definition: MathUtil.h:58
MEbinning ds_binning_
Definition: BPHMonitor.h:118
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:74
double p2[4]
Definition: TauolaWrapper.h:90
void analyze(edm::Event const &iEvent, edm::EventSetup const &iSetup) override
Definition: BPHMonitor.cc:347
double maxmass_
Definition: BPHMonitor.h:177
MEbinning prob_binning_
Definition: BPHMonitor.h:120
MEbinning cos_binning_
Definition: BPHMonitor.h:119
double dxdz() const
dxdz slope
Definition: BeamSpot.h:82
double ndof() const
Definition: Vertex.h:105
std::vector< std::string > hltpaths_den
Definition: BPHMonitor.h:197
StringCutObjectSelector< reco::Muon, true > muoSelection_
Definition: BPHMonitor.h:165
METME muz0_
Definition: BPHMonitor.h:126
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:639
METME DiMuProb_
Definition: BPHMonitor.h:153
MEbinning d0_binning_
Definition: BPHMonitor.h:113
int k[5][pyjets_maxn]
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: BPHMonitor.cc:219
std::vector< TriggerObject > TriggerObjectCollection
collection of trigger physics objects (e.g., all isolated muons)
Definition: TriggerObject.h:81
double minmassJpsiTk
Definition: BPHMonitor.h:186
Float e1
Definition: deltaR.h:20
const edm::InputTag filterTag(trigger::size_type index) const
Definition: TriggerEvent.h:103
edm::EDGetTokenT< reco::PhotonCollection > phToken_
Definition: BPHMonitor.h:108
edm::EDGetTokenT< reco::BeamSpot > bsToken_
Definition: BPHMonitor.h:106
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: BPHMonitor.cc:997
~BPHMonitor() override
Definition: BPHMonitor.cc:131
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:279
double minDS
Definition: BPHMonitor.h:193
StringCutObjectSelector< reco::Muon, true > muoSelection_ref
Definition: BPHMonitor.h:166
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:136
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
This class takes a vector of HLT paths and returns a weight based on their HLT and L1 prescales...
std::vector< size_type > Keys
T rerr(const GlobalPoint &aPoint) const
double minmassUpsilon
Definition: BPHMonitor.h:182
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
const T & get() const
Definition: EventSetup.h:59
static MEbinning getHistoPSet(edm::ParameterSet pset)
Definition: BPHMonitor.cc:138
double py() const final
y coordinate of momentum vector
void add(std::string const &label, ParameterSetDescription const &psetDescription)
double maxmassUpsilon
Definition: BPHMonitor.h:181
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:510
Float e2
Definition: deltaR.h:21
double minmassTkTk
Definition: BPHMonitor.h:184
fixed size matrix
METME muEta_
Definition: BPHMonitor.h:123
HLT enums.
double p1[4]
Definition: TauolaWrapper.h:89
METME DiMuMass_
Definition: BPHMonitor.h:156
METME DiMuPt_
Definition: BPHMonitor.h:151
MEbinning mass_binning_
Definition: BPHMonitor.h:116
METME mud0_
Definition: BPHMonitor.h:125
double y0() const
y coordinate
Definition: BeamSpot.h:66
METME mu1Pt_
Definition: BPHMonitor.h:130
METME mu3Phi_
Definition: BPHMonitor.h:138
std::vector< std::string > hltpaths_num
Definition: BPHMonitor.h:196
METME mu3Eta_
Definition: BPHMonitor.h:139
double minmassJpsi
Definition: BPHMonitor.h:180
void initRun(const edm::Run &run, const edm::EventSetup &setup)
To be called from beginRun() methods.
METME DiMuPVcos_
Definition: BPHMonitor.h:152
long double T
METME DiMudR_
Definition: BPHMonitor.h:157
T x() const
Definition: PV3DBase.h:62
double eta() const
momentum pseudorapidity
Definition: Particle.h:110
GenericTriggerEventFlag * num_genTriggerEventFlag_
Definition: BPHMonitor.h:162
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
bool isValid() const
METME muPhi_
Definition: BPHMonitor.h:122
StringCutObjectSelector< reco::Track, true > trSelection_ref
Definition: BPHMonitor.h:199
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:633
Definition: Run.h:43
double p3[4]
Definition: TauolaWrapper.h:91
double x0() const
x coordinate
Definition: BeamSpot.h:64
double min_dR
Definition: BPHMonitor.h:189