MakeGroupList.h

Go to the documentation of this file.

#ifndef MAKEGROUPLIST__H
#define MAKEGROUPLIST__H
 
#include "CLHEP/Vector/LorentzVector.h"
#include "Criteria.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
#include "Pi0Cut.h"
#include "VertexFit/KinematicFit.h"
#include <algorithm>
#include <list>
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);
  }*/
} // namespace Pi0
#endif