PndTrackingQA.cxx

Go to the documentation of this file.

/*
 * PndTrackingQA.cxx
 *
 *  Created on: Aug 23, 2013
 *      Author: stockman
 */

#include <PndTrackingQA.h>
#include "PndTrack.h"
#include "PndMCTrack.h"
#include "FairHit.h"

ClassImp(PndTrackingQA);

PndTrackingQA::PndTrackingQA (TString trackBranchName, TString idealTrackName, Bool_t pndTrackData):
       fTrackBranchName(trackBranchName), fIdealTrackName(idealTrackName), fPndTrackOrTrackCand(pndTrackData), fPossibleTrack(0), fCleanFunctor(kFALSE), fNGhosts(0), fUseCorrectedSkewedHits(kFALSE), fVerbose(0)
{
        if(fPossibleTrack == 0){
                std::cout << "-I- PndTrackingQA::PndTrackingQA no PossibleTrackFunctor given. Taking Standard!" << std::endl;
                if (trackBranchName == "MVDTrack" ){
                        fPossibleTrack = new RiemannMvdSttGemFunctor();
                } else if (trackBranchName == "CombiTrackCand" ) {
                        fPossibleTrack = new OnlySttFunctor();
                } else {
                        fPossibleTrack = new StandardTrackFunctor();
                }
                fCleanFunctor = kTRUE;
        }
}

PndTrackingQA::PndTrackingQA (TString trackBranchName, TString idealTrackName, PndTrackFunctor* posTrack, Bool_t pndTrackData):
        fTrackBranchName(trackBranchName), fIdealTrackName(idealTrackName), fPndTrackOrTrackCand(pndTrackData), fPossibleTrack(posTrack), fCleanFunctor(kFALSE), fNGhosts(0), fUseCorrectedSkewedHits(kFALSE), fVerbose(0)
{
        if(fPossibleTrack == 0){
                std::cout << "-I- PndTrackingQA::PndTrackingQA no PossibleTrackFunctor given. Taking Standard!" << std::endl;
                if (trackBranchName == "MVDTrack" ){
                        fPossibleTrack = new RiemannMvdSttGemFunctor();
                } else if (trackBranchName == "CombiTrackCand" ) {
                        fPossibleTrack = new OnlySttFunctor();
                } else {
                        fPossibleTrack = new StandardTrackFunctor();
                }
                fCleanFunctor = kTRUE;
        }
}

PndTrackingQA::~PndTrackingQA()
{
        if (fCleanFunctor)
                delete (fPossibleTrack);
}

void PndTrackingQA::Init()
{
        ioman = FairRootManager::Instance();
        if (!ioman) {
                std::cout << "-E- PndTrackingQualityTask::Init: "
                                << "RootManager not instantiated!" << std::endl;
                return;
        }

        fTrack = (TClonesArray*) ioman->GetObject(fTrackBranchName);
        fMCTrack = (TClonesArray*) ioman->GetObject("MCTrack");
        fIdealTrack = (TClonesArray*) ioman->GetObject(fIdealTrackName);


        if (fBranchNames.size() == 0){
                AddHitsBranchName("MVDHitsPixel");
                AddHitsBranchName("MVDHitsStrip");
                AddHitsBranchName("STTHit");
                AddHitsBranchName("GEMHit");
                AddHitsBranchName("FTSHit");

        //      if (FairRootManager::Instance()->GetBranchId("CorrectedSkewedHits")  > 0){
        //              AddHitsBranchName("CorrectedSkewedHits");
        //      }
        }
        if (fVerbose > 0){
                std::cout << "-I- PndTrackingQA::Init: PossibleTrackFunctor: ";
                fPossibleTrack->Print();
        }
}

void PndTrackingQA::AnalyseEvent(TClonesArray *recoTrackInfo)
{
        FillMapTrackQualifikation();
        if (fVerbose > 2){
                std::cout << "PndTrackingQA::AnalyseEvent() Track quality map before analysis: " << std::endl << std::endl;
                PrintTrackQualityMap(kTRUE);
        }

        for (Int_t i = 0; i < fTrack->GetEntriesFast(); i++){
                if (fVerbose > 2){
                        std::cout << "----------------------------------" << std::endl;
                        std::cout << "Analyse Track: " << i << std::endl;
                }
                std::map<TString, FairMultiLinkedData> trackInfo;

                if (fPndTrackOrTrackCand){
                        PndTrack* myTrack = (PndTrack*)fTrack->At(i);
                        trackInfo = AnalyseTrackCand(myTrack->GetTrackCandPtr());
                } else {
                        PndTrackCand* myTrack = (PndTrackCand*)fTrack->At(i);
                        trackInfo = AnalyseTrackCand(myTrack);
                }
                Int_t mostProbableTrack = AnalyseTrackInfo(trackInfo, i);

                if (fVerbose > 1){
                        std::cout << "PndTrackingQA::AnalyseEvent Analyse track: " << i << std::endl;
                        std::cout << "mostProbableTrack " << mostProbableTrack << std::endl;
                }

                if (mostProbableTrack == -1) continue;

                fTrackIdMCId[i] = mostProbableTrack;
                CalcEfficiencies(mostProbableTrack, trackInfo);

                if (fMapTrackQualification.count(mostProbableTrack) > 0)
                        fMCTrackFound[mostProbableTrack]++;

                PndTrackingQualityRecoInfo recoinfo = GetRecoInfoFromRecoTrack(i, mostProbableTrack);
                int nof_asso_mctracks = trackInfo["AllHits"].GetNLinks();
                recoinfo.SetNofMCTracks(nof_asso_mctracks);
                int size = recoTrackInfo->GetEntriesFast();
                new((*recoTrackInfo)[size]) PndTrackingQualityRecoInfo(recoinfo);
        }

        for (std::map<Int_t, Int_t>::iterator iter = fMapTrackQualification.begin(); iter != fMapTrackQualification.end(); iter++) {

                if (iter->first > -1 && iter->second > 0) {
                        PndMCTrack* mcTrack = (PndMCTrack*) fMCTrack->At(iter->first);

                        if (fPndTrackOrTrackCand) {
                                PndTrack* myTrack = (PndTrack*) fTrack->At(fMCIdTrackId[iter->first]);
                                TVector3 mom(myTrack->GetParamFirst().GetPx(),
                                                myTrack->GetParamFirst().GetPy(),
                                                myTrack->GetParamFirst().GetPz());
                                TVector3 McMom(mcTrack->GetMomentum());

                                fMapPResolution[iter->first] = (mom.Mag() - McMom.Mag());
                                fMapP[iter->first] = mom;
                                fMapPtResolution[iter->first] = (mom.Pt() - McMom.Pt());
                                fMapPt[iter->first] = mom.Pt();
                                fMapPlResolution[iter->first] = (mom.Pz() - McMom.Pz());
                                fMapPl[iter->first] = mom.Pz();
                                fMapPResolutionRel[iter->first] = (mom.Mag() - McMom.Mag())
                                                / McMom.Mag();
                                fMapPtResolutionRel[iter->first] = (mom.Pt() - McMom.Pt())
                                                / McMom.Pt();
                                fMapPlResolutionRel[iter->first] = (mom.Pz() - McMom.Pz())
                                                / McMom.Pz();
                        }
                }

        }
        if (fVerbose > 2)
                std::cout << "AnalyseEvent End" << std::endl;
}

FairMultiLinkedData PndTrackingQA::GetMCInfoForBranch(TString branchName, PndTrackCand* trackCand)
{
        FairMultiLinkedData result;
        result.SetInsertHistory(kFALSE);
        //std::cout << "GetMCInforForBranch: " << branchName << std::endl;
        FairMultiLinkedData linksOfType = trackCand->GetLinksWithType(ioman->GetBranchId(branchName));
//      std::cout << "GetMCInforForBranch: LinksOfType " << branchName << " : " << linksOfType << std::endl;

        for (int j = 0; j < linksOfType.GetNLinks(); j++){
                FairMultiLinkedData_Interface* linkData = (FairMultiLinkedData_Interface*)FairRootManager::Instance()->GetCloneOfLinkData(linksOfType.GetLink(j));
                //std::cout << "CloneOfLinkData: " << *linkData << std::endl;
                if (linkData != 0){
                        FairMultiLinkedData linkDataType = linkData->GetLinksWithType(FairRootManager::Instance()->GetBranchId("MCTrack"));
                        linkDataType.SetAllWeights(1.);
                        result.AddLinks(linkDataType);
                        linkData->Delete();
                }
        }
        return result;
}

std::map<TString, FairMultiLinkedData> PndTrackingQA::AnalyseTrackCand(PndTrackCand* trackCand)
{
        std::map<TString, FairMultiLinkedData> trackInfo;

        if (fVerbose > 2) {
                std::cout << "PndTrackingQualityData::AnalyseTrackCand: TrackInfo" << std::endl;
                std::cout << *trackCand << std::endl;
        }
        for (size_t branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
//              std::cout << "BranchName: " << fBranchNames[branchIndex] << std::endl;
                trackInfo[fBranchNames[branchIndex]] = GetMCInfoForBranch(fBranchNames[branchIndex], trackCand);
                trackInfo["AllHits"].AddLinks(trackInfo[fBranchNames[branchIndex]]);
        }
        if (fVerbose > 2)
                PrintTrackInfo(trackInfo);
        return trackInfo;
}

Int_t PndTrackingQA::AnalyseTrackInfo(std::map<TString, FairMultiLinkedData>& trackInfo, Int_t trackId)
{

        Int_t mostProbableTrack = -1;   //the MCTrack most FairLinks of a TrackCand are pointing to
        if (fVerbose > 2){
                std::cout << "PndTrackingQA::AnalyseTrackInfo: NMCLinks: " << trackInfo["AllHits"].GetNLinks() << std::endl;
                PrintTrackDataSummary(trackInfo["AllHits"]);
                std::cout << "PndTrackingQA::AnalyseTrackInfo " << trackInfo["AllHits"] << std::endl;
        }

        std::vector<FairLink> sortedMCTracks = trackInfo["AllHits"].GetSortedMCTracks();

        if (sortedMCTracks.size() == 0) return mostProbableTrack; //returns -1; no MCTracks

        if (sortedMCTracks.size() == 1){
                mostProbableTrack = sortedMCTracks[0].GetIndex();
                if (fMCIdIdealTrackId.count(mostProbableTrack) == 0){
                        std::cout << "-W- PndTrackingQA::AnalyseTrackInfo fMCIdIdealTrackId does not contain mostProbableTrack" << mostProbableTrack << std::endl;
                        return -1;
                }
                if (!(fMCIdIdealTrackId[mostProbableTrack] < fIdealTrack->GetEntriesFast())){
                        std::cout << "-W- PndTrackingQA::AnalyseTrackInfo fMCIdIdealTrackId the ideal track is not in fIdealTrack. mostProbableTrack:" << mostProbableTrack << " idealTrack " << fMCIdIdealTrackId[mostProbableTrack] << std::endl;
                        return -1;
                }
                PndTrackCand* myIdealTrack = ((PndTrack*)fIdealTrack->At(fMCIdIdealTrackId[mostProbableTrack]))->GetTrackCandPtr();
                Int_t nMCHits = GetSumOfAllValidMCHits(myIdealTrack->GetPointerToLinks());              //get the number of hits which should be found

                if (nMCHits == sortedMCTracks[0].GetWeight()){
                        fMapTrackQualification[sortedMCTracks[0].GetIndex()] = qualityNumbers::kFullyFound;
                        fMCIdTrackId[mostProbableTrack] = trackId;
                } else {
                        if (fMapTrackQualification[sortedMCTracks[0].GetIndex()] != qualityNumbers::kFullyFound){
                                fMapTrackQualification[sortedMCTracks[0].GetIndex()] = qualityNumbers::kPartiallyFound;
                                fMCIdTrackId[mostProbableTrack] = trackId;
                        }
                }
        } else {
                Int_t highestCount = sortedMCTracks[0].GetWeight();
                mostProbableTrack = sortedMCTracks[0].GetIndex();
                Int_t allCounts = 0;
                for (size_t i = 0; i < sortedMCTracks.size(); i++){
                        allCounts += sortedMCTracks[i].GetWeight();
//                      if (trackInfo["AllHits"].GetLink(i).GetWeight() > highestCount){
//                              highestCount = trackInfo["AllHits"].GetLink(i).GetWeight();
//                              mostProbableTrack = trackInfo["AllHits"].GetLink(i).GetIndex();
//                      }
                }


                if ((Double_t)highestCount/(Double_t)allCounts > 0.7){
                        if (fMapTrackQualification[mostProbableTrack] != qualityNumbers::kFullyFound
                                        && fMapTrackQualification[mostProbableTrack] != qualityNumbers::kPartiallyFound){
                                fMapTrackQualification[mostProbableTrack] = qualityNumbers::kSpuriousFound;
                                fMCIdTrackId[mostProbableTrack] = trackId;
                        }
                }
                else {
                        fNGhosts++;
                }
        }
        if (fVerbose > 0){
                for (size_t j = 0; j < sortedMCTracks.size(); j++){
                        FairLink myLink = sortedMCTracks[j];
            if (fMCIdIdealTrackId.count(myLink.GetIndex()) > 0){
                PndTrackCand* myIdealTrack = ((PndTrack*)fIdealTrack->At(fMCIdIdealTrackId[myLink.GetIndex()]))->GetTrackCandPtr();
                if (fVerbose > 1) std::cout << "Ideal Tracking: Track " << myLink.GetIndex() << ": ";
                if (fVerbose > 1) PrintTrackDataSummary(*myIdealTrack->GetPointerToLinks());
            } else {
                std::cout << "Ideal Tracking: Track " << myLink.GetIndex() << " not available" << std::endl;
            }

                }
//              std::cout << "MostProbableTrack: " << mostProbableTrack << " Quality: " << fMapTrackQualification[mostProbableTrack] << std::endl;
//              std::cout << std::endl;
        }

        return mostProbableTrack;
}

void PndTrackingQA::FillMapTrackQualifikation()
{
        fMapTrackQualification.clear();
        fMapTrackMCStatus.clear();
        fMCIdIdealTrackId.clear();
//      std::cout << " fIdealTrack.size() " << fIdealTrack->GetEntriesFast() << std::endl;
        for (int i = 0; i < fIdealTrack->GetEntriesFast(); i++){
                PndTrackCand* idealTrackCand = ((PndTrack*)fIdealTrack->At(i))->GetTrackCandPtr();
//              std::cout << i << " : PndTrackingQA::FillMapTrackQualifikation: " << *idealTrackCand << std::endl;

                PndMCTrack* mcTrack = (PndMCTrack*)fMCTrack->At(idealTrackCand->getMcTrackId());


                fMCIdIdealTrackId[idealTrackCand->getMcTrackId()] = i;
                Bool_t primaryTrack = (mcTrack->GetMotherID() < 0);
                Bool_t atLeastThreeHits = kFALSE;
                Int_t nHits = 0;
//                      for (int branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
//                              TString branchName = fBranchNames[branchIndex];
//                              nHits += GetNIdealHits(i, branchName);
//                      }
                nHits += GetNIdealHits(*(idealTrackCand->GetPointerToLinks()), "MVDHitsPixel");
                nHits += GetNIdealHits(*(idealTrackCand->GetPointerToLinks()), "MVDHitsStrip");
                nHits += GetNIdealHits(*(idealTrackCand->GetPointerToLinks()), "STTHit");
                nHits += GetNIdealHits(*(idealTrackCand->GetPointerToLinks()), "GEMHit");
                nHits += GetNIdealHits(*(idealTrackCand->GetPointerToLinks()), "FTSHit");
                //std::cout << "FillMapTrackQualifikation: NHits: " << nHits << std::endl;

                if (nHits > 2) atLeastThreeHits = kTRUE;

                if (atLeastThreeHits) {
                        if (primaryTrack) {
                                fMapTrackQualification[idealTrackCand->getMcTrackId()] = qualityNumbers::kAtLeastThreePrim;
                        } else {
                                fMapTrackQualification[idealTrackCand->getMcTrackId()] = qualityNumbers::kAtLeastThreeSec;
                        }
                //      std::cout << "AtLeastThreeHits: Track "<< i << " status: " << fMapTrackQualification[i] << std::endl;
                }

                else if (primaryTrack){  //No hits for primary track in tracking detectors
                        fMapTrackQualification[idealTrackCand->getMcTrackId()] = qualityNumbers::kLessThanThreePrim;
                }

                PndTrackCand* entry = ((PndTrack*)fIdealTrack->At(i))->GetTrackCandPtr();
                if ((*fPossibleTrack)((FairMultiLinkedData*)entry->GetPointerToLinks(), primaryTrack))
                {
                        if (primaryTrack) {
                                fMapTrackQualification[idealTrackCand->getMcTrackId()] = qualityNumbers::kPossiblePrim;
                        } else {
                                fMapTrackQualification[idealTrackCand->getMcTrackId()] = qualityNumbers::kPossibleSec;
                        }
                }
        }
        if (fVerbose > 1){
                std::cout << "-I- PndMCTestPatternRecoQuality::FillMapTrackQualifikation:" << std::endl;
//              PrintTrackQualityMap();
        }
        fMapTrackMCStatus = fMapTrackQualification;
}

//Bool_t PndTrackingQA::PossibleTrack(FairMultiLinkedData& mcForward)
//{
//      Bool_t possibleTrack = kFALSE;
//
//      for (int i = 0; i < fBranchNames.size(); i++){
//              if (fBranchNames[i] == "MVDHitsPixel"){
//                      possibleTrack = possibleTrack | (mcForward.GetLinksWithType(ioman->GetBranchId("MVDHitsPixel")).GetNLinks() +
//                                                                                   mcForward.GetLinksWithType(ioman->GetBranchId("MVDHitsStrip")).GetNLinks() > 3);
//
//              }
//              if (fBranchNames[i] == "STTHit"){
//                      possibleTrack = possibleTrack | mcForward.GetLinksWithType(ioman->GetBranchId("STTHit")).GetNLinks() > 5;
//              }
//      }
//
//      return possibleTrack;
//}

Int_t PndTrackingQA::GetSumOfAllValidMCHits(FairMultiLinkedData* trackData)
{
        Int_t result = 0;
        for (size_t branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
                if (fBranchNames[branchIndex] == "GEMHit"){
                        FairMultiLinkedData gemHits = trackData->GetLinksWithType(ioman->GetBranchId("GEMPoint"));
                        result += gemHits.GetNLinks();
                } else {
                        result += trackData->GetLinksWithType(ioman->GetBranchId(fBranchNames[branchIndex])).GetNLinks();
                }
        }
        return result;
}


void PndTrackingQA::CalcEfficiencies(Int_t mostProbableTrack, std::map<TString, FairMultiLinkedData>& trackInfo)
{
        if (mostProbableTrack < 0) return;
        if (fVerbose > 2)
                std::cout << "PndTrackingQA::CalcEfficiencies() for mostProbableTrack " << mostProbableTrack
                                        <<  " count " << fMCIdIdealTrackId.count(mostProbableTrack) << std::endl;
        for (size_t branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
                if (fMCIdIdealTrackId.count(mostProbableTrack) > 0){
                        PndTrackCand* trackCand = ((PndTrack*)fIdealTrack->At(fMCIdIdealTrackId[mostProbableTrack]))->GetTrackCandPtr();
                        if (trackCand == 0) return;
                        Int_t nMcHits = GetNIdealHits(*trackCand->GetPointerToLinks(), fBranchNames[branchIndex]);
                        FairMultiLinkedData foundHits = trackInfo[fBranchNames[branchIndex]];
                        for (int i = 0; i < foundHits.GetNLinks(); i++){
                                if (foundHits.GetLink(i).GetIndex() == mostProbableTrack){
                                        Double_t nFoundHits = foundHits.GetLink(i).GetWeight();
                                        if ((nFoundHits/nMcHits) > fMapEfficiencies[mostProbableTrack][fBranchNames[branchIndex]].first){
                                                std::pair<Double_t, Int_t> result(nFoundHits/nMcHits, nMcHits);
                                                fMapEfficiencies[mostProbableTrack][fBranchNames[branchIndex]]=result;
                                        }
                                }
                        }
                }
        }
}

//Int_t PndTrackingQA::GetNIdealHits(Int_t trackId, TString branchName)
//{
//      PndMCEntry idealTrack = fIdealTrackCand.GetEntry(trackId);
//      return GetNIdealHits(idealTrack, branchName);
//
//}


Int_t PndTrackingQA::GetNIdealHits(FairMultiLinkedData& track, TString branchName)
{
        Int_t numberGemHits = 0;
        if (branchName == "GEMHit"){
                numberGemHits = track.GetLinksWithType(ioman->GetBranchId("GEMPoint")).GetNLinks();
                return numberGemHits;
        }
        return track.GetLinksWithType(ioman->GetBranchId(branchName)).GetNLinks();
}

void PndTrackingQA::PrintTrackDataSummary(FairMultiLinkedData& trackData, Bool_t detailedInfo)
{
        if (detailedInfo == kTRUE) std::cout << std::endl;
        for (size_t branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
                TString branchName = fBranchNames[branchIndex];
                if (GetNIdealHits(trackData, branchName) > 0){
                        std::cout << branchName << " " << GetNIdealHits(trackData, branchName);
                        if (detailedInfo == kTRUE){
                                std::cout << " : ";
                                if (trackData.GetLinksWithType(ioman->GetBranchId(branchName)).GetNLinks() > 0)
                                         std::cout << trackData.GetLinksWithType(ioman->GetBranchId(branchName));
                                std::cout << std::endl;
                        } else {
                                std::cout << " | ";
                        }
                }
        }
        std::cout << std::endl;
}

void PndTrackingQA::PrintTrackDataSummaryCompare(FairMultiLinkedData& recoTrackData, FairMultiLinkedData& idealTrackData)
{
        for (size_t branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
                TString branchName = fBranchNames[branchIndex];
                std::cout << branchName << " " << GetNIdealHits(recoTrackData, branchName);
                std::cout << "/" << GetNIdealHits(idealTrackData, branchName);
                        std::cout << " | ";
        }
        std::cout << std::endl;
}


void PndTrackingQA::PrintTrackQualityMap(Bool_t detailedInfo)
{
        for (std::map<Int_t, Int_t>::iterator iter = fMapTrackQualification.begin(); iter != fMapTrackQualification.end(); iter++){
                std::cout << "TrackID: " << iter->first << " MCQuality: "  << qualityNumbers::QualityNumberToString(fMapTrackMCStatus[iter->first]) << " Quality: ";
                if (iter->second < 0)
                        std::cout << " NOT FOUND ";
                else
                        std::cout << qualityNumbers::QualityNumberToString(iter->second);
                std::cout << " NFound: " << fMCTrackFound[iter->first];
                std::cout << " MCData: ";
                if (fMCIdIdealTrackId.count(iter->first) > 0){
                        PndTrackCand* idealTrackCand = ((PndTrack*)fIdealTrack->At(fMCIdIdealTrackId[iter->first]))->GetTrackCandPtr();
                        if (iter->second < 0)
                                PrintTrackDataSummary(*idealTrackCand->GetPointerToLinks(), detailedInfo);
                        else {
                                PndTrackCand* recoTrackCand = ((PndTrack*)fTrack->At(fMCIdTrackId[iter->first]))->GetTrackCandPtr();
                                if (recoTrackCand != 0 && idealTrackCand != 0)
                                        PrintTrackDataSummaryCompare(*recoTrackCand->GetPointerToLinks(), *idealTrackCand->GetPointerToLinks());
                        }
                }
        }
        std::cout << std::endl;
}

void PndTrackingQA::PrintTrackMCStatusMap()
{
        std::cout << "PrintTrackMCStatusMap: " << std::endl;
        for (std::map<Int_t, Int_t>::iterator iter = fMapTrackMCStatus.begin(); iter != fMapTrackMCStatus.end(); iter++){
                std::cout << "TrackID: " << iter->first << " Quality: "  << qualityNumbers::QualityNumberToString(iter->second) << " Found: " << fMCTrackFound[iter->first];
                std::cout << std::endl;
        }
        std::cout << std::endl;
}


void PndTrackingQA::PrintTrackInfo(std::map<TString, FairMultiLinkedData> info)
{
        std::cout << "TrackInfo: (MC-ID/NHits) : ";
        for (std::map<TString, FairMultiLinkedData>::iterator iter = info.begin(); iter != info.end(); iter++){
                std::cout << iter->first;
                for (int i = 0; i < iter->second.GetNLinks(); i++){
                        std::cout << " : (" << iter->second.GetLink(i).GetIndex() << "/" << iter->second.GetLink(i).GetWeight() << ")";
                }
                std::cout << " || ";
        }
        std::cout << std::endl;
}


Bool_t PndTrackingQA::IsBetterTrackExisting(Int_t& mcIndex,  int quality)
{
        if (fMapTrackQualification.count(mcIndex) == 1){
                if (fMapTrackQualification[mcIndex] > 0 && fMapTrackQualification[mcIndex] > quality) return true;
        }
        return false;
}

PndTrackingQualityRecoInfo PndTrackingQA::GetRecoInfoFromRecoTrack(Int_t trackId, Int_t mctrackId)
{
  PndTrackingQualityRecoInfo recoinfo(trackId);

  // CHECK this can be modified in the future, for now it is ok
  // mvd pix / mvd str / stt paral / stt skew / gem

  if (fVerbose > 2)
          std::cout << "PndTrackingQA::GetRecoInfoFromRecoTrack()" << std::endl;

    std::map<int, int> noftruehits;
        std::map<int, int> noffakehits;
        std::map<int, int> nofmissinghits;

        for (size_t branchIndex = 0; branchIndex < fBranchNames.size();
                        branchIndex++) {
                noftruehits.clear();
                noffakehits.clear();
                nofmissinghits.clear();

                // get the reco track...
                PndTrack *track = (PndTrack*) fTrack->At(trackId);
                // .. and the track cand
                PndTrackCand* trackcand = track->GetTrackCandPtr();

                recoinfo.SetPositionFirst(track->GetParamFirst().GetPosition());
                recoinfo.SetMomentumFirst(track->GetParamFirst().GetMomentum());

                recoinfo.SetPositionLast(track->GetParamLast().GetPosition());
                recoinfo.SetMomentumLast(track->GetParamLast().GetMomentum());

                recoinfo.SetCharge(track->GetParamFirst().GetQ());
                recoinfo.SetFlag(track->GetFlag());

                // get links associated to the reco track
                FairMultiLinkedData ptrlink = *trackcand->GetPointerToLinks();
                // get links corresponding to the hits of the specific detector
                FairMultiLinkedData links = ptrlink.GetLinksWithType(
                                ioman->GetBranchId(fBranchNames[branchIndex]));
                // get their number
                Int_t nHits = links.GetNLinks();
                if (fVerbose > 2){
                        if (nHits > 0)
                                std::cout << "----- reco track " << trackId << " (mc track " << mctrackId<< ") has " << nHits << " from " << fBranchNames[branchIndex] << std::endl;
                }
                // get mc track info from each hit
                for (int ihit = 0; ihit < nHits; ihit++) {
                        FairLink link = links.GetLink(ihit);
                        int assomctrack = -1;
                        // std::cout << "ihit " << ihit << " " << link.GetIndex() << " " << link.GetType() << " " << link.GetWeight() << std::endl;
                        FairHit * hit = (FairHit*) links.GetData(link);
                        if (!hit) {
                                std::cout << "ihit " << ihit << " " << link
                                                << " is FAKE" << std::endl;
                                //    std::cout << "No Obj Hit" << std::endl;

                                //
                                if (noffakehits.count(branchIndex) > 0)
                                        noffakehits[branchIndex]++;
                                else
                                        noffakehits[branchIndex] = 1; // if not there
                                continue;
                        }

                        // get links of the hit
                        FairMultiLinkedData hitlink = *hit->GetPointerToLinks();
                        // get the links corresponding to the mc track associated to the hit
                        FairMultiLinkedData mclinks = hitlink.GetLinksWithType(
                                        ioman->GetBranchId("MCTrack"));
                        //        std::cout << "hit " << ihit  << " belongs to " << mclinks.GetNLinks() << " mc tracks" << std::endl;
                        Bool_t isgood = kFALSE;
                        FairMultiLinkedData mvdstrhits = links.GetLinksWithType(
                                        FairRootManager::Instance()->GetBranchId("MVDHitsStrip"));
                        FairMultiLinkedData gemhits = links.GetLinksWithType(
                                        FairRootManager::Instance()->GetBranchId("GEMHit"));
                        if ((gemhits.GetNLinks() > 0 || mvdstrhits.GetNLinks() > 0) && mclinks.GetNLinks() > 1) {
                                isgood = kFALSE;
                        } else {
                                for (int imctrk = 0; imctrk < mclinks.GetNLinks(); imctrk++) {
                                        FairLink mclink = mclinks.GetLink(imctrk);
                                        assomctrack = mclink.GetIndex();
                                        //           std::cout << "imctrk " << imctrk << " " << mclink.GetIndex() << " " << mclink.GetType() << " " << mclink.GetWeight() << std::endl;
//        std::cout << "ihit " << ihit << " (hitid " << link.GetIndex() << ") belongs to MC track " << assomctrack << std::endl;
                                        if (assomctrack == mctrackId)
                                                isgood = kTRUE;
                                }
                        }

                        // if true
                        if (isgood == kTRUE) {
                                // if true and already there
                                if (noftruehits.count(branchIndex) > 0)
                                        noftruehits[branchIndex]++;
                                else
                                        noftruehits[branchIndex] = 1; // if not there
                        } else { // if not
                                if (noffakehits.count(branchIndex) > 0)
                                        noffakehits[branchIndex]++;
                                else
                                        noffakehits[branchIndex] = 1; // if not there
                        }

                }

                // retrieve the ideal track cand associated
                if (fMCIdIdealTrackId.count(mctrackId) > 0) {
                        int idealtrackid = fMCIdIdealTrackId[mctrackId];
                        PndTrack *idealtrack = (PndTrack*) fIdealTrack->At(idealtrackid);
                        PndTrackCand* idealtrackcand = idealtrack->GetTrackCandPtr();

                        Int_t nMcHits = GetNIdealHits(*idealtrackcand->GetPointerToLinks(),
                                        fBranchNames[branchIndex]);

                        nofmissinghits[branchIndex] = nMcHits - noftruehits[branchIndex];

//      std::cout << "**** ideal track " << idealtrackid << " has " << nMcHits << " from " << fBranchNames[branchIndex] << std::endl;

                }

//    std::cout << "===> BRANCH " << fBranchNames[branchIndex] << " of track " << trackId << " (mc track " << mctrackId << ") has " << noftruehits[branchIndex] << " true, " << noffakehits[branchIndex] << " fake and " << nofmissinghits[branchIndex] << " missing hits" << std::endl;

                if (fBranchNames[branchIndex] == "MVDHitsPixel") {
                        recoinfo.SetNofMvdPixelTrueHits(noftruehits[branchIndex]);
                        recoinfo.SetNofMvdPixelFakeHits(noffakehits[branchIndex]);
                        recoinfo.SetNofMvdPixelMissingHits(nofmissinghits[branchIndex]);
                } else if (fBranchNames[branchIndex] == "MVDHitsStrip") {
                        recoinfo.SetNofMvdStripTrueHits(noftruehits[branchIndex]);
                        recoinfo.SetNofMvdStripFakeHits(noffakehits[branchIndex]);
                        recoinfo.SetNofMvdStripMissingHits(nofmissinghits[branchIndex]);
                } else if (fBranchNames[branchIndex] == "STTHit") {
                        recoinfo.SetNofSttTrueHits(noftruehits[branchIndex]);
                        recoinfo.SetNofSttFakeHits(noffakehits[branchIndex]);
                        recoinfo.SetNofSttMissingHits(nofmissinghits[branchIndex]);
                } else if (fBranchNames[branchIndex] == "GEMHit") {
                        recoinfo.SetNofGemTrueHits(noftruehits[branchIndex]);
                        recoinfo.SetNofGemFakeHits(noffakehits[branchIndex]);
                        recoinfo.SetNofGemMissingHits(nofmissinghits[branchIndex]);
                } else if (fBranchNames[branchIndex] == "FTSHit") {
                        recoinfo.SetNofFtsTrueHits(noftruehits[branchIndex]);
                        recoinfo.SetNofFtsFakeHits(noffakehits[branchIndex]);
                        recoinfo.SetNofFtsMissingHits(nofmissinghits[branchIndex]);
                }

  }

  recoinfo.SetMCTrackID(mctrackId);
  return recoinfo;
}