test
CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DiMuonHistograms.cc
Go to the documentation of this file.
1 /* This Class Header */
2 #include "DQMOffline/Muon/interface/DiMuonHistograms.h"
3 
4 /* Collaborating Class Header */
5 #include "FWCore/Framework/interface/MakerMacros.h"
6 #include "FWCore/Framework/interface/Frameworkfwd.h"
7 #include "FWCore/Framework/interface/Event.h"
8 #include "FWCore/ParameterSet/interface/ParameterSet.h"
9 #include "FWCore/Framework/interface/ESHandle.h"
10 #include "FWCore/Utilities/interface/Exception.h"
11 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
12 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
13 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
14 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
15 #include "DataFormats/TrackReco/interface/Track.h"
16 
17 #include "DQMServices/Core/interface/DQMStore.h"
18 #include "DQMServices/Core/interface/MonitorElement.h"
19 
20 #include "TLorentzVector.h"
21 #include "TFile.h"
22 #include <vector>
23 #include <cmath>
24 
25 /* C++ Headers */
26 #include <iostream>
27 #include <fstream>
28 #include <cmath>
29 using namespace std;
30 using namespace edm;
31 
32 DiMuonHistograms::DiMuonHistograms(const edm::ParameterSet& pSet){
33 
34  // initialise parameters:
35  parameters = pSet;
36 
37 
38  // declare consumes:
39  theMuonCollectionLabel_ = consumes<reco::MuonCollection> (parameters.getParameter<edm::InputTag>("MuonCollection"));
40  theBeamSpotLabel_ = mayConsume<reco::BeamSpot> (parameters.getParameter<edm::InputTag>("BeamSpotLabel"));
41  theVertexLabel_ = consumes<reco::VertexCollection>(parameters.getParameter<edm::InputTag>("VertexLabel"));
42 
43  etaBin = parameters.getParameter<int>("etaBin");
44  etaBBin = parameters.getParameter<int>("etaBBin");
45  etaEBin = parameters.getParameter<int>("etaEBin");
46 
47  etaBinLM = parameters.getParameter<int>("etaBinLM");
48  etaBBinLM = parameters.getParameter<int>("etaBBinLM");
49  etaEBinLM = parameters.getParameter<int>("etaEBinLM");
50 
51  etaBMin = parameters.getParameter<double>("etaBMin");
52  etaBMax = parameters.getParameter<double>("etaBMax");
53  etaECMin = parameters.getParameter<double>("etaECMin");
54  etaECMax = parameters.getParameter<double>("etaECMax");
55 
56  LowMassMin = parameters.getParameter<double>("LowMassMin");
57  LowMassMax = parameters.getParameter<double>("LowMassMax");
58  HighMassMin = parameters.getParameter<double>("HighMassMin");
59  HighMassMax = parameters.getParameter<double>("HighMassMax");
60 
61 }
62 
63 DiMuonHistograms::~DiMuonHistograms() { }
64 
65 void DiMuonHistograms::bookHistograms(DQMStore::IBooker & ibooker,
66  edm::Run const & /*iRun*/,
67  edm::EventSetup const & /* iSetup */){
68 
69  ibooker.cd();
70  ibooker.setCurrentFolder("Muons/DiMuonHistograms");
71 
72  int nBin = 0, nBinLM = 0;
73  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
74  if (iEtaRegion==0) { EtaName = ""; nBin = etaBin;}
75  if (iEtaRegion==1) { EtaName = "_Barrel"; nBin = etaBBin;}
76  if (iEtaRegion==2) { EtaName = "_EndCap"; nBin = etaEBin;}
77 
78 
79  if (etaBinLM == 0) { nBinLM = nBin;} //for HeavyIons
80  else{
81  if (iEtaRegion==0) { EtaName = ""; nBinLM = etaBinLM;}
82  if (iEtaRegion==1) { EtaName = "_Barrel"; nBinLM = etaBBinLM;}
83  if (iEtaRegion==2) { EtaName = "_EndCap"; nBinLM = etaEBinLM;}
84  }
85 
86  GlbGlbMuon_LM.push_back(ibooker.book1D("GlbGlbMuon_LM"+EtaName,"InvMass_{GLB,GLB}"+EtaName,nBinLM, LowMassMin, LowMassMax));
87  TrkTrkMuon_LM.push_back(ibooker.book1D("TrkTrkMuon_LM"+EtaName,"InvMass_{TRK,TRK}"+EtaName,nBinLM, LowMassMin, LowMassMax));
88  StaTrkMuon_LM.push_back(ibooker.book1D("StaTrkMuon_LM"+EtaName,"InvMass_{STA,TRK}"+EtaName,nBinLM, LowMassMin, LowMassMax));
89 
90  GlbGlbMuon_HM.push_back(ibooker.book1D("GlbGlbMuon_HM"+EtaName,"InvMass_{GLB,GLB}"+EtaName,nBin, HighMassMin, HighMassMax));
91  TrkTrkMuon_HM.push_back(ibooker.book1D("TrkTrkMuon_HM"+EtaName,"InvMass_{TRK,TRK}"+EtaName,nBin, HighMassMin, HighMassMax));
92  StaTrkMuon_HM.push_back(ibooker.book1D("StaTrkMuon_HM"+EtaName,"InvMass_{STA,TRK}"+EtaName,nBin, HighMassMin, HighMassMax));
93 
94  // arround the Z peak
95  TightTightMuon.push_back(ibooker.book1D("TightTightMuon"+EtaName,"InvMass_{Tight,Tight}"+EtaName,nBin, 55.0, 125.0));
96 
97  // low-mass resonances
98  SoftSoftMuon.push_back(ibooker.book1D("SoftSoftMuon"+EtaName,"InvMass_{Soft,Soft}"+EtaName,nBin, 5.0, 55.0));
99  }
100 }
101 
102 void DiMuonHistograms::analyze(const edm::Event & iEvent,const edm::EventSetup& iSetup){
103 
104  LogTrace(metname)<<"[DiMuonHistograms] Analyze the mu in different eta regions";
105  edm::Handle<reco::MuonCollection> muons;
106  iEvent.getByToken(theMuonCollectionLabel_, muons);
107 
108  // =================================================================================
109  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
110  reco::Vertex::Point posVtx;
111  reco::Vertex::Error errVtx;
112  unsigned int theIndexOfThePrimaryVertex = 999.;
113 
114  edm::Handle<reco::VertexCollection> vertex;
115  iEvent.getByToken(theVertexLabel_, vertex);
116  if (vertex.isValid()){
117  for (unsigned int ind=0; ind<vertex->size(); ++ind) {
118  if ( (*vertex)[ind].isValid() && !((*vertex)[ind].isFake()) ) {
119  theIndexOfThePrimaryVertex = ind;
120  break;
121  }
122  }
123  }
124 
125  if (theIndexOfThePrimaryVertex<100) {
126  posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
127  errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
128  }
129  else {
130  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
131 
132  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
133  iEvent.getByToken(theBeamSpotLabel_,recoBeamSpotHandle);
134  reco::BeamSpot bs = *recoBeamSpotHandle;
135 
136  posVtx = bs.position();
137  errVtx(0,0) = bs.BeamWidthX();
138  errVtx(1,1) = bs.BeamWidthY();
139  errVtx(2,2) = bs.sigmaZ();
140  }
141 
142  const reco::Vertex vtx(posVtx,errVtx);
143 
144  if(!muons.isValid()) return;
145 
146  // Loop on muon collection
147  TLorentzVector Mu1, Mu2;
148  float charge = 99.;
149  float InvMass = -99.;
150 
151  for (reco::MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1!=muons->end(); ++recoMu1) {
152  LogTrace(metname)<<"[DiMuonHistograms] loop over 1st muon"<<endl;
153 
154  // Loop on second muons to fill invariant mass plots
155  for (reco::MuonCollection::const_iterator recoMu2 = recoMu1; recoMu2!=muons->end(); ++recoMu2){
156  LogTrace(metname)<<"[DiMuonHistograms] loop over 2nd muon"<<endl;
157  if (recoMu1==recoMu2) continue;
158 
159  // Global-Global Muon
160  if (recoMu1->isGlobalMuon() && recoMu2->isGlobalMuon()) {
161  LogTrace(metname)<<"[DiMuonHistograms] Glb-Glb pair"<<endl;
162  reco::TrackRef recoCombinedGlbTrack1 = recoMu1->combinedMuon();
163  reco::TrackRef recoCombinedGlbTrack2 = recoMu2->combinedMuon();
164  Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(), recoCombinedGlbTrack1->py(),recoCombinedGlbTrack1->pz(), recoCombinedGlbTrack1->p());
165  Mu2.SetPxPyPzE(recoCombinedGlbTrack2->px(), recoCombinedGlbTrack2->py(),recoCombinedGlbTrack2->pz(), recoCombinedGlbTrack2->p());
166 
167  charge = recoCombinedGlbTrack1->charge()*recoCombinedGlbTrack2->charge();
168  if (charge < 0) {
169  InvMass = (Mu1+Mu2).M();
170  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
171  if (iEtaRegion==0) {EtaCutMin= 0.; EtaCutMax=2.4; }
172  if (iEtaRegion==1) {EtaCutMin= etaBMin; EtaCutMax=etaBMax; }
173  if (iEtaRegion==2) {EtaCutMin= etaECMin; EtaCutMax=etaECMax; }
174 
175  if(fabs(recoCombinedGlbTrack1->eta())>EtaCutMin && fabs(recoCombinedGlbTrack1->eta())<EtaCutMax &&
176  fabs(recoCombinedGlbTrack2->eta())>EtaCutMin && fabs(recoCombinedGlbTrack2->eta())<EtaCutMax){
177  if (InvMass < LowMassMax) GlbGlbMuon_LM[iEtaRegion]->Fill(InvMass);
178  if (InvMass > HighMassMin) GlbGlbMuon_HM[iEtaRegion]->Fill(InvMass);
179  }
180  }
181  }
182  // Also Tight-Tight Muon Selection
183 
184  if ( muon::isTightMuon(*recoMu1, vtx) &&
185  muon::isTightMuon(*recoMu2, vtx) ) {
186 
187  LogTrace(metname)<<"[DiMuonHistograms] Tight-Tight pair"<<endl;
188  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
189  if (iEtaRegion==0) {EtaCutMin= 0.; EtaCutMax=2.4; }
190  if (iEtaRegion==1) {EtaCutMin= etaBMin; EtaCutMax=etaBMax; }
191  if (iEtaRegion==2) {EtaCutMin= etaECMin; EtaCutMax=etaECMax; }
192 
193  if(fabs(recoCombinedGlbTrack1->eta())>EtaCutMin && fabs(recoCombinedGlbTrack1->eta())<EtaCutMax &&
194  fabs(recoCombinedGlbTrack2->eta())>EtaCutMin && fabs(recoCombinedGlbTrack2->eta())<EtaCutMax){
195  if (InvMass > 55. && InvMass < 125.) TightTightMuon[iEtaRegion]->Fill(InvMass);
196  }
197  }
198  }
199  }
200 
201  // Now check for STA-TRK
202  if (recoMu2->isStandAloneMuon() && recoMu1->isTrackerMuon()) {
203  LogTrace(metname)<<"[DiMuonHistograms] STA-Trk pair"<<endl;
204  reco::TrackRef recoStaTrack = recoMu2->standAloneMuon();
205  reco::TrackRef recoTrack = recoMu1->track();
206  Mu2.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(),recoStaTrack->pz(), recoStaTrack->p());
207  Mu1.SetPxPyPzE(recoTrack->px(), recoTrack->py(),recoTrack->pz(), recoTrack->p());
208 
209  charge = recoStaTrack->charge()*recoTrack->charge();
210  if (charge < 0) {
211  InvMass = (Mu1+Mu2).M();
212  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
213  if (iEtaRegion==0) {EtaCutMin= 0.; EtaCutMax=2.4; }
214  if (iEtaRegion==1) {EtaCutMin= etaBMin; EtaCutMax=etaBMax; }
215  if (iEtaRegion==2) {EtaCutMin= etaECMin; EtaCutMax=etaECMax; }
216 
217  if(fabs(recoStaTrack->eta())>EtaCutMin && fabs(recoStaTrack->eta())<EtaCutMax &&
218  fabs(recoTrack->eta())>EtaCutMin && fabs(recoTrack->eta())<EtaCutMax){
219  if (InvMass < LowMassMax) StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
220  if (InvMass > HighMassMin) StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
221  }
222  }
223  }
224  }
225  if (recoMu1->isStandAloneMuon() && recoMu2->isTrackerMuon()) {
226  LogTrace(metname)<<"[DiMuonHistograms] STA-Trk pair"<<endl;
227  reco::TrackRef recoStaTrack = recoMu1->standAloneMuon();
228  reco::TrackRef recoTrack = recoMu2->track();
229  Mu1.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(),recoStaTrack->pz(), recoStaTrack->p());
230  Mu2.SetPxPyPzE(recoTrack->px(), recoTrack->py(),recoTrack->pz(), recoTrack->p());
231 
232  charge = recoStaTrack->charge()*recoTrack->charge();
233  if (charge < 0) {
234  InvMass = (Mu1+Mu2).M();
235  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
236  if (iEtaRegion==0) {EtaCutMin= 0.; EtaCutMax=2.4; }
237  if (iEtaRegion==1) {EtaCutMin= etaBMin; EtaCutMax=etaBMax; }
238  if (iEtaRegion==2) {EtaCutMin= etaECMin; EtaCutMax=etaECMax; }
239 
240  if(fabs(recoStaTrack->eta())>EtaCutMin && fabs(recoStaTrack->eta())<EtaCutMax &&
241  fabs(recoTrack->eta())>EtaCutMin && fabs(recoTrack->eta())<EtaCutMax){
242  if (InvMass < LowMassMax) StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
243  if (InvMass > HighMassMin) StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
244  }
245  }
246  }
247  }
248 
249  // TRK-TRK dimuon
250  if (recoMu1->isTrackerMuon() && recoMu2->isTrackerMuon()) {
251  LogTrace(metname)<<"[DiMuonHistograms] Trk-Trk dimuon pair"<<endl;
252  reco::TrackRef recoTrack2 = recoMu2->track();
253  reco::TrackRef recoTrack1 = recoMu1->track();
254  Mu2.SetPxPyPzE(recoTrack2->px(), recoTrack2->py(),recoTrack2->pz(), recoTrack2->p());
255  Mu1.SetPxPyPzE(recoTrack1->px(), recoTrack1->py(),recoTrack1->pz(), recoTrack1->p());
256 
257  charge = recoTrack1->charge()*recoTrack2->charge();
258  if (charge < 0) {
259  InvMass = (Mu1+Mu2).M();
260  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
261  if (iEtaRegion==0) {EtaCutMin= 0.; EtaCutMax=2.4; }
262  if (iEtaRegion==1) {EtaCutMin= etaBMin; EtaCutMax=etaBMax; }
263  if (iEtaRegion==2) {EtaCutMin= etaECMin; EtaCutMax=etaECMax; }
264 
265  if(fabs(recoTrack1->eta())>EtaCutMin && fabs(recoTrack1->eta())<EtaCutMax &&
266  fabs(recoTrack2->eta())>EtaCutMin && fabs(recoTrack2->eta())<EtaCutMax){
267  if (InvMass < LowMassMax) TrkTrkMuon_LM[iEtaRegion]->Fill(InvMass);
268  if (InvMass > HighMassMin) TrkTrkMuon_HM[iEtaRegion]->Fill(InvMass);
269  }
270  }
271  }
272 
273 
274  LogTrace(metname)<<"[DiMuonHistograms] Soft-Soft pair"<<endl;
275 
276  if (muon::isSoftMuon(*recoMu1, vtx) &&
277  muon::isSoftMuon(*recoMu2, vtx) ) {
278 
279  if (charge < 0) {
280  InvMass = (Mu1+Mu2).M();
281  for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
282  if (iEtaRegion==0) {EtaCutMin= 0.; EtaCutMax=2.4; }
283  if (iEtaRegion==1) {EtaCutMin= etaBMin; EtaCutMax=etaBMax; }
284  if (iEtaRegion==2) {EtaCutMin= etaECMin; EtaCutMax=etaECMax; }
285 
286  if(fabs(recoTrack1->eta())>EtaCutMin && fabs(recoTrack1->eta())<EtaCutMax &&
287  fabs(recoTrack2->eta())>EtaCutMin && fabs(recoTrack2->eta())<EtaCutMax){
288  SoftSoftMuon[iEtaRegion]->Fill(InvMass);
289  }
290  }
291  }
292  }
293  }
294  } //muon2
295  } //Muon1
296 }
297 
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:44
DQMStore::IBooker::cd
void cd(void)
Definition: DQMStore.cc:268
reco::Vertex::Error
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:43
muon::isSoftMuon
bool isSoftMuon(const reco::Muon &, const reco::Vertex &)
Definition: MuonSelectors.cc:779
iEvent
int iEvent
Definition: GenABIO.cc:230
DQMStore::IBooker::book1D
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
reco::Vertex::Point
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
reco::BeamSpot::BeamWidthX
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:86
DiMuonHistograms::DiMuonHistograms
DiMuonHistograms(const edm::ParameterSet &pset)
Definition: DiMuonHistograms.cc:32
HLT_25ns10e33_v2_cff.parameters
tuple parameters
Definition: HLT_25ns10e33_v2_cff.py:3864
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:75
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:276
reco::BeamSpot::sigmaZ
double sigmaZ() const
sigma z
Definition: BeamSpot.h:80
DiMuonHistograms::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: DiMuonHistograms.cc:65
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:88
DiMuonHistograms::~DiMuonHistograms
virtual ~DiMuonHistograms()
Definition: DiMuonHistograms.cc:63
patZpeak.muons
tuple muons
Definition: patZpeak.py:38
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:62
muon::isTightMuon
bool isTightMuon(const reco::Muon &, const reco::Vertex &)
Definition: MuonSelectors.cc:745
DiMuonHistograms::analyze
void analyze(const edm::Event &, const edm::EventSetup &)
Definition: DiMuonHistograms.cc:102
edm::Run
Definition: Run.h:43
Generated for CMSSW Reference Manual by doxygen 1.8.5