Tau3MuReco.cc

Go to the documentation of this file.

#include <TMath.h>
 
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
#include "HeavyFlavorAnalysis/Skimming/interface/Combinatorics.h"
 
#include "HeavyFlavorAnalysis/Skimming/interface/Tau3MuReco.h"
 
Tau3MuReco::Tau3MuReco(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
    : m_kMatchingDeltaR(iConfig.getParameter<double>("RecoAnalysisMatchingDeltaR")),
      m_kMatchingPt(iConfig.getParameter<double>("RecoAnalysisMatchingPt")),
      m_kTauMassCut(iConfig.getParameter<double>("RecoAnalysisTauMassCut")),
      m_kTauMass(iConfig.getParameter<double>("RecoAnalysisTauMass")),
      m_kMuonMass(iConfig.getParameter<double>("RecoAnalysisMuonMass")),
      m_kMuonSourceToken(iC.consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("MuonSourceTag"))),
      m_kTrackSourceToken(iC.consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("TrackSourceTag"))) {}
 
Tau3MuReco::~Tau3MuReco() {}
 
bool Tau3MuReco::doTau3MuReco(const edm::Event& iEvent,
                              const edm::EventSetup& iSetup,
                              reco::MuonCollection* muonCollection,
                              reco::TrackCollection* trackCollection) {
  m_MuonCollection = muonCollection;
  m_TrackCollection = trackCollection;
 
  edm::Handle<reco::MuonCollection> muons;
  edm::Handle<reco::TrackCollection> tracks;
 
  reco::MuonCollection::const_iterator muon;
 
  iEvent.getByToken(m_kMuonSourceToken, muons);
  iEvent.getByToken(m_kTrackSourceToken, tracks);
 
  for (muon = muons->begin(); muon != muons->end(); ++muon) {
    m_TrackCollection->push_back(*(muon->track().get()));
    m_MuonCollection->push_back(*muon);
  }
 
  if (m_TrackCollection->size() > 3) {
    //find the right three ones coming from tau
    if (findCorrectPairing() == false) {  //maybe implement something like in ==3 (throw away ....)
      LogDebug("Tau3MuReco") << "Could not find correct combination!" << std::endl;
      return false;
    }
 
    return true;
  }
 
  if (m_TrackCollection->size() == 3) {
    if (fabs(getInvariantMass(m_TrackCollection, m_kMuonMass) - m_kTauMass) <= m_kTauMassCut)
      return true;
    else  //throw away one muon which don't match
      removeIncorrectMuon();
  }
 
  if (m_TrackCollection->size() == 2) {
    //search the third track
 
    //get 3rd muon canidate from tracks
    if (find3rdTrack(iEvent, iSetup, *(tracks.product())) == false) {
      LogDebug("Tau3MuReco") << "A 3rd Track can not be assigned!" << std::endl;
      return false;
    } else
      return true;
  }
 
  // cannot use this event, because less than 2 muons have been found
 
  LogDebug("Tau3MuReco") << "Not enough (" << m_TrackCollection->size() << ") muons found! Event skipped!" << std::endl;
 
  return false;
}
 
//private
bool Tau3MuReco::check4MuonTrack(const reco::Track& track) {
  reco::TrackCollection::const_iterator iter;
 
  for (iter = m_TrackCollection->begin(); iter != m_TrackCollection->end(); iter++) {
    //check if the track has the right charge
    //and check if dR is smaller than fMatchingDeltaR
    if ((*iter).charge() == track.charge() && getDeltaR(*iter, track) < m_kMatchingDeltaR &&
        fabs(((*iter).pt()) - (track.pt())) <= (track.pt() * m_kMatchingPt)) {
      LogDebug("Tau3MuReco") << "Found muon track in Tracks with DeltaR: " << getDeltaR(*iter, track)
                             << " Pt: " << track.pt() << " Muons Pt is: " << (*iter).pt() << std::endl;
      return true;
    }
  }
  return false;
}
 
bool Tau3MuReco::find3rdTrack(const edm::Event& iEvent,
                              const edm::EventSetup& iSetup,
                              const reco::TrackCollection& Tracks) {
  //check size of TrackVector (should be two!)
  if (m_TrackCollection->size() != 2)
    return false;
 
  //more then two tracks should be in the event
  if (Tracks.size() <= 2)
    return false;
 
  double SumDeltaR = 0;
 
  double MuonDeltaR = getDeltaR(m_TrackCollection->at(0), m_TrackCollection->at(1));
 
  //Loop overall tracks
 
  LogDebug("Tau3MuReco") << "Number of tracks found: " << Tracks.size() << std::endl;
 
  std::multimap<double, const reco::Track> TrackMultiMap;
 
  unsigned short muonCounter = 0;
 
  reco::TrackCollection::const_iterator track;
 
  for (track = Tracks.begin(); track != Tracks.end(); track++) {
    if (check4MuonTrack(*track)) {
      muonCounter++;
      continue;
    }
 
    SumDeltaR = MuonDeltaR;
 
    SumDeltaR += getDeltaR(m_TrackCollection->at(1), *track);
    SumDeltaR += getDeltaR(*track, m_TrackCollection->at(0));
 
    std::pair<double, const reco::Track> actTrack(SumDeltaR, *track);
 
    TrackMultiMap.insert(actTrack);
  }
 
  //two tracks should be clearly identified as muons by check4MuonTrack
  //else event is not useable
  if (muonCounter < 2) {
    LogDebug("Tau3MuReco") << "Not enough muons (" << muonCounter << ") assigned to a track! Event skipped!"
                           << std::endl;
    return false;
  }
 
  std::multimap<double, const reco::Track>::iterator it = TrackMultiMap.begin();
 
  if (it == TrackMultiMap.end()) {
    LogDebug("Tau3MuReco") << "Not enough tracks (0) left! Event skipped!" << std::endl;
    return false;
  }
 
  //get 2mu+track with minimal DeltaR Sum (MultiMaps are sorted)
  m_TrackCollection->push_back((*it).second);
 
  //and check charge of this track
  //and check invariant mass of this combination
  //and make a vertex fit
 
  char Charge = m_TrackCollection->at(0).charge() * m_TrackCollection->at(1).charge();
 
  unsigned int count = 0;
 
  //Charge > 0 means the two muons have same charge, so the third track has to have the opposit charge
  while ((Charge > 0 && ((*it).second).charge() == (m_TrackCollection->at(0)).charge()) ||
         fabs(getInvariantMass(m_TrackCollection) - m_kTauMass) > m_kTauMassCut) {
    count++;
 
    LogDebug("Tau3MuReco") << "Track canidate: " << count << std::endl;
 
    if (Charge > 0 && ((*it).second).charge() != (m_TrackCollection->at(0)).charge())
      LogDebug("Tau3MuReco") << "\tWrong charge!" << std::endl;
    LogDebug("Tau3MuReco") << "\tInvariant Mass deviation! " << fabs(getInvariantMass(m_TrackCollection) - m_kTauMass)
                           << std::endl;
    LogDebug("Tau3MuReco") << "\tTrack Pt: " << (*it).second.pt() << std::endl;
    LogDebug("Tau3MuReco") << "\tDelta R: " << (*it).first << std::endl;
    LogDebug("Tau3MuReco") << "\tChi2: " << ((*it).second).normalizedChi2() << std::endl;
 
    ++it;
 
    //was not the best canidate
    m_TrackCollection->pop_back();
 
    if (it == TrackMultiMap.end())
      return false;
 
    //get next to minimal (Delta R Sum) track
    m_TrackCollection->push_back((*it).second);
  }
 
  LogDebug("Tau3MuReco") << "Found corresponding 3rd track: " << std::endl;
  LogDebug("Tau3MuReco") << "Track canidate: " << count << std::endl;
  LogDebug("Tau3MuReco") << "\tInvariant Mass deviation! " << fabs(getInvariantMass(m_TrackCollection) - m_kTauMass)
                         << std::endl;
  LogDebug("Tau3MuReco") << "\tDelta R: " << (*it).first << std::endl;
  LogDebug("Tau3MuReco") << "\tNormChi2: " << ((*it).second).normalizedChi2() << std::endl;
 
  //choose this track, because it is the best canidate
  return true;
}
 
bool Tau3MuReco::findCorrectPairing() {
  Combinatorics myCombinatorics(m_TrackCollection->size(), 3);
 
  std::vector<std::vector<UInt_t> > CombinationVec = myCombinatorics.GetCombinations();
 
  std::vector<std::vector<UInt_t> >::iterator it = CombinationVec.begin();
 
  char Charge = 0;
 
  reco::TrackCollection tempTrackCollection;
  reco::MuonCollection tempMuonCollection;
 
  do {
    if (it == CombinationVec.end())
      return false;
 
    Charge = 0;
 
    tempMuonCollection.clear();
    tempTrackCollection.clear();
 
    for (UInt_t i = 0; i < (*it).size(); i++) {
      Charge += m_TrackCollection->at((*it).at(i)).charge();
      tempTrackCollection.push_back(m_TrackCollection->at((*it).at(i)));
      tempMuonCollection.push_back(m_MuonCollection->at((*it).at(i)));
    }
 
    LogDebug("Tau3MuReco") << "Charge is: " << (int)Charge << " Have to be -1 or 1!!!" << std::endl;
    LogDebug("Tau3MuReco") << "Invariant mass is: " << fabs(getInvariantMass(&tempTrackCollection) - m_kTauMass)
                           << " Have to be smaller than " << m_kTauMassCut << std::endl;
 
    it++;
  } while (abs(Charge) != 1 || fabs(getInvariantMass(&tempTrackCollection) - m_kTauMass) > m_kTauMassCut);
 
  *m_MuonCollection = tempMuonCollection;
  *m_TrackCollection = tempTrackCollection;
 
  return true;
}
 
bool Tau3MuReco::removeIncorrectMuon() {
  double deltaR12 = getDeltaR(m_TrackCollection->at(0), m_TrackCollection->at(1));
  double deltaR23 = getDeltaR(m_TrackCollection->at(1), m_TrackCollection->at(2));
  double deltaR31 = getDeltaR(m_TrackCollection->at(2), m_TrackCollection->at(0));
 
  //if DeltaR12 is the smallest, than the 3rd one seems to be wrong
  //if DeltaR23 is the smallest, than the 2nd one seems to be wrong
  //if DeltaR31 is the smallest, than the 1st one seems to be wrong
 
  unsigned char temp;
  double junk;
 
  deltaR12 < deltaR23 ? temp = 3 : temp = 1;
  deltaR12 < deltaR23 ? junk = deltaR12 : junk = deltaR23;
 
  if (deltaR31 < junk)
    temp = 2;
 
  m_TrackCollection->erase(m_TrackCollection->begin() + temp - 1);
 
  return true;
}
 
double Tau3MuReco::getInvariantMass(const reco::TrackCollection* tracks, const double MuonMass) {
  unsigned int numOfParticles = tracks->size();
 
  double SumPx = 0;
  double SumPy = 0;
  double SumPz = 0;
 
  double SumE = 0;
 
  for (unsigned int i = 0; i < numOfParticles; i++) {
    SumPx += tracks->at(i).px();
    SumPy += tracks->at(i).py();
    SumPz += tracks->at(i).pz();
 
    SumE += sqrt(pow(tracks->at(i).p(), 2) + pow(MuonMass, 2));
  }
 
  double invmass = sqrt(pow(SumE, 2) - pow(SumPx, 2) - pow(SumPy, 2) - pow(SumPz, 2));
 
  return invmass;
}
 
double Tau3MuReco::getDeltaR(const reco::Track& track1, const reco::Track& track2) {
  double dEta = track1.eta() - track2.eta();
  double dPhi = track1.phi() - track2.phi();
 
  while (dPhi >= TMath::Pi())
    dPhi -= (2.0 * TMath::Pi());
  while (dPhi < (-1.0 * TMath::Pi()))
    dPhi += (2.0 * TMath::Pi());
 
  return sqrt(pow(dEta, 2) + pow(dPhi, 2));
}