//--------------------------------------------------------------------------
// File and Version Information:
//      MyBToDlnuAnalysis.cc 2004/10/11
//      based on
//
// Description:
//      Semileptonic B -> Dlnu decay analysis
//
// Environment:
//      Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//      Kenji Hamano
//
// Copyright Information:
//      Copyright (C) 
//
//------------------------------------------------------------------------

#include "BaBar/BaBar.hh"

//-----------------------
// This Class's Header --
//-----------------------
#include "BToDlnuMiniUser/MyBToDlnuAnalysis.hh"
	
//-------------
// C Headers --
//-------------
#include <assert.h>

//---------------
// C++ Headers --
//---------------
//#include <iostream>
#include <math.h>

//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "AbsEnv/AbsEnv.hh"
#include "AbsEvent/AbsEvent.hh"
#include "Beta/BtaCandidate.hh"
#include "Beta/BtaMcTruth.hh"
#include "Beta/EventInfo.hh"
#include "BetaCoreTools/BtaOpMakeTree.hh"
#include "BetaCoreTools/BtaTreeNavigator.hh"
#include "VtxBase/BtaAbsVertexer.hh"
#include "VtxBase/VtxAbsAlgorithm.hh"
#include "VtxFitter/VtxFitterOper.hh"
//#include "Beta/BtaAbsVertex.hh"
//#include "FastVtx/BtaOpFastVtx.hh"
#include "ProbTools/probab.hh"
#include "BetaCoreTools/BtaMcAssoc.hh"
//#include "BetaCoreTools/BtaCandMap2.hh"
#include "BetaTools/BtaSemilepCand.hh"
#include "CLHEP/Alist/AList.h"
#include "CLHEP/Alist/AIterator.h"
#include "ErrLogger/ErrLog.hh"
#include "GenEnv/GenEnv.hh"
#include "HepTuple/Tuple.h"
//#include "HepTuple/HepHistID.hh"
#include "HepTuple/TupleManager.h"
#include "HepTuple/Histogram.h"
#include "HepTuple/HTValOrderedVector.h"
#include "PDT/Pdt.hh"
#include "PDT/PdtEntry.hh"
#include "PDT/PdtLund.hh"
#include "ProxyDict/IfdStrKey.hh"
#include "ProxyDict/IfdKey.hh"
#include "ProxyDict/Ifd.hh"
#include "ProbTools/probab.hh"

#include <string>
#include <map>
using std::endl;
using std::ostream;

//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------

//----------------
// Constructors --
//----------------

MyBToDlnuAnalysis::MyBToDlnuAnalysis( const char* const theName, 
				  const char* const theDescription )
  : AppFilterModule( theName, theDescription )
  , _eventInfoList(        "eventInfoList",        this, "Default"  )
  , _btaTruthList(         "btaTruthList",         this, "MCTruth" )
  , _chargedList(          "chargedList",          this, "ChargedTracks")
  , _GTVLList(             "GTLVList",             this, "GoodTracksVeryLoose")
  , _GTLList(              "GTLList",              this, "GoodTracksLoose")
  , _photonList(           "photonList",           this, "CalorNeutral")
  , _gammaConvList(        "gammaConvList",    this, "gammaConversionDefault")
  , _electronList(         "electronList",         this, "PidLHElectrons")
  , _muonList(             "muonList",             this, "muNNLoose")
  , _muonTightList(        "muonTightList",        this, "muNNTight")
  , _kaonList(             "kaonList",             this, "KMicroNotPionGTL")
  , _compDDstarLList(      "compDDstarLList",      this, "BToDlnuSkimmedList")
  , _btaTruthMap(          "btaTruthMap",          this, "GHit")
  //, _vtxFitAlgo(           "vtxFitAlgo", this, VtxAbsAlgorithm::Cascade)
  , _requireBhabhaVeto(    "requireBhabhaVeto",    this, true)
  , _requireDOnly(         "requireDOnly",         this, true)
  , _requireTightK(        "requireTightK",        this, true)
  , _requireDVtxCut(       "requireDVtxCut",       this, true)
  , _fitDprobMin(          "fitDprobMin",          this, 0.01)
  , _requireBVtxCut(       "requireBVtxCut",       this, true)
  , _fitBprobMin(          "fitBprobMin",          this, 0.01)
  , _requreThrustCut(      "requreThrustCut",      this, true)
  , _cosThrMax(            "cosThrMax",            this, 0.88)
  , _requireTightMu(       "requireTightMu",       this, false)
  , _requireDlLepHighPCut( "requireDlLepHighPCut", this, false)
  , _DlLepMomMax(          "DlLepMomMax",          this, 2.8)
  , _booster()
  , _numberEvt(0)
  , _numPassedEvt(0)
{
  commands()->append( &_eventInfoList );
  commands()->append( &_btaTruthList );
  commands()->append( &_chargedList);
  commands()->append( &_GTVLList);
  commands()->append( &_GTLList);
  commands()->append( &_photonList);
  commands()->append( &_gammaConvList);
  commands()->append( &_electronList);
  commands()->append( &_muonList);
  commands()->append( &_muonTightList);
  commands()->append( &_kaonList);
  commands()->append( &_compDDstarLList);
  commands()->append( &_btaTruthMap );
  //commands()->append( &_vtxFitAlgo );
  commands()->append( &_requireBhabhaVeto);
  commands()->append( &_requireDOnly);
  commands()->append( &_requireTightK);
  commands()->append( &_requireDVtxCut);
  commands()->append( &_fitDprobMin);
  commands()->append( &_requireBVtxCut);
  commands()->append( &_fitBprobMin);
  commands()->append( &_requreThrustCut);
  commands()->append( &_cosThrMax);
  commands()->append( &_requireTightMu);
  commands()->append( &_requireDlLepHighPCut);
  commands()->append( &_DlLepMomMax);
}

//--------------
// Destructor --
//--------------

MyBToDlnuAnalysis::~MyBToDlnuAnalysis( )
{}

//--------------
// Operations --
//--------------

AppResult
MyBToDlnuAnalysis::beginJob( AbsEvent* anEvent )
{
    ErrMsg(routine)<<"begin Job"<<endmsg;

    HepTupleManager* manager = gblEnv->getGen()->ntupleManager();
    assert(manager != 0);

    // create output tuples
    _analysisTuple = manager->ntuple("BToDlnu Analysis ntuple");

    // Names starting with underscores are defined in the
    // MyBToDlnuAnalysis.hh file as class variables. They live for the
    // life of the analysis object (this one), so can be used to communicate
    // between the begin, event and end functions. If you want to book another
    // here, make another entry in the obvious way.

    _nTrials = 0;
    _nSuccesses = 0;

    _cutBhabha = 0;
    _cutDOnly = 0;
    _cutTightK = 0;
    _cutDVtx = 0;
    _cutBVtx = 0;
    _cutThrust = 0;
    _cutTightMu = 0;
    _cutDlLepHighP = 0;
    _cutCandRequirement = 0;
    _cutTooManyDl = 0;

  return AppResult::OK;
}

//------------------------------------------------------------

AppResult
MyBToDlnuAnalysis::event( AbsEvent* anEvent )
{
  _numberEvt++;
  bool result = pass( anEvent );
  setPassed( result );

  return AppResult::OK;
}

//-----------------------------------------------------------
  
AppResult
MyBToDlnuAnalysis::endJob( AbsEvent* anEvent )
{
  ErrMsg(routine)<<" end Job" << endmsg;

  // Print statistics
  if ( ErrLogging(routine) ) {
    ErrMsg(routine) << "===========================================" << endmsg;
    ErrMsg(routine) << " Select candidates selection in " << name()  << endmsg;
    ErrMsg(routine) << "===========================================" << endmsg;
    ErrMsg(routine) << "---- " 
                    << _numberEvt  << "\tevents processed" << endmsg;
    ErrMsg(routine) << "---- " << _cutBhabha    
                    << "\tfailed Bhabha cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutCandRequirement    
                    << "\tfailed Dl cand requirements below"   << endmsg;
    ErrMsg(routine) << "---- " << _cutTooManyDl    
                    << "\tfailed  too many Dl candidates cut"   << endmsg;
    ErrMsg(routine) << "---- " << _numPassedEvt  
		    << "\tevents passed" << endmsg;
    if(_numberEvt>0) {
      ErrMsg(routine) << "---- " << (float)_numPassedEvt/(float)_numberEvt 
		      << "\tFraction of events" << endmsg;
    }
    ErrMsg(routine) << "---- " 
                    << _nTrials    << "\ttotal Candidates" << endmsg;
    ErrMsg(routine) << "---- " << _cutDOnly    
                    << "\tfailed D only (not Dstar) cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutTightK    
                    << "\tfailed tight K cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutDVtx    
                    << "\tfailed D vertex cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutBVtx    
                    << "\tfailed B vertex cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutThrust    
                    << "\tfailed thrust angle cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutTightMu    
                    << "\tfailed tight Dl Muon cut"   << endmsg;
    ErrMsg(routine) << "---- " << _cutDlLepHighP    
                    << "\tfailed Dlnu lepton high momentum cut"   << endmsg;
    ErrMsg(routine) << "---- " 
                    << _nSuccesses << "\tselected Candidates" << endmsg;
    if(_nTrials>0) {
      ErrMsg(routine) << "---- " << (float)_nSuccesses/(float)_nTrials 
		      << "\tFraction of candidates" << endmsg;
    }
  }

  return AppResult::OK;
}

//------------------------------------------------------------------
bool
MyBToDlnuAnalysis::pass( AbsEvent* anEvent )
{

    ErrMsg(trace) << "start MyBToDlnuAnalysis::pass " << endmsg;
    // Get all the lists to be used.
    const HepAList<EventInfo>* infoList=
      Ifd<HepAList<EventInfo> >::get(anEvent,_eventInfoList.value());
    if (infoList == 0) ErrMsg(fatal) << "No EventInfo !!" << endmsg;
    EventInfo* eventInfo = infoList->first();

    const HepAList<BtaCandidate>* truthList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _btaTruthList.value());
    int nTruth(0);
    if (truthList != NULL) nTruth = truthList->length();

    HepAList<BtaCandidate>* chargedList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _chargedList.value());
    int nChgTrk(0);
    if (chargedList != NULL) nChgTrk = chargedList->length();

    HepAList<BtaCandidate>* GTVLList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _GTVLList.value());
    int nGTVLTrk(0);
    if (GTVLList != NULL) nGTVLTrk = GTVLList->length();

    HepAList<BtaCandidate>* GTLList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _GTLList.value());
    int nGTLTrk(0);
    if (GTLList != NULL) nGTLTrk = GTLList->length();

    HepAList<BtaCandidate>* photonList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _photonList.value());
    int nPhoton(0);
    if (photonList != NULL) nPhoton = photonList->length();

    HepAList<BtaCandidate>* gammaConvList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _gammaConvList.value());
    int nGamConv(0);
    if (gammaConvList != NULL) nGamConv = gammaConvList->length();

    const HepAList<BtaCandidate>* electronList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _electronList.value());
    int nElectron(0);
    if (electronList != NULL) nElectron = electronList->length();
 
    const HepAList<BtaCandidate>* muonList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _muonList.value());
    int nMuon(0);
    if (muonList != NULL) nMuon = muonList->length();

    const HepAList<BtaCandidate>* muonTightList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _muonTightList.value());
    int nTightMu(0);
    if (muonTightList != NULL) nTightMu = muonTightList->length();

    const HepAList<BtaCandidate>* kaonList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _kaonList.value());
    int nKaon(0);
    if (kaonList != NULL) nKaon = kaonList->length();

    const HepAList<BtaCandidate>* compDDstarLList=
      Ifd<HepAList<BtaCandidate> >::get(anEvent, _compDDstarLList.value());
    int nDlCand(0);
    if (compDDstarLList != NULL) nDlCand = compDDstarLList->length();
    
    BtaMcAssoc* truthMap = 
      Ifd< BtaMcAssoc >::get( anEvent, _btaTruthMap.value());
    bool mcTruth(truthMap != 0);


    // Define and calculate some constants
    HepLorentzVector p4CM = eventInfo->cmFrame();
    _booster.setLorentzVector(p4CM);

    //double pdgUpsMass= Pdt::lookup("Upsilon(4S)")->mass();
    //double mCM = p4CM.m();
    //double mB = 5.280;
    //if( mCM/2 < mB ) mB = mB*mCM/pdgUpsMass; // off-peak data
    //    double pBmag2 = mCM*mCM/4 - mB*mB;



    //##########################################
    //### Radiative Bhabha veto ############
    //######################################
    ErrMsg(trace) << "BhaBha veto" << endmsg;
    int isABhabha(0);
    int nGamma(0);  // number of photon with E>0.8 GeV
    if (nPhoton > 0) {
      HepAListIterator<BtaCandidate> photonIter(*photonList);
      BtaCandidate* photoCand(0);
      while (photoCand = photonIter()) {
	float photonE = photoCand->p4().e();
	if (photonE > 0.080) nGamma++;
      }
    }
    int nTrk = nGTLTrk;
    int nGam = nGamma;
    int nConv = nGamConv;
    int nTrk2 = nTrk - 2*nConv;
    int nGam2 = nGam + nConv;
    int mult2 = nTrk2 + nGam2/2;
    if (nTrk < 4 || mult2 < 5)  isABhabha=1;
    ErrMsg(trace) << "nGTLTrk = " << nGTLTrk << endmsg;
    ErrMsg(trace) << "nPhoton = " << nPhoton << endmsg;
    ErrMsg(trace) << "nGamma = " << nGamma << endmsg;
    ErrMsg(trace) << "nConv = " << nConv << endmsg;
    if (_requireBhabhaVeto.value() && isABhabha==1) {
      _cutBhabha++;
      return false;
    }

    // ###########################################################
    // ### Dlnu block ##################
    // #################################
    ErrMsg(trace) << "Dlnu block" << endmsg;
    // Quantities to be stored in the ntuple.
    // MC truth info
    HTValOrderedVector<int>   trueDVec; 
    HTValOrderedVector<int>   trueDhadVec;   
    HTValOrderedVector<int>   trueDstarVec; 
    HTValOrderedVector<int>   trueDFromBVec;
    HTValOrderedVector<int>   trueDlnuVec;   
    HTValOrderedVector<int>   trueDhadFromBVec;
    HTValOrderedVector<int>   trueDstarFromBVec;
    HTValOrderedVector<int>   trueLepFromBVec; //Lepton originated from B,
                                               //including B->D->X+l, etc.
    HTValOrderedVector<int>   trueLepFromSameBVec; //lepton from same B as D
    HTValOrderedVector<int>   trueDhadlnuVec; 
    HTValOrderedVector<int>   trueDstarlnuVec;
    HTValOrderedVector<int>   directLeptVec; //Lepton came directly from a B
    HTValOrderedVector<int>   mcLeptIdVec; // lepton ID
    HTValOrderedVector<int>   mcLeptBIdVec; // Mother B ID
    HTValOrderedVector<int>   trueDecModVec; // B decay mode from MC truth
    HTValOrderedVector<int>   trueDDecModVec; // D decay mode from MC truth
    HTValOrderedVector<int>   noExtraKorPiVec; // D decay mode from MC truth
    HTValOrderedVector<int>   mcKPifromSameDVec;
    HTValOrderedVector<int>   mcKMissIDVec;
    HTValOrderedVector<int>   mcPiMissIDVec;
    HTValOrderedVector<int>   mcKIDVec;
    HTValOrderedVector<int>   mcPi1IDVec;
    HTValOrderedVector<int>   mcPi2IDVec;
    HTValOrderedVector<int>   mcPi3IDVec;
    HTValOrderedVector<int>   stepsToKMomVec;
    HTValOrderedVector<int>   stepsToPi1MomVec;
    HTValOrderedVector<int>   stepsToPi2MomVec;
    HTValOrderedVector<int>   stepsToPi3MomVec;
    // Reconstructed info
    HTValOrderedVector<int>   DhadTypeVec; //D* =>1, D =>0 
    HTValOrderedVector<int>   DhadModeVec; // Lund ID of Dhadron (D* or D)
                                           //   D_plus       => +411
                                           //   D_minus      => -411
                                           //   D0           => +421
                                           //   anti_D0      => -421
                                           //   D_star_plus  => +413
                                           //   D_star_minus => -413
                                           //   D_star0      => +423
                                           //   anti_D_star0 => -423
    HTValOrderedVector<int>   DDecayModeVec; // D0/D+ decay mode
                                             //   1 => D0ToKpi
                                             //   2 => D0ToK3pi 
                                             //   3 => D0ToKpipi0
                                             //   4 => D0ToKspipi
                                             //   5 => DchToKpipi
                                             //   6 => DchToKspipi0
                                             //   7 => DchToKspi
    HTValOrderedVector<float> DhadChrgVec; // Charge of D hadron (D* or D)
    HTValOrderedVector<float> DChrgVec; // Charge of D
    HTValOrderedVector<float> DhadMassVec;  // Dhadron (D* or D) mass
    HTValOrderedVector<float> DMassVec;  // D mass 
    HTValOrderedVector<float> DhadMomVec; // Dhadron momentum in CM
    HTValOrderedVector<float> DMomVec; // D momentum in CM
    HTValOrderedVector<float> DlLepChrgVec;  // Dlnu side lepton charge
    HTValOrderedVector<int>   DlLepTypeVec;  //electron -> 1, muon -> 2
    HTValOrderedVector<float> DlLeptMomVec; // Dlnu side lepton momentum in CM
    HTValOrderedVector<float> DLcosBYVec;  // cosTheta_ B-Dl
    HTValOrderedVector<float> DlThrustVec; // Thrust axis of Dl
    HTValOrderedVector<float> DlThrustAxVec;
    HTValOrderedVector<float> DlThrustAyVec;
    HTValOrderedVector<float> DlThrustAzVec;
    HTValOrderedVector<float> nonDlThrustVec; //Thrust axis of non Dl particles
    HTValOrderedVector<float> nonDlThrustAxVec;
    HTValOrderedVector<float> nonDlThrustAyVec;
    HTValOrderedVector<float> nonDlThrustAzVec;
    HTValOrderedVector<float> cosDlnonDlVec; // angle of Dl and non-Dl
    HTValOrderedVector<int>   numChargedDlVec; // #charged trk in this Dl cand
    HTValOrderedVector<int>   numNeutralDlVec; // #neutrals in this Dl cand
    HTValOrderedVector<int>   numChargednonDlVec; //#non Dl charged trk
    HTValOrderedVector<int>   numNeutralnonDlVec; //#non Dl neutrals
    HTValOrderedVector<double> fitDchiSqVec; //chisquare of the D fit
    HTValOrderedVector<int> fitDnDofVec; // number of degree of freedom
    HTValOrderedVector<int> fitDstatVec; // status of the D fit
    HTValOrderedVector<double> DprobVec; //chisquare of the B fit
    HTValOrderedVector<double> fitBchiSqVec; //chisquare of the B fit
    HTValOrderedVector<int> fitBnDofVec; // number of degree of freedom
    HTValOrderedVector<int> fitBstatVec; // status of the B fit
    HTValOrderedVector<double> BprobVec; //chisquare of the B fit

    // Loop over Dlnu candidates and store information into ntuple
    int nSelectedDlCand(0);
    for (size_t i=0; i<nDlCand; i++) {
      ErrMsg(trace) << "looping over Dlnu cand" << endmsg;
      _nTrials++;  //number of Dlnu candidates
      BtaCandidate* cand = (*compDDstarLList)[i]; //compDDstarLList candidate
      BtaSemilepCand dlCand(*cand, p4CM);  //BtaSemilepCand candidate

      bool isDstar(dlCand.theDstar()!=0);
      int dDecayMode(dlCand.getDDecayMode());
      const int DLund(dlCand.theD()->pdtEntry()->lundId());

      // D Only cut
      ErrMsg(trace) << "Donly cut" << endmsg;
      if (_requireDOnly.value() && isDstar) {_cutDOnly++; continue;}
	

      // Tight K cut
      ErrMsg(trace) << "Tight K cut" << endmsg;
      if (_requireTightK.value()) {
	// Loop over kaonList and check theK is in the list
	ErrMsg(trace) << "looping over tight K list" << endmsg;
	bool foundKaon(false);
	for ( size_t j=0; j<nKaon; j++ ) {
	  BtaCandidate* kaonCand = (*kaonList)[j];
	  if (kaonCand->overlaps(*cand)) {
	    ErrMsg(trace) << "found a overlapped kaon" << endmsg;
	    foundKaon = true;
	  }
	} //end of kaon loop
	if (!foundKaon) {_cutTightK++; continue;}
      }
	  

      // Get D vertex info
      ErrMsg(trace) << "Vertexing" << endmsg;
      BtaOpMakeTree comb;
      BtaCandidate* D(0);
      double fitDchiSq(-1);
      int fitDnDof(-1);
      int fitDstat(-1);
      double Dprob(-1);
      static const PdtEntry* dzero = Pdt::lookup(PdtLund::D0) ;
      static const PdtEntry* dzeroBar = Pdt::lookup(PdtLund::anti_D0) ;
      static const PdtEntry* dplus = Pdt::lookup(PdtLund::D_plus) ;
      static const PdtEntry* dminus = Pdt::lookup(PdtLund::D_minus) ;

      if (dlCand.thePi1()!=0 && dlCand.theK()!=0) {
	ErrMsg(trace) << "found Pi1 and K" << endmsg;
	if (dDecayMode==1) {
	  ErrMsg(trace) << "dDecayMode=1" << endmsg;
	  D = comb.create(*dlCand.thePi1(),*dlCand.theK());
	  if (DLund == 421)  D->setType(dzero);
	  else D->setType(dzeroBar);
	} 
	if (dDecayMode==7) {
	  ErrMsg(trace) << "dDecayMode=7" << endmsg;
	  D = comb.create(*dlCand.thePi1(), *dlCand.theK());
	  if (DLund == 411)  D->setType(dplus);
	  if (DLund == -411) D->setType(dminus);
	}
	if (dlCand.thePi2() != 0) {
	  ErrMsg(trace) << "found Pi2" << endmsg;
	  if (dDecayMode==3) {
	    ErrMsg(trace) << "dDecayMode=3" << endmsg;
	    D = comb.create(*dlCand.thePi1(),*dlCand.thePi2(),*dlCand.theK());
	    if (DLund == 421)  D->setType(dzero);
	    else D->setType(dzeroBar);
	  }
	  if (dDecayMode==4) {
	    ErrMsg(trace) << "dDecayMode=4" << endmsg;
	    D = comb.create(*dlCand.thePi1(),*dlCand.thePi2(),*dlCand.theK());
	    if (DLund == 421)  D->setType(dzero);
	    else D->setType(dzeroBar);
	  }
	  if (dDecayMode==5) {
	    ErrMsg(trace) << "dDecayMode=5" << endmsg;
	    D = comb.create(*dlCand.thePi1(),*dlCand.thePi2(), *dlCand.theK());
	    if (DLund == 411)  D->setType(dplus);
	    else  D->setType(dminus);
	  }
	  if (dlCand.thePi3()!=0) {
	    ErrMsg(trace) << "found Pi3" << endmsg;
	    if (dDecayMode==2) {
	      ErrMsg(trace) << "dDecayMode=2" << endmsg;
	      D = comb.create(*dlCand.thePi1(),*dlCand.thePi2(),
			      *dlCand.thePi3(),*dlCand.theK());
	      if (DLund == 421)  D->setType(dzero);
	      else D->setType(dzeroBar);
	    }
	  }
	}
      }
      if (D == 0) {
	ErrMsg(trace) << "D = 0" << endmsg;
      }
      if (D!=0) {
	//setGeoConstraint(*D);
	VtxFitterOper dFitter(*D, VtxFitterOper::TreeFitter);
	dFitter.fit();
	BtaCandidate fitD = dFitter.getFitted(*D);
	if (fitD!=0) {
	  fitDchiSq = fitD.decayVtx()->chiSquared();
	  fitDnDof = fitD.decayVtx()->nDof();
	  fitDstat = fitD.decayVtx()->status();
	  Dprob = probab( fitDnDof, fitDchiSq );
	}
      }

      if (_requireDVtxCut.value() && Dprob<_fitDprobMin.value()) {
	_cutDVtx++;
	continue;
      }

      // Fit to B decay
      BtaCandidate* B(0);
      double fitBchiSq(-1);
      int fitBnDof(-1);
      int fitBstat(-1);
      double Bprob(-1);
      if (dlCand.theD()!=0 && dlCand.theLepton()!=0) {
	float  Bcharge 
	  = dlCand.theD()->charge() + dlCand.theLepton()->charge();
	B = comb.create(*dlCand.theD(),*dlCand.theLepton());
	if (B == 0) {
	  ErrMsg(trace) << "B = 0" << endmsg;
	  _cutDVtx++;
	}
	//if (B == 0) {_cutDVtx++; continue;}
	if (B!=0) {
	  static const PdtEntry* bplus = Pdt::lookup("B+") ;
	  static const PdtEntry* bzero = Pdt::lookup("B0") ;
	  B->setType( Bcharge==0 ? bzero : 
		    (Bcharge>0 ? bplus : bplus->conjugate()));
	  VtxFitterOper bFitter(*B, VtxFitterOper::TreeFitter);
	  bFitter.fit();
	  BtaCandidate fitB = bFitter.getFitted(*B);
	  if (fitB!=0) {
	    fitBchiSq = fitB.decayVtx()->chiSquared();
	    fitBnDof = fitB.decayVtx()->nDof();
	    fitBstat = fitB.decayVtx()->status();
	    Bprob = probab( fitBnDof, fitBchiSq );
	  }
	}
      }

      if (_requireBVtxCut.value() && Bprob<_fitBprobMin.value()) {
	_cutBVtx++;
	continue;
      }

      

      // Check theLepton (Dl lepton) is electron or muon
      ErrMsg(trace) << "Lepton flavour and tight nu cut" << endmsg;
      // Loop over electron list and check theLepton is in the list
      bool foundElec(false);
      for ( size_t j=0; j<nElectron; j++ ) {
	BtaCandidate* electronCand = (*electronList)[j];
	if (electronCand->overlaps(*cand)) {
	  ErrMsg(trace) << "found a overlapped Dl electron" << endmsg;
	  foundElec = true;
	}
      } // end of electron loop
      // Loop over tight muon list and check theLepton is in the list
      bool foundTightMu(false);
      for ( size_t j=0; j<nTightMu; j++ ) {
	BtaCandidate* tightMuCand = (*muonTightList)[j];
	if (tightMuCand->overlaps(*cand)) {
	  ErrMsg(trace) << "found a overlapped tight muon" << endmsg;
	  foundTightMu = true;
	}
      } // end of tight muon loop
      int DlLeptonType(0);
      if(foundElec) DlLeptonType=1;
      else if (foundTightMu) DlLeptonType=2;
      
      // Tight Dl muon cut
      if (_requireTightMu.value()) {
	if (!foundTightMu && !foundElec) {
	  _cutTightMu++;
	  continue;
	}
      }

      // Dlnu side lepton high mementum cut
      ErrMsg(trace) << "Dl lep high mom cut" << endmsg;
      if(_requireDlLepHighPCut.value()) {
	if (dlCand.leptonPstar() > _DlLepMomMax.value()) {
	  _cutDlLepHighP++;
	  continue;
	}
      }

      // Get thrust axis for this Dl candidate
      ErrMsg(trace) << "geting thrust axis" << endmsg;
      // Create, DlCandList and nonDlCandList
      HepAList<BtaCandidate>* DlCandList = new HepAList<BtaCandidate>;
      if (!DlCandList) ErrMsg(fatal) << "No DlCandList" << endmsg;
      HepAList<BtaCandidate>* nonDlCandList = new HepAList<BtaCandidate>;
      if (!nonDlCandList) ErrMsg(fatal) << "No nonDlCandList" << endmsg;
      int nChargedDl(0), nChargednonDl(0);
      int nNeutralDl(0), nNeutralnonDl(0);
      ErrMsg(trace) << "looping over charged list" << endmsg;
      HepAListIterator<BtaCandidate> ChrgIter(*chargedList);
      BtaCandidate* Chrgcand;
      while (Chrgcand = ChrgIter()) {
	if (Chrgcand->overlaps(*cand)) {
	  nChargedDl++;
	  DlCandList->append(*Chrgcand);
	} else {
	  nChargednonDl++;
	  nonDlCandList->append(*Chrgcand);
	}
      }
      ErrMsg(trace) << "looping over neutral list" << endmsg;
      HepAListIterator<BtaCandidate> NeutIter(*photonList);
      BtaCandidate* Neutcand;
      while (Neutcand = NeutIter()) {
	if (Neutcand->overlaps(*cand)) {
	  nNeutralDl++;
	  DlCandList->append(*Neutcand);
	} else {
	  nNeutralnonDl++;
	  nonDlCandList->append(*Neutcand);
	}
      }

      double DlThu(0), DlThuX(0), DlThuY(0), DlThuZ(0), DlThuMag(0);
      double noDlThu(0), noDlThuX(0), noDlThuY(0), noDlThuZ(0), noDlThuMag(0);
      double denDlnonDl(0);
      if (DlCandList==0 || nonDlCandList==0) {
	ErrMsg(trace) << "No DlCandList or nonDlCandList" << endmsg;
      }
      if (DlCandList!=0 && nonDlCandList!=0) {
	// Create BtaThrust object and get thrust axis
	BtaThrust * DlThrust = 
	  new BtaThrust(*DlCandList,*eventInfo,BtaThrust::BTAllParticles); 
	DlThu = DlThrust->thrust();
	DlThuX = DlThrust->thrust_axis().x();
	DlThuY = DlThrust->thrust_axis().y();
	DlThuZ = DlThrust->thrust_axis().z();
	ErrMsg(trace) << "DlThu = " << DlThu << endmsg;
	ErrMsg(trace) << "DlThuX = " << DlThuX << endmsg;
	ErrMsg(trace) << "DlThuY = " << DlThuY << endmsg;
	ErrMsg(trace) << "DlThuZ = " << DlThuZ << endmsg;
	BtaThrust * nonDlThrust =
	  new BtaThrust(*nonDlCandList,*eventInfo,BtaThrust::BTAllParticles); 
	noDlThu = nonDlThrust->thrust();
	noDlThuX = nonDlThrust->thrust_axis().x();
	noDlThuY = nonDlThrust->thrust_axis().y();
	noDlThuZ = nonDlThrust->thrust_axis().z();
	ErrMsg(trace) << "noDlThu = " << noDlThu << endmsg;
	ErrMsg(trace) << "noDlThuX = " << noDlThuX << endmsg;
	ErrMsg(trace) << "noDlThuY = " << noDlThuY << endmsg;
	ErrMsg(trace) << "noDlThuZ = " << noDlThuZ << endmsg;
	delete DlThrust;
	delete nonDlThrust;

	// Calculate angles
	DlThuMag 
	  = sqrt(DlThuX*DlThuX + DlThuY*DlThuY + DlThuZ*DlThuZ);
	ErrMsg(trace) << "DlThuMag1 = " << DlThuMag << endmsg;
	if (DlThuMag < 0.1) DlThuMag=0;
	ErrMsg(trace) << "DlThuMag2 = " << DlThuMag << endmsg;
	noDlThuMag 
	  = sqrt(noDlThuX*noDlThuX + noDlThuY*noDlThuY + noDlThuZ*noDlThuZ);
	ErrMsg(trace) << "noDlThuMag1 = " << noDlThuMag << endmsg;
	if (noDlThuMag < 0.1) noDlThuMag=0;
	ErrMsg(trace) << "noDlThuMag2 = " << noDlThuMag << endmsg;
	denDlnonDl = DlThuMag*noDlThuMag;
	ErrMsg(trace) << "denDlnonDl = " << denDlnonDl << endmsg;
	
     } 

      float cosDlnonDl(-2);
      if (denDlnonDl>0) {
	cosDlnonDl =
	  (DlThuX*noDlThuX + DlThuY*noDlThuY + DlThuZ*noDlThuZ)/denDlnonDl;
      }
      ErrMsg(trace) << "cosDlnonDl = " << cosDlnonDl << endmsg;

      // Thrust angle cut
      if (_requreThrustCut.value() && 
	  fabs(cosDlnonDl)>_cosThrMax.value()) {
	_cutThrust++;
	continue;
      }


      // Fill variables for ntuple
      DhadTypeVec.append(isDstar);
      const int DhadLund(dlCand.theHadron()->pdtEntry()->lundId());
      //const int DstarLund(dlCand.theDstar()->pdtEntry()->lundId());
      DhadModeVec.append(DhadLund);
      DDecayModeVec.append(dDecayMode);
      DMassVec.append(dlCand.theD()->p4().m());
      DhadMassVec.append(dlCand.theHadron()->p4().m());
      const BtaCandidate 
	boostedD(_booster.boostTo(*dlCand.theD()) );
      DMomVec.append(boostedD.p4().rho());
      const BtaCandidate 
	boostedDhad(_booster.boostTo(*dlCand.theHadron()) );
      DhadMomVec.append(boostedDhad.p4().rho());
      DChrgVec.append(dlCand.theD()->charge());
      DhadChrgVec.append(dlCand.theHadron()->charge());
      DlLeptMomVec.append(dlCand.leptonPstar());
      DlLepChrgVec.append(dlCand.theLepton()->charge());
      DlLepTypeVec.append(DlLeptonType);
      DLcosBYVec.append(dlCand.cosThetaBY());
      fitDchiSqVec.append(fitDchiSq);
      fitDnDofVec.append(fitDnDof);
      fitDstatVec.append(fitDstat);
      DprobVec.append(Dprob);
      fitBchiSqVec.append(fitBchiSq);
      fitBnDofVec.append(fitBnDof);
      fitBstatVec.append(fitBstat);
      BprobVec.append(Bprob);
      
      numChargedDlVec.append(nChargedDl);
      numNeutralDlVec.append(nNeutralDl);
      numChargednonDlVec.append(nChargednonDl);
      numNeutralnonDlVec.append(nNeutralnonDl);
      DlThrustVec.append(DlThu);
      DlThrustAxVec.append(DlThuX);
      DlThrustAyVec.append(DlThuY);
      DlThrustAzVec.append(DlThuZ);
      nonDlThrustVec.append(noDlThu);
      nonDlThrustAxVec.append(noDlThuX);
      nonDlThrustAyVec.append(noDlThuY);
      nonDlThrustAzVec.append(noDlThuZ);
      cosDlnonDlVec.append(cosDlnonDl);


      // Get MC truth info
      if(mcTruth) {
	bool isMcTruth = dlCand.fillMcCache(truthMap);
	if (isMcTruth) {
	  // Fill MC truth values for ntuples
	  dlCand.resolveMcAssoc();
	  trueDVec.append(dlCand.realD());
	  trueDhadVec.append(dlCand.realHadron());
	  trueDstarVec.append(dlCand.realDstar());
	  trueDFromBVec.append(dlCand.dFromB());
	  trueDlnuVec.append(dlCand.dlnuExclusive());
	  trueDhadFromBVec.append(dlCand.hadronFromB());
	  trueDstarFromBVec.append(dlCand.dstarFromB());
	  trueLepFromBVec.append(dlCand.leptFromB());
	  trueLepFromSameBVec.append(dlCand.leptFromSameB());
	  trueDhadlnuVec.append(dlCand.leptFromSameB() && dlCand.hadronFromB());
	  trueDstarlnuVec.append(dlCand.leptFromSameB() && dlCand.dstarFromB());
	  directLeptVec.append(dlCand.directLept());
	  mcLeptIdVec.append(dlCand.mcLeptId());
	  mcLeptBIdVec.append(dlCand.mcLeptBId());
	  
	  // Get true B decay mode
	  ErrMsg(trace) << "getting trueDecMod" << endmsg;
	  int trueDecMod(0);
	  int theLeptId = dlCand.mcLeptId();
	  ErrMsg(trace) << "theLeptId = " << theLeptId << endmsg;
	  if (-15<theLeptId && theLeptId<15) { 
	    BtaCandidate* theMcLept = truthMap->mcFromReco(dlCand.theLepton());
	    if (theMcLept == 0) {
	      ErrMsg(routine) << "No true lepton" << endmsg;
	    }
	    if (theMcLept != 0) {
	      ErrMsg(trace) << "found a true lepton" << endmsg;
	      BtaCandidate* aMoth = theMcLept->theMother();
	      if (aMoth == 0) {
		ErrMsg(routine) << "No true B mother" << endmsg;
	      }
	      if (aMoth != 0) {
		int mothID1 = aMoth->pdtEntry()->lundId();
		//int mothID2 = aMoth->uid();
		ErrMsg(trace) << "mothID1 = " << mothID1 << endmsg;
		//ErrMsg(trace) << "mothID2 = " << mothID2 << endmsg;
		if ((-600<mothID1 && mothID1<-500) || (500<mothID1 && mothID1<600)) {
		  int digi1(5), digi2(0), digi3(0), digi4(0);
		  if (mothID1 == 521)  digi1 = 1; //Bp
		  if (mothID1 == -521)  digi1 = 2; //Bm
		  if (mothID1 == 511)  digi1 = 3; //B0
		  if (mothID1 == -511) digi1 = 4;  //B0bar
		  if (digi1<5) {
		    HepAListIterator<BtaCandidate> 
		      iterDaut(aMoth->daughterIterator());
		    BtaCandidate* daut(0);
		    while(daut = iterDaut()) {
		      ErrMsg(trace) << "looping over B daughters" << endmsg;
		      int daughtID = daut->pdtEntry()->lundId();
		      ErrMsg(trace) << "daughtID = " << daughtID << endmsg;
		      if (daughtID == 411)  {digi2 = 1; digi3 = 1;} //Dp
		      if (daughtID == -411)  {digi2 = 1; digi3 = 2;} //Dm
		      if (daughtID == 421)  {digi2 = 1; digi3 = 3;} //D0
		      if (daughtID == -421)  {digi2 = 1; digi3 = 4;} //D0bar
		      if (daughtID == 413)  {digi2 = 2; digi3 = 1;} //Dsrtarp
		      if (daughtID == -413)  {digi2 = 2; digi3 = 2;} //Dsrtarm
		      if (daughtID == 423)  {digi2 = 2; digi3 = 3;} //Dsrtar0
		      if (daughtID == -423)  {digi2 = 2; digi3 = 4;} //Dsrtar0bar
		      if (daughtID == 415)  {digi2 = 3; digi3 = 1;} //D2srtarp
		      if (daughtID == -415)  {digi2 = 3; digi3 = 2;} //D2srtarm
		      if (daughtID == 425)  {digi2 = 3; digi3 = 3;} //D2srtar0
		      if (daughtID == -425)  {digi2 = 3; digi3 = 4;} //D2srtar0bar
		      if (daughtID == 10411)  {digi2 = 4; digi3 = 1;} //D0srtarp
		      if (daughtID == -10411)  {digi2 = 4; digi3 = 2;} //D0srtarm
		      if (daughtID == 10421)  {digi2 = 4; digi3 = 3;} //D0srtar0
		      if (daughtID == -10421)  {digi2 = 4; digi3 = 4;} //D0srtar0bar
		      if (daughtID == 10413)  {digi2 = 5; digi3 = 1;} //D1(2420)p
		      if (daughtID == -10413)  {digi2 = 5; digi3 = 2;} //D1(2420)m
		      if (daughtID == 10423)  {digi2 = 5; digi3 = 3;} //D1(2420)0
		      if (daughtID == -10423)  {digi2 = 5; digi3 = 4;} //D1(2420)0bar
		      if (daughtID == 20413)  {digi2 = 6; digi3 = 1;} //D1(H)p
		      if (daughtID == -20413)  {digi2 = 6; digi3 = 2;} //D1(H)m
		      if (daughtID == 20423)  {digi2 = 6; digi3 = 3;} //D1(H)0
		      if (daughtID == -20423)  {digi2 = 6; digi3 = 4;} //D1(H)0bar
		      if (daughtID == 211)  digi4 = 1;  //pip
		      if (daughtID == -211)  digi4 = 2;  //pim
		      if (daughtID == 111)  digi4 = 3;  //pi0
		    }
		  }
		  trueDecMod = digi1*1000 + digi2*100 + digi3*10 + digi4;
		} else if ((-500<mothID1 && mothID1<-400) || 
			   (400<mothID1 && mothID1<500)) {
		  trueDecMod = 6000;
		} else  if ((-400<mothID1 && mothID1<-100) ||
			    (100<mothID1 && mothID1<400)) {
		  trueDecMod = 7000;
		} else {
		  trueDecMod = 10000;
		}
		ErrMsg(trace) << "trueDecMod before = " 
			      << trueDecMod << endmsg;
	      }
	    }
	  }
	      

	  // Get true D info
	  ErrMsg(trace) << "getting true D info" << endmsg;
	  int trueDDecMod(0), noExtraKPi(0);
	  int mcKID(0), kaon1(0);
	  int mcKMissId(0);
	  int mcDID(0);
	  int KMothID(0), stepsToKMom(0);
	  int mcPi1Id(0), mcPi2Id(0), mcPi3Id(0);
	  int stepsToPi1Mom(0), stepsToPi2Mom(0), stepsToPi3Mom(0);
	  int pi1(0), pi2(0), pi3(0);
	  int mcPiMissId(0);
	  int KPifromSameD(0);
	  int Pi1fromSameD(0), Pi2fromSameD(0), Pi3fromSameD(0);
	  BtaCandidate* theMcPi1 = truthMap->mcFromReco(dlCand.thePi1());
	  BtaCandidate* theMcPi2 = truthMap->mcFromReco(dlCand.thePi2());
	  BtaCandidate* theMcPi3 = truthMap->mcFromReco(dlCand.thePi3());
	  BtaCandidate* theMcK = truthMap->mcFromReco(dlCand.theK());
	  BtaCandidate* aKMoth(0);
	  
	  // Get KId and D mother of the K
	  if (theMcK == 0) {
	    ErrMsg(trace) << "No true kaon" << endmsg;
	  }
	  if (theMcK != 0) {
	    ErrMsg(trace) << "found a true kaon" << endmsg;
	    mcKID = theMcK->pdtEntry()->lundId();
	    ErrMsg(trace) << "mcKID = " << mcKID << endmsg;
	    if (mcKID == 321)  kaon1 = 1; //Kp
	    if (mcKID == -321)  kaon1 = 2;  //Km
	    if (mcKID == 310)  kaon1 = 3; //Ks0
	    if (mcKID == -310)  kaon1 = 4;  //Ks0bar
	    if (kaon1==0) mcKMissId = 1;
	    BtaCandidate* aCand(theMcK);
	    while ((aKMoth = aCand->theMother()) && aKMoth!= 0) {
	      stepsToKMom++;
	      KMothID = aKMoth->pdtEntry()->lundId();
	      ErrMsg(trace) << "KMothID in the loop = " << KMothID << endmsg;
	      if (KMothID == 411)  mcDID = 1; //Dp
	      if (KMothID == -411)  mcDID = 2; //Dm
	      if (KMothID == 421)  mcDID = 3; //D0
	      if (KMothID == -421)  mcDID = 4;  //D0bar
	      if (mcDID!=0) {
		break;
	      } else aCand= aKMoth;
	    }
	  }
	  ErrMsg(trace) << "KMothID after the loop = " << KMothID << endmsg;

	  // Get Pi ID and K Pi from same D or not
	  if (mcDID!=0 && aKMoth!=0) {
	    if (theMcPi1 == 0) {
	      ErrMsg(trace) << "No true Pi1" << endmsg;
	    }
	    if (theMcPi1 != 0) {
	      mcPi1Id = theMcPi1->pdtEntry()->lundId();
	      ErrMsg(trace) << "mcPi1Id = " << mcPi1Id << endmsg;
	      if (mcPi1Id == 211) pi1=1; //Pip
	      if (mcPi1Id == -211) pi1=2; //Pim
	      if (mcPi1Id == 111) pi1=3; //Pi0
	      if (pi1 == 0) mcPiMissId=1;
	      BtaCandidate* aPi1Moth(0);
	      BtaCandidate* aCand(theMcPi1);
	      while ((aPi1Moth = aCand->theMother()) && aPi1Moth!= 0) {
		stepsToPi1Mom++;
		if ( aPi1Moth->uid() == aKMoth->uid()) {
		  Pi1fromSameD=1;
		  break;
		} else aCand= aPi1Moth;
	      }
	      if (theMcPi2 == 0) {
		ErrMsg(trace) << "No true Pi2" << endmsg;
		if (Pi1fromSameD==1) KPifromSameD=1;
	      }
	      if (theMcPi2 != 0) {
		mcPi2Id = theMcPi2->pdtEntry()->lundId();
		ErrMsg(trace) << "mcPi2Id = " << mcPi2Id << endmsg;
		if (mcPi2Id == 211) pi2=1; //Pip
		if (mcPi2Id == -211) pi2=2; //Pim
		if (mcPi2Id == 111) pi2=3; //Pi0
		if (pi2 == 0) mcPiMissId=1;
		BtaCandidate* aPi2Moth(0);
		BtaCandidate* aCand(theMcPi2);
		while ((aPi2Moth = aCand->theMother()) && aPi2Moth!= 0) {
		  stepsToPi2Mom++;
		  if ( aPi2Moth->uid() == aKMoth->uid()) {
		    Pi2fromSameD=1;
		    break;
		  } else aCand= aPi2Moth;
		}
		if (theMcPi3 == 0) {
		  ErrMsg(trace) << "No true Pi3" << endmsg;
		  if (Pi1fromSameD==1 && Pi2fromSameD==1)  KPifromSameD=1;
		}
		if (theMcPi3 != 0) {
		  mcPi3Id = theMcPi3->pdtEntry()->lundId();
		  ErrMsg(trace) << "mcPi3Id = " << mcPi3Id << endmsg;
		  if (mcPi3Id == 211) pi3=1; //Pip
		  if (mcPi3Id == -211) pi3=2; //Pim
		  if (mcPi3Id == 111) pi3=3; //Pi0
		  if (pi3 == 0) mcPiMissId=1;
		  BtaCandidate* aPi3Moth(0);
		  BtaCandidate* aCand(theMcPi3);
		  while ((aPi3Moth = aCand->theMother()) && aPi3Moth!= 0) {
		    stepsToPi3Mom++;
		    if ( aPi3Moth->uid() == aKMoth->uid()) {
		      Pi3fromSameD=1;
		      break;
		    } else aCand= aPi3Moth;
		  }
		  if (Pi1fromSameD==1 && Pi2fromSameD==1 && Pi3fromSameD==1) {
		    KPifromSameD=1;
		  }
		}
	      }
	    }
	  }

	  trueDDecMod = 10000*mcDID + 1000*kaon1 + 100*pi1 + 10*pi2 + pi3;
	  ErrMsg(trace) << "trueDDecMod before = " << trueDDecMod << endmsg;

	  // Check if there are extra kaons or pis
	  int foundExtraK(0), foundExtraPi(0);
	  // Loop over MC truth candidates and store information into ntuple
	  if(nTruth>0 && mcDID!=0 && aKMoth!=0 && theMcK!=0 && theMcPi1!=0) {
	    BtaCandidate* MCcand(0);
	    for ( size_t j=0; j<nTruth; j++ ) {
	      ErrMsg(trace) << "looping over MC cand" << endmsg;
	      MCcand = (*truthList)[j];
	      int MCID = MCcand->pdtEntry()->lundId();
	      ErrMsg(trace) << "ID of this MC cand = " << MCID << endmsg;
	      if (MCID==321 || MCID==-321 || MCID==310 || MCID==-310) {
		if (MCcand->uid() == theMcK->uid()) {
		  ErrMsg(trace) << "found theMcK" << endmsg;
		  continue;
		}
		BtaCandidate* aExtraKMom(0);
		BtaCandidate* aCand(MCcand);
		while ((aExtraKMom = aCand->theMother()) && aExtraKMom!= 0) {
		  int ExtraKMomId = aExtraKMom->pdtEntry()->lundId();
		  ErrMsg(trace) << "ID of ExtKMom = " << ExtraKMomId << endmsg;
		  if ( aExtraKMom->uid() == aKMoth->uid()) {
		    ErrMsg(trace) << "found a extra K" << endmsg;
		    foundExtraK++;
		    break;
		  } else if (ExtraKMomId==11 || ExtraKMomId==-11     //e
			     || ExtraKMomId==13 || ExtraKMomId==-13  //mu
			     || ExtraKMomId==2212 || ExtraKMomId==-2212 //p
			     || ExtraKMomId==211 || ExtraKMomId==-211  //Pi
			     || ExtraKMomId==111                       //Pi0
			     || ExtraKMomId==321 || ExtraKMomId==-321 //K
			     || ExtraKMomId==310 || ExtraKMomId==310  //Ks
			     || ExtraKMomId==130 || ExtraKMomId==-130) {  //KL
		    break;
		  } else aCand= aExtraKMom;
		}
	      }
	      if (MCID==211 || MCID==-211 || MCID==111) {
		if (MCcand->uid() == theMcPi1->uid()) {
		  ErrMsg(trace) << "found thePi1" << endmsg;
		  continue;
		}
		if (theMcPi2!=0 && MCcand->uid()==theMcPi2->uid()) {
		  ErrMsg(trace) << "found thePi2" << endmsg;
		  continue; 
		}
		if (theMcPi3!=0 && MCcand->uid()==theMcPi3->uid()) {
		  ErrMsg(trace) << "found thePi3" << endmsg;
		  continue;
		} 
		BtaCandidate* aExtraPiMom(0);
		BtaCandidate* aCand(MCcand);
		while ((aExtraPiMom = aCand->theMother()) && aExtraPiMom!= 0) {
		  int ExtraPiMomId = aExtraPiMom->pdtEntry()->lundId();
		  ErrMsg(trace) <<"ID of ExtPiMom = "<< ExtraPiMomId << endmsg;
		  if ( aExtraPiMom->uid() == aKMoth->uid()) {
		    ErrMsg(trace) << "found a extra Pi" << endmsg;
		    foundExtraPi++;
		    break;
		  } else if (ExtraPiMomId==11 || ExtraPiMomId==-11     //e
			     || ExtraPiMomId==13 || ExtraPiMomId==-13  //mu
			     || ExtraPiMomId==2212 || ExtraPiMomId==-2212 //p
			     || ExtraPiMomId==211 || ExtraPiMomId==-211  //Pi
			     || ExtraPiMomId==111                       //Pi0
			     || ExtraPiMomId==321 || ExtraPiMomId==-321 //K
			     || ExtraPiMomId==310 || ExtraPiMomId==310  //Ks
			     || ExtraPiMomId==130 || ExtraPiMomId==-130) { //KL
		    break;
		  } else aCand= aExtraPiMom;
		}
	      }
	    } // end of MC truth loop
	    
	    ErrMsg(trace) << "number of extra K = " << foundExtraK << endmsg;
	    ErrMsg(trace) << "number of extra Pi = " << foundExtraPi << endmsg;
	    if (foundExtraK==0 && foundExtraPi==0) noExtraKPi=1;
	  }


	  // Fill values for ntuple
	  trueDecModVec.append(trueDecMod);
	  ErrMsg(trace) << "trueBDecMod = " << trueDecMod << endmsg;
	  trueDDecModVec.append(trueDDecMod);
	  ErrMsg(trace) << "trueDDecMod = " << trueDDecMod << endmsg;
	  noExtraKorPiVec.append(noExtraKPi);
	  mcKPifromSameDVec.append(KPifromSameD);
	  mcKMissIDVec.append(mcKMissId);
	  mcPiMissIDVec.append(mcPiMissId);
	  mcKIDVec.append(mcKID);
	  mcPi1IDVec.append(mcPi1Id);
	  mcPi2IDVec.append(mcPi2Id);
	  mcPi3IDVec.append(mcPi3Id);
	  stepsToKMomVec.append(stepsToKMom);
	  stepsToPi1MomVec.append(stepsToPi1Mom);
	  stepsToPi2MomVec.append(stepsToPi2Mom);
	  stepsToPi3MomVec.append(stepsToPi3Mom);

	} else {
	  trueDVec.append(0);
	  trueDhadVec.append(0);
	  trueDstarVec.append(0);
	  trueDFromBVec.append(0);
	  trueDlnuVec.append(0);
	  trueDhadFromBVec.append(0);
	  trueDstarFromBVec.append(0);
	  trueLepFromBVec.append(0);
	  trueLepFromSameBVec.append(0);
	  trueDhadlnuVec.append(0);
	  trueDstarlnuVec.append(0);
	  directLeptVec.append(0);
	  mcLeptIdVec.append(0);
	  mcLeptBIdVec.append(0);
	  trueDecModVec.append(0);
	  trueDDecModVec.append(0);
	  noExtraKorPiVec.append(0);
	  mcKPifromSameDVec.append(0);
	  mcKMissIDVec.append(0);
	  mcPiMissIDVec.append(0);
	  mcKIDVec.append(0);
	  mcPi1IDVec.append(0);
	  mcPi2IDVec.append(0);
	  mcPi3IDVec.append(0);
	  stepsToKMomVec.append(0);
	  stepsToPi1MomVec.append(0);
	  stepsToPi2MomVec.append(0);
	  stepsToPi2MomVec.append(0);
	}
      } else {
	trueDVec.append(0);
	trueDhadVec.append(0);
	trueDstarVec.append(0);
	trueDFromBVec.append(0);
	trueDlnuVec.append(0);
	trueDhadFromBVec.append(0);
	trueDstarFromBVec.append(0);
	trueLepFromBVec.append(0);
	trueLepFromSameBVec.append(0);
	trueDhadlnuVec.append(0);
	trueDstarlnuVec.append(0);
	directLeptVec.append(0);
	mcLeptIdVec.append(0);
	mcLeptBIdVec.append(0);
	noExtraKorPiVec.append(0);
	trueDecModVec.append(0);
	trueDDecModVec.append(0);
	mcKPifromSameDVec.append(0);
	mcKMissIDVec.append(0);
	mcPiMissIDVec.append(0);
	mcKIDVec.append(0);
	mcPi1IDVec.append(0);
	mcPi2IDVec.append(0);
	mcPi3IDVec.append(0);
	stepsToKMomVec.append(0);
	stepsToPi1MomVec.append(0);
	stepsToPi2MomVec.append(0);
	stepsToPi2MomVec.append(0);
      }

 
      _nSuccesses++;
      nSelectedDlCand++;
     
    } // end of Dlnu candidate loop 
    // #############################################################

    if (nSelectedDlCand==0) {
      _cutCandRequirement++;
      return false;
    }

    if (nSelectedDlCand>150) {
      _cutTooManyDl++;
      return false;
    }


    // ############################################################

    ErrMsg(trace) << "fill ntuple" << endmsg;
    {
      _analysisTuple->column("nMC", nTruth, 0, "Event");
      _analysisTuple->column("nChgTrk", nChgTrk, 0, "Event");
      _analysisTuple->column("nGTVL", nGTVLTrk, 0, "Event");
      _analysisTuple->column("nGTL", nGTLTrk, 0, "Event");
      _analysisTuple->column("nPhoton", nPhoton, 0, "Event");
      _analysisTuple->column("nGamma", nGamma, 0, "Event");
      _analysisTuple->column("isABhabha", isABhabha, 0, "Event");
      _analysisTuple->column("nKaon", nKaon, 0, "Event");
      _analysisTuple->column("nDlCand", nDlCand, 0, "Event");

      _analysisTuple
	->column("nDl", nSelectedDlCand, 0, "DlnuCand", HTRange<int>(0,150));
      _analysisTuple->column("trueD", trueDVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("trueDhad", trueDhadVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("trueDstar", trueDstarVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("trueDFromB", trueDFromBVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("trueDlnu", trueDlnuVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueDlhadFromB", trueDhadFromBVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueDlstarFromB", trueDstarFromBVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueLepFromB", trueLepFromBVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueLepFromSameB", trueLepFromSameBVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("trueDhadlnu", trueDhadlnuVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueDstarlnu", trueDstarlnuVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcDirectLepFromB", directLeptVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("mcLepID", mcLeptIdVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("mcLepBID", mcLeptBIdVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueBDecMod", trueDecModVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("trueDDecMod", trueDDecModVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("noExtraKorPi", noExtraKorPiVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcKPifromSameD", mcKPifromSameDVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcKaonMissID",mcKMissIDVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcPiMissID",mcPiMissIDVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcKaonID",mcKIDVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcPi1ID",mcPi1IDVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcPi2ID",mcPi2IDVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("mcPi3ID",mcPi3IDVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("stepsToKMoth",stepsToKMomVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("stepsToPi1Moth",stepsToPi1MomVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("stepsToPi2Moth",stepsToPi2MomVec , "nDl", 0, "Dlnu");
      _analysisTuple
	->column("stepsToPi3Moth",stepsToPi3MomVec , "nDl", 0, "Dlnu");

      _analysisTuple->column("DhadType", DhadTypeVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DhadMode", DhadModeVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DDecayMode", DDecayModeVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DhadChrg", DhadChrgVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DChrg", DChrgVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DhadMass", DhadMassVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DMass", DMassVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DhadMom", DhadMomVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DMom", DMomVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlLepChrg", DlLepChrgVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlLepType", DlLepTypeVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlLeptMom", DlLeptMomVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DLcosBY", DLcosBYVec, "nDl", 0, "Dlnu");

      _analysisTuple->column("nChrgDl", numChargedDlVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("nNeutDl", numNeutralDlVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("nChrgNonDl", numChargednonDlVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("nNeutNonDl", numNeutralnonDlVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlThrust", DlThrustVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlThrustAx", DlThrustAxVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlThrustAy", DlThrustAyVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("DlThrustAz", DlThrustAzVec, "nDl", 0, "Dlnu");
      _analysisTuple->column("nonDlThrust", nonDlThrustVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("nonDlThrustAx", nonDlThrustAxVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("nonDlThrustAy", nonDlThrustAyVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("nonDlThrustAz", nonDlThrustAzVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("cosDlnonDl", cosDlnonDlVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitDchiSq", fitDchiSqVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitDnDof", fitDnDofVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitDstat", fitDstatVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitDprob", DprobVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitBchiSq", fitBchiSqVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitBnDof", fitBnDofVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitBstat", fitBstatVec, "nDl", 0, "Dlnu");
      _analysisTuple
	->column("FitBprob", BprobVec, "nDl", 0, "Dlnu");
    }
    
    // dumping to ntuple
    _analysisTuple->dumpData();

    _numPassedEvt++;
    return true;
    
}