/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/Pi0RecAlg/Pi0RecAlg-00-00-06/Pi0RecAlg/MakeGroupList.h

Go to the documentation of this file.

#ifndef MAKEGROUPLIST__H
#define MAKEGROUPLIST__H

#include <list>
#include <algorithm>
#include "Criteria.h"
#include "VertexFit/KinematicFit.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
#include "Pi0Cut.h"
namespace Pi0{


        typedef std::list<Criteria> CriteriaList;
        typedef std::list<Criteria>::iterator CriteriaIterator;
        typedef std::list<Criteria>::const_iterator Const_CriteriaIteator;

        typedef std::list<GammaPair> Pi0List;
        typedef std::list<GammaPair>::iterator Pi0Iterator;
        typedef std::list<GammaPair>::const_iterator Const_Pi0Iterator;

        typedef std::list<EvtRecTrack*> GammaList;
        typedef std::list<EvtRecTrack*>::iterator GammaIterator;
        typedef std::list<EvtRecTrack*>::const_iterator Const_GammaIterator;

        CriteriaList default_criteria_list;
        Pi0List candidate_pi0_list;
        Pi0List filter_pi0_list;
        GammaList default_gamma_list;//one copy of all photons
//      GammaList update_gamma_list;//a subset of default_gamma_list 

        GammaList& GetDefaultGammaList() { return default_gamma_list;}
//        GammaList& GetUpdateGammaList() { return update_gamma_list;}
        Pi0List& GetCandidatePi0List() { return candidate_pi0_list;}
        Pi0List& GetFilterPi0List() { return filter_pi0_list;}

//      GammaList& make_gamma_list(EvtRecEvent* recEvt, EvtRecTrackCol* recTrkCol)
        GammaList& make_gamma_list(UserPi0Cut& cut)
        {
                EvtRecEvent* recEvt = UserPi0Cut::recEvt;
                EvtRecTrackCol* recTrkCol = UserPi0Cut::recTrkCol;

                default_gamma_list.clear();
//              update_gamma_list.clear();
                for(int i1 = recEvt->totalCharged(); i1 < (recEvt->totalTracks()); ++i1) 
                {

                        EvtRecTrack* gTrk = *(recTrkCol->begin() + i1);
                        if(cut.isGoodPhoton(gTrk))
                        {
                                default_gamma_list.push_back(gTrk);
//                              update_gamma_list.push_back(gTrk);
                        }
                }
                return default_gamma_list;
        }
/*      GammaList& remake_gamma_list_by_remove(const GammaPair& gp)//remove two gamma ref by gp from list
        {
                update_gamma_list.remove(gp.First);
                update_gamma_list.remove(gp.Second);
                return update_gamma_list;
        }*/
        void print_gamma_list(const GammaList& list)
        {
                if(list.size() == 0) {std::cout<<"GammaList->()"<<std::endl; return;}
                std::cout<<"GammaList->(";
                for(Const_GammaIterator it = list.begin(); it!= list.end(); it++)
                {
                        std::cout<<(*it)->trackId()<<", "; 
                }
                std::cout<<')'<<std::endl;
        }
        //=============
        Pi0List& make_pi0_list(const GammaList& gamma_list)
        {
                candidate_pi0_list.clear();
                Const_GammaIterator i_it = gamma_list.begin();
                Const_GammaIterator i_it_end = gamma_list.end();
                Const_GammaIterator j_it_end = gamma_list.end();

                --i_it_end;

                //              KinematicFit * kmfit = KinematicFit::instance();
                Const_GammaIterator j_it ;
                for(; i_it != i_it_end; ++i_it)
                {
                        for( j_it = i_it, ++j_it ; j_it != j_it_end; ++j_it)
                        {

                                EvtRecTrack* g1Trk = *i_it;
                                EvtRecTrack* g2Trk = *j_it;

                                RecEmcShower* g1Shower = g1Trk->emcShower();
                                RecEmcShower* g2Shower = g2Trk->emcShower();
                                double inv_m = (getP4(g1Shower) + getP4(g2Shower)).m();

                                /*                              kmfit->init();
                                                                kmfit->AddTrack(0, 0.0, g1Shower);
                                                                kmfit->AddTrack(1, 0.0, g2Shower);
                                                                kmfit->AddResonance(0, 0.1349766, 0, 1);
                                                                kmfit->Fit(0);*/

                                candidate_pi0_list.push_back(GammaPair(*i_it, *j_it, inv_m)); 
                        }
                }
                return candidate_pi0_list;
        }
        void print_pi0_list(const Pi0List& list)
        {
                std::cout<<"OK Pi0List->{";
                for(Const_Pi0Iterator it = list.begin(); it!= list.end(); it++)
                {
                        std::cout<<"("<<(*it).inv_m<<","<<(*it).First->trackId()<<","<<(*it).Second->trackId()<<"), ";       
                }
                std::cout<<'}'<<std::endl;
        }
        Pi0List& apply_criteria(const Criteria& cri)
        {
                Pi0Iterator it = candidate_pi0_list.begin();
                for(; it!=candidate_pi0_list.end();)
                {
                        if( cri.check(*it))
                                ++it;
                        else{
                                filter_pi0_list.push_back(*it);
                                it = candidate_pi0_list.erase(it);
                        }

                }
                return candidate_pi0_list;
        }
        void Pi0ListToTDS(const Pi0List& pi0list, EvtRecPi0Col* recPi0Col)
        {
                assert(recPi0Col);
                static double xmpi0= 0.1349766;
                KinematicFit * kmfit = KinematicFit::instance();

                for(Const_Pi0Iterator it = pi0list.begin(); it != pi0list.end(); ++it)
                {
                        EvtRecTrack* g1Trk = (*it).First;
                        EvtRecTrack* g2Trk = (*it).Second;
                        RecEmcShower* g1Shower = g1Trk->emcShower();
                        RecEmcShower* g2Shower = g2Trk->emcShower();

                        kmfit->init();
                        kmfit->AddTrack(0, 0.0, g1Shower);
                        kmfit->AddTrack(1, 0.0, g2Shower);
                        kmfit->AddResonance(0, xmpi0, 0, 1);
                        kmfit->Fit(0);

                        HepLorentzVector g1P4 = getP4(g1Shower);
                        HepLorentzVector g2P4 = getP4(g2Shower);
                        HepLorentzVector p2g = g1P4 + g2P4;
                        EvtRecPi0* recPi0 = new EvtRecPi0();

                        recPi0->setUnconMass(p2g.restMass());
                        recPi0->setChisq(kmfit->chisq(0));
                        if ( g1P4.e() >= g2P4.e() ) {
                                recPi0->setHiPfit(kmfit->pfit(0));
                                recPi0->setLoPfit(kmfit->pfit(1));
                                recPi0->setHiEnGamma(g1Trk);
                                recPi0->setLoEnGamma(g2Trk);
                        }
                        else {
                                recPi0->setHiPfit(kmfit->pfit(1));
                                recPi0->setLoPfit(kmfit->pfit(0));
                                recPi0->setHiEnGamma(g2Trk);
                                recPi0->setLoEnGamma(g1Trk);
                        }

                        recPi0Col->push_back(recPi0);
                }
        }
        //pi0 reconstruction strategy
        /*  void priority_method(int nPi0, EvtRecPi0Col* pi0_col) //reconstruct pi0 one by one
            {
            assert(pi0_col);
            int i = 0;
            Pi0List& _this = make_pi0_list(update_gamma_list);
            print_pi0_list(_this);
            Pi0List ret = the_most_optimized(_this);
            Pi0List2TDS(ret, pi0_col);
            ++i;
            while(i<nPi0){
            remake_gamma_list_by_remove(*(ret.begin()));
            Pi0List& _this = make_pi0_list(update_gamma_list);
            print_pi0_list(_this);
            Pi0List ret = the_most_optimized(_this);
            Pi0List2TDS(ret, pi0_col);
            };
            }*/
        /*      void general_method(EvtRecPi0Col* pi0_col)
                {
                make_pi0_list(default_gamma_list);
                Criteria cri_inv;
                apply_criteria(cri_inv);
                BasicCriteria cri_kfit;
                apply_criteria(cri_kfit);
                Pi0ListToTDS(candidate_pi0_list, pi0_col);
                }*/
        void apply_criterias()
        {
                for(CriteriaIterator it = default_criteria_list.begin(); it != default_criteria_list.end(); ++it)
                {
                        apply_criteria(*it);
                }
        }
        void add_cut(const Criteria& cri)
        {
                default_criteria_list.push_back(cri);
        }
        /*      void rec_pi0(EvtRecEvent* recEvt, EvtRecTrackCol* recTrkCol, EvtRecPi0Col* pi0_col, int nPi0=0)//nPi0, you want to reconstruct n Pi0 
                {
        //              assert(recEvt);assert(recTrkCol);
        GammaList& gamma_list = make_gamma_list(recEvt, recTrkCol);
        //              std::cout<<"gamma in total = "<<gamma_list.size()<<std::endl;
        //              std::cout<<"maximum pi0 N = "<<gamma_list.size()/2<<std::endl;
        //              if(gamma_list.size()/2 <nPi0) nPi0 = gamma_list.size()/2;
        //              print_gamma_list(gamma_list);
        //                priority_method(nPi0, pi0_col);
        make_pi0_list(gamma_list);
        apply_criterias();
        //              Criteria cri_inv;
        //              apply_criteria(cri_inv);
        //              KFitCriteria cri_kfit;
        //              apply_criteria(cri_kfit);
        candidate_pi0_list.sort(high_momentum());
        print_pi0_list(candidate_pi0_list);
        Pi0ListToTDS(candidate_pi0_list, pi0_col);
        }*/
}
#endif

---