PatMCMatchingExtended.cc

Go to the documentation of this file.

#include <map>
#include <string>

#include "TH1.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "Math/VectorUtil.h"

#include "DataFormats/PatCandidates/interface/Muon.h"


class PatMCMatchingExtended : public edm::EDAnalyzer {

public:
  explicit PatMCMatchingExtended(const edm::ParameterSet&);
  ~PatMCMatchingExtended();

private:

  virtual void beginJob() override ;
  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
  virtual void endJob() override ;

  // simple map to contain all histograms;
  // histograms are booked in the beginJob()
  // method
  std::map<std::string,TH1F*> histContainer_;

  // input tags
  edm::EDGetTokenT<edm::View<pat::Muon> > muonSrcToken_;

  //counts how often a genParticle with different charge gives a match
  unsigned int diffCharge;

  //how many muons have no match
  unsigned int noMatch;

  //how many muons have no status 1 or 3 match, but decay in flight
  unsigned int decayInFlight;

  //count the number of muons in all events
  unsigned int numberMuons;

};



PatMCMatchingExtended::PatMCMatchingExtended(const edm::ParameterSet& iConfig):
  histContainer_(),
  muonSrcToken_(consumes<edm::View<pat::Muon> >(iConfig.getUntrackedParameter<edm::InputTag>("muonSrc")))
{
}

PatMCMatchingExtended::~PatMCMatchingExtended()
{
}

void
PatMCMatchingExtended::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{

  // get muon collection
  edm::Handle<edm::View<pat::Muon> > muons;
  iEvent.getByToken(muonSrcToken_,muons);

  for(edm::View<pat::Muon>::const_iterator muon=muons->begin(); muon!=muons->end(); ++muon){
    if(muon->genParticleById(0,1).isNonnull() ){
      histContainer_["DR_status1Match"]->Fill( ROOT::Math::VectorUtil::DeltaR(muon->p4() , (muon->genParticleById(0,1) )->p4() ) );
      histContainer_["DPt_status1Match"]->Fill(muon->pt() - (muon->genParticleById(0,1) )->pt() );
    }
    if(muon->genParticleById(0,3).isNonnull() ){
      histContainer_["DR_status3Match"]->Fill( ROOT::Math::VectorUtil::DeltaR(muon->p4() , (muon->genParticleById(0,3) )->p4() ) );
      histContainer_["DPt_status3Match"]->Fill(muon->pt() - (muon->genParticleById(0,3) )->pt() );
    }
    if(muon->genParticleById(0,-1).isNonnull() ){
      histContainer_["DR_defaultMatch"]->Fill( ROOT::Math::VectorUtil::DeltaR(muon->p4() , (muon->genParticleById(0,-1) )->p4() ) );
      histContainer_["DPt_defaultMatch"]->Fill(muon->pt() - (muon->genParticleById(0,-1) )->pt() );
    }
    if(muon->genParticleById(0,1).isNull() && muon->genParticleById(0,3).isNull() && muon->genParticleById(0,-1).isNull()) noMatch++;
    if(muon->genParticleById(0,1).isNull() && muon->genParticleById(0,3).isNull() && muon->genParticleById(0,-1).isNonnull())decayInFlight++;



    if( muon->genParticleById(-13,0, 1).isNonnull() ){
      diffCharge++;
      std::cout<<" DIFF CHARGE!!! charge gen: "<< muon->genParticleById(-13,0, true)->charge()<< " charge reco: "<< muon->charge()<<std::endl;
    }
    numberMuons++;
  }

}

void
PatMCMatchingExtended::beginJob()
{
  // register to the TFileService
  edm::Service<TFileService> fs;

  // book histograms:
  //DR
  histContainer_["DR_defaultMatch"  ]=fs->make<TH1F>("DR_defaultMatch",   "DR_defaultMatch",     100, 0,  0.02);
  histContainer_["DR_status1Match"  ]=fs->make<TH1F>("DR_status1Match",   "DR_status1Match",     100, 0,  0.02);
  histContainer_["DR_status3Match"  ]=fs->make<TH1F>("DR_status3Match",   "DR_status3Match",     100, 0,  0.02);
  //DPT
  histContainer_["DPt_defaultMatch"  ]=fs->make<TH1F>("DPt_defaultMatch",   "DPt_defaultMatch",     10, 0,  1.2);
  histContainer_["DPt_status1Match"  ]=fs->make<TH1F>("DPt_status1Match",   "DPt_status1Match",     10, 0,  1.2);
  histContainer_["DPt_status3Match"  ]=fs->make<TH1F>("DPt_status3Match",   "DPt_status3Match",     10, 0,  1.2);
  //some counters
  diffCharge=0;
  noMatch=0;
  decayInFlight=0;
  numberMuons=0;

}

void
PatMCMatchingExtended::endJob()
{
    std::cout<<"diffcharge: "<< diffCharge <<std::endl;
    std::cout<<"noMatch: "<<  noMatch <<std::endl;
    std::cout<<"decayInFlight: "<< decayInFlight <<std::endl;
    std::cout<<"numberMuons: "<< numberMuons <<std::endl;
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(PatMCMatchingExtended);