CMSSW Reference Manual

00001 #include "FWCore/Framework/interface/EDAnalyzer.h"
00002 #include "FWCore/ParameterSet/interface/InputTag.h"
00003 #include "TH1.h"
00004 
00005 class ZHistogrammer : public edm::EDAnalyzer {
00006 public:
00007   ZHistogrammer(const edm::ParameterSet& pset); 
00008 private:
00009   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00010   edm::InputTag  z_, gen_, match_;
00011   size_t nbinsMass_, nbinsPt_, nbinsAng_, nbinsMassRes_;
00012   double massMax_, ptMax_, angMax_, massResMax_;
00013   TH1F *h_nZ_, *h_mZ_, *h_ptZ_, *h_phiZ_, *h_thetaZ_, *h_etaZ_, *h_rapidityZ_;
00014   TH1F *h_invmMuMu_;
00015   TH1F *h_nZMC_, *h_mZMC_, *h_ptZMC_, *h_phiZMC_, *h_thetaZMC_, *h_etaZMC_, *h_rapidityZMC_;
00016   TH1F *h_invmMuMuMC_;
00017   //TH1F *h_mZ2vs3MC_, *h_ptZ2vs3MC_, *h_phiZ2vs3MC_, *h_thetaZ2vs3MC_, *h_etaZ2vs3MC_, *h_rapidityZ2vs3MC_;
00018   TH1F *h_mResZ_, *h_ptResZ_, *h_phiResZ_, *h_thetaResZ_, *h_etaResZ_, *h_rapidityResZ_;
00019   TH1F *h_mResZMuMu_, *h_mRatioZMuMu_;
00020   TH1F *h_mResZMuMuMC_, *h_mRatioZMuMuMC_;
00021 };
00022 
00023 #include "DataFormats/Candidate/interface/Candidate.h"
00024 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00025 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00026 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00027 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "DataFormats/Common/interface/Handle.h"
00030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00031 #include "FWCore/ServiceRegistry/interface/Service.h"
00032 #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
00033 #include <cmath>
00034 #include <iostream>
00035 
00036 using namespace std;
00037 using namespace reco;
00038 using namespace edm;
00039 
00040 ZHistogrammer::ZHistogrammer(const ParameterSet& pset) :
00041   z_(pset.getParameter<InputTag>("z")),
00042   gen_(pset.getParameter<InputTag>("gen")), 
00043   match_(pset.getParameter<InputTag>("match")), 
00044   nbinsMass_(pset.getUntrackedParameter<size_t>("nbinsMass")),
00045   nbinsPt_(pset.getUntrackedParameter<size_t>("nbinsPt")),
00046   nbinsAng_(pset.getUntrackedParameter<size_t>("nbinsAng")),
00047   nbinsMassRes_(pset.getUntrackedParameter<size_t>("nbinsMassRes")),
00048   massMax_(pset.getUntrackedParameter<double>("massMax")),
00049   ptMax_(pset.getUntrackedParameter<double>("ptMax")),
00050   angMax_(pset.getUntrackedParameter<double>("angMax")), 
00051   massResMax_(pset.getUntrackedParameter<double>("massResMax")) { 
00052   cout << ">>> Z Histogrammer constructor" << endl;
00053   Service<TFileService> fs;
00054   TFileDirectory ZHisto = fs->mkdir( "ZRecoHisto" );
00055   TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
00056   TFileDirectory ZResHisto = fs->mkdir( "ZResolutionHisto" );
00057   //TFileDirectory Z2vs3MCHisto = fs->mkdir( "Z2vs3MCHisto" );
00058   h_nZ_ = ZHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
00059   h_mZ_ = ZHisto.make<TH1F>("ZMass", "Z mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
00060   h_ptZ_ = ZHisto.make<TH1F>("ZPt", "Z p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00061   h_phiZ_ = ZHisto.make<TH1F>("ZPhi", "Z #phi", nbinsAng_,  -angMax_, angMax_);
00062   h_thetaZ_ = ZHisto.make<TH1F>("Ztheta", "Z #theta", nbinsAng_,  0, angMax_);
00063   h_etaZ_ = ZHisto.make<TH1F>("ZEta", "Z #eta", nbinsAng_,  -angMax_, angMax_);
00064   h_rapidityZ_ = ZHisto.make<TH1F>("ZRapidity", "Z rapidity", nbinsAng_,  -angMax_, angMax_);
00065   h_invmMuMu_ = ZHisto.make<TH1F>("MuMuMass", "#mu #mu invariant mass", 
00066                                   nbinsMass_,  0, massMax_);
00067   h_nZMC_ = ZMCHisto.make<TH1F>("ZMCNumber", "number of Z MC particles", 11, -0.5, 10.5);
00068   h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
00069   h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00070   h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_,  -angMax_, angMax_);
00071   h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_,  0, angMax_);
00072   h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_,  -angMax_, angMax_);
00073   h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC rapidity", 
00074                                        nbinsAng_,  -angMax_, angMax_);
00075   h_invmMuMuMC_ = ZMCHisto.make<TH1F>("MuMuMCMass", "#mu #mu MC invariant mass", 
00076                                       nbinsMass_,  0, massMax_);
00077   /*
00078   h_mZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCMass", "Z MC st 2 vs st 3 mass (GeV/c^{2})", 
00079                                         nbinsMassRes_, -massResMax_, massResMax_);
00080   h_ptZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPt", "Z MC st 2 vs st 3 p_{t} (GeV/c)", 
00081                                          nbinsPt_, -ptMax_, ptMax_);
00082   h_phiZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPhi", "Z MC st 2 vs st 3 #phi", 
00083                                           nbinsAng_,  -angMax_, angMax_);
00084   h_thetaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCTheta", "Z MC st 2 vs st 3 #theta", 
00085                                             nbinsAng_,  -angMax_, angMax_);
00086   h_etaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCEta", "Z MC st 2 vs st 3 #eta", 
00087                                           nbinsAng_,  -angMax_, angMax_);
00088   h_rapidityZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCRapidity", "Z MC st 2 vs st 3 rapidity", 
00089                                        nbinsAng_,  -angMax_, angMax_);
00090   */
00091   h_mResZ_ = ZResHisto.make<TH1F>("ZMassResolution", "Z mass Resolution (GeV/c^{2})", 
00092                                   nbinsMassRes_, -massResMax_, massResMax_);
00093   h_ptResZ_ = ZResHisto.make<TH1F>("ZPtResolution", "Z p_{t} Resolution (GeV/c)", 
00094                                    nbinsPt_, -ptMax_, ptMax_);
00095   h_phiResZ_ = ZResHisto.make<TH1F>("ZPhiResolution", "Z #phi Resolution", 
00096                                     nbinsAng_,  -angMax_, angMax_);
00097   h_thetaResZ_ = ZResHisto.make<TH1F>("ZThetaResolution", "Z #theta Resolution", 
00098                                       nbinsAng_, -angMax_, angMax_);
00099   h_etaResZ_ = ZResHisto.make<TH1F>("ZEtaResolution", "Z #eta Resolution", 
00100                                     nbinsAng_,  -angMax_, angMax_);
00101   h_rapidityResZ_ = ZResHisto.make<TH1F>("ZRapidityResolution", "Z rapidity Resolution", 
00102                                          nbinsAng_,  -angMax_, angMax_);
00103   h_mResZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassResolution", 
00104                                       "Z Reco vs matched final state #mu #mu mass Difference (GeV/c^{2})", 
00105                                       nbinsMassRes_, -massResMax_, massResMax_);
00106   h_mRatioZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassRatio", 
00107                                         "Z Reco vs matched final state #mu #mu mass Ratio", 
00108                                         4000, 0, 2);
00109   h_mResZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassResolution", 
00110                                         "Z vs final state #mu #mu MC mass Difference (GeV/c^{2})", 
00111                                         nbinsMassRes_/2 + 1, -2*massResMax_/nbinsMassRes_, massResMax_);
00112   h_mRatioZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassRatio", 
00113                                           "Z vs final state #mu #mu MC mass Ratio", 
00114                                           2002, 0.999, 2);
00115 }
00116 
00117 void ZHistogrammer::analyze(const edm::Event& event, const edm::EventSetup& setup) { 
00118   cout << ">>> Z Histogrammer analyze" << endl;
00119   Handle<CandidateCollection> z;
00120   Handle<CandidateCollection> gen;
00121   Handle<CandMatchMap> match;
00122   event.getByLabel(z_, z);
00123   event.getByLabel(gen_, gen);
00124   event.getByLabel(match_, match);
00125   h_nZ_->Fill(z->size());
00126   for(size_t i = 0; i < z->size(); ++i) {
00127     const Candidate &zCand = (*z)[i];
00128     h_mZ_->Fill(zCand.mass());
00129     h_ptZ_->Fill(zCand.pt());
00130     h_phiZ_->Fill(zCand.phi());
00131     h_thetaZ_->Fill(zCand.theta());
00132     h_etaZ_->Fill(zCand.eta());
00133     h_rapidityZ_->Fill(zCand.rapidity());
00134     CandidateRef zCandRef(z, i);
00135     CandidateRef zMCMatch = (*match)[zCandRef];
00136     if(zMCMatch.isNonnull() && zMCMatch->pdgId()==23) {
00137       h_mResZ_->Fill(zCandRef->mass() - zMCMatch->mass());
00138       h_ptResZ_->Fill(zCandRef->pt() - zMCMatch->pt());
00139       h_phiResZ_->Fill(zCandRef->phi() - zMCMatch->phi());
00140       h_thetaResZ_->Fill(zCandRef->theta() - zMCMatch->theta());
00141       h_etaResZ_->Fill(zCandRef->eta() - zMCMatch->eta());
00142       h_rapidityResZ_->Fill(zCandRef->rapidity() - zMCMatch->rapidity());
00143       const Candidate * dau0 = zMCMatch->daughter(0);
00144       const Candidate * dau1 = zMCMatch->daughter(1);
00145       for(size_t i0 = 0; i0 < dau0->numberOfDaughters(); ++i0) {
00146         const Candidate * ddau0 = dau0->daughter(i0);
00147         if(abs(ddau0->pdgId())==13 && ddau0->status()==1) {
00148           dau0 = ddau0; break;
00149         }
00150       }
00151       for(size_t i1 = 0; i1 < dau1->numberOfDaughters(); ++i1) {
00152         const Candidate * ddau1 = dau1->daughter(i1);
00153         if(abs(ddau1->pdgId())==13 && ddau1->status()==1) {
00154           dau1 = ddau1; break;
00155         }
00156       }
00157       assert(abs(dau0->pdgId())==13 && dau0->status()==1);
00158       assert(abs(dau1->pdgId())==13 && dau1->status()==1);
00159       double invMass = (dau0->p4()+dau1->p4()).mass();
00160       h_invmMuMu_->Fill(invMass);
00161       h_mResZMuMu_->Fill(zCand.mass() - invMass);
00162       h_mRatioZMuMu_->Fill(zCand.mass()/invMass);
00163     }
00164   }
00165   h_nZMC_->Fill(gen->size());
00166   for(size_t i = 0; i < gen->size(); ++i) {
00167     const Candidate &genCand = (*gen)[i];
00168     if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0
00169       cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters() 
00170            << " daughters" << endl;
00171     if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0
00172       cout << ">>> Z0 found, with " << genCand.numberOfDaughters() 
00173            << " daughters" << endl;
00174       h_mZMC_->Fill(genCand.mass());
00175       h_ptZMC_->Fill(genCand.pt());
00176       h_phiZMC_->Fill(genCand.phi());
00177       h_thetaZMC_->Fill(genCand.theta());
00178       h_etaZMC_->Fill(genCand.eta());
00179       h_rapidityZMC_->Fill(genCand.rapidity());
00180       Particle::LorentzVector pZ(0, 0, 0, 0);
00181       int nMu = 0;
00182       for(size_t j = 0; j < genCand.numberOfDaughters(); ++j) { 
00183         const Candidate *dauGen = genCand.daughter(j);
00184         /*
00185         if((dauGen->pdgId() == 23) && (dauGen->status() == 2)) { 
00186           h_mZ2vs3MC_->Fill(genCand.mass() - dauGen->mass());
00187           h_ptZ2vs3MC_->Fill(genCand.pt() - dauGen->pt());
00188           h_phiZ2vs3MC_->Fill(genCand.phi() - dauGen->phi());
00189           h_thetaZ2vs3MC_->Fill(genCand.theta() - dauGen->theta());
00190           h_etaZ2vs3MC_->Fill(genCand.eta() - dauGen->eta());
00191           h_rapidityZ2vs3MC_->Fill(genCand.rapidity() - dauGen->rapidity());
00192         }
00193         */
00194         if((abs(dauGen->pdgId()) == 13) && (dauGen->numberOfDaughters() != 0)) {
00195           //we are looking for photons of final state radiation
00196           cout << ">>> The muon " << j 
00197                << " has " << dauGen->numberOfDaughters() << " daughters" <<endl;
00198           for(size_t k = 0; k < dauGen->numberOfDaughters(); ++k) {
00199             const Candidate * dauMuGen = dauGen->daughter(k);
00200             cout << ">>> Mu " << j 
00201                  << " daughter MC " << k 
00202                  << " PDG Id " << dauMuGen->pdgId() 
00203                  << ", status " << dauMuGen->status() 
00204                  << ", charge " << dauMuGen->charge() 
00205                  << endl;
00206             if(abs(dauMuGen->pdgId()) == 13 && dauMuGen->status() ==1) {
00207               pZ += dauMuGen->p4();
00208               nMu ++;
00209             }
00210           }
00211         }
00212       }
00213       assert(nMu == 2);
00214       double mZ = pZ.mass();
00215       h_invmMuMuMC_->Fill(mZ);
00216       h_mResZMuMuMC_->Fill(genCand.mass() - mZ);
00217       h_mRatioZMuMuMC_->Fill(genCand.mass()/mZ);
00218     }
00219   }
00220 }
00221 
00222 #include "FWCore/Framework/interface/MakerMacros.h"
00223 
00224 DEFINE_FWK_MODULE(ZHistogrammer);
00225
ZHistogrammer.cc
