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