// $Id:$
// ============================================================================
// CVS tag $Name:  $ , version $Revision:$
// ============================================================================
// $Log:$
// ============================================================================
// Include files 
// ============================================================================
// LoKi
// ============================================================================
#include "LoKi/LoKi.h"
// ============================================================================

/** @file
 *  Simple algorithm for Bs -> mu+ mu- analysis 
 *
 *  @author Vanya  BELYAEV   Ivan.Belyaev@lapp.in2p3.fr 
 *  @author Andrey GOLUTVIN  Andrey.Goloutvin@cern.ch 
 *
 *  @date 2005-07-21 
 */

LOKI_ALGORITHM( LoKiMuMu ) 
{
  using namespace LoKi       ;
  using namespace LoKi::Cuts ;
  using namespace LoKi::Fits ;
  
  const StatusCode OK = StatusCode::SUCCESS ;

  const EventHeader* header = 
    get<EventHeader> ( EventHeaderLocation::Default ) ;
  if ( 0 == header  ) { return StatusCode::FAILURE ; }


  VRange primaries = vselect ( "PVs" , Vertex::Primary == VTYPE ) ;
  if ( primaries.empty() ) 
  { return Warning ( "No primary vertices are found!", OK ) ; }
  Fun mips = MIPCHI2 ( primaries , geo() ) ;
  
  // select muons 
  Range muons   = select 
    ( "muons"           , 
      "mu+" == ABSID   &&    // muons 
      PT   > 1.0 * GeV &&    // large pt 
      PPBITMUON        &&    // has muon information 
      mips > 4.0        ) ;  // do not point to the primary 
  
  Range muplus  = select ( "mu+"   , muons , Q > 0  ) ;
  Range muminus = select ( "mu-"   , muons , Q < 0  ) ;
  
  if ( muplus.empty() || muminus.empty() ) { return StatusCode::SUCCESS  ; }
  
  // a bit of MC 
  MCMatch mc     = mctruth() ;
  
  // true decay 
  MCRange _sig   = mc->findDecays( " [  B_s0         -> mu+ mu- {, gamma } { , gamma } ]cc" ) ;
  
  // hadronic decays 
  MCRange _pipi  = mc->findDecays( " [ { B_s0 , B0 } -> pi+ pi- {, gamma } { , gamma } ]cc ") ;
  MCRange _k1pi  = mc->findDecays( " [ { B_s0 , B0 } -> K+  pi- {, gamma } { , gamma } ]cc ") ;
  MCRange _k2pi  = mc->findDecays( " [ { B_s0 , B0 } -> K-  pi+ {, gamma } { , gamma } ]cc ") ;
  MCRange _kk    = mc->findDecays( " [ { B_s0 , B0 } -> K+  K-  {, gamma } { , gamma } ]cc ") ;
  
  // muons from pi,K 
  MCRange _muH   = mc->findDecays( " [ {K- , pi-} -> ^mu- ... ]cc " ) ;
  
  // all muons 
  MCCut   mccut  = abs( MCVXFUN( MCVZ ) ) < (  60 *  cm )  ;
  MCRange _mu    = mcselect ( "mcmu" , "mu+" == MCABSID && mccut ) ;
  
  MCTRUTH _hhMC   = MCTRUTH ( mc , _pipi ) ;
  _hhMC.add( _k1pi ) ;
  _hhMC.add( _k2pi ) ;
  _hhMC.add( _kk   ) ;
  
  Cut mmMC = MCTRUTH ( mc , _sig ) ;
  Cut hhMC = _hhMC                 ;
  Cut muMC = MCTRUTH ( mc , _mu  ) ;
  Cut mhMC = MCTRUTH ( mc , _muH ) ;
  

  Tuple tuple   = nTuple ("MuMu N-Tuple") ;
  
  for ( Loop Bs = loop( "mu+ mu-" , "B_s0" ) ; Bs ; ++Bs ) 
  {
    const double mass = Bs -> mass(1,2) ;
    if ( mass < 2.8 * GeV ) { continue ; }
    
    if ( VCHI2 ( Bs ) < 0 || VCHI2 ( Bs ) > 16 ) { continue ; } // CONTINUE 
    
    const Particle* B   = Bs ;
    const Particle* mu1 = Bs(1) ;
    const Particle* mu2 = Bs(2) ;
    
    if ( 0 == B || 0 == mu1 || 0 == mu2 ) { continue ; } // CONTINUE 
    
    // check the mass 
    //if ( M( B ) < 4 * GeV               ) { continue ; } // CONTINUE 
    
    // ask Bs to point to the primary vertex 
    if ( mips( Bs ) > 25                ) { continue ; } // CONTINUE

    // find Primary Vertex:
    // # from upstream PVxs select PVx with
    //   the minimal B-impact parameter significance 
    VFun vipChi2 =  VIPCHI2 ( Bs , geo  () ) ;
    VCut vCut    = ( vipChi2 >= 0 ) && ( VPSD( Bs , geo  () )  > -0.5 * mm ) ;
    VRange::iterator iprim  = 
      LoKi::select_min ( primaries.begin ()   , // begin of vertices
                         primaries.end   ()   , // end of vertices
                         vipChi2              , // to be minimized 
                         vCut                 ) ; // constraints  
    if ( primaries.end() == iprim ) 
    { Warning ( "'Good' primary vertex is not found " )  ; continue ; }
    const Vertex* prim = *iprim ;
    if ( 0 == prim ) 
    { Warning ( "'Good' primary vertex points to NULL " ) ; continue ; }
    
    tuple << Tuples::Column( "" , header ) ;
    
    tuple -> column ( "nPV"   , primaries.size() ) ;
    
    tuple -> column ( "M"     , M  ( Bs  )               / GeV ) ;
    
    tuple -> column ( "PT1"   , PT ( mu1 )               / GeV ) ;
    tuple -> column ( "PT2"   , PT ( mu2 )               / GeV ) ;
    
    tuple -> column ( "MU1"   , mu1->momentum()          / GeV ) ;
    tuple -> column ( "MU2"   , mu2->momentum()          / GeV ) ;
    
    tuple -> column ( "VCHI2" , VCHI2  ( Bs  )                 ) ;
    
    tuple -> column ( "BV"    , Bs->vertex()->position() / mm  ) ;
    tuple -> column ( "PV"    , prim->position()         / mm  ) ;
    
    Fun dira = DIRA    ( prim   ) ;
    // direction angle 
    tuple -> column ( "DIRA"  , dira ( Bs ) ) ;
    tuple -> column ( "LV01"  , LV01 ( Bs ) ) ;
    
    tuple -> column ( "mipB" , mips( Bs  ) ) ;
    
    tuple -> column ( "mip1" , mips( mu1 ) ) ;
    tuple -> column ( "mip2" , mips( mu2 ) ) ;
    
    Fun vdot = VDDOT ( prim ) ;
    // distance between PV and Bs vertex  along B-momentum
    tuple -> column ( "cdot" , vdot    ( Bs ) ) ;
    
    Fun vddtime = VDDTIME ( prim ) ;
    // time distance between PV and Bs vertex 
    // along B-momentum in c*tau units 
    tuple -> column ( "vddt" , vddtime ( Bs ) ) ;
    
    Fun mvd    = MVDCHI2( primaries ) ;
    // minimal chi2 of the distance between Bs vertex 
    // and *ALL* primary vertices 
    tuple -> column ( "mvd" , mvd ( Bs ) ) ;
    
    // verted distance in chi2 
    VFun vchi2d = VCHI2D( prim ) ;
    tuple -> column ( "vdchi2" , vchi2d( Bs ) ) ;
    
    // DLL(mu-pi) for mu1 
    tuple -> column ( "dll1"  , PIDmu     ( mu1 ) ) ;
    // DLL(mu-pi) for mu2 
    tuple -> column ( "dll2"  , PIDmu     ( mu2 ) ) ;
    // has MUON information ?
    tuple -> column ( "mbit1" , PPBITMUON ( mu1 ) ) ;
    // has MUON information ?
    tuple -> column ( "mbit2" , PPBITMUON ( mu2 ) ) ;
    
    // Bs -> mu+ mu- ? `
    tuple -> column ( "mc1bs"  , mmMC ( Bs  ) ) ;
    tuple -> column ( "mc1mu1" , mmMC ( mu1 ) ) ;
    tuple -> column ( "mc1mu2" , mmMC ( mu2 ) ) ;
    // B -> HH ? 
    tuple -> column ( "mc2bs"  , hhMC ( Bs  ) ) ;
    tuple -> column ( "mc2mu1" , hhMC ( mu1 ) ) ;
    tuple -> column ( "mc2mu2" , hhMC ( mu2 ) ) ;
    // true muons ? 
    tuple -> column ( "mcmu1"  , muMC ( mu1 ) ) ;
    tuple -> column ( "mcmu2"  , muMC ( mu2 ) ) ;
    // decay-in-flight ?
    tuple -> column ( "mchmu1" , mhMC ( mu1 ) ) ;
    tuple -> column ( "mchmu2" , mhMC ( mu2 ) ) ;

    Bs->save("Bs") ;
    
    tuple->write() ;
    
  } ;
  
  if ( !selected("Bs").empty() ) { setFilterPassed ( true ) ; }
  
  return StatusCode::SUCCESS ;
};


// ============================================================================
// The END 
// ============================================================================