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