/*****************************************************************************
 * Project: BaBar detector at the SLAC PEP-II B-factory
 * Package: BetaTools
 *    File: $Id: BtaSemilepCand.hh,v 1.27 2004/12/21 22:39:54 chcheng Exp $
 * Authors:
 *   DK, David Kirkby, Stanford University, kirkby@hep.stanford.edu
 * History:
 *   24-Aug-2000 DK Created initial version
 *   08-Sep-2000 DK Improve member function names, add deltaM()
 *
 * Copyright (C) 2000 Stanford University
 *****************************************************************************/
#ifndef BTA_SEMILEP_CAND
#define BTA_SEMILEP_CAND

#include "AbsEvent/AbsEvent.hh"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Geometry/HepPoint.h"
#include "difAlgebra/DifComplex.hh"
#include "difAlgebra/DifNumber.hh"
#include "BetaCoreTools/BtaMcAssoc.hh"
#include <string>
#include <map>

class BtaCandidate;

// This is a utility class for calculating observables associated with
// composite candidates for the semileptonic B decay chain:
//
//  B0 / B+ -> D* l (nu X) , D* -> D pi/gamma , D -> K X
//  B+ -> D0 l (nu X)      , D0 -> K- pi+

class BtaSemilepCand {
public:
  BtaSemilepCand(const BtaCandidate& mother, const HepLorentzVector& p4CMS);
  virtual ~BtaSemilepCand();
  inline bool isAllowedDecayTree() const {
    return _bMode != UnusedBMode &&
      (_dstarMode != UnusedDstarMode && _dMode != UnusedDMode) ||
      (_dMode != UnusedDMode);

  }
  // fill event cashe (user can specify the name of list used for PID map)
  bool fillEventCache(AbsEvent *theEvent, const std::string& key= "ChargedTracks");
  // return this candidate's B decay mode
  enum BDecayMode { UnusedBMode, B0ToDstarlnu, BchToDstarlnu, BchToDlnu, 
                                                              B0ToDlnu };
  inline BDecayMode getBDecayMode() const { return _bMode; }
  // return this candidate's D* decay mode
  enum DstarDecayMode { UnusedDstarMode, DstarToD0pi,
			DstarToDchpi0, DstarToDchgamma,
			Dstar0ToD0pi0, Dstar0ToD0gamma};
  inline DstarDecayMode getDstarDecayMode() const { return _dstarMode; }
  // return this candidate's D decay mode
  enum DDecayMode { UnusedDMode, D0ToKpi, D0ToK3pi, D0ToKpipi0, D0ToKspipi, DchToKpipi, DchToKspipi0, DchToKspi };
  inline DDecayMode getDDecayMode() const { return _dMode; }
  // get a pointer to our D* candidate (could be D or D*)
  inline const BtaCandidate *theHadron() const { return _hadron; }
  // get a pointer to our D* candidate (could be a D*0, D*+, or D*-)
  inline const BtaCandidate *theDstar() const { return _dstar; }
  // get a pointer to our lepton candidate (it might think it is a pion)
  inline const BtaCandidate *theLepton() const { return _lepton; }
  // get a pointer to our D*'s D daughter (could be a D0, D+, or D-)
  inline const BtaCandidate *theD() const { return _theD; }
  // get a pointer to our D*'s soft (non-D) daughter (could be a pi0,
  // pi+, pi-, or gamma)
  inline const BtaCandidate *theSoftDaughter() const { return _theSoft; }
  // get a pointer to our D's Kaon daughter (could be a Ks, K+, or K-)
  inline const BtaCandidate *theK() const { return _theK; }
  // get a pointer to our D's Pion daughter
  inline const BtaCandidate *thePi1() const { return _thePi1; }
  inline const BtaCandidate *thePi2() const { return _thePi2; }
  inline const BtaCandidate *thePi3() const { return _thePi3; }  
  // check track quality
  bool isGoodTrack(const BtaCandidate *aTrack, const double maxDoca=1.0, const double maxAbsZ0=3.0);
  // count bad tracks (including lepton, soft pi and D daughters)
  int nBadTracks();
  // get Kaon ID (user can specify the name of PID bit map )
  int getKaonId(const std::string& key="KSelectorsMap");
  // get lepton ID
  int getElecId(const std::string& key="eSelectorsMap");
  int getMuonId(const std::string& key="muSelectorsMap");
  // does our D decay into a mode with at least one pi0?
  bool  theDcontainPi0() const;
  // return the lepton candidate's momentum in the event CMS frame
  double leptonPstar() const;
  // return the charm hadron candidate's momentum in the event CMS frame
  double hadronPstar() const;
  // return the D* candidate's momentum in the event CMS frame
  double dstarPstar() const;
  // return the difference the masses of our D* and D candidates
  double deltaM() const;
  // return the cosine of the angle between our D* and lepton candidates
  // in the event CMS frame
  double cosThetaStar() const;
  // return the cosine of the angle between our B candidate and the vector
  // sum of our charm hadron and lepton candidates in the event CMS frame
  double cosThetaBY() const;
  // return the same above with flipped lepton
  double cosThetaBYflip() const;
  // return the Dalitz weight for our D decay (if applicable). 
  //The maximum is 1.
  double dDalitzWgt() const;
  //   Now make functions for the four kin vars needed for the FF (form factor) fit:
  //  (note that "B frame" thus far only means the event CMS frame, actually)
  // Return the cosine of the angle between the D* candidate in the B frame and its
  // D daughter in the D* frame:
  double cosThetaV() const;
  // Return the cosine of the angle between the virtual W in the B frame  and the lepton
  //  daughter in the W frame:
  double cosThetaL() const;
  // return the angle between the D*-D-pi_s plane and the W-l-nu plane
  // in the event CMS frame:
  double angleChi() const;
  // Return the invariant mass of the virtual W, which is equivalent to the 4-momentum transfer squared between the B and the D*
  // momentumTransfer2 = (pl+pnu)^2 = (pB-pDstar)^2
  // in the B rest frame.
  double momentumTransfer2() const;
  bool  boostOkFlag;

  bool fillMcCache(BtaMcAssoc* theMcMap);
  bool resolveMcAssoc();
  bool realD();
  bool realHadron();
  bool realDstar();
  bool dFromB();
  bool dlnuExclusive();
  bool hadronFromB();
  bool dstarFromB();
  bool leptFromB();
  bool leptFromSameB();
  bool directLept();
  int mcLeptId();
  int mcLeptBId();

protected:
  AbsEvent *_event;
  HepAList<BtaCandidate>* _pidTrkList;
  std::map<BtaCandidate*, int>* _kaonBitMap;
  std::map<BtaCandidate*, int>* _elecBitMap;
  std::map<BtaCandidate*, int>* _muonBitMap;
  std::string _kaonMapKey, _elecMapKey, _muonMapKey;
  BtaMcAssoc* _truthMap;
  bool findMother(BtaCandidate *dau, BtaCandidate *moth);
  BtaCandidate* findMother(BtaCandidate *dau, const char *moth);
  void testMcLept();
  BtaCandidate* _theMcD, *_theMcDstar, *_theMcSoft, *_theMcLept, *_theMcDstB, 
    *_theMcLB, *_theMcDB, *_theMcHadron, *_theMcHadB;
  
  bool _mcResolved, _realD, _realDstar, _realHadron, _dstarFromB, _leptFromB, _leptFromSameB, _directLept, _dFromB, _hadronFromB, _isDstar, _dlnuExclusive;
  int  _mcLeptId, _mcLeptBId;

  void loadBeamSpot();
  DifComplex brWigAmp(double msq12, double msq23, double mP,
		      double m1, double m2, double m3,
		      double ampl, double theta, double gamma,
		      double bwm, int spin) const;
  double decayAngle(double msq12, double msq23, double mP,
		    double m1, double m2, double m3) const;
  int nFinalDaus(const BtaCandidate& c, std::string terminals, bool print);
  void fillnFinal( const BtaCandidate& c, std::string terminals, int* nFinal, bool print);
  const BtaCandidate *_mother, *_dstar, *_lepton, *_theD, *_theSoft, *_theK;
  const BtaCandidate *_thePi1, *_thePi2, *_thePi3, *_hadron;
  BtaCandidate *_leptonTrk, *_kaonTrk;
  HepLorentzVector _dstarP4CMS, _leptonP4CMS, _dP4CMS, _hadronP4CMS;
  HepLorentzVector _dP4DstarFrame, _leptonP4WFrame, _pW;
  HepPoint  _beamSpot;
  double _energyBCMS;
  BDecayMode _bMode;
  DstarDecayMode _dstarMode;
  DDecayMode _dMode;
};

#endif