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):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_kMuonSource(iConfig.getParameter<edm::InputTag>("MuonSourceTag")),
                                                         m_kTrackSource(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.getByLabel(m_kMuonSource, muons);
    iEvent.getByLabel(m_kTrackSource, 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));
}