hparticlepairmaker.cc

Go to the documentation of this file.

#include "hparticlepairmaker.h"
#include "hparticlepairdraw.h"
#include "hcategorymanager.h"
#include "hphysicsconstants.h"
#include "hades.h"
#include "hrecevent.h"


#include "TCanvas.h"
#include "TH1F.h"
#include "TLatex.h"
#include "TEllipse.h"
#include "TLine.h"

// ROOT's IO and RTTI stuff is added here
ClassImp(HParticlePairMaker)

//_HADES_CLASS_DESCRIPTION
////////////////////////////////////////////////////////////////////////////////
//
//
// HParticlePairMaker
//
// Build vectors of freference candidates, fothers candidates and fpairs
// for the current event. As reference candidates marked as kIsLepton
// will be used.
//
// nectEvent() to be called inside eventloop after HParticleSorter selection.
// 1. clear all internal structures
// 2. fill vector freference by candidates which are flagged kIsLepton ( or kIsUsed if it is setup)
//    All other candidates (disregarding kIsLepton or kIsUsed) are added to vector fothers
//    To avoid unneeded pairs the user can set a userFilter by set setUserFilter(Bool_t (*userfilter)(HParticleCand*))
//    In this case only candidates passing this function (return kTRUE) will be taken into account
//    for freference and fothers.
// 3. created all pair combinations between freference x freference (no double counting) and freference x fothers.
//       If setDoSkippedFullCandPairs(kTRUE) is  set (default) in addition all pairs of full reconstructed
//       particles ( inner+outer MDC + META) of fOthers are build (ffullrecoOthers x ffullrecoOthers +
//       ffullrecoOthers x fnofullrecoOthers ). All pairs with full reco candidates of the event are build
//       this way for investigation, even the ones which are skipped from reference by any criteria or double hit
//       rejection.
//    a. pass the selectPID1 and selectPID2 functions
//       will get pid1 or pid2 assigned. Otherwise pid is -2(fake+), -1(fake-) (see HPhysicsConstants).
//       if pid2 fails it will get pid assigned opposite to reference cand in pairCase10.
//       fakes for pid2 will be assigned by polarity or opposite to reference cand.
//    b. eClosePairSelect flags will be set for each pair (see hparticledef.h).
//       This flags can be used later to classify the pairs.
//    c. pairs are inserted in fpairs vector. With static HParticlePair::setDoMomentumCorrection(Bool_t doit)
//       the enegryloss correction can be swithced on/off
//    d. secondary vertex and vertex cut vars are calulated. Which  primary ebevnt vertex
//       should be used has to be set by the user (default = kVertexParticle, see eVertex
//       in hparticledef.h). The user can provide an own event vertex per event by calling setVertex().
//
//
// Bool_t  HParticleTool::evalPairsFlags(UInt_t flag,HParticlePair& pair)
// can be used to classify the pair according to the predefined pair cases
// in ePairCase + eClosePairSelect (hparticledef.h). There are dedicated filter
// functions provided to select from the pairs vector
//
//BEGIN_HTML
// <img src="http://web-docs.gsi.de/~halo/docs/hydra/classDocumentation/docu_pics/hydra2/particle/lepton_pair_cases_1.png" alt="Lepton Pair Cases" width="750" >
// <img src="http://web-docs.gsi.de/~halo/docs/hydra/classDocumentation/docu_pics/hydra2/particle/lepton_pair_cases_2.png" alt="Lepton Pair Cases" width="750" >
//END_HTML
//   Fig 1. All predefined pairCases for Leptons
//
//BEGIN_HTML
// <img src="http://web-docs.gsi.de/~halo/docs/hydra/classDocumentation/docu_pics/hydra2/particle/hadron_pair_cases.png" alt="Hadron Pair Cases" width="750" >
//END_HTML
//   Fig 2. All predefined pairCases for Hadrons
//
//BEGIN_HTML
// <img src="http://web-docs.gsi.de/~halo/docs/hydra/classDocumentation/docu_pics/hydra2/particle/lepton_hadron_pair_cases.png" alt="Mixed Lepton Hadron Pair Cases" width="525" >
//END_HTML
//   Fig 2. All predefined pairCases for mixed Lepton Hadrons
//
    //
    //-----------------------------------------------------------------------------
    // USAGE:
    // #include  "hparticlepairmaker.h"
    // #include  "hparticledef.h"
    // #include  "hparticletool.h"
    // #include  <vector>
    // using namespace std;
    // ....
    //
    //
    // {
    // // setup HLoop etc ....
    //
    // ....
    //  HEnergyLossCorrPar enLossCorr;
    //  enLossCorr.setDefaultPar("apr12"); //  "jan04","nov02","aug04","apr12"
    //
    // //---------------------------------------------
    // // HPARTICLEPAIRMAKER SETUP :
    // HParticlePairMaker pairmaker;
    // // pairmaker.setPIDs(2,3,51) // default
    // // pairmaker.setPIDsSelection(mySelectionPID1,mySelectionPID2) // set your own selection function for pid1+pid2 (signature :  Bool_t myfunction(HParticleCand*) )
    // // pairmaker.setUseLeptons(kTRUE);// kTRUE : use kIsLepton, kFALSE: kIsUsed for selection of reference  (default kTRUE)
    // // pairmake.setRequireRich(kFALSE); // do not ask for rich Hit in selection functions ( default = kTRUE)
    // // pairmaker.setVertexCase(Particle::kVertexParticle); // select which event vertex to use (default = kVertexParticle)
    // // HParticlePair::setDoMomentumCorrection(kTRUE); // default : kTRUE
    // //---------------------------------------------
    //
    // vector<HParticlePair*> pairs;
    //
    // for (Int_t i=0; i < entries; i++) {
    //    Int_t nbytes =  loop.nextEvent(i);             // get next event. categories will be cleared before
    //    if(nbytes <= 0) { cout<<nbytes<<endl; break; } // last event reached
    //
    //
    //    //--------------------------------------------------------------------------
    //    // HPARTICLETRACKSORTER SETUP
    //    sorter.cleanUp();
    //    sorter.resetFlags(kTRUE,kTRUE,kTRUE,kTRUE);     // reset all flags for flags (0-28) ,reject,used,lepton
    //    Int_t nCandLep     = sorter.fill(HParticleTrackSorter::selectLeptons);   // fill only good leptons
    //    Int_t nCandLepBest = sorter.selectBest(HParticleTrackSorter::kIsBestRKRKMETA,HParticleTrackSorter::kIsLepton);
    //    Int_t nCandHad     = sorter.fill(HParticleTrackSorter::selectHadrons);   // fill only good hadrons (already marked good leptons will be skipped)
    //    Int_t nCandHadBest = sorter.selectBest(HParticleTrackSorter::kIsBestRKRKMETA,HParticleTrackSorter::kIsHadron);
    //    //--------------------------------------------------------------------------
    //
    //
    //    //###########################################################################
    //    // HPARTICLEPAIRMAKER ACTION:
    //    pairmaker.nextEvent();  // call it after your HParticleSorter action!
    //
    //
    //
    //    //--------------------------------------------------------------------------
    //    // example  1:
    //    // fill list of pairs for pairCase1 and require
    //    // fitted inner+out+RK for second leg in addition
    //    pairmaker.filterPairsVector(pairs,kPairCase1|kFittedInnerMDC2|kFittedInnerMDC2|kRK2);
    //
    //    for(UInt_t j = 0; j < pairs.size(); j++ ){ // print the pairs
    //       pairs[j]->print();
    //    }
    //    //--------------------------------------------------------------------------
    //
    //
    //    //--------------------------------------------------------------------------
    //    // example 2:
    //    // get list of all pairs and loop over the list
    //    // plot different pairCases
    //
    //    vector<HParticlePair>& pairsVec =  pairmaker.getPairsVector();
    //    for(UInt_t l = 0; l < pairsVec.size(); l++ ){
    //        HParticlePair& pair = pairsVec[l];
    //
    //        if(HParticleTool::evalPairsFlags(kPairCase1,pair)){
    //          cout<<"case1"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase2,pair)){
    //          cout<<"case2"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase3,pair)){
    //          cout<<"case3"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase4,pair)){
    //          cout<<"case4"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase5,pair)){
    //          cout<<"case5"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase6,pair)){
    //          cout<<"case6"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase7,pair)){
    //          cout<<"case7"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase8,pair)){
    //          cout<<"case8"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase9,pair)){
    //          cout<<"case9"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase10,pair)){
    //          cout<<"case10"<<endl;
    //        } else if(HParticleTool::evalPairsFlags(kPairCase11,pair)){
    //          cout<<"case11"<<endl;
    //        }
    //    }
    //    //--------------------------------------------------------------------------
    //    //###########################################################################
    //   }  // end event loop
    // }
////////////////////////////////////////////////////////////////////////////////

Bool_t HParticlePairMaker::frequireRich = kTRUE;

HParticlePairMaker::HParticlePairMaker()
{
    fPID1      = HPhysicsConstants::pid("e+");
    fPID2      = HPhysicsConstants::pid("e-");
    fMotherPID = HPhysicsConstants::pid("dilepton");
    fselectPID1 = 0;
    fselectPID2 = 0;
    fuserFilter = 0;
    fuse_kIsLepton = kTRUE;
    fVertexCase = kVertexParticle ;
    fdoSkippedFullCandPairs = kTRUE;
    fpairs    .resize(10000);
    fothers   .resize(1000);
    ffullrecoOthers   .resize(100);
    fnofullrecoOthers .resize(1000);
    freference.resize(100);


    fCaseVec.push_back(kPairCase1);
    fCaseVec.push_back(kPairCase2);
    fCaseVec.push_back(kPairCase3);
    fCaseVec.push_back(kPairCase4);
    fCaseVec.push_back(kPairCase5);
    fCaseVec.push_back(kPairCase6);
    fCaseVec.push_back(kPairCase7);
    fCaseVec.push_back(kPairCase8);
    fCaseVec.push_back(kPairCase9);
    fCaseVec.push_back(kPairCase10);
    fCaseVec.push_back(kPairCase11);
    fCaseVec.push_back(kPairCase12);
    fCaseVec.push_back(kPairCase13);
    fCaseVec.push_back(kPairCase14);
    fCaseVec.push_back(kPairCase15);

    for(UInt_t i=0;i<fCaseVec.size();i++) fCaseCt.push_back(0);
    richCandCt = 0 ;
}

HParticlePairMaker::~HParticlePairMaker()
{
     clearVectors();
}

Bool_t HParticlePairMaker::selectPos(HParticleCand* cand)
{
    // standard function for positrons
    // require richmatch (broad window,  frequireRich =kTRUE)
    // Select positive polatity

    if(getRequireRich() && cand->getRichInd() == -1) return kFALSE;

    if( cand->getCharge() > 0 ) return kTRUE;
    else                        return kFALSE;
}

Bool_t HParticlePairMaker::selectNeg(HParticleCand* cand)
{
    // standard function for electrons
    // require richmatch (broad window,  frequireRich =kTRUE)
    // Select negative polatity

    if(getRequireRich() && cand->getRichInd() == -1) return kFALSE;

    if( cand->getCharge() < 0 ) return kTRUE;
    else                        return kFALSE;
}

void HParticlePairMaker::bookHits(HParticleCand* cand1)
{
    // add candidate hits to maps

    if(cand1->getRichInd() !=-1 ) {
        if(mRichtoCand.find(cand1->getRichInd()) == mRichtoCand.end()){
            vector<HParticleCand*> v;
            v.push_back(cand1);
            mRichtoCand[cand1->getRichInd()] = v;
        } else {
            mRichtoCand[cand1->getRichInd()].push_back(cand1);
        }
    }

    if(cand1->getInnerSegInd() !=-1 ) {
        if(mInnerMdctoCand.find(cand1->getInnerSegInd()) == mInnerMdctoCand.end()){
            vector<HParticleCand*> v;
            v.push_back(cand1);
            mInnerMdctoCand[cand1->getInnerSegInd()] = v;
        } else {
            mInnerMdctoCand[cand1->getInnerSegInd()].push_back(cand1);
        }
    }

    if(cand1->getOuterSegInd() !=-1 ) {
        if(mOuterMdctoCand.find(cand1->getOuterSegInd()) == mOuterMdctoCand.end()){
            vector<HParticleCand*> v;
            v.push_back(cand1);
            mOuterMdctoCand[cand1->getOuterSegInd()] = v;
        } else {
            mOuterMdctoCand[cand1->getOuterSegInd()].push_back(cand1);
        }
    }

    if(cand1->getSystemUsed() != -1){
        Int_t meta = cand1->getMetaHitInd();
        if(cand1->isTofHitUsed()){
            if(mTofHittoCand.find(meta) == mTofHittoCand.end()){
                vector<HParticleCand*> v;
                v.push_back(cand1);
                mTofHittoCand[meta] = v;
            } else {
                mTofHittoCand[meta].push_back(cand1);
            }
        } else if (cand1->isTofClstUsed()){
            if(mTofClsttoCand.find(meta) == mTofClsttoCand.end()){
                vector<HParticleCand*> v;
                v.push_back(cand1);
                mTofClsttoCand[meta] = v;
            } else {
                mTofClsttoCand[meta].push_back(cand1);
            }
        } else if (cand1->isRpcClstUsed()){
            if(mRpcClsttoCand.find(meta) == mRpcClsttoCand.end()){
                vector<HParticleCand*> v;
                v.push_back(cand1);
                mRpcClsttoCand[meta] = v;
            } else {
                mRpcClsttoCand[meta].push_back(cand1);
            }
        } else if (cand1->isShowerUsed()){
            if(mShowertoCand.find(meta) == mShowertoCand.end()){
                vector<HParticleCand*> v;
                v.push_back(cand1);
                mShowertoCand[meta] = v;
            } else {
                mShowertoCand[meta].push_back(cand1);
            }
        } else if (cand1->isEmcUsed()){
            if(mEmctoCand.find(meta) == mEmctoCand.end()){
                vector<HParticleCand*> v;
                v.push_back(cand1);
                mEmctoCand[meta] = v;
            } else {
                mEmctoCand[meta].push_back(cand1);
            }
        }
    }
}
void HParticlePairMaker::selectPID(HParticleCand* cand1,Int_t& pid1,Bool_t warn)
{
    if(       (*fselectPID1)(cand1) ) { pid1 = fPID1;}    // one of the 2 functions should set the PID1 , take care about others!
    else if  ((*fselectPID2)(cand1) ) { pid1 = fPID2;}
    else     {
        if(warn){
            Warning("nextEvent()","Reference particle not labeled with PID! Should not happen. Check your selection functions. frequireRich = kTRUE but kIsLepton is not used as referencce this might cause empty pids.");
            cand1->print();
        }
    }
}

void HParticlePairMaker::nextEvent()
{
    // to be called inside eventloop after HParticleSorter selection.
    // 1. clear all internal structures
    // 2. fill vector freference by candidates which are flagged kIsLepton
    //    All other candidates are added to vector fothers
    // 3. created all pair combinations between freference and fothers which
    //    a. pass the selectPID1 and selectPID2 functions
    //       will get pid1 or pid2 assigned. Otherwise pid is -2(fake+), -1(fake-) (see HPhysicsConstants).
    //       if pid2 fails it will get pid assigend opposite to reference cand in pairCase10.
    //       fakes for pid2 will be assigned by polarity or opposite to reference cand.
    //    b. eClosePairSelect flags will be set for each pair.
    //       This flags can be used later to classify the pairs.
    //    c. pairs are inserted in fpairs vector.


    clearVectors();

    if(!fselectPID1) fselectPID1 = HParticlePairMaker::selectPos;
    if(!fselectPID2) fselectPID2 = HParticlePairMaker::selectNeg;

    HRecEvent* event = (HRecEvent*)gHades->getCurrentEvent();
    if(event){
        if(fVertexCase != kVertexUser){
            fVertex = HParticleTool::getGlobalVertex(fVertexCase,kFALSE);
        }

        HCategory* candCat = (HCategory*) event->getCategory(catParticleCand);
        if(candCat){
            UInt_t n = candCat->getEntries();
            HParticleCand* cand1 = 0 ;
            HParticleCand* cand2 = 0 ;
            //-------------------------------------------------------
            // fill lists for reference and others
            for(UInt_t i=0; i < n; i++){
                cand1 = HCategoryManager::getObject(cand1,candCat,i);

                //-------------------------------------------------------
                //  book the hits -> candiate lists
                bookHits(cand1);
                //-------------------------------------------------------

                if(cand1){

                    if(fuserFilter && !(*fuserFilter)(cand1)) continue;


                    if(( fuse_kIsLepton && cand1->isFlagBit(Particle::kIsLepton)) ||
                       (!fuse_kIsLepton && cand1->isFlagBit(Particle::kIsUsed))
                      ){
                        freference.push_back(cand1);
                    } else {
                        fothers.push_back(cand1);
                    }
                }
            } // end fill list loop
            //-------------------------------------------------------


            //-------------------------------------------------------
            //  create pair combinations
            Int_t pid1,pid2;
            pid1 = pid2 = -10;

            //-------------------------------------------------------
            //  first do good pairs from reference x reference candidates

            n = freference.size();
            if(n > 1){
                for(UInt_t i = 0 ; i < n ; i++){
                    cand1 = freference[i];
                    pid1 = -10;

                    selectPID(cand1,pid1,kTRUE);

                    for(UInt_t j = i+1 ; j < n; j++){
                        cand2 = freference[j];
                        if(cand2->getRichInd()!=-1) richCandCt++;
                        UInt_t flag = 0;
                        HParticleTool::setPairFlags(flag,cand2,cand1);

                        if(!fuse_kIsLepton) flag=flag|kNoUseRICH;

                        pid2 = -10;

                        selectPID(cand2,pid2,kTRUE);

                        HParticlePair pair;
                        pair.setPair(cand1,pid1,cand2,pid2,fMotherPID,flag,fVertex);
                        fpairs.push_back(pair);

                    } // end loop others
                } // end loop reference x reference
            }
            //-------------------------------------------------------

            //-------------------------------------------------------
            // now do all pairs reference x other candidates
            for(UInt_t i = 0 ; i < freference.size(); i++){
                cand1 = freference[i];
                pid1 = -10;

                selectPID(cand1,pid1,kTRUE);

                for(UInt_t j = 0 ; j < fothers.size(); j++){
                    cand2 = fothers[j];
                    if(cand2->getRichInd()!=-1) richCandCt++;
                    UInt_t flag = 0;
                    HParticleTool::setPairFlags(flag,cand2,cand1);

                    if(!fuse_kIsLepton) flag=flag|kNoUseRICH;

                    pid2 = -10;
                    if(        pid1 == fPID2 && (*fselectPID1)(cand2) ) { pid2 = fPID1;}  // PID2 could fail no polarity is defined
                    else if  ( pid1 == fPID1 && (*fselectPID2)(cand2) ) { pid2 = fPID2;}

                    if(pid2 == -10){

                        if((flag&kNoOuterMDC2) == kNoOuterMDC2) {
                            // in this case polarity is not defined
                            // selectPID will fail, but we want to set
                            // pid2 in this case

                            if(pid1 == fPID1) pid2 = fPID2;
                            else              pid2 = fPID1;
                        }
                    }

                    if(pid2 == -10) { // all cases with no RICH
                        if(cand2->getCharge() != 0){
                            if(cand2->getCharge() > 0) pid2 =  HPhysicsConstants::pid("fake+");
                            else                       pid2 =  HPhysicsConstants::pid("fake-");
                        } else {
                            if(HPhysicsConstants::charge(pid1) < 0 ) pid2 =  HPhysicsConstants::pid("fake+");
                            else                                     pid2 =  HPhysicsConstants::pid("fake-");
                        }

                    }

                    HParticlePair pair;
                    pair.setPair(cand1,pid1,cand2,pid2,fMotherPID,flag,fVertex);
                    fpairs.push_back(pair);
                    //#define dbug
#ifdef  dbug
                    Int_t ct = 0 ;
                    for(UInt_t k = 0; k < fCaseVec.size(); k ++){   // fill statistics for pair cases
                        if(HParticleTool::evalPairsFlags(fCaseVec[k],pair) || cand2->getRichInd()==-1) {
                            good2=kTRUE;
                            if(cand2->getRichInd()!=-1) ct++;
                        }
                    }
                    if(!good2 || ct>1) pair.print();
#endif
                } // end loop others
            } // end loop reference
            //-------------------------------------------------------




            if(fdoSkippedFullCandPairs)
            {
                
                //-------------------------------------------------------
                // now do all pairs full reco of others x no full reco other candidates and fullreco others X full reco others
                for(UInt_t j = 0 ; j < fothers.size(); j++){
                    HParticleCand* c = fothers[j];
                    if(c->isFlagAND(4,                                 // Fully reconstructed
                                    Particle::kIsAcceptedHitInnerMDC,
                                    Particle::kIsAcceptedHitOuterMDC,
                                    Particle::kIsAcceptedHitMETA,
                                    Particle::kIsAcceptedRK )){
                          ffullrecoOthers.push_back(c);
                    } else {
                          fnofullrecoOthers.push_back(c);
                    }
                }

                //-------------------------------------------------------
                //fullreco others X full reco others
                for(UInt_t i = 0 ; i < ffullrecoOthers.size(); i++){
                    cand1 = ffullrecoOthers[i];
                    pid1 = -10;

                    selectPID(cand1,pid1,kTRUE);

                    for(UInt_t j = i+1 ; j < ffullrecoOthers.size(); j++){
                        cand2 = ffullrecoOthers[j];
                        UInt_t flag = 0;
                        HParticleTool::setPairFlags(flag,cand2,cand1);

                        if(!fuse_kIsLepton) flag=flag|kNoUseRICH;

                        pid2 = -10;
                        if(        pid1 == fPID2 && (*fselectPID1)(cand2) ) { pid2 = fPID1;}  // PID2 could fail no polarity is defined
                        else if  ( pid1 == fPID1 && (*fselectPID2)(cand2) ) { pid2 = fPID2;}

                        if(pid2 == -10){

                            if((flag&kNoOuterMDC2) == kNoOuterMDC2) {
                                // in this case polarity is not defined
                                // selectPID will fail, but we want to set
                                // pid2 in this case

                                if(pid1 == fPID1) pid2 = fPID2;
                                else              pid2 = fPID1;
                            }
                        }

                        if(pid2 == -10) { // all cases with no RICH
                            if(cand2->getCharge() != 0){
                                if(cand2->getCharge() > 0) pid2 =  HPhysicsConstants::pid("fake+");
                                else                       pid2 =  HPhysicsConstants::pid("fake-");
                            } else {
                                if(HPhysicsConstants::charge(pid1) < 0 ) pid2 =  HPhysicsConstants::pid("fake+");
                                else                                     pid2 =  HPhysicsConstants::pid("fake-");
                            }

                        }

                        HParticlePair pair;
                        pair.setPair(cand1,pid1,cand2,pid2,fMotherPID,flag,fVertex);
                        fpairs.push_back(pair);

                    } // end loop fullreco
                } // end loop fullreco
                //-------------------------------------------------------

                //-------------------------------------------------------
                //fullreco others X no full reco others
                for(UInt_t i = 0 ; i < ffullrecoOthers.size(); i++){
                    cand1 = ffullrecoOthers[i];
                    pid1 = -10;

                    selectPID(cand1,pid1,kTRUE);

                    for(UInt_t j = 0 ; j < fnofullrecoOthers.size(); j++){
                        cand2 = fnofullrecoOthers[j];
                        UInt_t flag = 0;
                        HParticleTool::setPairFlags(flag,cand2,cand1);

                        if(!fuse_kIsLepton) flag=flag|kNoUseRICH;

                        pid2 = -10;
                        if(        pid1 == fPID2 && (*fselectPID1)(cand2) ) { pid2 = fPID1;}  // PID2 could fail no polarity is defined
                        else if  ( pid1 == fPID1 && (*fselectPID2)(cand2) ) { pid2 = fPID2;}

                        if(pid2 == -10){

                            if((flag&kNoOuterMDC2) == kNoOuterMDC2) {
                                // in this case polarity is not defined
                                // selectPID will fail, but we want to set
                                // pid2 in this case

                                if(pid1 == fPID1) pid2 = fPID2;
                                else              pid2 = fPID1;
                            }
                        }

                        if(pid2 == -10) { // all cases with no RICH
                            if(cand2->getCharge() != 0){
                                if(cand2->getCharge() > 0) pid2 =  HPhysicsConstants::pid("fake+");
                                else                       pid2 =  HPhysicsConstants::pid("fake-");
                            } else {
                                if(HPhysicsConstants::charge(pid1) < 0 ) pid2 =  HPhysicsConstants::pid("fake+");
                                else                                     pid2 =  HPhysicsConstants::pid("fake-");
                            }

                        }

                        HParticlePair pair;
                        pair.setPair(cand1,pid1,cand2,pid2,fMotherPID,flag,fVertex);
                        fpairs.push_back(pair);

                    } // end loop fullnoreco
                } // end loop fullreco
                //-------------------------------------------------------

            }  // end fdoSkippedFullCandPairs


            //-------------------------------------------------------
            //  book the candidate -> pair lists
            for(UInt_t i = 0; i < fpairs.size(); i++){
                HParticlePair* pair = &(fpairs[i]);
                cand1 = dynamic_cast<HParticleCand*>(pair->getCand(0));
                cand2 = dynamic_cast<HParticleCand*>(pair->getCand(1));


                //-------------------------------------------------------
                // fill statistics for pair cases
                for(UInt_t k = 0; k < fCaseVec.size(); k ++){
                    if(HParticleTool::evalPairsFlags(fCaseVec[k],*pair)) {
                        fCaseCt[k]++;
                    }
                }
                //-------------------------------------------------------


                if(mCandtoPair.find(cand1) == mCandtoPair.end()){
                    vector<HParticlePair*> pairs;
                    pairs.push_back(pair);
                    mCandtoPair[cand1] = pairs;
                } else {
                    mCandtoPair[cand1].push_back(pair);
                }
                if(mCandtoPair.find(cand2) == mCandtoPair.end()){
                    vector<HParticlePair*> pairs;
                    pairs.push_back(pair);
                    mCandtoPair[cand2] = pairs;
                } else {
                    mCandtoPair[cand2].push_back(pair);
                }
            }
            //-------------------------------------------------------

        }

    }
}

void HParticlePairMaker::clearVectors()
{
    // clear all internal vectors etc
    freference     .clear();
    fothers        .clear();
    fpairs         .clear();
    ffullrecoOthers.clear();
    fnofullrecoOthers.clear();

    mTofHittoCand  .clear();
    mTofClsttoCand .clear();
    mRpcClsttoCand .clear();
    mShowertoCand  .clear();
    mEmctoCand     .clear();
    mInnerMdctoCand.clear();
    mOuterMdctoCand.clear();
    mRichtoCand    .clear();
    mCandtoPair    .clear();
}

void HParticlePairMaker::filterPairsVector(vector<HParticlePair*>& filterpairs,UInt_t flag)
{
    // fill all pairs which fullfill the flag to filterspairs.
    // filterpairs will be cleared automatically before filling.
    filterpairs.clear();
    for(UInt_t i = 0; i < fpairs.size(); i++){
        HParticlePair& pair = fpairs[i];
        if(HParticleTool::evalPairsFlags(flag,pair)){
            filterpairs.push_back(&pair);
        }
    }
}
void HParticlePairMaker::filterPairsVector(vector<HParticlePair*>& filterpairs,vector<UInt_t>& flags)
{
    // fill all pairs which fullfill the flags to filterspairs.
    // filterpairs will be cleared automatically before filling.
    filterpairs.clear();
    for(UInt_t i = 0; i < fpairs.size(); i++){
        HParticlePair& pair = fpairs[i];
        if(HParticleTool::evalPairsFlags(flags,pair)){
            filterpairs.push_back(&pair);
        }
    }
}

Int_t HParticlePairMaker::filterCandidates(HVirtualCand* cand,vector<HVirtualCand*>& candidates,UInt_t flag,Float_t oAngle)
{
    // vector candidates with all candidates which have been combined into pairs.
    // if flag == 0 no filtering will be done. Filtering takes place on pairs.
    // Flag can be used to filter for candidates which share 1 or more detector
    // hits with the candidate (see flags eClosePairSelect+ ePairCase (hparticledef.h)).
    // In addition a filter of opening Angle of the track can be applied (oAngle < 0 does no selection)
    // The function returns the number of found candidates or -1 if the candidate
    // has not been found at all (should not happen).

    candidates.clear();

    map<HVirtualCand*,vector<HParticlePair*> >::iterator it = mCandtoPair.find(cand);
    if(it != mCandtoPair.end()){
        vector<HParticlePair*>& v = it->second;

        for(UInt_t i = 0 ; i < v.size(); i++){
            HParticlePair& pair = *(v[i]);
            if(HParticleTool::evalPairsFlags(flag,pair) && (oAngle < 0 || pair.getOpeningAngle() < oAngle )){

                if(pair.getCand(0) == cand) candidates.push_back(pair.getCand(1));
                else                        candidates.push_back(pair.getCand(0));
            }
        }

       return candidates.size();
    }
    return -1;
}

Int_t HParticlePairMaker::filterCandidates (HVirtualCand* cand,vector<HParticlePair*>& filterpairs,UInt_t flag,Float_t oAngle)
{
    // Fills vector filterpairs with all pairs which share cand.
    // if flag == 0 no filtering will be done. Filtering takes place on pairs.
    // Flag can be used to filter for candidates which share 1 or more detector
    // hits with the candidate (see flags eClosePairSelect+ ePairCase (hparticledef.h)).
    // In addition a filter of opening Angle of the track can be applied (oAngle < 0 does no selection)
    // The function returns the number of found pairs or -1 if the candidate
    // has not been found at all (should not happen).

    filterpairs.clear();
    map<HVirtualCand*,vector<HParticlePair*> >::iterator it = mCandtoPair.find(cand);
    if(it != mCandtoPair.end()){
        vector<HParticlePair*>& v = it->second;
        for(UInt_t i = 0 ; i < v.size(); i++){

            if(HParticleTool::evalPairsFlags(flag,*(v[i])) && ( oAngle < 0 || v[i]->getOpeningAngle() < oAngle) ){
                filterpairs.push_back(v[i]);
            }
        }

    } else return -1;

    return filterpairs.size();

}


Int_t HParticlePairMaker::getSameRich(HParticleCand* cand,vector<HParticleCand*>& candidates,UInt_t flag,Bool_t isReference)
{
    // fills list of candidates which share the same RICH. If the RICH is not found
    // at all -1 is returned, other wise the number of candidates found for the hit.
    // The input candidate will be not included in candidates. Candidates can be filtered
    // by flag ( 0 == no filter, see flags eClosePairSelect+ ePairCase (hparticledef.h)).
    // The input candidate will be treated as refererence if isReference == kTRUE.

    candidates.clear();
    if(cand->getRichInd() !=-1 ){
        map<Int_t,vector<HParticleCand*> >::iterator it = mRichtoCand.find(cand->getRichInd());
        if(it != mRichtoCand.end()){
            vector<HParticleCand*>& v = it->second;
            for(UInt_t i = 0 ; i < v.size(); i++){

                if(v[i] != cand) {  // skip same cand
                    if(flag == 0 ) candidates.push_back(v[i]);
                    else {
                        Bool_t fl = kFALSE;
                        if(!isReference) fl=HParticleTool::evalPairsFlags(fl,cand,v[i]);
                        else             fl=HParticleTool::evalPairsFlags(fl,v[i],cand);
                        if(fl){
                            candidates.push_back(v[i]);
                        }
                    }
                }
            } //end loop
        } else return -1;
    }
    return candidates.size();
}

Int_t HParticlePairMaker::getSameInnerMdc(HParticleCand* cand,vector<HParticleCand*>& candidates,UInt_t flag,Bool_t isReference)
{
    // fills list of candidates which share the same inner MDC seg. If the seg is not found
    // at all -1 is returned, other wise the number of candidates found for the hit.
    // The input candidate will be not included in candidates. Candidates can be filtered
    // by flag ( 0 == no filter, see flags eClosePairSelect+ ePairCase (hparticledef.h)).
    // The input candidate will be treated as refererence if isReference == kTRUE.

    candidates.clear();
    if(cand->getInnerSegInd() !=-1 ){
        map<Int_t,vector<HParticleCand*> >::iterator it = mInnerMdctoCand.find(cand->getInnerSegInd());
        if(it != mInnerMdctoCand.end()){
            vector<HParticleCand*>& v = it->second;
            for(UInt_t i = 0 ; i < v.size(); i++){

                if(v[i] != cand) {  // skip same cand
                    if(flag == 0 ) candidates.push_back(v[i]);
                    else {
                        Bool_t fl = kFALSE;
                        if(!isReference) fl=HParticleTool::evalPairsFlags(fl,cand,v[i]);
                        else             fl=HParticleTool::evalPairsFlags(fl,v[i],cand);
                        if(fl){
                            candidates.push_back(v[i]);
                        }
                    }
                }
            } //end loop
        } else return -1;
    }
    return candidates.size();
}

Int_t HParticlePairMaker::getSameOuterMdc(HParticleCand* cand,vector<HParticleCand*>& candidates,UInt_t flag,Bool_t isReference)
{
    // fills list of candidates which share the same outer MDC seg. If the seg is not found
    // at all -1 is returned, other wise the number of candidates found for the hit.
    // The input candidate will be not included in candidates. Candidates can be filtered
    // by flag ( 0 == no filter, see flags eClosePairSelect+ ePairCase (hparticledef.h)).
    // The input candidate will be treated as refererence if isReference == kTRUE.

    candidates.clear();
    if(cand->getOuterSegInd() !=-1 ){
        map<Int_t,vector<HParticleCand*> >::iterator it = mOuterMdctoCand.find(cand->getOuterSegInd());
        if(it != mOuterMdctoCand.end()){
            vector<HParticleCand*>& v = it->second;
            for(UInt_t i = 0 ; i < v.size(); i++){

                if(v[i] != cand) {  // skip same cand
                    if(flag == 0 ) candidates.push_back(v[i]);
                    else {
                        Bool_t fl = kFALSE;
                        if(!isReference) fl=HParticleTool::evalPairsFlags(fl,cand,v[i]);
                        else             fl=HParticleTool::evalPairsFlags(fl,v[i],cand);
                        if(fl){
                            candidates.push_back(v[i]);
                        }
                    }
                }
            } //end loop
        } else return -1;
    }
    return candidates.size();
}

Int_t HParticlePairMaker::getSameMeta(HParticleCand* cand,vector<HParticleCand*>& candidates,UInt_t flag,Bool_t isReference)
{
    // fills list of candidates which share the same META hit. The number of candidates
    // found for the hit is returned or -1 if the hit is not found at all..
    // The input candidate will be not included in candidates. Candidates can be filtered
    // by flag ( 0 == no filter, see flags eClosePairSelect+ ePairCase (hparticledef.h)).
    // The input candidate will be treated as refererence if isReference == kTRUE.

    candidates.clear();

    Int_t metaind = cand->getMetaHitInd();
    Int_t sel     = cand->getSelectedMeta();

    if(metaind != -1 ){
        map<Int_t,vector<HParticleCand*> >* mp =0;

        if     (sel == kTofClst)           mp = &mTofClsttoCand;
        else if(sel == kTofHit1||sel==kTofHit2) mp = &mTofHittoCand;
        else if(sel == kRpcClst)           mp = &mRpcClsttoCand;
        else if(sel == kShowerHit)         mp = &mShowertoCand;
        else if(sel == kEmcClst)           mp = &mEmctoCand;
        else  {
            Error("getSameMeta()","unknown Meta case!"); return 0;
        }
        map<Int_t,vector<HParticleCand*> >::iterator it = mp->find(metaind);

        if(it != mp->end()){
            vector<HParticleCand*>& v = it->second;
            for(UInt_t i = 0 ; i < v.size(); i++){

                if(v[i] != cand) {  // skip same cand
                    if(flag == 0 ) candidates.push_back(v[i]);
                    else {
                        Bool_t fl = kFALSE;
                        if(!isReference) fl=HParticleTool::evalPairsFlags(fl,cand,v[i]);
                        else             fl=HParticleTool::evalPairsFlags(fl,v[i],cand);
                        if(fl){
                            candidates.push_back(v[i]);
                        }
                    }
                }
            } //end loop
        } else return -1;
    } 
    return candidates.size();
}

void  HParticlePairMaker::plotPairCaseStat()
{
    // draws the fraction of pair cases for
    // the default pairCases1...12.

    TCanvas* c = new TCanvas("pairCase","pairCase",1500,800);
    c->SetGridx();
    c->SetGridy();
    c->SetBottomMargin(0.15);

    UInt_t bins = fCaseVec.size();

    TH1F* h = new TH1F("h","",bins,0,bins);
    h->SetLineColor(kRed);
    h->SetLineWidth(3);
    h->SetYTitle("Fraction of pairs [%]");
    h->GetXaxis()->SetNdivisions(100+bins);
    h->GetXaxis()->SetLabelSize(0);

    h->GetYaxis()->SetRangeUser(0.,50.);
    Float_t sum = 0;
    for(UInt_t i=0;i<bins;i++)  sum+=fCaseCt[i];

    cout<<"Pairs Cases for Leptons :"<<endl;
    for(UInt_t i=0;i<bins;i++){
        Float_t frac = 0;
        if(sum != 0) frac = (fCaseCt[i]/sum)*100.;
        if(sum !=0 ) h->Fill(i,frac);
        cout<<Form("case%2i : ",i+1)<<setw(15)<< frac << " cts "<<setw(15)<<fCaseCt[i]<<endl;

    }
    h->Draw();

    Double_t offsetX  = 0.25;
    Double_t offsetY  =   -8;
    Double_t scaleX   = 1.;
    Double_t scaleY   = 1.;

    HParticlePairDraw pairdraw;

    for(UInt_t i=0; i < bins; i++){
        pairdraw.drawPair(offsetX+i,offsetY,scaleX,scaleY,fCaseVec[i],Form("Case%i",i+1),"");
    }
    cout<<"sum "<<sum<< " check " <<richCandCt<<endl;
}