Inheritance diagram for ZMCHistogrammer:

Public Member Functions
	ZMCHistogrammer (const edm::ParameterSet &pset)
Private Member Functions
virtual void	analyze (const edm::Event &event, const edm::EventSetup &setup)
Private Attributes
double	angMax_
edm::InputTag	gen_
TH1F *	h_etaResZ_
TH1F *	h_etaZ_
TH1F *	h_etaZMC_
TH1F *	h_invmMuMu_
TH1F *	h_invmMuMuMC_
TH1F *	h_mRatioZMuMu_
TH1F *	h_mRatioZMuMuMC_
TH1F *	h_mResZ_
TH1F *	h_mResZMuMu_
TH1F *	h_mResZMuMuMC_
TH1F *	h_mZ_
TH1F *	h_mZMC_
TH1F *	h_nZ_
TH1F *	h_nZMC_
TH1F *	h_phiResZ_
TH1F *	h_phiZ_
TH1F *	h_phiZMC_
TH1F *	h_ptResZ_
TH1F *	h_ptZ_
TH1F *	h_ptZMC_
TH1F *	h_rapidityResZ_
TH1F *	h_rapidityZ_
TH1F *	h_rapidityZMC_
TH1F *	h_thetaResZ_
TH1F *	h_thetaZ_
TH1F *	h_thetaZMC_
TH1F *	h_weight_histo
bool	isMCatNLO_
double	massMax_
double	massResMax_
edm::InputTag	match_
unsigned int	nbinsAng_
unsigned int	nbinsMass_
unsigned int	nbinsMassRes_
unsigned int	nbinsPt_
double	ptMax_
edm::InputTag	z_

Detailed Description

Definition at line 5 of file ZMCHistogrammer.cc.

Constructor & Destructor Documentation

ZMCHistogrammer::ZMCHistogrammer ( const edm::ParameterSet & pset )

Definition at line 45 of file ZMCHistogrammer.cc.

References angMax_, gather_cfg::cout, h_etaResZ_, h_etaZ_, h_etaZMC_, h_invmMuMu_, h_invmMuMuMC_, h_mRatioZMuMu_, h_mRatioZMuMuMC_, h_mResZ_, h_mResZMuMu_, h_mResZMuMuMC_, h_mZ_, h_mZMC_, h_nZ_, h_nZMC_, h_phiResZ_, h_phiZ_, h_phiZMC_, h_ptResZ_, h_ptZ_, h_ptZMC_, h_rapidityResZ_, h_rapidityZ_, h_rapidityZMC_, h_thetaResZ_, h_thetaZ_, h_thetaZMC_, h_weight_histo, TFileDirectory::make(), massMax_, massResMax_, nbinsAng_, nbinsMass_, nbinsMassRes_, nbinsPt_, and ptMax_.

                                                         :
  z_(pset.getParameter<InputTag>("z")),
  gen_(pset.getParameter<InputTag>("gen")), 
  match_(pset.getParameter<InputTag>("match")), 
  nbinsMass_(pset.getUntrackedParameter<unsigned int>("nbinsMass")),
  nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
  nbinsAng_(pset.getUntrackedParameter<unsigned int>("nbinsAng")),
  nbinsMassRes_(pset.getUntrackedParameter<unsigned int>("nbinsMassRes")),
  massMax_(pset.getUntrackedParameter<double>("massMax")),
  ptMax_(pset.getUntrackedParameter<double>("ptMax")),
  angMax_(pset.getUntrackedParameter<double>("angMax")), 
  massResMax_(pset.getUntrackedParameter<double>("massResMax")),
  isMCatNLO_(pset.getUntrackedParameter<bool>("isMCatNLO"))  { 
  cout << ">>> Z Histogrammer constructor" << endl;
  Service<TFileService> fs;
  TFileDirectory ZHisto = fs->mkdir( "ZRecoHisto" );
  TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
  TFileDirectory ZResHisto = fs->mkdir( "ZResolutionHisto" );
  //TFileDirectory Z2vs3MCHisto = fs->mkdir( "Z2vs3MCHisto" );
  h_nZ_ = ZHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
  h_mZ_ = ZHisto.make<TH1F>("ZMass", "Z mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
  h_ptZ_ = ZHisto.make<TH1F>("ZPt", "Z p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
  h_phiZ_ = ZHisto.make<TH1F>("ZPhi", "Z #phi", nbinsAng_,  -angMax_, angMax_);
  h_thetaZ_ = ZHisto.make<TH1F>("Ztheta", "Z #theta", nbinsAng_,  0, angMax_);
  h_etaZ_ = ZHisto.make<TH1F>("ZEta", "Z #eta", nbinsAng_,  -angMax_, angMax_);
  h_rapidityZ_ = ZHisto.make<TH1F>("ZRapidity", "Z rapidity", nbinsAng_,  -angMax_, angMax_);
  h_invmMuMu_ = ZHisto.make<TH1F>("MuMuMass", "#mu #mu invariant mass",nbinsMass_,  0, massMax_);
  h_weight_histo  = ZHisto.make<TH1F>("weight_histo","weight_histo",20,-10,10); 
  h_nZMC_ = ZMCHisto.make<TH1F>("ZMCNumber", "number of Z MC particles", 11, -0.5, 10.5);
  h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
  h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
  h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_,  -angMax_, angMax_);
  h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_,  0, angMax_);
  h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_,  -angMax_, angMax_);
  h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC rapidity",  nbinsAng_,  -angMax_, angMax_);
  h_invmMuMuMC_ = ZMCHisto.make<TH1F>("MuMuMCMass", "#mu #mu MC invariant mass", 
                                      nbinsMass_,  0, massMax_);
  /*
  h_mZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCMass", "Z MC st 2 vs st 3 mass (GeV/c^{2})", 
                                        nbinsMassRes_, -massResMax_, massResMax_);
  h_ptZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPt", "Z MC st 2 vs st 3 p_{t} (GeV/c)", 
                                         nbinsPt_, -ptMax_, ptMax_);
  h_phiZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPhi", "Z MC st 2 vs st 3 #phi", 
                                          nbinsAng_,  -angMax_, angMax_);
  h_thetaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCTheta", "Z MC st 2 vs st 3 #theta", 
                                            nbinsAng_,  -angMax_, angMax_);
  h_etaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCEta", "Z MC st 2 vs st 3 #eta", 
                                          nbinsAng_,  -angMax_, angMax_);
  h_rapidityZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCRapidity", "Z MC st 2 vs st 3 rapidity", 
                                       nbinsAng_,  -angMax_, angMax_);
  */
  h_mResZ_ = ZResHisto.make<TH1F>("ZMassResolution", "Z mass Resolution (GeV/c^{2})", 
                                  nbinsMassRes_, -massResMax_, massResMax_);
  h_ptResZ_ = ZResHisto.make<TH1F>("ZPtResolution", "Z p_{t} Resolution (GeV/c)", 
                                   nbinsPt_, -ptMax_, ptMax_);
  h_phiResZ_ = ZResHisto.make<TH1F>("ZPhiResolution", "Z #phi Resolution", 
                                    nbinsAng_,  -angMax_, angMax_);
  h_thetaResZ_ = ZResHisto.make<TH1F>("ZThetaResolution", "Z #theta Resolution", 
                                      nbinsAng_, -angMax_, angMax_);
  h_etaResZ_ = ZResHisto.make<TH1F>("ZEtaResolution", "Z #eta Resolution", 
                                    nbinsAng_,  -angMax_, angMax_);
  h_rapidityResZ_ = ZResHisto.make<TH1F>("ZRapidityResolution", "Z rapidity Resolution", 
                                         nbinsAng_,  -angMax_, angMax_);
  h_mResZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassResolution", 
                                      "Z Reco vs matched final state #mu #mu mass Difference (GeV/c^{2})", 
                                      nbinsMassRes_, -massResMax_, massResMax_);
  h_mRatioZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassRatio", 
                                        "Z Reco vs matched final state #mu #mu mass Ratio", 
                                        4000, 0, 2);
  h_mResZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassResolution", 
                                        "Z vs final state #mu #mu MC mass Difference (GeV/c^{2})", 
                                        nbinsMassRes_/2 + 1, -2*massResMax_/nbinsMassRes_, massResMax_);
  h_mRatioZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassRatio", 
                                          "Z vs final state #mu #mu MC mass Ratio", 
                                          2002, 0.999, 2);
}

Member Function Documentation

void ZMCHistogrammer::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	setup
	)		`[private, virtual]`

Implements edm::EDAnalyzer.

Definition at line 122 of file ZMCHistogrammer.cc.

                                                                               { 
  cout << ">>> Z Histogrammer analyze" << endl;
  Handle<CandidateView> z;
  Handle<CandidateView> gen;
  Handle<std::vector<GenParticleRef> > match;
  event.getByLabel(z_, z);
  event.getByLabel(gen_, gen);
  event.getByLabel(match_, match);
  h_nZ_->Fill(z->size());
 
// get HepMC::GenEvent ...
   Handle<HepMCProduct> evt_h;
   event.getByType(evt_h);
   HepMC::GenEvent* evt = new  HepMC::GenEvent(*(evt_h->GetEvent()));


   // get weight and fill it to histogram
  HepMC::WeightContainer weights = evt->weights();   
  double weight = weights.front();
  if(!weight) weight=1.;
  h_weight_histo->Fill(weight); 
  
  if(isMCatNLO_) {
    weight > 0 ?  weight=1. : weight=-1.;
  }

  for(unsigned int i = 0; i < z->size(); ++i) {
    const Candidate &zCand = (*z)[i];
    h_mZ_->Fill(zCand.mass(),weight );
    h_ptZ_->Fill(zCand.pt(),weight);
    h_phiZ_->Fill(zCand.phi(),weight);
    h_thetaZ_->Fill(zCand.theta(),weight);
    h_etaZ_->Fill(zCand.eta(),weight);
    h_rapidityZ_->Fill(zCand.rapidity(),weight);
    CandidateBaseRef zCandRef = z->refAt(i);   
   
 GenParticleRef zMCMatch = (*match)[i];
    if(zMCMatch.isNonnull() && zMCMatch->pdgId()==23) {
      h_mResZ_->Fill(zCandRef->mass() - zMCMatch->mass());
      h_ptResZ_->Fill(zCandRef->pt() - zMCMatch->pt());
      h_phiResZ_->Fill(zCandRef->phi() - zMCMatch->phi());
      h_thetaResZ_->Fill(zCandRef->theta() - zMCMatch->theta());
      h_etaResZ_->Fill(zCandRef->eta() - zMCMatch->eta());
      h_rapidityResZ_->Fill(zCandRef->rapidity() - zMCMatch->rapidity());
      const Candidate * dau0 = zMCMatch->daughter(0);
      const Candidate * dau1 = zMCMatch->daughter(1);
      for(unsigned int i0 = 0; i0 < dau0->numberOfDaughters(); ++i0) {
        const Candidate * ddau0 = dau0->daughter(i0);
        if(abs(ddau0->pdgId())==13 && ddau0->status()==1) {
          dau0 = ddau0; break;
        }
      }
      for(unsigned int i1 = 0; i1 < dau1->numberOfDaughters(); ++i1) {
        const Candidate * ddau1 = dau1->daughter(i1);
        if(abs(ddau1->pdgId())==13 && ddau1->status()==1) {
          dau1 = ddau1; break;
        }
      }
      assert(abs(dau0->pdgId())==13 && dau0->status()==1);
      assert(abs(dau1->pdgId())==13 && dau1->status()==1);
      double invMass = (dau0->p4()+dau1->p4()).mass();
      h_invmMuMu_->Fill(invMass,weight);
      h_mResZMuMu_->Fill(zCand.mass() - invMass);
      h_mRatioZMuMu_->Fill(zCand.mass()/invMass);
    }
  }
  h_nZMC_->Fill(gen->size());
  for(unsigned int i = 0; i < gen->size(); ++i) {
    const Candidate &genCand = (*gen)[i];
    if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0
      cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters() 
           << " daughters" << endl;
    if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0
      cout << ">>> Z0 found, with " << genCand.numberOfDaughters() 
           << " daughters" << endl;
      h_mZMC_->Fill(genCand.mass(),weight);
      h_ptZMC_->Fill(genCand.pt(),weight);
      h_phiZMC_->Fill(genCand.phi(),weight);
      h_thetaZMC_->Fill(genCand.theta(),weight);
      h_etaZMC_->Fill(genCand.eta(),weight);
      h_rapidityZMC_->Fill(genCand.rapidity(),weight);
      Particle::LorentzVector pZ(0, 0, 0, 0);
      int nMu = 0;
      for(unsigned int j = 0; j < genCand.numberOfDaughters(); ++j) { 
        const Candidate *dauGen = genCand.daughter(j);
        /*
        if((dauGen->pdgId() == 23) && (dauGen->status() == 2)) { 
          h_mZ2vs3MC_->Fill(genCand.mass() - dauGen->mass());
          h_ptZ2vs3MC_->Fill(genCand.pt() - dauGen->pt());
          h_phiZ2vs3MC_->Fill(genCand.phi() - dauGen->phi());
          h_thetaZ2vs3MC_->Fill(genCand.theta() - dauGen->theta());
          h_etaZ2vs3MC_->Fill(genCand.eta() - dauGen->eta());
          h_rapidityZ2vs3MC_->Fill(genCand.rapidity() - dauGen->rapidity());
        }
        */
        if((abs(dauGen->pdgId()) == 13) && (dauGen->numberOfDaughters() != 0)) {
          //we are looking for photons of final state radiation
          cout << ">>> The muon " << j 
               << " has " << dauGen->numberOfDaughters() << " daughters" <<endl;
          for(unsigned int k = 0; k < dauGen->numberOfDaughters(); ++k) {
            const Candidate * dauMuGen = dauGen->daughter(k);
            cout << ">>> Mu " << j 
                 << " daughter MC " << k 
                 << " PDG Id " << dauMuGen->pdgId() 
                 << ", status " << dauMuGen->status() 
                 << ", charge " << dauMuGen->charge() 
                 << endl;
            if(abs(dauMuGen->pdgId()) == 13 && dauMuGen->status() ==1) {
              pZ += dauMuGen->p4();
              nMu ++;
            }
          }
        }
      }
      assert(nMu == 2);
      double mZ = pZ.mass();
      h_invmMuMuMC_->Fill(mZ,weight);
      h_mResZMuMuMC_->Fill(genCand.mass() - mZ);
      h_mRatioZMuMuMC_->Fill(genCand.mass()/mZ);
    }
  }
}