Reference Manual

Go to the documentation of this file.
00001 /* This Class Header */
00002 #include "DQMOffline/Muon/src/DiMuonHistograms.h"
00003 
00004 /* Collaborating Class Header */
00005 #include "FWCore/Framework/interface/MakerMacros.h"
00006 #include "FWCore/Framework/interface/Frameworkfwd.h"
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 #include "FWCore/Framework/interface/ESHandle.h"
00010 #include "FWCore/Utilities/interface/Exception.h"
00011 #include "DataFormats/MuonReco/interface/Muon.h"
00012 #include "DataFormats/MuonReco/interface/MuonFwd.h" 
00013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00014 
00015 using namespace edm;
00016 
00017 #include "TLorentzVector.h"
00018 #include "TFile.h"
00019 #include <vector>
00020 #include "math.h"
00021 
00022 
00023 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00024 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00025 
00026 
00027 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00028 #include "DataFormats/TrackReco/interface/Track.h"
00029 
00030 /* C++ Headers */
00031 #include <iostream>
00032 #include <fstream>
00033 #include <cmath>
00034 using namespace std;
00035 using namespace edm;
00036 
00037 DiMuonHistograms::DiMuonHistograms(const edm::ParameterSet& pSet, MuonServiceProxy *theService):MuonAnalyzerBase(theService){ 
00038   parameters = pSet;
00039 }
00040 
00041 DiMuonHistograms::~DiMuonHistograms() { }
00042 
00043 void DiMuonHistograms::beginJob(DQMStore * dbe) {
00044   
00045   metname = "DiMuonhistograms";
00046   LogTrace(metname)<<"[DiMuonHistograms] Parameters initialization";
00047   dbe->setCurrentFolder("Muons/DiMuonHistograms");  
00048 
00049   theMuonCollectionLabel = parameters.getParameter<edm::InputTag>("MuonCollection");
00050   
00051   etaBin = parameters.getParameter<int>("etaBin");
00052   etaBBin = parameters.getParameter<int>("etaBBin");
00053   etaEBin = parameters.getParameter<int>("etaEBin");
00054   
00055   etaBMin = parameters.getParameter<double>("etaBMin");
00056   etaBMax = parameters.getParameter<double>("etaBMax");
00057   etaECMin = parameters.getParameter<double>("etaECMin");
00058   etaECMax = parameters.getParameter<double>("etaECMax");
00059   
00060   LowMassMin = parameters.getParameter<double>("LowMassMin");  
00061   LowMassMax = parameters.getParameter<double>("LowMassMax");
00062   HighMassMin = parameters.getParameter<double>("HighMassMin");
00063   HighMassMax = parameters.getParameter<double>("HighMassMax");
00064 
00065   int nBin = 0;
00066   for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00067     if (iEtaRegion==0) { EtaName = "";         nBin = etaBin;} 
00068     if (iEtaRegion==1) { EtaName = "_Barrel";  nBin = etaBBin;}
00069     if (iEtaRegion==2) { EtaName = "_EndCap";  nBin = etaEBin;}
00070     
00071     GlbGlbMuon_LM.push_back(dbe->book1D("GlbGlbMuon_LM"+EtaName,"InvMass_{GLB,GLB}"+EtaName,nBin, LowMassMin, LowMassMax));
00072     TrkTrkMuon_LM.push_back(dbe->book1D("TrkTrkMuon_LM"+EtaName,"InvMass_{TRK,TRK}"+EtaName,nBin, LowMassMin, LowMassMax));
00073     StaTrkMuon_LM.push_back(dbe->book1D("StaTrkMuon_LM"+EtaName,"InvMass_{STA,TRK}"+EtaName,nBin, LowMassMin, LowMassMax));
00074     
00075     GlbGlbMuon_HM.push_back(dbe->book1D("GlbGlbMuon_HM"+EtaName,"InvMass_{GLB,GLB}"+EtaName,nBin, HighMassMin, HighMassMax));
00076     TrkTrkMuon_HM.push_back(dbe->book1D("TrkTrkMuon_HM"+EtaName,"InvMass_{TRK,TRK}"+EtaName,nBin, HighMassMin, HighMassMax));
00077     StaTrkMuon_HM.push_back(dbe->book1D("StaTrkMuon_HM"+EtaName,"InvMass_{STA,TRK}"+EtaName,nBin, HighMassMin, HighMassMax));
00078     
00079     // arround the Z peak
00080     TightTightMuon.push_back(dbe->book1D("TightTightMuon"+EtaName,"InvMass_{Tight,Tight}"+EtaName,nBin, 55.0, 125.0));
00081 
00082     // low-mass resonances
00083     SoftSoftMuon.push_back(dbe->book1D("SoftSoftMuon"+EtaName,"InvMass_{Soft,Soft}"+EtaName,nBin, 5.0, 55.0));
00084   }
00085 }
00086 void DiMuonHistograms::analyze(const edm::Event & iEvent,const edm::EventSetup& iSetup) {
00087 
00088   LogTrace(metname)<<"[DiMuonHistograms] Analyze the mu in different eta regions";
00089   edm::Handle<reco::MuonCollection> muons;
00090   iEvent.getByLabel(theMuonCollectionLabel, muons);
00091 
00092   reco::BeamSpot beamSpot;
00093   Handle<reco::BeamSpot> beamSpotHandle;
00094   iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
00095   beamSpot = *beamSpotHandle;
00096 
00097   if(!muons.isValid()) return;
00098 
00099   // Loop on muon collection
00100   TLorentzVector Mu1, Mu2;
00101   float charge = 99.;
00102   float InvMass = -99.;
00103 
00104   for (reco::MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1!=muons->end(); ++recoMu1) {
00105     LogTrace(metname)<<"[DiMuonHistograms] loop over 1st muon"<<endl;
00106 
00107     // Loop on second muons to fill invariant mass plots
00108     for (reco::MuonCollection::const_iterator recoMu2 = recoMu1; recoMu2!=muons->end(); ++recoMu2){ 
00109       LogTrace(metname)<<"[DiMuonHistograms] loop over 2nd muon"<<endl;
00110       if (recoMu1==recoMu2) continue;
00111       
00112       // Global-Global Muon
00113       if (recoMu1->isGlobalMuon() && recoMu2->isGlobalMuon()) {
00114         LogTrace(metname)<<"[DiMuonHistograms] Glb-Glb pair"<<endl;
00115         reco::TrackRef recoCombinedGlbTrack1 = recoMu1->combinedMuon();
00116         reco::TrackRef recoCombinedGlbTrack2 = recoMu2->combinedMuon();
00117         Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(), recoCombinedGlbTrack1->py(),recoCombinedGlbTrack1->pz(), recoCombinedGlbTrack1->p());
00118         Mu2.SetPxPyPzE(recoCombinedGlbTrack2->px(), recoCombinedGlbTrack2->py(),recoCombinedGlbTrack2->pz(), recoCombinedGlbTrack2->p());
00119         
00120         charge  = recoCombinedGlbTrack1->charge()*recoCombinedGlbTrack2->charge();
00121         if (charge < 0) {
00122           InvMass = (Mu1+Mu2).M();
00123           for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00124             if (iEtaRegion==0) {EtaCutMin= 0.;         EtaCutMax=2.4;       }
00125             if (iEtaRegion==1) {EtaCutMin= etaBMin;    EtaCutMax=etaBMax;   } 
00126             if (iEtaRegion==2) {EtaCutMin= etaECMin;   EtaCutMax=etaECMax;  }
00127             
00128             if(fabs(recoCombinedGlbTrack1->eta())>EtaCutMin && fabs(recoCombinedGlbTrack1->eta())<EtaCutMax && 
00129                fabs(recoCombinedGlbTrack2->eta())>EtaCutMin && fabs(recoCombinedGlbTrack2->eta())<EtaCutMax){
00130               if (InvMass < LowMassMax)  GlbGlbMuon_LM[iEtaRegion]->Fill(InvMass);
00131               if (InvMass > HighMassMin) GlbGlbMuon_HM[iEtaRegion]->Fill(InvMass);
00132             }
00133           }
00134         }
00135         // Also Tight-Tight Muon Selection
00136         if (recoMu1->isGlobalMuon() && recoMu1->isTrackerMuon() && recoMu1->combinedMuon()->normalizedChi2()<10. 
00137             && recoMu1->combinedMuon()->hitPattern().numberOfValidMuonHits()>0 && fabs(recoMu1->combinedMuon()->dxy(beamSpot.position()))<0.2 
00138             && recoMu1->combinedMuon()->hitPattern().numberOfValidPixelHits()>0 && recoMu1->numberOfMatches() > 1   
00139             && recoMu2->isGlobalMuon() && recoMu2->isTrackerMuon() && recoMu2->combinedMuon()->normalizedChi2()<10. 
00140             && recoMu2->combinedMuon()->hitPattern().numberOfValidMuonHits()>0 && fabs(recoMu2->combinedMuon()->dxy(beamSpot.position()))<0.2 
00141             && recoMu2->combinedMuon()->hitPattern().numberOfValidPixelHits()>0 && recoMu2->numberOfMatches() > 1) {
00142           LogTrace(metname)<<"[DiMuonHistograms] Tight-Tight pair"<<endl;
00143           for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00144             if (iEtaRegion==0) {EtaCutMin= 0.;         EtaCutMax=2.4;       }
00145             if (iEtaRegion==1) {EtaCutMin= etaBMin;    EtaCutMax=etaBMax;   } 
00146             if (iEtaRegion==2) {EtaCutMin= etaECMin;   EtaCutMax=etaECMax;  }
00147 
00148             if(fabs(recoCombinedGlbTrack1->eta())>EtaCutMin && fabs(recoCombinedGlbTrack1->eta())<EtaCutMax && 
00149                fabs(recoCombinedGlbTrack2->eta())>EtaCutMin && fabs(recoCombinedGlbTrack2->eta())<EtaCutMax){
00150               if (InvMass > 55. && InvMass < 125.) TightTightMuon[iEtaRegion]->Fill(InvMass);
00151             }
00152           }
00153         }
00154       }
00155     
00156       // Now check for STA-TRK 
00157       if (recoMu2->isStandAloneMuon() && recoMu1->isTrackerMuon()) {
00158         LogTrace(metname)<<"[DiMuonHistograms] STA-Trk pair"<<endl;
00159         reco::TrackRef recoStaTrack = recoMu2->standAloneMuon();
00160         reco::TrackRef recoTrack    = recoMu1->track();
00161         Mu2.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(),recoStaTrack->pz(), recoStaTrack->p());
00162         Mu1.SetPxPyPzE(recoTrack->px(), recoTrack->py(),recoTrack->pz(), recoTrack->p());
00163 
00164         charge  = recoStaTrack->charge()*recoTrack->charge();
00165         if (charge < 0) {
00166           InvMass = (Mu1+Mu2).M();
00167           for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00168             if (iEtaRegion==0) {EtaCutMin= 0.;         EtaCutMax=2.4;       }
00169             if (iEtaRegion==1) {EtaCutMin= etaBMin;    EtaCutMax=etaBMax;   } 
00170             if (iEtaRegion==2) {EtaCutMin= etaECMin;   EtaCutMax=etaECMax;  }
00171             
00172             if(fabs(recoStaTrack->eta())>EtaCutMin && fabs(recoStaTrack->eta())<EtaCutMax &&
00173                fabs(recoTrack->eta())>EtaCutMin && fabs(recoTrack->eta())<EtaCutMax){
00174               if (InvMass < LowMassMax)  StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
00175               if (InvMass > HighMassMin) StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
00176             }
00177           }
00178         }
00179       }
00180       if (recoMu1->isStandAloneMuon() && recoMu2->isTrackerMuon()) {
00181         LogTrace(metname)<<"[DiMuonHistograms] STA-Trk pair"<<endl;
00182         reco::TrackRef recoStaTrack = recoMu1->standAloneMuon();
00183         reco::TrackRef recoTrack    = recoMu2->track();
00184         Mu1.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(),recoStaTrack->pz(), recoStaTrack->p());
00185         Mu2.SetPxPyPzE(recoTrack->px(), recoTrack->py(),recoTrack->pz(), recoTrack->p());
00186 
00187         charge  = recoStaTrack->charge()*recoTrack->charge();
00188         if (charge < 0) {
00189           InvMass = (Mu1+Mu2).M();
00190           for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00191             if (iEtaRegion==0) {EtaCutMin= 0.;         EtaCutMax=2.4;       }
00192             if (iEtaRegion==1) {EtaCutMin= etaBMin;    EtaCutMax=etaBMax;   } 
00193             if (iEtaRegion==2) {EtaCutMin= etaECMin;   EtaCutMax=etaECMax;  }
00194             
00195             if(fabs(recoStaTrack->eta())>EtaCutMin && fabs(recoStaTrack->eta())<EtaCutMax &&
00196                fabs(recoTrack->eta())>EtaCutMin && fabs(recoTrack->eta())<EtaCutMax){
00197               if (InvMass < LowMassMax)  StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
00198               if (InvMass > HighMassMin) StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
00199             }
00200           }
00201         }
00202       }
00203 
00204       // TRK-TRK dimuon 
00205       if (recoMu1->isTrackerMuon() && recoMu2->isTrackerMuon()) {
00206         LogTrace(metname)<<"[DiMuonHistograms] Trk-Trk dimuon pair"<<endl;
00207         reco::TrackRef recoTrack2 = recoMu2->track();
00208         reco::TrackRef recoTrack1 = recoMu1->track();
00209         Mu2.SetPxPyPzE(recoTrack2->px(), recoTrack2->py(),recoTrack2->pz(), recoTrack2->p());
00210         Mu1.SetPxPyPzE(recoTrack1->px(), recoTrack1->py(),recoTrack1->pz(), recoTrack1->p());
00211         
00212         charge  = recoTrack1->charge()*recoTrack2->charge();
00213         if (charge < 0) {
00214           InvMass = (Mu1+Mu2).M();
00215           for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00216             if (iEtaRegion==0) {EtaCutMin= 0.;         EtaCutMax=2.4;       }
00217             if (iEtaRegion==1) {EtaCutMin= etaBMin;    EtaCutMax=etaBMax;   } 
00218             if (iEtaRegion==2) {EtaCutMin= etaECMin;   EtaCutMax=etaECMax;  }
00219             
00220             if(fabs(recoTrack1->eta())>EtaCutMin && fabs(recoTrack1->eta())<EtaCutMax &&
00221                fabs(recoTrack2->eta())>EtaCutMin && fabs(recoTrack2->eta())<EtaCutMax){
00222               if (InvMass < LowMassMax)  TrkTrkMuon_LM[iEtaRegion]->Fill(InvMass);
00223               if (InvMass > HighMassMin) TrkTrkMuon_HM[iEtaRegion]->Fill(InvMass);
00224             }
00225           }
00226         }
00227         
00228         if (recoMu1->isTrackerMuon() && recoMu2->isTrackerMuon() &&
00229             recoMu1->innerTrack()->found() > 11 && recoMu2->innerTrack()->found() &&
00230             recoMu1->innerTrack()->chi2()/recoMu1->innerTrack()->ndof() < 4.0 && 
00231             recoMu2->innerTrack()->chi2()/recoMu2->innerTrack()->ndof() < 4.0 &&
00232             recoMu1->numberOfMatches() > 0 && recoMu2->numberOfMatches() > 0 && 
00233             recoMu1->innerTrack()->hitPattern().pixelLayersWithMeasurement() > 1 &&
00234             recoMu2->innerTrack()->hitPattern().pixelLayersWithMeasurement() > 1 &&
00235             fabs(recoMu1->innerTrack()->dxy()) < 3.0 && fabs(recoMu1->innerTrack()->dxy()) < 3.0 &&
00236             fabs(recoMu1->innerTrack()->dz()) < 15.0 && fabs(recoMu1->innerTrack()->dz()) < 15.0){
00237           
00238           if (charge < 0) {
00239             InvMass = (Mu1+Mu2).M();
00240             for (unsigned int iEtaRegion=0; iEtaRegion<3; iEtaRegion++){
00241               if (iEtaRegion==0) {EtaCutMin= 0.;         EtaCutMax=2.4; }
00242               if (iEtaRegion==1) {EtaCutMin= etaBMin;    EtaCutMax=etaBMax; }
00243               if (iEtaRegion==2) {EtaCutMin= etaECMin;   EtaCutMax=etaECMax; }
00244               
00245               if(fabs(recoTrack1->eta())>EtaCutMin && fabs(recoTrack1->eta())<EtaCutMax &&
00246                  fabs(recoTrack2->eta())>EtaCutMin && fabs(recoTrack2->eta())<EtaCutMax){
00247                 SoftSoftMuon[iEtaRegion]->Fill(InvMass); 
00248               }
00249             }
00250           }
00251         } 
00252       }
00253     } //muon2
00254   } //Muon1
00255 }
00256
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/DQMOffline/Muon/src/DiMuonHistograms.cc
