/*****************************************************************************
 * Project: BaBar detector at the SLAC PEP-II B-factory
 * Package: BetaTools
 *    File: $Id: BtaSemilepCand.cc,v 1.45 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
 *   06-Aug-2000 DK Fix typo which caused D0->Kspipi decays to be
 *                  considered invalid: "KS_0" -> "K_S0"
 *   21-Sep-2000 C.Cheng Implement Dalitz weighting
 *   05-Oct-2000 M.S. Gill -- Add new kin variables for FF fit
 *
 * Copyright (C) 2000 Stanford University
 *****************************************************************************/
#include "BaBar/BaBar.hh"
#include "BetaTools/BtaSemilepCand.hh"

#include "AbsEnv/AbsEnv.hh"
#include "AbsEvent/getTmpAList.hh"
#include "Beta/EventInfo.hh"  
#include "Beta/BtaCandidate.hh"
#include "CLHEP/Alist/AIterator.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "PDT/PdtPdg.hh"
#include "PDT/Pdt.hh"
#include "ErrLogger/ErrLog.hh"
#include "TrkFitter/TrkPocaXY.hh"
#include "TrkBase/TrkFit.hh"
#include "TrkBase/TrkRecoTrk.hh"
#include "PepEnv/PepEnv.hh"
#include "PepEnvData/PepBeamSpotCal.hh"
//#include "ProxyDict/IfdStrKey.hh"
#include "BetaMicroTruth/BtaMcAssocMicro.hh"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/String/Strings.h"
#include <string>
#include <sstream>
#include <map>
using std::cout;
using std::endl;
using std::string;

BtaSemilepCand::BtaSemilepCand(const BtaCandidate& mother,
			       const HepLorentzVector& p4CMS) :
  _event(0), _pidTrkList(0), _kaonBitMap(0), _elecBitMap(0), _muonBitMap(0),
  _kaonMapKey(""),_elecMapKey(""), _muonMapKey(""),
  _truthMap(0), _isDstar(false),
  _mother(&mother), _dstar(0), _lepton(0),
  _theD(0), _theSoft(0), _theK(0), _thePi1(0), _thePi2(0), _thePi3(0),  
  _hadron(0),  _leptonTrk(0), _kaonTrk(0),
  _beamSpot(-999.,-999.,-999.), _energyBCMS(0),
  _bMode(UnusedBMode), _dstarMode(UnusedDstarMode), _dMode(UnusedDMode)
{
  // we must have exactly 2 daughters
  if(mother.nDaughters() != 2) return;
  // store pointers to these daughters, assuming the first is the hadron
  _energyBCMS= p4CMS.m()/2;
  HepAListIterator<BtaCandidate> next= mother.daughterIterator();
  _hadron= next();
  _lepton= next();

  // check decay mode
  // for now, we only deal with hadron= D*+-, D*0, D+- or D0
  const int lund(_hadron->pdtEntry()->lundId());
  switch(lund) {
  case PdtLund::D_star_plus:
  case PdtLund::D_star_minus:
  case PdtLund::D_star0:
  case PdtLund::anti_D_star0:
    _isDstar = true;
    _dstar = _hadron;
  case PdtLund::D_plus:
  case PdtLund::D_minus:
  case PdtLund::D0:
  case PdtLund::anti_D0:
    break;
  default:
    return;
  }

  if (_isDstar){
    if (_dstar->nDaughters() != 2) return;
    if (mother.charge()!=0 && _dstar->charge() == 0) _bMode= BchToDstarlnu;
    else if (mother.charge()==0 && _dstar->charge()!= 0) _bMode= B0ToDstarlnu;
    else
      return;
  }
  else {
    _theD = _hadron;
        if (mother.charge()!=0 && _hadron->charge()== 0) _bMode= BchToDlnu;
  else if ( mother.charge()==0 && _hadron->charge()!= 0) _bMode= B0ToDlnu;
  else
    return;
  }

  BtaCandidate *theDau;
  if(_isDstar){
    // D*'s decay products must contain one D and a Pi or gamma
    next= _dstar->daughterIterator();
    while ( (theDau= next()) ) {
      const PdtEntry* dauPdt= theDau->pdtEntry();
      if (dauPdt == Pdt::lookup("D0") ||
	  dauPdt == Pdt::lookup("anti-D0") ||
	  dauPdt == Pdt::lookup("D+") ||
	  dauPdt == Pdt::lookup("D-") )  _theD= theDau;
      else if ( dauPdt == Pdt::lookup("pi+") ||
		dauPdt == Pdt::lookup("pi-") ||
		dauPdt == Pdt::lookup("pi0") ||
		dauPdt == Pdt::lookup("gamma") ) _theSoft= theDau;
    }
    if ( !_theD || !_theSoft) return;
    if (_dstar->charge()== 0) {
      if(_theD->charge() != 0) return;
      if(_theSoft->pdtEntry() == Pdt::lookup("pi0"))
	_dstarMode= Dstar0ToD0pi0;
      else if(_theSoft->pdtEntry()== Pdt::lookup("gamma")) 
	_dstarMode= Dstar0ToD0gamma;
    }
    else { //D*+
      if ( _theD->charge()==0) { 
	if (_theSoft->pdtEntry()== Pdt::lookup("pi+") || 
	    _theSoft->pdtEntry()== Pdt::lookup("pi-") ) _dstarMode= DstarToD0pi;
      } 
      else { 
	if (_theSoft->pdtEntry()== Pdt::lookup("pi0")) 
	  _dstarMode= DstarToDchpi0;
	else if (_theSoft->pdtEntry()== Pdt::lookup("gamma")) 
	  _dstarMode= DstarToDchgamma;
      }  
    } 
  }

  // D's daughters
  int nPi(0), nPi0(0);
  next= _theD->daughterIterator();
  while ( (theDau= next()) ) {
    const PdtEntry* dauPdt= theDau->pdtEntry();
    if (dauPdt == Pdt::lookup("K+") ||
	dauPdt == Pdt::lookup("K-") ||
	dauPdt == Pdt::lookup("K_S0") ) _theK= theDau;
    else if (dauPdt == Pdt::lookup("pi+") ||
	     dauPdt == Pdt::lookup("pi-") ) {
      nPi++;
      if ( !_thePi1) _thePi1= theDau;
      else if (!_thePi2) _thePi2= theDau;
      else if (!_thePi3) _thePi3= theDau;
    }
    else if (dauPdt == Pdt::lookup("pi0") ) {
      nPi0++;
      if ( _theD->nDaughters()==2) {
	if (!_thePi1) _thePi1= theDau;
      } else if ( _theD->nDaughters()==3) {
	if (!_thePi2) _thePi2= theDau;
      } else if ( _theD->nDaughters()==4) {
	if (!_thePi3) _thePi3= theDau;
      }
    }
  }
  if ( !_theK) return;
  if ( _theD->charge()!= 0) {  // Dch
    if ( _theK->pdtEntry()== Pdt::lookup("K_S0")) {
      if ( _theD->nDaughters()== 3 && nPi==1 && nPi0==1) _dMode= DchToKspipi0;
      if ( _theD->nDaughters()== 2 && nPi==1 && nPi0==0) _dMode= DchToKspi;
    }
    else 
      if ( _theD->nDaughters()== 3 && nPi==2 ) _dMode= DchToKpipi;  
  }
  else { // D0
    if ( _theK->pdtEntry()== Pdt::lookup("K_S0")) {
      if ( _theD->nDaughters()== 3 && nPi== 2) _dMode= D0ToKspipi;
    }
    else
      if ( _theD->nDaughters()== 2 && nPi== 1) _dMode= D0ToKpi;
      else if ( _theD->nDaughters()== 4 && nPi== 3) _dMode= D0ToK3pi;
      else if ( _theD->nDaughters()== 3 && nPi== 1 && nPi0== 1) _dMode= D0ToKpipi0;
  }

  // set flag saying that we haven't resolved mc association
  _mcResolved= false;

  // Calculate the daughter 4-momenta in the CMS frame  
  _hadronP4CMS= _hadron->p4();
  _hadronP4CMS.boost(-p4CMS.boostVector());
  
  if(_isDstar){
    _dstarP4CMS= _hadronP4CMS;
  }

  _leptonP4CMS= _lepton->p4();
  _leptonP4CMS.boost(-p4CMS.boostVector());

  _dP4CMS= _theD->p4();
  _dP4CMS.boost(-p4CMS.boostVector());

  // And the D, lepton 4-momenta in the D*, W frames

  if(_isDstar){
    _dP4DstarFrame = _dP4CMS;
    _dP4DstarFrame.boost(-_hadronP4CMS.boostVector());
  }
  
  _leptonP4WFrame = _leptonP4CMS;
  
  double massB0= Pdt::lookup("B0")->mass();
  
  _pW.setVect(-Hep3Vector(_hadronP4CMS));
  _pW.setE(massB0- _hadronP4CMS.e());
  
  if (_pW.e() >= _pW.rho() ) 
    {
      _leptonP4WFrame.boost(-_pW.boostVector());  
      boostOkFlag=true;
    }
  else boostOkFlag=false;
}

BtaSemilepCand::~BtaSemilepCand() { }

bool BtaSemilepCand::fillEventCache(AbsEvent* theEvent, const string& key) {
  // check if we have already seen this event
  if ( theEvent== _event) return true;
  
  // do we have a valid event?
  if (0 == theEvent){
    ErrMsg(error) << "unable to get a pointer to this event" << endmsg;
    return false;
  }
  
  // lookup this event's info
  HepAList<EventInfo>* eventInfoList;
  getTmpAList(theEvent, eventInfoList, IfdStrKey("Default") );
  if(0 == eventInfoList) {
    ErrMsg(error) << "unable to get the info list for this event" << endmsg;
    return false;
  }
  EventInfo* eventInfo = eventInfoList->first();
  if(0 == eventInfo) {
    ErrMsg(error) << "event has an empty info list" << endmsg;
    return false;
  }
  
  // get list for PID
  getTmpAList( theEvent, _pidTrkList, IfdStrKey(key.c_str()));
  if ( 0== _pidTrkList) {
    ErrMsg(error) << "unable to get the " << key  
		  << " list for this event" << endmsg;
    return false;
  }

  // find the lepton and the Kaon in the list
  HepAListIterator<BtaCandidate> trkIter(*_pidTrkList);
  BtaCandidate *trk;
  BtaGenUid leptonUid= _lepton->uid();
  BtaGenUid kaonUid= _theK->uid();

  while ((trk= trkIter()) && !_leptonTrk)
    if (trk->uid() == leptonUid ) _leptonTrk= trk;

  if ( _theK->charge() !=0) {
    trkIter.rewind();  
    while ((trk= trkIter()) &&  !_kaonTrk ) 
      if (trk->uid() == kaonUid ) _kaonTrk= trk;
  }

  if (!_leptonTrk) {
    ErrMsg(error) << "unable to find the lepton in the list "
		  << key << endmsg;
    return false;
  }

  if ( _theK->charge() !=0  && !_kaonTrk) {
    ErrMsg(error) << "unable to find the kaon in the list "
		  << key << endmsg;
    return false;
  }

  // remember this event
  _event= theEvent;
  return true;
}

bool BtaSemilepCand::isGoodTrack(const BtaCandidate *aTrack, const double maxDoca, const double maxAbsZ0)  {
  if (0==aTrack) return false;
  if (aTrack->charge() == 0) return false;

  if (_beamSpot.x()<-99.) loadBeamSpot();
  const TrkFit* fit= aTrack->recoTrk()->fitResult();
  if (!fit) return false;
  TrkPocaXY pocaxy(fit->traj(),0,_beamSpot);
  Hep3Vector poca= fit->position(pocaxy.fltl1())- _beamSpot;
  if (sqrt(poca.x()*poca.x()+poca.y()*poca.y()) > maxDoca) return false;
  if (fabs(poca.z()) > maxAbsZ0) return false;
  return true;
}

int BtaSemilepCand::nBadTracks() {
  int nBad=0;
  if ( _lepton && !isGoodTrack(_lepton)) nBad++;
  if ( _theSoft && _theSoft->charge()!=0 && !isGoodTrack(_theSoft)) nBad++;
  if ( _theK && _theK->charge()!=0 && !isGoodTrack(_theK)) nBad++;
  if ( _thePi1 && _thePi1->charge()!=0 && !isGoodTrack(_thePi1)) nBad++;
  if ( _thePi2 && _thePi2->charge()!=0 && !isGoodTrack(_thePi2)) nBad++;
  if ( _thePi3 && _thePi3->charge()!=0 && !isGoodTrack(_thePi3)) nBad++;

  return nBad;
}

bool BtaSemilepCand::theDcontainPi0() const {
  HepAListIterator<BtaCandidate> next= _theD->daughterIterator();
  BtaCandidate *thePi0(0), *dau;
  while( (dau= next()) ) {
    if(dau->pdtEntry() == Pdt::lookup("pi0") ) 
      thePi0= dau;
  }
  return (thePi0!=0);
}

double BtaSemilepCand::leptonPstar() const {
  return _leptonP4CMS.rho();
}

double BtaSemilepCand::hadronPstar() const {
  return _hadronP4CMS.rho();
}

double BtaSemilepCand::dstarPstar() const {
  return _dstarP4CMS.rho();
}

double BtaSemilepCand::cosThetaStar() const {
  return cos(_hadronP4CMS.angle(_leptonP4CMS.vect()));
}

double BtaSemilepCand::cosThetaBY() const {
  double massB0= Pdt::lookup("B0")->mass();
  // "Effective" B energy used as Reference Point
  // to rescale cos Theta BY for off-peak data. 
  // 10.578210106 GeV is the mean Sqrt(s)  for Run1-4
  static const double energyBCMS_OffPeak_RP = 10.578210106 / 2;

  // check if it's off-peak data
  if (_energyBCMS < massB0) {
    massB0= massB0* _energyBCMS / energyBCMS_OffPeak_RP;
  }
  double pB0= sqrt(_energyBCMS*_energyBCMS - massB0*massB0);

  HepLorentzVector hadLepP4CMS= _hadronP4CMS + _leptonP4CMS;
  double hadLepMass= hadLepP4CMS.m();
  double hadLepE= hadLepP4CMS.e();
  double hadLepP= hadLepP4CMS.rho();

  return -(massB0*massB0 + hadLepMass*hadLepMass - 2*_energyBCMS*hadLepE)/
    (2*pB0*hadLepP);
}

double BtaSemilepCand::cosThetaBYflip() const {
  double massB0= Pdt::lookup("B0")->mass();
  double pdgUpsMass= Pdt::lookup("Upsilon(4S)")->mass();
  // check if it's off-peak data
  if (_energyBCMS < massB0) {
    massB0= massB0* _energyBCMS*2 / pdgUpsMass;
  }
  double pB0= sqrt(_energyBCMS*_energyBCMS - massB0*massB0);

  HepLorentzVector leptonFlip(-_leptonP4CMS.x(),-_leptonP4CMS.y(),-_leptonP4CMS.z(),_leptonP4CMS.t());
  HepLorentzVector hadLepP4CMS= _dstarP4CMS + leptonFlip;
  double hadLepMass= hadLepP4CMS.m();
  double hadLepE= hadLepP4CMS.e();
  double hadLepP= hadLepP4CMS.rho();

  return -(massB0*massB0 + hadLepMass*hadLepMass - 2*_energyBCMS*hadLepE)/
    (2*pB0*hadLepP);

}

// Get now the 4 kin vars for the FF fit, (see .hh file for a description of each of them)

double BtaSemilepCand::cosThetaV() const {
  return cos(_dstarP4CMS.angle(_dP4DstarFrame.vect()));
}

double BtaSemilepCand::cosThetaL() const {
  if (boostOkFlag)
    {
      return cos((-_hadronP4CMS).angle(_leptonP4WFrame.vect()));
    }
  else return -2;
}

double BtaSemilepCand::angleChi() const {
  int flag = 0;
  double cthv = cos(_dstarP4CMS.angle(_dP4DstarFrame.vect()));
  double cthl =  cos((-_dstarP4CMS).angle(_leptonP4WFrame.vect()));

  Hep3Vector dP3Trans, leptonP3Trans;

  if (_dstarP4CMS.rho()!=0)  {
    flag=flag+1;
    dP3Trans= _dP4DstarFrame.vect() - cthv*_dP4DstarFrame.rho() * _dstarP4CMS.vect()*(1/_dstarP4CMS.rho())  ;
  }

  if (_pW.rho()!=0)  {
    flag=flag+1;
    leptonP3Trans= _leptonP4WFrame.vect() - cthl* _leptonP4WFrame.rho() * _pW.vect()* (1 / _pW.rho())  ;
  }

  double chi=-2;	
  if (flag==2)  chi = dP3Trans.angle(leptonP3Trans);
  return chi;

}

double BtaSemilepCand::momentumTransfer2() const {
  double massB0= Pdt::lookup("B0")->mass();
  double EHadron = _hadronP4CMS.e();
  double massHadron= _hadron->pdtEntry()->mass();
  return (massB0*massB0 + massHadron*massHadron - 2*massB0*EHadron);
}



double BtaSemilepCand::deltaM() const {
  if(0 == _dstar || 0 == _theD) return -1;
  return _dstar->mass() - _theD->mass();
}

double BtaSemilepCand::dDalitzWgt() const {
  if(0 == _theD) return -2;
  if ( _theD->nDaughters() != 3) return -1;
  if ( !_theK  || !_thePi1 || !_thePi2) return -1;

  double pi180inv= 3.141592653589793238/180.;
  //  double mP= _theD->mass();
  double mP= (_theK->p4()+ _thePi1->p4()+ _thePi2->p4()).m();
  //double mPPDG= _theD->pdtEntry()->mass();
  double mPi= Pdt::lookup("pi+")->mass();
  double mPi0= Pdt::lookup("pi0")->mass();
  // K is always #1
  //  double m1= _theK->pdtEntry()->mass();
  double m1= _theK->mass();
  double m2(0.0),m3(0.0),msq12(0.0),msq23(0.0),msq31(0.0);


  double wgt(0.0),maxWgt(1.0);
  DifComplex resSum,res1,res2,res3,res4,res5,res6;
  DifNumber resSumReal, resSumImag, resSumMagSq;

  switch (_dMode) {
  case (D0ToKpipi0):
    // charged pi should be #2
    m2= mPi;
    m3= mPi0;
    if (_thePi1->charge() != 0 ) {
      msq12= (_theK->p4() + _thePi1->p4()).m2();
      msq23= (_thePi1->p4() + _thePi2->p4()).m2();
      msq31= (_thePi2->p4() + _theK->p4()).m2();
    }
    else {
      msq12= (_theK->p4() + _thePi2->p4()).m2();
      msq23= (_thePi2->p4() + _thePi1->p4()).m2();
      msq31= (_thePi1->p4() + _theK->p4()).m2();
    }
    // K*0   (12)
    res1= brWigAmp(msq12,msq23,mP,m1,m2,m3,0.0862,-2.,0.0507,0.8961,1);
    // Rho    (23)
    res2= brWigAmp(msq23,msq31,mP,m2,m3,m1,0.1942,0.,0.1502,0.7693,1);
    // K*-    (13)
    res3= brWigAmp(msq31,msq23,mP,m1,m3,m2,0.0813,162,0.0508,0.8917,1);
    // nonresonance
    res4= 0.1784* DifComplex(cos(-122*pi180inv),sin(-122*pi180inv));
    resSum= res1+ res2+ res3+ res4;

    resSumReal = resSum.real();
    resSumImag = resSum.imag();
    resSumMagSq = resSumReal*resSumReal + resSumImag*resSumImag;
    wgt = resSumMagSq.number();
    maxWgt= 3.9;
    break;

  case (D0ToKspipi):
    m2= mPi;
    m3= mPi;
    // if D0, the order should be K_SO pi+ pi-
    // if anti-D0, it should be   K_S0 pi- pi+
    if (_theD->pdtEntry()== Pdt::lookup("D0"))
      if (_thePi1->charge() > 0 ) {
        msq12= (_theK->p4() + _thePi1->p4()).m2();
        msq23= (_thePi1->p4() + _thePi2->p4()).m2();
        msq31= (_thePi2->p4() + _theK->p4()).m2();
      } else {
        msq12= (_theK->p4() + _thePi2->p4()).m2();
        msq23= (_thePi2->p4() + _thePi1->p4()).m2();
        msq31= (_thePi1->p4() + _theK->p4()).m2();
      }
    else if (_theD->pdtEntry()== Pdt::lookup("anti-D0"))
      if (_thePi1->charge() < 0 ) {
        msq12= (_theK->p4() + _thePi1->p4()).m2();
        msq23= (_thePi1->p4() + _thePi2->p4()).m2();
        msq31= (_thePi2->p4() + _theK->p4()).m2();
      } else {
        msq12= (_theK->p4() + _thePi2->p4()).m2();
        msq23= (_thePi2->p4() + _thePi1->p4()).m2();
        msq31= (_thePi1->p4() + _theK->p4()).m2();
      }
    else {
      wgt= -1;
      maxWgt= 1;
      ErrMsg(trace) << "Wrong D0 mode! " << endmsg;
      break;
    }

    // K*(892)  (13)
    res1= brWigAmp(msq31,msq23,mP,m1,m3,m2,0.1703,0.,0.0508,0.8917,1);
    // K*0(1430)  (13)
    res2= brWigAmp(msq31,msq23,mP,m1,m3,m2,0.1222,-166.,0.294,1.412,0);
    // Rho  (23)
    res3= brWigAmp(msq23,msq31,mP,m2,m3,m1,0.1324,-136.,0.1502,0.7693,1);
    // f0(975)  (23)
    res4= brWigAmp(msq23,msq31,mP,m2,m3,m1,0.0283,38.,0.07,0.98,0);
    // f2(1270)  (23)
    res5= brWigAmp(msq23,msq31,mP,m2,m3,m1,0.6032,-174.,0.1851,1.2754,2);
    // f0(1400)  (23)
    res6= brWigAmp(msq23,msq31,mP,m2,m3,m1,0.1622,-45.,0.35,1.4,0);
    resSum= res1+ res2+ res3+ res4+ res5+ res6;

    resSumReal = resSum.real();
    resSumImag = resSum.imag();
    resSumMagSq = resSumReal*resSumReal + resSumImag*resSumImag;
    wgt = resSumMagSq.number();
    maxWgt= 8.1;
    break;

  case (DchToKpipi):
    if ( _thePi1->charge() != _thePi2->charge() ) {
      wgt= -1;
      maxWgt= 1;
      ErrMsg(trace) << "Wrong D0 mode! " << endmsg;
      break;
    }
    m1= mPi;
    m2= mPi;
    msq12= (_theK->p4() + _thePi1->p4()).m2();
    msq23= (_thePi1->p4() + _thePi2->p4()).m2();
    msq31= (_thePi2->p4() + _theK->p4()).m2();

    // K*(892)  (12+23)
    res1= brWigAmp(msq12,msq23,mP,m1,m2,m3,0.1131,48.,0.0507,0.8961,1);
    res2= brWigAmp(msq31,msq23,mP,m1,m3,m2,0.1131,48.,0.0507,0.8961,1);
    // K*(1430)
    res3= brWigAmp(msq12,msq23,mP,m1,m2,m3,0.2389,63.,0.294,1.412,0);
    res4= brWigAmp(msq31,msq23,mP,m1,m3,m2,0.2389,63.,0.294,1.412,0);
    // K*(1680)
    res5= brWigAmp(msq12,msq23,mP,m1,m2,m3,1.341,73.,0.322,1.717,1);
    res6= brWigAmp(msq31,msq23,mP,m1,m3,m2,1.341,73.,0.322,1.717,1);
    resSum= 0.5586*DifComplex(1.,0.)+
      0.5*(res1+ res2+ res3+ res4+ res5+ res6);

    resSumReal = resSum.real();
    resSumImag = resSum.imag();
    resSumMagSq = resSumReal*resSumReal + resSumImag*resSumImag;
    wgt = resSumMagSq.number();

    maxWgt= 2.6;
    break;

  case(DchToKspipi0):
    // charged pi should be #2
    m2= mPi;
    m3= mPi0;
    if (_thePi1->charge() != 0 ) {
      msq12= (_theK->p4() + _thePi1->p4()).m2();
      msq23= (_thePi1->p4() + _thePi2->p4()).m2();
      msq31= (_thePi2->p4() + _theK->p4()).m2();
    }
    else {
      msq12= (_theK->p4() + _thePi2->p4()).m2();
      msq23= (_thePi2->p4() + _thePi1->p4()).m2();
      msq31= (_thePi1->p4() + _theK->p4()).m2();
    }
    // Rho  (23)
    res1= brWigAmp(msq23,msq31,mP,m2,m3,m1,0.1818,0.,0.1502,0.7693,1); 
    // K*0  (31)
    res2= brWigAmp(msq31,msq12,mP,m3,m1,m2,0.0921,-137.,0.0507,0.8961,1); 
    resSum= 0.2018* DifComplex(cos(250*pi180inv),sin(250*pi180inv))+
      res1+ res2;

    resSumReal = resSum.real();
    resSumImag = resSum.imag();
    resSumMagSq = resSumReal*resSumReal + resSumImag*resSumImag;
    wgt = resSumMagSq.number();
    maxWgt= 3.6;
    break; 

  default:
    wgt= 1;
    maxWgt= 1;
    break;
  }
  return wgt/maxWgt;
}


  // copy most of the code from EvtGen/util/EvtResonance.cc
  // for (12) resonance
DifComplex BtaSemilepCand::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 pi180inv= 3.141592653589793238/180.;
  
  double m12= sqrt(msq12);
  double msqP= mP*mP;

  // 3-magnitude of p2 and p3 at 12 rest frame 
  double p2= sqrt((msq12-(m1+m2)*(m1+m2)) * (msq12-(m1-m2)*(m1-m2)))/(2.0*m12);
  double p3= sqrt((msqP-(m12+m3)*(m12+m3))*(msqP-(m12-m3)*(m12-m3)))/(2.0*m12);
  // p2 when m12= bwm
  double p2R = sqrt((bwm*bwm-(m1+m2)*(m1+m2))*(bwm*bwm-(m1-m2)*(m1-m2)))/(2.0*bwm);
  
  // angle 3 make with 2 in rest frame of (12)
  double cos23= decayAngle(msq12,msq23,mP,m1,m2,m3);
  //angle 3 makes with 1 in 12 is, of course, -cos23

  double gam(0), R(0);
  R= 1.5; // Blatt and Weiskopf Meson Radii (GeV-1)
  double p2rel= p2/p2R;
  double Rp2sq= R*p2*R*p2;
  double Rp2Rsq= R*p2R*R*p2R;

  // Blatt and Weiskopf Form factor squares
  double Fsq(0),FRsq(0);
  // angular dependence
  double S;

  switch (spin) {
  case 0:
    Fsq= 1;
    gam= gamma*(bwm/m12)*p2rel;
    S= 1;

    break;

  case 1:
    Fsq= 1.0/(1 + Rp2sq);
    FRsq= 1.0/(1 + Rp2Rsq);
    gam= gamma*(bwm/m12)* p2rel*p2rel*p2rel* Fsq/FRsq;
    S= 2.0* p2*p3*cos23;
    break;

  case 2:
    Fsq= 1.0/(9+ 3*Rp2sq+ Rp2sq*Rp2sq);
    FRsq= 1.0/(9+ 3*Rp2Rsq+ Rp2Rsq*Rp2Rsq);
    gam= gamma*(bwm/m12)*(p2/p2R)*(p2/p2R)*(p2/p2R)*(p2/p2R)*(p2/p2R)*Fsq/FRsq;
    S= 2.0*p2*p2*p3*p3*(3.0*cos23*cos23 -1) ;
    break;

  default:
    ErrMsg(warning) << "BtaSemilepCand: wrong spin ! "<< endmsg;
    S= 0.;
    break;         
  }

  DifComplex ampBW= ampl* DifComplex(cos(theta*pi180inv),sin(theta*pi180inv))*
    S* sqrt(Fsq)/(bwm*bwm- msq12 - DifComplex(0.0,1.0)*gam*bwm );

  return ampBW;

}

// cosTh(32) in (12) rest frame
double BtaSemilepCand::decayAngle(double msq12, double msq23, double mP,
                                  double m1, double m2, double m3) 
  const {
  double msqP= mP*mP;
  double msq1= m1*m1;
  double msq2= m2*m2;
  double msq3= m3*m3;
  double e2= (msq12+ msq2 -msq1)/2./sqrt(msq12);
  double e3= (msqP- msq12- msq3)/2./sqrt(msq12);
  double p2= sqrt(e2*e2- msq2);
  double p3= sqrt(e3*e3- msq3);
  double cost= (msq2+ msq3+ 2.*e2*e3- msq23)/2./p2/p3;
  return cost;
}

int BtaSemilepCand::getKaonId(const string& key) {
  if (!_event) {
    ErrMsg(error) << "Event cache is not set! Do yourDstarl.fillEventCache() first. " << endmsg;
    return -1;
  }

  if (_theK->charge() == 0) return -1;

  if ( !_kaonBitMap || _kaonMapKey != key) {
    //    ErrMsg(trace) << "Get a new kaonBitMap with key " << key << endmsg;
    _kaonBitMap= 
      Ifd< std::map<BtaCandidate*, int> >::get(_event,
							   IfdStrKey(key.c_str()));
  }
  if ( !_kaonBitMap) {
    ErrMsg(error) << "Config error: "<< key << " not available." << endmsg;
    return -1;
  } else 
    _kaonMapKey= key;

  // locate kaon in kaon bit map
  int bitmask=0;
  std::map<BtaCandidate*, int>::iterator iter = _kaonBitMap->find(_kaonTrk);
  if (iter != _kaonBitMap->end()) {
    bitmask = iter->second;
  }
 
  return bitmask;
}

int BtaSemilepCand::getElecId(const string& key) {
  if (!_event) {
    ErrMsg(error) << "Event cache is not set! Do yourDstarl.fillEventCache() first. " << endmsg;
    return -1;
  }

  if ( !_elecBitMap || _elecMapKey!= key) {
    //    ErrMsg(trace) << "Get a new elecBitMap with key " << key << endmsg;
    _elecBitMap= 
      Ifd< std::map<BtaCandidate*, int> >::get(_event,
							   IfdStrKey(key.c_str()));
  }
  if ( !_elecBitMap) {
    ErrMsg(error) << "Config error: " << key << " not available." << endmsg;
    return -1;
  } else
    _elecMapKey= key;

  // locate elec in elec bit map
  int bitmask=0;
  std::map<BtaCandidate*, int>::iterator iter = _elecBitMap->find(_leptonTrk);
  if (iter != _elecBitMap->end()) {
    bitmask = iter->second;
  }

 
  return bitmask;
}

int BtaSemilepCand::getMuonId(const string& key) {
  if (!_event) {
    ErrMsg(error) << "Event cache is not set! Do yourDstarl.fillEventCache() first. " << endmsg;
    return -1;
  }

  if ( !_muonBitMap || _muonMapKey!= key) {
    //    ErrMsg(trace) << "Get a new muonBitMap with key "<< key << endmsg;
    _muonBitMap= 
      Ifd< std::map<BtaCandidate*, int> >::get(_event,
							   IfdStrKey(key.c_str()));
  }
  if ( !_muonBitMap) {
    ErrMsg(error) << "Config error: " << key << " not available." << endmsg;
    return -1;
  } else 
    _muonMapKey= key;
  
  // locate elec in muon bit map
  int bitmask=0;
  std::map<BtaCandidate*, int>::iterator iter =  _muonBitMap->find(_leptonTrk);
  if (iter != _muonBitMap->end()) {
    bitmask = iter->second;
  }

  return bitmask;
}


void BtaSemilepCand::loadBeamSpot() {
  const BbrPointErr bs= gblEnv->getPep()->pepBeamSpotCal()->beamSpot();
  _beamSpot= HepPoint(bs.x(),bs.y(),bs.z());
}


bool BtaSemilepCand::fillMcCache(BtaMcAssoc* theMcMap) {
  if ( 0== theMcMap) {
    ErrMsg(fatal) << "no MC truth!" << endmsg;
    return false;
  }

  if ( _mcResolved) return true;
  if ( _truthMap) return true;
  _truthMap= theMcMap;
  return true;
}

bool BtaSemilepCand::resolveMcAssoc() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return true;

  // reset
  _realD= _realDstar= _dstarFromB= _leptFromB= _leptFromSameB= _directLept= _dFromB= _hadronFromB= _realHadron= _dlnuExclusive= false;
  _mcLeptId= _mcLeptBId= 0;
  _theMcD= _theMcDstar= _theMcSoft= _theMcLept= _theMcDstB= _theMcLB= _theMcDB=  _theMcHadron= _theMcHadB= 0;

  if ( 0== _truthMap ) {
    ErrMsg(fatal) << "truthMap is null !! Do fillMcCache(..)"
                  << endmsg; 
    return false;
  }
  _realD= realD();
  _realHadron= realHadron();
  _realDstar= realDstar();
  _dstarFromB= dstarFromB();
  _hadronFromB= hadronFromB();
  _dFromB= dFromB();
  _dlnuExclusive= dlnuExclusive();
  testMcLept();
  _mcLeptBId= mcLeptBId();

  _mcResolved= true;
  return _mcResolved;
}

bool BtaSemilepCand::findMother(BtaCandidate *dau, BtaCandidate *moth) {
  if ( !_truthMap) return false;

  BtaCandidate* aMoth;
  BtaCandidate* aCand(dau);
  while ((aMoth= aCand->theMother()) && aMoth!= 0) {
    if ( aMoth->uid() == moth->uid() ) return true;
    else aCand= aMoth;
  }

  return false;

}

BtaCandidate* BtaSemilepCand::findMother(BtaCandidate *dau, const char *moth) {
  if ( !_truthMap) return 0;

  BtaCandidate* aMoth(0);
  BtaCandidate* aCand(dau);
  bool found= false;
  while ((aMoth= aCand->theMother()) && aMoth!= 0) {
    if ( aMoth->pdtEntry() == Pdt::lookup(moth) ) {
      found= true;
      break;
    } else aCand= aMoth;
  }
  
  if (found) return aMoth;
  else return 0;

}

bool BtaSemilepCand::realHadron() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _realHadron;
  if ( _dstar) {
    return realDstar();
  } else {
    return realD();
  }
}

bool BtaSemilepCand::realD() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _realD;
  
  _realD= false;
  
  if ( 0== _theK || 0==_theD ) return false;

  BtaCandidate *theMcK= _truthMap->mcFromReco(_theK);
  if (0== theMcK) return false;
  if (theMcK->pdtEntry()!= _theK->pdtEntry()) return false;
  // find the Mc D from Mc K's mother
  _theMcD= findMother( theMcK, _theD->pdtEntry()->name());
  if ( 0== _theMcD) return false;
  if (_theMcD->pdtEntry()!=_theD->pdtEntry()) return false;

  bool match= false;
  if ( 0== _thePi1) return false;
  BtaCandidate *theMcPi1= _truthMap->mcFromReco(_thePi1);
  if (0== theMcPi1) return false;
  if (theMcPi1->pdtEntry()!= _thePi1->pdtEntry()) return false;
  match= findMother( theMcPi1, _theMcD);

  if ( _theD->nDaughters() >= 3) {
    if (0== _thePi2) return false;
    BtaCandidate *theMcPi2= _truthMap->mcFromReco(_thePi2);
    if (0== theMcPi2) return false;
    if (theMcPi2->pdtEntry()!= _thePi2->pdtEntry()) return false;
    match= findMother( theMcPi2, _theMcD);
  }
  if ( _theD->nDaughters() == 4) {
    if (0== _thePi3) return false;
    BtaCandidate *theMcPi3= _truthMap->mcFromReco(_thePi3);
    if (0== theMcPi3) return false;
    if (theMcPi3->pdtEntry()!= _thePi3->pdtEntry()) return false;
    match= findMother( theMcPi3, _theMcD);
  }

  // are the number of final particles the same 
  if (match) {
    bool print=false;
    int nFinal=nFinalDaus(*_theMcD,"pi+ pi0 K_S0 K+ e+ mu+ gamma",print);
    match = ( nFinal == _theD->nDaughters() );
  }
  _realD= match;

  return _realD;
}

bool BtaSemilepCand::realDstar() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _realDstar;

  _realDstar= false;

  if ( 0== _dstar || 0== _theSoft ) return false;

  // do we have an mc D?
  if ( 0==_theMcD ) {
    // test real D again
    _realD= realD();
    if ( 0==_theMcD || !_realD ) return false;    // still no mc D or D is not real
  }
  // real D has been tested
  if ( !_realD) return false;   

  _theMcSoft= _truthMap->mcFromReco(_theSoft);
  if ( 0==_theMcSoft) return false;

  _theMcDstar= findMother(_theMcD, _dstar->pdtEntry()->name());
  if ( 0==_theMcDstar) return false;

  _realDstar= findMother(_theMcSoft, _theMcDstar);

  return _realDstar;

}

bool BtaSemilepCand::dstarFromB() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _dstarFromB;

  _dstarFromB= false;
  if ( 0==_dstar) return false;

  // do we have an mc Dstar
  if ( 0== _theMcDstar) {
    // test real D* again
    _realDstar= realDstar();
    if ( 0== _theMcDstar || !_realDstar ) return false;  // still no mc D* or D* not real
  }
  // real D* has been tested
  if ( !_realDstar) return false;

  _theMcDstB= findMother(_theMcDstar, "B0");
  if ( 0== _theMcDstB) _theMcDstB= findMother(_theMcDstar, "anti-B0");
  if ( 0== _theMcDstB) _theMcDstB= findMother(_theMcDstar, "B+");
  if ( 0== _theMcDstB) _theMcDstB= findMother(_theMcDstar, "B-");
  
  if ( 0== _theMcDstB) _dstarFromB= false;
  else _dstarFromB= true;

  return _dstarFromB;

}

bool BtaSemilepCand::hadronFromB() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _hadronFromB;

  _hadronFromB= false;
  if ( 0==_hadron) return false;

  // do we have an mc Hadron?
  if ( _dstar) {
    _theMcHadron = _theMcDstar;
  } else {
    _theMcHadron = _theMcD;
  }

  if ( 0== _theMcHadron) {
    // test real D* again
    _realHadron= realHadron();
    if ( 0== _theMcHadron || !_realHadron ) return false;  // still no mc D* or D* not real
  }
  // real D* has been tested
  if ( !_realHadron) return false;

  _theMcHadB= findMother(_theMcHadron, "B0");
  if ( 0== _theMcHadB) _theMcHadB= findMother(_theMcHadron, "anti-B0");
  if ( 0== _theMcHadB) _theMcHadB= findMother(_theMcHadron, "B+");
  if ( 0== _theMcHadB) _theMcHadB= findMother(_theMcHadron, "B-");
  
  if ( 0== _theMcHadB) _hadronFromB= false;
  else _hadronFromB= true;

  return _hadronFromB;

}

bool BtaSemilepCand::dFromB() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _dFromB;

  _dFromB= false;
  if ( 0==_theD) return false;

  // do we have an mc D
  if ( 0== _theMcD) {
    // test real D again
    _realD= realD();
    if ( 0== _theMcD || !_realD ) return false;  // still no mc D or D not real
  }
  // real D has been tested
  if ( !_realD) return false;

  _theMcDB= findMother(_theMcD, "B0");
  if ( 0== _theMcDB) _theMcDB= findMother(_theMcD, "anti-B0");
  if ( 0== _theMcDB) _theMcDB= findMother(_theMcD, "B+");
  if ( 0== _theMcDB) _theMcDB= findMother(_theMcD, "B-");
  
  if ( 0== _theMcDB) _dFromB= false;
  else _dFromB= true;

  return _dFromB;

}

bool BtaSemilepCand::dlnuExclusive() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _dlnuExclusive;

  _dlnuExclusive= false;
  if ( 0==_theD) return false;

  // do we have an mc D
  if ( 0== _theMcD) {
    // test real D again
    _realD= realD();
    if ( 0== _theMcD || !_realD ) return false;  // still no mc D or D not real
  }
  // real D has been tested
  if ( !_realD) return false;

  if ( !_theMcD->theMother() || 
       (_theMcD->theMother()->pdtEntry() != Pdt::lookup("B0") &&
	_theMcD->theMother()->pdtEntry() != Pdt::lookup("anti-B0") &&
	_theMcD->theMother()->pdtEntry() != Pdt::lookup("B+") &&
	_theMcD->theMother()->pdtEntry() != Pdt::lookup("B-") ) ||
       _theMcD->theMother()->nDaughters()!=3 )
    return false;

  // check lepton
  
  _dlnuExclusive = leptFromSameB();
  return _dlnuExclusive;

}

bool BtaSemilepCand::leptFromB() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _leptFromB;
  else testMcLept();

  return _leptFromB;
}

bool BtaSemilepCand::leptFromSameB() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _leptFromSameB;
  else testMcLept();

  return _leptFromSameB;
}

bool BtaSemilepCand::directLept() {
  if ( !_truthMap) return false;
  if ( _mcResolved) return _directLept;
  else testMcLept();

  return _directLept;
}

int BtaSemilepCand::mcLeptId() {
  if ( !_truthMap) return 0;
  if ( _mcResolved) return _mcLeptId;
  else testMcLept();

  return _mcLeptId;
}
  
void BtaSemilepCand::testMcLept() {
  if ( !_truthMap) return;
  if ( _mcResolved) return;

  _leptFromB= false;
  _leptFromSameB= false;
  _directLept= false;
  _mcLeptId= 0;
  
  if ( 0== _lepton ) return;
  _theMcLept= _truthMap->mcFromReco(_lepton);
  if ( 0== _theMcLept) return;

  // check true lepton ID
  _mcLeptId= _theMcLept->pdtEntry()->pdgId();

  // check if lepton cand from a B
  _theMcLB= findMother(_theMcLept, "B0");
  if ( 0== _theMcLB ) _theMcLB= findMother(_theMcLept, "anti-B0");
  if ( 0== _theMcLB ) _theMcLB= findMother(_theMcLept, "B+");
  if ( 0== _theMcLB ) _theMcLB= findMother(_theMcLept, "B-");
  if ( 0!= _theMcLB ) _leptFromB= true;
  else return;

  // now lepton is from a B, check if it is a direct lepton
  if ( _theMcLept->theMother()->uid()== _theMcLB->uid()) _directLept= true;

  // check if this B is the same as B->D* or D
  if ( 0== _theMcHadB) hadronFromB();
  if ( 0== _theMcHadB || !_hadronFromB ) return;
  if ( _theMcHadB->uid() == _theMcLB->uid()) _leptFromSameB= true;

  return;

}

int BtaSemilepCand::mcLeptBId() {
  if ( !_truthMap) return 0;
  if ( _mcResolved) return _mcLeptBId;

  if ( 0== _theMcLB) _mcLeptBId=0;
  else _mcLeptBId=_theMcLB->pdtEntry()->pdgId();
  return _mcLeptBId;

}

int BtaSemilepCand::nFinalDaus(const BtaCandidate& c, std::string terminals, bool print) {
  int nFinal(0);
  fillnFinal(c,terminals,&nFinal, print);
  return nFinal;

}

void BtaSemilepCand::fillnFinal( const BtaCandidate& c, std::string terminals, int* nFinal, bool print){
  bool stable= (c.nDaughters()==0);

  const PdtEntry* cEntry= c.pdtEntry();
  int cPdgId(0);
  if (cEntry) cPdgId= cEntry->pdgId();
  
  if ( cPdgId != 0 && !stable) {
    std::istringstream terminalsStream(terminals);
    std::string terminal;

    while (terminalsStream >> terminal ){
      const PdtEntry *entry= Pdt::lookup(terminal);
      int pdgId(0);
      if (entry) pdgId= entry->pdgId();
      else ErrMsg(fatal) << "Bad terminal name: " << terminal << endmsg;

      if (abs(cPdgId)== abs(pdgId)) stable= true;
    }
  }

  if (ErrLogging(trace) || print)
    cout << "BtaSemilepCand.cc::fillnFinal: " << cEntry->name() << ": " ;

  if (!stable) {
    BtaCandidate* dau(0);
    HepAListIterator<BtaCandidate> iterDau=c.daughterIterator();
    if (ErrLogging(trace) || print) {
      while( (dau=iterDau()) ) 
        cout << dau->pdtEntry()->name() << " ";
      cout << endl;
      iterDau.rewind();
    }
    while( (dau=iterDau()) ) 
      fillnFinal( *dau ,terminals, nFinal, print);
  } else
    (*nFinal)++;
  
}