MuScleFitMuonSelector.h

#ifndef MUSCLEFITMUONSELECTOR
#define MUSCLEFITMUONSELECTOR

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
#include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/GenMuonPair.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/Muon.h"

#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"

#include "MuScleFitPlotter.h"

#include <vector>

/**
* This class is used to perform the selection of muon pairs in MuScleFit. <br>
* It receives the event and returns a muon collection containing all the selected muons.
*/
typedef reco::Particle::LorentzVector lorentzVector;

class MuScleFitMuonSelector
{
 public:
  MuScleFitMuonSelector(const edm::InputTag & muonLabel, const int muonType, const bool PATmuons,
			const std::vector<int> & resfind,
			const bool speedup, const std::string & genParticlesName,
			const bool compareToSimTracks, const edm::InputTag & simTracksCollectionName,
			const bool sherpa, const bool debug) :
    muonLabel_(muonLabel),
    muonType_(muonType),
    PATmuons_(PATmuons),
    resfind_(resfind),
    speedup_(speedup),
    genParticlesName_(genParticlesName),
    compareToSimTracks_(compareToSimTracks),
    simTracksCollectionName_(simTracksCollectionName),
    sherpa_(sherpa),
    debug_(debug)
  {}
  ~MuScleFitMuonSelector() {}

  //Method to get the muon after FSR (status 1 muon in PYTHIA6) starting from status 3 muon which is daughter of the Z
  const reco::Candidate* getStatus1Muon(const reco::Candidate* status3Muon);

  //Method to get the muon before FSR (status 3 muon in PYTHIA6) starting from status 3 muon which is daughter of the Z
  const reco::Candidate* getStatus3Muon(const reco::Candidate* status3Muon);
  
  /// Main method used to select muons of type specified by muonType_ from the collection specified by muonLabel_ and PATmuons_
  void selectMuons(const edm::Event & event, std::vector<MuScleFitMuon> & muons,
		   std::vector<GenMuonPair> & genPair,
		   std::vector<std::pair<lorentzVector,lorentzVector> > & simPair,
		   MuScleFitPlotter * plotter);


 protected:
  /// Apply the Onia cuts to select globalMuons
  bool selGlobalMuon(const pat::Muon* aMuon);
  /// Apply the Onia cuts to select trackerMuons
  bool selTrackerMuon(const pat::Muon* aMuon);

  // Generator and simulation level information
  GenMuonPair findGenMuFromRes( const reco::GenParticleCollection* genParticles);
  GenMuonPair findGenMuFromRes( const edm::HepMCProduct* evtMC );
  std::pair<lorentzVector, lorentzVector> findSimMuFromRes( const edm::Handle<edm::HepMCProduct> & evtMC,
							    const edm::Handle<edm::SimTrackContainer> & simTracks );
  void selectGeneratedMuons(const edm::Handle<pat::CompositeCandidateCollection> & collAll,
			    const std::vector<const pat::Muon*> & collMuSel,
			    std::vector<GenMuonPair> & genPair,
			    MuScleFitPlotter * plotter);
  void selectGenSimMuons(const edm::Event & event,
			 std::vector<GenMuonPair> & genPair,
			 std::vector<std::pair<lorentzVector,lorentzVector> > & simPair,
			 MuScleFitPlotter * plotter);
  // void selectGeneratedMuons(const edm::Event & event, std::vector<std::pair<lorentzVector,lorentzVector> > & genPair);
  void selectSimulatedMuons(const edm::Event & event,
			    const bool ifHepMC, edm::Handle<edm::HepMCProduct> evtMC,
			    std::vector<std::pair<lorentzVector,lorentzVector> > & simPair,
			    MuScleFitPlotter * plotter);

  /// Template function used to convert the muon collection to a vector of reco::LeafCandidate
  template<typename T>
  std::vector<MuScleFitMuon> fillMuonCollection( const std::vector<T>& tracks )
  {
    std::vector<MuScleFitMuon> muons;
    typename std::vector<T>::const_iterator track;
    for( track = tracks.begin(); track != tracks.end(); ++track ) {
      reco::Particle::LorentzVector mu;
      mu = reco::Particle::LorentzVector(track->px(),track->py(),track->pz(),
					 sqrt(track->p()*track->p() + mMu2));
      
      Double_t hitsTk(0), hitsMuon(0), ptError(0);
      if ( const reco::Muon* myMu = dynamic_cast<const reco::Muon*>(&(*track))  ){
	hitsTk =   myMu->innerTrack()->hitPattern().numberOfValidTrackerHits();
	hitsMuon = myMu->innerTrack()->hitPattern().numberOfValidMuonHits();
	ptError =  myMu->innerTrack()->ptError();
      }
      else if ( const pat::Muon* myMu = dynamic_cast<const pat::Muon*>(&(*track)) ) {
	hitsTk =   myMu->innerTrack()->hitPattern().numberOfValidTrackerHits();
	hitsMuon = myMu->innerTrack()->hitPattern().numberOfValidMuonHits();
	ptError =  myMu->innerTrack()->ptError();
      }
      else if (const reco::Track* myMu = dynamic_cast<const reco::Track*>(&(*track))){
	hitsTk =   myMu->hitPattern().numberOfValidTrackerHits();
	hitsMuon = myMu->hitPattern().numberOfValidMuonHits();
	ptError =  myMu->ptError();
      }
      
      MuScleFitMuon muon(mu,track->charge(),ptError,hitsTk,hitsMuon);

    if (debug_>0) {
      std::cout<<"[MuScleFitMuonSelector::fillMuonCollection] after MuScleFitMuon initialization"<<std::endl;
      std::cout<<"  muon = "<<muon<<std::endl;
    }

    muons.push_back(muon);
    }
    return muons;
  }

  /// Template function used to extract the selected muon type from the muon collection
  template<typename T>
  void takeSelectedMuonType(const T & muon, std::vector<reco::Track> & tracks)
  {
    // std::cout<<"muon "<<muon->isGlobalMuon()<<muon->isStandAloneMuon()<<muon->isTrackerMuon()<<std::endl;
    //NNBB: one muon can be of many kinds at once but with the muonType_ we are sure
    // to avoid double counting of the same muon
    if(muon->isGlobalMuon() && muonType_==1)
      tracks.push_back(*(muon->globalTrack()));
    else if(muon->isStandAloneMuon() && muonType_==2)
      tracks.push_back(*(muon->outerTrack()));
    else if(muon->isTrackerMuon() && muonType_==3)
      tracks.push_back(*(muon->innerTrack()));

    else if( muonType_ == 10 && !(muon->isStandAloneMuon()) ) //particular case!!
      tracks.push_back(*(muon->innerTrack()));
    else if( muonType_ == 11 && muon->isGlobalMuon() )
      tracks.push_back(*(muon->innerTrack()));
    else if( muonType_ == 13 && muon->isTrackerMuon() )
      tracks.push_back(*(muon->innerTrack()));
  }

  const edm::InputTag muonLabel_;
  const int muonType_;
  const bool PATmuons_;
  const std::vector<int> resfind_;
  const bool speedup_;
  const std::string genParticlesName_;
  const bool compareToSimTracks_;
  const edm::InputTag simTracksCollectionName_;
  const bool sherpa_;
  const bool debug_;
  static const double mMu2;
  static const unsigned int motherPdgIdArray[6];
};

#endif