Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2010/11/11 17:33:03 $
00005  *  $Revision: 1.6 $
00006  *  \author S. Bolognesi, Erik - CERN
00007  */
00008 
00009 #include "DQM/Physics/src/BPhysicsOniaDQM.h"
00010 
00011 #include "FWCore/Framework/interface/Event.h"
00012 #include "FWCore/Framework/interface/EventSetup.h"
00013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00014 
00015 #include "DQMServices/Core/interface/DQMStore.h"
00016 #include "DQMServices/Core/interface/MonitorElement.h"
00017 #include "FWCore/ServiceRegistry/interface/Service.h"
00018 
00019 #include "DataFormats/Common/interface/Handle.h"
00020 #include "DataFormats/MuonReco/interface/MuonFwd.h" 
00021 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00022 #include "DataFormats/TrackReco/interface/Track.h"
00023 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00024 #include "DataFormats/VertexReco/interface/Vertex.h"
00025 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00026 
00027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00028 
00029 using namespace std;
00030 using namespace edm;
00031 using namespace reco;
00032 
00033 BPhysicsOniaDQM::BPhysicsOniaDQM(const ParameterSet& parameters) {
00034   // Muon Collection Label
00035   theMuonCollectionLabel = parameters.getParameter<InputTag>("MuonCollection");
00036   vertex = parameters.getParameter<InputTag>("vertex");
00037 
00038   global_background = NULL;
00039   diMuonMass_global = NULL;
00040   tracker_background = NULL;
00041   diMuonMass_tracker = NULL;
00042   standalone_background = NULL;
00043   diMuonMass_standalone = NULL;
00044   
00045   glbSigCut = NULL;
00046   glbSigNoCut = NULL;
00047   glbBkgNoCut = NULL;
00048   staSigCut = NULL;
00049   staSigNoCut = NULL;
00050   staBkgNoCut = NULL;
00051   trkSigCut = NULL;
00052   trkSigNoCut = NULL;
00053   trkBkgNoCut = NULL;
00054 
00055   JPsiGlbYdLumi = NULL;
00056   JPsiStaYdLumi = NULL;
00057   JPsiTrkYdLumi = NULL;
00058 }
00059 
00060 BPhysicsOniaDQM::~BPhysicsOniaDQM() { 
00061 }
00062 
00063 void BPhysicsOniaDQM::beginJob() {
00064   // the services
00065   theDbe = Service<DQMStore>().operator->();
00066 
00067   metname = "oniaAnalyzer";
00068   LogTrace(metname)<<"[BPhysicsOniaDQM] Parameters initialization";
00069 
00070   if(theDbe!=NULL){
00071     theDbe->setCurrentFolder("Physics/BPhysics");  // Use folder with name of PAG
00072     global_background = theDbe->book1D("global_background", "Same-sign global-global dimuon mass", 750, 0, 15);
00073     diMuonMass_global = theDbe->book1D("diMuonMass_global", "Opposite-sign global-global dimuon mass", 750, 0, 15);
00074     tracker_background = theDbe->book1D("tracker_background", "Same-sign tracker-tracker (arbitrated) dimuon mass", 750, 0, 15);
00075     diMuonMass_tracker = theDbe->book1D("diMuonMass_tracker", "Opposite-sign tracker-tracker (arbitrated) dimuon mass", 750, 0, 15);
00076     standalone_background = theDbe->book1D("standalone_background", "Same-sign standalone-standalone dimuon mass", 500, 0, 15);
00077     diMuonMass_standalone = theDbe->book1D("diMuonMass_standalone", "Opposite-sign standalone-standalone dimuon mass", 500, 0, 15);
00078 
00079     glbSigCut = theDbe->book1D("glbSigCut", "Opposite-sign glb-glb dimuon mass", 650, 0, 130);
00080     glbSigNoCut = theDbe->book1D("glbSigNoCut", "Opposite-sign glb-glb dimuon mass (no cut)", 650, 0, 130);
00081     glbBkgNoCut = theDbe->book1D("glbBkgNoCut", "Same-sign glb-glb dimuon mass (no cut)", 650, 0, 130);
00082     staSigCut = theDbe->book1D("staSigCut", "Opposite-sign sta-sta dimuon mass", 430, 0, 129);
00083     staSigNoCut = theDbe->book1D("staSigNoCut", "Opposite-sign sta-sta dimuon mass (no cut)", 430, 0, 129);
00084     staBkgNoCut  = theDbe->book1D("staBkgNoCut", "Same-sign sta-sta dimuon mass (no cut)", 430, 0, 129);
00085     trkSigCut = theDbe->book1D("trkSigCut", "Opposite-sign trk-trk dimuon mass", 650, 0, 130);
00086     trkSigNoCut = theDbe->book1D("trkSigNoCut", "Opposite-sign trk-trk dimuon mass (no cut)", 650, 0, 130);
00087     trkBkgNoCut = theDbe->book1D("trkBkgNoCutt", "Same-sign trk-trk dimuon mass (no cut)", 650, 0, 130);
00088   }
00089 
00090 }
00091 
00092 void BPhysicsOniaDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
00093 
00094   LogTrace(metname)<<"[BPhysicsOniaDQM] Analysis of event # ";
00095   
00096   // Take the STA muon container
00097   Handle<MuonCollection> muons;
00098   iEvent.getByLabel(theMuonCollectionLabel,muons);
00099 
00100   Handle<reco::VertexCollection> privtxs;
00101   iEvent.getByLabel(vertex,privtxs);
00102   VertexCollection::const_iterator privtx;
00103 
00104   if(privtxs->begin() != privtxs->end()){
00105     privtx = privtxs->begin();
00106     RefVtx = privtx->position();
00107   } else {
00108     RefVtx.SetXYZ(0.,0.,0.);
00109   }
00110 
00111   if(muons.isValid()){
00112     for (MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1!=muons->end(); ++recoMu1){
00113 
00114       // only loop over the remaining muons if recoMu1 is one of the following
00115       if(recoMu1->isGlobalMuon() || recoMu1->isTrackerMuon() || recoMu1->isStandAloneMuon()){
00116         for (MuonCollection::const_iterator recoMu2 = recoMu1+1; recoMu2!=muons->end(); ++recoMu2){
00117 
00118           // fill the relevant histograms if recoMu2 satisfies one of the following
00119           if (recoMu1->isGlobalMuon() && recoMu2->isGlobalMuon()){
00120             math::XYZVector vec1 = recoMu1->globalTrack()->momentum();
00121             math::XYZVector vec2 = recoMu2->globalTrack()->momentum();
00122             float massJPsi = computeMass(vec1,vec2);
00123 
00124             // if opposite charges, fill glbSig, else fill glbBkg
00125             if (((*recoMu1).charge()*(*recoMu2).charge())<0) {
00126               if(diMuonMass_global!=NULL){  // BPhysicsOniaDQM original one
00127                 diMuonMass_global->Fill(massJPsi);
00128               }
00129 
00130               if(glbSigNoCut!=NULL){
00131                 glbSigNoCut->Fill(massJPsi);
00132                 if (selGlobalMuon(*recoMu1) && selGlobalMuon(*recoMu2)) {
00133                   if (glbSigCut!=NULL) glbSigCut->Fill(massJPsi);
00134                   if (massJPsi >= 3.0 && massJPsi <= 3.2) jpsiGlbSigPerLS++;
00135                 }
00136               }
00137             } else {
00138               if(global_background!=NULL){  // BPhysicsOniaDQM original one
00139                 global_background->Fill (massJPsi);
00140               }
00141 
00142               if(glbBkgNoCut!=NULL){
00143                 glbBkgNoCut->Fill(massJPsi);
00144               }
00145             }
00146           }
00147           
00148           if(recoMu1->isStandAloneMuon() && recoMu2->isStandAloneMuon() &&
00149             fabs(recoMu1->outerTrack()->d0()) < 5 && fabs(recoMu1->outerTrack()->dz()) < 30 &&
00150             fabs(recoMu2->outerTrack()->d0()) < 5 && fabs(recoMu2->outerTrack()->dz()) < 30){
00151             math::XYZVector vec1 = recoMu1->outerTrack()->momentum();
00152             math::XYZVector vec2 = recoMu2->outerTrack()->momentum();
00153             float massJPsi = computeMass(vec1,vec2);
00154 
00155             // if opposite charges, fill staSig, else fill staBkg
00156             if (((*recoMu1).charge()*(*recoMu2).charge())<0) {
00157               if(diMuonMass_standalone!=NULL){
00158                 diMuonMass_standalone->Fill(massJPsi);
00159               }
00160 
00161               if(staSigNoCut!=NULL){
00162                 staSigNoCut->Fill(massJPsi);
00163                 /*if (selStandaloneMuon(*recoMu1) && selStandaloneMuon(*recoMu2)) {
00164                   if (staSigCut!=NULL) staSigCut->Fill(massJPsi);
00165                   if (massJPsi >= 3.0 && massJPsi <= 3.2) jpsiStaSigPerLS++;
00166                 }*/
00167               }
00168             } else {
00169               if(standalone_background!=NULL){
00170                 standalone_background->Fill (massJPsi);
00171               }
00172 
00173               if(staBkgNoCut!=NULL){
00174                 staBkgNoCut->Fill(massJPsi);
00175               }
00176             }
00177           }
00178 
00179           if(recoMu1->isTrackerMuon() && recoMu2->isTrackerMuon() &&
00180             muon::isGoodMuon(*recoMu1, muon::TrackerMuonArbitrated) &&
00181             muon::isGoodMuon(*recoMu2, muon::TrackerMuonArbitrated)){
00182             math::XYZVector vec1 = recoMu1->innerTrack()->momentum();
00183             math::XYZVector vec2 = recoMu2->innerTrack()->momentum();
00184             float massJPsi = computeMass(vec1,vec2);
00185 
00186             // if opposite charges, fill trkSig, else fill trkBkg
00187             if (((*recoMu1).charge()*(*recoMu2).charge())<0) {
00188               if(diMuonMass_tracker!=NULL){
00189                 diMuonMass_tracker->Fill(massJPsi);
00190               }
00191 
00192               if(trkSigNoCut!=NULL){
00193                 trkSigNoCut->Fill(massJPsi);
00194                 if (selTrackerMuon(*recoMu1) && selTrackerMuon(*recoMu2)) {
00195                   if (trkSigCut!=NULL) trkSigCut->Fill(massJPsi);
00196                   if(massJPsi >= 3.0 && massJPsi <= 3.2) jpsiTrkSigPerLS++;
00197                 }
00198               }
00199             } else {
00200               if(tracker_background!=NULL){
00201                 tracker_background->Fill (massJPsi);
00202               }
00203 
00204               if(trkBkgNoCut!=NULL){
00205                 trkBkgNoCut->Fill(massJPsi);
00206               }
00207             }
00208           }
00209 
00210         }//end of 2nd MuonCollection
00211       }//end of GLB,STA,TRK muon check
00212     }//end of 1st MuonCollection
00213   }//Is this MuonCollection vaild?
00214 
00215 }
00216 
00217 void BPhysicsOniaDQM::endJob(void) {
00218   LogTrace(metname)<<"[BPhysicsOniaDQM] EndJob";
00219 }
00220 
00221 void BPhysicsOniaDQM::beginLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup)
00222 {
00223   LogTrace(metname)<<"[BPhysicsOniaDQM] Start of a LuminosityBlock";
00224   
00225   jpsiGlbSigPerLS = 0;
00226   jpsiStaSigPerLS = 0;
00227   jpsiTrkSigPerLS = 0;
00228 }
00229 
00230 void BPhysicsOniaDQM::endLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup)
00231 {
00232   LogTrace(metname)<<"[BPhysicsOniaDQM] Start of a LuminosityBlock";
00233 
00234   edm::Handle<LumiSummary> lumiSummary;
00235   lumiBlock.getByLabel("lumiProducer",lumiSummary);
00236 
00237   int LBlockNum = lumiBlock.id().luminosityBlock();
00238   
00239   jpsiGlbSig.insert( pair<int,int>(LBlockNum, jpsiGlbSigPerLS) );
00240   jpsiStaSig.insert( pair<int,int>(LBlockNum, jpsiStaSigPerLS) );
00241   jpsiTrkSig.insert( pair<int,int>(LBlockNum, jpsiTrkSigPerLS) );
00242 //  cout << "lumi: " << LBlockNum << "\t" << jpsiGlbSig[LBlockNum] << "\t" << jpsiStaSig[LBlockNum] << "\t" << jpsiTrkSig[LBlockNum] << endl;
00243 
00244   if (jpsiGlbSig.size()%5 != 0) return;
00245 
00246   theDbe->setCurrentFolder("Physics/BPhysics");
00247   if(JPsiGlbYdLumi!=NULL) {
00248     theDbe->removeElement("JPsiGlbYdLumi");   // Remove histograms from previous run
00249     theDbe->removeElement("JPsiStaYdLumi");
00250     theDbe->removeElement("JPsiTrkYdLumi");
00251   }
00252 
00253   int xmin = (*jpsiGlbSig.begin()).first;
00254   int xmax = (*jpsiGlbSig.rbegin()).first;
00255   int nx   = (xmax - xmin + 1)/5 + 1; // Merge 5 lumisections into 1 bin
00256 //  cout << "x-axis " << xmin << " " << xmax << endl;
00257 
00258   JPsiGlbYdLumi = theDbe->book1D("JPsiGlbYdLumi", "JPsi yield from global-global dimuon", nx, xmin, xmax);
00259   JPsiStaYdLumi = theDbe->book1D("JPsiStaYdLumi", "JPsi yield from standalone-standalone dimuon", nx, xmin, xmax);
00260   JPsiTrkYdLumi = theDbe->book1D("JPsiTrkYdLumi", "JPsi yield from tracker-tracker dimuon", nx, xmin, xmax);
00261 
00262   map<int,int>::iterator glb;
00263   map<int,int>::iterator sta;
00264   map<int,int>::iterator trk;
00265   for (glb = jpsiGlbSig.begin(); glb != jpsiGlbSig.end(); ++glb)
00266   {
00267     int bin = ((*glb).first - xmin + 1)/5 + 1;  //X-axis bin #
00268     sta = jpsiStaSig.find((*glb).first);
00269     trk = jpsiTrkSig.find((*glb).first);
00270     JPsiGlbYdLumi->setBinContent(bin,JPsiGlbYdLumi->getBinContent(bin)+(*glb).second);
00271     JPsiStaYdLumi->setBinContent(bin,JPsiStaYdLumi->getBinContent(bin)+(*sta).second);
00272     JPsiTrkYdLumi->setBinContent(bin,JPsiTrkYdLumi->getBinContent(bin)+(*trk).second);
00273 //    cout << "glb: " << bin << "\t" << (*glb).first << "\t" << (*glb).second << endl;
00274 //    cout << "sta: " << bin << "\t" << (*sta).first << "\t" << (*sta).second << endl;
00275 //    cout << "trk: " << bin << "\t" << (*trk).first << "\t" << (*trk).second << endl;
00276   }
00277 }
00278 
00279 void BPhysicsOniaDQM::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00280 {
00281   LogTrace(metname)<<"[BPhysicsOniaDQM] Start of a Run";
00282 }
00283 
00284 void BPhysicsOniaDQM::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00285 {
00286   LogTrace(metname)<<"[BPhysicsOniaDQM] End of a Run";
00287   
00288   if (!jpsiGlbSig.empty()) {
00289     jpsiGlbSig.clear();
00290     jpsiStaSig.clear();
00291     jpsiTrkSig.clear();
00292   }
00293 }
00294 
00295 float BPhysicsOniaDQM::computeMass(const math::XYZVector &vec1,const math::XYZVector &vec2){
00296   // mass of muon
00297   float massMu = 0.10566;
00298   float eMu1 = -999;
00299   if(massMu*massMu + vec1.Mag2()>0)
00300     eMu1 = sqrt(massMu*massMu + vec1.Mag2());
00301   float eMu2 = -999;
00302   if(massMu*massMu + vec2.Mag2()>0)
00303     eMu2 = sqrt(massMu*massMu + vec2.Mag2());
00304 
00305   float pJPsi = -999;
00306   if((vec1+vec2).Mag2()>0)
00307     pJPsi = sqrt((vec1+vec2).Mag2());
00308   float eJPsi = eMu1 + eMu2;
00309 
00310   float massJPsi = -999;
00311   if((eJPsi*eJPsi - pJPsi*pJPsi) > 0)
00312     massJPsi = sqrt(eJPsi*eJPsi - pJPsi*pJPsi);
00313  
00314  return massJPsi;
00315 }
00316 
00317 bool BPhysicsOniaDQM::isMuonInAccept(const reco::Muon &recoMu)
00318 {
00319   return (fabs(recoMu.eta()) < 2.4 &&
00320          ((fabs(recoMu.eta()) < 1.3 && recoMu.pt() > 3.3) ||
00321           (fabs(recoMu.eta()) > 1.3 && fabs(recoMu.eta()) < 2.2 && recoMu.p() > 2.9) ||
00322           (fabs(recoMu.eta()) > 2.2 && recoMu.pt() > 0.8)));
00323 }
00324 
00325 bool BPhysicsOniaDQM::selGlobalMuon(const reco::Muon &recoMu)
00326 {
00327   TrackRef iTrack = recoMu.innerTrack();
00328   const reco::HitPattern &p = iTrack->hitPattern();
00329   
00330   TrackRef gTrack = recoMu.globalTrack();
00331   const reco::HitPattern &q = gTrack->hitPattern();
00332 
00333   return (isMuonInAccept(recoMu) &&
00334           iTrack->found() > 11 &&
00335           gTrack->chi2()/gTrack->ndof() < 20.0 &&
00336           q.numberOfValidMuonHits() > 0 &&
00337           iTrack->chi2()/iTrack->ndof() < 4.0 &&
00338           //recoMu.muonID("TrackerMuonArbitrated") &&
00339           //recoMu.muonID("TMLastStationAngTight") &&
00340           p.pixelLayersWithMeasurement() > 1 &&
00341           fabs(iTrack->dxy(RefVtx)) < 3.0 &&
00342           fabs(iTrack->dz(RefVtx)) < 15.0 );
00343 }
00344 
00345 bool BPhysicsOniaDQM::selTrackerMuon(const reco::Muon &recoMu)
00346 {
00347   TrackRef iTrack = recoMu.innerTrack();
00348   const reco::HitPattern &p = iTrack->hitPattern();
00349 
00350   return (isMuonInAccept(recoMu) &&
00351           iTrack->found() > 11 &&
00352           iTrack->chi2()/iTrack->ndof() < 4.0 &&
00353           //recoMu.muonID("TrackerMuonArbitrated") &&
00354           //recoMu.muonID("TMLastStationAngTight") &&
00355           p.pixelLayersWithMeasurement() > 1 &&
00356           fabs(iTrack->dxy(RefVtx)) < 3.0 &&
00357           fabs(iTrack->dz(RefVtx)) < 15.0 );
00358 }
00359
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/DQM/Physics/src/BPhysicsOniaDQM.cc
